IV - 1D computations¶
The two following links might help you during the computation of 1D flames :
https://cantera.org/documentation/docs-3.0/sphinx/html/matlab/one-dim.html
https://cantera.org/science/flames.html
The first one is about the function that can be used for computing a one dimensional flame and the second one explains the equations that are computed.
The main difference between Cantera and a solver like AVBP is that equations are already one dimensional in Cantera.
1. ADIABATIC FLAME - A freely-propagating, premixed flat flame¶
This script will show you the creation of a premixed flame. First the initial solution is created and then, the calculation is performed before plotting the interesting results.
import cantera as ct
import numpy as np
from matplotlib.pylab import *
1.1 Initial solution¶
Import solution¶
gas = ct.Solution('gri30.yaml') # Import gas phases with mixture transport model
Set the parameters¶
# General
p = 1e5 # pressure
tin = 600.0 # unburned gas temperature
phi = 0.8 # equivalence ratio
fuel = {'CH4': 1} # Methane composition
oxidizer = {'O2': 1, 'N2': 3.76} # Oxygen composition
Set gas state to that of the unburned gas¶
gas.TP = tin, p
gas.set_equivalence_ratio(phi, fuel, oxidizer)
Create the flame and set inlet conditions¶
f = ct.FreeFlame(gas, width=0.02) # Create the free laminar premixed flame specifying the width of the grid
f.inlet.X = gas.X # Inlet condition on mass fraction
f.inlet.T = gas.T # Inlet condition on temperature
1.2 Program starts here¶
The idea here is to solve four flames to make the calculation converge. The first flame will be solved without equation energy whereas the others will continue the computation by enabling the calculation using the equation energy and by adding more point to the flame front.
First flame¶
f.energy_enabled = False # No energy for starters
tol_ss = [1.0e-5, 1.0e-9] # tolerance [rtol atol] for steady-state problem
tol_ts = [1.0e-5, 1.0e-9] # tolerance [rtol atol] for time stepping
f.flame.set_steady_tolerances(default=tol_ss)
f.flame.set_transient_tolerances(default=tol_ts)
# Set calculation parameters
f.set_max_jac_age(50, 50) # Maximum number of times the Jacobian will be used before it
# must be re-evaluated
f.set_time_step(1.0e-5, [2, 5, 10, 20, 80]) # Set time steps (in seconds) whenever Newton convergence fails
f.set_refine_criteria(ratio=10.0, slope=1, curve=1) # Refinement criteria
# Calculation
loglevel = 1 # amount of diagnostic output (0 to 5)
refine_grid = 'disabled' # 'refine' or 'remesh' to enable refinement
# 'disabled' to disable
f.solve(loglevel, refine_grid) # solve the flame on the grid
f.save('4-Output/ch4_adiabatic.yaml', 'no energy','solution with no energy') # save solution
Second flame¶
%%capture
f.energy_enabled = True # Energy equation enabled
refine_grid = 'refine' # Calculation and save of the results
f.set_refine_criteria(ratio=5.0, slope=0.5, curve=0.5) # Refinement criteria when energy equation is enabled
f.solve(loglevel, refine_grid)
f.save('4-Output/ch4_adiabatic.yaml', 'energy', 'solution with the energy equation enabled')
print('mixture-averaged flamespeed = {0:7f} m/s'.format(f.velocity[0])) # m/s
Third flame and so on ...¶
%%capture
f.set_refine_criteria(ratio=2.0, slope=0.05, curve=0.05) # Refinement criteria should be changed ...
f.solve(loglevel, refine_grid)
f.save('4-Output/ch4_adiabatic.yaml', 'energy continuation','solution with the energy equation enabled continuation')
points = f.flame.n_points
print('mixture-averaged flamespeed continuation = {0:7f} m/s'.format(f.velocity[0])) # m/s
print('mixture-averaged final T = {0:7f} K'.format(f.T[points - 1])) # K
Fourth flame and so on ...¶
%%capture
f.transport_model = 'Multi' # Switch transport model
f.solve(loglevel, refine_grid)
f.save('4-Output/ch4_adiabatic.yaml', 'energy_multi',
'solution with the multicomponent transport and energy equation enabled')
points = f.flame.n_points
print('multicomponent flamespeed = {0:7f} m/s'.format(f.velocity[0])) # m/s
print('multicomponent final T = {0:7f} K'.format(f.T[points - 1])) # K
You may notice that there are small differences on the flamespeed and on the final temperature between the two transport models, as it was predictable.
1.3 Manipulate the results¶
Save your results¶
f.save('4-Output/ch4_adiabatic.csv') # Write the velocity, temperature, density, and mole fractions
# to a CSV file
Plot your results (temperature, density, velocity, ...)¶
rcParams['figure.figsize'] = (14, 10)
# Get the different arrays
z = f.flame.grid
T = f.T
u = f.velocity
ifuel = gas.species_index('CH4')
fig=figure(1)
# create first subplot - adiabatic flame temperature
a=fig.add_subplot(221)
a.plot(z,T)
title(r'$T_{adiabatic}$ vs. Position', fontsize=15)
xlabel(r'Position [m]')
ylabel("Adiabatic Flame Temperature [K]")
a.xaxis.set_major_locator(MaxNLocator(10)) # this controls the number of tick marks on the axis
# create second subplot - velocity
b=fig.add_subplot(222)
b.plot(z,u)
title(r'Velocity vs. Position', fontsize=15)
xlabel(r'Position [m]')
ylabel("velocity [m/s]")
b.xaxis.set_major_locator(MaxNLocator(10))
# create third subplot - rho
c=fig.add_subplot(223)
p = zeros(f.flame.n_points,'d')
for n in range(f.flame.n_points):
f.set_gas_state(n)
p[n]= gas.density_mass
c.plot(z,p)
title(r'Rho vs. Position', fontsize=15)
xlabel(r'Position [m]')
ylabel("Rho [$kg/m^3$]")
c.xaxis.set_major_locator(MaxNLocator(10))
# create fourth subplot - specie CH4
d=fig.add_subplot(224)
ch4 = zeros(f.flame.n_points,'d')
for n in range(f.flame.n_points):
f.set_gas_state(n)
ch4[n]= gas.Y[ifuel]
d.plot(z,ch4)
title(r'CH4 vs. Position', fontsize=15)
xlabel(r'Position [m]')
ylabel("CH4 Mole Fraction")
d.xaxis.set_major_locator(MaxNLocator(10))
# Set title
fig.text(0.5,0.95,r'Adiabatic $CH_{4}$ + Air Free Flame at Phi = 0.8 Ti = 600K and P = 1atm',fontsize=22,horizontalalignment='center')
subplots_adjust(left=0.08, right=0.96, wspace=0.25, hspace=0.25)
2. Compute your own premixed flame¶
Please go on the website :
https://chemistry.cerfacs.fr/en/home/
You will find the reference website of cerfacs as far as chemistry is concerned. If you go to Chemical Database and Data table, you will notice that several mechanisms can be used for your computation.
- compute two premixed flames : one including the gri30.yaml and the other the Lu.yaml (under methane/skeletal on the website).
- the flames will be fuel/air at an equivalence ratio of 1.1, a cold gases temperature of 400K and a pressure of 2 bars.
- Try to post-process your data by plotting the flame speed of both flames.
The solution you are supposed to get is shown below.
NB : the transport used should be a mix one.
Here is where you should code your premixed flame¶
import cantera as ct
import numpy as np
from matplotlib.pylab import *
names = ['gri30.yaml', 'Mechanisms/Lu.yaml']
fuels = ['CH4','CH4']
fig=figure(1)
for i,n in enumerate(names):
gas = ct.Solution(n) # Import gas phases with mixture transport model
p = 2e5 # pressure
tin = 400.0 # unburned gas temperature
phi = 1.1 # equivalence ratio
fuel = {fuels[i]: 1} # Fuel composition
oxidizer = {'O2': 1, 'N2': 3.76} # Oxygen composition
gas.TP = tin, p
gas.set_equivalence_ratio(phi, fuel, oxidizer)
f = ct.FreeFlame(gas, width=0.02) # Create the free laminar premixed flame specifying the width of the grid
f.inlet.X = gas.X # Inlet condition on mass fraction
f.inlet.T = gas.T # Inlet condition on temperature
f.energy_enabled = False # No energy for starters
tol_ss = [1.0e-5, 1.0e-9] # tolerance [rtol atol] for steady-state problem
tol_ts = [1.0e-5, 1.0e-9] # tolerance [rtol atol] for time stepping
f.flame.set_steady_tolerances(default=tol_ss)
f.flame.set_transient_tolerances(default=tol_ts)
# Set calculation parameters
f.set_max_jac_age(50, 50) # Maximum number of times the Jacobian will be used before it
# must be re-evaluated
f.set_time_step(1.0e-5, [2, 5, 10, 20, 80]) # Set time steps (in seconds) whenever Newton convergence fails
f.set_refine_criteria(ratio=10.0, slope=1, curve=1) # Refinement criteria
# Calculation
loglevel = 1 # amount of diagnostic output (0 to 5)
refine_grid = 'disabled' # 'refine' or 'remesh' to enable refinement
# 'disabled' to disable
f.solve(loglevel, refine_grid) # solve the flame on the grid
f.energy_enabled = True # Energy equation enabled
refine_grid = 'refine' # Calculation and save of the results
f.set_refine_criteria(ratio=5.0, slope=0.5, curve=0.5) # Refinement criteria when energy equation is enabled
f.solve(loglevel, refine_grid)
f.set_refine_criteria(ratio=2.0, slope=0.05, curve=0.05) # Refinement criteria should be changed ...
f.solve(loglevel, refine_grid)
rcParams['figure.figsize'] = (14, 10)
# Get the different arrays
if n == 'gri30.yaml':
z_gri3 = f.flame.grid
u_gri3 = f.velocity
else:
z_lu = f.flame.grid
u_lu = f.velocity
# create second subplot - velocity
b=fig.add_subplot(111)
b.plot(z_gri3,u_gri3)
b.plot(z_lu, u_lu)
title(r'Velocity vs. Position', fontsize=25)
xlabel(r'Position [m]', fontsize=15)
ylabel("velocity [m/s]", fontsize=15)
b.xaxis.set_major_locator(MaxNLocator(10))
plt.show()
fig.savefig('4-Output/Comparison.png')
This introduces the notion of types of mechanisms. To represent your chemistry, you will have :
- detailed mechanisms such as the gri30 which are very complete but also containing a lot of species (so hard to integrate into solvers such as AVBP).
- global mechanisms that are very simple (usually 6 species and 2 equations) but that may be inaccurate.
- ARC mechanisms that are an inbetween solution (among others) and that is chosen to be used at CERFACS. The Lu mechanism used here is an intermediate step to access to an ARC mechanism and thus shows some differences with the detailed, as you might see on the graph.
Nevertheless, be always careful when you use this or this detailed mechanism as a reference for your computations. Indeed, it might be accurate with a certain error bar or developed for special conditions that are not yours.
3. A burner-stabilized rich premixed methane-oxygen flame at low pressure¶
This script will show you the creation of a burner. This is basically the same thing as a premixed flame, except that the flame is stabilized on a burner and not freely proapagating.
From R.J. Kee, M.E. Coltrin and P. Glarborg Chemically Reacting Flow
- the gri30 mechanism (CH4/air flame)
- with a pressure of 1 bar, a temperature of 373 K, an equivalence ratio of 1.3
# Import gas phases with mixture transport model
gas = ct.Solution('gri30.yaml')
# Parameter values :
# General
p = 1e5 # pressure
tin = 373 # unburned gas temperature
phi = 1.3
fuel = 'CH4: 1'
oxidizer = 'O2:1.0, N2:3.76'
# Set gas state to that of the unburned gas
gas.TP = tin, p
gas.set_equivalence_ratio(phi, fuel, oxidizer)
f = ct.BurnerFlame(gas, width=0.2)
f.burner.T = gas.T
f.burner.X = gas.X # Conditions
mdot = 0.04
f.burner.mdot = mdot
Here, the strategy is, as above, a four step calculation. For the purpose of this training, we have commented lines so that the result is only done in one calculation.
As you might notice, the calculation works and it converges, but it takes more time than doing it by running four calculations.
Therefore, the choice is entirely yours and dependent on the fuel you are using (heptane will take more time to converge than methane for instance).
It is possible that your flame does not converge in some simulations (sometimes, the Newton solver used to compute the flame is crashing if the chemistry is too sharp - meaning that the chemical timesteps are very small - or for other reasons). You have several parameters that you can play on to make it work :
- increase or decrease the tolerances for steady-state and time stepping (whether absolute or relative).
- increase your ratio in the refine criteria.
- take a slower time step with more steps.
NB : If you do not want the solver to plot the lines of computing, just uncomment the %%capture command.
NB2 : Note that two functions might be useful for your future developments :
- if you want to save a solution in a xml file, you can use the
save()function on the freeflame object. - if you want to restore a solution from a xml file, you can use the
restore()function. This is very useful if you want to analyse old results or if you want a flame to converge without starting from a linear initialization.
NB3 : To specify the grid refinement criteria, it is possible to set the function set_refine_criteria with correct values. This enables to have more control on the way the solver refine the curve. The different fields definitions are :
- ratio : additional points will be added if the ratio of the spacing on either side of a grid point exceeds this value
- slope : maximum difference in value between two adjacent points, scaled by the maximum difference in the profile (0.0 < slope < 1.0). Adds points in regions of high slope.
- curve : maximum difference in slope between two adjacent intervals, scaled by the maximum difference in the profile (0.0 < curve < 1.0). Adds points in regions of high curvature.
- prune : if the slope or curve criteria are satisfied to the level of prune, the grid point is assumed not to be needed and is removed. Set prune significantly smaller than slope and curve. Set to zero to disable pruning the grid.
For a first computation, it is best to put the ratio value to a sufficiently high value; and to let the 'slope' and 'curve' values to 1. Use 'prune' if you see that Cantera tries to refine your grid too much, and that the number of points starts to exceed a reasonable value (approximately 200 points is enough).
#%%capture
#################
#f.energy_enabled = False
tol_ss = [1.0e-5, 1.0e-9] # [rtol atol] for steady-state
tol_ts = [1.0e-5, 1.0e-4] # [rtol atol] for time stepping
f.flame.set_steady_tolerances(default=tol_ss)
f.flame.set_transient_tolerances(default=tol_ts)
f.set_refine_criteria(ratio=10.0, slope=1, curve=1)
f.set_max_jac_age(50, 50)
f.set_time_step(1.0e-5, [1, 2, 5, 10, 20])
loglevel = 1 # amount of diagnostic output (0 to 5)
#refine_grid = 'refine' # True to enable refinement, False to
#f.solve(loglevel, refine_grid)
#f.save('4-Output/ch4_burner_flame.xml', 'no_energy',
# 'solution with the energy equation disabled')
#################
f.energy_enabled = True
f.set_refine_criteria(ratio=3.0, slope=0.1, curve=0.2)
#f.solve(loglevel, refine_grid='refine')
#f.save('4-Output/ch4_burner_flame.xml', 'energy',
# 'solution with the energy equation enabled')
#################
f.transport_model = 'multicomponent'
#f.solve(loglevel, refine_grid='refine')
#f.save('4-Output/ch4_burner_flame.xml', 'energy_multi',
# 'solution with the energy equation enabled and multicomponent transport')
#################
f.soret_enabled = True
f.solve(loglevel, refine_grid='refine')
f.save('4-Output/ch4_burner_flame.yaml', 'energy_soret',
'solution with the energy equation enabled and multicomponent transport')
import matplotlib.pyplot as plt
plt.plot(f.grid, f.velocity)
plt.xlabel('grid (m)',fontsize=15)
plt.ylabel('Laminar flame speed ($m/s$)',fontsize=15)
plt.title('Laminar flame speed vs x-axis for an equivalence ratio of 1.3', fontsize=20)
plt.show()
f.save('4-Output/ch4_burner_flame.csv')
NB : This type of flame is not used in practice as most of the combustion chambers are based on the study of premixed flames configurations.. Nevertheless, it can be useful for some special configurations that you may want to study.
Next, we will see the counter-flow diffusion flame as this is an important configuration to study when studying a combustion chamber. However, if you need to look for other types, some more are available (see above).
4. An opposed-flow ethane/air diffusion flame¶
A counter-flow diffusion flame is basically a flame where the oxidizer and the fuel meet to burn (they are not mixed before burning, contrary to the two others types that we have seen before).
From Pfitzner, Michael & Müller, Hagen & Ferraro, Federica. (2013). Implementation of a Steady Laminar Flamelet Model for nonpremixed combustion in LES and RANS simulations.
First of all, the solution object is instanciated as usual. The mass flow rate of the oxidizer and the fuel have been fitted so that the flame is in stoichiometric conditions. In addition to the composition, the temperature and the pressure, the density and the mass flow rate are given. By dividing the second by the first, the velocity of each gases is given which leads to the calculation of the strain rate, an important parameter to counter-flow diffusion flame.
gas = ct.Solution('gri30.yaml', transport='mixture-averaged')
p = 1e5 # pressure
comp_f = 'C2H6:1' # fuel composition
tin_f = 300.0 # fuel inlet temperature
rho_f = p / (8.314 / 0.030 * tin_f) # fuel inlet density
mdot_f = 0.24 # fuel inlet mass flow rate (kg/m^2/s)
vel_f = mdot_f / rho_f # fuel inlet velocity
print('Velocity of the fuel : ' + str(vel_f))
comp_o = 'O2:0.21, N2:0.78, AR:0.01' # oxidizer composition
tin_o = 300.0 # oxidizer inlet temperature
rho_o = p / (8.314 / 0.029 * tin_o) # oxidizer inlet density
mdot_o = 0.72 # oxidizer inlet mass flow rate (kg/m^2/s)
vel_o = mdot_o / rho_o # oxidizer inlet velocity
print('Velocity of the oxidizer : ' + str(vel_o))
gas.TP = tin_o, p
width = 0.02
a = (vel_o + vel_f) / width # calculation of the strain rate (s^-1)
print('Strain rate of the diffusion flame : ' + str(a))
This will instanciate the counterflow diffusion flame. As you might notice, and this is valid for all the other cases, you can instanciate the initial grid yourself. This allows you to refine the grid when you want. This can be a good trick when your simulation is crashing and you want to compute it again.
As far as the inlets are concerned, the composition, temperature and mass flow rate needs to be specified.
initial_grid = np.linspace(0, width, 6)
f = ct.CounterflowDiffusionFlame(gas, initial_grid)
# Set the state of the two inlets
f.fuel_inlet.mdot = mdot_f
f.fuel_inlet.X = comp_f
f.fuel_inlet.T = tin_f
f.oxidizer_inlet.mdot = mdot_o
f.oxidizer_inlet.X = comp_o
f.oxidizer_inlet.T = tin_o
This starts the calculation of the counterflow diffusion flame. By experience, you will see that, out of this wonderful tutorial case, diffusion flames are hard to converge. You can either play with the size of the domain and every parameter said all along the tutorial (indeed, sometimes it converges better if the domain is smaller), or you can try a method that is developed here in CERFACS that resolves the diffusion flame in the z-space.
If you want some more information about that, come and ask the chemistry team of cerfacs.
#%%capture
# First flame:
# disable the energy equation
f.energy_enabled = False
# Set error tolerances
tol_ss = [1.0e-5, 1.0e-11] # [rtol, atol] for steady-state problem
tol_ts = [1.0e-5, 1.0e-11] # [rtol, atol] for time stepping
f.flame.set_steady_tolerances(default=tol_ss)
f.flame.set_transient_tolerances(default=tol_ts)
# and solve the problem without refining the grid
f.solve(loglevel=1, refine_grid='disabled')
#################
# Second flame:
# specify grid refinement criteria, turn on the energy equation,
f.energy_enabled = True
f.set_refine_criteria(ratio=3, slope=0.8, curve=0.8)
# and solve the problem again
f.solve(loglevel=1, refine_grid='refine')
# save your results
f.save('4-Output/c2h6_diffusion.yaml', 'energy')
#################
# Third flame:
# specify grid refinement criteria, turn on the energy equation,
f.energy_enabled = True
f.set_refine_criteria(ratio=2, slope=0.2, curve=0.2, prune=0.04)
# and solve the problem again
f.solve(loglevel=1, refine_grid='refine')
# save your results
f.save('4-Output/c2h6_diffusion.yaml', 'energy continuation')
#################################################################
# Save your results if needed
#################################################################
# write the velocity, temperature, and mole fractions to a CSV file
f.save('4-Output/c2h6_diffusion.csv')
These are interesting quantities to plot, such as the different mass fractions and the temperature of the flame. This is all done in the following graph (with the values being adimensionalised).
# Get interesting values
z = f.flame.grid
T = f.T
u = f.velocity
# Get interesting indices for computation of species
fuel_species = 'C2H6'
ifuel = gas.species_index(fuel_species)
io2 = gas.species_index('O2')
in2 = gas.species_index('N2')
# Initiate interesting vectors
c2h6 = np.zeros(f.flame.n_points,'d')
o2 = np.zeros(f.flame.n_points,'d')
hr = np.zeros(f.flame.n_points,'d')
# Computes interesting quantities for analyzing a counter-flow flame
for n in range(f.flame.n_points):
f.set_gas_state(n)
c2h6[n]= gas.Y[ifuel]
o2[n]= gas.Y[io2]
hr[n] = - np.dot(gas.net_production_rates, gas.partial_molar_enthalpies)
fig=figure(1)
a=fig.add_subplot(111)
a.plot(z,T/np.max(T),z,c2h6/np.max(c2h6),z,o2/np.max(o2))
plt.title(r'$T_{adiabatic}$ vs. Position',fontsize=25)
plt.xlabel(r'Position [m]', fontsize=15)
plt.ylabel('Normalized values of different quantities',fontsize=15)
plt.legend(['Temperature','$Y_{C_2H_6}$', '$Y_{O_2}$'],fontsize=15)
show()
