Merge branch 'master' of /home/dm4/OpenFOAM/OpenFOAM-dev

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/Allwclean

@@ -2,9 +2,9 @@
 cd ${0%/*} || exit 1    # run from this directory
 set -x
 
-wclean libso phaseModel
+wclean libso twoPhaseSystem
 wclean libso interfacialModels
-wclean libso kineticTheoryModels
+wclean libso phaseIncompressibleTurbulenceModels
 wclean
 
 # ----------------------------------------------------------------- end-of-file

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/Allwmake

@@ -2,9 +2,10 @@
 cd ${0%/*} || exit 1    # run from this directory
 set -x
 
-wmake libso phaseModel
+wmakeLnInclude interfacialModels
+wmake libso twoPhaseSystem
 wmake libso interfacialModels
-wmake libso kineticTheoryModels
+wmake libso phaseIncompressibleTurbulenceModels
 wmake
 
 # ----------------------------------------------------------------- end-of-file

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/EEqns.H

@@ -1,22 +1,12 @@
 {
-    volScalarField k1
-    (
-        "k1",
-        alpha1*(thermo1.alpha()/rho1 + sqr(Ct)*nut2*thermo1.CpByCpv()/Prt)
-    );
-
-    volScalarField k2
-    (
-        "k2",
-        alpha2*(thermo2.alpha()/rho2 + nut2*thermo2.CpByCpv()/Prt)
-    );
-
     volScalarField& he1 = thermo1.he();
     volScalarField& he2 = thermo2.he();
 
     volScalarField Cpv1(thermo1.Cpv());
     volScalarField Cpv2(thermo2.Cpv());
 
+    volScalarField heatTransferCoeff(fluid.heatTransferCoeff());
+
     fvScalarMatrix he1Eqn
     (
         fvm::ddt(alpha1, he1) + fvm::div(alphaPhi1, he1)
@@ -32,7 +22,8 @@
           : -alpha1*dpdt
         )/rho1
 
-      - fvm::laplacian(k1, he1)
+      //***HGW- fvm::laplacian(alpha1*turbulence1->alphaEff(), he1)
+      - fvm::laplacian(alpha1*turbulence1->nuEff(), he1)
      ==
         heatTransferCoeff*(thermo2.T() - thermo1.T())/rho1
       + heatTransferCoeff*he1/Cpv1/rho1
@@ -54,7 +45,8 @@
           : -alpha2*dpdt
         )/rho2
 
-      - fvm::laplacian(k2, he2)
+      //***HGW- fvm::laplacian(alpha2*turbulence2->alphaEff(), he2)
+      - fvm::laplacian(alpha2*turbulence2->nuEff(), he2)
      ==
         heatTransferCoeff*(thermo1.T() - thermo2.T())/rho2
       + heatTransferCoeff*he2/Cpv2/rho2

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/Make/options

@@ -2,18 +2,22 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/transportModels/incompressible/lnInclude \
-    -IturbulenceModel \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModel/lnInclude \
-    -IkineticTheoryModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/phaseIncompressible/lnInclude \
+    -IphaseIncompressibleTurbulenceModels/lnInclude \
     -IinterfacialModels/lnInclude \
-    -IphaseModel/lnInclude \
+    -ItwoPhaseSystem/lnInclude \
     -Iaveraging
 
 EXE_LIBS = \
     -lfluidThermophysicalModels \
     -lspecie \
+    -lturbulenceModels \
+    -lincompressibleTurbulenceModels \
+    -lphaseIncompressibleTurbulenceModels \
     -lincompressibleTransportModels \
-    -lcompressiblePhaseModel \
+    -lcompressibleTwoPhaseSystem \
     -lcompressibleEulerianInterfacialModels \
-    -lcompressibleKineticTheoryModel \
     -lfiniteVolume \
     -lmeshTools

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/UEqns.H

@@ -5,31 +5,12 @@ mrfZones.correctBoundaryVelocity(U);
 fvVectorMatrix U1Eqn(U1, U1.dimensions()*dimVol/dimTime);
 fvVectorMatrix U2Eqn(U2, U2.dimensions()*dimVol/dimTime);
 
+volScalarField dragCoeff(fluid.dragCoeff());
+
 {
+    volVectorField liftForce(fluid.liftForce(U));
+
     {
-        volTensorField gradU1T(fvc::grad(U1)().T());
-
-        if (kineticTheory.on())
-        {
-            kineticTheory.solve(gradU1T);
-            nuEff1 = kineticTheory.mu1()/rho1;
-        }
-        else // If not using kinetic theory is using Ct model
-        {
-            nuEff1 = sqr(Ct)*nut2 + thermo1.mu()/rho1;
-        }
-
-        volTensorField Rc1
-        (
-            "Rc",
-            (((2.0/3.0)*I)*nuEff1)*tr(gradU1T) - nuEff1*gradU1T
-        );
-
-        if (kineticTheory.on())
-        {
-            Rc1 -= ((kineticTheory.lambda()/rho1)*tr(gradU1T))*tensor(I);
-        }
-
         U1Eqn =
         (
             fvm::ddt(alpha1, U1)
@@ -38,31 +19,23 @@ fvVectorMatrix U2Eqn(U2, U2.dimensions()*dimVol/dimTime);
             // Compressibity correction
           - fvm::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), U1)
 
-          + Cvm*rho2*alpha1*alpha2/rho1*
+          + fluid.Cvm()*rho2*alpha1*alpha2/rho1*
             (
                 fvm::ddt(U1)
               + fvm::div(phi1, U1)
               - fvm::Sp(fvc::div(phi1), U1)
             )
 
-          - fvm::laplacian(alpha1*nuEff1, U1)
+          + turbulence1->divDevReff(U1)
-          + fvc::div(alpha1*Rc1)
          ==
           - fvm::Sp(dragCoeff/rho1, U1)
-          - alpha1*alpha2/rho1*(liftForce - Cvm*rho2*DDtU2)
+          - alpha1*alpha2/rho1*(liftForce - fluid.Cvm()*rho2*DDtU2)
         );
-        mrfZones.addCoriolis(alpha1*(1 + Cvm*rho2*alpha2/rho1), U1Eqn);
+        mrfZones.addCoriolis(alpha1*(1 + fluid.Cvm()*rho2*alpha2/rho1), U1Eqn);
         U1Eqn.relax();
     }
 
     {
-        volTensorField gradU2T(fvc::grad(U2)().T());
-        volTensorField Rc2
-        (
-            "Rc",
-            (((2.0/3.0)*I)*nuEff2)*tr(gradU2T) - nuEff2*gradU2T
-        );
-
         U2Eqn =
         (
             fvm::ddt(alpha2, U2)
@@ -71,20 +44,18 @@ fvVectorMatrix U2Eqn(U2, U2.dimensions()*dimVol/dimTime);
             // Compressibity correction
           - fvm::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), U2)
 
-          + Cvm*rho2*alpha1*alpha2/rho2*
+          + fluid.Cvm()*rho2*alpha1*alpha2/rho2*
             (
                 fvm::ddt(U2)
               + fvm::div(phi2, U2)
               - fvm::Sp(fvc::div(phi2), U2)
             )
-
+          + turbulence2->divDevReff(U2)
-          - fvm::laplacian(alpha2*nuEff2, U2)
-          + fvc::div(alpha2*Rc2)
          ==
           - fvm::Sp(dragCoeff/rho2, U2)
-          + alpha1*alpha2/rho2*(liftForce + Cvm*rho2*DDtU1)
+          + alpha1*alpha2/rho2*(liftForce + fluid.Cvm()*rho2*DDtU1)
         );
-        mrfZones.addCoriolis(alpha2*(1 + Cvm*rho2*alpha1/rho2), U2Eqn);
+        mrfZones.addCoriolis(alpha2*(1 + fluid.Cvm()*rho2*alpha1/rho2), U2Eqn);
         U2Eqn.relax();
     }
 }

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/alphaEqn.H

@@ -1,6 +1,3 @@
-surfaceScalarField alphaPhi1("alphaPhi" + phase1Name, phi1);
-surfaceScalarField alphaPhi2("alphaPhi" + phase2Name, phi2);
-
 {
     word alphaScheme("div(phi," + alpha1.name() + ')');
     word alpharScheme("div(phir," + alpha1.name() + ')');
@@ -8,12 +5,23 @@ surfaceScalarField alphaPhi2("alphaPhi" + phase2Name, phi2);
     surfaceScalarField phic("phic", phi);
     surfaceScalarField phir("phir", phi1 - phi2);
 
-    if (g0.value() > 0.0)
+    surfaceScalarField alpha1f(fvc::interpolate(max(alpha1, 0.0)));
+
+    tmp<surfaceScalarField> pPrimeByA;
+
+    if (implicitPhasePressure)
     {
-        surfaceScalarField alpha1f(fvc::interpolate(alpha1));
+        pPrimeByA =
-        surfaceScalarField phipp(ppMagf*fvc::snGrad(alpha1)*mesh.magSf());
+            fvc::interpolate((1.0/rho1)*rAU1*turbulence1().pPrime())
-        phir += phipp;
+          + fvc::interpolate((1.0/rho2)*rAU2*turbulence2().pPrime());
-        phic += alpha1f*phipp;
+
+        surfaceScalarField phiP
+        (
+            pPrimeByA()*fvc::snGrad(alpha1, "bounded")*mesh.magSf()
+        );
+
+        phic += alpha1f*phiP;
+        phir += phiP;
     }
 
     for (int acorr=0; acorr<nAlphaCorr; acorr++)
@@ -92,7 +100,7 @@ surfaceScalarField alphaPhi2("alphaPhi" + phase2Name, phi2);
                 alphaPhic1,
                 Sp,
                 Su,
-                (g0.value() > 0 ? alphaMax : 1),
+                1,
                 0
             );
 
@@ -104,61 +112,19 @@ surfaceScalarField alphaPhi2("alphaPhi" + phase2Name, phi2);
             {
                 alphaPhi1 = alphaPhic1;
             }
-
-            /*
-            // Legacy semi-implicit and potentially unbounded form
-            fvScalarMatrix alpha1Eqn
-            (
-                fvm::ddt(alpha1)
-              + fvm::div(phic, alpha1, alphaScheme)
-              + fvm::div
-                (
-                    -fvc::flux(-phir, alpha2, alpharScheme),
-                    alpha1,
-                    alpharScheme
-                )
-             ==
-               fvm::Sp(Sp, alpha1) + Su
-            );
-
-            alpha1Eqn.relax();
-            alpha1Eqn.solve();
-
-            if (nAlphaSubCycles > 1)
-            {
-                alphaPhi1 += (runTime.deltaT()/totalDeltaT)*alpha1Eqn.flux();
-            }
-            else
-            {
-                alphaPhi1 = alpha1Eqn.flux();
-            }
-            */
         }
 
-        if (g0.value() > 0.0)
+        if (implicitPhasePressure)
         {
-            surfaceScalarField alpha1f(fvc::interpolate(alpha1));
-
-            ppMagf =
-                fvc::interpolate((1.0/rho1)*rAU1)
-               *g0*min(exp(preAlphaExp*(alpha1f - alphaMax)), expMax);
-
             fvScalarMatrix alpha1Eqn
             (
                 fvm::ddt(alpha1) - fvc::ddt(alpha1)
-              - fvm::laplacian
+              - fvm::laplacian(alpha1f*pPrimeByA, alpha1, "bounded")
-                (
-                    alpha1f*ppMagf,
-                    alpha1,
-                    "laplacian(alpha1PpMag,alpha1)"
-                )
             );
 
             alpha1Eqn.relax();
             alpha1Eqn.solve();
 
-            #include "packingLimiter.H"
-
             alphaPhi1 += alpha1Eqn.flux();
         }
 
@@ -173,4 +139,4 @@ surfaceScalarField alphaPhi2("alphaPhi" + phase2Name, phi2);
     }
 }
 
-rho = alpha1*rho1 + alpha2*rho2;
+rho = fluid.rho();

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/compressibleTwoPhaseEulerFoam.C

@@ -34,21 +34,11 @@ Description
 #include "MULES.H"
 #include "subCycle.H"
 #include "rhoThermo.H"
-#include "nearWallDist.H"
+#include "twoPhaseSystem.H"
-#include "wallFvPatch.H"
-#include "fixedValueFvsPatchFields.H"
-#include "Switch.H"
-
-#include "IFstream.H"
-#include "OFstream.H"
-
-#include "phaseModel.H"
 #include "dragModel.H"
 #include "heatTransferModel.H"
-#include "kineticTheoryModel.H"
+#include "PhaseIncompressibleTurbulenceModel.H"
-
 #include "pimpleControl.H"
-
 #include "IOMRFZoneList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -62,7 +52,6 @@ int main(int argc, char *argv[])
     #include "readGravitationalAcceleration.H"
     #include "createFields.H"
     #include "createMRFZones.H"
-    #include "readPPProperties.H"
     #include "initContinuityErrs.H"
     #include "readTimeControls.H"
     #include "CourantNos.H"
@@ -87,8 +76,6 @@ int main(int argc, char *argv[])
         while (pimple.loop())
         {
             #include "alphaEqn.H"
-            #include "kEpsilon.H"
-            #include "interfacialCoeffs.H"
             #include "EEqns.H"
             #include "UEqns.H"
 
@@ -99,6 +86,12 @@ int main(int argc, char *argv[])
             }
 
             #include "DDtU.H"
+
+            if (pimple.turbCorr())
+            {
+                turbulence1->correct();
+                turbulence2->correct();
+            }
         }
 
         #include "write.H"

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/createFields.H

@@ -1,64 +1,38 @@
-    Info<< "Reading transportProperties\n" << endl;
+    Info<< "Creating twoPhaseSystem\n" << endl;
 
-    IOdictionary transportProperties
+    twoPhaseSystem fluid(mesh);
+
+    phaseModel& phase1 = fluid.phase1();
+    phaseModel& phase2 = fluid.phase2();
+
+    volScalarField& alpha1 = phase1;
+    volScalarField& alpha2 = phase2;
+
+    volVectorField& U1 = phase1.U();
+    surfaceScalarField& phi1 = phase1.phi();
+    surfaceScalarField alphaPhi1
     (
-        IOobject
+        IOobject::groupName("alphaPhi", phase1.name()),
-        (
+        fvc::interpolate(alpha1)*phi1
-            "transportProperties",
-            runTime.constant(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE
-        )
     );
 
-    word phase1Name
+    volVectorField& U2 = phase2.U();
+    surfaceScalarField& phi2 = phase2.phi();
+    surfaceScalarField alphaPhi2
     (
-        transportProperties.found("phases")
+        IOobject::groupName("alphaPhi", phase2.name()),
-      ? wordList(transportProperties.lookup("phases"))[0]
+        fvc::interpolate(alpha2)*phi2
-      : "1"
     );
 
-    word phase2Name
-    (
-        transportProperties.found("phases")
-      ? wordList(transportProperties.lookup("phases"))[1]
-      : "2"
-    );
-
-    autoPtr<phaseModel> phase1 = phaseModel::New
-    (
-        mesh,
-        transportProperties,
-        phase1Name
-    );
-
-    autoPtr<phaseModel> phase2 = phaseModel::New
-    (
-        mesh,
-        transportProperties,
-        phase2Name
-    );
-
-    volScalarField& alpha1 = phase1();
-    volScalarField& alpha2 = phase2();
-    alpha2 = scalar(1) - alpha1;
-
-    volVectorField& U1 = phase1->U();
-    surfaceScalarField& phi1 = phase1->phi();
-
-    volVectorField& U2 = phase2->U();
-    surfaceScalarField& phi2 = phase2->phi();
-
     dimensionedScalar pMin
     (
         "pMin",
         dimPressure,
-        transportProperties.lookup("pMin")
+        fluid.lookup("pMin")
     );
 
-    rhoThermo& thermo1 = phase1->thermo();
+    rhoThermo& thermo1 = phase1.thermo();
-    rhoThermo& thermo2 = phase2->thermo();
+    rhoThermo& thermo2 = phase2.thermo();
 
     volScalarField& p = thermo1.p();
 
@@ -78,7 +52,7 @@
             IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
-        alpha1*U1 + alpha2*U2
+        fluid.U()
     );
 
     surfaceScalarField phi
@@ -91,7 +65,7 @@
             IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
-        fvc::interpolate(alpha1)*phi1 + fvc::interpolate(alpha2)*phi2
+        fluid.phi()
     );
 
     volScalarField rho
@@ -104,7 +78,7 @@
             IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
-        alpha1*rho1 + alpha2*rho2
+        fluid.rho()
     );
 
     Info<< "Calculating field DDtU1 and DDtU2\n" << endl;
@@ -127,119 +101,11 @@
     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
