Unphysical temperature reduction when mixing different density gases

[aansari2]

I tried the depthCharge2D example with compressibleInterFlow. I replaced steam with helium gas.

    specie
    {
        molWeight   4.002;
    }
    thermodynamics
    {
        Cv          5190.0;
        Hf          0;
    }
    transport
    {
        mu          1.96e-05;
        Pr          0.71;
    }

At the gas-gas interface, a nearly 100 kelvin temperature reduction is seen. Here is the following result at t = 0.01. My guess is this likely has to do with some confusion between mass-fraction and volume-fraction.

[HenningScheufler]

depthCharge2D.zip

The behavior is nonphysical. The attached testcase makes it more obvious by setting a uniform alpha1( =0), pressure and temperature field.

If you run the testcase above with the current solver and compare it to the results with following lines

//   + (
//         fvc::div(fvc::absolute(phi, U), p)()() // - contErr/rho*p
//        + (fvc::ddt(rho, K) + fvc::div(rhoPhi, K))()() - contErr*K
//     )
//    *(
//        alpha1()/mixture.thermo1().Cv()()
//      + alpha2()/mixture.thermo2().Cv()()
//     )

commented out in TEqn.H. This indicates a problem in the temperature equation.

[HenningScheufler]

The error is a result of a discretization error in the energy/temperature equation. bouyant.zip

The energy equation of the solver rhoReactingFoam changes based on the thermodynamic model:

       + (
             he.name() == "e"
           ? fvc::div
             (
                 fvc::absolute(phi/fvc::interpolate(rho), U),
                 p,
                 "div(phiv,p)"
             )
           : -dpdt
         )

With the thermodynamic model energy sensibleEnthalpy; we get a realistic species distribution (red has a higher density heavier)

With the thermodynamic model energy sensibleInternalEnergy; the simulation crashes after 0.01 seconds and the temperature field is nonpysical

[aansari2]

Thanks for the prompt response. The reason I said it might be mass-fraction/volume fraction problem in the energy equation is because I noticed inconsistency between gas mixture density and gas concentrations. Here pure air is ~1.2 kg/m^3 and pure steam is 0.16 kg/m^3. The density should be equal to 1.2 air.air + 0.16 steam.air but that doesn't seem to be the case.

[HenningScheufler]

The solvers is based on the mass fraction Y and the formulation is similar to the reaction based solver (rhoReactionFoam ....)

The basic idea is to have a concentration field for every species and phase, so with three species and two phases we would have 6 fields in total.

https://www.openfoam.com/documentation/guides/v2112/api/reactingEuler_2multiphaseSystem_2phaseModel_2MultiComponentPhaseModel_2MultiComponentPhaseModel_8C_source.html

If you want volume-based concentration, you need to create a new class e.g. MultiComponentPhaseModelVolume (start by copying MultiComponentPhaseModel and rename it). In order to use the class you need to compile it with the macro #include "makePhaseTypes.H" (examples can be found in phaseModels.C)

An example of volume based concentration is here.

https://www.openfoam.com/documentation/guides/v2112/api/multiphaseInter_2phasesSystem_2phaseModel_2MultiComponentPhaseModel_2MultiComponentPhaseModel_8C_source.html