Solving the flamelet equations with FiPy

[marvinpom]

Hello,

I have some issues with applying FiPy to my problem and would like to ask for your assistance.

Physical problem description:

I would like to use FiPy to solve the steady-state counterflow diffusion flame problem at constant prescribed background pressure, i.e. the classical "flamelet" equations in mixture fraction coordinates (1D mesh) for combustion modeling. The coupled set of equations is composed of the species transport equations as well as the energy equation:

$\frac{\partial Y_{i}}{\partial t} = 0 = \frac{\chi}{2}\frac{\partial^{2} Y_{i}}{\partial Z^{2}} + \frac{\dot{\omega_{i}}}{\rho}$ where $i = 1,...,N_{\mathrm{species}}$
$\frac{\partial T}{\partial t} = 0 = \frac{\chi}{2}\frac{\partial^{2} T}{\partial Z^{2}} - \frac{1}{\rho c_{p}} \sum_{i}{h_{i} \dot{\omega}_{i}}$

$Y_{i}$: mass fraction of chemical species $i$ [kg kg^-1]
$T$: temperature [K]
$\chi$: scalar dissipation rate [s^-1]
$Z$: mixture fraction [-]
$\dot{\omega_{i}}$: chemical source term of chemical species $i$ [kg m^-3 s^-1]
$\rho$: mass-specific density [kg m^-3]
$c_{p}$: mass-specific isobaric heat capacity [J kg^-1 K^-1]
$h_{i}$: mass-specific (absolute) enthalpy of chemical species $i$ [J kg^-1]

In the end, I am interested in the chemical mixture composition in terms of mass fractions as a function of mixture fraction $Y_{i}(Z)$ as well as in the temperature profile $T(Z)$. The mixture fraction represents the local fuel content in the mixture and is restricted to the interval [0, 1]. The scalar dissipation rate, which acts as the diffusion coefficient of the problem under consideration ($D = \chi/2$), depends on the mesh (i.e. the mixture fraction) and results from the analytical profile

$\chi(Z) = \chi_{\mathrm{st}}\mathrm{exp}(2((\mathrm{erfc}^{-1}(2 Z_{\mathrm{st}}))^{2}-(\mathrm{erfc}^{-1}(2 Z))^{2}))$

where the index "st" represents the value at stoichiometric mixing conditions. The value of $\chi_{\mathrm{st}}$ is prescribed and kept constant. The value of $Z_{\mathrm{st}}$ is also constant for a fixed propellant combination and thus known a-priori.

The remaining variables are functions of pressure $p$ and the solution variables $Y_{i}, T$, i.e. $\dot{\omega_{i}}, \rho, c_{p}, h_{i} = f(p, Y_{i}, T)$. I can determine their values with the help of another class of my python program, as soon as I know (the pressure,) the mass fraction vector and the temperature.

I would like to apply Dirichlet boundary conditions for the propellant compositions and injection temperatures:

@ oxidizer inlet (corresponds to mixture fraction $Z = 0$ ~ no fuel in mixture)
if i == ox_index:
$Y_{i}(Z=0) = 1$
else:
$Y_{i}(Z=0) = 0$

$T(Z=0) = T_{\mathrm{ox}}$

@ fuel inlet (corresponds to mixture fraction $Z = 1$ ~ pure fuel in mixture)
if i == fu_index:
$Y_{i}(Z=1) = 1$
else:
$Y_{i}(Z=1) = 0$

$T(Z=1) = T_{\mathrm{fu}}$

Furthermore, resonable initial guesses for the profiles of the solution variables are available.

I would like to employ the Generalized Minimum Residual (GMRES) solver method from PETSc.

Code:

I attach the code snippets of my attempt to solve the Flamelet equations with FiPy as .txt file:

solver.txt

Issues:

1. With this code, I observe a strange solver behavior, which may be due to the fact that both the respective source terms in the species transport equations
   $S_{\mathrm{species}} = \frac{\dot{\omega_{i}}}{\rho}$
   as well as the source term in the energy equation
   $S_{\mathrm{energy}} = -\frac{1}{\rho c_{p}} \sum_{i}{h_{i} \dot{\omega}_{i}}$
   are evaluated only once at the beginning and are not updated again according to the current values of the solution variables in each iteration. How can I ensure that the source terms always correspond to the current values of the solution variables and are therefore consistent?

