Counter flow diffusion flame 1D simulation

[Alan Kong]

Hi all,

Sorry for all the recent questions. This time, it has something to do with counter flow diffusion flame simulations.

I am intending to run a 1D flame simulations with my customized detailed mechanism and later my reduced mechanism.

Firstly, I run successfully the example python file for counter flow diffusion flame, diffusion_flame.py, noticing that it takes only a few tries before Cantera is successful with the Newton time stepping solver.

Next, I tried testing the same code with my mechanism (207 species, 2329 reactions) and applying the same grid points and refinement setting as with the example files.

The only change I made was to comment out this line set_boundary_emissivities as Cantera could not recognized this attribute of CounterFlowDiffusionFlames.

I am using Pure CH4 as fuel and Pure O2 as oxidizer, performing at pressure of 60 bar and inlet temp of 300k. Grid size is the same as diffusion_flame.py.

However, the Newton time stepping solver keep running for a long time (2 hours) and fail to find a proper solutions to solve the 1D flame calculations.

I do not know if my grid is too coarse or that I need to change the tolerance setting for the steady state problems and time stepping. I am also confused as to how the size of the mechanism affects the solver.

Any help on this matter will be great!

Thanks in advance!

Alan

[Thomas Fiala]

Hi Alan,

I once faced a similar problem. Counterflow flame simulations are quite stiff and they are very prone to diverge unless you have a good initial solution.
The best way to get to your desired operating conditions is to slowly iterate successful counterflow flame simulations to it. Ideally, you should adapt the grid in the process to avoid time stepping failures and to speed up the process. I found out a set of optimal parameters for doing so (which are explained and published here), and created a Cantera tutorial, which is in the repository for Cantera 2.2. You can already copy it from here.
The script uses a strain rate computation also introduced in the development version. If you don't need that, you can comment it out and run it in Cantera 2.1.

The procedure generally is:

1. Create a good solution with your desired fuel/oxidizer combination at a low strain rate and 1 bar. If Cantera already diverges while computing that, try to first use the fuel/oxidizer combination from a working example (like `diffusion_flame.py`) and slowly adjust your species composition and temperature. Also, ensure that the entire flame is inside the domain and not to close to either inlet (i.e. you see parabolic profiles of the axial velocity u at both oxidizer and fuel side). Additionally, the flame should be far from extinction, i.e. the peak temperature is close to the adiabatic flame temperature.

2. Use the pressure loop from `diffusion_flame_batch.py` to get to your desired pressure. If wanted, you can then apply the strain rate loop to compute counterflow flames at various strain rates.

I hope this will help you. If you have any questions, feel free to ask!

Thomas

[Alan Kong]

Hi Thomas,

Is Cantera 2.2 already released? I am using the latest version of 2.1.2.

I notice that the time stepping fails most of the time. For initial calculation, the solver only finds the correct time stepping after 6-7 attempts at 20 time step and step size of 1.e-05. (I am in the process of twerking this setting)

I am actually quite lost right now with what to setting to change first, grid points etc. but I will give you methods a try.

Does changing the composition and temperature affect the results greatly? I am using CH4:1 and O2:1, both using the same kinetic mechanisms. I was going to try stoichiometric ratio of CH4:1 and O2:2.

Is the results only possible for certain inlet temperature for fuel and oxidizer? In this sense, I am quite limited to a range of temperature.

The reason I am asking this is that I am using both the mole fraction profile of certain species and temperature profile as some of the criteria for validation with experiment data and these constraints may lead to result that do not match those data, thus if possible, I will like to know what kind of changes it may do to the results

I am not sure what you meant by ' slowly iterate successful counterflow flame simulations to it', is it about iterating and redo the simulations with successive parameters and setting that works?

I am new to the area of combustion and especially in the field of diffusion flame and the type of reactors and experiment conducted for diffusion flames (counterflow, co flow etc.)

My final aim, if time permits, is to modify the code and implement a co flow reactor which is the most closest to my operating conditions so if there are any informations, website, books as such that can help me understand, that will be great.

And thanks so much Thomas for your help!

Cheers

Alan

[Thomas Fiala]

Hi Alan,

No, Cantera 2.2 is not yet released. However, you can check out the development version and compile it.

For literature, I can recommend the following books and review paper:
R. J. Kee, M. E. Coltrin, and P. Glarborg. Chemically Reacting Flow.
C. K. Law. Combustion Physics.
H. Tsuji. Counterflow diffusion flames.
S. R. Turns. An introduction to combustion: concepts and applications.

If I understand you correctly, you want to know whether it makes a difference if you set `'O2:1'` or `'O2:2'` as inlet compositions. The answer to that question is no, because mole and mass fractions will always be normalized and unless you specify some other species (as, for example, N2) the oxidizer composition will always be just oxygen. You can specify the 'global equivalence ratio' by the mass fluxes at the boundaries. However, keep in mind that diffusion flames always burn around the stoichiometric surface, and the 'global equivalence ratio' is essentially meaningless. The ratio between the fuel and oxidizer mass fluxes only controls whether the flame burns closer to the fuel or the oxidizer side of the domain. You specify the counterflow flame by the (normalized) fuel and oxidizer compositions and temperatures, pressure, and strain rate (for the latter, several definitions exist, see also the new function in Cantera 2.2 or my paper).

