Initial Condition in Cantera

749 views
Skip to first unread message

Anidro

unread,
Feb 14, 2012, 5:56:39 AM2/14/12
to Cantera User's Group
Several methods require restart from previous computation. This is
true especially if the interest is in approaches not directly
programmed in Cantera, for instance for optimization routines or
continuation methods or even coupling with CFD codes, for which
Cantera would be adopted as a library instead of a separate
application. So I write to ask if in Cantera is possible to import a
previous saved solution as initial condition for advance integrator
method.
I write a Matlab Code to perform a continuation method that use a C++
Cantera Code as integrator but i have some problem to set the initial
state for C++ code.
I hope that somebody could Help me.
Best Regards
Luigi Acampora
PH.D. Student
Seconda Università degli Studi di Napoli

Ray Speth

unread,
Feb 14, 2012, 12:44:45 PM2/14/12
to canter...@googlegroups.com
Luigi,

Yes, in general it is possible to reinitialize a reactor network with a new initial condition and start integrating from there. If you can show your existing code, it would be easier to suggest the appropriate modifications. I'm unclear from your question whether you're using the C++ or Matlab interface to Cantera, and the answer will depend on which Cantera classes you're working with.

Ray

Anidro

unread,
Feb 14, 2012, 1:36:53 PM2/14/12
to Cantera User's Group
#include <cantera/Cantera.h>
#include <cantera/zerodim.h>
#include <cantera/IdealGasMix.h>
#include <iomanip>
#include "wsr.h"

using namespace CanteraZeroD;

OutputStructure wsr(InputStructure input, int i, int j) {
OutputStructure output;
double tign = 0;

// Reaction mechanism
IdealGasMix gas(input.KinMech + ".cti", input.KinMech);
int nsp = gas.nSpecies();

// Create a reservoir for the fuel inlet, and set to pure methane.
Reservoir fuel_in;
gas.setState_TPX(300.0, OneAtm, input.MolecularFormula + ":1.0");
fuel_in.insert(gas);
double fuel_mw = gas.meanMolecularWeight();

// Create a reservoir for the air inlet
Reservoir air_in;
IdealGasMix air("air.cti");
gas.setState_TPX(300.0, OneAtm, "N2:0.78, O2:0.21, AR:0.01");
double air_mw = air.meanMolecularWeight();
air_in.insert(gas);

// To ignite the fuel/air mixture, we'll introduce heat through a
wall.
// We'll create a reservoir to represent the environment and the
wall
// that separates the reactor from the environment.
gas.setState_TPX(300.0, OneAtm, "N2:0.78, O2:0.21, AR:0.01");
Reservoir env;
env.insert(gas);

// Create a reservoir for the exhaust. The initial composition
// doesn't matter.
Reservoir exhaust;
exhaust.insert(gas);

// Create the combustor, and fill it in initially with N2
double* MoleFArray = &input.MoleF[0]; //Converstion Vector to
Array
// input.MoleF is read from csv file containing 53 species mole
fraction evalueted in previous run of software
gas.setState_TPX(input.Temperature, input.Pressure, MoleFArray);
Reactor combustor;
combustor.insert(gas);
combustor.setInitialVolume(input.CombustorVolume);

// Create the wall through which passes the heat flux that allows
the ignition of the mixture
// The wall will use a Guassian 'functor' object to specify the
time-dependent heat flux.
Gaussian heat(input.A, input.t_0, input.FWHM);
Wall w;
w.install(env, combustor);
w.setArea(1.0);
w.setHeatFlux(&heat);

// Compute fuel and air mass flow rates
air_mdot ....
fuel_mdot ....

// Create and install the mass flow controllers. Controllers
// m1 and m2 provide constant mass flow rates
// Fuel mass flow controller.
MassFlowController m1;
m1.install(fuel_in, combustor);
m1.setFunction(&fuel_mdot);

// Now create the air mass flow controller.
MassFlowController m2;
m2.install(air_in, combustor);
m2.setFunction(&air_mdot);

// Put a valve on the exhaust line to regulate the pressure
PressureController v;
v.install(combustor, exhaust);
double Kv = 1e1;
v.setParameters(1, &Kv);
v.setMaster(&m2);

// Define the simulation time.
double tnow = input.ISimTime;
double ttot = 0.0;
double tres, tresnom; // Residence Time and Nominal Residence Time
double phi; // Equivalence Ratio
double Xf_in, Xair_in; // Fuel and Air Inlet Molar Fraction

// The simulation only contains one reactor
ReactorNet sim;
sim.addReactor(&combustor);
sim.setMaxTimeStep(0.5); // Set Max Time Step
sim.initialize(tnow); // Use this to assign the desired initial
time
//sim.setInitialTime(tnow);


// Simulation
double dt = 0.1; // Time step
double nsteps = (input.FSimTime - input.ISimTime)/dt;
for (int n = 1; n <= nsteps; n++) {
tnow += dt;
sim.advance(tnow); // Time to the current step
....
}

return output;

Anidro

unread,
Feb 14, 2012, 1:43:29 PM2/14/12
to Cantera User's Group
I use a Matlab Code that use a .exe written in c++ and based on
Cantera.

Ray Speth

unread,
Feb 15, 2012, 4:57:57 PM2/15/12
to canter...@googlegroups.com
To start an existing reactor network with a new initial condition,
what you should do is set the state on each of the Reactors, using the
Reactor.updateState method.  Its syntax is a little awkward. You pass
an array consisting of:

    (1) the total internal energy (in J)
    (2) the total volume (in m^3)
    (3..N+2) the masses (in kg) of each species

Then call ReactorNet.initialize followed by the advance method.

Hope this helps,
Ray

Anidro

unread,
Feb 16, 2012, 7:29:01 AM2/16/12
to Cantera User's Group
Thanks for your useful reply. I have another question. When I recreate
the reactor how should I set the state of gases to restart the
simulation from previous solution??
Now i use this code:
// input.MoleF is read from .csv file containing 53 species mole
fraction evaluated in previous run of software
gas.setState_TPX(input.Temperature, input.Pressure, MoleFArray);
Reactor combustor;
combustor.insert(gas);

Best Regards
Luigi

On 15 Feb, 22:57, Ray Speth <yarm...@gmail.com> wrote:

Ray Speth

unread,
Feb 16, 2012, 11:46:47 AM2/16/12
to canter...@googlegroups.com
Here is an example that shows how to extract the solution at a certain point in time and use it to restart a reactor network:

    https://gist.github.com/1846315

The key part is that you need to set the state of each Reactor using the 'insert' method, then call the 'initialize' method on the ReactorNet object. I think this is significantly simpler than my previous suggestion.

Hope this helps,
Ray

Anidro

unread,
Feb 17, 2012, 12:59:16 PM2/17/12
to Cantera User's Group
Thanks very much. Your example has been very useful.
Eventually I discovered my mistake.
Best Regards
Luigi
Reply all
Reply to author
Forward
0 new messages