[ChrisA]

I originally posted this issue here: http://www.cfd-online.com/Forums/ope...tml#post414308 but after reading this post: http://www.cfd-online.com/Forums/ope...t-mixture.html and this bug report: http://www.openfoam.org/mantisbt/view.php?id=327 I'm beginning to think the issue is inherent in the code and therefor not suitable for the running and solving section. I'd like to confirm this fact and figure out how to stop the behavior...

Anyway, what appears to be happening is that the viscosity calculation for a mixture is behaving in a strange way. Changing the parameters for a gas that is not present in the mixture for a given cell is changing the viscosity of that cell. It seems like the behavior is similar to that described in the Janaf bug report above, where there's some summing of coefficients prior to calculating the viscosity.

This suspicion is further reinforced by the fact that all the sutherland coefficients for all of the tutorials are identical. This would effectively make the mixture viscosity independent of composition... this seems odd given that the cellMixture functions exist in the mxitureThermos from combustionThermo.

I've tried playing with the operators in the viscosity library I'm using but I'm afraid that's had no effect and I'm not sure that it is the correct approach to this issue. Can anyone offer any advice regarding this? I've managed to get my case solving properly by writing the viscosity equations into the solver code, doing the mixing there and not using the transport libraries but this is obviously a non-ideal solution.

[dkxls]

Hi Chris!

After I ran into the same problem with wrong viscosities I looked into it a big deeper. I am quite confident that it is a bug in OpenFOAM, more specific in the sutherlandTransport class.
As you pointed out already, the mixture averaging in OpenFOAM is done by averaging the polynomial or whatever coefficients and then computing the properties only once with these averaged coefficients.
The averaging is done in multiComponentMixture.C and the important lines are:

    mixture_ = Y_[0][celli]/speciesData_[0].W()*speciesData_[0];

    for (label n=1; n<Y_.size(); n++)
    {
        mixture_ += Y_[n][celli]/speciesData_[n].W()*speciesData_[n];
    }

Now in the cases of sutherlandTransport the += operator is not defined. I actually don't know why this does not lead to an error during compilation, but what it apparently leads to is a viscosity computation that only uses the coefficients for the first species!

I have quickly coded a small utility that computes the species data for

sutherlandTransport<specieThermo<janafThermo<perfectGas> > >

in two ways.
1) explicitly mass averaged
2) the OpenFOAM way (using the code above)

If I replace "+=" by "mixture_ = mixture_ + ...", I get sane values for the viscosity. The tool is attached, feel free to test it yourself.

Now the whole problem only shows up if one uses different coefficients for the sutherlandTransport, which is not done in any tutorial.
As species data are often needed for combustion simulation, the thermo/ transport coefficients come usually from the CHEMKIN reader in OpenFOAM, which sets the Sutherland coefficients for all species to values that should resemble nitrogen. This is also a quite questionable practice, but also most likely the reason why this bug didn't show up earlier.

Best,
Armin

[dkxls]

The bug got resolved in 2.1.x and 2.2.x. Here is the link to the bug-report:
http://www.openfoam.org/mantisbt/view.php?id=799

[ChrisA]

Thanks Armin, for looking into it and confirming my suspicions, also for making the bug report, it looks like the revisions havn't hit a source pack yet, at least from what I can tell.

This leads me to another issue, the viscosity model I'd like to use (kinetic theory based) is non linear with respect to the coefficients. As a result I cannot use molar averaging of the coefficients and hope to get real results. I imagine to get around this I need to re-define my += operator in such a way that the viscosity, kappa and alpha are all explicitly calculated for each specie then mass averaged in multiComponentMixture.C. Is it possible to do this? If so, any idea on how? This level of coding is going a bit over my head.

I guess my final option, if the above isn't possible and which I'd like to avoid, would be using polynomial transport with curve fits for my functions...

Edit:

I should note that I'm using OF 2.1.1.

Edit2:

Follow up question regarding Kappa and Alpha calculations, when Cp, Cv and R are called in sutherlandTransport (or any transport), does it receive molar averaged values or per-specie values?

Also, in polynomialTransport... the coefficients are being divided by MW before being passed on, this implies molar coefficients? but if the mixture.C also divides by MW does this not lead to incorrect values?

[dkxls]

I guess it still takes some time until the point release (2.2.1) comes, so you are left up to a git installation if you want to have the fix now.

Originally Posted by ChrisA

I imagine to get around this I need to re-define my += operator in such a way that the viscosity, kappa and alpha are all explicitly calculated for each specie then mass averaged in multiComponentMixture.C. Is it possible to do this?

I don't see how you can do that in the += operator. I'm no C++ expert, but I would say it's not possible.

Originally Posted by ChrisA

I guess my final option, if the above isn't possible and which I'd like to avoid, would be using polynomial transport with curve fits for my functions...

I don't think that polynomial fits solve your problem with the mixture viscosity. At least if you use different polynomial coefficients for each species you still have to combine them somehow and mass averaging is rather a rough approximation.
This leaves you up with writing a mixture viscosity model, e.g. according to Wilke's formula.
Of course, you can also use a polynomial fit for the mixture viscosity, thus use only one polynomial. But I don't think that's what you have in mind as it is probably an even rougher approximation.

Originally Posted by ChrisA

Follow up question regarding Kappa and Alpha calculations, when Cp, Cv and R are called in sutherlandTransport (or any transport), does it receive molar averaged values or per-specie values?

Yes, given you use the multiComponentMixture class.

Originally Posted by ChrisA

Also, in polynomialTransport... the coefficients are being divided by MW before being passed on, this implies molar coefficients? but if the mixture.C also divides by MW does this not lead to incorrect values?

No, I don't see there anything wrong in the way how the coefficients are calculated.

Hope you have an idea of how to go on, and as Henry Weller pointed out in the bug report, the new thermo structure in 2.2.x facilitates the implementation of different mixture models.

So things got easier in that perspective.

-Armin

[ChrisA]

Thanks again for taking the time to answer my questions. The mixing in OF certainly is a tough nut to crack with many intricacies.

Originally Posted by dkxls

I don't see how you can do that in the += operator. I'm no C++ expert, but I would say it's not possible.

I'm probably missing one of the finer details in the multiComponentMixture.C cellMixture function but it seems to me that when it runs through the loop for each specie, the mixture can be calculated as (using mu for an example)

mixture mu += mu(T),molarbased*Y[i]/MW_i

where += just adds the previous mixture mu quantity to the new mixture mu quantity, ie. functions as the original definition of += from back in C++ 101. And the mu function just gets called each time for the per specie coefficients, as if it were being called by a pure mixture thermotype like psiThermo. I think it does function this for rho in perfectGas? but if i try to copy that format I still get incorrect results.

I don't think that polynomial fits solve your problem with the mixture viscosity. At least if you use different polynomial coefficients for each species you still have to combine them somehow and mass averaging is rather a rough approximation.
This leaves you up with writing a mixture viscosity model, e.g. according to Wilke's formula.
Of course, you can also use a polynomial fit for the mixture viscosity, thus use only one polynomial. But I don't think that's what you have in mind as it is probably an even rougher approximation.

I have done surface plot fits for viscosity before using PVT data for a fluid who's composition wasn't constant but that wasn't in OF and defiantly isn't the solution I'm looking for with this problem :P

The idea behind the mass averaging is I'm trying to match up to previous work done with FLUENT before moving forward with a more accurate model. You know, show that OF can achieve the baseline complexity and then go from there. The FLUENT work used a kinetic theory user defined function for the individual viscosity then mass weighted them together.

The other issue is I'm simulating somewhat outside of the sutherland range of applicability, as least as far as my thesis supervisor is concerned.

Yes, given you use the multiComponentMixture class.

I'm guessing this was yes to the molar averaged constants, but just confirming since it was a double ended question. This further exacerbates the issue as I'm using a thermal conductivity model that does not lend itself to that either.