By 'slowly iterating counterflow flame simulations' I basically meant:

1. Compute a counterflow flame simulation which you know converges, like from one of the examples. E.g. you use the one from `diffusion_flame.py`, the volumetric composition of the oxidizer is about 21% O2 and 79% N2. After computing it, keep your Python terminal open (or use a script for the following steps).
   
2. Now you want to have pure oxygen as oxidizer. To get to it, set `f.oxidizer_inlet.X` to `'O2:30,N2:70'`. Then, again compute your simulation with `f.solve()`.
3. Keep going until you reach `'O2:1,N2:0'`.
4. In the mean time, check that your flame does not touch the boundaries. I use matplotlib to plot the temperature and velocity profiles to assure that this does not happen.
   

I hope that this is clearer now.

Best regards,

Thomas

[Alan Kong]

Hi Thomas,

I have tried your method with relatively good success although the scripts takes pretty long to compute.

For my mechanism, (207 species and 2329 reactions), it took 8 hours to compute my results and I have still not done the full pressure range as defined in the code.

That said, I have so question.

1) How different it is to use a linear space pressure range instead of a log based pressure range, does it affect the computation time?

2) For this case, are we doing the simulation in 1 instant of time or does this following some rules or constraints ?

3) The result (from some modification of your code) came out well but I will like to introduce my own initial grid system after taking into consideration where the flame is likely to be from the initial guess. Well this help in the computation time?

4) For the iteration, can I implement the change in composition in a loop like what u did for pressure? However, I ran the code once and got some result with the CH4:1, O2:1 composition for fuel and oxidizer. Are these results correct?

5) I am intending to change the values of mass flow velocity and the grid distance to fit some experimental setup, will this affect result greatly? (keeping in mind of using my original initial solution obtained as a guide to set the boundary condition)

Lastly, it has something to do with how the Counterflowdiffusionflame() solver works. How does the result get calculated? 

For example, in the case of solving for one pressure, does the solver do several iteration for one pressure value before saving the result? What about the initial newton time stepping, do we use those time step to calculate the end result?

In the case of the premixed flame simulations, I can specify the reaction to advance 1 step but I dont know if I can do the same with the counterflowdiffusionflame module.

I asked my supervisor who only worked with CHEMKIN and he told me that the solver maybe calculating the propagation of the flame as it meets, until the temperature profile reaches a certain threshold before saving the result in f.T.

All in all, I received great help from you and I will be using parts of your code and literature for my thesis. How do I properly cite your work in my code.

Really appreciate your help in this!

Alan

[Thomas Fiala]

Hi Alan,

I'm glad that you found my suggestions helpful.

Your long computation times can be due to an unfavorable initial grid (I sometimes had problems with ultra-low strain rates), but I think in your case it probably because of the huge mechanism you use. It might be worth using a simpler mechanism like GRI to get to your operating conditions, and then use this solution as the initial solution for your very detailed mechanism.

Regarding your questions:

1. Sure, you can pick any range you like. As many values scale exponentially, a logarithmic scale might be the fastest to reach a certain pressure. However, I have not tested what is the best step size. You could go ahead and try coarser steps. If the solution fails, just reduce the step size.
2. I don't really get what you mean by this question - all simulation results are ideally steady-state results.
3. It definitely would be best to use initial grid and mass flux settings optimized for your composition. This probably will help to speed up computation time.
4. Sure, that should be the best way to do it. The question whether your solution is "correct" is not so easy to answer. If Cantera computes a solution, you should have a physical solution for your boundary conditions (with respect to the underlying physical models). However, if you set an unfavorable combination of grid and boundary conditions, you might have a result different to what you expect it to be due to interaction of the flame with the boundary conditions.
5. Again, this depends: If your flame is fully contained in both your numerical and experimental domain, i.e. the boundaries are so far away that you effectively have potential flow boundary conditions, the flame structure should be identical. If this is the case, then you can also "inflate" your computational domain by extrapolating the mass fluxes from your solution to meet the experimental values.

I'm not confident enough to explain in detail how the solver works. It's probably best to have a look at the C++ source code yourself. Also, I think the book by Kee, Coltrin, and Glarborg mentions some schemes which should be similar.

If you found the batch procedure and the underlying scaling rules useful, it surely would be nice if you mentioned the corresponding paper;)

Best regards,

Thomas

[Alan Kong]

Hi All,

Been really successful with the batch calculation and using it for later iteration of calculations.

However, recently, I have been seeing this error popping out during solving

alloc: invalid block: 000000000730A9D0: 80 7

and another time it was

alloc: invalid block: 0000000005767300: 50 4

and then python.exe crash. (This only happens at the end of the calculations)

