compressibleTwoPhaseEulerFoam: Updated thermodynamics

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/Make/options

@@ -1,4 +1,5 @@
 EXE_INC = \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/transportModels/incompressible/lnInclude \
     -IturbulenceModel \
@@ -8,6 +9,8 @@ EXE_INC = \
     -Iaveraging
 
 EXE_LIBS = \
+    -lfluidThermophysicalModels \
+    -lspecie \
     -lincompressibleTransportModels \
     -lcompressiblePhaseModel \
     -lcompressibleEulerianInterfacialModels \

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/TEqns.H

@@ -1,36 +0,0 @@
-{
-    volScalarField kByCp1("kByCp1", alpha1*(k1/Cp1/rho1 + sqr(Ct)*nut2/Prt));
-    volScalarField kByCp2("kByCp2", alpha2*(k2/Cp2/rho2 + nut2/Prt));
-
-    fvScalarMatrix T1Eqn
-    (
-        fvm::ddt(alpha1, T1)
-      + fvm::div(alphaPhi1, T1)
-      - fvm::laplacian(kByCp1, T1)
-     ==
-        heatTransferCoeff*T2/Cp1/rho1
-      - fvm::Sp(heatTransferCoeff/Cp1/rho1, T1)
-      + alpha1*(dpdt/rho1 - (fvc::ddt(K1) + fvc::div(phi1, K1)))/Cp1
-    );
-
-    fvScalarMatrix T2Eqn
-    (
-        fvm::ddt(alpha2, T2)
-      + fvm::div(alphaPhi2, T2)
-      - fvm::laplacian(kByCp2, T2)
-     ==
-        heatTransferCoeff*T1/Cp2/rho2
-      - fvm::Sp(heatTransferCoeff/Cp2/rho2, T2)
-      + alpha2*(dpdt/rho2 - (fvc::ddt(K2) + fvc::div(phi2, K2)))/Cp2
-    );
-
-    T1Eqn.relax();
-    T1Eqn.solve();
-
-    T2Eqn.relax();
-    T2Eqn.solve();
-
-    // Update compressibilities
-    psi1 = 1.0/(R1*T1);
-    psi2 = 1.0/(R2*T2);
-}

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/UEqns.H

@@ -12,12 +12,12 @@ fvVectorMatrix U2Eqn(U2, U2.dimensions()*dimVol/dimTime);
         }
         else // If not using kinetic theory is using Ct model
         {
-            nuEff1 = sqr(Ct)*nut2 + nu1;
+            nuEff1 = sqr(Ct)*nut2 + thermo1.mu()/rho1;
         }
 
         volTensorField Rc1
         (
-            "Rc1",
+            "Rc",
             (((2.0/3.0)*I)*nuEff1)*tr(gradU1T) - nuEff1*gradU1T
         );
 
@@ -52,7 +52,7 @@ fvVectorMatrix U2Eqn(U2, U2.dimensions()*dimVol/dimTime);
         volTensorField gradU2T(fvc::grad(U2)().T());
         volTensorField Rc2
         (
-            "Rc2",
+            "Rc",
             (((2.0/3.0)*I)*nuEff2)*tr(gradU2T) - nuEff2*gradU2T
         );
 

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/alphaEqn.H

@@ -1,9 +1,9 @@
-surfaceScalarField alphaPhi1("alphaPhi1", phi1);
+surfaceScalarField alphaPhi1("alphaPhi", phi1);
-surfaceScalarField alphaPhi2("alphaPhi2", phi2);
+surfaceScalarField alphaPhi2("alphaPhi", phi2);
 
 {
-    word scheme("div(phi,alpha1)");
+    word scheme("div(phi,alpha)");
-    word schemer("div(phir,alpha1)");
+    word schemer("div(phir,alpha)");
 
     surfaceScalarField phic("phic", phi);
     surfaceScalarField phir("phir", phi1 - phi2);
@@ -13,7 +13,7 @@ surfaceScalarField alphaPhi2("alphaPhi2", phi2);
         surfaceScalarField alpha1f(fvc::interpolate(alpha1));
         surfaceScalarField phipp(ppMagf*fvc::snGrad(alpha1)*mesh.magSf());
         phir += phipp;
-        phic += fvc::interpolate(alpha1)*phipp;
+        phic += alpha1f*phipp;
     }
 
     for (int acorr=0; acorr<nAlphaCorr; acorr++)
