0D reactor IdealGasConstPressure

[Alan Kong]

Hi all,

Recently, I am trying to generate a series of Reaction Path Diagram for a 0D IdealGasConstPressure reactor for my mechanism of around 200 species and 2000 plus reactions.

My settings are CH4:1, O2:2, OH:0.01 with initial temp 300K and I did the same loop as the code reaction_path.py, stopping the temperature at steady state combustion. (i.e about 3644K for my simulation)

I have been running into trouble trying to introduce an igniter for the problem (a small 1% stream of OH radical into the reactor at 300K). CVODES kept having ths:

I am not sure if the integrator of CVODES is actually working after this as my program still runs and still generate a reaction path diagram. However, upon checking the.dot files, I released it gave me some really strange diagram where my fuel CH4 is not even included in the starting node of the diagram. I have started printing the temperature for each interaction and found that this occurred around 315.0913896821292K. The diagram that I have now looks really different from the simulation I did where the starting temp is near ignition temp (from 520K and stop at 3644, no OH radical added) which shows me a more plausible path diagram. Is there anyone who can help me.

Another strange thing that I discover is that I was about to generate a proper reaction path diagram with the setting CH4:1, O2:2, initial temp 300K, stop at 1000K few weeks back without even having to introduce the OH radical as a igniter (kind of strange but it works at that time)

I was able to generate a diagram for C, O and H element. These calculation was done in Cantera 2.1.2 while the calculation above is done in Cantera .2.2. I have check that it is not the problem with the version of Cantera as I tested the code with both version and I get the same problem.

I had also tried with putting the setting  CH4:1, O2:2, H:0.01 for the inital composition but still, I get the CVODES error. 

It appears that I need another way of introducing the igniter radical.

If someone here encountered the same problem, please enlighten me! Any help will be great,

Cheers

Alan

[Ray Speth]

Alan,

It's hard to know exactly what's going on without seeing your code. In addition, providing the mechanism file as well would allow others to try reproducing the behavior you see.

If you're just simulating a closed reactor as in the reaction_path.py example, at equilibrium you of course will have very little CH4 left in your reactor, so it's not surprisng for it to be excluded from the diagram.

To simulate a steady reactor with an inflow and outflow, I would recommend initializing this with the reactor contents at the constant pressure equilibrium, rather than using a stream of radicals. I know there is a different example that uses this concept, but it's problematic numerically, since it introduces some very short timescales because the radicals are so far from their equilibrium concentrations when introduced into a gas that's at low temperature.

While you should be concerned about the CVODES warnings printed above, note that they are not errors, and since the integrator eventually returns, you are getting a solution to the problem, although it's possible there's some loss of accuracy associated with those warnings.

Thank you for testing the new version of Cantera. I especially appreciate reports of any unexpected differences between the versions.

Regards,

Ray

[Alan Kong]

Hi Ray,

I have attached my code and mechanism files for you reference.

Something strange happened again. Using the same code, I tried running it but did not get any CVODES warnings this time. This happened after I did a test run with the step() function going to a large time (180s).

I also tried what you suggested by initializing the igniter with a separate reactor at 1500K and composition of OH:0.01 and joining the 2 resulting reactor together. However this doesnt work as there was no change to the initial temp.

I have switched back to Cantera 2.2 and I am using python 3.4. The resulting diagram I get when taking the reaction path at 1000K is definitely wrong as it only shows 1 reaction path.

I had also run the code with initial temp at 500K which works fine without the igniter, however I do know that if the initial temp is different, the resulting net flux in the diagram will be affected greatly, therefore, I have to be extra careful with that.

Cheers

Alan

I have attached also, the diagram that I produced successfully without any igniter and at initial temp of 300K. Mole fraction CH4:1, O2:2. Snapshot at 1000K

I dont know if consecutive code testing with the command prompt open will affect the result? 

[Alan Kong]

[Nick Curtis]