My guess right now is that there is an overload on the memory used by Cantera, causing python to crash. However, I am not very confident that this is the true issue.

There was also times where the solution is unable to converge in the middle of the batch calculations for 1 pressure setting even though it was successfully for other pressures (which I find weird as I have quite a good initial guess from previous calculations which work find)

CanteraError thrown by OneDim::timeStep:

Time integration failed

I will be very grateful if some wise men here can elaborate on what has happened here.

Cheers

Alan

[Thomas Fiala]

Hi Alan,

Do you execute the matplotlib routines at the end? That might be the source of your 'alloc' error. Matplotlib might fail while handling huge sets of data. However, I haven't experienced that error before, so I'm also not sure.

Cantera is very sensitive to all boundary and operating conditions, and sometimes the solution fails for one particular constellation. For me, it typically helped to just change one property by a tiny amount. Employing the scaling rules solved this problem for me completely. However, your large reaction mechanism might need a little more care.

Thomas

[Alan Kong]

Hi Thomas,

The scaling rule work just fine for the my mechanism. I will try to write a separate scripts for plotting to check on the alloc error. Thanks for the heads up. If the error continue to appear, I will switch to matlab for plotting purposes.

I am also wondering if its possible to get a time dependent solutions, similar to the one used in sim1D.advance where we can specify the duration of time for the simulation with the counterflow flame. I dont know if this is helpful but I am trying to plot the evolution of the Fuel and exhaust product profile with time (mole fractions of CH4, O2, H2O and CO2). I see that there is another code by Ray Speth, 'Ember' which is capable of doing that but just wondering if its possible to do it in Cantera too.

Ultimate goal of all this is the generate a reduce mechanisms and I have started my first round of reduction and trying out a few mechanisms already (RAMEC, GRI MECH3.0, etc) which all work really well with your scaling law.

Cheers

Alan

[Thomas Fiala]

Hi Alan,

Good to hear that most things work fine.

I never really looked into the solution algorithm, but I assume that you won't be able to get the transient behavior out of Cantera, because Cantera tries to solve a steady-state problem and the pseudo time stepping is not necessarily physical time stepping (Ray or anyone, please correct me if I'm wrong). Ember is especially designed for unsteady simulations, so it is probably the way to go.
If you just want to monitor the solution behavior, you could consider increasing the log level, or writing data at interrupts. An example for the latter would be:

def interrupt_extinction(t):
    plt.plot(f.grid, f.T)
    return 0.
f.set_interrupt(interrupt_extinction)

(This probably is a bad example as it will surely blow up matplotlib with a huge number of lines)

Best regards,

Thomas

[Alan Kong]

Hi Thomas,

Just a quick question. The scaling law proposed, does it keep the profile of T and mass fraction constant while changing the grid, u profile, V profile, mass flow rate and lambda profile according to the desired pressure and strain rate?

Does it make sense to apply the scaling only to the mass flow rate but retain the grid size and re-evaluating the u, V and lambda profile?

Cheers

Alan

[Thomas Fiala]

Hi Alan,

Yes, the profiles of temperature and mass fractions are held constant with respect to the relative grid position. As the grid size is varied, the "width" of these profiles are scaled, whereas their peak values remain untouched.

Keep in mind that the scaling is there just for the initialization of the new solution. Therefore, after solving the problem, the solution should be identical for constant operating and boundary conditions (that is, pressure, strain rate, inlet temperature/composition, tolerances).

The scaling rules are set up to quickly compute a batch of counterflow flame simulations (up to extinction), and also to iterate the solution to operating conditions which fail to converge with wrong initial conditions. Depending on your situation, it certainly might make sense to only vary some parameters, e.g. to match experimental boundary conditions. It just depends on what you want to do;)

Best regards,

Thomas

[Alan Kong]

Hi Thomas,

I have a slight problem. I used the scaling law to compute a solution up to the pressure of 60 bar, keeping strain rate constant and using the a ~ p^-3/2

However, I had arrive at a grid which 0.1 mm wide which is a little too small for my case of 20mm. (The inlet temp, composition are at my operating condition)

I also find that the mass flow rate are too large for my application (more then 10 times my initial input for mdot f and mdot o)

When I tried to move towards my condition, by just changing the mass flow by 0.1, the solution fail to converge. 

I will try putting all my experimental setting at the start and try to see if I move closer to my experimental setting. 

Thanks for the tips

Alan

[Thomas Fiala]

Hi Alan,

From the solution you computed, you can try to reduce the strain rate like in the strain rate loop of the example (just with a factor < 1). This will give you a wider flame and lower inlet mass fluxes.

The relevant part of the flame (with respect to the strain rate) will be captured by the simulation in any case. Left and right of the reacting zone, there will be just fuel and oxidizer, respectively. Their axial velocity profiles are parabolic.
If you want to model this entire zone of your experiment, you can manually extend the grid on either side and fill it up with fuel or oxidizer. Don't forget to adjust your mass flow rates accordingly.

Best regards,

Thomas