@@ -68,16 +68,22 @@ surfaceScalarField alphaPhi2("alphaPhi2", phi2);
 
         if (g0.value() > 0.0)
         {
+            surfaceScalarField alpha1f(fvc::interpolate(alpha1));
+
             ppMagf =
                 fvc::interpolate((1.0/rho1)*rAU1)
-               *fvc::interpolate
+               *g0*min(exp(preAlphaExp*(alpha1f - alphaMax)), expMax);
-                (
+
-                    g0*min(exp(preAlphaExp*(alpha1 - alphaMax)), expMax)
+            // ppMagf =
-                );
+            //     fvc::interpolate((1.0/rho1)*rAU1)
+            //    *fvc::interpolate
+            //     (
+            //         g0*min(exp(preAlphaExp*(alpha1 - alphaMax)), expMax)
+            //     );
 
             alpha1Eqn -= fvm::laplacian
             (
-                (fvc::interpolate(alpha1) + scalar(0.0001))*ppMagf,
+                alpha1f*ppMagf,
                 alpha1,
                 "laplacian(alphaPpMag,alpha1)"
             );

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/compressibleTwoPhaseEulerFoam.C

@@ -31,6 +31,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
+#include "rhoThermo.H"
 #include "nearWallDist.H"
 #include "wallFvPatch.H"
 #include "fixedValueFvsPatchFields.H"
@@ -86,7 +87,7 @@ int main(int argc, char *argv[])
             #include "alphaEqn.H"
             #include "kEpsilon.H"
             #include "interfacialCoeffs.H"
-            #include "TEqns.H"
+            #include "EEqns.H"
             #include "UEqns.H"
 
             // --- Pressure corrector loop

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/createFields.H

@@ -12,55 +12,43 @@
         )
     );
 
+    word phase1Name
+    (
+        transportProperties.found("phases")
+      ? wordList(transportProperties.lookup("phases"))[0]
+      : "phase1"
+    );
+
+    word phase2Name
+    (
+        transportProperties.found("phases")
+      ? wordList(transportProperties.lookup("phases"))[1]
+      : "phase2"
+    );
+
     autoPtr<phaseModel> phase1 = phaseModel::New
     (
         mesh,
         transportProperties,
-        "1"
+        phase1Name
     );
 
     autoPtr<phaseModel> phase2 = phaseModel::New
     (
         mesh,
         transportProperties,
-        "2"
+        phase2Name
     );
 
+    volScalarField& alpha1 = phase1();
+    volScalarField& alpha2 = phase2();
+    alpha2 = scalar(1) - alpha1;
+
     volVectorField& U1 = phase1->U();
     surfaceScalarField& phi1 = phase1->phi();
-    const dimensionedScalar& nu1 = phase1->nu();
-    const dimensionedScalar& k1 = phase1->kappa();
-    const dimensionedScalar& Cp1 = phase1->Cp();
-    dimensionedScalar rho10
-    (
-        "rho0",
-        dimDensity,
-        transportProperties.subDict("phase1").lookup("rho0")
-    );
-    dimensionedScalar R1
-    (
-        "R",
-        dimensionSet(0, 2, -2, -1, 0),
-        transportProperties.subDict("phase1").lookup("R")
-    );
 
     volVectorField& U2 = phase2->U();
     surfaceScalarField& phi2 = phase2->phi();
-    const dimensionedScalar& nu2 = phase2->nu();
-    const dimensionedScalar& k2 = phase2->kappa();
-    const dimensionedScalar& Cp2 = phase2->Cp();
-    dimensionedScalar rho20
-    (
-        "rho0",
-        dimDensity,
-        transportProperties.subDict("phase2").lookup("rho0")
-    );
-    dimensionedScalar R2
-    (
-        "R",
-        dimensionSet(0, 2, -2, -1, 0),
-        transportProperties.subDict("phase2").lookup("R")
-    );
 
     dimensionedScalar pMin
     (
@@ -69,102 +57,16 @@
         transportProperties.lookup("pMin")
     );
 
+    rhoThermo& thermo1 = phase1->thermo();
+    rhoThermo& thermo2 = phase2->thermo();
 
