Porting, compilation

Allwmake

@@ -1,19 +1,27 @@
 #!/bin/sh
-set -x
+cd ${0%/*} || exit 1    # run from this directory
 
-# run from this directory only
-cd ${0%/*} || exit 1
+if [ "$PWD" != "$WM_PROJECT_DIR" ]
+then
+    echo "Error: Current directory is not \$WM_PROJECT_DIR"
+    echo "    The environment variable are not consistent with the installation."
+    echo "    Check the OpenFOAM entries in your dot-files and source them."
+    exit 1
+fi
 
 # wmake is required for subsequent targets
 ( cd wmake/src && make )
 
+# build ThirdParty sources
 ( cd $WM_THIRD_PARTY_DIR && ./Allwmake )
 
-(cd src && ./Allwmake)
-
-(cd applications && ./Allwmake)
+# build OpenFOAM libraries and applications
+src/Allwmake
+applications/Allwmake
 
 if [ "$1" = doc ]
 then
-   (cd doc && ./Allwmake)
+    doc/Allwmake
 fi
+
+# ----------------------------------------------------------------- end-of-file

applications/Allwmake

@@ -1,5 +1,17 @@
 #!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
+
+if [ "$PWD" != "$WM_PROJECT_DIR/applications" ]
+then
+    echo "Error: Current directory in not \$WM_PROJECT_DIR/applications"
+    echo "    The environment variable are not consistent with the installation."
+    echo "    Check the OpenFOAM entries in your dot-files and source them."
+    exit 1
+fi
+
 set -x
 
-( cd solvers && wmake all )
-( cd utilities && wmake all )
+wmake all solvers
+wmake all utilities
+
+# ----------------------------------------------------------------- end-of-file

applications/solvers/combustion/dieselEngineFoam/dieselEngineFoam.C

@@ -77,12 +77,6 @@ int main(int argc, char *argv[])
         Info<< "Crank angle = " << runTime.theta() << " CA-deg" << endl;
 
         mesh.move();
-        Check that this is unnecessary. HJ
-        // 1.6.x merge.  HJ, 26/Aug/2010
-        const_cast<volPointInterpolation&>
-        (
-            volPointInterpolation::New(mesh)
-        ).updateMesh();
 
         dieselSpray.evolve();
 

applications/solvers/engine/sonicTurbDyMEngineFoam/Make/options

@@ -4,7 +4,7 @@ EXE_INC = \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
-    -I$(LIB_SRC)/turbulenceModels/RAS
+    -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel
 
 EXE_LIBS = \
     -lengine \
@@ -12,6 +12,7 @@ EXE_LIBS = \
     -ldynamicMesh \
     -ltopoChangerFvMesh \
     -lcompressibleRASModels \
+    -lcompressibleLESModels \
     -lbasicThermophysicalModels \
     -lspecie \
     -lmeshTools \

applications/solvers/engine/sonicTurbDyMEngineFoam/createFields.H

@@ -2,10 +2,11 @@
 
     Info<< "Reading thermophysical properties\n" << endl;
 
-    autoPtr<basicThermo> thermo
+    autoPtr<basicPsiThermo> pThermo
     (
-        basicThermo::New(mesh)
+        basicPsiThermo::New(mesh)
     );
+    basicPsiThermo& thermo = pThermo();
 
     volScalarField rho
     (
@@ -17,16 +18,16 @@
             IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
-        thermo->rho()
+        thermo.rho()
     );
     rho.oldTime();
 
-    volScalarField& p = thermo->p();
+    volScalarField& p = thermo.p();
     p.oldTime();
 
-    const volScalarField& psi = thermo->psi();
-    const volScalarField& T = thermo->T();
-    volScalarField& h = thermo->h();
+    const volScalarField& psi = thermo.psi();
+    const volScalarField& T = thermo.T();
+    volScalarField& h = thermo.h();
 
 
     Info<< "\nReading field U\n" << endl;
@@ -47,14 +48,14 @@
 
 
     Info<< "Creating turbulence model\n" << endl;
-    autoPtr<compressible::RASModel> turbulence
+    autoPtr<compressible::turbulenceModel> turbulence
     (
-        compressible::RASModel::New
+        compressible::turbulenceModel::New
         (
             rho,
             U,
             phi,
-            thermo()
+            thermo
         )
     );
 

applications/solvers/engine/sonicTurbDyMEngineFoam/hEqn.H

@@ -11,5 +11,5 @@
 
     hEqn.solve();
 
-    thermo->correct();
+    thermo.correct();
 }

applications/solvers/engine/sonicTurbDyMEngineFoam/pEqn.H

@@ -1,5 +1,5 @@
 {
-    rho = thermo->rho();
+    rho = thermo.rho();
 
     rUA = 1.0/UEqn.A();
     H = UEqn.H();
@@ -35,8 +35,8 @@
     // rho does not carry working boundary conditions and needs to be updated
     // strictly according to the thermodynamics package
     // HJ, 22/Aug/2007
-    thermo->correct();
-    rho = thermo->rho();
+    thermo.correct();
+    rho = thermo.rho();
 
     DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
 

applications/solvers/engine/sonicTurbDyMEngineFoam/sonicTurbDyMEngineFoam.C

@@ -36,8 +36,8 @@ Description
 #include "engineTime.H"
 #include "dynamicFvMesh.H"
 #include "engineTopoChangerMesh.H"
-#include "basicThermo.H"
-#include "compressible/RASModel/RASModel.H"
+#include "basicPsiThermo.H"
+#include "turbulenceModel.H"
 #include "Switch.H"
 #include "OFstream.H"
 
@@ -61,7 +61,7 @@ int main(int argc, char *argv[])
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-    thermo->correct();
+    thermo.correct();
 
     Info << "\nStarting time loop\n" << endl;
 
@@ -83,7 +83,7 @@ int main(int argc, char *argv[])
 
         if(meshChanged)
         {
-            thermo->correct();
+            thermo.correct();
 
 #           include "checkTotalVolume.H"
 #           include "compressibleCorrectPhi.H"
@@ -112,7 +112,7 @@ int main(int argc, char *argv[])
 
 #       include "logSummary.H"
 
-        rho = thermo->rho();
+        rho = thermo.rho();
 
         runTime.write();
 

applications/solvers/engine/turbDyMEngineFoam/Make/options

@@ -4,8 +4,9 @@ EXE_INC = \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/topoChangerFvMesh/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/turbulenceModels/RAS \
     -I$(LIB_SRC)/transportModels \
+    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
+    -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
     -I$(LIB_SRC)/finiteVolume/lnInclude
 
 EXE_LIBS = \
@@ -15,6 +16,7 @@ EXE_LIBS = \
     -ltopoChangerFvMesh \
     -lmeshTools \
     -lincompressibleRASModels \
+    -lincompressibleLESModels \
     -lincompressibleTransportModels \
     -lfiniteVolume \
     $(WM_DECOMP_LIBS)

applications/solvers/engine/turbDyMEngineFoam/createFields.H

@@ -39,9 +39,9 @@
 
     singlePhaseTransportModel laminarTransport(U, phi);
 
-    autoPtr<incompressible::RASModel> turbulence
+    autoPtr<incompressible::turbulenceModel> turbulence
     (
-        incompressible::RASModel::New(U, phi, laminarTransport)
+        incompressible::turbulenceModel::New(U, phi, laminarTransport)
     );
 
     Info<< "Reading field rUA if present\n" << endl;

applications/solvers/engine/turbDyMEngineFoam/turbDyMEngineFoam.C

@@ -32,8 +32,8 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "incompressible/singlePhaseTransportModel/singlePhaseTransportModel.H"
-#include "incompressible/RASModel/RASModel.H"
+#include "singlePhaseTransportModel.H"
+#include "turbulenceModel.H"
 #include "dynamicFvMesh.H"
 #include "engineTime.H"
 

applications/solvers/heatTransfer/boussinesqBuoyantFoam/boussinesqBuoyantFoam.C

@@ -43,7 +43,7 @@ int main(int argc, char *argv[])
 #   include "createTime.H"
 #   include "createMesh.H"
 #   include "readTransportProperties.H"
-#   include "readEnvironmentalProperties.H"
+#   include "readGravitationalAcceleration.H"
 #   include "createFields.H"
 #   include "initContinuityErrs.H"
 

applications/solvers/multiphase/bubbleFoam/bubbleFoam.C

@@ -40,12 +40,11 @@ Description
 
 int main(int argc, char *argv[])
 {
-
 #   include "setRootCase.H"
 
 #   include "createTime.H"
 #   include "createMesh.H"
-#   include "readEnvironmentalProperties.H"
+#   include "readGravitationalAcceleration.H"
 #   include "createFields.H"
 #   include "initContinuityErrs.H"
 
@@ -53,7 +52,7 @@ int main(int argc, char *argv[])
 
     Info<< "\nStarting time loop\n" << endl;
 
-    for (runTime++; !runTime.end(); runTime++)
+    while (runTime.loop())
     {
         Info<< "Time = " << runTime.timeName() << nl << endl;
 

applications/solvers/multiphase/cavitatingFoam/Make/files

@@ -0,0 +1,3 @@
+cavitatingFoam.C
+
+EXE = $(FOAM_APPBIN)/cavitatingFoam

applications/solvers/multiphase/cavitatingFoam/Make/options

@@ -3,12 +3,12 @@ EXE_INC = \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/transportModels/incompressible/lnInclude \
     -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
-    -I$(LIB_SRC)/turbulenceModels/RAS \
+    -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
     -I$(LIB_SRC)/thermophysicalModels/barotropicCompressibilityModel/lnInclude
 
 EXE_LIBS = \
     -lincompressibleTransportModels \
     -lincompressibleRASModels \
+    -lincompressibleLESModels \
     -lfiniteVolume \
     -lbarotropicCompressibilityModel
-

applications/solvers/multiphase/cavitatingFoam/UEqn.H

@@ -14,7 +14,11 @@
       - fvc::div(muEff*(fvc::interpolate(dev(fvc::grad(U))) & mesh.Sf()))
     );
 
+    UEqn.relax();
+
     if (momentumPredictor)
     {
         solve(UEqn == -fvc::grad(p));
     }
+
+    Info<< "max(U) " << max(mag(U)).value() << endl;

applications/solvers/multiphase/cavitatingFoam/cavitatingFoam.C

@@ -23,23 +23,25 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Application
-    rasCavitatingFoam
+    cavitatingFoam
 
 Description
-    Transient cavitation code with RAS turbulence.
+    Transient cavitation code based on the homogeneous equilibrium model
+    from which the compressibility of the liquid/vapour "mixture" is obtained.
+
+    Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
 
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
 #include "barotropicCompressibilityModel.H"
 #include "twoPhaseMixture.H"
-#include "incompressible/RASModel/RASModel.H"
+#include "turbulenceModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 int main(int argc, char *argv[])
 {
-
     #include "setRootCase.H"
 
     #include "createTime.H"
@@ -51,7 +53,7 @@ int main(int argc, char *argv[])
     #include "compressibleCourantNo.H"
     #include "setInitialDeltaT.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 
@@ -87,7 +89,7 @@ int main(int argc, char *argv[])
 
     Info<< "\n end \n";
 
-    return(0);
+    return 0;
 }
 
 

applications/solvers/multiphase/cavitatingFoam/createFields.H

@@ -78,8 +78,8 @@
 
     twoPhaseMixture twoPhaseProperties(U, phiv, "gamma");
 
-    // Create RAS turbulence model
-    autoPtr<incompressible::RASModel> turbulence
+    // Create incompressible turbulence model
+    autoPtr<incompressible::turbulenceModel> turbulence
     (
-        incompressible::RASModel::New(U, phiv, twoPhaseProperties)
+        incompressible::turbulenceModel::New(U, phiv, twoPhaseProperties)
     );

applications/solvers/multiphase/cavitatingFoam/pEqn.H

@@ -37,7 +37,14 @@
           - fvm::laplacian(rUAf, p)
         );
 
-        pEqn.solve();
+        if (corr == nCorr-1 && nonOrth == nNonOrthCorr)
+        {
+            pEqn.solve(mesh.solver(p.name() + "Final"));
+        }
+        else
+        {
+            pEqn.solve(mesh.solver(p.name()));
+        }
 
         if (nonOrth == nNonOrthCorr)
         {
@@ -45,9 +52,32 @@
         }
     }
 
-    Info<< "max-min p: " << max(p).value()
+    Info<< "Predicted p max-min : " << max(p).value()
         << " " << min(p).value() << endl;
 
+    rho == max
+    (
+        psi*p
+      + (1.0 - gamma)*rhol0
+      + ((gamma*psiv + (1.0 - gamma)*psil) - psi)*pSat,
+        rhoMin
+    );
+
+    #include "gammaPsi.H"
+
+    p =
+    (
+        rho
+      - (1.0 - gamma)*rhol0
+      - ((gamma*psiv + (1.0 - gamma)*psil) - psi)*pSat
+    )/psi;
+
+    p.correctBoundaryConditions();
+
+    Info<< "Phase-change corrected p max-min : " << max(p).value()
+        << " " << min(p).value() << endl;
+
+    // Correct velocity
 
     U = HbyA - rUA*fvc::grad(p);
 
@@ -63,18 +93,4 @@
     U.correctBoundaryConditions();
 
     Info<< "max(U) " << max(mag(U)).value() << endl;
-
-    rho == max
-    (
-        psi*p
-      + (1.0 - gamma)*rhol0
-      + ((gamma*psiv + (1.0 - gamma)*psil) - psi)*pSat,
-        rhoMin
-    );
-
-    Info<< "max-min rho: " << max(rho).value()
-        << " " << min(rho).value() << endl;
-
-#   include "gammaPsi.H"
-
 }

applications/solvers/multiphase/compressibleInterDyMFoam/Make/files

@@ -0,0 +1,3 @@
+compressibleInterDyMFoam.C
+
+EXE = $(FOAM_APPBIN)/compressibleInterDyMFoam

applications/solvers/multiphase/compressibleInterDyMFoam/Make/options

@@ -1,16 +1,20 @@
 EXE_INC = \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/transportModels/incompressible/lnInclude \
     -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
-    -I$(LIB_SRC)/turbulenceModels/LES \
-    -I$(LIB_SRC)/turbulenceModels/LES/incompressible/lnInclude \
-    -I$(LIB_SRC)/turbulenceModels/LES/LESdeltas/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/barotropicCompressibilityModel/lnInclude
+    -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/dynamicMesh/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/dynamicFvMesh/lnInclude
 
 EXE_LIBS = \
+    -linterfaceProperties \
     -lincompressibleTransportModels \
+    -lincompressibleRASModels \
     -lincompressibleLESModels \
     -lfiniteVolume \
-    -lbarotropicCompressibilityModel
+    -ldynamicMesh \
+    -lmeshTools \
+    -ldynamicFvMesh
 

applications/solvers/multiphase/compressibleInterDyMFoam/UEqn.H

@@ -2,7 +2,7 @@
     (
         "muEff",
         twoPhaseProperties.muf()
-      + fvc::interpolate(rho*turbulence->nuSgs())
+      + fvc::interpolate(rho*turbulence->nut())
     );
 
     fvVectorMatrix UEqn
@@ -11,9 +11,11 @@
       + fvm::div(rhoPhi, U)
       - fvm::laplacian(muEff, U)
       - (fvc::grad(U) & fvc::grad(muEff))
-    //- fvc::div(muEff*(fvc::interpolate(dev(fvc::grad(U))) & mesh.Sf()))
+    //- fvc::div(muf*(mesh.Sf() & fvc::interpolate(fvc::grad(U)().T())))
     );
 
+    UEqn.relax();
+
     if (momentumPredictor)
     {
         solve
@@ -22,10 +24,10 @@
          ==
             fvc::reconstruct
             (
-                (
-                    fvc::interpolate(interface.sigmaK())*fvc::snGrad(gamma)
-                  - ghf*fvc::snGrad(rho)
-                  - fvc::snGrad(pd)
+                fvc::interpolate(rho)*(g & mesh.Sf())
+              + (
+                    fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
+                  - fvc::snGrad(p)
                 ) * mesh.magSf()
             )
         );

applications/solvers/multiphase/compressibleInterDyMFoam/alphaEqns.H

@@ -59,7 +59,7 @@
                 alpharScheme
             );
 
-        MULES::explicitSolve(oneField(), alpha1, phi, phiAlpha1, Sp, Su, 1, 0);
+        MULES::explicitSolve(geometricOneField(), alpha1, phi, phiAlpha1, Sp, Su, 1, 0);
 
         surfaceScalarField rho1f = fvc::interpolate(rho1);
         surfaceScalarField rho2f = fvc::interpolate(rho2);

applications/solvers/multiphase/compressibleInterDyMFoam/alphaEqnsSubCycle.H

@@ -10,9 +10,11 @@
     );
 
     surfaceScalarField phic = mag(phi/mesh.magSf());
-    phic = min(interface.cGamma()*phic, max(phic));
+    phic = min(interface.cAlpha()*phic, max(phic));
 
+    fvc::makeAbsolute(phi, U);
     volScalarField divU = fvc::div(phi);
+    fvc::makeRelative(phi, U);
 
     if (nAlphaSubCycles > 1)
     {

applications/solvers/multiphase/compressibleInterDyMFoam/compressibleInterDyMFoam.C

@@ -0,0 +1,142 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright held by original author
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Application
+    compressibleLesInterFoam
+
+Description
+    Solver for 2 compressible, isothermal immiscible fluids using a VOF
+    (volume of fluid) phase-fraction based interface capturing approach,
+    with optional mesh motion and mesh topology changes including adaptive
+    re-meshing.
+
+    The momentum and other fluid properties are of the "mixture" and a
+    single momentum equation is solved.  Turbulence modelling is generic,
+    i.e. laminar, RAS or LES may be selected.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "dynamicFvMesh.H"
+#include "MULES.H"
+#include "subCycle.H"
+#include "interfaceProperties.H"
+#include "twoPhaseMixture.H"
+#include "turbulenceModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createDynamicFvMesh.H"
+    #include "readGravitationalAcceleration.H"
+    #include "readControls.H"
+    #include "initContinuityErrs.H"
+    #include "createFields.H"
+    #include "CourantNo.H"
+    #include "setInitialDeltaT.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.run())
+    {
+        #include "readControls.H"
+        #include "CourantNo.H"
+
+        // Make the fluxes absolute
+        fvc::makeAbsolute(phi, U);
+
+        #include "setDeltaT.H"
+
+        runTime++;
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        {
+            // Store divU from the previous mesh for the correctPhi
+            volScalarField divU = fvc::div(phi);
+
+            scalar timeBeforeMeshUpdate = runTime.elapsedCpuTime();
+
+            // Do any mesh changes
+            mesh.update();
+
+            if (mesh.changing())
+            {
+                Info<< "Execution time for mesh.update() = "
+                    << runTime.elapsedCpuTime() - timeBeforeMeshUpdate
+                    << " s" << endl;
+            }
+
+            if (mesh.changing() && correctPhi)
+            {
+                #include "correctPhi.H"
+            }
+        }
+
+        // Make the fluxes relative to the mesh motion
+        fvc::makeRelative(phi, U);
+
+        if (mesh.changing() && checkMeshCourantNo)
+        {
+            #include "meshCourantNo.H"
+        }
+
+        turbulence->correct();
+
+        // --- Outer-corrector loop
+        for (int oCorr=0; oCorr<nOuterCorr; oCorr++)
+        {
+            #include "alphaEqnsSubCycle.H"
+
+            solve(fvm::ddt(rho) + fvc::div(rhoPhi));
+
+            #include "UEqn.H"
+
+            // --- PISO loop
+            for (int corr=0; corr<nCorr; corr++)
+            {
+                #include "pEqn.H"
+            }
+        }
+
+        rho = alpha1*rho1 + alpha2*rho2;
+
+        runTime.write();
+
+        Info<< "ExecutionTime = "
+            << runTime.elapsedCpuTime()
+            << " s\n\n" << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/multiphase/compressibleInterDyMFoam/correctPhi.H

@@ -0,0 +1,61 @@
+{
+    if (mesh.changing())
+    {
+        forAll(U.boundaryField(), patchi)
+        {
+            if (U.boundaryField()[patchi].fixesValue())
+            {
+                U.boundaryField()[patchi].initEvaluate();
+            }
+        }
+
+        forAll(U.boundaryField(), patchi)
+        {
+            if (U.boundaryField()[patchi].fixesValue())
+            {
+                U.boundaryField()[patchi].evaluate();
+
+                phi.boundaryField()[patchi] =
+                U.boundaryField()[patchi] & mesh.Sf().boundaryField()[patchi];
+            }
+        }
+    }
+
+    #include "continuityErrs.H"
+
+    volScalarField pcorr
+    (
+        IOobject
+        (
+            "pcorr",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("pcorr", p.dimensions(), 0.0),
+        pcorrTypes
+    );
+
+    dimensionedScalar rAUf("(1|A(U))", dimTime/rho.dimensions(), 1.0);
+
+    adjustPhi(phi, U, pcorr);
+
+    for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    {
+        fvScalarMatrix pcorrEqn
+        (
+            fvm::laplacian(rAUf, pcorr) == fvc::div(phi) - divU
+        );
+
+        pcorrEqn.solve();
+
+        if (nonOrth == nNonOrthCorr)
+        {
+            phi -= pcorrEqn.flux();
+        }
+    }
+
+    #include "continuityErrs.H"
+}

applications/solvers/multiphase/compressibleInterDyMFoam/createFields.H

@@ -1,9 +1,9 @@
-    Info<< "Reading field pd\n" << endl;
-    volScalarField pd
+    Info<< "Reading field p\n" << endl;
+    volScalarField p
     (
         IOobject
         (
-            "pd",
+            "p",
             runTime.timeName(),
             mesh,
             IOobject::MUST_READ,
@@ -46,11 +46,6 @@
     #include "createPhi.H"
 
 
-    Info<< "Calculating field g.h\n" << endl;
-    volScalarField gh("gh", g & mesh.C());
-    surfaceScalarField ghf("ghf", g & mesh.Cf());
-
-
     Info<< "Reading transportProperties\n" << endl;
     twoPhaseMixture twoPhaseProperties(U, phi);
 
@@ -88,24 +83,6 @@
 
     dimensionedScalar pMin(twoPhaseProperties.lookup("pMin"));
 
-    volScalarField p
-    (
-        IOobject
-        (
-            "p",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        max
-        (
-            (pd + gh*(alpha1*rho10 + alpha2*rho20))
-           /(1.0 - gh*(alpha1*psi1 + alpha2*psi2)),
-            pMin
-        )
-    );
-
     volScalarField rho1 = rho10 + psi1*p;
     volScalarField rho2 = rho20 + psi2*p;
 
@@ -145,8 +122,23 @@
     // Construct interface from alpha1 distribution
     interfaceProperties interface(alpha1, U, twoPhaseProperties);
 
-    // Construct LES model
-    autoPtr<incompressible::LESModel> turbulence
+    // Construct incompressible turbulence model
+    autoPtr<incompressible::turbulenceModel> turbulence
     (
-        incompressible::LESModel::New(U, phi, twoPhaseProperties)
+        incompressible::turbulenceModel::New(U, phi, twoPhaseProperties)
     );
+
+
+    wordList pcorrTypes
+    (
+        p.boundaryField().size(),
+        zeroGradientFvPatchScalarField::typeName
+    );
+
+    for (label i=0; i<p.boundaryField().size(); i++)
+    {
+        if (p.boundaryField()[i].fixesValue())
+        {
+            pcorrTypes[i] = fixedValueFvPatchScalarField::typeName;
+        }
+    }

applications/solvers/multiphase/compressibleInterDyMFoam/pEqn.H

@@ -0,0 +1,94 @@
+{
+    volScalarField rUA = 1.0/UEqn.A();
+    surfaceScalarField rUAf = fvc::interpolate(rUA);
+
+    tmp<fvScalarMatrix> pEqnComp;
+
+    if (transonic)
+    {
+        pEqnComp =
+            (fvm::ddt(p) + fvm::div(phi, p) - fvm::Sp(fvc::div(phi), p));
+    }
+    else
+    {
+        pEqnComp =
+            (fvm::ddt(p) + fvc::div(phi, p) - fvc::Sp(fvc::div(phi), p));
+    }
+
+
+    U = rUA*UEqn.H();
+
+    surfaceScalarField phiU
+    (
+        "phiU",
+        (fvc::interpolate(U) & mesh.Sf())
+    );
+
+    phi = phiU +
+        (
+            fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)*mesh.magSf()
+          + fvc::interpolate(rho)*(g & mesh.Sf())
+        )*rUAf;
+
+    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    {
+        fvScalarMatrix pEqnIncomp
+        (
+            fvc::div(phi)
+          - fvm::laplacian(rUAf, p)
+        );
+
+        if
+        (
+            oCorr == nOuterCorr-1
+            && corr == nCorr-1
+            && nonOrth == nNonOrthCorr
+        )
+        {
+            solve
+            (
+                (
+                    max(alpha1, scalar(0))*(psi1/rho1)
+                  + max(alpha2, scalar(0))*(psi2/rho2)
+                )
+               *pEqnComp()
+              + pEqnIncomp,
+                mesh.solver(p.name() + "Final")
+            );
+        }
+        else
+        {
+            solve
+            (
+                (
+                    max(alpha1, scalar(0))*(psi1/rho1)
+                  + max(alpha2, scalar(0))*(psi2/rho2)
+                )
+               *pEqnComp()
+              + pEqnIncomp
+            );
+        }
+
+        if (nonOrth == nNonOrthCorr)
+        {
+            dgdt =
+                (pos(alpha2)*(psi2/rho2) - pos(alpha1)*(psi1/rho1))
+               *(pEqnComp & p);
+            phi += pEqnIncomp.flux();
+        }
+    }
+
+    U += rUA*fvc::reconstruct((phi - phiU)/rUAf);
+    U.correctBoundaryConditions();
+
+    p.max(pMin);
+
+    rho1 = rho10 + psi1*p;
+    rho2 = rho20 + psi2*p;
+
+    Info<< "max(U) " << max(mag(U)).value() << endl;
+    Info<< "min(p) " << min(p).value() << endl;
+
+    // Make the fluxes relative to the mesh motion
+    fvc::makeRelative(phi, U);
+}

applications/solvers/multiphase/compressibleInterDyMFoam/readControls.H

@@ -0,0 +1,32 @@
+   #include "readPISOControls.H"
+   #include "readTimeControls.H"
+
+    label nAlphaCorr
+    (
+        readLabel(piso.lookup("nAlphaCorr"))
+    );
+
+    label nAlphaSubCycles
+    (
+        readLabel(piso.lookup("nAlphaSubCycles"))
+    );
+
+    if (nAlphaSubCycles > 1 && nOuterCorr != 1)
+    {
+        FatalErrorIn(args.executable())
+            << "Sub-cycling alpha is only allowed for PISO, "
+               "i.e. when the number of outer-correctors = 1"
+            << exit(FatalError);
+    }
+
+    bool correctPhi = true;
+    if (piso.found("correctPhi"))
+    {
+        correctPhi = Switch(piso.lookup("correctPhi"));
+    }
+
+    bool checkMeshCourantNo = false;
+    if (piso.found("checkMeshCourantNo"))
+    {
+        checkMeshCourantNo = Switch(piso.lookup("checkMeshCourantNo"));
+    }

applications/solvers/multiphase/compressibleInterFoam/Make/files

@@ -0,0 +1,3 @@
+compressibleInterFoam.C
+
+EXE = $(FOAM_APPBIN)/compressibleInterFoam

applications/solvers/multiphase/compressibleInterFoam/Make/options

@@ -1,14 +1,13 @@
 EXE_INC = \
-    -I../interFoam \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/transportModels/incompressible/lnInclude \
     -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
-    -I$(LIB_SRC)/turbulenceModels/RAS \
+    -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
     -I$(LIB_SRC)/finiteVolume/lnInclude
 
 EXE_LIBS = \
     -linterfaceProperties \
     -lincompressibleTransportModels \
     -lincompressibleRASModels \
-    -lfiniteVolume \
-    -llduSolvers
+    -lincompressibleLESModels \
+    -lfiniteVolume

applications/solvers/multiphase/compressibleInterFoam/UEqn.H

@@ -24,10 +24,10 @@
          ==
             fvc::reconstruct
             (
-                (
-                    fvc::interpolate(interface.sigmaK())*fvc::snGrad(gamma)
-                  - ghf*fvc::snGrad(rho)
-                  - fvc::snGrad(pd)
+                fvc::interpolate(rho)*(g & mesh.Sf())
+              + (
+                    fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
+                  - fvc::snGrad(p)
                 ) * mesh.magSf()
             )
         );

applications/solvers/multiphase/compressibleInterFoam/alphaEqns.H

@@ -0,0 +1,76 @@
+{
+    word alphaScheme("div(phi,alpha)");
+    word alpharScheme("div(phirb,alpha)");
+
+    surfaceScalarField phir = phic*interface.nHatf();
+
+    for (int gCorr=0; gCorr<nAlphaCorr; gCorr++)
+    {
+        volScalarField::DimensionedInternalField Sp
+        (
+            IOobject
+            (
+                "Sp",
+                runTime.timeName(),
+                mesh
+            ),
+            mesh,
+            dimensionedScalar("Sp", dgdt.dimensions(), 0.0)
+        );
+
+        volScalarField::DimensionedInternalField Su
+        (
+            IOobject
+            (
+                "Su",
+                runTime.timeName(),
+                mesh
+            ),
+            // Divergence term is handled explicitly to be
+            // consistent with the explicit transport solution
+            divU*min(alpha1, scalar(1))
+        );
+
+        forAll(dgdt, celli)
+        {
+            if (dgdt[celli] > 0.0 && alpha1[celli] > 0.0)
+            {
+                Sp[celli] -= dgdt[celli]*alpha1[celli];
+                Su[celli] += dgdt[celli]*alpha1[celli];
+            }
+            else if (dgdt[celli] < 0.0 && alpha1[celli] < 1.0)
+            {
+                Sp[celli] += dgdt[celli]*(1.0 - alpha1[celli]);
+            }
+        }
+
+
+        surfaceScalarField phiAlpha1 =
+            fvc::flux
+            (
+                phi,
+                alpha1,
+                alphaScheme
+            )
+          + fvc::flux
+            (
+                -fvc::flux(-phir, alpha2, alpharScheme),
+                alpha1,
+                alpharScheme
+            );
+
+        MULES::explicitSolve(geometricOneField(), alpha1, phi, phiAlpha1, Sp, Su, 1, 0);
+
+        surfaceScalarField rho1f = fvc::interpolate(rho1);
+        surfaceScalarField rho2f = fvc::interpolate(rho2);
+        rhoPhi = phiAlpha1*(rho1f - rho2f) + phi*rho2f;
+
+        alpha2 = scalar(1) - alpha1;
+    }
+
+    Info<< "Liquid phase volume fraction = "
+        << alpha1.weightedAverage(mesh.V()).value()
+        << "  Min(alpha1) = " << min(alpha1).value()
+        << "  Min(alpha2) = " << min(alpha2).value()
+        << endl;
+}

applications/solvers/multiphase/compressibleInterFoam/alphaEqnsSubCycle.H

@@ -0,0 +1,43 @@
+{
+    label nAlphaCorr
+    (
+        readLabel(piso.lookup("nAlphaCorr"))
+    );
+
+    label nAlphaSubCycles
+    (
+        readLabel(piso.lookup("nAlphaSubCycles"))
+    );
+
+    surfaceScalarField phic = mag(phi/mesh.magSf());
+    phic = min(interface.cAlpha()*phic, max(phic));
+
+    volScalarField divU = fvc::div(phi);
+
+    if (nAlphaSubCycles > 1)
+    {
+        dimensionedScalar totalDeltaT = runTime.deltaT();
+        surfaceScalarField rhoPhiSum = 0.0*rhoPhi;
+
+        for
+        (
+            subCycle<volScalarField> alphaSubCycle(alpha1, nAlphaSubCycles);
+            !(++alphaSubCycle).end();
+        )
+        {
+            #include "alphaEqns.H"
+            rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi;
+        }
+
+        rhoPhi = rhoPhiSum;
+    }
+    else
+    {
+        #include "alphaEqns.H"
+    }
+
+    if (oCorr == 0)
+    {
+        interface.correct();
+    }
+}

applications/solvers/multiphase/compressibleInterFoam/compressibleInterFoam.C

@@ -23,14 +23,16 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Application
-    compressibleLesInterFoam
+    compressibleInterFoam
 
 Description
     Solver for 2 compressible, isothermal immiscible fluids using a VOF
     (volume of fluid) phase-fraction based interface capturing approach.
+
     The momentum and other fluid properties are of the "mixture" and a single
-    momentum equation is solved.  Turbulence is modelled using a run-time
-    selectable incompressible LES model.
+    momentum equation is solved.
+
+    Turbulence modelling is generic, i.e.  laminar, RAS or LES may be selected.
 
 \*---------------------------------------------------------------------------*/
 
@@ -39,7 +41,7 @@ Description
 #include "subCycle.H"
 #include "interfaceProperties.H"
 #include "twoPhaseMixture.H"
-#include "incompressible/LESModel/LESModel.H"
+#include "turbulenceModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -48,7 +50,7 @@ int main(int argc, char *argv[])
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
-    #include "readEnvironmentalProperties.H"
+    #include "readGravitationalAcceleration.H"
     #include "readControls.H"
     #include "initContinuityErrs.H"
     #include "createFields.H"
@@ -69,8 +71,6 @@ int main(int argc, char *argv[])
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        turbulence->correct();
-
         // --- Outer-corrector loop
         for (int oCorr=0; oCorr<nOuterCorr; oCorr++)
         {
@@ -89,6 +89,8 @@ int main(int argc, char *argv[])
 
         rho = alpha1*rho1 + alpha2*rho2;
 
+        turbulence->correct();
+
         runTime.write();
 
         Info<< "ExecutionTime = "
@@ -98,7 +100,7 @@ int main(int argc, char *argv[])
 
     Info<< "End\n" << endl;
 
-    return(0);
+    return 0;
 }
 
 

applications/solvers/multiphase/compressibleInterFoam/createFields.H

@@ -0,0 +1,129 @@
+    Info<< "Reading field p\n" << endl;
+    volScalarField p
+    (
+        IOobject
+        (
+            "p",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< "Reading field alpha1\n" << endl;
+    volScalarField alpha1
+    (
+        IOobject
+        (
+            "alpha1",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< "Calculating field alpha1\n" << endl;
+    volScalarField alpha2("alpha2", scalar(1) - alpha1);
+
+    Info<< "Reading field U\n" << endl;
+    volVectorField U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    #include "createPhi.H"
+
+
+    Info<< "Reading transportProperties\n" << endl;
+    twoPhaseMixture twoPhaseProperties(U, phi);
+
+    dimensionedScalar rho10
+    (
+        twoPhaseProperties.subDict
+        (
+            twoPhaseProperties.phase1Name()
+        ).lookup("rho0")
+    );
+
+    dimensionedScalar rho20
+    (
+        twoPhaseProperties.subDict
+        (
+            twoPhaseProperties.phase2Name()
+        ).lookup("rho0")
+    );
+
+    dimensionedScalar psi1
+    (
+        twoPhaseProperties.subDict
+        (
+            twoPhaseProperties.phase1Name()
+        ).lookup("psi")
+    );
+
+    dimensionedScalar psi2
+    (
+        twoPhaseProperties.subDict
+        (
+            twoPhaseProperties.phase2Name()
+        ).lookup("psi")
+    );
+
+    dimensionedScalar pMin(twoPhaseProperties.lookup("pMin"));
+
+    volScalarField rho1 = rho10 + psi1*p;
+    volScalarField rho2 = rho20 + psi2*p;
+
+    volScalarField rho
+    (
+        IOobject
+        (
+            "rho",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        alpha1*rho1 + alpha2*rho2
+    );
+
+
+    // Mass flux
+    // Initialisation does not matter because rhoPhi is reset after the
+    // alpha1 solution before it is used in the U equation.
+    surfaceScalarField rhoPhi
+    (
+        IOobject
+        (
+            "rho*phi",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        fvc::interpolate(rho)*phi
+    );
+
+    volScalarField dgdt =
+        pos(alpha2)*fvc::div(phi)/max(alpha2, scalar(0.0001));
+
+    // Construct interface from alpha1 distribution
+    interfaceProperties interface(alpha1, U, twoPhaseProperties);
+
+    // Construct incompressible turbulence model
+    autoPtr<incompressible::turbulenceModel> turbulence
+    (
+        incompressible::turbulenceModel::New(U, phi, twoPhaseProperties)
+    );

applications/solvers/multiphase/compressibleInterFoam/pEqn.H

@@ -2,17 +2,17 @@
     volScalarField rUA = 1.0/UEqn.A();
     surfaceScalarField rUAf = fvc::interpolate(rUA);
 
-    tmp<fvScalarMatrix> pdEqnComp;
+    tmp<fvScalarMatrix> pEqnComp;
 
     if (transonic)
     {
-        pdEqnComp =
-            (fvm::ddt(pd) + fvm::div(phi, pd) - fvm::Sp(fvc::div(phi), pd));
+        pEqnComp =
+            (fvm::ddt(p) + fvm::div(phi, p) - fvm::Sp(fvc::div(phi), p));
     }
     else
     {
-        pdEqnComp =
-            (fvm::ddt(pd) + fvc::div(phi, pd) - fvc::Sp(fvc::div(phi), pd));
+        pEqnComp =
+            (fvm::ddt(p) + fvc::div(phi, p) - fvc::Sp(fvc::div(phi), p));
     }
 
 
@@ -26,16 +26,16 @@
 
     phi = phiU +
         (
-            fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
-          - ghf*fvc::snGrad(rho)
-        )*rUAf*mesh.magSf();
+            fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)*mesh.magSf()
+          + fvc::interpolate(rho)*(g & mesh.Sf())
+        )*rUAf;
 
     for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
-        fvScalarMatrix pdEqnIncomp
+        fvScalarMatrix pEqnIncomp
         (
             fvc::div(phi)
-          - fvm::laplacian(rUAf, pd)
+          - fvm::laplacian(rUAf, p)
         );
 
         solve
@@ -44,31 +44,27 @@
                 max(alpha1, scalar(0))*(psi1/rho1)
               + max(alpha2, scalar(0))*(psi2/rho2)
             )
-           *pdEqnComp()
-          + pdEqnIncomp
+           *pEqnComp()
+          + pEqnIncomp
         );
 
         if (nonOrth == nNonOrthCorr)
         {
             dgdt =
                 (pos(alpha2)*(psi2/rho2) - pos(alpha1)*(psi1/rho1))
-               *(pdEqnComp & pd);
-            phi += pdEqnIncomp.flux();
+               *(pEqnComp & p);
+            phi += pEqnIncomp.flux();
         }
     }
 
     U += rUA*fvc::reconstruct((phi - phiU)/rUAf);
     U.correctBoundaryConditions();
 
-    p = max
-        (
-            (pd + gh*(alpha1*rho10 + alpha2*rho20))/(1.0 - gh*(alpha1*psi1 + alpha2*psi2)),
-            pMin
-        );
+    p.max(pMin);
 
     rho1 = rho10 + psi1*p;
     rho2 = rho20 + psi2*p;
 
     Info<< "max(U) " << max(mag(U)).value() << endl;
-    Info<< "min(pd) " << min(pd).value() << endl;
+    Info<< "min(p) " << min(p).value() << endl;
 }

applications/solvers/multiphase/compressibleLesInterFoam/Make/files

@@ -1,3 +0,0 @@
-compressibleLesInterFoam.C
-
-EXE = $(FOAM_APPBIN)/compressibleLesInterFoam

applications/solvers/multiphase/compressibleLesInterFoam/UEqn.H

@@ -1,29 +0,0 @@
-    surfaceScalarField muf =
-        twoPhaseProperties.muf()
-      + fvc::interpolate(rho*turbulence->nuSgs());
-
-    fvVectorMatrix UEqn
-    (
-        fvm::ddt(rho, U)
-      + fvm::div(rhoPhi, U)
-      - fvm::laplacian(muf, U)
-      - (fvc::grad(U) & fvc::grad(muf))
-    //- fvc::div(muf*(mesh.Sf() & fvc::interpolate(fvc::grad(U)().T())))
-    );
-
-    if (momentumPredictor)
-    {
-        solve
-        (
-            UEqn
-         ==
-            fvc::reconstruct
-            (
-                (
-                    fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
-                  - ghf*fvc::snGrad(rho)
-                  - fvc::snGrad(pd)
-                ) * mesh.magSf()
-            )
-        );
-    }

applications/solvers/multiphase/interDyMFoam/Make/options

@@ -1,9 +1,9 @@
 EXE_INC = \
+    -I../interFoam \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/transportModels/incompressible/lnInclude \
     -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
-    -I$(LIB_SRC)/turbulenceModels/RAS \
-    -I$(LIB_SRC)/turbulenceModels/RAS/incompressible/lnInclude \
+    -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
@@ -13,8 +13,10 @@ EXE_LIBS = \
     -linterfaceProperties \
     -lincompressibleTransportModels \
     -lincompressibleRASModels \
+    -lincompressibleLESModels \
     -lfiniteVolume \
     -ldynamicMesh \
     -lmeshTools \
-    -ldynamicFvMesh
-
+    -ldynamicFvMesh \
+    -ltopoChangerFvMesh \
+    -llduSolvers

applications/solvers/multiphase/interDyMFoam/UEqn.H

@@ -1,29 +0,0 @@
-    surfaceScalarField muf = twoPhaseProperties.muf();
-
-    fvVectorMatrix UEqn
-    (
-        fvm::ddt(rho, U)
-      + fvm::div(rhoPhi, U)
-      - fvm::laplacian(muf, U)
-      - (fvc::grad(U) & fvc::grad(muf))
-    //- fvc::div(muf*(fvc::interpolate(dev(fvc::grad(U))) & mesh.Sf()))
-    );
-
-    UEqn.relax();
-
-    if (momentumPredictor)
-    {
-        solve
-        (
-            UEqn
-         ==
-            fvc::reconstruct
-            (
-                (
-                    fvc::interpolate(interface.sigmaK())*fvc::snGrad(gamma)
-                  - ghf*fvc::snGrad(rho)
-                  - fvc::snGrad(pd)
-                ) * mesh.magSf()
-            )
-        );
-    }

applications/solvers/multiphase/interDyMFoam/correctPhi.H

@@ -1,4 +1,26 @@
 {
+    if (mesh.changing())
+    {
+        forAll(U.boundaryField(), patchi)
+        {
+            if (U.boundaryField()[patchi].fixesValue())
+            {
+                U.boundaryField()[patchi].initEvaluate();
+            }
+        }
+
+        forAll(U.boundaryField(), patchi)
+        {
+            if (U.boundaryField()[patchi].fixesValue())
+            {
+                U.boundaryField()[patchi].evaluate();
+
+                phi.boundaryField()[patchi] =
+                U.boundaryField()[patchi] & mesh.Sf().boundaryField()[patchi];
+            }
+        }
+    }
+
 #   include "continuityErrs.H"
 
     volScalarField pcorr

applications/solvers/multiphase/interDyMFoam/createFields.H

@@ -12,12 +12,12 @@
         mesh
     );
 
-    Info<< "Reading field gamma\n" << endl;
-    volScalarField gamma
+    Info<< "Reading field alpha1\n" << endl;
+    volScalarField alpha1
     (
         IOobject
         (
-            "gamma",
+            "alpha1",
             runTime.timeName(),
             mesh,
             IOobject::MUST_READ,
@@ -44,7 +44,7 @@
 
 
     Info<< "Reading transportProperties\n" << endl;
-    twoPhaseMixture twoPhaseProperties(U, phi, "gamma");
+    twoPhaseMixture twoPhaseProperties(U, phi, "alpha1");
 
     const dimensionedScalar& rho1 = twoPhaseProperties.rho1();
     const dimensionedScalar& rho2 = twoPhaseProperties.rho2();
@@ -60,15 +60,15 @@
             mesh,
             IOobject::READ_IF_PRESENT
         ),
-        gamma*rho1 + (scalar(1) - gamma)*rho2,
-        gamma.boundaryField().types()
+        alpha1*rho1 + (scalar(1) - alpha1)*rho2,
+        alpha1.boundaryField().types()
     );
     rho.oldTime();
 
 
     // Mass flux
     // Initialisation does not matter because rhoPhi is reset after the
-    // gamma solution before it is used in the U equation.
+    // alpha1 solution before it is used in the U equation.
     surfaceScalarField rhoPhi
     (
         IOobject
@@ -83,16 +83,20 @@
     );
 
 
-    // Construct interface from gamma distribution
-    interfaceProperties interface(gamma, U, twoPhaseProperties);
+    // Construct interface from alpha1 distribution
+    interfaceProperties interface(alpha1, U, twoPhaseProperties);
 
-    // Construct incompressible RAS model
-    autoPtr<incompressible::RASModel> turbulence
+    // Construct incompressible turbulence model
+    autoPtr<incompressible::turbulenceModel> turbulence
     (
-        incompressible::RASModel::New(U, phi, twoPhaseProperties)
+        incompressible::turbulenceModel::New(U, phi, twoPhaseProperties)
     );
 
-    wordList pcorrTypes(pd.boundaryField().types());
+    wordList pcorrTypes
+    (
+        pd.boundaryField().size(),
+        zeroGradientFvPatchScalarField::typeName
+    );
 
     for (label i = 0; i < pd.boundaryField().size(); i++)
     {

applications/solvers/multiphase/interDyMFoam/gammaEqn.H

@@ -1,35 +0,0 @@
-{
-    word gammaScheme("div(phi,gamma)");
-    word gammarScheme("div(phirb,gamma)");
-
-    surfaceScalarField phic = mag(phi/mesh.magSf());
-    phic = min(interface.cGamma()*phic, max(phic));
-    surfaceScalarField phir = phic*interface.nHatf();
-
-    for (int gCorr=0; gCorr<nGammaCorr; gCorr++)
-    {
-        surfaceScalarField phiGamma = 
-            fvc::flux
-            (
-                phi,
-                gamma,
-                gammaScheme
-            )
-          + fvc::flux
-            (
-                -fvc::flux(-phir, scalar(1) - gamma, gammarScheme),
-                gamma,
-                gammarScheme
-            );
-
-        MULES::explicitSolve(gamma, phi, phiGamma, 1, 0);
-
-        rhoPhi = phiGamma*(rho1 - rho2) + phi*rho2;
-    }
-
-    Info<< "Liquid phase volume fraction = "
-        << gamma.weightedAverage(mesh.V()).value()
-        << "  Min(gamma) = " << min(gamma).value()
-        << "  Max(gamma) = " << max(gamma).value()
-        << endl;
-}

applications/solvers/multiphase/interDyMFoam/gammaEqnSubCycle.H

@@ -1,35 +0,0 @@
-label nGammaCorr
-(
-    readLabel(piso.lookup("nGammaCorr"))
-);
-
-label nGammaSubCycles
-(
-    readLabel(piso.lookup("nGammaSubCycles"))
-);
-
-if (nGammaSubCycles > 1)
-{
-    dimensionedScalar totalDeltaT = runTime.deltaT();
-    surfaceScalarField rhoPhiSum = 0.0*rhoPhi;
-
-    for
-    (
-        subCycle<volScalarField> gammaSubCycle(gamma, nGammaSubCycles);
-        !(++gammaSubCycle).end();
-    )
-    {
-#       include "gammaEqn.H"
-        rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi;
-    }
-
-    rhoPhi = rhoPhiSum;
-}
-else
-{
-#       include "gammaEqn.H"
-}
-
-interface.correct();
-
-rho == gamma*rho1 + (scalar(1) - gamma)*rho2;

applications/solvers/multiphase/interDyMFoam/interDyMFoam.C

@@ -39,8 +39,7 @@ Description
 #include "subCycle.H"
 #include "interfaceProperties.H"
 #include "twoPhaseMixture.H"
-#include "incompressible/RASModel/RASModel.H"
-#include "EulerDdtScheme.H"
+#include "turbulenceModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -49,7 +48,7 @@ int main(int argc, char *argv[])
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createDynamicFvMesh.H"
-    #include "readEnvironmentalProperties.H"
+    #include "readGravitationalAcceleration.H"
     #include "readPISOControls.H"
     #include "initContinuityErrs.H"
     #include "createFields.H"
@@ -106,7 +105,7 @@ int main(int argc, char *argv[])
 
         twoPhaseProperties.correct();
 
-        #include "gammaEqnSubCycle.H"
+        #include "alphaEqnSubCycle.H"
 
         #include "UEqn.H"
 
@@ -139,7 +138,7 @@ int main(int argc, char *argv[])
 
     Info<< "End\n" << endl;
 
-    return(0);
+    return 0;
 }
 
 

applications/solvers/multiphase/interDyMFoam/pEqn.H

@@ -5,17 +5,18 @@
     U = rAU*UEqn.H();
     surfaceScalarField phiU("phiU", (fvc::interpolate(U) & mesh.Sf()));
 
-    phi = phiU +
-    (
-        fvc::interpolate(interface.sigmaK())*fvc::snGrad(gamma)
-      - ghf*fvc::snGrad(rho)
-    )*rAUf*mesh.magSf();
-
     if (pd.needReference())
     {
         adjustPhi(phi, U, pd);
     }
 
+    phi = phiU +
+    (
+        fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
+      - ghf*fvc::snGrad(rho)
+    )*rAUf*mesh.magSf();
+
+
     for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {
         fvScalarMatrix pdEqn

applications/solvers/multiphase/interFoam/Make/options

@@ -2,10 +2,13 @@ EXE_INC = \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/transportModels/incompressible/lnInclude \
     -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
+    -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
     -I$(LIB_SRC)/finiteVolume/lnInclude
 
 EXE_LIBS = \
     -linterfaceProperties \
     -lincompressibleTransportModels \
+    -lincompressibleRASModels \
+    -lincompressibleLESModels \
     -lfiniteVolume \
     -llduSolvers

applications/solvers/multiphase/interFoam/UBlendingFactor.H

@@ -1,6 +0,0 @@
-    surfaceScalarField gammaf = fvc::interpolate(gamma);
-    surfaceScalarField UBlendingFactor
-    (
-        "UBlendingFactor",
-        sqrt(max(min(4*gammaf*(1.0 - gammaf), 1.0), 0.0))
-    );

applications/solvers/multiphase/interFoam/UEqn.H

@@ -1,12 +1,17 @@
-    surfaceScalarField muf = twoPhaseProperties.muf();
+    surfaceScalarField muEff
+    (
+        "muEff",
+        twoPhaseProperties.muf()
+      + fvc::interpolate(rho*turbulence->nut())
+    );
 
     fvVectorMatrix UEqn
     (
         fvm::ddt(rho, U)
       + fvm::div(rhoPhi, U)
-      - fvm::laplacian(muf, U)
-      - (fvc::grad(U) & fvc::grad(muf))
-    //- fvc::div(muf*(fvc::interpolate(dev(fvc::grad(U))) & mesh.Sf()))
+      - fvm::laplacian(muEff, U)
+      - (fvc::grad(U) & fvc::grad(muEff))
+    //- fvc::div(muEff*(fvc::interpolate(dev(fvc::grad(U))) & mesh.Sf()))
     );
 
     UEqn.relax();
@@ -20,7 +25,7 @@
             fvc::reconstruct
             (
                 (
-                    fvc::interpolate(interface.sigmaK())*fvc::snGrad(gamma)
+                    fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
                   - ghf*fvc::snGrad(rho)
                   - fvc::snGrad(pd)
                 ) * mesh.magSf()

applications/solvers/multiphase/interFoam/alphaEqn.H

@@ -0,0 +1,35 @@
+{
+    word alphaScheme("div(phi,alpha)");
+    word alpharScheme("div(phirb,alpha)");
+
+    surfaceScalarField phic = mag(phi/mesh.magSf());
+    phic = min(interface.cAlpha()*phic, max(phic));
+    surfaceScalarField phir = phic*interface.nHatf();
+
+    for (int aCorr=0; aCorr<nAlphaCorr; aCorr++)
+    {
+        surfaceScalarField phiAlpha =
+            fvc::flux
+            (
+                phi,
+                alpha1,
+                alphaScheme
+            )
+          + fvc::flux
+            (
+                -fvc::flux(-phir, scalar(1) - alpha1, alpharScheme),
+                alpha1,
+                alpharScheme
+            );
+
+        MULES::explicitSolve(alpha1, phi, phiAlpha, 1, 0);
+
+        rhoPhi = phiAlpha*(rho1 - rho2) + phi*rho2;
+    }
+
+    Info<< "Liquid phase volume fraction = "
+        << alpha1.weightedAverage(mesh.V()).value()
+        << "  Min(alpha1) = " << min(alpha1).value()
+        << "  Max(alpha1) = " << max(alpha1).value()
+        << endl;
+}

applications/solvers/multiphase/interFoam/alphaEqnSubCycle.H

@@ -0,0 +1,35 @@
+label nAlphaCorr
+(
+    readLabel(piso.lookup("nAlphaCorr"))
+);
+
+label nAlphaSubCycles
+(
+    readLabel(piso.lookup("nAlphaSubCycles"))
+);
+
+if (nAlphaSubCycles > 1)
+{
+    dimensionedScalar totalDeltaT = runTime.deltaT();
+    surfaceScalarField rhoPhiSum = 0.0*rhoPhi;
+
+    for
+    (
+        subCycle<volScalarField> alphaSubCycle(alpha1, nAlphaSubCycles);
+        !(++alphaSubCycle).end();
+    )
+    {
+#       include "alphaEqn.H"
+        rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi;
+    }
+
+    rhoPhi = rhoPhiSum;
+}
+else
+{
+#       include "alphaEqn.H"
+}
+
+interface.correct();
+
+rho == alpha1*rho1 + (scalar(1) - alpha1)*rho2;

applications/solvers/multiphase/interFoam/correctPhi.H

@@ -1,7 +1,11 @@
 {
 #   include "continuityErrs.H"
 
-    wordList pcorrTypes(pd.boundaryField().types());
+    wordList pcorrTypes
+    (
+        pd.boundaryField().size(),
+        zeroGradientFvPatchScalarField::typeName
+    );
 
     for (label i=0; i<pd.boundaryField().size(); i++)
     {

applications/solvers/multiphase/interFoam/createFields.H

@@ -12,12 +12,12 @@
         mesh
     );
 
-    Info<< "Reading field gamma\n" << endl;
-    volScalarField gamma
+    Info<< "Reading field alpha1\n" << endl;
+    volScalarField alpha1
     (
         IOobject
         (
-            "gamma",
+            "alpha1",
             runTime.timeName(),
             mesh,
             IOobject::MUST_READ,
@@ -44,7 +44,7 @@
 
 
     Info<< "Reading transportProperties\n" << endl;
-    twoPhaseMixture twoPhaseProperties(U, phi, "gamma");
+    twoPhaseMixture twoPhaseProperties(U, phi, "alpha1");
 
     const dimensionedScalar& rho1 = twoPhaseProperties.rho1();
     const dimensionedScalar& rho2 = twoPhaseProperties.rho2();
@@ -60,15 +60,15 @@
             mesh,
             IOobject::READ_IF_PRESENT
         ),
-        gamma*rho1 + (scalar(1) - gamma)*rho2,
-        gamma.boundaryField().types()
+        alpha1*rho1 + (scalar(1) - alpha1)*rho2,
+        alpha1.boundaryField().types()
     );
     rho.oldTime();
 
 
     // Mass flux
     // Initialisation does not matter because rhoPhi is reset after the
-    // gamma solution before it is used in the U equation.
+    // alpha1 solution before it is used in the U equation.
     surfaceScalarField rhoPhi
     (
         IOobject
@@ -87,7 +87,6 @@
     volScalarField gh("gh", g & mesh.C());
     surfaceScalarField ghf("gh", g & mesh.Cf());
 
-
     volScalarField p
     (
         IOobject
@@ -123,5 +122,11 @@
         );
     }
 
-    // Construct interface from gamma distribution
-    interfaceProperties interface(gamma, U, twoPhaseProperties);
+    // Construct interface from alpha1 distribution
+    interfaceProperties interface(alpha1, U, twoPhaseProperties);
+
+     // Construct incompressible turbulence model
+    autoPtr<incompressible::turbulenceModel> turbulence
+    (
+        incompressible::turbulenceModel::New(U, phi, twoPhaseProperties)
+    );

applications/solvers/multiphase/interFoam/gammaEqn.H

@@ -1,35 +0,0 @@
-{
-    word gammaScheme("div(phi,gamma)");
-    word gammarScheme("div(phirb,gamma)");
-
-    surfaceScalarField phic = mag(phi/mesh.magSf());
-    phic = min(interface.cGamma()*phic, max(phic));
-    surfaceScalarField phir = phic*interface.nHatf();
-
-    for (int gCorr=0; gCorr<nGammaCorr; gCorr++)
-    {
-        surfaceScalarField phiGamma = 
-            fvc::flux
-            (
-                phi,
-                gamma,
-                gammaScheme
-            )
-          + fvc::flux
-            (
-                -fvc::flux(-phir, scalar(1) - gamma, gammarScheme),
-                gamma,
-                gammarScheme
-            );
-
-        MULES::explicitSolve(gamma, phi, phiGamma, 1, 0);
-
-        rhoPhi = phiGamma*(rho1 - rho2) + phi*rho2;
-    }
-
-    Info<< "Liquid phase volume fraction = "
-        << gamma.weightedAverage(mesh.V()).value()
-        << "  Min(gamma) = " << min(gamma).value()
-        << "  Max(gamma) = " << max(gamma).value()
-        << endl;
-}

applications/solvers/multiphase/interFoam/gammaEqnSubCycle.H

@@ -1,35 +0,0 @@
-label nGammaCorr
-(
-    readLabel(piso.lookup("nGammaCorr"))
-);
-
-label nGammaSubCycles
-(
-    readLabel(piso.lookup("nGammaSubCycles"))
-);
-
-if (nGammaSubCycles > 1)
-{
-    dimensionedScalar totalDeltaT = runTime.deltaT();
-    surfaceScalarField rhoPhiSum = 0.0*rhoPhi;
-
-    for
-    (
-        subCycle<volScalarField> gammaSubCycle(gamma, nGammaSubCycles);
-        !(++gammaSubCycle).end();
-    )
-    {
-#       include "gammaEqn.H"
-        rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi;
-    }
-
-    rhoPhi = rhoPhiSum;
-}
-else
-{
-#       include "gammaEqn.H"
-}
-
-interface.correct();
-
-rho == gamma*rho1 + (scalar(1) - gamma)*rho2;

applications/solvers/multiphase/interFoam/interFoam.C

@@ -28,9 +28,12 @@ Application
 Description
     Solver for 2 incompressible, isothermal immiscible fluids using a VOF
     (volume of fluid) phase-fraction based interface capturing approach.
+
     The momentum and other fluid properties are of the "mixture" and a single
     momentum equation is solved.
 
+    Turbulence modelling is generic, i.e.  laminar, RAS or LES may be selected.
+
     For a two-fluid approach see twoPhaseEulerFoam.
 
 \*---------------------------------------------------------------------------*/
@@ -40,6 +43,7 @@ Description
 #include "subCycle.H"
 #include "interfaceProperties.H"
 #include "twoPhaseMixture.H"
+#include "turbulenceModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -48,7 +52,7 @@ int main(int argc, char *argv[])
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
-    #include "readEnvironmentalProperties.H"
+    #include "readGravitationalAcceleration.H"
     #include "readPISOControls.H"
     #include "initContinuityErrs.H"
     #include "createFields.H"
@@ -74,7 +78,7 @@ int main(int argc, char *argv[])
 
         twoPhaseProperties.correct();
 
-        #include "gammaEqnSubCycle.H"
+        #include "alphaEqnSubCycle.H"
 
         #include "UEqn.H"
 
@@ -107,7 +111,7 @@ int main(int argc, char *argv[])
 
     Info<< "End\n" << endl;
 
-    return(0);
+    return 0;
 }
 
 

applications/solvers/multiphase/interFoam/pEqn.H

@@ -7,16 +7,18 @@
     surfaceScalarField phiU
     (
         "phiU",
-        (fvc::interpolate(U) & mesh.Sf()) + fvc::ddtPhiCorr(rUA, rho, U, phi)
+        (fvc::interpolate(U) & mesh.Sf())
+      + fvc::ddtPhiCorr(rUA, rho, U, phi)
     );
 
+    adjustPhi(phiU, U, p);
+
     phi = phiU +
         (
-            fvc::interpolate(interface.sigmaK())*fvc::snGrad(gamma)
+            fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
           - ghf*fvc::snGrad(rho)
         )*rUAf*mesh.magSf();
 
-    adjustPhi(phi, U, pd);
 
     for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
     {

applications/solvers/multiphase/interMixingFoam/Make/files

@@ -0,0 +1,6 @@
+incompressibleThreePhaseMixture/threePhaseMixture.C
+threePhaseInterfaceProperties/threePhaseInterfaceProperties.C
+interMixingFoam.C
+
+EXE = $(FOAM_APPBIN)/interMixingFoam
+

applications/solvers/multiphase/interMixingFoam/Make/options

@@ -0,0 +1,17 @@
+INTERFOAM = $(FOAM_SOLVERS)/multiphase/interFoam
+
+EXE_INC = \
+    -I$(INTERFOAM) \
+    -IincompressibleThreePhaseMixture \
+    -IthreePhaseInterfaceProperties \
+    -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
+    -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/transportModels
+
+EXE_LIBS = \
+    -linterfaceProperties \
+    -lincompressibleTransportModels \
+    -lincompressibleRASModels \
+    -lincompressibleLESModels \
+    -lfiniteVolume

applications/solvers/multiphase/interMixingFoam/alphaEqns.H

@@ -0,0 +1,164 @@
+{
+    word alphaScheme("div(phi,alpha)");
+    word alpharScheme("div(phirb,alpha)");
+
+    surfaceScalarField phir
+    (
+        IOobject
+        (
+            "phir",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        interface.cAlpha()*mag(phi/mesh.magSf())*interface.nHatf()
+    );
+
+    for (int gCorr=0; gCorr<nAlphaCorr; gCorr++)
+    {
+        // Create the limiter to be used for all phase-fractions
+        scalarField allLambda(mesh.nFaces(), 1.0);
+
+        // Split the limiter into a surfaceScalarField
+        slicedSurfaceScalarField lambda
+        (
+            IOobject
+            (
+                "lambda",
+                mesh.time().timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            mesh,
+            dimless,
+            allLambda
+        );
+
+
+        // Create the complete convection flux for alpha1
+        surfaceScalarField phiAlpha1 =
+            fvc::flux
+            (
+                phi,
+                alpha1,
+                alphaScheme
+            )
+          + fvc::flux
+            (
+                -fvc::flux(-phir, alpha2, alpharScheme),
+                alpha1,
+                alpharScheme
+            )
+          + fvc::flux
+            (
+                -fvc::flux(-phir, alpha3, alpharScheme),
+                alpha1,
+                alpharScheme
+            );
+
+        // Create the bounded (upwind) flux for alpha1
+        surfaceScalarField phiAlpha1BD =
+            upwind<scalar>(mesh, phi).flux(alpha1);
+
+        // Calculate the flux correction for alpha1
+        phiAlpha1 -= phiAlpha1BD;
+
+        // Calculate the limiter for alpha1
+        MULES::limiter
+        (
+            allLambda,
+            geometricOneField(),
+            alpha1,
+            phiAlpha1BD,
+            phiAlpha1,
+            zeroField(),
+            zeroField(),
+            1,
+            0,
+            3
+        );
+
+        // Create the complete flux for alpha2
+        surfaceScalarField phiAlpha2 =
+            fvc::flux
+            (
+                phi,
+                alpha2,
+                alphaScheme
+            )
+          + fvc::flux
+            (
+                -fvc::flux(phir, alpha1, alpharScheme),
+                alpha2,
+                alpharScheme
+            );
+
+        // Create the bounded (upwind) flux for alpha2
+        surfaceScalarField phiAlpha2BD =
+            upwind<scalar>(mesh, phi).flux(alpha2);
+
+        // Calculate the flux correction for alpha2
+        phiAlpha2 -= phiAlpha2BD;
+
+        // Further limit the limiter for alpha2
+        MULES::limiter
+        (
+            allLambda,
+            geometricOneField(),
+            alpha2,
+            phiAlpha2BD,
+            phiAlpha2,
+            zeroField(),
+            zeroField(),
+            1,
+            0,
+            3
+        );
+
+        // Construct the limited fluxes
+        phiAlpha1 = phiAlpha1BD + lambda*phiAlpha1;
+        phiAlpha2 = phiAlpha2BD + lambda*phiAlpha2;
+
+        // Solve for alpha1
+        solve(fvm::ddt(alpha1) + fvc::div(phiAlpha1));
+
+        // Create the diffusion coefficients for alpha2<->alpha3
+        volScalarField Dc23 = D23*max(alpha3, scalar(0))*pos(alpha2);
+        volScalarField Dc32 = D23*max(alpha2, scalar(0))*pos(alpha3);
+
+        // Add the diffusive flux for alpha3->alpha2
+        phiAlpha2 -= fvc::interpolate(Dc32)*mesh.magSf()*fvc::snGrad(alpha1);
+
+        // Solve for alpha2
+        fvScalarMatrix alpha2Eqn
+        (
+            fvm::ddt(alpha2)
+          + fvc::div(phiAlpha2)
+          - fvm::laplacian(Dc23 + Dc32, alpha2)
+        );
+        alpha2Eqn.solve();
+
+        // Construct the complete mass flux
+        rhoPhi =
+              phiAlpha1*(rho1 - rho3)
+            + (phiAlpha2 + alpha2Eqn.flux())*(rho2 - rho3)
+            + phi*rho3;
+
+        alpha3 = 1.0 - alpha1 - alpha2;
+    }
+
+    Info<< "Air phase volume fraction = "
+        << alpha1.weightedAverage(mesh.V()).value()
+        << "  Min(alpha1) = " << min(alpha1).value()
+        << "  Max(alpha1) = " << max(alpha1).value()
+        << endl;
+
+    Info<< "Liquid phase volume fraction = "
+        << alpha2.weightedAverage(mesh.V()).value()
+        << "  Min(alpha2) = " << min(alpha2).value()
+        << "  Max(alpha2) = " << max(alpha2).value()
+        << endl;
+}

applications/solvers/multiphase/interMixingFoam/alphaEqnsSubCycle.H

@@ -0,0 +1,43 @@
+label nAlphaCorr
+(
+    readLabel(piso.lookup("nAlphaCorr"))
+);
+
+label nAlphaSubCycles
+(
+    readLabel(piso.lookup("nAlphaSubCycles"))
+);
+
+if (nAlphaSubCycles > 1)
+{
+    surfaceScalarField rhoPhiSum = 0.0*rhoPhi;
+    dimensionedScalar totalDeltaT = runTime.deltaT();
+
+    for
+    (
+        subCycle<volScalarField> alphaSubCycle(alpha1, nAlphaSubCycles);
+        !(++alphaSubCycle).end();
+    )
+    {
+#       include "alphaEqns.H"
+        rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi;
+    }
+
+    rhoPhi = rhoPhiSum;
+}
+else
+{
+#       include "alphaEqns.H"
+}
+
+interface.correct();
+
+{
+    volScalarField rhoNew = alpha1*rho1 + alpha2*rho2 + alpha3*rho3;
+
+    //solve(fvm::ddt(rho) + fvc::div(rhoPhi));
+    //Info<< "density error = "
+    //    << max((mag(rho - rhoNew)/mag(rhoNew))().internalField()) << endl;
+
+    rho == rhoNew;
+}

applications/solvers/multiphase/interMixingFoam/createFields.H

@@ -1,4 +1,4 @@
-    Info<< "Reading field pd\n" << endl;
+    Info<< "Reading field p\n" << endl;
     volScalarField pd
     (
         IOobject
@@ -12,12 +12,12 @@
         mesh
     );
 
-    Info<< "Reading field gamma\n" << endl;
-    volScalarField gamma
+    Info<< "Reading field alpha1\n" << endl;
+    volScalarField alpha1
     (
         IOobject
         (
-            "gamma",
+            "alpha1",
             runTime.timeName(),
             mesh,
             IOobject::MUST_READ,
@@ -26,6 +26,39 @@
         mesh
     );
 
+
+    Info<< "Reading field alpha2\n" << endl;
+    volScalarField alpha2
+    (
+        IOobject
+        (
+            "alpha2",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+
+    Info<< "Reading field alpha3\n" << endl;
+    volScalarField alpha3
+    (
+        IOobject
+        (
+            "alpha3",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    alpha3 == 1.0 - alpha1 - alpha2;
+
+
     Info<< "Reading field U\n" << endl;
     volVectorField U
     (
@@ -42,12 +75,13 @@
 
 #   include "createPhi.H"
 
-    Info<< "Reading transportProperties\n" << endl;
-    twoPhaseMixture twoPhaseProperties(U, phi, "gamma");
+    threePhaseMixture threePhaseProperties(U, phi);
 
-    const dimensionedScalar& rho1 = twoPhaseProperties.rho1();
-    const dimensionedScalar& rho2 = twoPhaseProperties.rho2();
+    const dimensionedScalar& rho1 = threePhaseProperties.rho1();
+    const dimensionedScalar& rho2 = threePhaseProperties.rho2();
+    const dimensionedScalar& rho3 = threePhaseProperties.rho3();
 
+    dimensionedScalar D23(threePhaseProperties.lookup("D23"));
 
     // Need to store rho for ddt(rho, U)
     volScalarField rho
@@ -59,14 +93,15 @@
             mesh,
             IOobject::READ_IF_PRESENT
         ),
-        gamma*rho1 + (scalar(1) - gamma)*rho2,
-        gamma.boundaryField().types()
+        alpha1*rho1 + alpha2*rho2 + alpha3*rho3,
+        alpha1.boundaryField().types()
     );
     rho.oldTime();
 
+
     // Mass flux
     // Initialisation does not matter because rhoPhi is reset after the
-    // gamma solution before it is used in the U equation.
+    // alpha solution before it is used in the U equation.
     surfaceScalarField rhoPhi
     (
         IOobject
@@ -80,7 +115,6 @@
         rho1*phi
     );
 
-
     Info<< "Calculating field g.h\n" << endl;
     volScalarField gh("gh", g & mesh.C());
     surfaceScalarField ghf("gh", g & mesh.Cf());
@@ -99,7 +133,6 @@
         pd + rho*gh
     );
 
-
     label pdRefCell = 0;
     scalar pdRefValue = 0.0;
     setRefCell(pd, mesh.solutionDict().subDict("PISO"), pdRefCell, pdRefValue);
@@ -121,11 +154,12 @@
         );
     }
 
-    // Construct interface from gamma distribution
-    interfaceProperties interface(gamma, U, twoPhaseProperties);
+    // Construct interface from alpha distribution
+    threePhaseInterfaceProperties interface(threePhaseProperties);
 
-    // Construct LES model
-    autoPtr<incompressible::LESModel> turbulence
+
+    // Construct incompressible turbulence model
+    autoPtr<incompressible::turbulenceModel> turbulence
     (
-        incompressible::LESModel::New(U, phi, twoPhaseProperties)
+        incompressible::turbulenceModel::New(U, phi, threePhaseProperties)
     );

applications/solvers/multiphase/interMixingFoam/incompressibleThreePhaseMixture/threePhaseMixture.C

@@ -0,0 +1,204 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright held by original author
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+
+Class
+    threePhaseMixture
+
+\*---------------------------------------------------------------------------*/
+
+#include "threePhaseMixture.H"
+#include "addToRunTimeSelectionTable.H"
+#include "surfaceFields.H"
+#include "fvc.H"
+
+// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
+
+//- Calculate and return the laminar viscosity
+void Foam::threePhaseMixture::calcNu()
+{
+    // Average kinematic viscosity calculated from dynamic viscosity
+    nu_ = mu()/(alpha1_*rho1_ + alpha2_*rho2_ + alpha3_*rho3_);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::threePhaseMixture::threePhaseMixture
+(
+    const volVectorField& U,
+    const surfaceScalarField& phi
+)
+:
+    transportModel(U, phi),
+
+    phase1Name_("phase1"),
+    phase2Name_("phase2"),
+    phase3Name_("phase3"),
+
+    nuModel1_
+    (
+        viscosityModel::New
+        (
+            "nu1",
+            subDict(phase1Name_),
+            U,
+            phi
+        )
+    ),
+    nuModel2_
+    (
+        viscosityModel::New
+        (
+            "nu2",
+            subDict(phase2Name_),
+            U,
+            phi
+        )
+    ),
+    nuModel3_
+    (
+        viscosityModel::New
+        (
+            "nu3",
+            subDict(phase2Name_),
+            U,
+            phi
+        )
+    ),
+
+    rho1_(nuModel1_->viscosityProperties().lookup("rho")),
+    rho2_(nuModel2_->viscosityProperties().lookup("rho")),
+    rho3_(nuModel3_->viscosityProperties().lookup("rho")),
+
+    U_(U),
+    phi_(phi),
+
+    alpha1_(U_.db().lookupObject<const volScalarField> ("alpha1")),
+    alpha2_(U_.db().lookupObject<const volScalarField> ("alpha2")),
+    alpha3_(U_.db().lookupObject<const volScalarField> ("alpha3")),
+
+    nu_
+    (
+        IOobject
+        (
+            "nu",
+            U_.time().timeName(),
+            U_.db()
+        ),
+        U_.mesh(),
+        dimensionedScalar("nu", dimensionSet(0, 2, -1, 0, 0), 0),
+        calculatedFvPatchScalarField::typeName
+    )
+{
+    calcNu();
+}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField> Foam::threePhaseMixture::mu() const
+{
+    return tmp<volScalarField>
+    (
+        new volScalarField
+        (
+            "mu",
+            alpha1_*rho1_*nuModel1_->nu()
+          + alpha2_*rho2_*nuModel2_->nu()
+          + alpha3_*rho3_*nuModel3_->nu()
+        )
+    );
+}
+
+
+Foam::tmp<Foam::surfaceScalarField> Foam::threePhaseMixture::muf() const
+{
+    surfaceScalarField alpha1f = fvc::interpolate(alpha1_);
+    surfaceScalarField alpha2f = fvc::interpolate(alpha2_);
+    surfaceScalarField alpha3f = fvc::interpolate(alpha3_);
+
+    return tmp<surfaceScalarField>
+    (
+        new surfaceScalarField
+        (
+            "mu",
+            alpha1f*rho1_*fvc::interpolate(nuModel1_->nu())
+          + alpha2f*rho2_*fvc::interpolate(nuModel2_->nu())
+          + alpha3f*rho3_*fvc::interpolate(nuModel3_->nu())
+        )
+    );
+}
+
+
+Foam::tmp<Foam::surfaceScalarField> Foam::threePhaseMixture::nuf() const
+{
+    surfaceScalarField alpha1f = fvc::interpolate(alpha1_);
+    surfaceScalarField alpha2f = fvc::interpolate(alpha2_);
+    surfaceScalarField alpha3f = fvc::interpolate(alpha3_);
+
+    return tmp<surfaceScalarField>
+    (
+        new surfaceScalarField
+        (
+            "nu",
+            (
+                alpha1f*rho1_*fvc::interpolate(nuModel1_->nu())
+              + alpha2f*rho2_*fvc::interpolate(nuModel2_->nu())
+              + alpha3f*rho3_*fvc::interpolate(nuModel3_->nu())
+            )/(alpha1f*rho1_ + alpha2f*rho2_ + alpha3f*rho3_)
+        )
+    );
+}
+
+
+bool Foam::threePhaseMixture::read()
+{
+    if (transportModel::read())
+    {
+        if
+        (
+            nuModel1_().read(*this)
+         && nuModel2_().read(*this)
+         && nuModel3_().read(*this)
+        )
+        {
+            nuModel1_->viscosityProperties().lookup("rho") >> rho1_;
+            nuModel2_->viscosityProperties().lookup("rho") >> rho2_;
+            nuModel3_->viscosityProperties().lookup("rho") >> rho3_;
+
+            return true;
+        }
+        else
+        {
+            return false;
+        }
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+// ************************************************************************* //

applications/solvers/multiphase/interMixingFoam/incompressibleThreePhaseMixture/threePhaseMixture.H

@@ -0,0 +1,197 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright held by original author
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+
+Class
+    threePhaseMixture
+
+Description
+
+SourceFiles
+    threePhaseMixture.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef threePhaseMixture_H
+#define threePhaseMixture_H
+
+#include "incompressible/transportModel/transportModel.H"
+#include "incompressible/viscosityModels/viscosityModel/viscosityModel.H"
+#include "dimensionedScalar.H"
+#include "volFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class threePhaseMixture Declaration
+\*---------------------------------------------------------------------------*/
+
+class threePhaseMixture
+:
+    public transportModel
+{
+    // Private data
+
+        word phase1Name_;
+        word phase2Name_;
+        word phase3Name_;
+
+        autoPtr<viscosityModel> nuModel1_;
+        autoPtr<viscosityModel> nuModel2_;
+        autoPtr<viscosityModel> nuModel3_;
+
+        dimensionedScalar rho1_;
+        dimensionedScalar rho2_;
+        dimensionedScalar rho3_;
+
+        const volVectorField& U_;
+        const surfaceScalarField& phi_;
+
+        const volScalarField& alpha1_;
+        const volScalarField& alpha2_;
+        const volScalarField& alpha3_;
+
+        volScalarField nu_;
+
+
+    // Private Member Functions
+
+        //- Calculate and return the laminar viscosity
+        void calcNu();
+
+
+public:
+
+    // Constructors
+
+        //- Construct from components
+        threePhaseMixture
+        (
+            const volVectorField& U,
+            const surfaceScalarField& phi
+        );
+
+
+    // Destructor
+
+        ~threePhaseMixture()
+        {}
+
+
+    // Member Functions
+
+        //- Return const-access to phase1 viscosityModel
+        const viscosityModel& nuModel1() const
+        {
+            return nuModel1_();
+        }
+
+        //- Return const-access to phase2 viscosityModel
+        const viscosityModel& nuModel2() const
+        {
+            return nuModel2_();
+        }
+
+        //- Return const-access to phase3 viscosityModel
+        const viscosityModel& nuModel3() const
+        {
+            return nuModel3_();
+        }
+
+        //- Return const-access to phase1 density
+        const dimensionedScalar& rho1() const
+        {
+            return rho1_;
+        }
+
+        //- Return const-access to phase2 density
+        const dimensionedScalar& rho2() const
+        {
+            return rho2_;
+        };
+
+        //- Return const-access to phase3 density
+        const dimensionedScalar& rho3() const
+        {
+            return rho3_;
+        };
+
+        const volScalarField& alpha1() const
+        {
+            return alpha1_;
+        }
+
+        const volScalarField& alpha2() const
+        {
+            return alpha2_;
+        }
+
+        const volScalarField& alpha3() const
+        {
+            return alpha3_;
+        }
+
+        //- Return the velocity
+        const volVectorField& U() const
+        {
+            return U_;
+        }
+
+        //- Return the dynamic laminar viscosity
+        tmp<volScalarField> mu() const;
+
+        //- Return the face-interpolated dynamic laminar viscosity
+        tmp<surfaceScalarField> muf() const;
+
+        //- Return the kinematic laminar viscosity
+        const volScalarField& nu() const
+        {
+            return nu_;
+        }
+
+        //- Return the face-interpolated dynamic laminar viscosity
+        tmp<surfaceScalarField> nuf() const;
+
+        //- Correct the laminar viscosity
+        void correct()
+        {
+            calcNu();
+        }
+
+        //- Read base transportProperties dictionary
+        bool read();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/multiphase/interMixingFoam/interMixingFoam.C

@@ -20,26 +20,23 @@ License
 
     You should have received a copy of the GNU General Public License
     along with OpenFOAM; if not, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+    Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
 
 Application
-    lesInterFoam
+    interMixingFoam
 
 Description
-    Solver for 2 incompressible, isothermal immiscible fluids using a VOF
-    (volume of fluid) phase-fraction based interface capturing approach.
-    The momentum and other fluid properties are of the "mixture" and a single
-    momentum equation is solved.  Turbulence is modelled using a run-time
-    selectable incompressible LES model.
+    Solver for 3 incompressible fluids, two of which are miscible,
+    using a VOF method to capture the interface.
 
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
 #include "MULES.H"
 #include "subCycle.H"
-#include "interfaceProperties.H"
-#include "twoPhaseMixture.H"
-#include "incompressible/LESModel/LESModel.H"
+#include "threePhaseMixture.H"
+#include "threePhaseInterfaceProperties.H"
+#include "turbulenceModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -48,16 +45,16 @@ int main(int argc, char *argv[])
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
-    #include "readEnvironmentalProperties.H"
+    #include "readGravitationalAcceleration.H"
     #include "readPISOControls.H"
     #include "initContinuityErrs.H"
     #include "createFields.H"
     #include "readTimeControls.H"
-    #include "correctPhi.H"
     #include "CourantNo.H"
     #include "setInitialDeltaT.H"
+    #include "correctPhi.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 
@@ -69,12 +66,14 @@ int main(int argc, char *argv[])
         #include "setDeltaT.H"
 
         runTime++;
+
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        #include "gammaEqnSubCycle.H"
+        threePhaseProperties.correct();
 
-        turbulence->correct();
+        #include "alphaEqnsSubCycle.H"
 
+        #define twoPhaseProperties threePhaseProperties
         #include "UEqn.H"
 
         // --- PISO loop
@@ -85,17 +84,7 @@ int main(int argc, char *argv[])
 
         #include "continuityErrs.H"
 
-        p = pd + rho*gh;
-
-        if (pd.needReference())
-        {
-            p += dimensionedScalar
-            (
-                "p",
-                p.dimensions(),
-                pRefValue - getRefCellValue(p, pdRefCell)
-            );
-        }
+        turbulence->correct();
 
         runTime.write();
 
@@ -104,7 +93,7 @@ int main(int argc, char *argv[])
             << nl << endl;
     }
 
-    Info<< "End\n" << endl;
+    Info<< "\n end \n";
 
     return(0);
 }

applications/solvers/multiphase/interMixingFoam/threePhaseInterfaceProperties/threePhaseInterfaceProperties.C

@@ -0,0 +1,209 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright held by original author
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+
+Application
+    threePhaseInterfaceProperties
+
+Description
+    Properties to aid interFoam :
+    1. Correct the alpha boundary condition for dynamic contact angle.
+    2. Calculate interface curvature.
+
+\*---------------------------------------------------------------------------*/
+
+#include "threePhaseInterfaceProperties.H"
+#include "alphaContactAngleFvPatchScalarField.H"
+#include "mathematicalConstants.H"
+#include "surfaceInterpolate.H"
+#include "fvcDiv.H"
+#include "fvcGrad.H"
+
+// * * * * * * * * * * * * * * * Static Member Data  * * * * * * * * * * * * //
+
+const Foam::scalar Foam::threePhaseInterfaceProperties::convertToRad =
+    Foam::mathematicalConstant::pi/180.0;
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+// Correction for the boundary condition on the unit normal nHat on
+// walls to produce the correct contact angle.
+
+// The dynamic contact angle is calculated from the component of the
+// velocity on the direction of the interface, parallel to the wall.
+
+void Foam::threePhaseInterfaceProperties::correctContactAngle
+(
+    surfaceVectorField::GeometricBoundaryField& nHatb
+) const
+{
+    const volScalarField::GeometricBoundaryField& alpha1 =
+        mixture_.alpha1().boundaryField();
+    const volScalarField::GeometricBoundaryField& alpha2 =
+        mixture_.alpha2().boundaryField();
+    const volScalarField::GeometricBoundaryField& alpha3 =
+        mixture_.alpha3().boundaryField();
+    const volVectorField::GeometricBoundaryField& U =
+        mixture_.U().boundaryField();
+
+    const fvMesh& mesh = mixture_.U().mesh();
+    const fvBoundaryMesh& boundary = mesh.boundary();
+
+    forAll(boundary, patchi)
+    {
+        if (isA<alphaContactAngleFvPatchScalarField>(alpha1[patchi]))
+        {
+            const alphaContactAngleFvPatchScalarField& a2cap =
+                refCast<const alphaContactAngleFvPatchScalarField>
+                (alpha2[patchi]);
+
+            const alphaContactAngleFvPatchScalarField& a3cap =
+                refCast<const alphaContactAngleFvPatchScalarField>
+                (alpha3[patchi]);
+
+            scalarField twoPhaseAlpha2 = max(a2cap, scalar(0));
+            scalarField twoPhaseAlpha3 = max(a3cap, scalar(0));
+
+            scalarField sumTwoPhaseAlpha =
+                twoPhaseAlpha2 + twoPhaseAlpha3 + SMALL;
+
+            twoPhaseAlpha2 /= sumTwoPhaseAlpha;
+            twoPhaseAlpha3 /= sumTwoPhaseAlpha;
+
+            fvsPatchVectorField& nHatp = nHatb[patchi];
+
+            scalarField theta =
+                convertToRad
+               *(
+                   twoPhaseAlpha2*(180 - a2cap.theta(U[patchi], nHatp))
+                 + twoPhaseAlpha3*(180 - a3cap.theta(U[patchi], nHatp))
+               );
+
+            vectorField nf = boundary[patchi].nf();
+
+            // Reset nHatPatch to correspond to the contact angle
+
+            scalarField a12 = nHatp & nf;
+
+            scalarField b1 = cos(theta);
+
+            scalarField b2(nHatp.size());
+
+            forAll(b2, facei)
+            {
+                b2[facei] = cos(acos(a12[facei]) - theta[facei]);
+            }
+
+            scalarField det = 1.0 - a12*a12;
+
+            scalarField a = (b1 - a12*b2)/det;
+            scalarField b = (b2 - a12*b1)/det;
+
+            nHatp = a*nf + b*nHatp;
+
+            nHatp /= (mag(nHatp) + deltaN_.value());
+        }
+    }
+}
+
+
+void Foam::threePhaseInterfaceProperties::calculateK()
+{
+    const volScalarField& alpha1 = mixture_.alpha1();
+
+    const fvMesh& mesh = alpha1.mesh();
+    const surfaceVectorField& Sf = mesh.Sf();
+
+    // Cell gradient of alpha
+    volVectorField gradAlpha = fvc::grad(alpha1);
+
+    // Interpolated face-gradient of alpha
+    surfaceVectorField gradAlphaf = fvc::interpolate(gradAlpha);
+
+    // Face unit interface normal
+    surfaceVectorField nHatfv = gradAlphaf/(mag(gradAlphaf) + deltaN_);
+    correctContactAngle(nHatfv.boundaryField());
+
+    // Face unit interface normal flux
+    nHatf_ = nHatfv & Sf;
+
+    // Simple expression for curvature
+    K_ = -fvc::div(nHatf_);
+
+    // Complex expression for curvature.
+    // Correction is formally zero but numerically non-zero.
+    //volVectorField nHat = gradAlpha/(mag(gradAlpha) + deltaN_);
+    //nHat.boundaryField() = nHatfv.boundaryField();
+    //K_ = -fvc::div(nHatf_) + (nHat & fvc::grad(nHatfv) & nHat);
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::threePhaseInterfaceProperties::threePhaseInterfaceProperties
+(
+    const threePhaseMixture& mixture
+)
+:
+    mixture_(mixture),
+    cAlpha_
+    (
+        readScalar
+        (
+            mixture.U().mesh().solutionDict().subDict("PISO").lookup("cAlpha")
+        )
+    ),
+    sigma12_(mixture.lookup("sigma12")),
+    sigma13_(mixture.lookup("sigma13")),
+
+    deltaN_
+    (
+        "deltaN",
+        1e-8/pow(average(mixture.U().mesh().V()), 1.0/3.0)
+    ),
+
+    nHatf_
+    (
+        (
+            fvc::interpolate(fvc::grad(mixture.alpha1()))
+           /(mag(fvc::interpolate(fvc::grad(mixture.alpha1()))) + deltaN_)
+        ) & mixture.alpha1().mesh().Sf()
+    ),
+
+    K_
+    (
+        IOobject
+        (
+            "K",
+            mixture.alpha1().time().timeName(),
+            mixture.alpha1().mesh()
+        ),
+        -fvc::div(nHatf_)
+    )
+{
+    calculateK();
+}
+
+
+// ************************************************************************* //

applications/solvers/multiphase/interMixingFoam/threePhaseInterfaceProperties/threePhaseInterfaceProperties.H

@@ -0,0 +1,157 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright held by original author
+     \\/     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 2 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, write to the Free Software Foundation,
+    Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
+
+Class
+    threePhaseInterfaceProperties
+
+Description
+    Properties to aid interFoam :
+    1. Correct the alpha boundary condition for dynamic contact angle.
+    2. Calculate interface curvature.
+
+SourceFiles
+    threePhaseInterfaceProperties.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef threePhaseInterfaceProperties_H
+#define threePhaseInterfaceProperties_H
+
+#include "threePhaseMixture.H"
+#include "surfaceFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+              Class threePhaseInterfaceProperties Declaration
+\*---------------------------------------------------------------------------*/
+
+class threePhaseInterfaceProperties
+{
+    // Private data
+
+        const threePhaseMixture& mixture_;
+
+        //- Compression coefficient
+        scalar cAlpha_;
+
+        //- Surface tension 1-2
+        dimensionedScalar sigma12_;
+
+        //- Surface tension 1-3
+        dimensionedScalar sigma13_;
+
+        //- Stabilisation for normalisation of the interface normal
+        const dimensionedScalar deltaN_;
+
+        surfaceScalarField nHatf_;
+        volScalarField K_;
+
+
+    // Private Member Functions
+
+        //- Disallow default bitwise copy construct and assignment
+        threePhaseInterfaceProperties(const threePhaseInterfaceProperties&);
+        void operator=(const threePhaseInterfaceProperties&);
+
+        //- Correction for the boundary condition on the unit normal nHat on
+        //  walls to produce the correct contact dynamic angle
+        //  calculated from the component of U parallel to the wall
+        void correctContactAngle
+        (
+            surfaceVectorField::GeometricBoundaryField& nHat
+        ) const;
+
+        //- Re-calculate the interface curvature
+        void calculateK();
+
+
+public:
+
+    //- Conversion factor for degrees into radians
+    static const scalar convertToRad;
+
+
+    // Constructors
+
+        //- Construct from volume fraction field alpha and IOdictionary
+        threePhaseInterfaceProperties(const threePhaseMixture& mixture);
+
+
+    // Member Functions
+
+        scalar cAlpha() const
+        {
+            return cAlpha_;
+        }
+
+        const dimensionedScalar& deltaN() const
+        {
+            return deltaN_;
+        }
+
+        const surfaceScalarField& nHatf() const
+        {
+            return nHatf_;
+        }
+
+        const volScalarField& K() const
+        {
+            return K_;
+        }
+
+        tmp<volScalarField> sigma() const
+        {
+            volScalarField limitedAlpha2 = max(mixture_.alpha2(), scalar(0));
+            volScalarField limitedAlpha3 = max(mixture_.alpha3(), scalar(0));
+
+            return
+                (limitedAlpha2*sigma12_ + limitedAlpha3*sigma13_)
+               /(limitedAlpha2 + limitedAlpha3 + SMALL);
+        }
+
+        tmp<volScalarField> sigmaK() const
+        {
+            return sigma()*K_;
+        }
+
+        void correct()
+        {
+            calculateK();
+        }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/multiphase/interPhaseChangeFoam/Make/options

@@ -2,13 +2,14 @@ EXE_INC = \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/transportModels/incompressible/lnInclude \
     -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
-    -I$(LIB_SRC)/turbulenceModels/LES \
-    -I$(LIB_SRC)/turbulenceModels/LES/LESdeltas/lnInclude \
+    -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
     -IphaseChangeTwoPhaseMixtures/phaseChangeTwoPhaseMixture \
     -I$(LIB_SRC)/finiteVolume/lnInclude
 
 EXE_LIBS = \
     -linterfaceProperties \
     -lincompressibleTransportModels \
+    -lincompressibleRASModels \
     -lincompressibleLESModels \
-    -lfiniteVolume
+    -lfiniteVolume \
+    -llduSolvers

applications/solvers/multiphase/interPhaseChangeFoam/UEqn.H

@@ -1,15 +1,18 @@
-    surfaceScalarField muf =
+    surfaceScalarField muEff
+    (
+        "muEff",
         twoPhaseProperties->muf()
-      + fvc::interpolate(rho*turbulence->nuSgs());
+      + fvc::interpolate(rho*turbulence->nut())
+    );
 
     fvVectorMatrix UEqn
     (
         fvm::ddt(rho, U)
       + fvm::div(rhoPhi, U)
       - fvm::Sp(fvc::ddt(rho) + fvc::div(rhoPhi), U)
-      - fvm::laplacian(muf, U)
-      - (fvc::grad(U) & fvc::grad(muf))
-    //- fvc::div(muf*(fvc::interpolate(dev2(fvc::grad(U))) & mesh.Sf()))
+      - fvm::laplacian(muEff, U)
+      - (fvc::grad(U) & fvc::grad(muEff))
+    //- fvc::div(muEff*(fvc::interpolate(dev2(fvc::grad(U))) & mesh.Sf()))
     );
 
     UEqn.relax();
@@ -23,7 +26,7 @@
             fvc::reconstruct
             (
                 (
-                    fvc::interpolate(interface.sigmaK())*fvc::snGrad(gamma)
+                    fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
                   - ghf*fvc::snGrad(rho)
                   - fvc::snGrad(pd)
                 ) * mesh.magSf()

applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H

@@ -1,23 +1,23 @@
 {
-    word gammaScheme("div(phi,gamma)");
-    word gammarScheme("div(phirb,gamma)");
+    word alphaScheme("div(phi,alpha)");
+    word alpharScheme("div(phirb,alpha)");
 
     surfaceScalarField phir("phir", phic*interface.nHatf());
 
-    for (int gCorr=0; gCorr<nGammaCorr; gCorr++)
+    for (int gCorr=0; gCorr<nAlphaCorr; gCorr++)
     {
-        surfaceScalarField phiGamma =
+        surfaceScalarField phiAlpha =
             fvc::flux
             (
                 phi,
-                gamma,
-                gammaScheme
+                alpha1,
+                alphaScheme
             )
           + fvc::flux
             (
-                -fvc::flux(-phir, scalar(1) - gamma, gammarScheme),
-                gamma,
-                gammarScheme
+                -fvc::flux(-phir, scalar(1) - alpha1, alpharScheme),
+                alpha1,
+                alpharScheme
             );
 
         Pair<tmp<volScalarField> > vDotAlphal =
@@ -46,22 +46,22 @@
             ),
             // Divergence term is handled explicitly to be
             // consistent with the explicit transport solution
-            divU*gamma
+            divU*alpha1
           + vDotcAlphal
         );
 
-        //MULES::explicitSolve(gamma, phi, phiGamma, 1, 0);
-        //MULES::explicitSolve(oneField(), gamma, phi, phiGamma, Sp, Su, 1, 0);
-        MULES::implicitSolve(oneField(), gamma, phi, phiGamma, Sp, Su, 1, 0);
+        //MULES::explicitSolve(alpha1, phi, phiAlpha, 1, 0);
+        //MULES::explicitSolve(geometricOneField(), alpha1, phi, phiAlpha, Sp, Su, 1, 0);
+        MULES::implicitSolve(geometricOneField(), alpha1, phi, phiAlpha, Sp, Su, 1, 0);
 
         rhoPhi +=
             (runTime.deltaT()/totalDeltaT)
-           *(phiGamma*(rho1 - rho2) + phi*rho2);
+           *(phiAlpha*(rho1 - rho2) + phi*rho2);
     }
 
     Info<< "Liquid phase volume fraction = "
-        << gamma.weightedAverage(mesh.V()).value()
-        << "  Min(gamma) = " << min(gamma).value()
-        << "  Max(gamma) = " << max(gamma).value()
+        << alpha1.weightedAverage(mesh.V()).value()
+        << "  Min(alpha1) = " << min(alpha1).value()
+        << "  Max(alpha1) = " << max(alpha1).value()
         << endl;
 }

applications/solvers/multiphase/interPhaseChangeFoam/alphaEqnSubCycle.H

@@ -11,37 +11,37 @@ surfaceScalarField rhoPhi
 );
 
 {
-    label nGammaCorr
+    label nAlphaCorr
     (
-        readLabel(piso.lookup("nGammaCorr"))
+        readLabel(piso.lookup("nAlphaCorr"))
     );
 
-    label nGammaSubCycles
+    label nAlphaSubCycles
     (
-        readLabel(piso.lookup("nGammaSubCycles"))
+        readLabel(piso.lookup("nAlphaSubCycles"))
     );
 
     surfaceScalarField phic = mag(phi/mesh.magSf());
-    phic = min(interface.cGamma()*phic, max(phic));
+    phic = min(interface.cAlpha()*phic, max(phic));
 
     volScalarField divU = fvc::div(phi);
 
     dimensionedScalar totalDeltaT = runTime.deltaT();
 
-    if (nGammaSubCycles > 1)
+    if (nAlphaSubCycles > 1)
     {
         for
         (
-            subCycle<volScalarField> gammaSubCycle(gamma, nGammaSubCycles);
-            !(++gammaSubCycle).end();
+            subCycle<volScalarField> alphaSubCycle(alpha1, nAlphaSubCycles);
+            !(++alphaSubCycle).end();
         )
         {
-#           include "gammaEqn.H"
+#           include "alphaEqn.H"
         }
     }
     else
     {
-#       include "gammaEqn.H"
+#       include "alphaEqn.H"
     }
 
     if (nOuterCorr == 1)
@@ -49,5 +49,5 @@ surfaceScalarField rhoPhi
         interface.correct();
     }
 
-    rho == gamma*rho1 + (scalar(1) - gamma)*rho2;
+    rho == alpha1*rho1 + (scalar(1) - alpha1)*rho2;
 }

applications/solvers/multiphase/interPhaseChangeFoam/correctPhi.H

@@ -1,7 +1,11 @@
 {
 #   include "continuityErrs.H"
 
-    wordList pcorrTypes(pd.boundaryField().types());
+    wordList pcorrTypes
+    (
+        pd.boundaryField().size(),
+        zeroGradientFvPatchScalarField::typeName
+    );
 
     for (label i=0; i<pd.boundaryField().size(); i++)
     {

applications/solvers/multiphase/interPhaseChangeFoam/createFields.H

@@ -12,12 +12,12 @@
         mesh
     );
 
-    Info<< "Reading field gamma\n" << endl;
-    volScalarField gamma
+    Info<< "Reading field alpha1\n" << endl;
+    volScalarField alpha1
     (
         IOobject
         (
-            "gamma",
+            "alpha1",
             runTime.timeName(),
             mesh,
             IOobject::MUST_READ,
@@ -44,7 +44,7 @@
 
     Info<< "Creating phaseChangeTwoPhaseMixture\n" << endl;
     autoPtr<phaseChangeTwoPhaseMixture> twoPhaseProperties =
-        phaseChangeTwoPhaseMixture::New(U, phi, "gamma");
+        phaseChangeTwoPhaseMixture::New(U, phi, "alpha1");
 
     const dimensionedScalar& rho1 = twoPhaseProperties->rho1();
     const dimensionedScalar& rho2 = twoPhaseProperties->rho2();
@@ -60,19 +60,32 @@
             mesh,
             IOobject::READ_IF_PRESENT
         ),
-        gamma*rho1 + (scalar(1) - gamma)*rho2,
-        gamma.boundaryField().types()
+        alpha1*rho1 + (scalar(1) - alpha1)*rho2,
+        alpha1.boundaryField().types()
     );
     rho.oldTime();
 
 
-    label pdRefCell = 0;
-    scalar pdRefValue = 0.0;
-    setRefCell(pd, mesh.solutionDict().subDict("PISO"), pdRefCell, pdRefValue);
+    // Mass flux
+    // Initialisation does not matter because rhoPhi is reset after the
+    // alpha1 solution before it is used in the U equation.
+    surfaceScalarField rhoPhi
+    (
+        IOobject
+        (
+            "rho*phi",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        rho1*phi
+    );
 
-    Info<< "Calculating field g.h" << endl;
+
+    Info<< "Calculating field g.h\n" << endl;
     volScalarField gh("gh", g & mesh.C());
-    surfaceScalarField ghf("ghf", g & mesh.Cf());
+    surfaceScalarField ghf("gh", g & mesh.Cf());
 
     volScalarField p
     (
@@ -88,11 +101,32 @@
     );
 
 
-    // Construct interface from gamma distribution
-    interfaceProperties interface(gamma, U, twoPhaseProperties());
+    label pdRefCell = 0;
+    scalar pdRefValue = 0.0;
+    setRefCell(pd, mesh.solutionDict().subDict("PISO"), pdRefCell, pdRefValue);
 
-    // Construct LES model
-    autoPtr<incompressible::LESModel> turbulence
+    scalar pRefValue = 0.0;
+
+    if (pd.needReference())
+    {
+        pRefValue = readScalar
         (
-        incompressible::LESModel::New(U, phi, twoPhaseProperties())
+            mesh.solutionDict().subDict("PISO").lookup("pRefValue")
+        );
+
+        p += dimensionedScalar
+        (
+            "p",
+            p.dimensions(),
+            pRefValue - getRefCellValue(p, pdRefCell)
+        );
+    }
+
+    // Construct interface from alpha1 distribution
+    interfaceProperties interface(alpha1, U, twoPhaseProperties());
+
+    // Construct incompressible turbulence model
+    autoPtr<incompressible::turbulenceModel> turbulence
+    (
+        incompressible::turbulenceModel::New(U, phi, twoPhaseProperties())
     );

applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C

@@ -28,13 +28,17 @@ Application
 Description
     Solver for 2 incompressible, isothermal immiscible fluids with phase-change
     (e.g. cavitation).  Uses a VOF (volume of fluid) phase-fraction based
-    interface capturing approach.  The momentum and other fluid properties are
-    of the "mixture" and a single momentum equation is solved.
+    interface capturing approach.
+
+    The momentum and other fluid properties are of the "mixture" and a
+    single momentum equation is solved.
 
     The set of phase-change models provided are designed to simulate cavitation
     but other mechanisms of phase-change are supported within this solver
     framework.
 
+    Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
+
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
@@ -42,7 +46,7 @@ Description
 #include "subCycle.H"
 #include "interfaceProperties.H"
 #include "phaseChangeTwoPhaseMixture.H"
-#include "incompressible/LESModel/LESModel.H"
+#include "turbulenceModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -51,7 +55,7 @@ int main(int argc, char *argv[])
 #   include "setRootCase.H"
 #   include "createTime.H"
 #   include "createMesh.H"
-    #include "readEnvironmentalProperties.H"
+#   include "readGravitationalAcceleration.H"
 #   include "readPISOControls.H"
 #   include "initContinuityErrs.H"
 #   include "createFields.H"
@@ -75,9 +79,7 @@ int main(int argc, char *argv[])
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        twoPhaseProperties->correct();
-
-        #include "gammaEqnSubCycle.H"
+#       include "alphaEqnSubCycle.H"
 
         turbulence->correct();
 
@@ -95,6 +97,8 @@ int main(int argc, char *argv[])
 #            include "continuityErrs.H"
         }
 
+        twoPhaseProperties->correct();
+
         runTime.write();
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
@@ -104,7 +108,7 @@ int main(int argc, char *argv[])
 
     Info<< "End\n" << endl;
 
-    return(0);
+    return 0;
 }
 
 

applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H

@@ -11,13 +11,14 @@
       + fvc::ddtPhiCorr(rUA, rho, U, phi)
     );
 
+    adjustPhi(phiU, U, p);
+
     phi = phiU +
         (
-            fvc::interpolate(interface.sigmaK())*fvc::snGrad(gamma)
+            fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
           - ghf*fvc::snGrad(rho)
         )*rUAf*mesh.magSf();
 
-    adjustPhi(phi, U, pd);
 
     Pair<tmp<volScalarField> > vDotP = twoPhaseProperties->vDotP();
     const volScalarField& vDotcP = vDotP[0]();

applications/solvers/multiphase/lesCavitatingFoam/Make/files

@@ -1,5 +0,0 @@
-lesCavitatingFoam.C
-
-devOneEqEddy/devOneEqEddy.C
-
-EXE = $(FOAM_APPBIN)/lesCavitatingFoam

applications/solvers/multiphase/lesCavitatingFoam/UEqn.H

@@ -1,20 +0,0 @@
-    surfaceScalarField muEff
-    (
-        "muEff",
-        twoPhaseProperties.muf()
-      + fvc::interpolate(rho*turbulence->nuSgs())
-    );
-
-    fvVectorMatrix UEqn
-    (
-        fvm::ddt(rho, U)
-      + fvm::div(phi, U)
-      - fvm::laplacian(muEff, U)
-    //- (fvc::grad(U) & fvc::grad(muf))
-      - fvc::div(muEff*(fvc::interpolate(dev(fvc::grad(U))) & mesh.Sf()))
-    );
-
-    if (momentumPredictor)
-    {
-        solve(UEqn == -fvc::grad(p));
-    }

applications/solvers/multiphase/lesCavitatingFoam/createFields.H

@@ -1,85 +0,0 @@
-    Info<< "Reading field p\n" << endl;
-    volScalarField p
-    (
-        IOobject
-        (
-            "p",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-    volScalarField rho
-    (
-        IOobject
-        (
-            "rho",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-    volScalarField gamma
-    (
-        IOobject
-        (
-            "gamma",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        max(min((rho - rholSat)/(rhovSat - rholSat), scalar(1)), scalar(0))
-    );
-    gamma.oldTime();
-
-    Info<< "Creating compressibilityModel\n" << endl;
-    autoPtr<barotropicCompressibilityModel> psiModel = 
-        barotropicCompressibilityModel::New
-        (
-            thermodynamicProperties,
-            gamma
-        );
-
-    const volScalarField& psi = psiModel->psi();
-
-    rho == max
-    (
-        psi*p
-      + (1.0 - gamma)*rhol0
-      + ((gamma*psiv + (1.0 - gamma)*psil) - psi)*pSat,
-        rhoMin
-    );
-
-    Info<< "Reading field U\n" << endl;
-    volVectorField U
-    (
-        IOobject
-        (
-            "U",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-#   include "createPhiv.H"
-#   include "compressibleCreatePhi.H"
-
-    Info<< "Reading transportProperties\n" << endl;
-
-    twoPhaseMixture twoPhaseProperties(U, phiv, "gamma");
-
-    // Create LES model
-    autoPtr<incompressible::LESModel> turbulence
-    (
-        incompressible::LESModel::New(U, phiv, twoPhaseProperties)
-    );

applications/solvers/multiphase/lesCavitatingFoam/devOneEqEddy/devOneEqEddy.C

@@ -1,130 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright held by original author
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-\*---------------------------------------------------------------------------*/
-
-#include "devOneEqEddy.H"
-#include "addToRunTimeSelectionTable.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace incompressible
-{
-namespace LESModels
-{
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-defineTypeNameAndDebug(devOneEqEddy, 0);
-addToRunTimeSelectionTable(LESModel, devOneEqEddy, dictionary);
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-devOneEqEddy::devOneEqEddy
-(
-    const volVectorField& U,
-    const surfaceScalarField& phi,
-    transportModel& transport
-)
-:
-    LESModel(typeName, U, phi, transport),
-    GenEddyVisc(U, phi, transport),
-
-    k_
-    (
-        IOobject
-        (
-            "k",
-            runTime_.timeName(),
-            mesh_,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh_
-    ),
-    ck_
-    (
-        dimensioned<scalar>::lookupOrAddToDict
-        (
-            "ck",
-            coeffDict(),
-            0.07
-        )
-    )
-{
-    printCoeffs();
-}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-void devOneEqEddy::correct(const tmp<volTensorField>& gradU)
-{
-    GenEddyVisc::correct(gradU);
-
-    //volScalarField G = 2*nuSgs_*magSqr(symm(gradU));
-    volScalarField G = 2*nuSgs_*(gradU() && dev(symm(gradU())));
-
-    solve
-    (
-       fvm::ddt(k_)
-     + fvm::div(phi(), k_)
-     - fvm::Sp(fvc::div(phi()), k_)
-     - fvm::laplacian(DkEff(), k_)
-    ==
-       G
-     - fvm::Sp(ce_*sqrt(k_)/delta(), k_)
-    );
-
-    bound(k_, k0());
-
-    nuSgs_ = ck_*sqrt(k_)*delta();
-    nuSgs_.correctBoundaryConditions();
-}
-
-
-bool devOneEqEddy::read()
-{
-    if (GenEddyVisc::read())
-    {
-        ck_.readIfPresent(coeffDict());
-
-        return true;
-    }
-    else
-    {
-        return false;
-    }
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace LESModels
-} // End namespace incompressible
-} // End namespace Foam
-
-// ************************************************************************* //

applications/solvers/multiphase/lesCavitatingFoam/lesCavitatingFoam.C

@@ -1,94 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright held by original author
-     \\/     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 2 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-Application
-    lesCavitatingFoam
-
-Description
-    Transient cavitation code with LES turbulence.
-
-\*---------------------------------------------------------------------------*/
-
-#include "fvCFD.H"
-#include "barotropicCompressibilityModel.H"
-#include "twoPhaseMixture.H"
-#include "incompressible/LESModel/LESModel.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-int main(int argc, char *argv[])
-{
-
-#   include "setRootCase.H"
-
-#   include "createTime.H"
-#   include "createMesh.H"
-#   include "readThermodynamicProperties.H"
-#   include "readControls.H"
-#   include "createFields.H"
-#   include "initContinuityErrs.H"
-#   include "compressibleCourantNo.H"
-#   include "setInitialDeltaT.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-    Info<< "\nStarting time loop\n" << endl;
-
-    while (runTime.run())
-    {
-#       include "readControls.H"
-#       include "CourantNo.H"
-#       include "setDeltaT.H"
-
-        runTime++;
-        Info<< "Time = " << runTime.timeName() << nl << endl;
-
-        turbulence->correct();
-
-        for (int outerCorr=0; outerCorr<nOuterCorr; outerCorr++)
-        {
-#           include "rhoEqn.H"
-#           include "gammaPsi.H"
-#           include "UEqn.H"
-
-            for (int corr=0; corr<nCorr; corr++)
-            {
-#               include "pEqn.H"
-            }
-        }
-
-        runTime.write();
-
-        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
-            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
-            << nl << endl;
-    }
-
-    Info<< "\n end \n";
-
-    return(0);
-}
-
-
-// ************************************************************************* //

applications/solvers/multiphase/lesCavitatingFoam/pEqn.H

@@ -1,80 +0,0 @@
-{
-    if (nOuterCorr == 1)
-    {
-        p = 
-        (
-            rho
-          - (1.0 - gamma)*rhol0
-          - ((gamma*psiv + (1.0 - gamma)*psil) - psi)*pSat
-        )/psi;
-    }
-
-    surfaceScalarField rhof = fvc::interpolate(rho, "rhof");
-
-    volScalarField rUA = 1.0/UEqn.A();
-    surfaceScalarField rUAf("rUAf", rhof*fvc::interpolate(rUA));
-    volVectorField HbyA = rUA*UEqn.H();
-
-    phiv = (fvc::interpolate(HbyA) & mesh.Sf())
-         + fvc::ddtPhiCorr(rUA, rho, U, phiv);
-
-    p.boundaryField().updateCoeffs();
-
-    surfaceScalarField phiGradp = rUAf*mesh.magSf()*fvc::snGrad(p);
-
-    phiv -= phiGradp/rhof;
-
-#   include "resetPhivPatches.H"
-
-    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
-    {
-        fvScalarMatrix pEqn
-        (
-            fvm::ddt(psi, p)
-          - (rhol0 + (psil - psiv)*pSat)*fvc::ddt(gamma) - pSat*fvc::ddt(psi)
-          + fvc::div(phiv, rho)
-          + fvc::div(phiGradp)
-          - fvm::laplacian(rUAf, p)
-        );
-
-        pEqn.solve();
-
-        if (nonOrth == nNonOrthCorr)
-        {
-            phiv += (phiGradp + pEqn.flux())/rhof;
-        }
-    }
-
-    Info<< "max-min p: " << max(p).value()
-        << " " << min(p).value() << endl;
-
-
-    U = HbyA - rUA*fvc::grad(p);
-
-    // Remove the swirl component of velocity for "wedge" cases
-    if (piso.found("removeSwirl"))
-    {
-        label swirlCmpt(readLabel(piso.lookup("removeSwirl")));
-
-        Info<< "Removing swirl component-" << swirlCmpt << " of U" << endl;
-        U.field().replace(swirlCmpt, 0.0);
-    }
-
-    U.correctBoundaryConditions();
-
-    Info<< "max(U) " << max(mag(U)).value() << endl;
-
-    rho == max
-    (
-        psi*p
-      + (1.0 - gamma)*rhol0
-      + ((gamma*psiv + (1.0 - gamma)*psil) - psi)*pSat,
-        rhoMin
-    );
-
-    Info<< "max-min rho: " << max(rho).value()
-        << " " << min(rho).value() << endl;
-
-#   include "gammaPsi.H"
-
-}

applications/solvers/multiphase/lesInterFoam/Make/files

@@ -1,3 +0,0 @@
-lesInterFoam.C
-
-EXE = $(FOAM_APPBIN)/lesInterFoam

applications/solvers/multiphase/lesInterFoam/Make/options

@@ -1,16 +0,0 @@
-EXE_INC = \
-    -Iaveraging \
-    -I../interFoam \
-    -I$(LIB_SRC)/transportModels \
-    -I$(LIB_SRC)/transportModels/incompressible/lnInclude \
-    -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
-    -I$(LIB_SRC)/turbulenceModels/LES \
-    -I$(LIB_SRC)/turbulenceModels/LES/LESdeltas/lnInclude \
-    -I$(LIB_SRC)/finiteVolume/lnInclude
-
-EXE_LIBS = \
-    -linterfaceProperties \
-    -lincompressibleTransportModels \
-    -lincompressibleLESModels \
-    -lfiniteVolume \
-    -llduSolvers

applications/solvers/multiphase/multiphaseInterFoam/Make/options

@@ -6,10 +6,14 @@ EXE_INC = \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/transportModels/incompressible/lnInclude \
     -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
+    -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
     -I$(LIB_SRC)/finiteVolume/lnInclude
 
 EXE_LIBS = \
     -linterfaceProperties \
     -lincompressibleTransportModels \
+    -lincompressibleTransportModels \
+    -lincompressibleRASModels \
+    -lincompressibleLESModels \
     -lfiniteVolume \
     -llduSolvers

applications/solvers/multiphase/multiphaseInterFoam/createFields.H

@@ -46,9 +46,46 @@
 
 
     Info<< "Calculating field g.h\n" << endl;
+    volScalarField gh("gh", g & mesh.C());
     surfaceScalarField ghf("gh", g & mesh.Cf());
 
+    volScalarField p
+    (
+        IOobject
+        (
+            "p",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        pd + rho*gh
+    );
+
 
     label pdRefCell = 0;
     scalar pdRefValue = 0.0;
     setRefCell(pd, mesh.solutionDict().subDict("PISO"), pdRefCell, pdRefValue);
+
+    scalar pRefValue = 0.0;
+
+    if (pd.needReference())
+    {
+        pRefValue = readScalar
+        (
+            mesh.solutionDict().subDict("PISO").lookup("pRefValue")
+        );
+
+        p += dimensionedScalar
+        (
+            "p",
+            p.dimensions(),
+            pRefValue - getRefCellValue(p, pdRefCell)
+        );
+    }
+
+    // Construct incompressible turbulence model
+    autoPtr<incompressible::turbulenceModel> turbulence
+    (
+        incompressible::turbulenceModel::New(U, phi, mixture)
+    );

applications/solvers/multiphase/multiphaseInterFoam/multiphaseInterFoam.C

@@ -27,22 +27,24 @@ Application
 
 Description
     Solver for n incompressible fluids which captures the interfaces and
-    includes surface-tension and contact-angle effects for each.
+    includes surface-tension and contact-angle effects for each phase.
+
+    Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
 
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
 #include "multiphaseMixture.H"
+#include "turbulenceModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 int main(int argc, char *argv[])
 {
-
 #   include "setRootCase.H"
 #   include "createTime.H"
 #   include "createMesh.H"
-#   include "readEnvironmentalProperties.H"
+    #include "readGravitationalAcceleration.H"
 #   include "readPISOControls.H"
 #   include "initContinuityErrs.H"
 #   include "createFields.H"
@@ -79,6 +81,8 @@ int main(int argc, char *argv[])
 
 #       include "continuityErrs.H"
 
+        turbulence->correct();
+
         runTime.write();
 
         Info<< "ExecutionTime = "
@@ -88,7 +92,7 @@ int main(int argc, char *argv[])
 
     Info<< "\n end \n";
 
-    return(0);
+    return 0;
 }
 
 

applications/solvers/multiphase/porousInterFoam/Make/files

@@ -0,0 +1,3 @@
+porousInterFoam.C
+
+EXE = $(FOAM_APPBIN)/porousInterFoam

applications/solvers/multiphase/porousInterFoam/Make/options

@@ -1,15 +1,19 @@
 INTERFOAM = $(FOAM_SOLVERS)/multiphase/interFoam
 
 EXE_INC = \
+    -I$(INTERFOAM) \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/transportModels/incompressible/lnInclude \
     -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
-    -I$(LIB_SRC)/turbulenceModels/LES \
-    -I$(LIB_SRC)/turbulenceModels/LES/LESdeltas/lnInclude \
-    -I$(LIB_SRC)/finiteVolume/lnInclude
+    -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude
 
 EXE_LIBS = \
     -linterfaceProperties \
     -lincompressibleTransportModels \
+    -lincompressibleRASModels \
     -lincompressibleLESModels \
-    -lfiniteVolume
+    -lfiniteVolume \
+    -lmeshTools
+