No, I don't see there anything wrong in the way how the coefficients are calculated.

Alright, thanks for the confirmation, at any rate I've decided to fit a thermal conductivity polynomial and continue to use my mu definition directly in the solver. Once I start matching I'll look into using the existing code with a better mixing law.

Hope you have an idea of how to go on, and as Henry Weller pointed out in the bug report, the new thermo structure in 2.2.x facilitates the implementation of different mixture models.

So things got easier in that perspective.

I'm certainly going to investigate the 2.2.x structure, but right now there's a pressing paper deadline and not enough time to learn how they implemented that. From what I can tell though they still insist on using the molar averaged coefficients at the bottom level for the thermodynamic quantities for mixtures, forcing the user into having coefficients that are strictly linear operators in the function. Or am I wrong about that?

[dkxls]

Originally Posted by ChrisA

I'm probably missing one of the finer details in the multiComponentMixture.C cellMixture function but it seems to me that when it runs through the loop for each specie, the mixture can be calculated as (using mu for an example)

mixture mu += mu(T),molarbased*Y[i]/MW_i

where += just adds the previous mixture mu quantity to the new mixture mu quantity, ie. functions as the original definition of += from back in C++ 101. And the mu function just gets called each time for the per specie coefficients, as if it were being called by a pure mixture thermotype like psiThermo. I think it does function this for rho in perfectGas? but if i try to copy that format I still get incorrect results.

