mutationpp / Mutationpp

The MUlticomponent Thermodynamic And Transport library for IONized gases in C++
GNU Lesser General Public License v3.0
103 stars 58 forks source link

RamshawDiffMat header & averageDiffusionCoeffs #151

Open ggange opened 3 years ago

ggange commented 3 years ago

Hello everybody,

I am currently working with the RamshawDiffMat file and I noticed a disparity with what it is written in the header:

 D_{ij} = \frac{(\delta_{ij}-y_j)}{x_j}\frac{(1-y_j)}{(1-x_j)}D_{jm}

and what it is in the source code:

for (int j = 0; j < ns; ++j) { m_Dij.col(j).fill(-Y[j]/X[j]*(1.-Y[j])*Dim(j)); m_Dij(j,j) -= m_Dij(j,j)/Y[j]; } which, to my understanding, performs:

 D_{ij} = (\delta_{ij}-y_j)\frac{(1-y_j)}{(x_j)}D_{jm}

Would it be possible to have the reference where this expression are from? In the Ramshaw's papers I found they are slightly different...

Another thing: the averageDiffusionCoeffs computed with the routine Dim() in CollisionDB are only valid for the single temperature case (Ramshaw,1993); is there already an implementation for the MultiTemperature case?

Thank you! Giuseppe

jbscoggi commented 3 years ago

Hi @ggange, thanks for the question! Indeed this is a little confusing...

I believe the formula is taken from this classical paper from Sutton and Gnoffo, Eq. 18, with the "Ramshaw" correction. Based on the Sutton and Gnoffo, we can write an approximate diffusion velocity as

V_i = -\frac{1}{x_i}\frac{1-y_i}{1-xi} D{im} \nabla x_i.

We can correct this using the Ramshaw correction, such that sum of the diffusion fluxes is zero: \sum_i y_iV_i^R = 0, where

V_i^R = V_i - \sum_j y_j V_j.

Working this out and putting it in the form V_i = \sumj D{ij} \nabla x_j will then lead to the formula for the RamshawDiffMat header and this is indeed what is implemented. Note that

D_{im} = \frac{1-x_i}{n \sum_j xj / (nA{ij})},

where I use A{ij} to be the binary diffusion coefficients (just not to confuse with the D{ij}). Note that there is a (1-xi) in the numerator which will cancel with the denominator in the expression for D{ij}. Because of this, there is an option in the function call to D_{im} that simply ignores this term for efficiency, instead of computing it twice and canceling it out.

https://github.com/mutationpp/Mutationpp/blob/acc42b991d8c7eb7e433dd8b5a46a66ca8736202/src/transport/RamshawDiffMat.cpp#L68

As for the second question, I think there might be a little confusion. This matrix is the approximate multicomponent diffusion matrix, which is different from the Ramshaw approximation of the Stefan-Maxwell equations. In terms of multiple temperatures, I believe that this matrix is already in a good form. For the Ramshaw approximation for the SM equations, I believe this is already what is actually used in the stefanMaxwell() function and is already working for 2T. You may need to check on that to be sure.

ggange commented 3 years ago

Hi @jbscoggi, thank you for the complete and clear reply. Indeed I did not know the paper by Sutton and Gnoffo, thanks for sharing it!

The first part is very clear but I still have some doubts on the second one: in the paper I attach here (In ramshaw1993.pdf ) Ramshaw defines (eq. 61) the effective binary diffusivities as:

D_i = (1 - xi) (\sum{j \neq i} \frac{zj}{D{ij}})^{-1}

where I substituted the weighting factors with x_i and z_j = p_j / p is the pressure fraction. I am quoting form the last paragraph "... In the case of equal temperatures (i.e., T_i = T for all i ), z_i reduces to x_i ..." (as in eq.4 by Sutton and Gnoffo).

I took a look at both Dim() and stefanMaxwell() and I see that only x_i is used so I assumed they were valid only for thermal equilibrium! In the stefanMaxwell() function I see a factor T_e / T_h multiplying the electron subsystem but z_e ~ x_e * T_e/T_h only if you could assume T_h as the temperature of the entire mixture, right? In my case I have a fairly low ionization degree so this reasoning holds due to the large number of neutrals, but I don't know if it can be always valid...

Am I missing something?

Thank you again for your help! Giuseppe