[Transport] Add unity Lewis number transport model

[awehrfritz]

Please fill in the issue number this pull request is fixing:

Fixes #507

Changes proposed in this pull request:

• Add unity Lewis number transport model

Supersedes pull requests #427 and #508 . The current implementation requires only minimal code changes and has been tested. The extension of the unity Lewis number transport model to a per-species Lewis number transport model is rather straight forward, but should be done in an own class.

[awehrfritz]

I am not sure why the OSx test is failing, but it appears unrelated to the current pull request, since the Linux test passed and the code changes should be agnostic to the operating system.

@speth you have reviewed and commented in detail the earlier pull request (#427). I believe I have considered most (all) your comments in this one and would appreciate if you could review this.

[bryanwweber]

The error in the macOS build is because of Homebrew ☹️ I'll look into fixing it

Can you add some tests for this new functionality? Tests for the Python interface go in the interfaces/cython/cantera/test directory, and tests for C++ go in the tests directory.

[awehrfritz]

@bryanwweber: I added a test with some basic sanity checks. The test passes and I think this is good to go on that end.

[speth]

I think this implementation of mixDiffCoeffs does not preserve one of the key properties of a Le=1 model, namely that Dh/Dt=0 when ∇h=0 (where h is specific enthalpy and D/Dt is the total derivative).

If we start with Eq. 3.203 from Chemically Reacting Flow (Kee, Coltrin, Glarborg; 1st ed.), and neglect the pressure and viscous dissipation terms, we have:
ρ Dh/dt = ∇·(k∇T) - ∑∇·(hk jk)
where hk and jk are the specific enthalpy and diffusive mass flux for species k, respectively, and k is the thermal conductivity.
In order to satisfy Dh/Dt = 0, then we must have
k ∇T = ∑ hk jk

If we define jk = - ρ Dkm ∇ Yk, (Eq. 12.177) which corresponds to the diffusion coefficients returned by the mixDiffCoeffsMass function, then we have
k ∇T = ∑ hk (-ρ Dkm ∇Yk)
If we assume that there is a single value for Dkm=Dm, we can pull it out of the summation and write
k ∇T = -ρ Dm ∑ hk∇Yk
Now, considering the condition ∇h=0, we can write
∇h=∇(∑hk Yk) = ∑ (hk ∇Yk + ∇hk Yk) = ∑ (hk ∇Yk + cpk Yk ∇T)=0
Substituting this result into the previous expression, we obtain
k ∇T = ρ Dm ∑ cpk Yk ∇T = ρ Dm ∇T ∑ cpk Yk = ρ Dm ∇T cp
And we find that Dm = k / (ρ cp), the expected result corresponding to Le=1. So implementing mixDiffCoeffsMass in this way would be fine.

On the other hand, if we define jk = - ρ Yk/Xk Dkm' ∇Xk (from Eq. 12.179), the definition which corresponds to mixDiffCoeffs, things turn out differently. Substituting this in to the enthalpy equation yields
k ∇T = ∑ hk (-ρ Dkm' Wk/Wmx ∇Xk)
where Wk is the molecular weight of species k and Wmx is the mixture molecular weight (using the identity Wk/Wmx = Yk/Xk)
If we again assume that Dkm' = Dm', some constant value, we can write
k ∇T = - ρ Dm' / Wmx ∑ hk Wk ∇Xk
The total enthalpy h can be written in terms of Xk as
h = ∑(hk Wk Xk / Wmx)
Then
0 = ∇h = ∑(∇hk Wk Xk / Wmx + hk Wk ∇Xk / Wmx - hk Wk Xk ∇Wmx / Wmx^2)
= ∑(cpk ∇T Wk Xk / Wmx + hk Wk ∇Xk / Wmx - hk Wk Xk ∇Wmx / Wmx^2)
Solving for the term in the enthalpy equation we wish to replace,
∑ hk Wk ∇Xk = ∑( - cpk Wk Xk ∇T + hk Wk Xk ∇Wmx / Wmx)
In this case, the appearance of a term proportional to ∇Wmx means that substituting this in does not let us cancel out ∇T and find a value for Dm' which is independent of these gradients.

I believe the same issue will occur with mixDiffCoeffsMole, where the mole fraction gradients again make an appearance.

[wandadars]

@speth From your analysis, do you think that a consistent expression for the diffusion coefficient under the Le=1 assumption can be derived for each of the 3 cases?

[speth]

What I think I have shown is that there is no definition of a single value (for all species) which could be returned by mixDiffCoeffs or mixDiffCoeffsMole that would actually give you the expected unity-Lewis number behavior. It doesn't necessarily show that there isn't a set of values (different for each species, presumably a function of molecular weight) that would work. Figuring out what those formulas might be (and demonstrating that they are correct) is left as an exercise for someone who actually wants to use such a model. 😂

[wandadars]

@speth Under what conditions would the different calls to MixDiffCoeffsMole and MixDiffCoeffss be made? Would it make sense to restrict this transport option to code which calls only the mixDiffCoeffsMass function? Something like having a 'feature not implemented' error upon calling the other functions?

[speth]

Well, you could do that, but the 1D flame code in Cantera uses the MixDiffCoeffs method, so I don't think it would be very useful.

[awehrfritz]

Apologies for the late reply on this issue, I could not find my original notes on the derivation of the unity Lewis number model anymore and then got caught up in other tasks.

@speth: In your derivation above, you did not consider the correction velocity which is required to ensure \sum j_k = 0. If you consider the correction velocity as defined in Kee et al. (2003), Eq. 12.183, the extra terms (containing the mean molecular weight) will drop out. Because the correction velocity is also implemented in Cantera I thus believe that my implementation of the unity Lewis number model is correct, and yields the properties you outlined above.

I have attached a document in which I derive the respective expressions:
UnityLewisNumberModelDerivation.pdf
I have not done the derivation for the equations you used above (i.e. starting from Eq. 12.177 in Kee et al.) and thus I prefer to not implement the respective functions (also considering that mixDiffCoeffsMass or mixDiffCoeffsMole are not used anywhere in Cantera).

Regarding the implementation, as you noted above, Cantera's 1D flame solver uses mixDiffCoeffs and thus also computes the diffusive mass flux from mole fraction gradients (including the correction term to ensure \sum j_k = 0). This is not the most efficient way for the unity Lewis number assumption, as one could simplify the equations if one was to use mass fraction gradients (outlined in the attached document). This would however require at least one additional if-then-else statement in the code. Either way, I would only implement mixDiffCoeffs for reasons mentioned above.

[speth]

@awehrfritz Thanks, that derivation is very helpful. The fact that applying that form of the correction velocity fixes the issue is not very intuitive (at least to me), but it does indeed work out as you have shown.

In the case of mixDiffCoeffsMass, things end up being pretty simple. The correction mass flux is ∑ ρ D_k ∇ Y_k, which just ends up being zero in the case where all of the D_k are the same. So, I think it would be reasonable to implement this method as well.

[speth]

I think this is just about good to go, pending some minor style cleanup.

[speth]

Please delete the commented out writelog statement.

[speth]

Please use double instead of doublereal (as noted in CONTRIBUTING.md)

[speth]

Dkm2 is assigned but never used

[speth]

Testing that the viscosity and binary diffusion coefficients are the same as returned by the Mix model is probably unnecessary.

[speth]

This should still be D^\prime_{km} if we're being consistent with Kee's notation, right?

[speth]

No need to define the constructor if it doesn't do anything.

[speth]

No need to explicitly reference this-> (here or anywhere else).

I think it would actually be more clear to just write this as

double Dm = thermalConductivity() / (m_thermo->density() * m_thermo->cp_mass());

[speth]

The @ingroup tranprops directive can be on just another line starting with //! rather than also starting a /*-style block.

[codecov]

Codecov Report

Merging #510 into master will increase coverage by <.01%.
The diff coverage is 90%.

@@            Coverage Diff             @@
##           master     #510      +/-   ##
==========================================
+ Coverage   64.74%   64.75%   +<.01%     
==========================================
  Files         383      384       +1     
  Lines       40689    40699      +10     
==========================================
+ Hits        26345    26354       +9     
- Misses      14344    14345       +1

Impacted Files	Coverage Δ
src/transport/TransportFactory.cpp	45.23% <100%> (+0.16%)	⬆️
include/cantera/transport/UnityLewisTransport.h	88.88% <88.88%> (ø)

---

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 3c978cd...ad31fd7. Read the comment docs.