2. The documentation (https://www.ctcms.nist.gov/fipy/documentation/SOLVERS.html) states that no preconditioners are implemented for PETSc. However, in the generated log file it is stated that the incomplete LU-decomposition ("ilu") is applied. Is the documentation not up to date or is it just a "dummy" entry in the log-file? Since the system of equations is very stiff due to the chemical reaction terms, a suitable preconditioning is probably indispensable. Since Trilinos and PySparse are only available for python 2 and the documentation states that there are no preconditioners available for either PETSc or scipy.sparse solvers as well, I wonder what I could do. Is PyAMG in combination with a scipy solver the only alternative?

I would like to thank you for your help in advance.

Best regards

[guyer]

1. FiPy supports lazy evaluation and can operate on whole fields at once, so, rather than calculating a field and assigning it to the value of a CellVariable, you should write the expression you want and us it as a CellVariable. Here's my attempt to do so with compute_species_source_term():

   def compute_diffusion_coefficient(
           Z_st,
           SDR_st,
           mesh
   ):
       Z = mesh.faceCenters[0]
   
       return SDR_st * np.exp(2. * (np.square(erfcinv(2. * Z_st)) - np.square(erfcinv(2. * Z)))) / 2.
   
   def compute_species_source_term(
           phase,
           phi_vec,
           mesh,
           species_index_of_interest
   ):
       # update thermophysical state in terms of mass fractions and temperature
       phase.Y_vec = phi_vec[:len(phase.species_name_vec)]
       phase.T = phi_vec[-1]
   
       return phase.net_mass_production_rate_vec()[species_index_of_interest] / phase.rho()
   
   def compute_energy_source_term(
           phase,
           phi_vec,
           mesh
   ):
       """
       Description: Source term in energy conservation equation.
       """
       # update thermophysical state in terms of mass fractions and temperature
       phase.Y_vec = phi_vec[:len(phase.species_name_vec)]
       phase.T = phi_vec[-1]
   
       # compute auxiliary [J m^-3 s^-1] vector
       auxiliary_vec = phase.ha_vec().dot(phase.net_mass_production_rate_vec())
   
       # store respective value of energy source term in each cell center
       return -1. / (phase.rho() * phase.cp()) * auxiliary_vec

   Similar changes are likely needed in
   phase.net_mass_production_rate_vec().

2. The documentation does not say that no preconditioners are implemented for PETSc. It says

   FiPy does not implement any precoditioner objects for PETSc. Simply pass one of the PCType strings in the precon= argument when declaring the solver.

   So, you can, for example, call LinearGMRESSolver(..., precon="jacobi") or LinearGMRESSolver(..., precon="rowscalingviennacl") to get the corresponding preconditioner from https://petsc.org/release/manualpages/PC/PCType/. (No, I have no idea what the latter one is).

[marvinpom]

Thank you very much for the quick assistance.

Regarding 1.:
Thank you for pointing that out. I understand what you mean. I have adjusted the evaluation of the diffusion coefficient in the associated function accordingly. However, for the source terms, such an adjustment is not as straightforward.

To explain why this is the case, I would like to provide an executable program: flamelet_solver_FiPyWithCantera.zip. I have still calculated the source terms and assigned them to a CellVariable to keep it executable for the time being. To run the code, you only need a Cantera installation, since the abstract 'phase' object is now a Cantera object, which is used to get the thermochemical quantities of interest. However, the principle is exactly the same. To be able to (re-)compute the thermochemical quantities $\dot{\omega_{i}}, \rho, c_{p}$ and $h_{i}$ from Cantera, for instance, by

rho = phase.density_mass

and eventually determine the source terms $S_{\mathrm{species}}$ and $S_{\mathrm{energy}}$ as a function of the solution variables $Y_{i}$ and $T$, the state of the Cantera 'phase' object must first be updated by the current values of the solution variables temperature $T$ and mass fraction vector $Y_{i}$ (note: operating pressure $p$ is constant). This is accomplished with the following Cantera 'phase' object method:

phase.TPY = T, p, Y_vec

At this point, Cantera expects a scalar value for the temperature and pressure, and a 1D mass fraction array (containing the mass fractions of each species involved in mixture). It is therefore not possible to pass the entire temperature field and a 2D mass fraction array over all species and mixture fractions (i.e. over the entire mesh). It is only possible to evaluate one grid point, i.e. one mixture fraction at a time.

Still, I hope there is a workaround that allows updating the source terms in each solution step according to the current values of the solution variables. Do you have any idea how I could handle this problem? I would really like to continue working with FiPy.

Final remark:
As you can see in the plots, the steady-state solution determined with FiPy represents pure inert mixing. No combustion processes take place, which are expressed by the source terms. Theoretically, the flamelet solution at the specified very low stoichiometric scalar dissipation rate should not deviate significantly from the initial solution (--> represents the chemical equilibrium solution which would theoretically result at $\chi_{\mathrm{st}} = 0$).

Regarding 2.:
Thanks for the clarification!

[guyer]

Those sources definitely complicate things. That's an unfortunate design on Cantera's part, but they're hardly the only C programmers who write bad (by which I mean non-vectorizable) Python interfaces.

Here's how I would start to approach things, to at least get the sources to update every sweep (right now, you only get the initial value of the sources, which might be why you're not observing combustion):

