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

Changes from 1 commit
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
26 changes: 19 additions & 7 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 @@ -579,12 +579,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 +610,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