Gas combustion: What is reaction kinetic model (flame stability)

Questions and remarks about code_saturne usage
Forum rules
Please read the forum usage recommendations before posting.
Post Reply
Antech
Posts: 197
Joined: Wed Jun 10, 2015 10:02 am

Gas combustion: What is reaction kinetic model (flame stability)

Post by Antech »

Hello. I have a specific task to check the flame stability of industrial boiler burner. Our standard tool for this is Fluent, but I'd like to have another result to compare so I tried Saturne 7.0.2. As I understand, the only choice with reaction kinetic in Saturne is Libby-Williams model. Generally, it works and it's stable on test case, but I get non-physical result. Although fuel and oxydant temperatures at inlet are both 300 K (<30*C), gas is ignited. In reality, ~500...700*C is needed for gaseous fuel (usually CH4) to ignite. Is it model limitation or I set some parameters wrong? I tried different statistic submodels (2-peaks ... 3-peaks) and Libby-Williams model parameter tweaking in cs_user_parameters.f90, but the mixture always ignite despite of very low inlet temperature.
I attached a complete test case. You only need to run calculation for some minutes and check results with ParaView, there is a PVSM file for convenience. It runs fast because of small mesh. Please see if something can be done to tame the reaction kinetics. Thanks for your attention!
Attachments
CasCombCase.zip
(1.48 MiB) Downloaded 54 times
Antech
Posts: 197
Joined: Wed Jun 10, 2015 10:02 am

Re: Gas combustion: What is reaction kinetic model (flame stability)

Post by Antech »

I tried to introduce kinetic dependence into Saturne 7.0.2 model but without any positive results. As I understand, from three gas combustion models available, only Libby-Williams one deals with kinetics, but it does not work correctly. I investigated the code in extent I can and found the following.

Main part of Libby-Williams model in Saturne is in pdflwc.f90 file. In this file component mass fraction, mixture properties (including temperature) and kinetic reaction rate (form of Arrhenius law) is calculated with or without PDF (statistical distribution function) depending on concentrations statistical dispersions (to save calculation time if dispersion is very low). Common algorithm in this file, applied to each mesh cell, is:
1. Initialize, set pointers...
2. Get mixture enthalpy h in cell:

Code: Select all