diff --git a/flamelet_solver_FiPyWithCantera.py b/flamelet_solver_FiPyWithCantera.py
index 26673cf..669a6fb 100644
--- a/flamelet_solver_FiPyWithCantera.py
+++ b/flamelet_solver_FiPyWithCantera.py
@@ -80,15 +80,21 @@ def solve_and_plot(
     # declare empty list for the equations
     eq_vec = []
 
+    source_vec = []
+
     # 1. define species transport equations
     # loop over species
     for species_index, species_name in enumerate(phase.species_names):
         # append equation
-        eq_vec.append(0 == DiffusionTerm(coeff=compute_D(Z_st, SDR_st, mesh), var=phi_vec[species_index]) + compute_S_species(phase, p, phi_vec, mesh, species_index))
+        S_species = CellVariable(mesh=mesh, name=f"S_{species_name}")
+        source_vec.append(S_species)
+        eq_vec.append(0 == DiffusionTerm(coeff=compute_D(Z_st, SDR_st, mesh), var=phi_vec[species_index]) + S_species)
 
     # 2. define temperature equation
     # append equation
-    eq_vec.append(0 == DiffusionTerm(coeff=compute_D(Z_st, SDR_st, mesh), var=phi_vec[-1]) + compute_S_energy(phase, p, phi_vec, mesh))
+    S_energy = CellVariable(mesh=mesh, name="S_energy")
+    source_vec.append(S_energy)
+    eq_vec.append(0 == DiffusionTerm(coeff=compute_D(Z_st, SDR_st, mesh), var=phi_vec[-1]) + S_energy)
 
     # 3. couple equations
     # loop over all equations
@@ -132,6 +138,11 @@ def solve_and_plot(
             print("STOP: Maximum number of iterations exceeded.")
             break
 
+        for species_index, (species_name, S_species) in enumerate(zip(phase.species_names, source_vec[:-1])):
+            S_species.value = compute_S_species(phase, p, phi_vec, mesh, species_index)
+
+        S_energy.value = compute_S_energy(phase, p, phi_vec, mesh)
+
         # sweep until convergence criterion met
         res = coupled_eq.sweep(dt=time_step, solver=solver)
 
@@ -226,8 +237,7 @@ def compute_S_species(
     # set mixture fraction vector
     Z_vec = mesh.cellCenters[0]
 
-    # create CellVariable
-    S_species = CellVariable(mesh=mesh)
+    S_species = []
 
     # declare array
     Y_vec = numerix.zeros(len(phase.species_names))
@@ -253,7 +263,7 @@ def compute_S_species(
         omega_vec = phase.net_production_rates * phase.molecular_weights
 
         # store value of species source term at current grid point
-        S_species.put(Z_index, omega_vec[species_index_of_interest] / rho)
+        S_species.append(omega_vec[species_index_of_interest] / rho)
 
     return S_species
 
@@ -275,7 +285,7 @@ def compute_S_energy(
     Z_vec = mesh.cellCenters[0]
 
     # create CellVariable
-    S_energy = CellVariable(mesh=mesh)
+    S_energy = []
 
     # declare array
     Y_vec = numerix.zeros(len(phase.species_names))
@@ -307,7 +317,7 @@ def compute_S_energy(
         h_vec = phase.partial_molar_enthalpies / phase.molecular_weights
 
         # store value of energy source term at current grid point
-        S_energy.put(Z_index, -numerix.dot(h_vec, omega_vec) / (rho * cp))
+        S_energy.append(-numerix.dot(h_vec, omega_vec) / (rho * cp))
 
     return S_energy

Unfortunately, this blows up after three sweeps with

*******************************************************************************
CanteraError thrown by Phase::setTemperature:
temperature must be positive. T = -1524829600791612.5
*******************************************************************************

To stabilize things, you may need to employ relaxation (add a TransientTerm and evolve the solution toward steady state).

Once this is working, I can show you how to make custom CellVariable objects that will update whenever their dependencies change.

Note: Since your equations don't involve TransientTerm, you don't need hasOld=True or .updateOld(). If you do end up using relaxation, then you should only call .updateOld() for each time step, not for each sweep.

[marvinpom]

Once again, thank you very much for the detailed guidance!

I have implemented the workaround to calculate the source terms accordingly. However, I did not achieve convergent solver behavior even including the transient terms... Instead, I realized two previously neglected modifications and was able to achieve convergence w/o transient terms, at least for large stoichiometric scalar dissipation rates:

flamelet_solver_FiPyWithCantera.zip

Details:

1. Modification:
The consideration of a density correction factor in the calculation procedure of the diffusion coefficient of species transport and temperature equations. The scalar dissipation rate is now determined from:

$\chi = \chi_{\mathrm{st}} \frac{\Theta}{\Theta_{\mathrm{st}}} \mathrm{exp} \left[2 \left((\mathrm{erfc}^{-1} (2 Z_{\mathrm{st}}))^{2} - (\mathrm{erfc}^{-1} (2 Z))^{2}\right)\right]$

where $\Theta = \frac{3}{4} \frac{\left(\sqrt{\frac{\rho_{\mathrm{ox}}}{\rho}} + 1\right)^{2}}{2 \sqrt{\frac{\rho_{\mathrm{ox}}}{\rho}} + 1}$

It now depends on the density, which in turn depends on the solution variables. As a consequence the diffusion coefficient must now also be updated continuously during the solution process

2. Modification:
I have included the convective term of the temperature equation which is usually neglected:

$0=\frac{\chi}{2 c_{p}} \left(\frac{\partial c_{p}}{\partial Z} + \sum_{i} c_{p,i} \frac{\partial Y_{i}}{\partial Z} \right) \frac{\partial T}{\partial Z} + \frac{\chi}{2} \frac{\partial^{2} T}{\partial Z^{2}} - \frac{1}{\rho c_{p}} \sum_{i} h_{i} \dot{\omega}_{i}$

However, the high stoichiometric scalar dissipation rate (~ high flame strechting) probably causes the flame to extinguish. Finally, I obtain again only the inert mixing solution w/o combustion for the given input parameters. As soon as I reduce its value further or choose unequal injection temperatures, I run into numerical instabilities until the solution diverges or negative temperature values result.

Questions:

1. With the workaround you suggest, can I really rule out the possibility that the stability problems and the fact that I can only reach the inert mixing solution are caused by inconsistent (i.e. not updated in intermediate solution steps) coefficients and source terms? Maybe updating them only after each sweep is not sufficient... Would the custom CellVariable you mentioned be an improvement in this context?

2. I tried to compute the coefficient of the convective temperature term in the same manner as for the diffusive term and the source terms. However, to make the gradient of the isobaric heat capacity accessible by the internal FiPy routines, I have defined a CellVariable within the corresponding function. In total, I need two loops over the computational grid. Is there a more elegant way to (re-)compute the value of the convection coefficient? Is the .grad method even the one I should use for $c_{p}$ and $Y_{i}$? How could I theoretically represent the second derivative of the species mass fraction vector (i.e. $\partial^{2} Y_{i} / \partial Z^{2}$)?

3. Do you have any other ideas on how I could stabilize the solution process of the coupled equations in FiPy? What about solving the uncoupled equations or introducing underrelaxation, for example? How do I specify it correctly? Specifying a scalar value (e.g. coupled_eq.sweep(dt=time_step, solver=solver, underRelaxation=0.1)) results in an error message.

I am very looking forward to your answers. Thank you in advance!

[guyer]

1. With the workaround you suggest, can I really rule out the possibility that the stability problems and the fact that I can only reach the inert mixing solution are caused by inconsistent (i.e. not updated in intermediate solution steps) coefficients and source terms? Maybe updating them only after each sweep is not sufficient... Would the custom CellVariable you mentioned be an improvement in this context?

You already update everything every sweep. The only way to make things update more often is to decouple and update after the sweep of each successive equation. There are some cases where that can help, e.g., semiconductor drift-diffusion equations, SIMPLE method for fluid flow, ... Usually, though (and definitely in those cases) there's a bunch of other algorithmic stuff or fancy damping going on for each partial solve. Don't get into that mess if you don't have to.

Generally speaking, making things more complicated seldom helps with convergence or stability.

For troubleshooting, I recommend some combination of:

• Remove the reaction sources (or use much simpler proxies).
• Use simpler expressions for your coefficients (D and C).
• In short, break the connection with Cantera until you have the basic equations working in FiPy.
• Remove the convection term (as I'll post below, do this anyway; this term is inappropriate in FiPy).
• Decouple your equations and solve them sequentially. They may require different damping or different solvers.
• Use relaxation (add a TransientTerm). Steady-state solutions often don't work without help.
• (As you have) use under relaxation.

Once you start getting converged solutions, you can start gradually reintroducing complexity.

[guyer]

2. I tried to compute the coefficient of the convective temperature term in the same manner as for the diffusive term and the source terms. However, to make the gradient of the isobaric heat capacity accessible by the internal FiPy routines, I have defined a CellVariable within the corresponding function. In total, I need two loops over the computational grid. Is there a more elegant way to (re-)compute the value of the convection coefficient? Is the .grad method even the one I should use for cp and Yi? How could I theoretically represent the second derivative of the species mass fraction vector (i.e. ∂2Yi/∂Z2)?

$(\partial c_p/\partial Z)(\partial T/\partial Z)$ is not a FiPy ConvectionTerm; $\partial(T \vec {u})/\partial Z$ is a ConvectionTerm.
For that matter, $(\chi/2)(\partial^2 T/\partial Z^2)$ isn't a FiPy DiffusionTerm; $\partial [(\chi/2) (\partial T/\partial Z)] / \partial Z$ is.

In case you are unaware, what you call the "convective term of the temperature equation" is a consequence of writing:

$$ \begin{aligned} \frac{\partial T}{\partial t} &= \frac{\partial}{\partial Z}\left(\frac{\chi}{2}\frac{\partial T}{\partial Z}\right) + \cdots \\ &= \frac{\partial}{\partial Z}\left(\frac{\Xi}{\rho c_p}\frac{\partial T}{\partial Z}\right) + \cdots \end{aligned} $$

Following the spatial heat equation, here $\Xi$ is some combination of the conventional thermal conductivity $k$ and a conversion like $(z/Z)^2$ between spatial and reaction coordinates. Running the chain rule:

$$ \begin{aligned} \frac{\partial T}{\partial t} &= \frac{\Xi}{\rho c_p}\frac{\partial^2 T}{\partial Z^2} +\frac{\partial}{\partial Z}\left(\frac{\Xi}{\rho c_p}\right) \frac{\partial T}{\partial Z} + \cdots \\ &= \frac{\Xi}{\rho c_p}\frac{\partial^2 T}{\partial Z^2} + \frac{1}{\rho c_p}\frac{\partial \Xi}{\partial Z} \frac{\partial T}{\partial Z} - \frac{\Xi}{\rho^2 c_p}\frac{\partial \rho}{\partial Z} \frac{\partial T}{\partial Z} - \frac{\Xi}{\rho c_p^2}\frac{\partial c_p}{\partial Z} \frac{\partial T}{\partial Z} + \cdots \\ &= \frac{\chi}{2}\frac{\partial^2 T}{\partial Z^2} + \frac{\chi}{2\Xi}\frac{\partial \Xi}{\partial Z} \frac{\partial T}{\partial Z} - \frac{\chi}{2\rho}\frac{\partial\rho}{\partial Z}\frac{\partial T}{\partial Z} - \frac{\chi}{2 c_p}\frac{\partial c_p}{\partial Z}\frac{\partial T}{\partial Z} + \cdots \end{aligned} $$

Since, evidently, $c_p = c_p(Z, Y_i)$, the fourth term on the right-hand side expands to

$$ \begin{aligned} \frac{\partial T}{\partial t} &= \frac{\chi}{2}\frac{\partial^2 T}{\partial Z^2} + \frac{\chi}{2\Xi}\frac{\partial \Xi}{\partial Z} \frac{\partial T}{\partial Z} - \frac{\chi}{2\rho}\frac{\partial\rho}{\partial Z}\frac{\partial T}{\partial Z} - \frac{\chi}{2 c_p}\left( \left.\frac{\partial c_p}{\partial Z}\right\vert_{Y_i} + \sum_i \left.\frac{\partial c_p}{\partial Y_i}\right\vert_{Z,Y_{i\ne j}} \left.\frac{\partial Y_i}{\partial Z}\right\vert_{Y_{i\ne j}} \right)\frac{\partial T}{\partial Z} + \cdots \end{aligned} $$

The expressions $\partial \Xi/\partial Z$ and $\partial\rho/\partial Z$ will have similar expansions.

Having done all that, however, since $\chi$ depends on $Z$, $\rho$, $c_p$, $Y_i$, $\ldots$ via your calls to Cantera, it is inappropriate to include these terms when writing your PDE in FiPy; their variation in space (or in the reaction coordinate) is already accounted for by the DiffusionTerm. You should solve

$$ \frac{\partial T}{\partial t} = 0 = \frac{\partial}{\partial Z}\left(\frac{\chi}{2}\frac{\partial T}{\partial Z}\right) + \cdots $$

and write it as

0 == DiffusionTerm(coeff=D, var=phi_vec[-1]) + source_vec[-1]

like you did originally.

If, despite all this, you're committed to including such terms, then, yes, $\partial c_p/\partial Z$ would be Cp.grad, $\partial Y_i/\partial Z$ would be phi_vec[species_index].grad, and $(\partial c_p/\partial Z) (\partial T/\partial Z)$ would be Cp.grad.dot(T.grad). Such an expression would be first-order in the reaction coordinate and would be evaluated explicitly, which are two more reasons not to do this. $\partial^2 Y_i/\partial Z^2$ would be phi_vec[species_index].faceGrad.divergence, which is explicit, but at least it's second-order.

If you really, really want to go down this route, I can show you how to remove some of the redundant looping in compute_C(), but I don't think you should be calling that function at all.

---

note: Many treatments of the heat equation (like on Wikipedia) put $\rho c_p$ on the left-hand side, with the time derivative, and outside the divergence:

$$ \rho c_p \frac{\partial T}{\partial t} = \nabla\cdot\left(k \nabla T\right) + \cdots $$

If this is the case, your terms arising from $\partial c_p(Z,Y_i)/\partial Z$, etc. would not arise at all. Further, a FiPy TransientTerm represents $\partial (\rho c_p T)/ \partial t$, not $\rho c_p \partial T/\partial t$, so further care would need to be taken to account for $\partial \rho/\partial t$ and $\partial c_p/\partial t$. I don't know what the right formulation is; I haven't researched it deeply, but most of the treatments I've seen of the heat equation are pretty sloppy about whether coefficients vary in time and/or space. If I were doing this for my own research, I'd want to go back to a statement of energy conservation (or entropy production) and be very careful and explicit about dependent variables.

[guyer]

3. Do you have any other ideas on how I could stabilize the solution process of the coupled equations in FiPy? What about solving the uncoupled equations or introducing underrelaxation, for example? How do I specify it correctly? Specifying a scalar value (e.g. coupled_eq.sweep(dt=time_step, solver=solver, underRelaxation=0.1)) results in an error message.

There's evidently a problem with underRelaxation when using the PETSc, Trilinos, and PySparse solvers on coupled equations (#971). It works with the scipy solvers. You are calling it correctly.