-    Info<< "Reading field T1" << endl;
+    volScalarField& p = thermo1.p();
-    volScalarField T1
-    (
-        IOobject
-        (
-            "T1",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
 
-    Info<< "Reading field T2" << endl;
+    volScalarField rho1("rho" + phase1Name, thermo1.rho());
-    volScalarField T2
+    const volScalarField& psi1 = thermo1.psi();
-    (
-        IOobject
-        (
-            "T2",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
 
-    volScalarField psi1
+    volScalarField rho2("rho" + phase2Name, thermo2.rho());
-    (
+    const volScalarField& psi2 = thermo2.psi();
-        IOobject
-        (
-            "psi1",
-            runTime.timeName(),
-            mesh
-        ),
-        1.0/(R1*T1)
-    );
-
-    volScalarField psi2
-    (
-        IOobject
-        (
-            "psi2",
-            runTime.timeName(),
-            mesh
-        ),
-        1.0/(R2*T2)
-    );
-    psi2.oldTime();
-
-    Info<< "Reading field p\n" << endl;
-    volScalarField p
-    (
-        IOobject
-        (
-            "p",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-    volScalarField rho1("rho1", rho10 + psi1*p);
-    volScalarField rho2("rho2", rho20 + psi2*p);
-
-
-    Info<< "Reading field alpha1\n" << endl;
-    volScalarField alpha1
-    (
-        IOobject
-        (
-            "alpha1",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-    volScalarField alpha2
-    (
-        IOobject
-        (
-            "alpha2",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        scalar(1) - alpha1
-    );
 
     volVectorField U
     (
@@ -179,27 +81,6 @@
         alpha1*U1 + alpha2*U2
     );
 
-    dimensionedScalar Cvm
-    (
-        "Cvm",
-        dimless,
-        transportProperties.lookup("Cvm")
-    );
-
-    dimensionedScalar Cl
-    (
-        "Cl",
-        dimless,
-        transportProperties.lookup("Cl")
-    );
-
-    dimensionedScalar Ct
-    (
-        "Ct",
-        dimless,
-        transportProperties.lookup("Ct")
-    );
-
     surfaceScalarField phi
     (
         IOobject
@@ -226,8 +107,6 @@
         alpha1*rho1 + alpha2*rho2
     );
 
-    #include "createRASTurbulence.H"
-
     Info<< "Calculating field DDtU1 and DDtU2\n" << endl;
 
     volVectorField DDtU1
@@ -248,6 +127,29 @@
     Info<< "Calculating field g.h\n" << endl;
     volScalarField gh("gh", g & mesh.C());
 
+    dimensionedScalar Cvm
+    (
+        "Cvm",
+        dimless,
+        transportProperties.lookup("Cvm")
+    );
+
+    dimensionedScalar Cl
+    (
+        "Cl",
+        dimless,
+        transportProperties.lookup("Cl")
+    );
+
+    dimensionedScalar Ct
+    (
+        "Ct",
+        dimless,
+        transportProperties.lookup("Ct")
+    );
+
+    #include "createRASTurbulence.H"
+
     IOdictionary interfacialProperties
     (
         IOobject
@@ -297,8 +199,8 @@
     if
     (
         !(
-            dispersedPhase == "1"
+            dispersedPhase == phase1Name
-         || dispersedPhase == "2"
+         || dispersedPhase == phase2Name
          || dispersedPhase == "both"
         )
     )
@@ -337,7 +239,7 @@
     (
         IOobject
         (
-            "rAU1",
+            "rAU" + phase1Name,
             runTime.timeName(),
             mesh,
             IOobject::NO_READ,
@@ -387,5 +289,5 @@
     );
 
     Info<< "Creating field kinetic energy K\n" << endl;
-    volScalarField K1("K1", 0.5*magSqr(U1));
+    volScalarField K1("K" + phase1Name, 0.5*magSqr(U1));
-    volScalarField K2("K2", 0.5*magSqr(U2));
+    volScalarField K2("K" + phase2Name, 0.5*magSqr(U2));

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/createRASTurbulence.H

@@ -151,7 +151,7 @@
     (
         IOobject
         (
-            "nut2",
+            "nut" + phase2Name,
             runTime.timeName(),
             mesh,
             IOobject::NO_READ,
@@ -165,13 +165,13 @@
     (
         IOobject
         (
-            "nuEff1",
+            "nuEff" + phase1Name,
             runTime.timeName(),
             mesh,
             IOobject::NO_READ,
             IOobject::NO_WRITE
         ),
-        sqr(Ct)*nut2 + nu1
+        sqr(Ct)*nut2 + thermo1.mu()/rho1
     );
 
     Info<< "Calculating field nuEff2\n" << endl;
@@ -179,11 +179,11 @@
     (
         IOobject
         (
-            "nuEff2",
+            "nuEff" + phase2Name,
             runTime.timeName(),
             mesh,
             IOobject::NO_READ,
             IOobject::NO_WRITE
         ),
-        nut2 + nu2
+        nut2 + thermo2.mu()/rho2
     );

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialCoeffs.H

@@ -38,12 +38,12 @@ volScalarField heatTransferCoeff
     volVectorField Ur(U1 - U2);
     volScalarField magUr(mag(Ur) + residualSlip);
 
-    if (dispersedPhase == "1")
+    if (dispersedPhase == phase1Name)
     {
         dragCoeff = drag1->K(magUr);
         heatTransferCoeff = heatTransfer1->K(magUr);
     }
-    else if (dispersedPhase == "2")
+    else if (dispersedPhase == phase2Name)
     {
         dragCoeff = drag2->K(magUr);
         heatTransferCoeff = heatTransfer2->K(magUr);

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/Make/options

@@ -1,7 +1,9 @@
 EXE_INC = \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I../phaseModel/lnInclude
 
 LIB_LIBS = \
-    -lcompressiblePhaseModel
+    -lcompressiblePhaseModel \
-
+    -lfluidThermophysicalModels \
+    -lspecie

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/heatTransferModels/RanzMarshall/RanzMarshall.C

@@ -72,8 +72,7 @@ Foam::tmp<Foam::volScalarField> Foam::heatTransferModels::RanzMarshall::K
 ) const
 {
     volScalarField Re(max(Ur*phase1_.d()/phase2_.nu(), scalar(1.0e-3)));
-    dimensionedScalar Prb =
+    volScalarField Prb(phase2_.rho()*phase2_.nu()*phase2_.Cp()/phase2_.kappa());
-        phase2_.rho()*phase2_.nu()*phase2_.Cp()/phase2_.kappa();
     volScalarField Nu(scalar(2) + 0.6*sqrt(Re)*cbrt(Prb));
 
     return 6.0*phase2_.kappa()*Nu/sqr(phase1_.d());

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/Make/options

@@ -1,5 +1,6 @@
 EXE_INC = \
     -I$(LIB_SRC)/foam/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I../phaseModel/lnInclude \
     -I../interfacialModels/lnInclude

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C

@@ -45,8 +45,8 @@ Foam::kineticTheoryModel::kineticTheoryModel
     phi1_(phase1.phi()),
     draga_(draga),
 
-    rho1_(phase1.rho()),
+    rho1_(phase1.rho()[0]), //***HGW
-    nu1_(phase1.nu()),
+    nu1_(phase1.nu()()[0]), //***HGW
 
     kineticTheoryProperties_
     (
@@ -120,7 +120,7 @@ Foam::kineticTheoryModel::kineticTheoryModel
     (
         IOobject
         (
-            "mu1",
+            "mu" + phase1.name(),
             U1_.time().timeName(),
             U1_.mesh(),
             IOobject::NO_READ,

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H

@@ -1,6 +1,6 @@
 {
-    rho1 = rho10 + psi1*p;
+    rho1 = thermo1.rho();
-    rho2 = rho20 + psi2*p;
+    rho2 = thermo2.rho();
 
     surfaceScalarField alpha1f(fvc::interpolate(alpha1));
     surfaceScalarField alpha2f(scalar(1) - alpha1f);
@@ -11,10 +11,10 @@
     surfaceScalarField rAlphaAU1f(fvc::interpolate(alpha1*rAU1));
     surfaceScalarField rAlphaAU2f(fvc::interpolate(alpha2*rAU2));
 
-    volVectorField HbyA1("HbyA1", U1);
+    volVectorField HbyA1("HbyA" + phase1Name, U1);
     HbyA1 = rAU1*U1Eqn.H();
 
-    volVectorField HbyA2("HbyA2", U2);
+    volVectorField HbyA2("HbyA" + phase2Name, U2);
     HbyA2 = rAU2*U2Eqn.H();
 
     mrfZones.absoluteFlux(phi1.oldTime());
@@ -38,14 +38,14 @@
 
     surfaceScalarField phiHbyA1
     (
-        "phiHbyA1",
+        "phiHbyA" + phase1Name,
         (fvc::interpolate(HbyA1) & mesh.Sf())
       + fvc::ddtPhiCorr(rAU1, alpha1, U1, phi1)
     );
 
     surfaceScalarField phiHbyA2
     (
-        "phiHbyA2",
+        "phiHbyA" + phase2Name,
         (fvc::interpolate(HbyA2) & mesh.Sf())
       + fvc::ddtPhiCorr(rAU2, alpha2, U2, phi2)
     );
@@ -96,8 +96,16 @@
     //}
     //else
     {
-        surfaceScalarField phid1("phid1", fvc::interpolate(psi1)*phi1);
+        surfaceScalarField phid1
-        surfaceScalarField phid2("phid2", fvc::interpolate(psi2)*phi2);
+        (
+            "phid" + phase1Name,
+            fvc::interpolate(psi1)*phi1
+        );
+        surfaceScalarField phid2
+        (
+            "phid" + phase2Name,
+            fvc::interpolate(psi2)*phi2
+        );
 
         pEqnComp1 =
             fvc::ddt(rho1) + psi1*correction(fvm::ddt(p))
@@ -174,13 +182,15 @@
 
     p = max(p, pMin);
 
-    rho1 = rho10 + psi1*p;
+    thermo1.correct();
-    rho2 = rho20 + psi2*p;
+    thermo2.correct();
+    rho1 = thermo1.rho();
+    rho2 = thermo2.rho();
 
     K1 = 0.5*magSqr(U1);
     K2 = 0.5*magSqr(U2);
 
-    //***HGW if (thermo.dpdt())
+    if (thermo1.dpdt())
     {
         dpdt = fvc::ddt(p);
     }

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseModel/Make/options

@@ -1,6 +1,9 @@
 EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/transportModels/incompressible/lnInclude
 
 LIB_LIBS = \
-    -lincompressibleTransportModels
+    -lincompressibleTransportModels \
+    -lfluidThermophysicalModels \
+    -lspecie

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseModel/phaseModel/phaseModel.C

@@ -37,35 +37,22 @@ Foam::phaseModel::phaseModel
     const word& phaseName
 )
 :
-    dict_
+    volScalarField
     (
-        transportProperties.subDict("phase" + phaseName)
+        IOobject
+        (
+            "alpha" + phaseName,
+            mesh.time().timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("alpha", dimless, 0)
     ),
     name_(phaseName),
-    nu_
+    phaseDict_(transportProperties.subDict(phaseName)),
-    (
+    thermo_(rhoThermo::New(mesh, phaseName)),
-        "nu",
-        dimensionSet(0, 2, -1, 0, 0),
-        dict_.lookup("nu")
-    ),
-    kappa_
-    (
-        "kappa",
-        dimensionSet(1, 1, -3, -1, 0),
-        dict_.lookup("kappa")
-    ),
-    Cp_
-    (
-        "Cp",
-        dimensionSet(0, 2, -2, -1, 0),
-        dict_.lookup("Cp")
-    ),
-    rho_
-    (
-        "rho",
-        dimDensity,
-        dict_.lookup("rho")
-    ),
     U_
     (
         IOobject
@@ -79,6 +66,8 @@ Foam::phaseModel::phaseModel
         mesh
     )
 {
+    thermo_->validate("phaseModel " + phaseName, "h", "e");
+
     const word phiName = "phi" + phaseName;
 
     IOobject phiHeader
@@ -147,7 +136,7 @@ Foam::phaseModel::phaseModel
 
     dPtr_ = diameterModel::New
     (
-        dict_,
+        phaseDict_,
         *this
     );
 }

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseModel/phaseModel/phaseModel.H

@@ -36,6 +36,7 @@ SourceFiles
 #include "dimensionedScalar.H"
 #include "volFields.H"
 #include "surfaceFields.H"
+#include "rhoThermo.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -50,25 +51,18 @@ class diameterModel;
 \*---------------------------------------------------------------------------*/
 
 class phaseModel
+:
+    public volScalarField
 {
     // Private data
 
-        dictionary dict_;
-
         //- Name of phase
         word name_;
 
-        //- Kinematic viscosity
+        dictionary phaseDict_;
-        dimensionedScalar nu_;
 
-        //- Thermal conductivity
+        //- Thermophysical properties
-        dimensionedScalar kappa_;
+        autoPtr<rhoThermo> thermo_;
-
-        //- Heat capacity
-        dimensionedScalar Cp_;
-
-        //- Density
-        dimensionedScalar rho_;
 
         //- Velocity
         volVectorField U_;
@@ -116,24 +110,34 @@ public:
 
         tmp<volScalarField> d() const;
 
-        const dimensionedScalar& nu() const
+        tmp<volScalarField> nu() const
         {
-            return nu_;
+            return thermo_->mu()/thermo_->rho();
         }
 
-        const dimensionedScalar& kappa() const
+        tmp<volScalarField> kappa() const
         {
-            return kappa_;
+            return thermo_->kappa();
         }
 
-        const dimensionedScalar& Cp() const
+        tmp<volScalarField> Cp() const
         {
-            return Cp_;
+            return thermo_->Cp();
         }
 
-        const dimensionedScalar& rho() const
+        const volScalarField& rho() const
         {
-            return rho_;
+            return thermo_->rho();
+        }
+
+        const rhoThermo& thermo() const
+        {
+            return thermo_();
+        }
+
+        rhoThermo& thermo()
+        {
+            return thermo_();
         }
 
         const volVectorField& U() const

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/turbulenceModel/kEpsilon.H

@@ -61,4 +61,4 @@ if (turbulence)
     #include "wallViscosity.H"
 }
 
-nuEff2 = nut2 + nu2;
+nuEff2 = nut2 + thermo2.mu()/rho2;

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/turbulenceModel/wallFunctions.H

@@ -4,7 +4,6 @@
     scalar Cmu25 = ::pow(Cmu.value(), 0.25);
     scalar Cmu75 = ::pow(Cmu.value(), 0.75);
     scalar kappa_ = kappa.value();
-    scalar nu2_ = nu2.value();
 
     const fvPatchList& patches = mesh.boundary();
 
@@ -30,6 +29,8 @@
     forAll(patches, patchi)
     {
         const fvPatch& currPatch = patches[patchi];
+        const scalarField& mu2_ = thermo2.mu().boundaryField()[patchi];
+        const scalarField& rho2_ = rho2.boundaryField()[patchi];
 
         if (isA<wallFvPatch>(currPatch))
         {
@@ -52,7 +53,8 @@
                     /(kappa_*y[patchi][facei]);
 
                 G[faceCelli] +=
-                    (nut2w[facei] + nu2_)*magFaceGradU[facei]
+                    (nut2w[facei] + mu2_[facei]/rho2_[facei])
+                   *magFaceGradU[facei]
                    *Cmu25*::sqrt(k[faceCelli])
                    /(kappa_*y[patchi][facei]);
             }

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/turbulenceModel/wallViscosity.H

@@ -2,13 +2,14 @@
     scalar Cmu25 = ::pow(Cmu.value(), 0.25);
     scalar kappa_ = kappa.value();
     scalar E_ = E.value();
-    scalar nu2_ = nu2.value();
 
     const fvPatchList& patches = mesh.boundary();
 
     forAll(patches, patchi)
     {
         const fvPatch& currPatch = patches[patchi];
+        const scalarField& mu2_ = thermo2.mu().boundaryField()[patchi];
+        const scalarField& rho2_ = rho2.boundaryField()[patchi];
 
         if (isA<wallFvPatch>(currPatch))
         {
@@ -20,11 +21,14 @@
 
                 // calculate yPlus
                 scalar yPlus =
-                    Cmu25*y[patchi][facei]*::sqrt(k[faceCelli])/nu2_;
+                    Cmu25*y[patchi][facei]*::sqrt(k[faceCelli])
+                   /(mu2_[facei]/rho2_[facei]);
 
                 if (yPlus > 11.6)
                 {
-                    nutw[facei] = nu2_*(yPlus*kappa_/::log(E_*yPlus) - 1);
+                    nutw[facei] =
+                        (mu2_[facei]/rho2_[facei])
+                       *(yPlus*kappa_/::log(E_*yPlus) - 1);
                 }
                 else
                 {