Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Fixing some errors calculating the Diffusion coefficients and thermal conductivity for kinetics theory. #267

Open
wants to merge 4 commits into
base: master
Choose a base branch
from
Open

Fixing some errors calculating the Diffusion coefficients and thermal conductivity for kinetics theory. #267

Show file tree
Hide file tree
Changes from 3 commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
49eff2e
Changing how we calculate the mixture averaged diffusion coefficients…
klbudzin Aug 20, 2018
3ad5a4a
Changing the way we calculate the mixture thermal conductivity from t…
klbudzin Aug 24, 2018
13f6179
Fixing the unit test that was tripped when altering the calculation o…
klbudzin Aug 24, 2018
aedef6f
Instead of clipping the mole fractions, implemented a method that wil…
klbudzin Aug 28, 2018
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
33 changes: 24 additions & 9 deletions src/transport/include/antioch/mixture_averaged_transport_evaluator.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -344,6 +344,8 @@ namespace Antioch

mu_mix = zero_clone(T);
k_mix = zero_clone(T);
StateType k_sum1 = zero_clone(T);
StateType k_sum2 = zero_clone(T);
D_vec = zero_clone(mass_fractions);

VectorStateType mu = zero_clone(mass_fractions);
Expand Down Expand Up @@ -398,11 +400,12 @@ namespace Antioch
mu[s] );
}

k_mix += k[s]*chi[s]/phi_s;
k_sum1 += k[s]*chi[s];
k_sum2 += chi[s]/k[s];
mu_mix += mu[s]*chi[s]/phi_s;

}

k_mix = .5*(k_sum1+1/k_sum2);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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).

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok, works for me then. Still curious why it would be an improvement for conductivity but not viscosity, though.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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...

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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. :(



if( DiffusionTraits<Diff>::is_species_diffusion )
Expand Down Expand Up @@ -579,12 +582,23 @@ namespace Antioch
}
case(MASS_FLUX_MASS_FRACTION):
{
VectorStateType molar_fractions = zero_clone(mass_fractions);

mixture.X(mixture.M(mass_fractions),mass_fractions,molar_fractions);


VectorStateType molar_fractions = zero_clone(mass_fractions);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What happened to the indentation here?

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@klbudzin add the following to your .emacs file in your home directory:

(setq-default indent-tabs-mode nil)
(add-hook 'before-save-hook 'delete-trailing-whitespace)

or if you're using VIM, ask @roystgnr for the appropriate settings for .vimrc. :)

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

mixture.X(mixture.M(mass_fractions),mass_fractions,molar_fractions);

//Clipping unrealistic values, and calculating MW_Mixture through the new mole fractions
typename value_type<VectorStateType>::type MW_Mixture = zero_clone(mass_fractions[0]);
for(unsigned int s=0; s < D_vec.size(); s++)
{
if(molar_fractions[s] > 1)
molar_fractions[s] = 1;
if(molar_fractions[s] < 1e-16)
molar_fractions[s] = 1e-16;
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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).

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Oh, I love the idea of a clip_values = false default argument.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

MW_Mixture += molar_fractions[s]*mixture.M(s);
}

typename value_type<VectorStateType>::type one = constant_clone(mass_fractions[0],1);

// term1 term2
// 1/D_s = (sum_{j\ne s} X_j/D_{s,j}) + X_s/(1-Y_s)\sum_{j\ne s} Y_j/D_{s,j}
for(unsigned int s = 0; s < D_vec.size(); s++)
Expand All @@ -599,15 +613,16 @@ namespace Antioch

term1 += molar_fractions[j]/D_mat[s][j];

term2 += mass_fractions[j]/D_mat[s][j];
term2 += molar_fractions[j]*mixture.M(j)/D_mat[s][j];
}

term2 *= molar_fractions[s]/(one - mass_fractions[s]);
term2 *= molar_fractions[s]/(MW_Mixture - mixture.M(s)*molar_fractions[s]);

D_vec[s] = one/(term1+term2);
}
break;
}

default:
{
antioch_error_msg("ERROR: Invalid DiffusivityType in MixtureAveragedTransportEvaluator::diffusion_mixing_rule");
Expand Down
2 changes: 1 addition & 1 deletion test/wilke_transport_unit.C
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -276,7 +276,7 @@ int tester()
*/

const Scalar mu_kt_long_double = 4.49877527305932602332e-05L;
const Scalar k_kt_long_double = 8.22050332419571328635e-02L;
const Scalar k_kt_long_double = 8.27884568936580333168e-02L;

std::vector<Scalar> D_kt_long_double(5,0);
D_kt_long_double[0] = 1.95418749838889089562e-04L;
Expand Down