Alan,
I think what Ray was suggesting was to have a flow reactor that is initialized at the HP equilibrium state, and then have stoichiometric inflow, e.g.

gas = ct.Solution('C1C4.xml')
gas.TPX = 300.0, 60.0e5 , 'CH4:1, O2:2, OH:0.01'  #stoichiometric condition
gas.equilibrate('HP')
r = ct.IdealGasConstPressureReactor(gas)


gas2 = ct.Solution('C1C4.xml')
gas.TPX = 300.0, 60.0e5 , 'CH4:1, O2:2, OH:0.01'  #stoichiometric condition
res = ct.Reservoir(gas)

env = ct.Reservoir(gas)

mdot = ct.MassFlowController(res, r, mdot = 0.01)

pres_regulator = ct.PressureController(r, env, master=mdot)

net = ct.ReactorNet([r])


Then simulate and do your reaction path diagram

This way, your stoichiometric flow will be ignited by the equilibrium products (which should be closer to the final state than the radial stream)

The reason your separate igniter didn't work is
1) that it was downstream from the base reactor, there is no upstream flow in Cantera to my knowledge
2) It doesn't actually appear to be connected through a flow device (I'm not sure what this implies in Cantera)

Nick


On Monday, June 1, 2015 at 8:59:36 AM UTC-4, Alan Kong wrote:

[Ray Speth]

Alan,

Rather than jumping straight to a reaction path diagram, you might want to try looking at the species concentrations that are actually coming out of that simulation. The result that I see is that at the point that you reach T=1000, the solution contains significant negative mass fractions for the species C3H5-S and H. This can happen when both species have small negative mass fractions (which are acceptable within the tolerances of the numerical integrator) but the reaction mechanism contains reactions between these species (which your mechanism does) resulting in negative creation rates for these species, and a resulting positive feedback loop. You can mitigate this problem by tightening the integrator tolerances. I found the following to be sufficient in this case:

net.rtol = 1e-11
net.atol = 1e-20

If you do this, what you will see is the behavior that I would expect -- the mixture will essentially never ignite. What happens is that the radical species you introduced just immediately recombine, increasing the temperature of the mixture slightly, but not nearly enough to reach a temperature at which the mixture will autoignite in any reasonable amount of time (I stopped the simulation after 1e7 seconds). I'm not sure I understand why you are trying to use high radical concentrations as a way of igniting the mixture, instead of by increasing the temperature, which is the physically-meaningful approach in all of the cases I'm aware of.

Regards,

Ray

[Alan Kong]

Hi Ray,

Thanks for taking the time to read the codes.

I am trying to introduce the radical at 1500K and at an amount of 1% to the fuel. (H or OH) This is to simulate an igniter spark in my case of application where I am trying to see how the reaction proceed from cryogenic condition. (which seem right now not plausible. This was what I initially expected when I first run the initial code. Seem odd that I was even able to generate the diagram in the first place as seen in the attached diagram.)

In fact, I am also trying to get the ignition delay time for the simulation at 0D and hence I am also plotting the time vs temperature profile which I am able to produce by taking fixed time step using the step() function in a loop but I start at 450K inital temperature.

I cant quite get my head around how the initial temperature affects the result in the reaction path diagram, maybe I had to try with different plot at different starting temp to see the effect. (I have tried with starting temp 500k and 550k, both produced similar diagram with identical major pathway, but the value in net flux different slightly.)

All in All, my main objective is to identify the ignition delay time from the temp profile (time to point where gradient dT/dt is largest) and produce reaction path diagram at 100-200K interval from the starting temp to steady state ignition to see how the reaction pathway change before, after and during ignition. (This will eventually lead to a short slideshow in my presentation) 

I had tried with a temp near the ignition temp (found from the time vs temp plot), and I am able to produce really accurate reaction pathways which will then aid me in my reduction of the full mechanism. (at least I think I do)

