Ideal Gas Reactors in Series

[cangt]

Hello All,

I am having a problem with multiple ideal gas reactors in series.

If I run one ideal gas reactor as:
        In Reservoir  -------->   mfc --------> ideal gas reactor --------> out Reservoir                (note: mfc = mass flow controller)

no matter which reaction mechanism I use, it works great. Technically, it is a bit slower or faster depending on the initial conditions that I use (in terms of T, p, etc). This makes sense to me-- It will take time for radicals to build up and the internal solver will need to take smaller time steps accordingly due to stiffness.


The problem comes in when I start running with chains of reactors linked together, as in the example below:

        In Reservoir  -------->   mfc --------> ideal gas reactor --------> mfc / wall --------> ideal gas reactor --------> mfc / wall --------> ideal gas reactor -------->    ...     --------> mfc / wall --------> ideal gas reactor --------> out Reservoir

things go much slower. Again, I believe this makes sense-- since each reactor is tied together in one network, thing will mix (idealized reactor) and send small amounts of other species to other reactors down stream.

The problem here is that things get slower than I would have thought they would a priori. Yes, as they should be slower (the solver will take smaller steps and now must act on all ideal gas reactors, not just one), but it is obnoxiously slow. Even when I test with very small mechanisms (2 steps/5 species) they are so slow that I can't even believe it takes this long. For example, with four ideal gas reactors and four mass flow controllers at the same rate, I usually just stop after minutes of execution since the network time doesn't even pass ~1.0e-4 seconds.

With this being said, here is the code. This is a simplfied version-- it is in C++, and I made it very simple (everything is hard coded).

#include <cantera/thermo.h>
#include "cantera/IdealGasMix.h"
#include "cantera/zerodim.h"
#include <iostream>
#include <fstream>
using namespace Cantera;

int runSimplePSR() {

  double mDot = 1.0;              //The current mass flow rate [kg/s]
  double Vol;                      //Volume of the PSR
  double Tfinal = 1.0;

  IdealGasMix gas("Mechanism.cti");
  IdealGasMix gas2("Mechanism.cti");
  IdealGasMix gas3("Mechanism.cti");
  IdealGasMix gas4("Mechanism.cti");
  IdealGasMix gasIn("Mechanism.cti");

  gas.setState_TPX(1000.0, 30.0*OneAtm, "KERO_LUCHE:1, O2:15.0");
  gasIn.setState_TPX(800.0, 30.0*OneAtm, "KERO_LUCHE:1, O2:15.0");

  IdealGasReactor PSR; IdealGasReactor PSR2; IdealGasReactor PSR3; IdealGasReactor PSR4;

  PSR.insert(gas); PSR2.insert(gas2); PSR3.insert(gas3); PSR4.insert(gas4);

  Vol = 0.01;
  PSR.setInitialVolume(Vol); PSR2.setInitialVolume(Vol); PSR3.setInitialVolume(Vol); PSR4.setInitialVolume(Vol);

  Reservoir envIn; Reservoir envOut;

  envIn.insert(gasIn); envOut.insert(gas);
 
  Wall wall; Wall wall2; Wall wall3; Wall wall4;

  wall.install(PSR, envIn); wall2.install(PSR, PSR2); wall3.install(PSR2, PSR3); wall4.install(PSR3, PSR4);

  wall.setArea(1.0);
  wall.setExpansionRateCoeff(1.0e6);
  wall2.setArea(1.0);
  wall2.setExpansionRateCoeff(1.0e6);
  wall3.setArea(1.0);
  wall3.setExpansionRateCoeff(1.0e6);
  wall4.setArea(1.0);
  wall4.setExpansionRateCoeff(1.0e6);

  MassFlowController mfcIn;
  MassFlowController mfcOut;
  MassFlowController mfc1;
  MassFlowController mfc2;
  MassFlowController mfc3;

  mfcIn.install(envIn, PSR);
  mfc1.install(PSR, PSR2);
  mfc2.install(PSR2, PSR3);
  mfc3.install(PSR3, PSR4);
  mfcOut.install(PSR, envOut);

  mfcIn.setMassFlowRate(mDot);
  mfc1.setMassFlowRate(mDot);
  mfc2.setMassFlowRate(mDot);
  mfc3.setMassFlowRate(mDot);
  mfcOut.setMassFlowRate(mDot);
 
  ReactorNet net;

  /* Insert the PSR into the network */
  net.addReactor(PSR);
  net.addReactor(PSR2);
  net.addReactor(PSR3);
  net.addReactor(PSR4);

  std::ofstream reactorStream;
  reactorStream.open("CanteraResults.dat");
 
  double t  = 0.0;
  double dt = 1.0e-3;

  reactorStream << "time      temperature      pressure     density    volume     MFR     resTime    ";
  for (int i=0; i<gas.nSpecies(); i++) reactorStream << gas.speciesName(i) << "    ";
  reactorStream << std::endl;

  while (net.time() < Tfinal) {   //Time less than 1 seconds

    ThermoPhase& conts = PSR.contents();       
    std::cout << net.time() << "  " << PSR.temperature() << "  " << PSR.pressure() << "  " << PSR.density() << "  " << PSR.volume() << "   " << "   " << mDot << "   " << PSR.density()*PSR.volume()/mDot  << std::endl;      
    reactorStream << net.time() << "  " << PSR.temperature() << "  " << PSR.pressure() << "  " << PSR.density() << "  " << PSR.volume() << "   " << "   " << mDot << "   " << PSR.density()*PSR.volume()/mDot << "   " ;       
    for (int i=0; i<gas.nSpecies(); i++) reactorStream << conts.moleFraction(i) << "  ";
    reactorStream << std::endl;
    net.step(dt);                
    t += dt;                
  }
   
  return 0;

}

