Slight inconsistency in 1D diffusion flux computation

[g3bk47]

For 1D flames, the diffusion coefficients are computed based on thermo-pyhsical properties at the midpoint between grid points, e.g. here: https://github.com/Cantera/cantera/blob/main/src/oneD/StFlow.cpp#L630

For the multi-component diffusion model, the density is evaluated based on these midpoint properties as well (https://github.com/Cantera/cantera/blob/main/src/oneD/StFlow.cpp#L614), which (to my understanding) is the correct approach.

For the mixture-averaged (and unity Lewis number) model, the density is evaluated at the grid point instead of the midpoint: https://github.com/Cantera/cantera/blob/main/src/oneD/StFlow.cpp#L675. density(j) refers to the grid point, not the midpoint: https://github.com/Cantera/cantera/blob/main/include/cantera/oneD/StFlow.h#L486. But since the density is part of the diffusion coefficient from the point of view of the Laplacian term, I think it should be evaluated at the midpoint as well. (This also affects IonFlow (https://github.com/Cantera/cantera/blob/baf28afdaf28e4998cdcdda9b727168f57e42e8a/src/oneD/IonFlow.cpp#L129).)

This change probably does not affect the results by much, but it might be more consistent. It looks like m_diff is not used outside of the 1D C++ code, so one easy fix would be to store the density together with the diffusion coefficient as m_diff[i] = rho*m_diff[i];, as currently done for the multi-component model.

Let me know what you think.

[g3bk47]

I did a test implementation of this in https://github.com/g3bk47/cantera/tree/midpointDensity.

For the standard C++ flame speed example (using methane with the GRI3.0 mechanism in samples/cxx/flamespeed/flamespeed.cpp), the difference in mixture-averaged flame speed is about 0.7%:

mixture-averaged flame speed = 0.316060775353 m/s (original)
mixture-averaged flame speed = 0.318441082012 m/s (with midpoint density)

For the standard python example for H2/O2 combustion with the h2o2.yaml mechanism (in samples/python/onedim/adiabatic_flame.py), the difference is about 0.8%:

mixture-averaged flame speed = 0.711305 m/s (original)
mixture-averaged flame speed = 0.705816 m/s (with midpoint density)

When the tolerances are chosen much tighter, i.e.

• ratio = 2.0, slope=0.001, curve=0.001
• steady tolerances: (1e-9, 1e-15)
• transient tolerances: (1e-9, 1e-15),

which results in flames with more than 6000 grid points, the difference in flame speed reduces to 0.005%.

mixture-averaged flame speed = 0.742162 m/s (original)
mixture-averaged flame speed = 0.742128 m/s (with midpoint density)

I now also included the mean molecular weight used in the flux expression at the midpoint as the same reasoning applies for this quantity.

[speth]

Thanks for looking into this. I think this is clearly an improvement in the discretization, even if the existing form is technically consistent (in that they converge to the original equation in the limit as $\Delta z \to 0$).

I was curious if this helped with a quirk that has been observed a few times, which is that there is some non-conservation of enthalpy, even when using the Le=1 transport model, although the degree of non-conservation does decrease as a finer grid is used. I tested this, and found that your changes definitely help in this case, at least within the flame, although the overall enthalpy difference across the flame is not affected by much. This test was based on the adiabatic_flame.py example.

[g3bk47]

This is interesting, I was not aware of this. I tried a few things to see if they affect the enthalpy balance. But they all had negligible effects:

• Interpolating the diffusion coefficients (and heat conductivity) evaluated on the grid points to the midpoint (i.e. D[midpoint] = 0.5*(D(i]+D[i-1])), instead of evaluating the thermo-properties at the midpoint
• formulating the convective and temporal term in a conservative way (cp*d(rho*T)/dt + cp*div(rho*u*T))
• Adding kinetic energy
• using the simplified Fickian flux expression for unity Lewis number as j=-rho*D*grad(Y) without correction velocity (which is mathematically equivalent to the current formulation)

All elemental mass fractions and rho*u seem to be conserved quite well, so it looks like the problem lies in the energy equation formulation? One idea would be to change to a (sensible) enthalpy equation, but this would make the fixed temperature point approach more complicated and would probably result in some performance penalty when converting from enthalpy back to temperature at each grid point.

For now, I will prepare a PR with the changes I have done so far.

[speth]

The other place I'd had a suspicion about this was in the computation of the "enthalpy flux" term, $\sum j_k dh_k/dz$, written for the ideal gas as $\sum jk c{p,k} dT/dz$. The update to support real gases was a relatively recent update (https://github.com/Cantera/cantera/pull/1079), but it preserved a feature that I thought might be part of the current issue. To evaluate this term, we compute the flux as 0.5*(m_flux(k,j-1) + m_flux(k,j), which depends on points $j-1$, $j$, and $j+1$ but then use an upwinded evaluation of $dh_k/dz$. Previously, this was the same, only with the temperature, and PR 1079 just preserved the existing stencil.

[g3bk47]

I tried replacing

for (size_t k = 0; k < m_nsp; k++) {
    double flxk = 0.5*(m_flux(k,j-1) + m_flux(k,j));
    sum += flxk * m_dhk_dz(k,j) / m_wt[k];
}

with

setGasAtMidpoint(x,j-1); // option 1
//setGas(x,j); // option 2
std::vector<double> D(m_nsp);
m_trans->getMixDiffCoeffs(D.data());
double rho = m_thermo->density();

for (size_t k = 0; k < m_nsp; k++) {
    // only works for this unity Lewis case
    sum += -rho*D[k]*dYdz(x,k,j) * m_dhk_dz(k,j) / m_wt[k];
 }

When running the adiabatic_flame.py example with UnityLewis and ratio=2, slope=0.02, curve=0.02, the enthalpy discrepancy only improves marginally:

Enthalpy in:  1305 J/kg/K
Enthalpy out: 1213 J/kg/K (using my PR and 0.5*(m_flux(k,j-1) + m_flux(k,j)))
Enthalpy out: 1229 J/kg/K (using my PR and -rho*D[k]*dYdz(x,k,j) with option 1)
Enthalpy out: 1216 J/kg/K (using my PR and -rho*D[k]*dYdz(x,k,j) with option 2)
Enthalpy out: 1240 J/kg/K (using my PR and -rho*D[k]*dYdz(x,k,j) with option 2, and also reverting back to cp_k dT/dz from dhk/dz)

[speth]

Resolved by https://github.com/Cantera/cantera/pull/1626