libantioch / antioch

C++ Chemical Kinetics, Thermodynaimics, and Transport Library
https://libantioch.github.io/
Other
23 stars 17 forks source link

Correct density passed to KineticsTheoryThermalConductivity? #146

Open pbauman opened 9 years ago

pbauman commented 9 years ago

@SylvainPlessis

Here

You are passing c*M_s, i.e. molar density times species molar mass. What is that? Cross checking the documentation with the thermal conductivity implementation, it looks like you want species mass density, i.e. \rho_s. That is not c*M_s; that is, for example, rho* Y_s (or whatever symbol we're using for mass fraction). I actually put rho*Y_s in and the thermal conductivity gold value failed in wilke_evaluator_unit, but using c*M_s, as you have it now, the test passes.

So what should it be? I'm going to press forward with rho*Y_s, but I wanted to get confirmation from you here.

Once I have confirmation, I'll put in a PR with the correction and rebase my transport-refactor branch.

SylvainPlessis commented 9 years ago

Yes there is some thing here about this density, and of course, documenting it is still on the TODO list... Sorry about that. So the density thinggy is that is you use the species density (which is the documentation version), we should pass c_s*M_s and the comparison with EG lib fail. Now if you says you need the total density as if it was a pure mixture of species s, thus passing c*M_s, you find the EG lib values. To make things even more easier, the EG lib values are polynomials fits, so we don't have access to the theoretical implemented equations...

It requires investigations on the theoretical side, I didn't have the time to do it yet...

variscarey commented 9 years ago

The discrepancy arises from the different pressures used to compute D_{ss} and \rhos. The binomial diffusion coeff. D{jk} has a pressure dependence of 1/P. The conductivity $\lambdas$ is not pressure-dependent(nor does it depend on mixture ratios). TRANSPORT computes D{jk} at a reference pressure and then fits the non-temp. part with a polynomial. It looks like(at least when we just fit the collision integrals) that we compute the whole thing, using ideal gas for the mixture to compute P and modify the Hirschfielder formula. This gives rise to a mixture molar density c term in the denominator for D_{ij}. If we pass in \rho_s to the conductivity formula we are using the "wrong pressure", namely the partial pressure for species s(P_s). So why does "c * M_s" work for \rho_s. There is a \mu/M_s sitting outside the conductivity formula. I might be looking at the wrong source, but it looks like we expand the rot and vib terms, so the \mu's cancel and we are left with \rho/Ms * D{ss}. If we use c * M_s as \rho, the corresponding mixture and pressure contributions cancel(the 1/c in D_ss). Don't ask where the M_s^{-1} goes in the trans contribution, I haven't tried to track it down. This is my best guess at the problem and why our "hack" solution works.

I suspect that we may want to switch to a reference pressure setup(ala TRANSPORT)(at least for computing \lambda_s). Then the corresponding term would look like (after canceling M_s and P) (R_0 T)^(-1) D_ss(the pressure ind part) * other stuff.

Varis

On Sun, May 10, 2015 at 10:32 AM, Sylvain Plessis notifications@github.com wrote:

Yes there is some thing here about this density, and of course, documenting it is still on the TODO list... Sorry about that. So the density thinggy is that is you use the species density (which is the documentation version), we should pass c_s_M_s''' and the comparison with EG lib fail. Now if you says you need the total density as if it was a pure mixture of speciess''', thus passing ```c_M_s''', you find the EG lib values. To make things even more easier, the EG lib values are polynomials fits, so we don't have access to the theoretical implemented equations...

It requires investigations on the theoretical side, I didn't have the time to do it yet...

— Reply to this email directly or view it on GitHub https://github.com/libantioch/antioch/issues/146#issuecomment-100655973.

variscarey commented 9 years ago

An update/clarification:

The discrepancy is indeed caused specifically by the use of the pressure P in the bimolecular diff coeff.
Right now the pure species thermal conductivity as written depends on rho_s and T (through diff, visc, and the formula). If we want the dependence on c_s or rho_s to drop out we need the relevant terms to cancel when we compute rho_s * D_ss. What do I have to multiply to compensate for the fact that I used P rather than P_s to compute D_ss? P/P_s. But rho_s_P/P_s= (P_s M_s)/(R_0 T) * (c_tot * R_0 * T) /P_s. This all cancels to give c_tot_M_s, which is exactly what Sylvain wrote in the code. So you shouldn't use rho Y_s but this silly expression.

This is of course bad practice. We should compute pressure ind. bi molecular diffusion coeffs, and then use (for ps D{ss} in the thermal conductivity formula) M_s_D_ss(pressure ind.)/(R_0_T). This will yield a temperature only thermal conductivity which we should be able to spline lambda in temperature, which will be a speedup in the reacting flow calculation. We should have splines in temperature for D_ij, lam_i, and visc_i.
2) Of course if we go the pressure ind. route then we will need an interface to divide through by 1/P. Some of "the competition" also have interfaces to return rho_m D_m and other convenient terms for R-NS and we will probably want to add something similar.

I'm willing to start a branch with these changes but someone else will probably need to clean it up as my meta-template-fu is white belt at best.

Until then, stop worrying and learn to love c_tot * M_s...

Varis

Related Issues: 1) Do we spline the D_ij in Temperature? I get lost in the templates but it looks like we do. This is safe in the low-mach version as P is constant, but if someone is using this in a compressible code then they are using the wrong D_ij !

pbauman commented 9 years ago

@variscarey Thanks for your explanation. I'm slowly catching up, but I'll comment on more once I've merged the transport refactoring into master --- I believe this will make things a little clearer.