int main() {

  /* Try to run the PSR */
  try {
    runSimplePSR();
  }
  /* This *should* handle ANY exception thrown by Cantera */
  catch (CanteraError& err) {
    std::cout << err.what() << std::endl;
  }

  return 0;
 
}
...

I have tried it numerous ways:

1. With one gas in all of the reactors (didn't make a difference)
2. Not outputting the data from the first reactor (didn't make a difference)
3. Changing properties so that the residence time (tau = r * vol / rho) is higher / lower (made minimal difference)

What I am wondering is this: I read that a Newton solver is the method used. I'm assuming when adding more reactors to the system, the Newton solver now has the algebraic equations for all of the reactors. Is there a way to change the solver type? If I used net.advance() as opposed to net.step(), would that make a difference? Basically, what I am looking for is any way to speed this up! It doesn't look like Cantera supports running in parallel (ie, one reactor on each core or processor) since there is a data dependency. Additionally, I'm not sure how much time it would save by doing this... each gas only contains five species and two possible reactions in my case.

Any help would be appreciated!!!!!!!! Thanks everybody for any advice!!!!!!

[Nick Curtis]

Hey,

One thing you should fix, replace:

net.step(dt)

with:

net.step(Tfinal)

the step function takes a single (internal) integration step towards the time passed in.

Since you are passing in dt, the network time will never pass dt

Alternatively if you want regular output you can do:

double t = 0;
//while loop etc.
net.advance(t + dt);
t += dt;

Which will integrate in steps of dt.

However, in your case the integration should still be going to 1e-3s before stopping.

I think it's likely that the expansion coefficient of the wall is inducing huge amounts of stiffness; a 1 Pa delta P would result in a wall velocity of 10^6 m/s! and the volume of the reactors is 1m^3!

Are you sure you didn't mean something a bit more reasonable there?

Best,

Nick

[cangt]

Hi Nick,

Thanks so much for the response, I'm banging my head around-- this is much appreciated!

As for imposing areas:

As I understood it, I set the area to a wall area of 1 m^3. I know my reactor volumes are much lower, but I have no feel for this sort of thing! Where does the delta P come in? I was just trying to use the wall as a means to hold the pressure constant. I tried valves but couldn't really get them to work well. Is there a better way do do this? Is it just a matter of changing wall areas/coefficients?

What would you suggest are good wall areas and good expansion coefficients? If this can help the stiffness of the system, that would be great! Like I said before, I would do anything to get it to run faster. I didn't realize I was even setting a deltaP-- I thoguht it would just stay constant since the wall coeffiencts were so high... there may be something I'm missing!

[Ray Speth]

Hi,

I think the first thing you need to fix is how your MassFlowControllers are connected: mfcOut should presumably connect PSR4 to envOut, not PSR to envOut. You also probably want to set the initial state of the other IdealGasMix objects so that the initial pressure is the same in each reactor.

The issue with the expansion coefficients isn’t that you’re setting a delta-P, it’s that a delta-P always exists during the transient (since the pressure is an independent variable for each reactor) and that there is a time constant associated with equalizing that that is inversely proportional to the expansion rate coefficient.

That said, I think you want to avoid walls entirely here, since using them means that your steady state will end up with each reactor having some arbitrary volume based on how the walls move during the initial transient. Instead, what you should use is a MassFlowController connecting the inlet to the first PSR, then use PressureController to connect each successive PSR. For each PressureController, you need to call the “setMaster” to associate it with the mass flow controller, and setParameters to set the value of the pressure coefficient. For example:

PressureController pc1;
pc1.setMaster(&mfcIn);
double pressureCoeff = 0.1; // kg/s/Pa
pc1.setParameters(1, &pressureCoeff);
pc1.install(PSR, PSR2);

Regards,
Ray

[cangt]

Hi Ray,

Thanks for your input, it is much appreciated!

I haven't fooled around with pressure controllers, I always used mass flow controllers in combination with walls. As I understand it, it will do the same thing. Point taken about the wall changing during transients, however the coefficient is so large that I didn't think it would make a difference... What do you thin, Ray??

Something interesting that you pointed out though: the outlet reservoir was connected to the first PSR, not the last, so it wasn't truly in series. When I changed it (to the correct version), it ran much, MUCH faster. That brings me to an interesting point: The code I wrote above was written very quickly and not as I wrote the original code-- the code above is completely hardcoded. When I wrote the original, it actually contained a matrix of reactors, gases, etc, not just a single row. Additionally, I use numerous pointers in the original code so that memory is allocated on the fly based on the number of reservoirs, reactors, etc. that I want at the point of program execution. I am wondering now if a pointer got messed up somewhere during memory allocation and is pointing somewhere else (to the wrong place). If so, it would have the same affect as I had before, where reactors are connected incorrectly and causing the time to increase by an order of magnitude (or more). This is something that I need to look into!

If anybody has any general comments on how to get things to run faster, please send them my way!! I thank you both, Nick and Ray for the insight!!

[cangt]

In case anybody is following, there does NOT seem to be any memory issues! Still not sure what the issue is!

[Ray Speth]

Hi,

Using walls will have the same effect as the flow controllers in that it will give you a constant pressure, but it is not the same because the change in the reactor volumes when using walls means that you aren't solving the same system, and so your solution will be different. You might be misinterpreting the meaning of the coefficient -- the larger the value, the faster the wall moves in response to a pressure differential. In either case, the advice about setting the coefficients to smaller values to avoid introducing very fast time scales applies.

I would suggest printing out as diagnostics the mass and volume of each reactor after each timestep and making sure they are behaving as you expect.

If you are trying to simulate a significantly larger network of reactors, it's worth noting that the overall solve time for the system scales as the number of variables cubed, which is set by the need to factorize the Jacobian.

Regards,

Ray

[cangt]

Hey Ray,

Thanks for all of the insight!

I may try switching from walls to pressure controllers, thanks! I actually do print out everything on each network solver step-- the pressure stays constant and volume changes (as the density changes due to temperature changes from reactions). This all made sense to me. It seems like it works, so I never bothered changing it. However, if the wall coefficients introduce numerical stiffness that the pressure controllers won't introduce, it is definitely worth changing to pressure controllers to decrease runtime.

As for the solver-- I am indeed looking at large networks, which is why the runtime is so important to me. If moving from walls to pressure controllers will decrease runtime, then it is worth a code overhaul. As I understand, all of the PSRs go into one network. The network is solved using Newton's method (Newton's method, which has quadratic convergence but requires O(n^3) operations per iteration). Do you know if there are any other solvers available? Is there any way to switch? Any other tips or tricks to achieving faster runtime?

[Ray Speth]

Hi,

High values for the pressure controller coefficients will also introduce short time scales, so this is an issue no matter what method you use for keeping the pressure equal in each reactor.

Are you interested in the transient solution or just the steady state? In the steady case, you can simulate the chain of reactors one at a time, using the steady-state solution for each reactor as the upstream reservoir state for the next reactor. This is the approach taken in the example surf_pfr.py, which you could adapt for your purposes. For the transient case, you could do something similar, but it would be more complex -- you would need to save the output mixture state as a function of time, and then in the next integration, use that to update the state of the upstream mixture after each integrator timestep.

The only solver available in Cantera is the dense, fully implicit mode of CVODES, which requires the O(n^3) Jacobian factorization.

Regards,

Ray