Fixing some errors calculating the Diffusion coefficients and thermal conductivity for kinetics theory. #267
Conversation
… by using mole fractions instead of mass fractions, like cantera, which keeps the mixture molecular weight more bearable. As well as clipping mole fractions that were negative or over 1.
…he species conductivities. The formula used here is taken from "Chemically Reacting Flow" by Kee, Coltrin, and Glarborg, page 519.
…f the thermal conductivity
| mixture.X(mixture.M(mass_fractions),mass_fractions,molar_fractions); | ||
|
|
||
|
|
||
| VectorStateType molar_fractions = zero_clone(mass_fractions); |
There was a problem hiding this comment.
What happened to the indentation here?
There was a problem hiding this comment.
I guess I didn't check everything all the way through. My bad, I will ad that to my .emacs file, so hopefully that doesn't happen again.
There was a problem hiding this comment.
Secondly I converted how we calculate the diffusing mixing rule to use mole fractions instead of mass fractions. For some reason, using mass fractions gives us a different mixture molecular weight than the mole fractions. From my testing, the diffusion coefficients are very sensitive to the mixture molecular weight causing our values to be highly different from the expected values for near pure species mixture concentrations.
The mole fraction route is the way to go: I asked Sylvain to implement this but he never did(not sure why). For further stability, Ern & G. recommend the following treatment:
For transport computations, define X_i^tr = X_i + \epsilon(number of species)^{-1} \sum_{species} X_j - X_i), where epsilon is something like machine(double) precision. This helps with stability in the pure species limit. It also does not change the sum of the mole fractions(up to rounding, and produces mole fractions between 0 and 1. You then compute an adjusted mass fraction Y_i^tr and transport coefficients using X_i^tr and Y_i^tr.
There was a problem hiding this comment.
I believe @SylvainPlessis did this in a PR that never got merged and is now stale (#132) :(. @klbudzin can you have a look at what is in #132 and see if it corresponds to what you have here. If not, we should get that integrated.
| if(molar_fractions[s] > 1) | ||
| molar_fractions[s] = 1; | ||
| if(molar_fractions[s] < 1e-16) | ||
| molar_fractions[s] = 1e-16; |
There was a problem hiding this comment.
Should we shout a warning here if we hit these cases? This seems at first glance like the best thing to do, but imagine e.g. someone trying to take a derivative by central finite differencing around 0; adding clipping under the hood could turn the relative error from O(epsilon^2) to O(.5).
There was a problem hiding this comment.
I wonder if we should do this at all and leave it up to the user to truncate themselves before passing in? Perhaps add a default argument to the function?
There was a problem hiding this comment.
Oh, I love the idea of a clip_values = false default argument.
There was a problem hiding this comment.
Yeah, most of this was just a quick fix to check that this was what was causing the discrepancies seen between antioch and cantera. The default argument would definitely be better.
| mu_mix += mu[s]*chi[s]/phi_s; | ||
|
|
||
| } | ||
|
|
||
| k_mix = .5*(k_sum1+1/k_sum2); |
There was a problem hiding this comment.
We need this better commented. Of course, we needed the previous version much better commented...
Why is this new formula appropriate for k but not mu?
There was a problem hiding this comment.
He does need to document it more, but I'll chime in. It seems this is the more "accepted" mixing rule for k. Kee et al use it in their book and it seems this is also what others use, including Cantera. I'm doubtful one is better than the other (they both look like hacks to me), but the is more consistent with the community (and has been helpful for code-to-code comparisons with Cantera that @klbudzin has been doing).
There was a problem hiding this comment.
Ok, works for me then. Still curious why it would be an improvement for conductivity but not viscosity, though.
There was a problem hiding this comment.
I'm not sure really, the reacting flow book I obtained it from, basically just stated this formula for the conductivity ontop of the normal viscosity one. I'll look into it more to see if I can obtain a better reasoning.
There was a problem hiding this comment.
There are a class of empirical averaging formulas(the Kee formula is a version of one of these) that are less accurate(but much cheaper) than the Wilke formula. There is a long discussion of this in
A. Ern and V. Giovangigli, Multicomponent Transport Algorithms, Lecture Notes in Physics, New Series Monographs, m 24, (Springer-Verlag, Heidelberg, 1994).
I don't know why Sylvain implemented the Wilke rule for thermal conductivity(I told him to use the Kee formula ala the Chemkin manual for consistency). Both Chemkin and Cantera use the Kee empirical rule.
There was a problem hiding this comment.
Hello all,
Well I won't be able to remember the whole situation, but it seems to me that there was lots of considerations in this part. Isn't it also the part where there was an implementation given by a french fortran code? I remember using their tricks to deal with Titan's high atmosphere (like huge amount of N2, little of CH4 and traces of what we want to follow), and theirs was the only one that would not give the numerical equivalent of kicking the ground in fury, then sit down while breaking in tears... But I cannot certify it was this exact part of the code.
Might worth a try to look for that code, they might have left a better documentation than me...
Now there is also the argument that I totally prefer rules (like Wilke) than empirical, but that might be because DFT looks sometimes like a bad magic trick to me (Chris can very probably explain what that means). In that case, my bad, I'm sorry. There is also the possibility that it was just before I left and I did not do it because I did not find the time before leaving...
There was a problem hiding this comment.
It's entirely possible @SylvainPlessis did the Kee formula originally and I changed it because I didn't know any better. (I don't have time to troll through the history). Anyway, I don't want @SylvainPlessis to feel like he's getting dumped on because that's not my intention.
Besides I haven't gotten to speak to either @SylvainPlessis or @variscarey in forever! :(
The French code is EGLib (Ern, Giovangigli) is still available. It's in Fortran. :(
|
Is this going to trip a bunch of GRINS regression tests too? |
At least one. I told him to do this first then we'll update GRINS down stream. |
…l calculate perturbed mass fractions that can be used to calculate the transport, this will make it easier to implement a switch for perturbing the mass fractions. Also added a little more documentation, but can still use more here.
|
I've updated a few of the concerns here. I added a separate method to calculate the perturbed mass fractions using the EGlib method, and call it in the mu_k_and_d method before calculating any of the transport properties. I have still yet to add a default argument for the clipping, but I could easily add one to the mu_k_and_d method now. I also added some more documentation in. |
When comparing Transport values with cantera I noticed a few problems that arose in the calculation of the mixture averaged Diffusion coefficients. First, we are not clipping non-realistic values for mass fractions, so if the user doesn't do this themselves, were throwing out non-physical values. Secondly I converted how we calculate the diffusing mixing rule to use mole fractions instead of mass fractions. For some reason, using mass fractions gives us a different mixture molecular weight than the mole fractions. From my testing, the diffusion coefficients are very sensitive to the mixture molecular weight causing our values to be highly different from the expected values for near pure species mixture concentrations.
I also changed how we calculate the thermal conductivity to match how it is calculated in ("Chemically Reacting Flow" by Kee, Coltrin, and Glarborg, page 519). These values now seem to match to canteras calculated values to at least the third digit. Finally there is also a changed to the unit test that was tripped with the thermal conductivity change.