h = ((hmax-hmin)*f + hmin*fmax - hmax*fmin)/(fmax-fmin)
where f is a sort of reaction progress variable or reacted fuel fraction (sorry, the code is very difficult to understand, most variables are not commented including doxygen and has not "clear" names like fuelMassPart or so, they named like "f", "y" etc, it's not easy for the user to get what they mean from the code). I suppressed array indexes to make expressions more suitable to get their meaning, the problem is not in PDF / Dirac's function etc, I tried without PDF.
3. Calculate components (fuel, oxydant and products) mass fractions (yfuel etc) in cell:

Code: Select all

yfuel = y
yoxyd = 1 - (coeff3+1)*f + coeff3*y
yprod = 1 - yfuel - yoxyd
yo2  = coeff1 - (coeff1 + coeff2) * f + coeff2 * y
where coeff3 and coeff1 are stoichiometric coefficients or something like this, the problem is not in stoichiometry.
4. Get mixture temperature and density in cell by enthalpy (cothht routine with mode=1).
5. Calculate kinetic reaction rate w in terms of the fuel consumption per second. It's a form of standard Arrhenius law + multiplication of fuel and oxygen (not oxydant) mass fractions, it's very common and OK, the problem is not in kinetic rate itself.
6. Set mass fractions, temperature, density and kinetic reaction rate (fuel mass source) in program-wide arrays.
7. Free memory/pointers etc.

Other important part of the model is in lwctss.f90 file. In this file, fuel concentration field sources are calculated. It's very important that the fuel field is associated with kinetic rate w, but other component fields are not.

The problem, in my understanding, is that product concentration (yprod), enthalpy and temperature are governed not by the reaction rate but by reaction progress variable and, maybe fuel concentration (y and f variables, sorry, it's not clear what they mean). If you disable lwctss subroutine (place return in it's beginning), you will see that reaction progress variable and fuel concentration are solved as transported scalars, but oxydant and product mass fractions are not, they are just a state functions in mesh cells. So even if initial/inlet temperature is very low (30*C), model fuel reacts fine! You cannot affect it by any kinetic parameters in input data. In the same cell (in Paraview), reaction rate (chemical source) tends to zero (10^-13), but temperature increases. The model behaves like it's high kinetic rate while it's zero (anyone can see it in my complete example in first post, just run it). As we know from reaction basics, if diffusion or reaction rate is zero, summary rate must also be zero, but in calculation, in opposite, it's very high.
The mechanism of this, as I can see, is follows. The program calculates progress variable and fuel concentration like transported scalars. Due to turbulent dispersion/mixing we have some intermediate (not 0 or 1) values of these scalars in cells. Then it calculates oxydant and fuel concentrations that are, roughly, arbitrary due to turbulent dispersion (or some balance error arises from reactant mixing/turbulence). Products concentration is calculated by balance: yprod=1-yfuel-yoxyd. So we have some product concentration (example in cell near inlet: ymfuel->0, ymoxyd=0.894, ymprod=0.106, T=447K, w->0; at inlet: Tin=300 K). Then, with reagents mixing downstream, product concentration increases up to unity in some cells and temperature gets higher up to adiabatic (in simple tests I found that enthalpy and temperature depend on yprod). Thus, we have fast reaction with kinetic-diffusion model at 30*C initial temperature that is non-physical. Kinetic part only affects fuel field source and cannot stop this. Even if you set w(idirac)=0 directly in the code you change quite few in test results.

I tried to correct this by setting yprod to zero at low temperatures etc but without luck. There are some artifacts and no correct behavior. The problem is that model based not on complete diffusion-kinetic rate and corresponding component sources like in standard approach but on implicit variables like "progress variable" not governed directly by the real reaction rate. So I didn't find an option to fix this. From my point of view, complete new kinetic-diffusion gas combustion model is required that operates with kinetc-diffusion rate of a given set of reactions and real component fields like CH4 and O2, not progress variable or products mass fractions. Such approach works well in CFX and Fluent, we also used kinetic-diffusion model for coal (not gas) combustion in our engineering programs for decades. An astronomic level is CHEMKIN import and complex solver for detailed mechanisms like GRI-Mech-3.0, but it's extremely hard to implement (Fluent can do it but with problems, not for every task). That's why an optimal solution is simple diffusion-kinetic algorithm with real components fields (CH4, CO, O2, CO2, N2, H2O). Maybe you will introduce it in future? For now, because I can be wrong, please, approve that current model cant take into account kinetic (flammability dependence from temperature) and, maybe, there is a workaround that I didn't find? I cannot introduce complete new combustion model in Saturne, it's outside of my knowledge... I can implement some engineering method described by algebraic relations but not complex things like numerical turbulent diffusion solving, PDFs, flamelets or so...
Li Ma
Posts: 7
Joined: Thu Nov 17, 2022 3:54 pm

Re: Gas combustion: What is reaction kinetic model (flame stability)

Post by Li Ma »

Hello,

Sorry for the response a little bit late. I tested this case rapidly. It is true that the ignite process is not observed whatever the inlet temperature is set. I got some information that may help you from one researcher who used Libby Williams model (maybe the only one at EDF R&D for now).

In fact the model should capture the ignition. However, the model's implementation was done decades ago and has not been maintained. And it aimed at the combustion within gas turbines in thermal power plants at that time. So some configurations and model parameters could be defined only for stable flames for the reason of convenience.

Some info:
f is mixture fraction 0 for oxidizer side, 1 for fuel side
yf is fuel mass fraction, equivalent to progress variable in this case.
And other properties are determined by some state functions.
If you choose 3 peak adiabatic condition for Libby Willams, the pdf used is coded in pdfpp3 instead of pdflwc

One "bug" that I've found is in lwcini which is the initialisation of transported variables

Code: Select all

cvar_yfm(iel) = fmax

cvar_fm(iel)   = fmax
They were initialized as fmax because the gas at the beginning was supposed to be hot for the gas turbine combustion.
But correcting them did not really succeed at capturing ignition.
I will look for further information in the source code and see if it is possible to fix this.

For more complex model that may be help you, the steady laminar flamelets model is implemented in v7.2. But using it is not quite easy, because one should provide flamlete base. And it is only used for diffusion flames. Otherwise, it not possible for premixed flames using the implemented flamelet approach.
We have some works using it for LES simulations with or without progress variables. Flamelet bases are created using FlameMaster with GriMEch 3. If you are interested by this, you can see this article
"Modelling extinction/re-ignition processes in fire plumes under oxygen-diluted conditions using flamelet tabulation approaches"
Antech
Posts: 197
Joined: Wed Jun 10, 2015 10:02 am

Re: Gas combustion: What is reaction kinetic model (flame stability)

Post by Antech »

Thanks for the info! So nostalgic to hear that the model is decades old... The world was different that times.
Although flamelet/PDF models are used by others, we cannot use diffusion-only models because we are interested in boiler burner simulations where ignition never takes place until hot gases get mixed with reactants.
I also noticed that the volume initially filled with fuel but the problem appears in fuel/oxydant boundary layer, not at outer (oxydant/ambient) boundary layer, so it's not the reason.
In my example I use 2-peaks PDF for simplicity, so I experimented with pdflwc.f90. The program reacts on changes but I didn't find the mean to fix ignition issue correctly. I thing about user scalars for gas species + user source according to one-reaction kinetic mechanism but it will be very simplified and time consuming for me because I'm not familiar with the code.
Hope that you will find what causes the problem. If current model can be fixed to capture ignition it would be nice!
Antech
Posts: 197
Joined: Wed Jun 10, 2015 10:02 am

Re: Gas combustion: What is reaction kinetic model (flame stability)

Post by Antech »

Hello.
I experimented with yf (ifmp) setting depending on temperature to make the model sensitive to ignition temperature. I created base results with original model so temperatures was high enough for ignition. Then I introduced the following simple thing in pdflwc.f90 around line 214:

Code: Select all

! ORIGINAL CODE
fmp    =  max(min(fmax,fm(iel)),fmin)
  yfmp   =  max(min(yfm(iel),fmp),                                &
                    zero,(fmp-fs(1))/(1.d0-fs(1)))
! MY ADDITION
  if (cpro_temp(iel).lt.(273.15+650)) then
   if (ipass.gt.1) yfmp=min(1.d0,fmp)
  end if
This suppresses reaction for temperatures under 650*C that is near practical values for CH4 ignition (there is C3H8 in Saturne default stoichiometric file but it's not the problem here). Because yf is a progress variable with starting value of 1.0, we keep the progress in the beginning under low temperatures. Condition with ipass is needed to avoid temperature initialization back to 300 K after calculation restart. If you start from initialization, all things look good: mass fractions and temperature correspond to "cold flow".
I restarted from diffusion results with high temperatures, discretization is SOLU, 2-peak PDF with enthalpy source term. The tweak had an effect but not enough and not as expected. First, high temperatures are still observed in fuel-oxydand internal boundary layer (near the inlet, they may only appear at the outer oxydant jet boundary layer in this case). Second, there is no flow-induced enthalpy transport. At the outer boundary layer of external (oxydant) jet there is a sharp temperature edge instead of the gradient that is characteristic for ejection (attachment) of external hot gases. Also, in internal fuel-oxydant boundary layer, there are no traces after "hot" cells that must be there due to convection of heat (enthalpy). So this tweak is not enough. I attached the temperature profile [in Kelvins]. It looks like, if the cell is hot, it will remain hot independent of convective flow, temperature field doesn't interact with convection despite of that enthalpy source is enabled via gas combustion option.

I tried also to skip cell processing in pdflwc if temperature is low:

Code: Select all

  if (cpro_temp(iel).lt.(273.15+650)) then
   if (ipass.gt.1) return
  end if
but, in this case, it turns back to pure diffusion behavior as with original code.

Please take a look at the code if you have time/intention. Maybe you will find the solution. Thanks for attention to my question.
Attachments
TemperatureWithTweak.png
Antech
Posts: 197
Joined: Wed Jun 10, 2015 10:02 am

Re: Gas combustion: What is reaction kinetic model (flame stability)

Post by Antech »

Hello.
Because my experiments with kinetics in built-in gas combustion model was not successful, I introduced my own standard model that is usually called "Finite Rate Chemistry/Eddy Dissipation". I made it for the company so I cannot share sources here, sorry, but I describe main features.
To get correct and reliable model behaviour under various concentrational and thermal conditions, we need both kinetic and diffusion dependency directly in reaction rate or fuel reacting flow expression. Also, we need to limit amount of reacting fuel to not exceed the fuel mass in mesh cell at current time step so satisfy mass conservation. Technically this can be done via user C functions (extra operation + sources). It's better to store model settings and initial data in separate header file so main files are user sources, user extra operations and initial data header. Steps to implement the model are:
1. Enable species transport and thermal scalar in GUI. Add CH4 (or other fuel), H2O, CO2, N2 and O2 species. Then these species mass fractions can be set on BCs. Add reaction rate (diffusion/kinetic/real) and reacting fuel flow as user arrays (it can be done in Saturne-7 GUI).
2. Create your "user extra operations" function and set all required pointers to current process fields.
3. Calculate maximum available fuel reacting flow in each cell as maxRate=FuelMass/LocalTimeStep. All these quantities are per cell (so calculation is done in cell loop for current process).
4. Calculate diffusion reaction rate with Eddy Dissipation model. It's very simple, just one equation that consist of stoichiometric multiplier, turbulent diffusion rate and component molar concentration [mol/m3] with stoichiometric coefficient. Physically it gives available diffusion flow of reactants (formula with coefficient A) and products (formula with coefficient B) on sub-cell level. Turbulent diffusion between cells are already calculated by solver as species transport. Taking into account product limiter is questionable (I disabled it as it done by default in CFX). Eddy Dissipation model description is easily accessible on the Internet. Don't forget to limit reaction rate as positive (difRate=fmax(rate,0)). I also suggest to always limit all divisors at some minimum value like 10^-50...10^-15 to make model absolutely insensitive to division by zero.
5. Calculate kinetic reaction rate coefficient with Arrhenius equation. Simple model assumes that there is only one global reaction (Fuel+Oxygen=ProductA+ProductB) that requires one kinetic rate. You may use any form of Arrhenius rate from simplest (kinRateCoef=K0*exp(Ea/(R*T))) to complicated like in CHEMKIN/Fluent (manuals are accessible on the Internet). I used the simplest form. Compute kinetic reaction rate as:
kinRate=(FuelMolarConcentration^n)*(OxygenMolarConcentration^m)*kinRateCoef
where n and m are reaction orders of fuel and oxygen. Kinetic parameters and reaction orders are given in famous Westbrook and Dryer article or you can use other sources.
Limit it with positive values as for diffusion one (it's important, otherwise your model will not work at all). Set kinetic rate to zero if temperature is lower than ignition temperature that can be taken as 600...700*C for CH4 (set this temperature in initial data header).
6. Calculate real or resulting reaction rate as minimum from available, diffusion and kinetic (rate=fmin(maxRate,difRate,kinRate)). Please note that if you consider rate as (moles of reactions)/s or as (kg of fuel)/s you need to adjust your reaction rated with corresponding multipliers. If you work with (kg of fuel)/s, maxRate is as is, difRate also has fuel stoichimetric coefficient and molar mass in it's formula to obtain (kg of fuel)/s, but kinRate must be added with FuelMolarMass*FuelStoichiometricCoefficient (kinetic reaction rate is initially in (mols of reactions)/s). Also, pay great attention to pre-exponential factor (K0) dimension! it may be for different dimentional basis (like cm3-mol-s or m3-kmol-s). If needed, compute K0 in desired units in a separate script like Mathcad or in combustion model routine.
7. Apply you reaction rate in user source function. I recommend to store it as user array so you can set it in extra-operations function and then use in user source function. You have reaction rate in some form, so you easily get any component sources with stoichiometric coefficients and molar masses, as usual.
8. Introduce statistics printing like volume part with high temperatures and species mass fractions check for usability. For Saturne, our species are just abstract transported scalars, it will not check if: specie is lower than zero, specie is greater than 1.0, species are not sum to 1.0. So automatic error test makes the model more usable!

Model advantages are:
1. Strict kinetic limiter. If mixture is not heated, it is not ignited, as in reality. It's important for boiler burner simulations.
2. Diffusion limiting with large-scale (inter-cell) and sub-cell diffusion.
3. Simplicity.

Model disadvantages are:
1. Just one global reaction. Usable only for basic gas combustion simulation. For complex chemistry you need implicit multi-reaction solver with reversable reactions, speed-up features and Eddy Dissipation Concept diffusion model. And convenient GUI to set up chemistry part of your simulation. Just text files are absolutely not plausible in this case. It's many times more complex than this simple algorithm.
2. Cannot incorporate precise ignition criteria. To get really precise tool you need indeed complex chemical mechanism like GRI-Mech and some approach to take into account active particles / molecules. Such solver will be very slow even with speed-up features (I tried with Fluent).
Post Reply