I would like to know how do we go about increasing the temperature (by using r.T and using a 'for' loop to increase the T in uniform interval until the ignition is self sustaining (stop after 550K)) 

Also, I am more interested in the case of non premixed co flow and counterflow flame and I know that it is currently not possible to conduct sensitivity analysis for non premixed flame in Cantera.

Sorry for the long post but currently am a little bit lost with the software and trying to see what can be done and what cannot be achieve at this stage.

Thanks for the great help so far.

Cheers

Alan

[Ray Speth]

Alan,

I think the idea of trying to represent a spark igniter is not consistent with a homogeneous reactor. If what you're interested in is the ignition process when starting at a particular condition, then you need to simulate that condition directly. If you add radicals or heat to the system, you aren't solving the same problem anymore. 

You're solving a system of nonlinear differential equations involving many variables. It's not sufficient to think about just the temperature -- there is important chemistry taking place during the initial induction period, even though the temperature does not increase significantly. There will certainly be similarities in the reaction paths extracted at the point where T=1000 (or any other value) for different starting temperatures, but each of these cases is a unique IVP, and you should expect them to be different.

I'm not sure why you think sensitivity analysis can't be done for the nonpremixed flame. The procedure outlined in the flamespeed_sensitivity.py example should work just fine, with the exception that you need to pick a different output to take the sensitivity of, such as the total heat release rate, since the flame speed isn't a meaningful concept in this case.

Regards,

Ray

[Alan Kong]

Thanks Ray for the help and advice.

I think my next step may be to run the simulation with fixed time step from 300K and then pick the point in which the temp gradient is the highest.

Am I right to say that there is still chemistry going on at 300K between the fuel and oxidizer? 

Will it be better if I simulate my case as a combustor with mass flow controller for  fuel and oxidizer and see if it ignites?

Lastly, I believe I have ask this question before but not sure it had been answered, can I use ReactionPathDiagram on 1D nonpremixed flames?

Cheers

Alan

[Alan Kong]

And also a short question on equilibrate()

This function shows the end of the chemical reaction right? Meaning, the amount of CH4 is reduce to a level that the forward and back reaction rate is equal and there is no consumption of CH4 anymore.

Cheers

Alan

[Ray Speth]

Alan,

According to the model, there is still chemistry happening at 300K, although it is extremely slow. You will probably run into numerical problems trying to simulate this case however, since it involves such a wide range of timescales that it is difficult to represent accurately using finite-precision arithmetic. At the point where this is a problem, though, the model isn't realistic anyway -- the approximation of an adiabatic reactor might be reasonable over a period of a few seconds, but not when the computed ignition delay time is measurable in years.

A reactor with an inlet and outlet flow is a different model. Whether it is "better" or not depends on what you're trying to model. This is a more complex system to analyze, as this is a system that for some boundary conditions (i.e. inlet composition and flow rate) will have multiple steady-state solutions, and which one you reach depends on the initial conditions (i.e. the initial state of the reactor volume).

The ReactionPathDiagram class takes as input a Solution object, not a flame. You can compute a reaction path diagram for any point in the flame, but not (directly) the integral over the entire flame domain. You could do the integration yourself by generating ReactionPathDiagrams for each point in the flame and parsing the output of ReactionPathDiagram.get_data(). This would give you all of the integral fluxes, but you'd have to build the graph yourself, unfortunately.

Equilibrate will show you the end state of the reactions if all of your reactions are reversible and your mechanism contains pathways that connect all of the species. I noticed in your mechanism that many reversible reactions are written as two irreversible reactions in opposite directions. Whether or not this leads to the same end state as equilibrate will depend on whether the rate constants are consistent with the species thermodynamic data, since specifying the reverse rate constant explicitly like this is actually overconstrained.

Regards,

Ray

[Alan Kong]

Hi Ray,

I am still confused with the method u suggest by initializing the reactor contents at the constant pressure equilibrium.