-        (
-            "interfacialProperties",
-            runTime.constant(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE
-        )
-    );
-
-    autoPtr<dragModel> drag1 = dragModel::New
-    (
-        interfacialProperties,
-        alpha1,
-        phase1,
-        phase2
-    );
-
-    autoPtr<dragModel> drag2 = dragModel::New
-    (
-        interfacialProperties,
-        alpha2,
-        phase2,
-        phase1
-    );
-
-    autoPtr<heatTransferModel> heatTransfer1 = heatTransferModel::New
-    (
-        interfacialProperties,
-        alpha1,
-        phase1,
-        phase2
-    );
-
-    autoPtr<heatTransferModel> heatTransfer2 = heatTransferModel::New
-    (
-        interfacialProperties,
-        alpha2,
-        phase2,
-        phase1
-    );
-
-    word dispersedPhase(interfacialProperties.lookup("dispersedPhase"));
-
-    if
-    (
-        !(
-            dispersedPhase == phase1Name
-         || dispersedPhase == phase2Name
-         || dispersedPhase == "both"
-        )
-    )
-    {
-        FatalErrorIn(args.executable())
-            << "invalid dispersedPhase " << dispersedPhase
-            << exit(FatalError);
-    }
-
-    Info << "dispersedPhase is " << dispersedPhase << endl;
-
-    scalar residualPhaseFraction
-    (
-        readScalar
-        (
-            interfacialProperties.lookup("residualPhaseFraction")
-        )
-    );
-
-    dimensionedScalar residualSlip
-    (
-        "residualSlip",
-        dimVelocity,
-        interfacialProperties.lookup("residualSlip")
-    );
-
-    kineticTheoryModel kineticTheory
-    (
-        phase1,
-        U2,
-        alpha1,
-        drag1
-    );
-
     volScalarField rAU1
     (
         IOobject
         (
-            "rAU" + phase1Name,
+            IOobject::groupName("rAU", phase1.name()),
             runTime.timeName(),
             mesh,
             IOobject::NO_READ,
@@ -249,21 +115,20 @@
         dimensionedScalar("zero", dimensionSet(0, 0, 1, 0, 0), 0.0)
     );
 
-    surfaceScalarField ppMagf
+    volScalarField rAU2
     (
         IOobject
         (
-            "ppMagf",
+            IOobject::groupName("rAU", phase2.name()),
             runTime.timeName(),
             mesh,
             IOobject::NO_READ,
             IOobject::NO_WRITE
         ),
         mesh,
-        dimensionedScalar("zero", dimensionSet(0, 2, -1, 0, 0), 0.0)
+        dimensionedScalar("zero", dimensionSet(0, 0, 1, 0, 0), 0.0)
     );
 
-
     label pRefCell = 0;
     scalar pRefValue = 0.0;
     setRefCell(p, mesh.solutionDict().subDict("PIMPLE"), pRefCell, pRefValue);
@@ -290,5 +155,31 @@
 
 
     Info<< "Creating field kinetic energy K\n" << endl;
-    volScalarField K1("K" + phase1Name, 0.5*magSqr(U1));
+    volScalarField K1(IOobject::groupName("K", phase1.name()), 0.5*magSqr(U1));
-    volScalarField K2("K" + phase2Name, 0.5*magSqr(U2));
+    volScalarField K2(IOobject::groupName("K", phase2.name()), 0.5*magSqr(U2));
+
+    autoPtr<PhaseIncompressibleTurbulenceModel<phaseModel> >
+    turbulence1
+    (
+        PhaseIncompressibleTurbulenceModel<phaseModel>::New
+        (
+            alpha1,
+            U1,
+            alphaPhi1,
+            phi1,
+            phase1
+        )
+    );
+
+    autoPtr<PhaseIncompressibleTurbulenceModel<phaseModel> >
+    turbulence2
+    (
+        PhaseIncompressibleTurbulenceModel<phaseModel>::New
+        (
+            alpha2,
+            U2,
+            alphaPhi2,
+            phi2,
+            phase2
+        )
+    );

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/createRASTurbulence.H

@@ -1,189 +0,0 @@
-    IOdictionary RASProperties
-    (
-        IOobject
-        (
-            "RASProperties",
-            runTime.constant(),
-            mesh,
-            //IOobject::MUST_READ_IF_MODIFIED,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE
-        )
-    );
-
-
-    Switch turbulence
-    (
-        RASProperties.lookup("turbulence")
-    );
-
-    dictionary kEpsilonDict
-    (
-        RASProperties.subDictPtr("kEpsilonCoeffs")
-    );
-
-    dimensionedScalar Cmu
-    (
-        dimensionedScalar::lookupOrAddToDict
-        (
-            "Cmu",
-            kEpsilonDict,
-            0.09
-        )
-    );
-
-    dimensionedScalar C1
-    (
-        dimensionedScalar::lookupOrAddToDict
-        (
-            "C1",
-            kEpsilonDict,
-            1.44
-        )
-    );
-
-    dimensionedScalar C2
-    (
-        dimensionedScalar::lookupOrAddToDict
-        (
-            "C2",
-            kEpsilonDict,
-            1.92
-        )
-    );
-
-    dimensionedScalar alpha1k
-    (
-        dimensionedScalar::lookupOrAddToDict
-        (
-            "alpha1k",
-            kEpsilonDict,
-            1.0
-        )
-    );
-
-    dimensionedScalar alpha1Eps
-    (
-        dimensionedScalar::lookupOrAddToDict
-        (
-            "alpha1Eps",
-            kEpsilonDict,
-            0.76923
-        )
-    );
-
-    dimensionedScalar Prt
-    (
-        dimensioned<scalar>::lookupOrAddToDict
-        (
-            "Prt",
-            kEpsilonDict,
-            1.0
-        )
-    );
-
-    dictionary wallFunctionDict
-    (
-        RASProperties.subDictPtr("wallFunctionCoeffs")
-    );
-
-    dimensionedScalar kappa
-    (
-        dimensionedScalar::lookupOrAddToDict
-        (
-            "kappa",
-            wallFunctionDict,
-            0.41
-        )
-    );
-
-    dimensionedScalar E
-    (
-        dimensionedScalar::lookupOrAddToDict
-        (
-            "E",
-            wallFunctionDict,
-            9.8
-        )
-    );
-
-    if (RASProperties.lookupOrDefault("printCoeffs", false))
-    {
-        Info<< "kEpsilonCoeffs" << kEpsilonDict << nl
-            << "wallFunctionCoeffs" << wallFunctionDict << endl;
-    }
-
-
-    nearWallDist y(mesh);
-
-
-    Info<< "Reading field k\n" << endl;
-    volScalarField k
-    (
-        IOobject
-        (
-            "k",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-    Info<< "Reading field epsilon\n" << endl;
-    volScalarField epsilon
-    (
-        IOobject
-        (
-            "epsilon",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-
-    Info<< "Calculating field nut2\n" << endl;
-    volScalarField nut2
-    (
-        IOobject
-        (
-            "nut" + phase2Name,
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        Cmu*sqr(k)/epsilon
-    );
-
-    Info<< "Calculating field nuEff1\n" << endl;
-    volScalarField nuEff1
-    (
-        IOobject
-        (
-            "nuEff" + phase1Name,
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        sqr(Ct)*nut2 + thermo1.mu()/rho1
-    );
-
-    Info<< "Calculating field nuEff2\n" << endl;
-    volScalarField nuEff2
-    (
-        IOobject
-        (
-            "nuEff" + phase2Name,
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        nut2 + thermo2.mu()/rho2
-    );

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialCoeffs.H

@@ -1,88 +0,0 @@
-volScalarField dragCoeff
-(
-    IOobject
-    (
-        "dragCoeff",
-        runTime.timeName(),
-        mesh
-    ),
-    mesh,
-    dimensionedScalar("dragCoeff", dimensionSet(1, -3, -1, 0, 0), 0)
-);
-
-volVectorField liftForce
-(
-    IOobject
-    (
-        "liftForce",
-        runTime.timeName(),
-        mesh
-    ),
-    mesh,
-    dimensionedVector("liftForce", dimensionSet(1, -2, -2, 0, 0), vector::zero)
-);
-
-volScalarField heatTransferCoeff
-(
-    IOobject
-    (
-        "heatTransferCoeff",
-        runTime.timeName(),
-        mesh
-    ),
-    mesh,
-    dimensionedScalar("heatTransferCoeff", dimensionSet(1, -1, -3, -1, 0), 0)
-);
-
-{
-    volVectorField Ur(U1 - U2);
-    volScalarField magUr(mag(Ur) + residualSlip);
-
-    if (dispersedPhase == phase1Name)
-    {
-        dragCoeff = drag1->K(magUr);
-        heatTransferCoeff = heatTransfer1->K(magUr);
-    }
-    else if (dispersedPhase == phase2Name)
-    {
-        dragCoeff = drag2->K(magUr);
-        heatTransferCoeff = heatTransfer2->K(magUr);
-    }
-    else if (dispersedPhase == "both")
-    {
-        dragCoeff =
-        (
-            alpha2*drag1->K(magUr)
-          + alpha1*drag2->K(magUr)
-        );
-
-        heatTransferCoeff =
-        (
-            alpha2*heatTransfer1->K(magUr)
-          + alpha1*heatTransfer2->K(magUr)
-        );
-    }
-    else
-    {
-        FatalErrorIn(args.executable())
-            << "dispersedPhase: " << dispersedPhase << " is incorrect"
-            << exit(FatalError);
-    }
-
-    volScalarField alphaCoeff(max(alpha1*alpha2, residualPhaseFraction));
-    dragCoeff *= alphaCoeff;
-    heatTransferCoeff *= alphaCoeff;
-
-    liftForce = Cl*(alpha1*rho1 + alpha2*rho2)*(Ur ^ fvc::curl(U));
-
-    // Remove lift, drag and phase heat-transfer at fixed-flux boundaries
-    forAll(phi1.boundaryField(), patchi)
-    {
-        if (isA<fixedValueFvsPatchScalarField>(phi1.boundaryField()[patchi]))
-        {
-            dragCoeff.boundaryField()[patchi] = 0.0;
-            heatTransferCoeff.boundaryField()[patchi] = 0.0;
-            liftForce.boundaryField()[patchi] = vector::zero;
-        }
-    }
-}

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

@@ -1,9 +1,10 @@
 EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/transportModels/incompressible/transportModel \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I../phaseModel/lnInclude
+    -I../twoPhaseSystem/lnInclude
 
 LIB_LIBS = \
-    -lcompressiblePhaseModel \
+    -lcompressibleTwoPhaseSystem \
     -lfluidThermophysicalModels \
     -lspecie

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/dragModels/dragModel/newDragModel.C

@@ -37,7 +37,7 @@ Foam::autoPtr<Foam::dragModel> Foam::dragModel::New
 {
     word dragModelType
     (
-        interfaceDict.lookup("dragModel" + phase1.name())
+        interfaceDict.lookup(phase1.name())
     );
 
     Info << "Selecting dragModel for phase "

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/interfacialModels/heatTransferModels/heatTransferModel/newHeatTransferModel.C

@@ -37,7 +37,7 @@ Foam::autoPtr<Foam::heatTransferModel> Foam::heatTransferModel::New
 {
     word heatTransferModelType
     (
-        interfaceDict.lookup("heatTransferModel" + phase1.name())
+        interfaceDict.lookup(phase1.name())
     );
 
     Info<< "Selecting heatTransferModel for phase "

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/kineticTheoryModels/Make/files

@@ -1,32 +0,0 @@
-kineticTheoryModel/kineticTheoryModel.C
-
-viscosityModel/viscosityModel/viscosityModel.C
-viscosityModel/viscosityModel/newViscosityModel.C
-viscosityModel/Gidaspow/GidaspowViscosity.C
-viscosityModel/Syamlal/SyamlalViscosity.C
-viscosityModel/HrenyaSinclair/HrenyaSinclairViscosity.C
-viscosityModel/none/noneViscosity.C
-
-conductivityModel/conductivityModel/conductivityModel.C
-conductivityModel/conductivityModel/newConductivityModel.C
-conductivityModel/Gidaspow/GidaspowConductivity.C
-conductivityModel/Syamlal/SyamlalConductivity.C
-conductivityModel/HrenyaSinclair/HrenyaSinclairConductivity.C
-
-radialModel/radialModel/radialModel.C
-radialModel/radialModel/newRadialModel.C
-radialModel/CarnahanStarling/CarnahanStarlingRadial.C
-radialModel/LunSavage/LunSavageRadial.C
-radialModel/SinclairJackson/SinclairJacksonRadial.C
-
-granularPressureModel/granularPressureModel/granularPressureModel.C
-granularPressureModel/granularPressureModel/newGranularPressureModel.C
-granularPressureModel/Lun/LunPressure.C
-granularPressureModel/SyamlalRogersOBrien/SyamlalRogersOBrienPressure.C
-
-frictionalStressModel/frictionalStressModel/frictionalStressModel.C
-frictionalStressModel/frictionalStressModel/newFrictionalStressModel.C
-frictionalStressModel/JohnsonJackson/JohnsonJacksonFrictionalStress.C
-frictionalStressModel/Schaeffer/SchaefferFrictionalStress.C
-
-LIB = $(FOAM_LIBBIN)/libcompressibleKineticTheoryModel

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

@@ -1,6 +0,0 @@
-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

@@ -1,389 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "kineticTheoryModel.H"
-#include "surfaceInterpolate.H"
-#include "mathematicalConstants.H"
-#include "fvCFD.H"
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::kineticTheoryModel::kineticTheoryModel
-(
-    const Foam::phaseModel& phase1,
-    const Foam::volVectorField& U2,
-    const Foam::volScalarField& alpha1,
-    const Foam::dragModel& draga
-)
-:
-    phase1_(phase1),
-    U1_(phase1.U()),
-    U2_(U2),
-    alpha1_(alpha1),
-    phi1_(phase1.phi()),
-    draga_(draga),
-
-    rho1_(phase1.rho()),
-
-    kineticTheoryProperties_
-    (
-        IOobject
-        (
-            "kineticTheoryProperties",
-            U1_.time().constant(),
-            U1_.mesh(),
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE
-        )
-    ),
-    kineticTheory_(kineticTheoryProperties_.lookup("kineticTheory")),
-    equilibrium_(kineticTheoryProperties_.lookup("equilibrium")),
-
-    viscosityModel_
-    (
-        kineticTheoryModels::viscosityModel::New
-        (
-            kineticTheoryProperties_
-        )
-    ),
-    conductivityModel_
-    (
-        kineticTheoryModels::conductivityModel::New
-        (
-            kineticTheoryProperties_
-        )
-    ),
-    radialModel_
-    (
-        kineticTheoryModels::radialModel::New
-        (
-            kineticTheoryProperties_
-        )
-    ),
-    granularPressureModel_
-    (
-        kineticTheoryModels::granularPressureModel::New
-        (
-            kineticTheoryProperties_
-        )
-    ),
-    frictionalStressModel_
-    (
-        kineticTheoryModels::frictionalStressModel::New
-        (
-            kineticTheoryProperties_
-        )
-    ),
-    e_(kineticTheoryProperties_.lookup("e")),
-    alphaMax_(kineticTheoryProperties_.lookup("alphaMax")),
-    alphaMinFriction_(kineticTheoryProperties_.lookup("alphaMinFriction")),
-    Fr_(kineticTheoryProperties_.lookup("Fr")),
-    eta_(kineticTheoryProperties_.lookup("eta")),
-    p_(kineticTheoryProperties_.lookup("p")),
-    phi_(dimensionedScalar(kineticTheoryProperties_.lookup("phi"))*M_PI/180.0),
-    Theta_
-    (
-        IOobject
-        (
-            "Theta",
-            U1_.time().timeName(),
-            U1_.mesh(),
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        U1_.mesh()
-    ),
-    mu1_
-    (
-        IOobject
-        (
-            "mu" + phase1.name(),
-            U1_.time().timeName(),
-            U1_.mesh(),
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        U1_.mesh(),
-        dimensionedScalar("zero", dimensionSet(1, -1, -1, 0, 0), 0.0)
-    ),
-    lambda_
-    (
-        IOobject
-        (
-            "lambda",
-            U1_.time().timeName(),
-            U1_.mesh(),
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        U1_.mesh(),
-        dimensionedScalar("zero", dimensionSet(1, -1, -1, 0, 0), 0.0)
-    ),
-    pa_
-    (
-        IOobject
-        (
-            "pa",
-            U1_.time().timeName(),
-            U1_.mesh(),
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        U1_.mesh(),
-        dimensionedScalar("zero", dimensionSet(1, -1, -2, 0, 0), 0.0)
-    ),
-    kappa_
-    (
-        IOobject
-        (
-            "kappa",
-            U1_.time().timeName(),
-            U1_.mesh(),
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        U1_.mesh(),
-        dimensionedScalar("zero", dimensionSet(1, -1, -1, 0, 0), 0.0)
-    ),
-    gs0_
-    (
-        IOobject
-        (
-            "gs0",
-            U1_.time().timeName(),
-            U1_.mesh(),
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        U1_.mesh(),
-        dimensionedScalar("zero", dimensionSet(0, 0, 0, 0, 0), 1.0)
-    )
-{}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::kineticTheoryModel::~kineticTheoryModel()
-{}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-void Foam::kineticTheoryModel::solve(const volTensorField& gradU1t)
-{
-    if (!kineticTheory_)
-    {
-        return;
-    }
-
-    const scalar sqrtPi = sqrt(constant::mathematical::pi);
-
-    volScalarField da_(phase1_.d());
-
-    surfaceScalarField phi(1.5*phi1_*fvc::interpolate(rho1_*alpha1_));
-
-    volTensorField dU(gradU1t.T());
-    volSymmTensorField D(symm(dU));
-
-    // NB, drag = K*alpha1*alpha2,
-    // (the alpha1 and alpha2 has been extracted from the drag function for
-    // numerical reasons)
-    volScalarField Ur(mag(U1_ - U2_));
-    volScalarField alpha2Prim(alpha1_*(1.0 - alpha1_)*draga_.K(Ur));
-
-    // Calculating the radial distribution function (solid volume fraction is
-    //  limited close to the packing limit, but this needs improvements)
-    //  The solution is higly unstable close to the packing limit.
-    gs0_ = radialModel_->g0
-    (
-        min(max(alpha1_, scalar(1e-6)), alphaMax_ - 0.01),
-        alphaMax_
-    );
-
-    // particle pressure - coefficient in front of Theta (Eq. 3.22, p. 45)
-    volScalarField PsCoeff
-    (
-        granularPressureModel_->granularPressureCoeff
-        (
-            alpha1_,
-            gs0_,
-            rho1_,
-            e_
-        )
-    );
-
-    // 'thermal' conductivity (Table 3.3, p. 49)
-    kappa_ = conductivityModel_->kappa(alpha1_, Theta_, gs0_, rho1_, da_, e_);
-
-    // particle viscosity (Table 3.2, p.47)
-    mu1_ = viscosityModel_->mu1(alpha1_, Theta_, gs0_, rho1_, da_, e_);
-
-    dimensionedScalar Tsmall
-    (
-        "small",
-        dimensionSet(0 , 2 ,-2 ,0 , 0, 0, 0),
-        1.0e-6
-    );
-
-    dimensionedScalar TsmallSqrt = sqrt(Tsmall);
-    volScalarField ThetaSqrt(sqrt(Theta_));
-
-    // dissipation (Eq. 3.24, p.50)
-    volScalarField gammaCoeff
-    (
-        12.0*(1.0 - sqr(e_))*sqr(alpha1_)*rho1_*gs0_*(1.0/da_)*ThetaSqrt/sqrtPi
-    );
-
-    // Eq. 3.25, p. 50 Js = J1 - J2
-    volScalarField J1(3.0*alpha2Prim);
-    volScalarField J2
-    (
-        0.25*sqr(alpha2Prim)*da_*sqr(Ur)
-       /(max(alpha1_, scalar(1e-6))*rho1_*sqrtPi*(ThetaSqrt + TsmallSqrt))
-    );
-
-    // bulk viscosity  p. 45 (Lun et al. 1984).
-    lambda_ = (4.0/3.0)*sqr(alpha1_)*rho1_*da_*gs0_*(1.0+e_)*ThetaSqrt/sqrtPi;
-
-    // stress tensor, Definitions, Table 3.1, p. 43
-    volSymmTensorField tau(2.0*mu1_*D + (lambda_ - (2.0/3.0)*mu1_)*tr(D)*I);
-
-    if (!equilibrium_)
-    {
-        // construct the granular temperature equation (Eq. 3.20, p. 44)
-        // NB. note that there are two typos in Eq. 3.20
-        // no grad infront of Ps
-        // wrong sign infront of laplacian
-        fvScalarMatrix ThetaEqn
-        (
-            fvm::ddt(1.5*alpha1_*rho1_, Theta_)
-          + fvm::div(phi, Theta_, "div(phi,Theta)")
-         ==
-            fvm::SuSp(-((PsCoeff*I) && dU), Theta_)
-          + (tau && dU)
-          + fvm::laplacian(kappa_, Theta_, "laplacian(kappa,Theta)")
-          + fvm::Sp(-gammaCoeff, Theta_)
-          + fvm::Sp(-J1, Theta_)
-          + fvm::Sp(J2/(Theta_ + Tsmall), Theta_)
-        );
-
-        ThetaEqn.relax();
-        ThetaEqn.solve();
-    }
-    else
-    {
-        // equilibrium => dissipation == production
-        // Eq. 4.14, p.82
-        volScalarField K1(2.0*(1.0 + e_)*rho1_*gs0_);
-        volScalarField K3
-        (
-            0.5*da_*rho1_*
-            (
-                (sqrtPi/(3.0*(3.0-e_)))
-               *(1.0 + 0.4*(1.0 + e_)*(3.0*e_ - 1.0)*alpha1_*gs0_)
-               +1.6*alpha1_*gs0_*(1.0 + e_)/sqrtPi
-            )
-        );
-
-        volScalarField K2
-        (
-            4.0*da_*rho1_*(1.0 + e_)*alpha1_*gs0_/(3.0*sqrtPi) - 2.0*K3/3.0
-        );
-
-        volScalarField K4(12.0*(1.0 - sqr(e_))*rho1_*gs0_/(da_*sqrtPi));
-
-        volScalarField trD(tr(D));
-        volScalarField tr2D(sqr(trD));
-        volScalarField trD2(tr(D & D));
-
-        volScalarField t1(K1*alpha1_ + rho1_);
-        volScalarField l1(-t1*trD);
-        volScalarField l2(sqr(t1)*tr2D);
-        volScalarField l3
-        (
-            4.0
-           *K4
-           *max(alpha1_, scalar(1e-6))
-           *(2.0*K3*trD2 + K2*tr2D)
-        );
-
-        Theta_ = sqr((l1 + sqrt(l2 + l3))/(2.0*(alpha1_ + 1.0e-4)*K4));
-    }
-
-    Theta_.max(1.0e-15);
-    Theta_.min(1.0e+3);
-
-    volScalarField pf
-    (
-        frictionalStressModel_->frictionalPressure
-        (
-            alpha1_,
-            alphaMinFriction_,
-            alphaMax_,
-            Fr_,
-            eta_,
-            p_
-        )
-    );
-
-    PsCoeff += pf/(Theta_+Tsmall);
-
-    PsCoeff.min(1.0e+10);
-    PsCoeff.max(-1.0e+10);
-
-    // update particle pressure
-    pa_ = PsCoeff*Theta_;
-
-    // frictional shear stress, Eq. 3.30, p. 52
-    volScalarField muf
-    (
-        frictionalStressModel_->muf
-        (
-            alpha1_,
-            alphaMax_,
-            pf,
-            D,
-            phi_
-        )
-    );
-
-    // add frictional stress
-    mu1_ += muf;
-    mu1_.min(1.0e+2);
-    mu1_.max(0.0);
-
-    Info<< "kinTheory: max(Theta) = " << max(Theta_).value() << endl;
-
-    volScalarField ktn(mu1_/rho1_);
-
-    Info<< "kinTheory: min(nu1) = " << min(ktn).value()
-        << ", max(nu1) = " << max(ktn).value() << endl;
-
-    Info<< "kinTheory: min(pa) = " << min(pa_).value()
-        << ", max(pa) = " << max(pa_).value() << endl;
-}
-
-
-// ************************************************************************* //

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

@@ -1,196 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
-
-Class
-    Foam::kineticTheoryModel
-
-Description
-
-SourceFiles
-    kineticTheoryModel.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef kineticTheoryModel_H
-#define kineticTheoryModel_H
-
-#include "dragModel.H"
-#include "phaseModel.H"
-#include "autoPtr.H"
-#include "viscosityModel.H"
-#include "conductivityModel.H"
-#include "radialModel.H"
-#include "granularPressureModel.H"
-#include "frictionalStressModel.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-                           Class kineticTheoryModel Declaration
-\*---------------------------------------------------------------------------*/
-
-class kineticTheoryModel
-{
-    // Private data
-
-        const phaseModel& phase1_;
-        const volVectorField& U1_;
-        const volVectorField& U2_;
-        const volScalarField& alpha1_;
-        const surfaceScalarField& phi1_;
-
-        const dragModel& draga_;
-
-        const volScalarField& rho1_;
-
-        //- dictionary holding the modeling info
-        IOdictionary kineticTheoryProperties_;
-
-        //- use kinetic theory or not.
-        Switch kineticTheory_;
-
-        //- use generation == dissipation
-        Switch equilibrium_;
-
-        autoPtr<kineticTheoryModels::viscosityModel> viscosityModel_;
-
-        autoPtr<kineticTheoryModels::conductivityModel> conductivityModel_;
-
-        autoPtr<kineticTheoryModels::radialModel> radialModel_;
-
-        autoPtr<kineticTheoryModels::granularPressureModel>
-            granularPressureModel_;
-
-        autoPtr<kineticTheoryModels::frictionalStressModel>
-            frictionalStressModel_;
-
-        //- coefficient of restitution
-        const dimensionedScalar e_;
-
-        //- maximum packing
-        const dimensionedScalar alphaMax_;
-
-        //- min value for which the frictional stresses are zero
-        const dimensionedScalar alphaMinFriction_;
-
-        //- material constant for frictional normal stress
-        const dimensionedScalar Fr_;
-
-        //- material constant for frictional normal stress
-        const dimensionedScalar eta_;
-
-        //- material constant for frictional normal stress
-        const dimensionedScalar p_;
-
-        //- angle of internal friction
-        const dimensionedScalar phi_;
-
-        //- The granular energy/temperature
-        volScalarField Theta_;
-
-        //- The granular viscosity
-        volScalarField mu1_;
-
-        //- The granular bulk viscosity
-        volScalarField lambda_;
-
-        //- The granular pressure
-        volScalarField pa_;
-
-        //- The granular temperature conductivity
-        volScalarField kappa_;
-
-        //- The radial distribution function
-        volScalarField gs0_;
-
-
-    // Private Member Functions
-
-        //- Disallow default bitwise copy construct
-        kineticTheoryModel(const kineticTheoryModel&);
-
-        //- Disallow default bitwise assignment
-        void operator=(const kineticTheoryModel&);
-
-
-public:
-
-    // Constructors
-
-        //- Construct from components
-        kineticTheoryModel
-        (
-            const phaseModel& phase1,
-            const volVectorField& U2,
-            const volScalarField& alpha1,
-            const dragModel& draga
-        );
-
-
-    //- Destructor
-    virtual ~kineticTheoryModel();
-
-
-    // Member Functions
-
-        void solve(const volTensorField& gradU1t);
-
-        bool on() const
-        {
-            return kineticTheory_;
-        }
-
-        const volScalarField& mu1() const
-        {
-            return mu1_;
-        }
-
-        const volScalarField& pa() const
-        {
-            return pa_;
-        }
-
-        const volScalarField& lambda() const
-        {
-            return lambda_;
-        }
-
-        const volScalarField& kappa() const
-        {
-            return kappa_;
-        }
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H

@@ -3,15 +3,23 @@
     surfaceScalarField alpha2f(scalar(1) - alpha1f);
 
     rAU1 = 1.0/U1Eqn.A();
-    volScalarField rAU2(1.0/U2Eqn.A());
+    rAU2 = 1.0/U2Eqn.A();
 
     surfaceScalarField rAlphaAU1f(fvc::interpolate(alpha1*rAU1));
     surfaceScalarField rAlphaAU2f(fvc::interpolate(alpha2*rAU2));
 
-    volVectorField HbyA1("HbyA" + phase1Name, U1);
+    volVectorField HbyA1
+    (
+        IOobject::groupName("HbyA", phase1.name()),
+        U1
+    );
     HbyA1 = rAU1*U1Eqn.H();
 
-    volVectorField HbyA2("HbyA" + phase2Name, U2);
+    volVectorField HbyA2
+    (
+        IOobject::groupName("HbyA", phase2.name()),
+        U2
+    );
     HbyA2 = rAU2*U2Eqn.H();
 
     mrfZones.absoluteFlux(phi1.oldTime());
@@ -19,30 +27,32 @@
     mrfZones.absoluteFlux(phi2.oldTime());
     mrfZones.absoluteFlux(phi2);
 
-    surfaceScalarField ppDrag("ppDrag", 0.0*phi1);
+    // Phase-1 pressure flux (e.g. due to particle-particle pressure)
+    surfaceScalarField phiP1
+    (
+        "phiP1",
+        fvc::interpolate((1.0/rho1)*rAU1*turbulence1().pPrime())
+       *fvc::snGrad(alpha1)*mesh.magSf()
+    );
 
-    if (g0.value() > 0.0)
+    // Phase-2 pressure flux (e.g. due to particle-particle pressure)
-    {
+    surfaceScalarField phiP2
-        ppDrag -= ppMagf*fvc::snGrad(alpha1)*mesh.magSf();
+    (
-    }
+        "phiP2",
-
+        fvc::interpolate((1.0/rho2)*rAU2*turbulence2().pPrime())
-    if (kineticTheory.on())
+       *fvc::snGrad(alpha2)*mesh.magSf()
-    {
+    );
-        ppDrag -=
-            fvc::interpolate(1.0/rho1)*rAlphaAU1f
-           *fvc::snGrad(kineticTheory.pa())*mesh.magSf();
-    }
 
     surfaceScalarField phiHbyA1
     (
-        "phiHbyA" + phase1Name,
+        IOobject::groupName("phiHbyA", phase1.name()),
         (fvc::interpolate(HbyA1) & mesh.Sf())
       + fvc::ddtPhiCorr(rAU1, alpha1, U1, phi1)
     );
 
     surfaceScalarField phiHbyA2
     (
-        "phiHbyA" + phase2Name,
+        IOobject::groupName("phiHbyA", phase2.name()),
         (fvc::interpolate(HbyA2) & mesh.Sf())
       + fvc::ddtPhiCorr(rAU2, alpha2, U2, phi2)
     );
@@ -53,7 +63,7 @@
     phiHbyA1 +=
     (
         fvc::interpolate((1.0/rho1)*rAU1*dragCoeff)*phi2
-      + ppDrag
+      - phiP1
       + rAlphaAU1f*(g & mesh.Sf())
     );
     mrfZones.relativeFlux(phiHbyA1);
@@ -61,6 +71,7 @@
     phiHbyA2 +=
     (
         fvc::interpolate((1.0/rho2)*rAU2*dragCoeff)*phi1
+      - phiP2
       + rAlphaAU2f*(g & mesh.Sf())
     );
     mrfZones.relativeFlux(phiHbyA2);
@@ -92,12 +103,12 @@
     {
         surfaceScalarField phid1
         (
-            "phid" + phase1Name,
+            IOobject::groupName("phid", phase1.name()),
             fvc::interpolate(psi1)*phi1
         );
         surfaceScalarField phid2
         (
-            "phid" + phase2Name,
+            IOobject::groupName("phid", phase2.name()),
             fvc::interpolate(psi2)*phi2
         );
 
@@ -173,29 +184,30 @@
             mSfGradp = pEqnIncomp.flux()/Dp;
 
             U1 = HbyA1
-               + fvc::reconstruct
+              + fvc::reconstruct
-                 (
+                (
-                     ppDrag
+                    rAlphaAU1f
-                   + rAlphaAU1f
+                   *(
-                    *(
+                        (g & mesh.Sf())
-                         (g & mesh.Sf())
+                      + mSfGradp/fvc::interpolate(rho1)
-                       + mSfGradp/fvc::interpolate(rho1)
+                    )
-                     )
+                  - phiP1
-                 );
+                );
             U1.correctBoundaryConditions();
 
             U2 = HbyA2
-               + fvc::reconstruct
+              + fvc::reconstruct
-                 (
+                (
-                     rAlphaAU2f
+                    rAlphaAU2f
-                    *(
+                   *(
                         (g & mesh.Sf())
                       + mSfGradp/fvc::interpolate(rho2)
                     )
-                 );
+                  - phiP2
+                );
             U2.correctBoundaryConditions();
 
-            U = alpha1*U1 + alpha2*U2;
+            U = fluid.U();
         }
     }
 

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/packingLimiter.H

@@ -1,36 +0,0 @@
-    if (packingLimiter)
-    {
-        // Calculating exceeding volume fractions
-        volScalarField alpha1Ex(max(alpha1 - alphaMax, scalar(0)));
-
-        // Finding neighbouring cells of the whole domain
-        labelListList neighbour = mesh.cellCells();
-        scalarField cellVolumes(mesh.cellVolumes());
-
-        forAll (alpha1Ex, celli)
-        {
-            // Finding the labels of the neighbouring cells
-            labelList neighbourCell = neighbour[celli];
-
-            // Initializing neighbouring cells contribution
-            scalar neighboursEx = 0.0;
-
-            forAll (neighbourCell, cellj)
-            {
-                labelList neighboursNeighbour = neighbour[neighbourCell[cellj]];
-                scalar neighboursNeighbourCellVolumes = 0.0;
-
-                forAll (neighboursNeighbour, cellk)
-                {
-                    neighboursNeighbourCellVolumes +=
-                        cellVolumes[neighboursNeighbour[cellk]];
-                }
-
-                neighboursEx +=
-                    alpha1Ex[neighbourCell[cellj]]*cellVolumes[celli]
-                   /neighboursNeighbourCellVolumes;
-            }
-
-            alpha1[celli] += neighboursEx - alpha1Ex[celli];
-        }
-    }

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/Make/files

@@ -0,0 +1,35 @@
+phaseIncompressibleTurbulenceModels.C
+phasePressureModel/phasePressureModel.C
+
+kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.C
+
+kineticTheoryModels/viscosityModel/viscosityModel/viscosityModel.C
+kineticTheoryModels/viscosityModel/viscosityModel/newViscosityModel.C
+kineticTheoryModels/viscosityModel/Gidaspow/GidaspowViscosity.C
+kineticTheoryModels/viscosityModel/Syamlal/SyamlalViscosity.C
+kineticTheoryModels/viscosityModel/HrenyaSinclair/HrenyaSinclairViscosity.C
+kineticTheoryModels/viscosityModel/none/noneViscosity.C
+
+kineticTheoryModels/conductivityModel/conductivityModel/conductivityModel.C
+kineticTheoryModels/conductivityModel/conductivityModel/newConductivityModel.C
+kineticTheoryModels/conductivityModel/Gidaspow/GidaspowConductivity.C
+kineticTheoryModels/conductivityModel/Syamlal/SyamlalConductivity.C
+kineticTheoryModels/conductivityModel/HrenyaSinclair/HrenyaSinclairConductivity.C
+
+kineticTheoryModels/radialModel/radialModel/radialModel.C
+kineticTheoryModels/radialModel/radialModel/newRadialModel.C
+kineticTheoryModels/radialModel/CarnahanStarling/CarnahanStarlingRadial.C
+kineticTheoryModels/radialModel/LunSavage/LunSavageRadial.C
+kineticTheoryModels/radialModel/SinclairJackson/SinclairJacksonRadial.C
+
+kineticTheoryModels/granularPressureModel/granularPressureModel/granularPressureModel.C
+kineticTheoryModels/granularPressureModel/granularPressureModel/newGranularPressureModel.C
+kineticTheoryModels/granularPressureModel/Lun/LunPressure.C
+kineticTheoryModels/granularPressureModel/SyamlalRogersOBrien/SyamlalRogersOBrienPressure.C
+
+kineticTheoryModels/frictionalStressModel/frictionalStressModel/frictionalStressModel.C
+kineticTheoryModels/frictionalStressModel/frictionalStressModel/newFrictionalStressModel.C
+kineticTheoryModels/frictionalStressModel/JohnsonJackson/JohnsonJacksonFrictionalStress.C
+kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.C
+
+LIB = $(FOAM_LIBBIN)/libphaseIncompressibleTurbulenceModels

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/Make/options

@@ -0,0 +1,10 @@
+EXE_INC = \
+    -I$(LIB_SRC)/foam/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/transportModels/incompressible/transportModel \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModel/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/phaseIncompressible/lnInclude \
+    -I../twoPhaseSystem/lnInclude \
+    -I../interfacialModels/lnInclude

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/kineticTheoryModels/conductivityModel/Gidaspow/GidaspowConductivity.C

@@ -80,7 +80,7 @@ Foam::kineticTheoryModels::conductivityModels::Gidaspow::kappa
 {
     const scalar sqrtPi = sqrt(constant::mathematical::pi);
 
-    return rho1*da*sqrt(Theta)*
+    return da*sqrt(Theta)*
     (
         2.0*sqr(alpha1)*g0*(1.0 + e)/sqrtPi
       + (9.0/8.0)*sqrtPi*g0*0.5*(1.0 + e)*sqr(alpha1)

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/kineticTheoryModels/conductivityModel/HrenyaSinclair/HrenyaSinclairConductivity.C

@@ -56,8 +56,8 @@ Foam::kineticTheoryModels::conductivityModels::HrenyaSinclair::HrenyaSinclair
 )
 :
     conductivityModel(dict),
-    coeffsDict_(dict.subDict(typeName + "Coeffs")),
+    coeffDict_(dict.subDict(typeName + "Coeffs")),
-    L_(coeffsDict_.lookup("L"))
+    L_("L", dimensionSet(0, 1, 0, 0, 0), coeffDict_.lookup("L"))
 {}
 
 
@@ -88,7 +88,7 @@ Foam::kineticTheoryModels::conductivityModels::HrenyaSinclair::kappa
         scalar(1) + da/(6.0*sqrt(2.0)*(alpha1 + scalar(1.0e-5)))/L_
     );
 
-    return rho1*da*sqrt(Theta)*
+    return da*sqrt(Theta)*
     (
         2.0*sqr(alpha1)*g0*(1.0 + e)/sqrtPi
       + (9.0/8.0)*sqrtPi*0.25*sqr(1.0 + e)*(2.0*e - 1.0)*sqr(alpha1)
@@ -101,4 +101,14 @@ Foam::kineticTheoryModels::conductivityModels::HrenyaSinclair::kappa
 }
 
 
+bool Foam::kineticTheoryModels::conductivityModels::HrenyaSinclair::read()
+{
+    coeffDict_ <<= dict_.subDict(typeName + "Coeffs");
+
+    L_.readIfPresent(coeffDict_);
+
+    return true;
+}
+
+
 // ************************************************************************* //

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/kineticTheoryModels/conductivityModel/HrenyaSinclair/HrenyaSinclairConductivity.H

@@ -53,12 +53,12 @@ class HrenyaSinclair
 :
     public conductivityModel
 {
-
+        dictionary coeffDict_;
-        dictionary coeffsDict_;
 
         //- characteristic length of geometry
         dimensionedScalar L_;
 
+
 public:
 
         //- Runtime type information
@@ -86,6 +86,8 @@ public:
             const volScalarField& da,
             const dimensionedScalar& e
         ) const;
+
+        virtual bool read();
 };
 
 

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/kineticTheoryModels/conductivityModel/Syamlal/SyamlalConductivity.C

@@ -80,7 +80,7 @@ Foam::kineticTheoryModels::conductivityModels::Syamlal::kappa
 {
     const scalar sqrtPi = sqrt(constant::mathematical::pi);
 
-    return rho1*da*sqrt(Theta)*
+    return da*sqrt(Theta)*
     (
         2.0*sqr(alpha1)*g0*(1.0 + e)/sqrtPi
       + (9.0/8.0)*sqrtPi*g0*0.25*sqr(1.0 + e)*(2.0*e - 1.0)*sqr(alpha1)

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/kineticTheoryModels/conductivityModel/conductivityModel/conductivityModel.H

@@ -113,6 +113,11 @@ public:
             const volScalarField& da,
             const dimensionedScalar& e
         ) const = 0;
+
+        virtual bool read()
+        {
+            return true;
+        }
 };
 
 

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/kineticTheoryModels/frictionalStressModel/JohnsonJackson/JohnsonJacksonFrictionalStress.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,6 +25,7 @@ License
 
 #include "JohnsonJacksonFrictionalStress.H"
 #include "addToRunTimeSelectionTable.H"
+#include "mathematicalConstants.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -55,8 +56,15 @@ JohnsonJackson
     const dictionary& dict
 )
 :
-    frictionalStressModel(dict)
+    frictionalStressModel(dict),
-{}
+    coeffDict_(dict.subDict(typeName + "Coeffs")),
+    Fr_("Fr", dimensionSet(1, -1, -2, 0, 0), coeffDict_.lookup("Fr")),
+    eta_("eta", dimless, coeffDict_.lookup("eta")),
+    p_("p", dimless, coeffDict_.lookup("p")),
+    phi_("phi", dimless, coeffDict_.lookup("phi"))
+{
+    phi_ *= constant::mathematical::pi/180.0;
+}
 
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
@@ -74,16 +82,13 @@ frictionalPressure
 (
     const volScalarField& alpha1,
     const dimensionedScalar& alphaMinFriction,
-    const dimensionedScalar& alphaMax,
+    const dimensionedScalar& alphaMax
-    const dimensionedScalar& Fr,
-    const dimensionedScalar& eta,
-    const dimensionedScalar& p
 ) const
 {
 
     return
-        Fr*pow(max(alpha1 - alphaMinFriction, scalar(0)), eta)
+        Fr_*pow(max(alpha1 - alphaMinFriction, scalar(0)), eta_)
-       /pow(max(alphaMax - alpha1, scalar(5.0e-2)), p);
+       /pow(max(alphaMax - alpha1, scalar(5.0e-2)), p_);
 }
 
 
@@ -93,32 +98,41 @@ frictionalPressurePrime
 (
     const volScalarField& alpha1,
     const dimensionedScalar& alphaMinFriction,
-    const dimensionedScalar& alphaMax,
+    const dimensionedScalar& alphaMax
-    const dimensionedScalar& Fr,
-    const dimensionedScalar& eta,
-    const dimensionedScalar& p
 ) const
 {
-    return Fr*
+    return Fr_*
     (
-        eta*pow(max(alpha1 - alphaMinFriction, scalar(0)), eta - 1.0)
+        eta_*pow(max(alpha1 - alphaMinFriction, scalar(0)), eta_ - 1.0)
        *(alphaMax-alpha1)
-      + p*pow(max(alpha1 - alphaMinFriction, scalar(0)), eta)
+      + p_*pow(max(alpha1 - alphaMinFriction, scalar(0)), eta_)
-    )/pow(max(alphaMax - alpha1, scalar(5.0e-2)), p + 1.0);
+    )/pow(max(alphaMax - alpha1, scalar(5.0e-2)), p_ + 1.0);
 }
 
 
 Foam::tmp<Foam::volScalarField>
-Foam::kineticTheoryModels::frictionalStressModels::JohnsonJackson::muf
+Foam::kineticTheoryModels::frictionalStressModels::JohnsonJackson::nu
 (
     const volScalarField& alpha1,
     const dimensionedScalar& alphaMax,
     const volScalarField& pf,
-    const volSymmTensorField& D,
+    const volSymmTensorField& D
-    const dimensionedScalar& phi
 ) const
 {
-    return dimensionedScalar("0.5", dimTime, 0.5)*pf*sin(phi);
+    return dimensionedScalar("0.5", dimTime, 0.5)*pf*sin(phi_);
+}
+
+
+bool Foam::kineticTheoryModels::frictionalStressModels::JohnsonJackson::read()
+{
+    coeffDict_ <<= dict_.subDict(typeName + "Coeffs");
+
+    Fr_.readIfPresent(coeffDict_);
+    eta_.readIfPresent(coeffDict_);
+    p_.readIfPresent(coeffDict_);
+    phi_.readIfPresent(coeffDict_);
+
+    return true;
 }
 
 

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/kineticTheoryModels/frictionalStressModel/JohnsonJackson/JohnsonJacksonFrictionalStress.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -53,6 +53,22 @@ class JohnsonJackson
 :
     public frictionalStressModel
 {
+    // Private data
+
+        dictionary coeffDict_;
+
+        //- Material constant for frictional normal stress
+        dimensionedScalar Fr_;
+
+        //- Material constant for frictional normal stress
+        dimensionedScalar eta_;
+
+        //- Material constant for frictional normal stress
+        dimensionedScalar p_;
+
+        //- Angle of internal friction
+        dimensionedScalar phi_;
+
 
 public:
 
@@ -76,30 +92,25 @@ public:
         (
             const volScalarField& alpha1,
             const dimensionedScalar& alphaMinFriction,
-            const dimensionedScalar& alphaMax,
+            const dimensionedScalar& alphaMax
-            const dimensionedScalar& Fr,
-            const dimensionedScalar& eta,
-            const dimensionedScalar& p
         ) const;
 
         virtual tmp<volScalarField> frictionalPressurePrime
         (
             const volScalarField& alpha1,
             const dimensionedScalar& alphaMinFriction,
-            const dimensionedScalar& alphaMax,
+            const dimensionedScalar& alphaMax
-            const dimensionedScalar& Fr,
-            const dimensionedScalar& eta,
-            const dimensionedScalar& p
         ) const;
 
-        virtual tmp<volScalarField> muf
+        virtual tmp<volScalarField> nu
         (
             const volScalarField& alpha1,
             const dimensionedScalar& alphaMax,
             const volScalarField& pf,
-            const volSymmTensorField& D,
+            const volSymmTensorField& D
-            const dimensionedScalar& phi
         ) const;
+
+        virtual bool read();
 };
 
 

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -54,7 +54,9 @@ Foam::kineticTheoryModels::frictionalStressModels::Schaeffer::Schaeffer
     const dictionary& dict
 )
 :
-    frictionalStressModel(dict)
+    frictionalStressModel(dict),
+    coeffDict_(dict.subDict(typeName + "Coeffs")),
+    phi_("phi", dimless, coeffDict_.lookup("phi"))
 {}
 
 
@@ -72,10 +74,7 @@ frictionalPressure
 (
     const volScalarField& alpha1,
     const dimensionedScalar& alphaMinFriction,
-    const dimensionedScalar& alphaMax,
+    const dimensionedScalar& alphaMax
-    const dimensionedScalar& Fr,
-    const dimensionedScalar& eta,
-    const dimensionedScalar& p
 ) const
 {
     return
@@ -90,10 +89,7 @@ frictionalPressurePrime
 (
     const volScalarField& alpha1,
     const dimensionedScalar& alphaMinFriction,
-    const dimensionedScalar& alphaMax,
+    const dimensionedScalar& alphaMax
-    const dimensionedScalar& Fr,
-    const dimensionedScalar& eta,
-    const dimensionedScalar& p
 ) const
 {
     return
@@ -103,42 +99,44 @@ frictionalPressurePrime
 
 
 Foam::tmp<Foam::volScalarField>
-Foam::kineticTheoryModels::frictionalStressModels::Schaeffer::muf
+Foam::kineticTheoryModels::frictionalStressModels::Schaeffer::nu
 (
     const volScalarField& alpha1,
     const dimensionedScalar& alphaMax,
     const volScalarField& pf,
-    const volSymmTensorField& D,
+    const volSymmTensorField& D
-    const dimensionedScalar& phi
 ) const
 {
     const scalar I2Dsmall = 1.0e-15;
 
-    // Creating muf assuming it should be 0 on the boundary which may not be
+    // Creating nu assuming it should be 0 on the boundary which may not be
     // true
-    tmp<volScalarField> tmuf
+    tmp<volScalarField> tnu
     (
         new volScalarField
         (
             IOobject
             (
-                "muf",
+                "Schaeffer:nu",
                 alpha1.mesh().time().timeName(),
-                alpha1.mesh()
+                alpha1.mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
             ),
             alpha1.mesh(),
-            dimensionedScalar("muf", dimensionSet(1, -1, -1, 0, 0), 0.0)
+            dimensionedScalar("nu", dimensionSet(0, 2, -1, 0, 0), 0.0)
         )
     );
 
-    volScalarField& muff = tmuf();
+    volScalarField& nuf = tnu();
 
     forAll (D, celli)
     {
         if (alpha1[celli] > alphaMax.value() - 5e-2)
         {
-            muff[celli] =
+            nuf[celli] =
-                0.5*pf[celli]*sin(phi.value())
+                0.5*pf[celli]*sin(phi_.value())
                /(
                     sqrt(1.0/6.0*(sqr(D[celli].xx() - D[celli].yy())
                   + sqr(D[celli].yy() - D[celli].zz())
@@ -149,9 +147,20 @@ Foam::kineticTheoryModels::frictionalStressModels::Schaeffer::muf
         }
     }
 
-    muff.correctBoundaryConditions();
+    // Correct coupled BCs
+    nuf.correctBoundaryConditions();
 
-    return tmuf;
+    return tnu;
+}
+
+
+bool Foam::kineticTheoryModels::frictionalStressModels::Schaeffer::read()
+{
+    coeffDict_ <<= dict_.subDict(typeName + "Coeffs");
+
+    phi_.readIfPresent(coeffDict_);
+
+    return true;
 }
 
 

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/kineticTheoryModels/frictionalStressModel/Schaeffer/SchaefferFrictionalStress.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -53,6 +53,13 @@ class Schaeffer
 :
     public frictionalStressModel
 {
+    // Private data
+
+        dictionary coeffDict_;
+
+        //- Angle of internal friction
+        dimensionedScalar phi_;
+
 
 public:
 
@@ -76,30 +83,25 @@ public:
         (
             const volScalarField& alpha1,
             const dimensionedScalar& alphaMinFriction,
-            const dimensionedScalar& alphaMax,
+            const dimensionedScalar& alphaMax
-            const dimensionedScalar& Fr,
-            const dimensionedScalar& eta,
-            const dimensionedScalar& p
         ) const;
 
         virtual tmp<volScalarField> frictionalPressurePrime
         (
             const volScalarField& alpha1,
             const dimensionedScalar& alphaMinFriction,
-            const dimensionedScalar& alphaMax,
+            const dimensionedScalar& alphaMax
-            const dimensionedScalar& Fr,
-            const dimensionedScalar& n,
-            const dimensionedScalar& p
         ) const;
 
-        virtual tmp<volScalarField> muf
+        virtual tmp<volScalarField> nu
         (
             const volScalarField& alpha1,
             const dimensionedScalar& alphaMax,
             const volScalarField& pf,
-            const volSymmTensorField& D,
+            const volSymmTensorField& D
-            const dimensionedScalar& phi
         ) const;
+
+        virtual bool read();
 };
 
 

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/kineticTheoryModels/frictionalStressModel/frictionalStressModel/frictionalStressModel.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/kineticTheoryModels/frictionalStressModel/frictionalStressModel/frictionalStressModel.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -63,6 +63,7 @@ protected:
 
     // Protected data
 
+        //- Reference to higher-level dictionary for re-read
         const dictionary& dict_;
 
 
@@ -108,30 +109,25 @@ public:
         (
             const volScalarField& alpha1,
             const dimensionedScalar& alphaMinFriction,
-            const dimensionedScalar& alphaMax,
+            const dimensionedScalar& alphaMax
-            const dimensionedScalar& Fr,
-            const dimensionedScalar& eta,
-            const dimensionedScalar& p
         ) const = 0;
 
         virtual tmp<volScalarField> frictionalPressurePrime
         (
             const volScalarField& alpha1f,
             const dimensionedScalar& alphaMinFriction,
-            const dimensionedScalar& alphaMax,
+            const dimensionedScalar& alphaMax
-            const dimensionedScalar& Fr,
-            const dimensionedScalar& eta,
-            const dimensionedScalar& p
         ) const = 0;
 
-        virtual tmp<volScalarField> muf
+        virtual tmp<volScalarField> nu
         (
             const volScalarField& alpha1,
             const dimensionedScalar& alphaMax,
             const volScalarField& pf,
-            const volSymmTensorField& D,
+            const volSymmTensorField& D
-            const dimensionedScalar& phi
         ) const = 0;
+
+        virtual bool read() = 0;
 };
 
 

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/kineticTheoryModels/granularPressureModel/granularPressureModel/granularPressureModel.H

@@ -122,6 +122,11 @@ public:
             const volScalarField& rho1,
             const dimensionedScalar& e
         ) const = 0;
+
+        virtual bool read()
+        {
+            return true;
+        }
 };
 
 

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

@@ -0,0 +1,546 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "kineticTheoryModel.H"
+#include "surfaceInterpolate.H"
+#include "mathematicalConstants.H"
+#include "twoPhaseSystem.H"
+#include "fvcDiv.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::kineticTheoryModel::kineticTheoryModel
+(
+    const volScalarField& alpha,
+    const geometricOneField& rho,
+    const volVectorField& U,
+    const surfaceScalarField& alphaPhi,
+    const surfaceScalarField& phi,
+    const transportModel& phase,
+    const word& propertiesName,
+    const word& type
+)
+:
+    RASModels::eddyViscosity<PhaseIncompressibleTurbulenceModel<phaseModel> >
+    (
+        type,
+        alpha,
+        rho,
+        U,
+        alphaPhi,
+        phi,
+        phase,
+        propertiesName
+    ),
+
+    phase_(phase),
+
+    draga_(phase.fluid().drag1()),
+
+    viscosityModel_
+    (
+        kineticTheoryModels::viscosityModel::New
+        (
+            this->coeffDict_
+        )
+    ),
+    conductivityModel_
+    (
+        kineticTheoryModels::conductivityModel::New
+        (
+            this->coeffDict_
+        )
+    ),
+    radialModel_
+    (
+        kineticTheoryModels::radialModel::New
+        (
+            this->coeffDict_
+        )
+    ),
+    granularPressureModel_
+    (
+        kineticTheoryModels::granularPressureModel::New
+        (
+            this->coeffDict_
+        )
+    ),
+    frictionalStressModel_
+    (
+        kineticTheoryModels::frictionalStressModel::New
+        (
+            this->coeffDict_
+        )
+    ),
+
+    equilibrium_(this->coeffDict_.lookup("equilibrium")),
+    e_("e", dimless, this->coeffDict_.lookup("e")),
+    alphaMax_("alphaMax", dimless, this->coeffDict_.lookup("alphaMax")),
+    alphaMinFriction_
+    (
+        "alphaMinFriction",
+        dimless,
+        this->coeffDict_.lookup("alphaMinFriction")
+    ),
+
+    Theta_
+    (
+        IOobject
+        (
+            IOobject::groupName("Theta", phase.name()),
+            U.time().timeName(),
+            U.mesh(),
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        U.mesh()
+    ),
+
+    lambda_
+    (
+        IOobject
+        (
+            IOobject::groupName("lambda", phase.name()),
+            U.time().timeName(),
+            U.mesh(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        U.mesh(),
+        dimensionedScalar("zero", dimensionSet(0, 2, -1, 0, 0), 0.0)
+    )
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::kineticTheoryModel::~kineticTheoryModel()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::kineticTheoryModel::correct()
+{
+    // Local references
+    volScalarField alpha(max(this->alpha_, 0.0));
+    const volScalarField& rho = phase_.rho();
+    const surfaceScalarField& alphaPhi = this->alphaPhi_;
+    const volVectorField& U = this->U_;
+    const volVectorField& Uc_ = phase_.fluid().otherPhase(phase_).U();
+
+    const scalar sqrtPi = sqrt(constant::mathematical::pi);
+    dimensionedScalar ThetaSmall("ThetaSmall", Theta_.dimensions(), 1.0e-6);
+    dimensionedScalar ThetaSmallSqrt(sqrt(ThetaSmall));
+
+    tmp<volScalarField> tda(phase_.d());
+    const volScalarField& da = tda();
+
+    tmp<volTensorField> tgradU(fvc::grad(this->U_));
+    const volTensorField& gradU(tgradU());
+    volSymmTensorField D(symm(gradU));
+
+    // Calculating the radial distribution function
+    volScalarField gs0(radialModel_->g0(alpha, alphaMinFriction_, alphaMax_));
+
+    if (!equilibrium_)
+    {
+        // particle viscosity (Table 3.2, p.47)
+        nut_ = viscosityModel_->nu(alpha, Theta_, gs0, rho, da, e_);
+
+        volScalarField ThetaSqrt(sqrt(Theta_));
+
+        // Bulk viscosity  p. 45 (Lun et al. 1984).
+        lambda_ = (4.0/3.0)*sqr(alpha)*da*gs0*(1.0 + e_)*ThetaSqrt/sqrtPi;
+
+        // Stress tensor, Definitions, Table 3.1, p. 43
+        volSymmTensorField tau(2.0*nut_*D + (lambda_ - (2.0/3.0)*nut_)*tr(D)*I);
+
+        // Dissipation (Eq. 3.24, p.50)
+        volScalarField gammaCoeff
+        (
+            12.0*(1.0 - sqr(e_))
+           *max(sqr(alpha), phase_.fluid().residualPhaseFraction())
+           *gs0*(1.0/da)*ThetaSqrt/sqrtPi
+        );
+
+        // NB, drag = K*alpha*alpha2,
+        // (the alpha and alpha2 has been extracted from the drag function for
+        // numerical reasons)
+        volScalarField magUr(mag(U - Uc_));
+
+        volScalarField alpha2Prim
+        (
+            max
+            (
+                alpha*(1.0 - alpha),
+                phase_.fluid().residualPhaseFraction()
+            )*draga_.K(magUr + phase_.fluid().residualSlip())/rho
+        );
+
+        // Eq. 3.25, p. 50 Js = J1 - J2
+        volScalarField J1(3.0*alpha2Prim);
+        volScalarField J2
+        (
+            0.25*sqr(alpha2Prim)*da*sqr(magUr)
+           /(
+               max(alpha, phase_.fluid().residualPhaseFraction())
+              *sqrtPi*(ThetaSqrt + ThetaSmallSqrt)
+            )
+        );
+
+        // particle pressure - coefficient in front of Theta (Eq. 3.22, p. 45)
+        volScalarField PsCoeff
+        (
+            granularPressureModel_->granularPressureCoeff
+            (
+                alpha,
+                gs0,
+                rho,
+                e_
+            )/rho
+        );
+
+        // 'thermal' conductivity (Table 3.3, p. 49)
+        volScalarField kappa
+        (
+            conductivityModel_->kappa(alpha, Theta_, gs0, rho, da, e_)
+        );
+
+        // Construct the granular temperature equation (Eq. 3.20, p. 44)
+        // NB. note that there are two typos in Eq. 3.20:
+        //     Ps should be without grad
+        //     the laplacian has the wrong sign
+        fvScalarMatrix ThetaEqn
+        (
+            1.5*
+            (
+                fvm::ddt(alpha, Theta_)
+              + fvm::div(alphaPhi, Theta_)
+              - fvc::Sp(fvc::ddt(alpha) + fvc::div(alphaPhi), Theta_)
+            )
+          - fvm::laplacian(kappa, Theta_, "laplacian(kappa, Theta)")
+         ==
+            fvm::SuSp(-((PsCoeff*I) && gradU), Theta_)
+          + (tau && gradU)
+          + fvm::Sp(-gammaCoeff, Theta_)
+          + fvm::Sp(-J1, Theta_)
+          + fvm::Sp(J2/(Theta_ + ThetaSmall), Theta_)
+        );
+
+        ThetaEqn.relax();
+        ThetaEqn.solve();
+    }
+    else
+    {
+        // Equilibrium => dissipation == production
+        // Eq. 4.14, p.82
+        volScalarField K1(2.0*(1.0 + e_)*rho*gs0);
+        volScalarField K3
+        (
+            0.5*da*rho*
+            (
+                (sqrtPi/(3.0*(3.0-e_)))
+               *(1.0 + 0.4*(1.0 + e_)*(3.0*e_ - 1.0)*alpha*gs0)
+               +1.6*alpha*gs0*(1.0 + e_)/sqrtPi
+            )
+        );
+
+        volScalarField K2
+        (
+            4.0*da*rho*(1.0 + e_)*alpha*gs0/(3.0*sqrtPi) - 2.0*K3/3.0
+        );
+
+        volScalarField K4(12.0*(1.0 - sqr(e_))*rho*gs0/(da*sqrtPi));
+
+        volScalarField trD
+        (
+            alpha/(alpha + phase_.fluid().residualPhaseFraction())
+           *fvc::div(this->phi_)
+        );
+        volScalarField tr2D(sqr(trD));
+        volScalarField trD2(tr(D & D));
+
+        volScalarField t1(K1*alpha + rho);
+        volScalarField l1(-t1*trD);
+        volScalarField l2(sqr(t1)*tr2D);
+        volScalarField l3
+        (
+            4.0
+           *K4
+           *alpha
+           *(2.0*K3*trD2 + K2*tr2D)
+        );
+
+        Theta_ = sqr
+        (
+            (l1 + sqrt(l2 + l3))
+           /(2.0*max(alpha, phase_.fluid().residualPhaseFraction())*K4)
+        );
+    }
+
+    Theta_.max(0);
+    Theta_.min(100);
+
+    {
+        // particle viscosity (Table 3.2, p.47)
+        nut_ = viscosityModel_->nu(alpha, Theta_, gs0, rho, da, e_);
+
+        volScalarField ThetaSqrt(sqrt(Theta_));
+
+        // Bulk viscosity  p. 45 (Lun et al. 1984).
+        lambda_ = (4.0/3.0)*sqr(alpha)*da*gs0*(1.0 + e_)*ThetaSqrt/sqrtPi;
+
+        // Frictional pressure
+        volScalarField pf
+        (
+            frictionalStressModel_->frictionalPressure
+            (
+                alpha,
+                alphaMinFriction_,
+                alphaMax_
+            )
+        );
+
+        // Add frictional shear viscosity, Eq. 3.30, p. 52
+        nut_ += frictionalStressModel_->nu
+        (
+            alpha,
+            alphaMax_,
+            pf/rho,
+            D
+        );
+
+        // Limit viscosity
+        nut_.min(100);
+    }
+
+    if (debug)
+    {
+        Info<< typeName << ':' << nl
+            << "    max(Theta) = " << max(Theta_).value() << nl
+            << "    max(nut) = " << max(nut_).value() << endl;
+    }
+}
+
+
+/*
+Foam::tmp<Foam::volScalarField> Foam::kineticTheoryModel::pp() const
+{
+
+    // Particle pressure coefficient
+    // Coefficient in front of Theta (Eq. 3.22, p. 45)
+    volScalarField PsCoeff
+    (
+        granularPressureModel_->granularPressureCoeff
+        (
+            alpha,
+            gs0,
+            rho,
+            e_
+        )
+    );
+
+    // Frictional pressure
+    volScalarField pf
+    (
+        frictionalStressModel_->frictionalPressure
+        (
+            alpha,
+            alphaMinFriction_,
+            alphaMax_
+        )
+    );
+
+    // Return total particle pressure
+    return PsCoeff*Theta_ + pf;
+}
+*/
+
+
+Foam::tmp<Foam::volScalarField> Foam::kineticTheoryModel::pPrime() const
+{
+    // Local references
+    const volScalarField& alpha = this->alpha_;
+    const volScalarField& rho = phase_.rho();
+
+    return
+    (
+        Theta_
+       *granularPressureModel_->granularPressureCoeffPrime
+        (
+            alpha,
+            radialModel_->g0(alpha, alphaMinFriction_, alphaMax_),
+            radialModel_->g0prime(alpha, alphaMinFriction_, alphaMax_),
+            rho,
+            e_
+        )
+     +  frictionalStressModel_->frictionalPressurePrime
+        (
+            alpha,
+            alphaMinFriction_,
+            alphaMax_
+        )
+    );
+}
+
+
+Foam::tmp<Foam::surfaceScalarField> Foam::kineticTheoryModel::pPrimef() const
+{
+    // Local references
+    const volScalarField& alpha = this->alpha_;
+    const volScalarField& rho = phase_.rho();
+
+    return fvc::interpolate
+    (
+        Theta_
+       *granularPressureModel_->granularPressureCoeffPrime
+        (
+            alpha,
+            radialModel_->g0(alpha, alphaMinFriction_, alphaMax_),
+            radialModel_->g0prime(alpha, alphaMinFriction_, alphaMax_),
+            rho,
+            e_
+        )
+     +  frictionalStressModel_->frictionalPressurePrime
+        (
+            alpha,
+            alphaMinFriction_,
+            alphaMax_
+        )
+    );
+}
+
+
+void Foam::kineticTheoryModel::correctNut()
+{}
+
+
+Foam::tmp<Foam::volScalarField> Foam::kineticTheoryModel::k() const
+{
+    notImplemented("kineticTheoryModel::k()");
+    return nut_;
+}
+
+
+Foam::tmp<Foam::volScalarField> Foam::kineticTheoryModel::epsilon() const
+{
+    notImplemented("kineticTheoryModel::epsilon()");
+    return nut_;
+}
+
+
+Foam::tmp<Foam::volSymmTensorField> Foam::kineticTheoryModel::R() const
+{
+    return tmp<volSymmTensorField>
+    (
+        new volSymmTensorField
+        (
+            IOobject
+            (
+                IOobject::groupName("R", this->U_.group()),
+                this->runTime_.timeName(),
+                this->mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+          - (this->nut_)*dev(twoSymm(fvc::grad(this->U_)))
+          - (lambda_*fvc::div(this->phi_))*symmTensor::I
+        )
+    );
+}
+
+
+Foam::tmp<Foam::volSymmTensorField> Foam::kineticTheoryModel::devRhoReff() const
+{
+    return tmp<volSymmTensorField>
+    (
+        new volSymmTensorField
+        (
+            IOobject
+            (
+                IOobject::groupName("devRhoReff", this->U_.group()),
+                this->runTime_.timeName(),
+                this->mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+          - (this->rho_*this->nut_)
+           *dev(twoSymm(fvc::grad(this->U_)))
+          - ((this->rho_*lambda_)*fvc::div(this->phi_))*symmTensor::I
+        )
+    );
+}
+
+
+Foam::tmp<Foam::fvVectorMatrix> Foam::kineticTheoryModel::divDevRhoReff
+(
+    volVectorField& U
+) const
+{
+    return
+    (
+      - fvm::laplacian(this->rho_*this->nut_, U)
+      - fvc::div
+        (
+            (this->rho_*this->nut_)*dev2(T(fvc::grad(U)))
+          + ((this->rho_*lambda_)*fvc::div(this->phi_))
+           *dimensioned<symmTensor>("I", dimless, symmTensor::I)
+        )
+    );
+}
+
+
+bool Foam::kineticTheoryModel::read()
+{
+    if
+    (
+        RASModels::eddyViscosity
+        <
+            PhaseIncompressibleTurbulenceModel<phaseModel>
+        >::read()
+    )
+    {
+        this->coeffDict().lookup("equilibrium") >> equilibrium_;
+        e_.readIfPresent(this->coeffDict());
+        alphaMax_.readIfPresent(this->coeffDict());
+        alphaMinFriction_.readIfPresent(this->coeffDict());
+
+        viscosityModel_->read();
+        conductivityModel_->read();
+        radialModel_->read();
+        granularPressureModel_->read();
+        frictionalStressModel_->read();
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+// ************************************************************************* //

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/kineticTheoryModels/kineticTheoryModel/kineticTheoryModel.H

@@ -0,0 +1,214 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::kineticTheoryModel
+
+Description
+
+SourceFiles
+    kineticTheoryModel.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef kineticTheoryModel_H
+#define kineticTheoryModel_H
+
+#include "eddyViscosity.H"
+#include "PhaseIncompressibleTurbulenceModel.H"
+#include "dragModel.H"
+#include "phaseModel.H"
+#include "autoPtr.H"
+#include "viscosityModel.H"
+#include "conductivityModel.H"
+#include "radialModel.H"
+#include "granularPressureModel.H"
+#include "frictionalStressModel.H"
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class kineticTheoryModel Declaration
+\*---------------------------------------------------------------------------*/
+
+class kineticTheoryModel
+:
+    public RASModels::eddyViscosity
+    <
+        PhaseIncompressibleTurbulenceModel<phaseModel>
+    >
+{
+    // Private data
+
+        // Input Fields
+
+            const phaseModel& phase_;
+
+
+        //- Drag model
+        const dragModel& draga_;
+
+        // Sub-models
+
+            //- Run-time selected viscosity model
+            autoPtr<kineticTheoryModels::viscosityModel> viscosityModel_;
+
+            //- Run-time selected conductivity model
+            autoPtr<kineticTheoryModels::conductivityModel> conductivityModel_;
+
+            //- Run-time selected radial distribution model
+            autoPtr<kineticTheoryModels::radialModel> radialModel_;
+
+            //- Run-time selected granular pressure model
+            autoPtr<kineticTheoryModels::granularPressureModel>
+                granularPressureModel_;
+
+            //- Run-time selected frictional stress model
+            autoPtr<kineticTheoryModels::frictionalStressModel>
+                frictionalStressModel_;
+
+
+        // Kinetic Theory Model coefficients
+
+            //- Use equilibrium approximation: generation == dissipation
+            Switch equilibrium_;
+
+            //- Coefficient of restitution
+            dimensionedScalar e_;
+
+            //- Maximum packing phase-fraction
+            dimensionedScalar alphaMax_;
+
+            //- Min value for which the frictional stresses are zero
+            dimensionedScalar alphaMinFriction_;
+
+
+        // Kinetic Theory Model Fields
+
+            //- The granular energy/temperature
+            volScalarField Theta_;
+
+            //- The granular bulk viscosity
+            volScalarField lambda_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        kineticTheoryModel(const kineticTheoryModel&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const kineticTheoryModel&);
+
+
+protected:
+
+    // Protected member functions
+
+        virtual void correctNut();
+
+
+public:
+
+    //- Runtime type information
+    TypeName("kineticTheory");
+
+
+    // Constructors
+
+        //- Construct from components
+        kineticTheoryModel
+        (
+            const volScalarField& alpha,
+            const geometricOneField& rho,
+            const volVectorField& U,
+            const surfaceScalarField& alphaPhi,
+            const surfaceScalarField& phi,
+            const phaseModel& transport,
+            const word& propertiesName = turbulenceModel::propertiesName,
+            const word& type = typeName
+        );
+
+
+    //- Destructor
+    virtual ~kineticTheoryModel();
+
+
+    // Member Functions
+
+        //- Return the effective viscosity
+        virtual tmp<volScalarField> nuEff() const
+        {
+            return this->nut();
+        }
+
+        //- Return the effective viscosity on patch
+        virtual tmp<scalarField> nuEff(const label patchi) const
+        {
+            return this->nut(patchi);
+        }
+
+        //- Return the turbulence kinetic energy
+        virtual tmp<volScalarField> k() const;
+
+        //- Return the turbulence kinetic energy dissipation rate
+        virtual tmp<volScalarField> epsilon() const;
+
+        //- Return the Reynolds stress tensor
+        virtual tmp<volSymmTensorField> R() const;
+
+        //- Return the phase-pressure'
+        // (derivative of phase-pressure w.r.t. phase-fraction)
+        virtual tmp<volScalarField> pPrime() const;
+
+        //- Return the face-phase-pressure'
+        // (derivative of phase-pressure w.r.t. phase-fraction)
+        virtual tmp<surfaceScalarField> pPrimef() const;
+
+        //- Return the effective stress tensor
+        virtual tmp<volSymmTensorField> devRhoReff() const;
+
+        //- Return the source term for the momentum equation
+        virtual tmp<fvVectorMatrix> divDevRhoReff(volVectorField& U) const;
+
+        //- Solve the kinetic theory equations and correct the viscosity
+        virtual void correct();
+
+        //- Re-read model coefficients if they have changed
+        virtual bool read();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/kineticTheoryModels/radialModel/CarnahanStarling/CarnahanStarlingRadial.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -70,6 +70,7 @@ Foam::tmp<Foam::volScalarField>
 Foam::kineticTheoryModels::radialModels::CarnahanStarling::g0
 (
     const volScalarField& alpha,
+    const dimensionedScalar& alphaMinFriction,
     const dimensionedScalar& alphaMax
 ) const
 {
@@ -77,7 +78,7 @@ Foam::kineticTheoryModels::radialModels::CarnahanStarling::g0
     return
         1.0/(1.0 - alpha)
       + 3.0*alpha/(2.0*sqr(1.0 - alpha))
-      + sqr(alpha)/(2.0*pow(1.0 - alpha, 3));
+      + sqr(alpha)/(2.0*pow3(1.0 - alpha));
 }
 
 
@@ -85,13 +86,14 @@ Foam::tmp<Foam::volScalarField>
 Foam::kineticTheoryModels::radialModels::CarnahanStarling::g0prime
 (
     const volScalarField& alpha,
+    const dimensionedScalar& alphaMinFriction,
     const dimensionedScalar& alphaMax
 ) const
 {
     return
         2.5/sqr(1.0 - alpha)
-      + 4.0*alpha/pow(1.0 - alpha, 3.0)
+      + 4.0*alpha/pow3(1.0 - alpha)
-      + 1.5*sqr(alpha)/pow(1.0 - alpha, 4.0);
+      + 1.5*sqr(alpha)/pow4(1.0 - alpha);
 }
 
 

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/kineticTheoryModels/radialModel/CarnahanStarling/CarnahanStarlingRadial.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -76,12 +76,14 @@ public:
         tmp<volScalarField> g0
         (
             const volScalarField& alpha,
+            const dimensionedScalar& alphaMinFriction,
             const dimensionedScalar& alphaMax
         ) const;
 
         tmp<volScalarField> g0prime
         (
             const volScalarField& alpha,
+            const dimensionedScalar& alphaMinFriction,
             const dimensionedScalar& alphaMax
         ) const;
 };

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/kineticTheoryModels/radialModel/LunSavage/LunSavageRadial.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -70,6 +70,7 @@ Foam::tmp<Foam::volScalarField>
 Foam::kineticTheoryModels::radialModels::LunSavage::g0
 (
     const volScalarField& alpha,
+    const dimensionedScalar& alphaMinFriction,
     const dimensionedScalar& alphaMax
 ) const
 {
@@ -82,10 +83,11 @@ Foam::tmp<Foam::volScalarField>
 Foam::kineticTheoryModels::radialModels::LunSavage::g0prime
 (
     const volScalarField& alpha,
+    const dimensionedScalar& alphaMinFriction,
     const dimensionedScalar& alphaMax
 ) const
 {
-    return 2.5*pow(1.0 - alpha/alphaMax, -1.0 - 2.5*alphaMax);
+    return 2.5*pow(1.0 - alpha/alphaMax, -2.5*alphaMax - 1);
 }
 
 

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/kineticTheoryModels/radialModel/LunSavage/LunSavageRadial.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -75,12 +75,14 @@ public:
         tmp<volScalarField> g0
         (
             const volScalarField& alpha,
+            const dimensionedScalar& alphaMinFriction,
             const dimensionedScalar& alphaMax
         ) const;
 
         tmp<volScalarField> g0prime
         (
             const volScalarField& alpha,
+            const dimensionedScalar& alphaMinFriction,
             const dimensionedScalar& alphaMax
         ) const;
 };

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/kineticTheoryModels/radialModel/SinclairJackson/SinclairJacksonRadial.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -70,10 +70,11 @@ Foam::tmp<Foam::volScalarField>
 Foam::kineticTheoryModels::radialModels::SinclairJackson::g0
 (
     const volScalarField& alpha,
+    const dimensionedScalar& alphaMinFriction,
     const dimensionedScalar& alphaMax
 ) const
 {
-    return 1.0/(1.0 - pow(alpha/alphaMax, 1.0/3.0));
+    return 1.0/(1.0 - cbrt(min(alpha, alphaMinFriction)/alphaMax));
 }
 
 
@@ -81,12 +82,16 @@ Foam::tmp<Foam::volScalarField>
 Foam::kineticTheoryModels::radialModels::SinclairJackson::g0prime
 (
     const volScalarField& alpha,
+    const dimensionedScalar& alphaMinFriction,
     const dimensionedScalar& alphaMax
 ) const
 {
-    return
+    volScalarField aByaMax
-       (1.0/3.0)*pow(max(alpha, scalar(1.0e-6))/alphaMax, -2.0/3.0)
+    (
-      /(alphaMax*sqr(1.0 - pow(alpha/alphaMax, 1.0/3.0)));
+        cbrt(min(max(alpha, 1e-3), alphaMinFriction)/alphaMax)
+    );
+
+    return (1.0/(3*alphaMax))/sqr(aByaMax - sqr(aByaMax));
 }
 
 

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/kineticTheoryModels/radialModel/SinclairJackson/SinclairJacksonRadial.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -75,12 +75,14 @@ public:
         tmp<volScalarField> g0
         (
             const volScalarField& alpha,
+            const dimensionedScalar& alphaMinFriction,
             const dimensionedScalar& alphaMax
         ) const;
 
         tmp<volScalarField> g0prime
         (
             const volScalarField& alpha,
+            const dimensionedScalar& alphaMinFriction,
             const dimensionedScalar& alphaMax
         ) const;
 };

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/kineticTheoryModels/radialModel/radialModel/radialModel.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/kineticTheoryModels/radialModel/radialModel/radialModel.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -108,6 +108,7 @@ public:
         virtual tmp<volScalarField> g0
         (
             const volScalarField& alpha,
+            const dimensionedScalar& alphaMinFriction,
             const dimensionedScalar& alphaMax
         ) const = 0;
 
@@ -115,8 +116,14 @@ public:
         virtual tmp<volScalarField> g0prime
         (
             const volScalarField& alpha,
+            const dimensionedScalar& alphaMinFriction,
             const dimensionedScalar& alphaMax
         ) const = 0;
+
+        virtual bool read()
+        {
+            return true;
+        }
 };
 
 

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/kineticTheoryModels/viscosityModel/Gidaspow/GidaspowViscosity.C

@@ -62,7 +62,7 @@ Foam::kineticTheoryModels::viscosityModels::Gidaspow::~Gidaspow()
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 Foam::tmp<Foam::volScalarField>
-Foam::kineticTheoryModels::viscosityModels::Gidaspow::mu1
+Foam::kineticTheoryModels::viscosityModels::Gidaspow::nu
 (
     const volScalarField& alpha1,
     const volScalarField& Theta,
@@ -74,7 +74,7 @@ Foam::kineticTheoryModels::viscosityModels::Gidaspow::mu1
 {
     const scalar sqrtPi = sqrt(constant::mathematical::pi);
 
-    return rho1*da*sqrt(Theta)*
+    return da*sqrt(Theta)*
     (
         (4.0/5.0)*sqr(alpha1)*g0*(1.0 + e)/sqrtPi
       + (1.0/15.0)*sqrtPi*g0*(1.0 + e)*sqr(alpha1)

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/kineticTheoryModels/viscosityModel/Gidaspow/GidaspowViscosity.H

@@ -71,7 +71,7 @@ public:
 
     // Member functions
 
-        tmp<volScalarField> mu1
+        tmp<volScalarField> nu
         (
             const volScalarField& alpha1,
             const volScalarField& Theta,

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/kineticTheoryModels/viscosityModel/HrenyaSinclair/HrenyaSinclairViscosity.C

@@ -56,8 +56,8 @@ Foam::kineticTheoryModels::viscosityModels::HrenyaSinclair::HrenyaSinclair
 )
 :
     viscosityModel(dict),
-    coeffsDict_(dict.subDict(typeName + "Coeffs")),
+    coeffDict_(dict.subDict(typeName + "Coeffs")),
-    L_(coeffsDict_.lookup("L"))
+    L_("L", dimensionSet(0, 1, 0, 0, 0), coeffDict_.lookup("L"))
 {}
 
 
@@ -70,7 +70,7 @@ Foam::kineticTheoryModels::viscosityModels::HrenyaSinclair::~HrenyaSinclair()
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 Foam::tmp<Foam::volScalarField>
-Foam::kineticTheoryModels::viscosityModels::HrenyaSinclair::mu1
+Foam::kineticTheoryModels::viscosityModels::HrenyaSinclair::nu
 (
     const volScalarField& alpha1,
     const volScalarField& Theta,
@@ -87,7 +87,7 @@ Foam::kineticTheoryModels::viscosityModels::HrenyaSinclair::mu1
         scalar(1) + da/(6.0*sqrt(2.0)*(alpha1 + scalar(1.0e-5)))/L_
     );
 
-    return rho1*da*sqrt(Theta)*
+    return da*sqrt(Theta)*
     (
         (4.0/5.0)*sqr(alpha1)*g0*(1.0 + e)/sqrtPi
       + (1.0/15.0)*sqrtPi*g0*(1.0 + e)*(3.0*e - 1)*sqr(alpha1)/(3.0-e)
@@ -98,4 +98,14 @@ Foam::kineticTheoryModels::viscosityModels::HrenyaSinclair::mu1
 }
 
 
+bool Foam::kineticTheoryModels::viscosityModels::HrenyaSinclair::read()
+{
+    coeffDict_ <<= dict_.subDict(typeName + "Coeffs");
+
+    L_.readIfPresent(coeffDict_);
+
+    return true;
+}
+
+
 // ************************************************************************* //

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/kineticTheoryModels/viscosityModel/HrenyaSinclair/HrenyaSinclairViscosity.H

@@ -55,9 +55,9 @@ class HrenyaSinclair
 {
     // Private data
 
-        dictionary coeffsDict_;
+        dictionary coeffDict_;
 
-        //- characteristic length of geometry
+        //- Characteristic length of geometry
         dimensionedScalar L_;
 
 
@@ -79,7 +79,7 @@ public:
 
     // Member functions
 
-        tmp<volScalarField> mu1
+        tmp<volScalarField> nu
         (
             const volScalarField& alpha1,
             const volScalarField& Theta,
@@ -88,6 +88,8 @@ public:
             const volScalarField& da,
             const dimensionedScalar& e
         ) const;
+
+        virtual bool read();
 };
 
 

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/kineticTheoryModels/viscosityModel/Syamlal/SyamlalViscosity.C

@@ -62,7 +62,7 @@ Foam::kineticTheoryModels::viscosityModels::Syamlal::~Syamlal()
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 Foam::tmp<Foam::volScalarField>
-Foam::kineticTheoryModels::viscosityModels::Syamlal::mu1
+Foam::kineticTheoryModels::viscosityModels::Syamlal::nu
 (
     const volScalarField& alpha1,
     const volScalarField& Theta,
@@ -74,7 +74,7 @@ Foam::kineticTheoryModels::viscosityModels::Syamlal::mu1
 {
     const scalar sqrtPi = sqrt(constant::mathematical::pi);
 
-    return rho1*da*sqrt(Theta)*
+    return da*sqrt(Theta)*
     (
         (4.0/5.0)*sqr(alpha1)*g0*(1.0 + e)/sqrtPi
       + (1.0/15.0)*sqrtPi*g0*(1.0 + e)*(3.0*e - 1.0)*sqr(alpha1)/(3.0 - e)

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/kineticTheoryModels/viscosityModel/Syamlal/SyamlalViscosity.H

@@ -72,7 +72,7 @@ public:
 
     // Member functions
 
-        tmp<volScalarField> mu1
+        tmp<volScalarField> nu
         (
             const volScalarField& alpha1,
             const volScalarField& Theta,

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/kineticTheoryModels/viscosityModel/none/noneViscosity.C

@@ -53,7 +53,7 @@ Foam::kineticTheoryModels::noneViscosity::~noneViscosity()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-Foam::tmp<Foam::volScalarField> Foam::kineticTheoryModels::noneViscosity::mu1
+Foam::tmp<Foam::volScalarField> Foam::kineticTheoryModels::noneViscosity::nu
 (
     const volScalarField& alpha1,
     const volScalarField& Theta,
@@ -66,7 +66,7 @@ Foam::tmp<Foam::volScalarField> Foam::kineticTheoryModels::noneViscosity::mu1
     return dimensionedScalar
     (
         "0",
-        dimensionSet(1, -1, -1, 0, 0, 0, 0),
+        dimensionSet(0, 2, -1, 0, 0, 0, 0),
         0.0
     )*alpha1;
 }

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/kineticTheoryModels/viscosityModel/none/noneViscosity.H

@@ -70,7 +70,7 @@ public:
 
     // Member functions
 
-        tmp<volScalarField> mu1
+        tmp<volScalarField> nu
         (
             const volScalarField& alpha1,
             const volScalarField& Theta,

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/kineticTheoryModels/viscosityModel/viscosityModel/viscosityModel.H

@@ -106,7 +106,7 @@ public:
 
     // Member Functions
 
-        virtual tmp<volScalarField> mu1
+        virtual tmp<volScalarField> nu
         (
             const volScalarField& alpha1,
             const volScalarField& Theta,
@@ -115,6 +115,11 @@ public:
             const volScalarField& da,
             const dimensionedScalar& e
         ) const = 0;
+
+        virtual bool read()
+        {
+            return true;
+        }
 };
 
 

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/phaseIncompressibleTurbulenceModels.C

@@ -0,0 +1,172 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "PhaseIncompressibleTurbulenceModel.H"
+#include "laminar.H"
+#include "RASModel.H"
+#include "kEpsilon.H"
+#include "LaheyKEpsilon.H"
+#include "continuousGasKEpsilon.H"
+#include "kineticTheoryModel.H"
+#include "phasePressureModel.H"
+#include "phaseModel.H"
+#include "addToRunTimeSelectionTable.H"
+
+namespace Foam
+{
+    typedef TurbulenceModel
+    <
+        volScalarField,
+        geometricOneField,
+        incompressibleTurbulenceModel,
+        phaseModel
+    > basePhaseIncompressibleTransportTurbulenceModel;
+
+    defineTemplateRunTimeSelectionTable
+    (
+        basePhaseIncompressibleTransportTurbulenceModel,
+        dictionary
+    );
+
+    typedef PhaseIncompressibleTurbulenceModel<phaseModel>
+        incompressibleTransportTurbulenceModel;
+
+    typedef laminar<incompressibleTransportTurbulenceModel>
+        incompressibleLaminar;
+
+    defineNamedTemplateTypeNameAndDebug(incompressibleLaminar, 0);
+
+    addToRunTimeSelectionTable
+    (
+        basePhaseIncompressibleTransportTurbulenceModel,
+        incompressibleLaminar,
+        dictionary
+    );
+
+
+    typedef RASModel<incompressibleTransportTurbulenceModel>
+        incompressibleRASModel;
+
+    defineNamedTemplateTypeNameAndDebug(incompressibleRASModel, 0);
+
+    defineTemplateRunTimeSelectionTable(incompressibleRASModel, dictionary);
+
+    addToRunTimeSelectionTable
+    (
+        basePhaseIncompressibleTransportTurbulenceModel,
+        incompressibleRASModel,
+        dictionary
+    );
+
+    namespace RASModels
+    {
+        typedef kEpsilon<incompressibleTransportTurbulenceModel>
+            incompressiblekEpsilon;
+
+        defineNamedTemplateTypeNameAndDebug(incompressiblekEpsilon, 0);
+
+        addToRunTimeSelectionTable
+        (
+            incompressibleRASModel,
+            incompressiblekEpsilon,
+            dictionary
+        );
+    }
+
+    namespace RASModels
+    {
+        typedef LaheyKEpsilon<incompressibleTransportTurbulenceModel>
+            incompressibleLaheyKEpsilon;
+
+        defineNamedTemplateTypeNameAndDebug(incompressibleLaheyKEpsilon, 0);
+
+        addToRunTimeSelectionTable
+        (
+            incompressibleRASModel,
+            incompressibleLaheyKEpsilon,
+            dictionary
+        );
+    }
+
+    namespace RASModels
+    {
+        typedef continuousGasKEpsilon<incompressibleTransportTurbulenceModel>
+            incompressiblecontinuousGasKEpsilon;
+
+        defineNamedTemplateTypeNameAndDebug
+        (
+            incompressiblecontinuousGasKEpsilon,
+            0
+        );
+
+        addToRunTimeSelectionTable
+        (
+            incompressibleRASModel,
+            incompressiblecontinuousGasKEpsilon,
+            dictionary
+        );
+    }
+}
+
+
+namespace Foam
+{
+    typedef PhaseIncompressibleTurbulenceModel<phaseModel>
+        incompressibleTransportTurbulenceModel;
+
+    typedef RASModel<incompressibleTransportTurbulenceModel>
+        incompressibleRASModel;
+
+    defineTypeNameAndDebug(kineticTheoryModel, 0);
+
+    addToRunTimeSelectionTable
+    (
+        incompressibleRASModel,
+        kineticTheoryModel,
+        dictionary
+    );
+}
+
+
+namespace Foam
+{
+    typedef PhaseIncompressibleTurbulenceModel<phaseModel>
+        incompressibleTransportTurbulenceModel;
+
+    typedef RASModel<incompressibleTransportTurbulenceModel>
+        incompressibleRASModel;
+
+    defineTypeNameAndDebug(phasePressureModel, 0);
+
+    addToRunTimeSelectionTable
+    (
+        incompressibleRASModel,
+        phasePressureModel,
+        dictionary
+    );
+}
+
+
+// ************************************************************************* //

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/phasePressureModel/phasePressureModel.C

@@ -0,0 +1,219 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "phasePressureModel.H"
+#include "surfaceInterpolate.H"
+#include "twoPhaseSystem.H"
+#include "dimensionedType.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::phasePressureModel::phasePressureModel
+(
+    const volScalarField& alpha,
+    const geometricOneField& rho,
+    const volVectorField& U,
+    const surfaceScalarField& alphaPhi,
+    const surfaceScalarField& phi,
+    const transportModel& phase,
+    const word& propertiesName,
+    const word& type
+)
+:
+    RASModels::eddyViscosity<PhaseIncompressibleTurbulenceModel<phaseModel> >
+    (
+        type,
+        alpha,
+        rho,
+        U,
+        alphaPhi,
+        phi,
+        phase,
+        propertiesName
+    ),
+
+    phase_(phase),
+
+    alphaMax_(readScalar(this->coeffDict_.lookup("alphaMax"))),
+    preAlphaExp_(readScalar(this->coeffDict_.lookup("preAlphaExp"))),
+    expMax_(readScalar(this->coeffDict_.lookup("expMax"))),
+    g0_
+    (
+        "g0",
+        dimensionSet(1, -1, -2, 0, 0),
+        this->coeffDict_.lookup("g0")
+    )
+{
+    this->nut_ == dimensionedScalar("zero", this->nut_.dimensions(), 0.0);
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::phasePressureModel::~phasePressureModel()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::phasePressureModel::correct()
+{}
+
+
+void Foam::phasePressureModel::correctNut()
+{}
+
+
+Foam::tmp<Foam::volScalarField> Foam::phasePressureModel::pPrime() const
+{
+    return
+        g0_
+       *min
+        (
+            exp(preAlphaExp_*(this->alpha_ - alphaMax_)),
+            expMax_
+        );
+}
+
+
+Foam::tmp<Foam::surfaceScalarField> Foam::phasePressureModel::pPrimef() const
+{
+    return
+        g0_
+       *min
+        (
+            exp(preAlphaExp_*(fvc::interpolate(this->alpha_) - alphaMax_)),
+            expMax_
+        );
+}
+
+
+Foam::tmp<Foam::volScalarField> Foam::phasePressureModel::k() const
+{
+    notImplemented("phasePressureModel::k()");
+    return nut_;
+}
+
+
+Foam::tmp<Foam::volScalarField> Foam::phasePressureModel::epsilon() const
+{
+    notImplemented("phasePressureModel::epsilon()");
+    return nut_;
+}
+
+
+Foam::tmp<Foam::volSymmTensorField> Foam::phasePressureModel::R() const
+{
+    return tmp<volSymmTensorField>
+    (
+        new volSymmTensorField
+        (
+            IOobject
+            (
+                IOobject::groupName("R", this->U_.group()),
+                this->runTime_.timeName(),
+                this->mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            this->mesh_,
+            dimensioned<symmTensor>
+            (
+                "R",
+                dimensionSet(0, 2, -2, 0, 0),
+                symmTensor::zero
+            )
+        )
+    );
+}
+
+
+Foam::tmp<Foam::volSymmTensorField> Foam::phasePressureModel::devRhoReff() const
+{
+    return tmp<volSymmTensorField>
+    (
+        new volSymmTensorField
+        (
+            IOobject
+            (
+                IOobject::groupName("devRhoReff", this->U_.group()),
+                this->runTime_.timeName(),
+                this->mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            this->mesh_,
+            dimensioned<symmTensor>
+            (
+                "R",
+                this->rho_.dimensions()*dimensionSet(0, 2, -2, 0, 0),
+                symmTensor::zero
+            )
+        )
+    );
+}
+
+
+Foam::tmp<Foam::fvVectorMatrix> Foam::phasePressureModel::divDevRhoReff
+(
+    volVectorField& U
+) const
+{
+    return tmp<fvVectorMatrix>
+    (
+        new fvVectorMatrix
+        (
+            U,
+            this->rho_.dimensions()*dimensionSet(0, 4, -2, 0, 0)
+        )
+    );
+}
+
+
+bool Foam::phasePressureModel::read()
+{
+    if
+    (
+        RASModels::eddyViscosity
+        <
+            PhaseIncompressibleTurbulenceModel<phaseModel>
+        >::read()
+    )
+    {
+        this->coeffDict().lookup("alphaMax") >> alphaMax_;
+        this->coeffDict().lookup("preAlphaExp") >> preAlphaExp_;
+        this->coeffDict().lookup("expMax") >> expMax_;
+        g0_.readIfPresent(this->coeffDict());
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+// ************************************************************************* //

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/phaseIncompressibleTurbulenceModels/phasePressureModel/phasePressureModel.H

@@ -0,0 +1,182 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::phasePressureModel
+
+Description
+
+SourceFiles
+    phasePressureModel.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef phasePressureModel_H
+#define phasePressureModel_H
+
+#include "eddyViscosity.H"
+#include "PhaseIncompressibleTurbulenceModel.H"
+#include "dragModel.H"
+#include "phaseModel.H"
+#include "autoPtr.H"
+#include "viscosityModel.H"
+#include "conductivityModel.H"
+#include "radialModel.H"
+#include "granularPressureModel.H"
+#include "frictionalStressModel.H"
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class phasePressureModel Declaration
+\*---------------------------------------------------------------------------*/
+
+class phasePressureModel
+:
+    public RASModels::eddyViscosity
+    <
+        PhaseIncompressibleTurbulenceModel<phaseModel>
+    >
+{
+    // Private data
+
+        // Input Fields
+
+            const phaseModel& phase_;
+
+
+        // Kinetic Theory Model coefficients
+
+            //- Maximum packing phase-fraction
+            scalar alphaMax_;
+
+            //- Pre-exponential factor
+            scalar preAlphaExp_;
+
+            //- Maximum limit of the exponential
+            scalar expMax_;
+
+            //- g0
+            dimensionedScalar g0_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct
+        phasePressureModel(const phasePressureModel&);
+
+        //- Disallow default bitwise assignment
+        void operator=(const phasePressureModel&);
+
+
+protected:
+
+    // Protected member functions
+
+        virtual void correctNut();
+
+
+public:
+
+    //- Runtime type information
+    TypeName("phasePressure");
+
+
+    // Constructors
+
+        //- Construct from components
+        phasePressureModel
+        (
+            const volScalarField& alpha,
+            const geometricOneField& rho,
+            const volVectorField& U,
+            const surfaceScalarField& alphaPhi,
+            const surfaceScalarField& phi,
+            const phaseModel& transport,
+            const word& propertiesName = turbulenceModel::propertiesName,
+            const word& type = typeName
+        );
+
+
+    //- Destructor
+    virtual ~phasePressureModel();
+
+
+    // Member Functions
+
+        //- Return the effective viscosity
+        virtual tmp<volScalarField> nuEff() const
+        {
+            return this->nut();
+        }
+
+        //- Return the effective viscosity on patch
+        virtual tmp<scalarField> nuEff(const label patchi) const
+        {
+            return this->nut(patchi);
+        }
+
+        //- Return the turbulence kinetic energy
+        virtual tmp<volScalarField> k() const;
+
+        //- Return the turbulence kinetic energy dissipation rate
+        virtual tmp<volScalarField> epsilon() const;
+
+        //- Return the Reynolds stress tensor
+        virtual tmp<volSymmTensorField> R() const;
+
+        //- Return the phase-pressure'
+        // (derivative of phase-pressure w.r.t. phase-fraction)
+        virtual tmp<volScalarField> pPrime() const;
+
+        //- Return the face-phase-pressure'
+        // (derivative of phase-pressure w.r.t. phase-fraction)
+        virtual tmp<surfaceScalarField> pPrimef() const;
+
+        //- Return the effective stress tensor
+        virtual tmp<volSymmTensorField> devRhoReff() const;
+
+        //- Return the source term for the momentum equation
+        virtual tmp<fvVectorMatrix> divDevRhoReff(volVectorField& U) const;
+
+        //- Solve the kinetic theory equations and correct the viscosity
+        virtual void correct();
+
+        //- Re-read model coefficients if they have changed
+        virtual bool read();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/readPPProperties.H

@@ -29,8 +29,3 @@
     (
         ppProperties.lookup("g0")
     );
-
-    Switch packingLimiter
-    (
-        ppProperties.lookup("packingLimiter")
-    );

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/readTwoPhaseEulerFoamControls.H

@@ -1,2 +1,7 @@
     #include "readTimeControls.H"
     #include "alphaControls.H"
+
+    Switch implicitPhasePressure
+    (
+        alphaControls.lookupOrDefault<Switch>("implicitPhasePressure", false)
+    );

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

@@ -1,72 +0,0 @@
-if (turbulence)
-{
-    if (mesh.changing())
-    {
-        y.correct();
-    }
-
-    tmp<volTensorField> tgradU2(fvc::grad(U2));
-    volScalarField G(2*nut2*(tgradU2() && dev(symm(tgradU2()))));
-    tgradU2.clear();
-
-    #include "wallFunctions.H"
-
-    // Dissipation equation
-    fvScalarMatrix epsEqn
-    (
-        fvm::ddt(alpha2, epsilon)
-      + fvm::div(alphaPhi2, epsilon)
-
-        // Compressibity correction
-      - fvm::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), epsilon)
-
-      - fvm::laplacian
-        (
-            alpha1Eps*nuEff2, epsilon,
-            "laplacian(DepsilonEff,epsilon)"
-        )
-      ==
-         C1*alpha2*G*epsilon/k
-       - fvm::Sp(C2*alpha2*epsilon/k, epsilon)
-    );
-
-    #include "wallDissipation.H"
-
-
-
-    epsEqn.relax();
-    epsEqn.solve();
-
-    epsilon.max(dimensionedScalar("zero", epsilon.dimensions(), 1.0e-15));
-
-
-    // Turbulent kinetic energy equation
-    fvScalarMatrix kEqn
-    (
-        fvm::ddt(alpha2, k)
-      + fvm::div(alphaPhi2, k)
-
-        // Compressibity correction
-      - fvm::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), k)
-
-      - fvm::laplacian
-        (
-            alpha1k*nuEff2, k,
-            "laplacian(DkEff,k)"
-        )
-      ==
-        alpha2*G
-      - fvm::Sp(alpha2*epsilon/k, k)
-    );
-    kEqn.relax();
-    kEqn.solve();
-
-    k.max(dimensionedScalar("zero", k.dimensions(), 1.0e-8));
-
-    //- Re-calculate turbulence viscosity
-    nut2 = Cmu*sqr(k)/epsilon;
-
-    #include "wallViscosity.H"
-}
-
-nuEff2 = nut2 + thermo2.mu()/rho2;

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

@@ -1,83 +0,0 @@
-{
-    labelList cellBoundaryFaceCount(epsilon.size(), 0);
-
-    scalar Cmu25 = ::pow(Cmu.value(), 0.25);
-    scalar Cmu75 = ::pow(Cmu.value(), 0.75);
-    scalar kappa_ = kappa.value();
-
-    const fvPatchList& patches = mesh.boundary();
-
-    //- Initialise the near-wall P field to zero
-    forAll(patches, patchi)
-    {
-        const fvPatch& currPatch = patches[patchi];
-
-        if (isA<wallFvPatch>(currPatch))
-        {
-            forAll(currPatch, facei)
-            {
-                label faceCelli = currPatch.faceCells()[facei];
-
-                epsilon[faceCelli] = 0.0;
-                G[faceCelli] = 0.0;
-            }
-        }
-    }
-
-    //- Accumulate the wall face contributions to epsilon and G
-    //  Increment cellBoundaryFaceCount for each face for averaging
-    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))
-        {
-            const scalarField& nut2w = nut2.boundaryField()[patchi];
-
-            scalarField magFaceGradU(mag(U2.boundaryField()[patchi].snGrad()));
-
-            forAll(currPatch, facei)
-            {
-                label faceCelli = currPatch.faceCells()[facei];
-
-                // For corner cells (with two boundary or more faces),
-                // epsilon and G in the near-wall cell are calculated
-                // as an average
-
-                cellBoundaryFaceCount[faceCelli]++;
-
-                epsilon[faceCelli] +=
-                     Cmu75*::pow(k[faceCelli], 1.5)
-                    /(kappa_*y[patchi][facei]);
-
-                G[faceCelli] +=
-                    (nut2w[facei] + mu2_[facei]/rho2_[facei])
-                   *magFaceGradU[facei]
-                   *Cmu25*::sqrt(k[faceCelli])
-                   /(kappa_*y[patchi][facei]);
-            }
-        }
-    }
-
-
-    // perform the averaging
-
-    forAll(patches, patchi)
-    {
-        const fvPatch& curPatch = patches[patchi];
-
-        if (isA<wallFvPatch>(curPatch))
-        {
-            forAll(curPatch, facei)
-            {
-                label faceCelli = curPatch.faceCells()[facei];
-
-                epsilon[faceCelli] /= cellBoundaryFaceCount[faceCelli];
-                G[faceCelli] /= cellBoundaryFaceCount[faceCelli];
-                cellBoundaryFaceCount[faceCelli] = 1;
-            }
-        }
-    }
-}

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

@@ -1,40 +0,0 @@
-{
-    scalar Cmu25 = ::pow(Cmu.value(), 0.25);
-    scalar kappa_ = kappa.value();
-    scalar E_ = E.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))
-        {
-            scalarField& nutw = nut2.boundaryField()[patchi];
-
-            forAll(currPatch, facei)
-            {
-                label faceCelli = currPatch.faceCells()[facei];
-
-                // calculate yPlus
-                scalar yPlus =
-                    Cmu25*y[patchi][facei]*::sqrt(k[faceCelli])
-                   /(mu2_[facei]/rho2_[facei]);
-
-                if (yPlus > 11.6)
-                {
-                    nutw[facei] =
-                        (mu2_[facei]/rho2_[facei])
-                       *(yPlus*kappa_/::log(E_*yPlus) - 1);
-                }
-                else
-                {
-                    nutw[facei] = 0.0;
-                }
-            }
-        }
-    }
-}

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/twoPhaseSystem/Make/files

@@ -4,4 +4,6 @@ diameterModels/diameterModel/newDiameterModel.C
 diameterModels/constantDiameter/constantDiameter.C
 diameterModels/isothermalDiameter/isothermalDiameter.C
 
-LIB = $(FOAM_LIBBIN)/libcompressiblePhaseModel
+twoPhaseSystem.C
+
+LIB = $(FOAM_LIBBIN)/libcompressibleTwoPhaseSystem

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/twoPhaseSystem/Make/options

@@ -1,4 +1,6 @@
 EXE_INC = \
+    -I../twoPhaseSystem \
+    -I../interfacialModels/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/transportModels/incompressible/lnInclude

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

@@ -24,7 +24,10 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "phaseModel.H"
+#include "twoPhaseSystem.H"
 #include "diameterModel.H"
+#include "dragModel.H"
+#include "heatTransferModel.H"
 #include "fixedValueFvPatchFields.H"
 #include "slipFvPatchFields.H"
 #include "surfaceInterpolate.H"
@@ -33,8 +36,8 @@ License
 
 Foam::phaseModel::phaseModel
 (
-    const fvMesh& mesh,
+    const twoPhaseSystem& fluid,
-    const dictionary& transportProperties,
+    const dictionary& phaseProperties,
     const word& phaseName
 )
 :
@@ -42,48 +45,44 @@ Foam::phaseModel::phaseModel
     (
         IOobject
         (
-            "alpha" + phaseName,
+            IOobject::groupName("alpha", phaseName),
-            mesh.time().timeName(),
+            fluid.mesh().time().timeName(),
-            mesh,
+            fluid.mesh(),
             IOobject::READ_IF_PRESENT,
             IOobject::AUTO_WRITE
         ),
-        mesh,
+        fluid.mesh(),
         dimensionedScalar("alpha", dimless, 0)
     ),
+    fluid_(fluid),
     name_(phaseName),
     phaseDict_
     (
-        transportProperties.subDict
+        phaseProperties.subDict(name_)
-        (
-            phaseName == "1" || phaseName == "2"
-          ? "phase" + phaseName
-          : word(phaseName)
-        )
     ),
-    thermo_(rhoThermo::New(mesh, phaseName)),
+    thermo_(rhoThermo::New(fluid.mesh(), name_)),
     U_
     (
         IOobject
         (
-            "U" + phaseName,
+            IOobject::groupName("U", name_),
-            mesh.time().timeName(),
+            fluid.mesh().time().timeName(),
-            mesh,
+            fluid.mesh(),
             IOobject::MUST_READ,
             IOobject::AUTO_WRITE
         ),
-        mesh
+        fluid.mesh()
     )
 {
-    thermo_->validate("phaseModel " + phaseName, "h", "e");
+    thermo_->validate("phaseModel " + name_, "h", "e");
 
-    const word phiName = "phi" + phaseName;
+    const word phiName = IOobject::groupName("phi", name_);
 
     IOobject phiHeader
     (
         phiName,
-        mesh.time().timeName(),
+        fluid_.mesh().time().timeName(),
-        mesh,
+        fluid_.mesh(),
         IOobject::NO_READ
     );
 
@@ -98,12 +97,12 @@ Foam::phaseModel::phaseModel
                 IOobject
                 (
                     phiName,
-                    mesh.time().timeName(),
+                    fluid_.mesh().time().timeName(),
-                    mesh,
+                    fluid_.mesh(),
                     IOobject::MUST_READ,
                     IOobject::AUTO_WRITE
                 ),
-                mesh
+                fluid_.mesh()
             )
         );
     }
@@ -136,12 +135,12 @@ Foam::phaseModel::phaseModel
                 IOobject
                 (
                     phiName,
-                    mesh.time().timeName(),
+                    fluid_.mesh().time().timeName(),
-                    mesh,
+                    fluid_.mesh(),
                     IOobject::NO_READ,
                     IOobject::AUTO_WRITE
                 ),
-                fvc::interpolate(U_) & mesh.Sf(),
+                fvc::interpolate(U_) & fluid_.mesh().Sf(),
                 phiTypes
             )
         );
@@ -155,20 +154,6 @@ Foam::phaseModel::phaseModel
 }
 
 
-Foam::autoPtr<Foam::phaseModel> Foam::phaseModel::New
-(
-    const fvMesh& mesh,
-    const dictionary& transportProperties,
-    const word& phaseName
-)
-{
-    return autoPtr<phaseModel>
-    (
-        new phaseModel(mesh, transportProperties, phaseName)
-    );
-}
-
-
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 Foam::phaseModel::~phaseModel()

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

@@ -36,6 +36,7 @@ SourceFiles
 #include "dimensionedScalar.H"
 #include "volFields.H"
 #include "surfaceFields.H"
+#include "transportModel.H"
 #include "rhoThermo.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -44,18 +45,24 @@ namespace Foam
 {
 
 // Forward declarations
+class twoPhaseSystem;
 class diameterModel;
 
+
 /*---------------------------------------------------------------------------*\
                            Class phaseModel Declaration
 \*---------------------------------------------------------------------------*/
 
 class phaseModel
 :
-    public volScalarField
+    public volScalarField,
+    public transportModel
 {
     // Private data
 
+        //- Reference to the twoPhaseSystem to which this phase belongs
+        const twoPhaseSystem& fluid_;
+
         //- Name of phase
         word name_;
 
@@ -80,19 +87,8 @@ public:
 
         phaseModel
         (
-            const fvMesh& mesh,
+            const twoPhaseSystem& fluid,
-            const dictionary& transportProperties,
+            const dictionary& phaseProperties,
-            const word& phaseName
-        );
-
-
-    // Selectors
-
-        //- Return a reference to the selected turbulence model
-        static autoPtr<phaseModel> New
-        (
-            const fvMesh& mesh,
-            const dictionary& transportProperties,
             const word& phaseName
         );
 
@@ -103,6 +99,12 @@ public:
 
     // Member Functions
 
+        //- Return the twoPhaseSystem to which this phase belongs
+        const twoPhaseSystem& fluid() const
+        {
+            return fluid_;
+        }
+
         const word& name() const
         {
             return name_;
@@ -165,6 +167,16 @@ public:
         {
             return phiPtr_();
         }
+
+        //- Dummy correct
+        void correct()
+        {}
+
+        //- Dummy read
+        bool read()
+        {
+            return true;
+        }
 };
 
 

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C

@@ -0,0 +1,390 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "twoPhaseSystem.H"
+#include "surfaceInterpolate.H"
+#include "fixedValueFvsPatchFields.H"
+#include "fvcCurl.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::twoPhaseSystem::twoPhaseSystem
+(
+    const fvMesh& mesh
+)
+:
+    IOdictionary
+    (
+        IOobject
+        (
+            "phaseProperties",
+            mesh.time().constant(),
+            mesh,
+            IOobject::MUST_READ_IF_MODIFIED,
+            IOobject::NO_WRITE
+        )
+    ),
+
+    mesh_(mesh),
+
+    phase1_
+    (
+        *this,
+        *this,
+        wordList(lookup("phases"))[0]
+    ),
+
+    phase2_
+    (
+        *this,
+        *this,
+        wordList(lookup("phases"))[1]
+    ),
+
+    Cvm_
+    (
+        "Cvm",
+        dimless,
+        lookup("Cvm")
+    ),
+
+    Cl_
+    (
+        "Cl",
+        dimless,
+        lookup("Cl")
+    ),
+
+    drag1_
+    (
+        dragModel::New
+        (
+            subDict("drag"),
+            phase1_,
+            phase1_,
+            phase2_
+        )
+    ),
+
+    drag2_
+    (
+        dragModel::New
+        (
+            subDict("drag"),
+            phase2_,
+            phase2_,
+            phase1_
+        )
+    ),
+
+    heatTransfer1_
+    (
+        heatTransferModel::New
+        (
+            subDict("heatTransfer"),
+            phase1_,
+            phase1_,
+            phase2_
+        )
+    ),
+
+    heatTransfer2_
+    (
+        heatTransferModel::New
+        (
+            subDict("heatTransfer"),
+            phase2_,
+            phase2_,
+            phase1_
+        )
+    ),
+
+    dispersedPhase_(lookup("dispersedPhase")),
+
+    residualPhaseFraction_
+    (
+        readScalar(lookup("residualPhaseFraction"))
+    ),
+
+    residualSlip_
+    (
+        "residualSlip",
+        dimVelocity,
+        lookup("residualSlip")
+    )
+{
+    if
+    (
+        !(
+            dispersedPhase_ == phase1_.name()
+         || dispersedPhase_ == phase2_.name()
+         || dispersedPhase_ == "both"
+        )
+    )
+    {
+        FatalErrorIn("twoPhaseSystem::twoPhaseSystem(const fvMesh& mesh)")
+            << "invalid dispersedPhase " << dispersedPhase_
+            << exit(FatalError);
+    }
+
+    Info << "dispersedPhase is " << dispersedPhase_ << endl;
+
+    // Ensure the phase-fractions sum to 1
+    phase2_.volScalarField::operator=(scalar(1) - phase1_);
+}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField> Foam::twoPhaseSystem::rho() const
+{
+    return phase1_*phase1_.thermo().rho() + phase2_*phase2_.thermo().rho();
+}
+
+
+Foam::tmp<Foam::volVectorField> Foam::twoPhaseSystem::U() const
+{
+    return phase1_*phase1_.U() + phase2_*phase2_.U();
+}
+
+
+Foam::tmp<Foam::surfaceScalarField> Foam::twoPhaseSystem::phi() const
+{
+    return
+        fvc::interpolate(phase1_)*phase1_.phi()
+      + fvc::interpolate(phase2_)*phase2_.phi();
+}
+
+
+Foam::tmp<Foam::volScalarField> Foam::twoPhaseSystem::dragCoeff() const
+{
+    tmp<volScalarField> tdragCoeff
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "dragCoeff",
+                mesh_.time().timeName(),
+                mesh_
+            ),
+            mesh_,
+            dimensionedScalar("dragCoeff", dimensionSet(1, -3, -1, 0, 0), 0)
+        )
+    );
+    volScalarField& dragCoeff = tdragCoeff();
+
+    volVectorField Ur(phase1_.U() - phase2_.U());
+    volScalarField magUr(mag(Ur) + residualSlip_);
+
+    if (dispersedPhase_ == phase1_.name())
+    {
+        dragCoeff = drag1().K(magUr);
+    }
+    else if (dispersedPhase_ == phase2_.name())
+    {
+        dragCoeff = drag2().K(magUr);
+    }
+    else if (dispersedPhase_ == "both")
+    {
+        dragCoeff =
+        (
+            phase2_*drag1().K(magUr)
+          + phase1_*drag2().K(magUr)
+        );
+    }
+    else
+    {
+        FatalErrorIn("twoPhaseSystem::dragCoeff()")
+            << "dispersedPhase: " << dispersedPhase_ << " is incorrect"
+            << exit(FatalError);
+    }
+
+    volScalarField alphaCoeff(max(phase1_*phase2_, residualPhaseFraction_));
+    dragCoeff *= alphaCoeff;
+
+    // Remove drag at fixed-flux boundaries
+    forAll(phase1_.phi().boundaryField(), patchi)
+    {
+        if
+        (
+            isA<fixedValueFvsPatchScalarField>
+            (
+                phase1_.phi().boundaryField()[patchi]
+            )
+        )
+        {
+            dragCoeff.boundaryField()[patchi] = 0.0;
+        }
+    }
+
+    return tdragCoeff;
+}
+
+
+Foam::tmp<Foam::volVectorField> Foam::twoPhaseSystem::liftForce
+(
+    const volVectorField& U
+) const
+{
+    tmp<volVectorField> tliftForce
+    (
+        new volVectorField
+        (
+            IOobject
+            (
+                "liftForce",
+                mesh_.time().timeName(),
+                mesh_
+            ),
+            mesh_,
+            dimensionedVector
+            (
+                "liftForce",
+                dimensionSet(1, -2, -2, 0, 0),
+                vector::zero
+            )
+        )
+    );
+    volVectorField& liftForce = tliftForce();
+
+    volVectorField Ur(phase1_.U() - phase2_.U());
+
+    liftForce =
+        Cl_*(phase1_*phase1_.rho() + phase2_*phase2_.rho())
+       *(Ur ^ fvc::curl(U));
+
+    // Remove lift at fixed-flux boundaries
+    forAll(phase1_.phi().boundaryField(), patchi)
+    {
+        if
+        (
+            isA<fixedValueFvsPatchScalarField>
+            (
+                phase1_.phi().boundaryField()[patchi]
+            )
+        )
+        {
+            liftForce.boundaryField()[patchi] = vector::zero;
+        }
+    }
+
+    return tliftForce;
+}
+
+
+Foam::tmp<Foam::volScalarField> Foam::twoPhaseSystem::heatTransferCoeff() const
+{
+    tmp<volScalarField> theatTransferCoeff
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "heatTransferCoeff",
+                mesh_.time().timeName(),
+                mesh_
+            ),
+            mesh_,
+            dimensionedScalar
+            (
+                "heatTransferCoeff",
+                dimensionSet(1, -1, -3, -1, 0),
+                0
+            )
+        )
+    );
+    volScalarField& heatTransferCoeff = theatTransferCoeff();
+
+    volVectorField Ur(phase1_.U() - phase2_.U());
+    volScalarField magUr(mag(Ur) + residualSlip_);
+
+    if (dispersedPhase_ == phase1_.name())
+    {
+        heatTransferCoeff = heatTransfer1().K(magUr);
+    }
+    else if (dispersedPhase_ == phase2_.name())
+    {
+        heatTransferCoeff = heatTransfer2().K(magUr);
+    }
+    else if (dispersedPhase_ == "both")
+    {
+        heatTransferCoeff =
+        (
+            phase2_*heatTransfer1().K(magUr)
+          + phase1_*heatTransfer2().K(magUr)
+        );
+    }
+    else
+    {
+        FatalErrorIn("twoPhaseSystem::heatTransferCoeff()")
+            << "dispersedPhase: " << dispersedPhase_ << " is incorrect"
+            << exit(FatalError);
+    }
+
+    volScalarField alphaCoeff(max(phase1_*phase2_, residualPhaseFraction_));
+    heatTransferCoeff *= alphaCoeff;
+
+    // Remove heatTransfer at fixed-flux boundaries
+    forAll(phase1_.phi().boundaryField(), patchi)
+    {
+        if
+        (
+            isA<fixedValueFvsPatchScalarField>
+            (
+                phase1_.phi().boundaryField()[patchi]
+            )
+        )
+        {
+            heatTransferCoeff.boundaryField()[patchi] = 0.0;
+        }
+    }
+
+    return theatTransferCoeff;
+}
+
+
+bool Foam::twoPhaseSystem::read()
+{
+    if (regIOobject::read())
+    {
+        bool readOK = true;
+
+        readOK &= phase1_.read();
+        readOK &= phase2_.read();
+
+        lookup("Cvm") >> Cvm_;
+        lookup("Cl") >> Cl_;
+
+        return readOK;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+// ************************************************************************* //

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.H

@@ -0,0 +1,228 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::twoPhaseSystem
+
+Description
+    Incompressible multi-phase mixture with built in solution for the
+    phase fractions with interface compression for interface-capturing.
+
+    Derived from transportModel so that it can be unsed in conjunction with
+    the incompressible turbulence models.
+
+    Surface tension and contact-angle is handled for the interface
+    between each phase-pair.
+
+SourceFiles
+    twoPhaseSystem.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef twoPhaseSystem_H
+#define twoPhaseSystem_H
+
+#include "IOdictionary.H"
+#include "phaseModel.H"
+#include "dragModel.H"
+#include "heatTransferModel.H"
+#include "volFields.H"
+#include "surfaceFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward declarations
+class dragModel;
+class heatTransferModel;
+
+
+/*---------------------------------------------------------------------------*\
+                      Class twoPhaseSystem Declaration
+\*---------------------------------------------------------------------------*/
+
+class twoPhaseSystem
+:
+    public IOdictionary
+{
+    // Private data
+
+        const fvMesh& mesh_;
+
+        phaseModel phase1_;
+        phaseModel phase2_;
+
+        dimensionedScalar Cvm_;
+        dimensionedScalar Cl_;
+
+        autoPtr<dragModel> drag1_;
+        autoPtr<dragModel> drag2_;
+
+        autoPtr<heatTransferModel> heatTransfer1_;
+        autoPtr<heatTransferModel> heatTransfer2_;
+
+        word dispersedPhase_;
+        scalar residualPhaseFraction_;
+        dimensionedScalar residualSlip_;
+
+
+public:
+
+    // Constructors
+
+        //- Construct from fvMesh
+        twoPhaseSystem(const fvMesh&);
+
+
+    //- Destructor
+    virtual ~twoPhaseSystem()
+    {}
+
+
+    // Member Functions
+
+        const fvMesh& mesh() const
+        {
+            return mesh_;
+        }
+
+        const phaseModel& phase1() const
+        {
+            return phase1_;
+        }
+
+        const phaseModel& phase2() const
+        {
+            return phase2_;
+        }
+
+        const phaseModel& otherPhase(const phaseModel& phase) const
+        {
+            if (&phase == &phase1_)
+            {
+                return phase2_;
+            }
+            else
+            {
+                return phase1_;
+            }
+        }
+
+        phaseModel& phase1()
+        {
+            return phase1_;
+        }
+
+        phaseModel& phase2()
+        {
+            return phase2_;
+        }
+
+        const dragModel& drag1() const
+        {
+            return drag1_();
+        }
+
+        const dragModel& drag2() const
+        {
+            return drag2_();
+        }
+
+        const dragModel& drag(const phaseModel& phase) const
+        {
+            if (&phase == &phase1_)
+            {
+                return drag1_();
+            }
+            else
+            {
+                return drag2_();
+            }
+        }
+
+        scalar residualPhaseFraction() const
+        {
+            return residualPhaseFraction_;
+        }
+
+        const dimensionedScalar& residualSlip() const
+        {
+            return residualSlip_;
+        }
+
+        tmp<volScalarField> dragCoeff() const;
+        tmp<volVectorField> liftForce(const volVectorField& U) const;
+
+        const heatTransferModel& heatTransfer1() const
+        {
+            return heatTransfer1_();
+        }
+
+        const heatTransferModel& heatTransfer2() const
+        {
+            return heatTransfer2_();
+        }
+
+        tmp<volScalarField> heatTransferCoeff() const;
+
+        //- Return the mixture density
+        tmp<volScalarField> rho() const;
+
+        //- Return the mixture velocity
+        tmp<volVectorField> U() const;
+
+        //- Return the mixture flux
+        tmp<surfaceScalarField> phi() const;
+
+        //- Return the virtual-mass coefficient
+        dimensionedScalar Cvm() const
+        {
+            return Cvm_;
+        }
+
+        //- Return the lift coefficient
+        dimensionedScalar Cl() const
+        {
+            return Cl_;
+        }
+
+        //- Dummy correct
+        void correct()
+        {}
+
+        //- Read base phaseProperties dictionary
+        bool read();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/test/PisoFoam/Make/options

@@ -1,13 +1,13 @@
 EXE_INC = \
-    -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModel/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
     -I$(LIB_SRC)/finiteVolume/lnInclude
 
 EXE_LIBS = \
-    -lincompressibleTurbulenceModel \
+    -lturbulenceModels \
-    -lincompressibleRASModels \
+    -lincompressibleTurbulenceModels \
-    -lincompressibleLESModels \
     -lincompressibleTransportModels \
     -lfiniteVolume \
     -lmeshTools

applications/test/PisoFoam/PisoFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -33,12 +33,7 @@ Description
 
 #include "fvCFD.H"
 #include "singlePhaseTransportModel.H"
-#include "turbulenceModel.H"
+#include "IncompressibleTurbulenceModel.H"
-
-#include "LduMatrix.H"
-#include "diagTensorField.H"
-
-typedef LduMatrix<vector, scalar, scalar> lduVectorMatrix;
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -73,50 +68,12 @@ int main(int argc, char *argv[])
               + turbulence->divDevReff(U)
             );
 
-            //UEqn.relax();
+            UEqn.relax();
 
-            fvVectorMatrix UEqnp(UEqn == -fvc::grad(p));
+            if (momentumPredictor)
-
-            lduVectorMatrix U3Eqnp(mesh);
-            U3Eqnp.diag() = UEqnp.diag();
-            U3Eqnp.upper() = UEqnp.upper();
-            U3Eqnp.lower() = UEqnp.lower();
-            U3Eqnp.source() = UEqnp.source();
-
-            UEqnp.addBoundaryDiag(U3Eqnp.diag(), 0);
-            UEqnp.addBoundarySource(U3Eqnp.source(), false);
-
-            U3Eqnp.interfaces() = U.boundaryField().interfaces();
-            U3Eqnp.interfacesUpper() = UEqnp.boundaryCoeffs().component(0);
-            U3Eqnp.interfacesLower() = UEqnp.internalCoeffs().component(0);
-
-            autoPtr<lduVectorMatrix::solver> U3EqnpSolver =
-            lduVectorMatrix::solver::New
-            (
-                U.name(),
-                U3Eqnp,
-                dictionary
-                (
-                    IStringStream
-                    (
-                        "{"
-                        "    /*solver          SmoothSolver;*/"
-                        "    smoother        GaussSeidel;"
-                        "    solver           PBiCCCG;"
-                        "    preconditioner   none;"
-                        "    tolerance        (1e-7 1e-7 1);"
-                        "    relTol           (0 0 0);"
-                        "}"
-                    )()
-                )
-            );
-
-            //for (int i=0; i<3; i++)
             {
-                U3EqnpSolver->solve(U).print(Info);
+                solve(UEqn == -fvc::grad(p));
-                U.correctBoundaryConditions();
             }
-            //solve(UEqnp);
 
             // --- PISO loop
 
@@ -173,6 +130,7 @@ int main(int argc, char *argv[])
             }
         }
 
+        laminarTransport.correct();
         turbulence->correct();
 
         runTime.write();

applications/test/PisoFoam/createFields.H

@@ -26,7 +26,7 @@
         mesh
     );
 
-#   include "createPhi.H"
+    #include "createPhi.H"
 
 
     label pRefCell = 0;
@@ -36,7 +36,13 @@
 
     singlePhaseTransportModel laminarTransport(U, phi);
 
-    autoPtr<incompressible::turbulenceModel> turbulence
+    autoPtr<IncompressibleTurbulenceModel<transportModel> >
+    turbulence
     (
-        incompressible::turbulenceModel::New(U, phi, laminarTransport)
+        IncompressibleTurbulenceModel<transportModel>::New
+        (
+            U,
+            phi,
+            laminarTransport
+        )
     );

applications/test/RhoPimpleFoam/EEqn.H

@@ -0,0 +1,33 @@
+{
+    volScalarField& he = thermo.he();
+
+    fvScalarMatrix EEqn
+    (
+        fvm::ddt(rho, he) + fvm::div(phi, he)
+      + fvc::ddt(rho, K) + fvc::div(phi, K)
+      + (
+            he.name() == "e"
+          ? fvc::div
+            (
+                fvc::absolute(phi/fvc::interpolate(rho), U),
+                p,
+                "div(phiv,p)"
+            )
+          : -dpdt
+        )
+      //- fvm::laplacian(turbulence->alphaEff(), he)
+      - fvm::laplacian(turbulence->muEff(), he)
+     ==
+        fvOptions(rho, he)
+    );
+
+    EEqn.relax();
+
+    fvOptions.constrain(EEqn);
+
+    EEqn.solve();
+
+    fvOptions.correct(he);
+
+    thermo.correct();
+}

applications/test/RhoPimpleFoam/Make/files

@@ -0,0 +1,3 @@
+rhoPimpleFoam.C
+
+EXE = $(FOAM_USER_APPBIN)/RhoPimpleFoam

applications/test/RhoPimpleFoam/Make/options

@@ -0,0 +1,20 @@
+EXE_INC = \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModel/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/cfdTools \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/fvOptions/lnInclude
+
+EXE_LIBS = \
+    -lfluidThermophysicalModels \
+    -lspecie \
+    -lturbulenceModels \
+    -lcompressibleTurbulenceModels \
+    -lincompressibleTransportModels \
+    -lfiniteVolume \
+    -lmeshTools \
+    -lsampling \
+    -lfvOptions

applications/test/RhoPimpleFoam/UEqn.H

@@ -0,0 +1,22 @@
+// Solve the Momentum equation
+
+tmp<fvVectorMatrix> UEqn
+(
+    fvm::ddt(rho, U)
+  + fvm::div(phi, U)
+  + turbulence->divDevRhoReff(U)
+ ==
+    fvOptions(rho, U)
+);
+
+UEqn().relax();
+
+fvOptions.constrain(UEqn());
+
+if (pimple.momentumPredictor())
+{
+    solve(UEqn() == -fvc::grad(p));
+
+    fvOptions.correct(U);
+    K = 0.5*magSqr(U);
+}

applications/test/RhoPimpleFoam/createFields.H

@@ -0,0 +1,73 @@
+    Info<< "Reading thermophysical properties\n" << endl;
+
+    autoPtr<psiThermo> pThermo
+    (
+        psiThermo::New(mesh)
+    );
+    psiThermo& thermo = pThermo();
+    thermo.validate(args.executable(), "h", "e");
+
+    volScalarField& p = thermo.p();
+    const volScalarField& psi = thermo.psi();
+
+    volScalarField rho
+    (
+        IOobject
+        (
+            "rho",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        thermo.rho()
+    );
+
+    Info<< "Reading field U\n" << endl;
+    volVectorField U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    #include "compressibleCreatePhi.H"
+
+    dimensionedScalar rhoMax(pimple.dict().lookup("rhoMax"));
+    dimensionedScalar rhoMin(pimple.dict().lookup("rhoMin"));
+
+    Info<< "Creating turbulence model\n" << endl;
+
+    autoPtr<CompressibleTurbulenceModel<fluidThermo> >
+    turbulence
+    (
+        CompressibleTurbulenceModel<fluidThermo>::New
+        (
+            rho,
+            U,
+            phi,
+            thermo
+        )
+    );
+
+    Info<< "Creating field dpdt\n" << endl;
+    volScalarField dpdt
+    (
+        IOobject
+        (
+            "dpdt",
+            runTime.timeName(),
+            mesh
+        ),
+        mesh,
+        dimensionedScalar("dpdt", p.dimensions()/dimTime, 0)
+    );
+
+    Info<< "Creating field kinetic energy K\n" << endl;
+    volScalarField K("K", 0.5*magSqr(U));

applications/test/RhoPimpleFoam/pEqn.H

@@ -0,0 +1,113 @@
+rho = thermo.rho();
+rho = max(rho, rhoMin);
+rho = min(rho, rhoMax);
+rho.relax();
+
+volScalarField rAU(1.0/UEqn().A());
+volVectorField HbyA("HbyA", U);
+HbyA = rAU*UEqn().H();
+
+if (pimple.nCorrPISO() <= 1)
+{
+    UEqn.clear();
+}
+
+if (pimple.transonic())
+{
+    surfaceScalarField phid
+    (
+        "phid",
+        fvc::interpolate(psi)
+       *(
+            (fvc::interpolate(HbyA) & mesh.Sf())
+          + fvc::ddtPhiCorr(rAU, rho, U, phi)
+        )
+    );
+
+    fvOptions.relativeFlux(fvc::interpolate(psi), phid);
+
+    volScalarField Dp("Dp", rho*rAU);
+
+    while (pimple.correctNonOrthogonal())
+    {
+        fvScalarMatrix pEqn
+        (
+            fvm::ddt(psi, p)
+          + fvm::div(phid, p)
+          - fvm::laplacian(Dp, p)
+          ==
+            fvOptions(psi, p, rho.name())
+        );
+
+        fvOptions.constrain(pEqn);
+
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
+
+        if (pimple.finalNonOrthogonalIter())
+        {
+            phi == pEqn.flux();
+        }
+    }
+}
+else
+{
+    surfaceScalarField phiHbyA
+    (
+        "phiHbyA",
+        fvc::interpolate(rho)
+       *(
+            (fvc::interpolate(HbyA) & mesh.Sf())
+          + fvc::ddtPhiCorr(rAU, rho, U, phi)
+        )
+    );
+
+    fvOptions.relativeFlux(fvc::interpolate(rho), phiHbyA);
+
+    volScalarField Dp("Dp", rho*rAU);
+
+    while (pimple.correctNonOrthogonal())
+    {
+        // Pressure corrector
+        fvScalarMatrix pEqn
+        (
+            fvm::ddt(psi, p)
+          + fvc::div(phiHbyA)
+          - fvm::laplacian(Dp, p)
+          ==
+            fvOptions(psi, p, rho.name())
+        );
+
+        fvOptions.constrain(pEqn);
+
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
+
+        if (pimple.finalNonOrthogonalIter())
+        {
+            phi = phiHbyA + pEqn.flux();
+        }
+    }
+}
+
+#include "rhoEqn.H"
+#include "compressibleContinuityErrs.H"
+
+// Explicitly relax pressure for momentum corrector
+p.relax();
+
+// Recalculate density from the relaxed pressure
+rho = thermo.rho();
+rho = max(rho, rhoMin);
+rho = min(rho, rhoMax);
+rho.relax();
+Info<< "rho max/min : " << max(rho).value()
+    << " " << min(rho).value() << endl;
+
+U = HbyA - rAU*fvc::grad(p);
+U.correctBoundaryConditions();
+fvOptions.correct(U);
+K = 0.5*magSqr(U);
+
+if (thermo.dpdt())
+{
+    dpdt = fvc::ddt(p);
+}

applications/test/RhoPimpleFoam/rhoPimpleFoam.C

@@ -0,0 +1,107 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    rhoPimpleFoam
+
+Description
+    Transient solver for laminar or turbulent flow of compressible fluids
+    for HVAC and similar applications.
+
+    Uses the flexible PIMPLE (PISO-SIMPLE) solution for time-resolved and
+    pseudo-transient simulations.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "psiThermo.H"
+#include "CompressibleTurbulenceModel.H"
+#include "bound.H"
+#include "pimpleControl.H"
+#include "fvIOoptionList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+
+    pimpleControl pimple(mesh);
+
+    #include "createFields.H"
+    #include "createFvOptions.H"
+    #include "initContinuityErrs.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.run())
+    {
+        #include "readTimeControls.H"
+        #include "compressibleCourantNo.H"
+        #include "setDeltaT.H"
+
+        runTime++;
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        if (pimple.nCorrPIMPLE() <= 1)
+        {
+            #include "rhoEqn.H"
+        }
+
+        // --- Pressure-velocity PIMPLE corrector loop
+        while (pimple.loop())
+        {
+            #include "UEqn.H"
+            #include "EEqn.H"
+
+            // --- Pressure corrector loop
+            while (pimple.correct())
+            {
+                #include "pEqn.H"
+            }
+
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
+            }
+        }
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

src/Allwmake

@@ -64,6 +64,7 @@ wmake $makeType randomProcesses
 thermophysicalModels/Allwmake $*
 transportModels/Allwmake $*
 turbulenceModels/Allwmake $*
+TurbulenceModels/Allwmake $*
 wmake $makeType combustionModels
 regionModels/Allwmake $*
 lagrangian/Allwmake $*

src/TurbulenceModels/Allwmake

@@ -0,0 +1,12 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
+makeType=${1:-libso}
+set -x
+
+wmake libso turbulenceModel
+wmake libso incompressible
+wmake libso compressible
+wmakeLnInclude phaseIncompressible
+wmakeLnInclude phaseCompressible
+
+# ----------------------------------------------------------------- end-of-file

src/TurbulenceModels/compressible/CompressibleTurbulenceModel/CompressibleTurbulenceModel.C

@@ -0,0 +1,98 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2013 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "CompressibleTurbulenceModel.H"
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class TransportModel>
+Foam::CompressibleTurbulenceModel<TransportModel>::
+CompressibleTurbulenceModel
+(
+    const geometricOneField& alpha,
+    const volScalarField& rho,
+    const volVectorField& U,
+    const surfaceScalarField& alphaPhi,
+    const surfaceScalarField& phi,
+    const transportModel& transport,
+    const word& propertiesName
+)
+:
+    TurbulenceModel
+    <
+        geometricOneField,
+        volScalarField,
+        compressibleTurbulenceModel,
+        transportModel
+    >
+    (
+        alpha,
+        rho,
+        U,
+        alphaPhi,
+        phi,
+        transport,
+        propertiesName
+    )
+{}
+
+
+// * * * * * * * * * * * * * * * * * Selectors * * * * * * * * * * * * * * * //
+
+template<class TransportModel>
+Foam::autoPtr<Foam::CompressibleTurbulenceModel<TransportModel> >
+Foam::CompressibleTurbulenceModel<TransportModel>::New
+(
+    const volScalarField& rho,
+    const volVectorField& U,
+    const surfaceScalarField& phi,
+    const transportModel& transport,
+    const word& propertiesName
+)
+{
+    return autoPtr<CompressibleTurbulenceModel>
+    (
+        static_cast<CompressibleTurbulenceModel*>(
+        TurbulenceModel
+        <
+            geometricOneField,
+            volScalarField,
+            compressibleTurbulenceModel,
+            transportModel
+        >::New
+        (
+            geometricOneField(),
+            rho,
+            U,
+            phi,
+            phi,
+            transport,
+            propertiesName
+        ).ptr())
+    );
+}
+
+
+// ************************************************************************* //