Indeed you are missing there the finer details. It is not done the way you describe it. As Niklas explained it here
http://www.cfd-online.com/Forums/ope...tml#post220990 ,
the for-loop in the multiComponentMixture class does not mass average viscosities, heat capacities or whatever. What essentially is done, is the construction of a new pseudo-species for each cell (called "mixture_"), which has the properties of the cell mixture. These properties are obtained by summing up the coefficients on a molar fraction basis (e.g. JANAF/NASA polynomial coeffs for Cp, H, S; Sutherland's coeffs for mu/alpha; ...).
The actual summation is done by the += operator in the respective classes, e.g. janafThermoI.H, sutherlandTransportI.H, ... .

Originally Posted by ChrisA

The other issue is I'm simulating somewhat outside of the sutherland range of applicability, as least as far as my thesis supervisor is concerned.

Well, with Sutherland's approximation you are kind of limited to a not-too-broad temperature range, as the error becomes rather big at either end (or in the middle part). You can minimize the error by optimizing your coefficients for a given temperature range, but accurate results for a broad temperature range you will only get with polynomial fits.

Originally Posted by ChrisA

I'm guessing this was yes to the molar averaged constants, but just confirming since it was a double ended question. This further exacerbates the issue as I'm using a thermal conductivity model that does not lend itself to that either.

Sorry for the confusion, but I guess the answer is given now properly above.

[ChrisA]

Hmm, well thank you very much for taking the time to clarify everything for me. The whole thing is somewhat disappointing to hear, the only real proper solution seems to be writing my own mixture class in which explicit per specie calculations and mass/whatever else weighting occurred. Oh well, polynomials it is!

On that topic, thank you very much for the application as well, defiantly made per-specie polynomial fits easier to make... since the formulations had already been coded by summer students no sense not using their contributions in some manner

[dkxls]

Another thing to keep in mind is the thermal diffusivity for enthalpy. Up to OF 2.1.x some Cpbar is used to compute alpha, i.e. alpha = kappa/Cpbar.

See also here: http://www.cfd-online.com/Forums/ope...tml#post399357

This Cpbar is in my opinion just wrong for sutherlandTransport or polynomialTransport and it also got replaced by a Cp from the thermo class in OF 2.2.x. So one more reason to switch to the newest version.

However, one should keep in mind that the thermal diffusivity for enthalpy is a mixture quantity and hence alphah = kappaMix/CpMix. (alphah got newly introduced in OF 2.2.x, previous versions call it just alpha).
In the present implementation, alphah is calculated in sutherlandTransport or polynomialTransport. If you implement your own mixture model you should keep in mind to re-implement the alphah calculation, using your new conductivity mixture model.

[ChrisA]

Good point, I compared k/Cp to k/Cpbar for what I'm doing and the Cpbar method results in an upwards of 4% difference. It also cuts down on the calculations as I only have to fit kmix and then calculate alpha from kmix and Cpmix like you said.

The difference is of course zero if Cp is constant.

[ziemowitzima]

Dear Foamers,

Can anybody tell me how in OF the viscosity of the mixture is calculated ?
Is it just linear superposition ? If yes then it seems to be wrong...

Where is it defined/calculated in the source code ?

best and thanks

ZMM

[dkxls]

Originally Posted by ziemowitzima

Can anybody tell me how in OF the viscosity of the mixture is calculated ?
Is it just linear superposition ? If yes then it seems to be wrong...

Where is it defined/calculated in the source code ?

Please read my posts above:
http://www.cfd-online.com/Forums/ope...tml#post416223
http://www.cfd-online.com/Forums/ope...tml#post416415

What you will essentially end up with is not even a strictly mass averaged mixture viscosity, but something else (i.e. a "Sutherland coefficient averaged mixture").
Either way, it's a rough approximation, but as long as you are doing RANS, there is probably not much to worry about that.
(I use the mixture models due to Wilke/Mathur for my LES.)


P.S.: I hope you are using either 2.2.x or 2.3.0, to benefit from the bug-fixes in the viscosity classes.

[ziemowitzima]

Thank you for your fast answer !
I read the whole thread more carefully, and I see that the bug you reported was fixed in OF2.2.1.

Although, if I understand correctly, in OF 2.2.1 viscosity and thermal conductivity is still calculated as linear averaging ?

In your "bug report"you said that you implemented your own thermoclass which uses
Wilke's formula (viscosity) and Mathur (thermal conductivity). Could you share it ?

Thanks
ZMM

[epi_c]

Hello All,

I have a question about sutherland coefficients in OpenFOAM. In tutorials like reactingFoam/counterFlowFlame2D, dynamic viscosity is calculated according to sutherland law, and the foamChemistryReader is used, so the sutherland coefficients are explicitly set in constant/thermo.compressibleGas as:

O2
{
    specie
    {
        nMoles          1;
        molWeight       31.9988;
    }
    thermodynamics
    {
        Tlow            200;
        Thigh           5000;
        Tcommon         1000;
        highCpCoeffs    ( 3.69758 0.00061352 -1.25884e-07 1.77528e-11 -1.13644e-15 -1233.93 3.18917 );
        lowCpCoeffs     ( 3.21294 0.00112749 -5.75615e-07 1.31388e-09 -8.76855e-13 -1005.25 6.03474 );
    }
    transport
    {
        As              1.67212e-06;
        Ts              170.672;
    }
}

But in sprayFoam/aachenBomb, the chemkinReader is used to read reaction mechanism and thermal data specified as:

CHEMKINFile     "$FOAM_CASE/chemkin/chem.inp";

CHEMKINThermoFile "$FOAM_CASE/chemkin/therm.dat";

Unfortunately, I can't find the corresponding coefficients As and Ts for each species like in tutorial reactingFoam/counterFlowFlame2D.

Could anybody tell me why? Thanks very much.

Best regards

Francis

[dkxls]

For the OpenFOAM CHEMKIN reader the sutherland coefficients are hard-coded into the chemkin lexer. I believe the coefficients are somewhat close to the ones of nitrogen.
The code is here:

src/thermophysicalModels/reactionThermo/chemistryReaders/chemkinReader/chemkinLexer.L

If you want to use coefficients per species, you have to first convert the CHEMKIN files to OpenFOAM format (chemkinToFoam) and then use the foamChemistryReader.
This way you can specify the Sutherland coefficients in the thermophysicalProperties, as done in the counterFlowFlame2D tutorial.

Good luck!

-Armin