Do u mean by using a Reservoir with the reactor contents at constant pressure equilibrium and introduce a flow inlet (mdot = 1 to ensure stoichiometric condition) to my IdealGasConstPressure Reactor which is at 300K. Subsequently, I will generate the ReactionPathDiagram in the Reactor.

Lastly, can I generate an empty reactor or is it mandatory to include a Solution object in the reactor.

Cheers

Alan

[Alan Kong]

Hi Ray,

The C++ documentation and the Python Documentation are not the same I presume? So if I find anything in C++ documentation, I may not be able to use it in Python.

Cheer

Alan

[Ray Speth]

Alan,

My suggested about initializing the contents at constant pressure equilibrium was based on the assumption that you were trying to simulate a well stirred reactor with inlet and outlet flows, since you initially mentioned a "stream" of radicals. That suggestion does not apply to simulating a closed reactor with no inputs and outputs.

You can create a Reactor without initially giving it a Solution object. However, you must then call Reactor.insert(gas_object) before trying to integrate the reactor network.

Regards,

Ray

[Alan Kong]

Hi Ray,

In the code combustor.py, what does the factor = 0.1 imply? 

Cheers

Alan

[Ray Speth]

Alan,

As you can see from the immediately following lines, it's just a scaling factor applied to the fuel and oxidizer mass flow rates:

factor = 0.1
air_mdot = factor * 9.52 * air_mw
fuel_mdot = factor * equiv_ratio * fuel_mw

The idea is that you can change 'factor' to change the residence time while holding the equivalence ratio constant.

Regards,

Ray

[Alan Kong]

Hi Ray,

I am pretty new to this but what is this 'residence time'? is it similar to the induction time, the time delay before ignition?

Cheers

Alan

[Alan Kong]

Hi Ray,

I have something more of an abstract question regarding the ignition of methane.

It is harder to ignite a CH4/O2 then CH4/air? I did a similar reaction like the one in combustor.py with pure CH4 and O2 but it seem that I have to increase the amplitude of the gaussian mass flow rate of the igniter  quite a bit to achieve ignition.

cheers

Alan

[Steven C. Decaluwe]

For flow reactors, the residence time is a measure of the ratio of the reactor volume to the volumetric flow rate.  This isn’t exactly right because the volumetric flow rate can change inside the reactor, but in general the name is actually pretty descriptive - it is a measure of the average time that a unit of mass spends inside the reactor, before flowing out.  Higher residence time, obviously, means more time for reactions to proceed.

Hope this is at least somewhat helpful,

Steven

----------------------------------------------------------
Steven DeCaluwe, Ph.D.
Assistant Professor of Mechanical Engineering
Colorado School of Mines
G. Brown Hall, W410 B
Golden, CO, 80401

email: deca...@mines.edu
ph: (303) 273-3666
Twitter: @DrDeCaluwe

--
You received this message because you are subscribed to the Google Groups "Cantera Users' Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to cantera-user...@googlegroups.com.
To post to this group, send email to canter...@googlegroups.com.
Visit this group at http://groups.google.com/group/cantera-users.
For more options, visit https://groups.google.com/d/optout.

[Ray Speth]

Alan,

The residence time is a property of the reactor system that you set up, not an intrinsic property of the mixture. In the combustor.py example, you have a well-stirred reactor with an inflow and an outflow. The residence time is the typical time that a unit of gas spends in the reactor. It can be computed as the mass of the gas in the reactor divided by the mass flow rate through the reactor. If you have a book on combustion that covers the well-stirred reactor model, it should talk about this idea in more detail.

Regards,

Ray

[Alan Kong]

Hi Ray,

Just a question on the previous suggested method of using the equilibrate gas content to ignite the gas at low temperature in a 0D IdealGasConstPressure reactor.

It is save to assume that the equilibrate gas is gas that are ignited by some external means and by introducing these 'ignited gas', it simulate the condition where part of the ignited gas is used to heat up the unheated gas in the combustion chamber?

Cheers

Alan