Merge commit 'OpenCFD/master' into olesenm

applications/solvers/combustion/PDRFoam/PDRModels/dragModels/basic/basic.C

@@ -117,7 +117,7 @@ Foam::tmp<Foam::volSymmTensorField> Foam::PDRDragModels::basic::Dcu() const
 {
     const volScalarField& betav = U_.db().lookupObject<volScalarField>("betav");
 
-    return rho_*CR_*mag(U_) + (Csu*I)*betav*turbulence_.muEff()*Aw2_;
+    return (0.5*rho_)*CR_*mag(U_) + (Csu*I)*betav*turbulence_.muEff()*Aw2_;
 }
 
 
@@ -125,8 +125,8 @@ Foam::tmp<Foam::volScalarField> Foam::PDRDragModels::basic::Gk() const
 {
     const volScalarField& betav = U_.db().lookupObject<volScalarField>("betav");
 
-    return 
+    return
-        rho_*mag(U_)*(U_ & CT_ & U_)
+        (0.5*rho_)*mag(U_)*(U_ & CT_ & U_)
       + Csk*betav*turbulence_.muEff()*Aw2_*magSqr(U_);
 }
 

applications/solvers/incompressible/boundaryFoam/Make/options

@@ -1,7 +1,9 @@
 EXE_INC = \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
-    -I$(LIB_SRC)/turbulenceModels/RAS \
+    -I$(LIB_SRC)/turbulenceModels/incompressible/RAS/RASModel \
     -I$(LIB_SRC)/transportModels \
+    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/sampling/lnInclude
 
 EXE_LIBS = \

applications/solvers/incompressible/boundaryFoam/boundaryFoam.C

@@ -37,8 +37,8 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "incompressible/singlePhaseTransportModel/singlePhaseTransportModel.H"
+#include "singlePhaseTransportModel.H"
-#include "incompressible/RASModel/RASModel.H"
+#include "RASModel.H"
 #include "wallFvPatch.H"
 #include "makeGraph.H"
 

applications/solvers/incompressible/channelFoam/Make/files

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

applications/solvers/incompressible/channelFoam/Make/options

@@ -1,9 +1,10 @@
 EXE_INC = \
-    -I$(LIB_SRC)/turbulenceModels/LES \
+    -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
+    -I$(LIB_SRC)/turbulenceModels/incompressible/LES/LESModel \
     -I$(LIB_SRC)/turbulenceModels/LES/LESdeltas/lnInclude \
     -I$(LIB_SRC)/transportModels \
+    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/sampling/lnInclude
 
 EXE_LIBS = \

applications/solvers/incompressible/channelFoam/channelFoam.C

@@ -23,7 +23,7 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Application
-    oodles
+    channelFoam
 
 Description
     Incompressible LES solver for flow in a channel.
@@ -31,8 +31,8 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "incompressible/singlePhaseTransportModel/singlePhaseTransportModel.H"
+#include "singlePhaseTransportModel.H"
-#include "incompressible/LESModel/LESModel.H"
+#include "LESModel.H"
 #include "IFstream.H"
 #include "OFstream.H"
 #include "Random.H"

applications/solvers/incompressible/channelOodles/Make/files

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

applications/solvers/incompressible/channelOodles/Make/options

@@ -1,13 +0,0 @@
-EXE_INC = \
-    -I$(LIB_SRC)/turbulenceModels/LES \
-    -I$(LIB_SRC)/turbulenceModels/LES/LESdeltas/lnInclude \
-    -I$(LIB_SRC)/transportModels \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/sampling/lnInclude \
-    -I../oodles
-
-EXE_LIBS = \
-    -lincompressibleLESModels \
-    -lincompressibleTransportModels \
-    -lfiniteVolume \
-    -lmeshTools

applications/solvers/incompressible/icoDyMFoam/Make/files

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

applications/solvers/incompressible/icoDyMFoam/Make/options

@@ -1,11 +0,0 @@
-EXE_INC = \
-    -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
-    -I$(LIB_SRC)/dynamicMesh/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/finiteVolume/lnInclude
-
-EXE_LIBS = \
-    -ldynamicFvMesh \
-    -ldynamicMesh \
-    -lmeshTools \
-    -lfiniteVolume

applications/solvers/incompressible/icoDyMFoam/UEqn.H

@@ -1,16 +0,0 @@
-    fvVectorMatrix UEqn
-    (
-        fvm::ddt(U)
-      + fvm::div(phi, U)
-      - fvm::laplacian(nu, U)
-    );
-
-    if (ocorr != nOuterCorr-1)
-    {
-        UEqn.relax();
-    }
-
-    if (momentumPredictor)
-    {
-        solve(UEqn == -fvc::grad(p));
-    }

applications/solvers/incompressible/icoDyMFoam/checkTotalVolume.H

@@ -1,6 +0,0 @@
-    scalar newTotalVolume = sum(mesh.V());
-    scalar totalVolRatio = newTotalVolume/totalVolume;
-
-    Info << "Total volume change: " << totalVolRatio - 1 << endl;
-
-    totalVolume = newTotalVolume;

applications/solvers/incompressible/icoDyMFoam/createFields.H

@@ -1,72 +0,0 @@
-    Info<< "Reading transportProperties\n" << endl;
-
-    IOdictionary transportProperties
-    (
-        IOobject
-        (
-            "transportProperties",
-            runTime.constant(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE
-        )
-    );
-
-    dimensionedScalar nu
-    (
-        transportProperties.lookup("nu")
-    );
-
-
-    Info<< "Reading field p\n" << endl;
-    volScalarField p
-    (
-        IOobject
-        (
-            "p",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-
-    Info<< "Reading field U\n" << endl;
-    volVectorField U
-    (
-        IOobject
-        (
-            "U",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-#   include "createPhi.H"
-
-
-    label pRefCell = 0;
-    scalar pRefValue = 0.0;
-    setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);
-
-
-    Info<< "Reading field rAU if present\n" << endl;
-    volScalarField rAU
-    (
-        IOobject
-        (
-            "rAU",
-            runTime.timeName(),
-            mesh,
-            IOobject::READ_IF_PRESENT,
-            IOobject::AUTO_WRITE
-        ),
-        mesh,
-        runTime.deltaT(),
-        zeroGradientFvPatchScalarField::typeName
-    );

applications/solvers/incompressible/nonNewtonianIcoFoam/Make/options

@@ -1,6 +1,7 @@
 EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/transportModels
+    -I$(LIB_SRC)/transportModels \
+    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel
 
 EXE_LIBS = \
     -lfiniteVolume \

applications/solvers/incompressible/nonNewtonianIcoFoam/nonNewtonianIcoFoam.C

@@ -31,7 +31,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "incompressible/singlePhaseTransportModel/singlePhaseTransportModel.H"
+#include "singlePhaseTransportModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

applications/solvers/incompressible/oodles/Make/files

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

applications/solvers/incompressible/oodles/oodles.C

@@ -1,128 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
-     \\/     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
-    oodles
-
-Description
-    Incompressible LES solver.
-
-\*---------------------------------------------------------------------------*/
-
-#include "fvCFD.H"
-#include "incompressible/singlePhaseTransportModel/singlePhaseTransportModel.H"
-#include "incompressible/transportModel/transportModel.H"
-#include "incompressible/LESModel/LESModel.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-int main(int argc, char *argv[])
-{
-    #include "setRootCase.H"
-    #include "createTime.H"
-    #include "createMeshNoClear.H"
-    #include "createFields.H"
-    #include "initContinuityErrs.H"
-
-    Info<< "\nStarting time loop\n" << endl;
-
-    for (runTime++; !runTime.end(); runTime++)
-    {
-        Info<< "Time = " << runTime.timeName() << nl << endl;
-
-        #include "readPISOControls.H"
-        #include "CourantNo.H"
-
-        sgsModel->correct();
-
-        fvVectorMatrix UEqn
-        (
-            fvm::ddt(U)
-          + fvm::div(phi, U)
-          + sgsModel->divDevBeff(U)
-        );
-
-        // Optionally ensure diagonal-dominance of the momentum matrix
-        UEqn.relax();
-
-        if (momentumPredictor)
-        {
-            solve(UEqn == -fvc::grad(p));
-        }
-
-        // --- PISO loop
-        for (int corr=0; corr<nCorr; corr++)
-        {
-            volScalarField rUA = 1.0/UEqn.A();
-
-            U = rUA*UEqn.H();
-            phi = (fvc::interpolate(U) & mesh.Sf())
-                + fvc::ddtPhiCorr(rUA, U, phi);
-
-            adjustPhi(phi, U, p);
-
-            for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
-            {
-                fvScalarMatrix pEqn
-                (
-                    fvm::laplacian(rUA, p) == fvc::div(phi)
-                );
-
-                pEqn.setReference(pRefCell, pRefValue);
-
-                if (corr == nCorr-1 && nonOrth == nNonOrthCorr)
-                {
-                    pEqn.solve(mesh.solver(p.name() + "Final"));
-                }
-                else
-                {
-                    pEqn.solve(mesh.solver(p.name()));
-                }
-
-                if (nonOrth == nNonOrthCorr)
-                {
-                    phi -= pEqn.flux();
-                }
-            }
-
-            #include "continuityErrs.H"
-
-            U -= rUA*fvc::grad(p);
-            U.correctBoundaryConditions();
-        }
-
-        runTime.write();
-
-        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
-            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
-            << nl << endl;
-    }
-
-    Info<< "End\n" << endl;
-
-    return(0);
-}
-
-
-// ************************************************************************* //

applications/solvers/incompressible/pimpleDyMFoam/Make/files

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

applications/solvers/incompressible/pimpleDyMFoam/Make/options

@@ -2,14 +2,16 @@ EXE_INC = \
     -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/turbulenceModels/RAS \
+    -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
     -I$(LIB_SRC)/transportModels \
+    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
     -I$(LIB_SRC)/finiteVolume/lnInclude
 
 EXE_LIBS = \
     -ldynamicFvMesh \
     -ldynamicMesh \
     -lmeshTools \
-    -lincompressibleRASModels \
     -lincompressibleTransportModels \
+    -lincompressibleRASModels \
+    -lincompressibleLESModels \
     -lfiniteVolume

applications/solvers/incompressible/pimpleDyMFoam/createFields.H

@@ -37,9 +37,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 rAU if present\n" << endl;

applications/solvers/incompressible/pimpleDyMFoam/pimpleDyMFoam.C

@@ -23,15 +23,19 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Application
-    icoDyMFoam
+    turbDyMFoam
 
 Description
-    Transient solver for incompressible, laminar flow of Newtonian fluids
+    Transient solver for incompressible, flow of Newtonian fluids
-    with moving mesh.
+    on a moving mesh using the PIMPLE (merged PISO-SIMPLE) algorithm.
+
+    Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
 
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
+#include "singlePhaseTransportModel.H"
+#include "turbulenceModel.H"
 #include "dynamicFvMesh.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -84,7 +88,10 @@ int main(int argc, char *argv[])
         // --- PIMPLE loop
         for (int ocorr=0; ocorr<nOuterCorr; ocorr++)
         {
-            p.storePrevIter();
+            if (nOuterCorr != 1)
+            {
+                p.storePrevIter();
+            }
 
 #           include "UEqn.H"
 
@@ -112,7 +119,11 @@ int main(int argc, char *argv[])
 
                     pEqn.setReference(pRefCell, pRefValue);
 
-                    if (corr == nCorr-1 && nonOrth == nNonOrthCorr)
+                    if
+                    (
+                        ocorr == nOuterCorr-1
+                     && corr == nCorr-1
+                     && nonOrth == nNonOrthCorr)
                     {
                         pEqn.solve(mesh.solver(p.name() + "Final"));
                     }
@@ -143,6 +154,8 @@ int main(int argc, char *argv[])
             }
         }
 
+        turbulence->correct();
+
         runTime.write();
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"

applications/solvers/incompressible/pimpleFoam/Make/options

@@ -1,10 +1,12 @@
 EXE_INC = \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
-    -I$(LIB_SRC)/turbulenceModels/RAS \
+    -I$(LIB_SRC)/transportModels \
-    -I$(LIB_SRC)/transportModels
+    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
+    -I$(LIB_SRC)/finiteVolume/lnInclude
 
 EXE_LIBS = \
-    -lincompressibleRASModels \
     -lincompressibleTransportModels \
+    -lincompressibleRASModels \
+    -lincompressibleLESModels \
     -lfiniteVolume \
     -lmeshTools

applications/solvers/incompressible/pimpleFoam/createFields.H

@@ -36,7 +36,7 @@ setRefCell(p, mesh.solutionDict().subDict("PIMPLE"), pRefCell, pRefValue);
 
 singlePhaseTransportModel laminarTransport(U, phi);
 
-autoPtr<incompressible::RASModel> turbulence
+autoPtr<incompressible::turbulenceModel> turbulence
 (
-    incompressible::RASModel::New(U, phi, laminarTransport)
+    incompressible::turbulenceModel::New(U, phi, laminarTransport)
 );

applications/solvers/incompressible/pimpleFoam/pimpleFoam.C

@@ -26,14 +26,16 @@ Application
     pimpleFoam
 
 Description
-    Large time-step transient solver for incompressible, turbulent flow using
+    Large time-step transient solver for incompressible, flow using the PIMPLE
-    the PIMPLE (merged PISO-SIMPLE) algorithm.
+    (merged PISO-SIMPLE) algorithm.
+
+    Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
 
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "incompressible/singlePhaseTransportModel/singlePhaseTransportModel.H"
+#include "singlePhaseTransportModel.H"
-#include "incompressible/RASModel/RASModel.H"
+#include "turbulenceModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -58,14 +60,14 @@ int main(int argc, char *argv[])
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        if (nOuterCorr != 1)
-        {
-            p.storePrevIter();
-        }
-
         // --- Pressure-velocity PIMPLE corrector loop
         for (int oCorr=0; oCorr<nOuterCorr; oCorr++)
         {
+            if (nOuterCorr != 1)
+            {
+                p.storePrevIter();
+            }
+
             #include "UEqn.H"
 
             // --- PISO loop

applications/solvers/incompressible/pisoFoam/Make/files

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

applications/solvers/incompressible/pisoFoam/Make/options

@@ -1,10 +1,12 @@
 EXE_INC = \
-    -I$(LIB_SRC)/turbulenceModels/RAS \
+    -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
     -I$(LIB_SRC)/transportModels \
+    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
     -I$(LIB_SRC)/finiteVolume/lnInclude
 
 EXE_LIBS = \
     -lincompressibleRASModels \
+    -lincompressibleLESModels \
     -lincompressibleTransportModels \
     -lfiniteVolume \
     -lmeshTools

applications/solvers/incompressible/pisoFoam/createFields.H

@@ -12,7 +12,6 @@
         mesh
     );
 
-
     Info<< "Reading field U\n" << endl;
     volVectorField U
     (
@@ -37,7 +36,7 @@
 
     singlePhaseTransportModel laminarTransport(U, phi);
 
-    autoPtr<incompressible::LESModel> sgsModel
+    autoPtr<incompressible::turbulenceModel> turbulence
     (
-        incompressible::LESModel::New(U, phi, laminarTransport)
+        incompressible::turbulenceModel::New(U, phi, laminarTransport)
     );

applications/solvers/incompressible/pisoFoam/pisoFoam.C

@@ -26,13 +26,15 @@ Application
     turbFoam
 
 Description
-    Transient solver for incompressible, turbulent flow.
+    Transient solver for incompressible flow.
+
+    Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
 
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "incompressible/singlePhaseTransportModel/singlePhaseTransportModel.H"
+#include "singlePhaseTransportModel.H"
-#include "incompressible/RASModel/RASModel.H"
+#include "turbulenceModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -80,7 +82,7 @@ int main(int argc, char *argv[])
                 volScalarField rUA = 1.0/UEqn.A();
 
                 U = rUA*UEqn.H();
-                phi = (fvc::interpolate(U) & mesh.Sf()) 
+                phi = (fvc::interpolate(U) & mesh.Sf())
                     + fvc::ddtPhiCorr(rUA, U, phi);
 
                 adjustPhi(phi, U, p);

applications/solvers/incompressible/simpleFoam/Make/options

@@ -1,7 +1,9 @@
 EXE_INC = \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
-    -I$(LIB_SRC)/turbulenceModels/RAS \
+    -I$(LIB_SRC)/turbulenceModels/incompressible/RAS/RASModel \
-    -I$(LIB_SRC)/transportModels
+    -I$(LIB_SRC)/transportModels \
+    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
+    -I$(LIB_SRC)/finiteVolume/lnInclude
 
 EXE_LIBS = \
     -lincompressibleRASModels \

applications/solvers/incompressible/simpleFoam/simpleFoam.C

@@ -31,8 +31,8 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "incompressible/singlePhaseTransportModel/singlePhaseTransportModel.H"
+#include "singlePhaseTransportModel.H"
-#include "incompressible/RASModel/RASModel.H"
+#include "RASModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

applications/solvers/incompressible/turbDyMFoam/Make/files

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

applications/solvers/incompressible/turbDyMFoam/correctPhi.H

@@ -1,44 +0,0 @@
-{
-    wordList pcorrTypes(p.boundaryField().types());
-
-    for (label i=0; i<p.boundaryField().size(); i++)
-    {
-        if(p.boundaryField()[i].fixesValue())
-        {
-            pcorrTypes[i] = fixedValueFvPatchScalarField::typeName;
-        }
-    }
-
-    volScalarField pcorr
-    (
-        IOobject
-        (
-            "pcorr",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh,
-        dimensionedScalar("pcorr", p.dimensions(), 0.0),
-        pcorrTypes
-    );
-
-    for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
-    {
-        fvScalarMatrix pcorrEqn
-        (
-            fvm::laplacian(rAU, pcorr) == fvc::div(phi)
-        );
-
-        pcorrEqn.setReference(pRefCell, pRefValue);
-        pcorrEqn.solve();
-
-        if (nonOrth == nNonOrthCorr)
-        {
-            phi -= pcorrEqn.flux();
-        }
-    }
-}
-
-#include "continuityErrs.H"

applications/solvers/incompressible/turbDyMFoam/readControls.H

@@ -1,14 +0,0 @@
-#   include "readTimeControls.H"
-#   include "readPISOControls.H"
-
-    bool correctPhi = false;
-    if (piso.found("correctPhi"))
-    {
-        correctPhi = Switch(piso.lookup("correctPhi"));
-    }
-
-    bool checkMeshCourantNo = false;
-    if (piso.found("checkMeshCourantNo"))
-    {
-        checkMeshCourantNo = Switch(piso.lookup("checkMeshCourantNo"));
-    }

applications/solvers/incompressible/turbDyMFoam/turbDyMFoam.C

@@ -1,163 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
-     \\/     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
-    turbDyMFoam
-
-Description
-    Transient solver for incompressible, turbulent flow of Newtonian fluids
-    with moving mesh.
-
-\*---------------------------------------------------------------------------*/
-
-#include "fvCFD.H"
-#include "incompressible/singlePhaseTransportModel/singlePhaseTransportModel.H"
-#include "incompressible/RASModel/RASModel.H"
-#include "dynamicFvMesh.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-int main(int argc, char *argv[])
-{
-
-#   include "setRootCase.H"
-
-#   include "createTime.H"
-#   include "createDynamicFvMesh.H"
-#   include "readPISOControls.H"
-#   include "initContinuityErrs.H"
-#   include "createFields.H"
-#   include "readTimeControls.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-    Info<< "\nStarting time loop\n" << endl;
-
-    while (runTime.run())
-    {
-#       include "readControls.H"
-#       include "CourantNo.H"
-
-        p.storePrevIter();
-
-        // Make the fluxes absolute
-        fvc::makeAbsolute(phi, U);
-
-#       include "setDeltaT.H"
-
-        runTime++;
-
-        Info<< "Time = " << runTime.timeName() << nl << endl;
-
-        bool meshChanged = mesh.update();
-
-        if (correctPhi && meshChanged)
-        {
-#           include "correctPhi.H"
-        }
-
-        // Make the fluxes relative to the mesh motion
-        fvc::makeRelative(phi, U);
-
-        if (meshChanged && checkMeshCourantNo)
-        {
-#           include "meshCourantNo.H"
-        }
-
-        // --- PIMPLE loop
-        for (int ocorr=0; ocorr<nOuterCorr; ocorr++)
-        {
-#           include "UEqn.H"
-
-            // --- PISO loop
-            for (int corr=0; corr<nCorr; corr++)
-            {
-                rAU = 1.0/UEqn.A();
-
-                U = rAU*UEqn.H();
-                phi = (fvc::interpolate(U) & mesh.Sf());
-                  //+ fvc::ddtPhiCorr(rAU, U, phi);
-
-                adjustPhi(phi, U, p);
-
-                for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
-                {
-                    fvScalarMatrix pEqn
-                    (
-                        fvm::laplacian(rAU, p) == fvc::div(phi)
-                    );
-
-                    pEqn.setReference(pRefCell, pRefValue);
-
-                    if
-                    (
-                        ocorr == nOuterCorr-1
-                     && corr == nCorr-1
-                     && nonOrth == nNonOrthCorr)
-                    {
-                        pEqn.solve(mesh.solver(p.name() + "Final"));
-                    }
-                    else
-                    {
-                        pEqn.solve(mesh.solver(p.name()));
-                    }
-                
-                    if (nonOrth == nNonOrthCorr)
-                    {
-                        phi -= pEqn.flux();
-                    }
-                }
-
-#               include "continuityErrs.H"
-
-                // Explicitly relax pressure for momentum corrector
-                if (ocorr != nOuterCorr-1)
-                {
-                    p.relax();
-                }
-
-                // Make the fluxes relative to the mesh motion
-                fvc::makeRelative(phi, U);
-
-                U -= rAU*fvc::grad(p);
-                U.correctBoundaryConditions();
-            }
-        }
-
-        turbulence->correct();
-
-        runTime.write();
-
-        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
-            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
-            << nl << endl;
-    }
-
-    Info<< "End\n" << endl;
-
-    return(0);
-}
-
-
-// ************************************************************************* //

applications/solvers/incompressible/turbFoam/Make/files

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

applications/solvers/incompressible/turbFoam/createFields.H

@@ -1,42 +0,0 @@
-    Info<< "Reading field p\n" << endl;
-    volScalarField p
-    (
-        IOobject
-        (
-            "p",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-    Info<< "Reading field U\n" << endl;
-    volVectorField U
-    (
-        IOobject
-        (
-            "U",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-#   include "createPhi.H"
-
-
-    label pRefCell = 0;
-    scalar pRefValue = 0.0;
-    setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);
-
-
-    singlePhaseTransportModel laminarTransport(U, phi);
-
-    autoPtr<incompressible::RASModel> turbulence
-    (
-        incompressible::RASModel::New(U, phi, laminarTransport)
-    );

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,6 +14,8 @@
       - fvc::div(muEff*(fvc::interpolate(dev(fvc::grad(U))) & mesh.Sf()))
     );
 
+    UEqn.relax();
+
     if (momentumPredictor)
     {
         solve(UEqn == -fvc::grad(p));

applications/solvers/multiphase/cavitatingFoam/cavitatingFoam.C

@@ -23,17 +23,19 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Application
-    lesCavitatingFoam
+    cavitatingFoam
 
 Description
-    Transient cavitation code with LES turbulence.
+    Transient cavitation code based on the barotropic equation of state.
+
+    Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
 
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
 #include "barotropicCompressibilityModel.H"
 #include "twoPhaseMixture.H"
-#include "incompressible/LESModel/LESModel.H"
+#include "turbulenceModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -64,8 +66,6 @@ int main(int argc, char *argv[])
         runTime++;
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        turbulence->correct();
-
         for (int outerCorr=0; outerCorr<nOuterCorr; outerCorr++)
         {
 #           include "rhoEqn.H"
@@ -78,6 +78,8 @@ int main(int argc, char *argv[])
             }
         }
 
+        turbulence->correct();
+
         runTime.write();
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"

applications/solvers/multiphase/cavitatingFoam/createFields.H

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

applications/solvers/multiphase/compressibleInterDyMFoam/Make/files

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

applications/solvers/multiphase/compressibleInterDyMFoam/Make/options

@@ -0,0 +1,22 @@
+INTERFOAM = $(FOAM_SOLVERS)/multiphase/interFoam
+
+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 \
+    -I$(LIB_SRC)/dynamicMesh/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/dynamicFvMesh/lnInclude
+
+EXE_LIBS = \
+    -linterfaceProperties \
+    -lincompressibleTransportModels \
+    -lincompressibleRASModels \
+    -lincompressibleLESModels \
+    -lfiniteVolume \
+    -ldynamicMesh \
+    -lmeshTools \
+    -ldynamicFvMesh
+

applications/solvers/multiphase/compressibleInterDyMFoam/UEqn.H

@@ -1,6 +1,6 @@
     surfaceScalarField muf =
         twoPhaseProperties.muf()
-      + fvc::interpolate(rho*turbulence->nuSgs());
+      + fvc::interpolate(rho*turbulence->nut());
 
     fvVectorMatrix UEqn
     (
@@ -11,6 +11,8 @@
     //- fvc::div(muf*(mesh.Sf() & fvc::interpolate(fvc::grad(U)().T())))
     );
 
+    UEqn.relax();
+
     if (momentumPredictor)
     {
         solve

applications/solvers/multiphase/compressibleInterDyMFoam/alphaEqnsSubCycle.H

@@ -10,7 +10,7 @@
     );
 
     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);
 

applications/solvers/multiphase/compressibleInterDyMFoam/compressibleInterDyMFoam.C

@@ -23,23 +23,25 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Application
-    rasInterFoam
+    compressibleLesInterFoam
 
 Description
-    Solver for 2 incompressible, isothermal immiscible fluids using a VOF
+    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
+    momentum equation is solved.
-    selectable incompressible RAS model.
+
+    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 "incompressible/RASModel/RASModel.H"
+#include "turbulenceModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -47,52 +49,84 @@ int main(int argc, char *argv[])
 {
     #include "setRootCase.H"
     #include "createTime.H"
-    #include "createMesh.H"
+    #include "createDynamicFvMesh.H"
     #include "readEnvironmentalProperties.H"
-    #include "readPISOControls.H"
+    #include "readControls.H"
     #include "initContinuityErrs.H"
     #include "createFields.H"
-    #include "readTimeControls.H"
-    #include "correctPhi.H"
     #include "CourantNo.H"
     #include "setInitialDeltaT.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
     Info<< "\nStarting time loop\n" << endl;
 
     while (runTime.run())
     {
-        #include "readPISOControls.H"
+        #include "readControls.H"
-        #include "readTimeControls.H"
         #include "CourantNo.H"
+
+        // Make the fluxes absolute
+        fvc::makeAbsolute(phi, U);
+
         #include "setDeltaT.H"
 
         runTime++;
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        #include "gammaEqnSubCycle.H"
+        scalar timeBeforeMeshUpdate = runTime.elapsedCpuTime();
 
-        #include "UEqn.H"
+        // Do any mesh changes
+        mesh.update();
 
-        // --- PISO loop
+        if (mesh.changing())
-        for (int corr=0; corr < nCorr; corr++)
         {
-            #include "pEqn.H"
+            Info<< "Execution time for mesh.update() = "
+                << runTime.elapsedCpuTime() - timeBeforeMeshUpdate
+                << " s" << endl;
+
+            gh = g & mesh.C();
+            ghf = g & mesh.Cf();
         }
 
-        #include "continuityErrs.H"
+        if (mesh.changing() && correctPhi)
+        {
+            //***HGW#include "correctPhi.H"
+        }
 
-        p = pd + rho*gh;
+        // 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"
+        Info<< "ExecutionTime = "
-            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << runTime.elapsedCpuTime()
-            << nl << endl;
+            << " s\n\n" << endl;
     }
 
     Info<< "End\n" << endl;

applications/solvers/multiphase/compressibleInterDyMFoam/createFields.H

@@ -145,8 +145,8 @@
     // Construct interface from alpha1 distribution
     interfaceProperties interface(alpha1, U, twoPhaseProperties);
 
-    // Construct LES model
+    // Construct incompressible turbulence model
-    autoPtr<incompressible::LESModel> turbulence
+    autoPtr<incompressible::turbulenceModel> turbulence
     (
-        incompressible::LESModel::New(U, phi, twoPhaseProperties)
+        incompressible::turbulenceModel::New(U, phi, twoPhaseProperties)
     );

applications/solvers/multiphase/compressibleInterDyMFoam/pEqn.H

@@ -0,0 +1,77 @@
+{
+    volScalarField rUA = 1.0/UEqn.A();
+    surfaceScalarField rUAf = fvc::interpolate(rUA);
+
+    tmp<fvScalarMatrix> pdEqnComp;
+
+    if (transonic)
+    {
+        pdEqnComp =
+            (fvm::ddt(pd) + fvm::div(phi, pd) - fvm::Sp(fvc::div(phi), pd));
+    }
+    else
+    {
+        pdEqnComp =
+            (fvm::ddt(pd) + fvc::div(phi, pd) - fvc::Sp(fvc::div(phi), pd));
+    }
+
+
+    U = rUA*UEqn.H();
+
+    surfaceScalarField phiU
+    (
+        "phiU",
+        (fvc::interpolate(U) & mesh.Sf()) + fvc::ddtPhiCorr(rUA, rho, U, phi)
+    );
+
+    phi = phiU +
+        (
+            fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
+          - ghf*fvc::snGrad(rho)
+        )*rUAf*mesh.magSf();
+
+    for(int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    {
+        fvScalarMatrix pdEqnIncomp
+        (
+            fvc::div(phi)
+          - fvm::laplacian(rUAf, pd)
+        );
+
+        solve
+        (
+            (
+                max(alpha1, scalar(0))*(psi1/rho1)
+              + max(alpha2, scalar(0))*(psi2/rho2)
+            )
+           *pdEqnComp()
+          + pdEqnIncomp
+        );
+
+        if (nonOrth == nNonOrthCorr)
+        {
+            dgdt =
+                (pos(alpha2)*(psi2/rho2) - pos(alpha1)*(psi1/rho1))
+               *(pdEqnComp & pd);
+            phi += pdEqnIncomp.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
+        );
+
+    rho1 = rho10 + psi1*p;
+    rho2 = rho20 + psi2*p;
+
+    Info<< "max(U) " << max(mag(U)).value() << endl;
+    Info<< "min(pd) " << min(pd).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

@@ -4,12 +4,12 @@ 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/incompressible/turbulenceModel \
-    -I$(LIB_SRC)/turbulenceModels/LES/LESdeltas/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude
 
 EXE_LIBS = \
     -linterfaceProperties \
     -lincompressibleTransportModels \
+    -lincompressibleRASModels \
     -lincompressibleLESModels \
     -lfiniteVolume

applications/solvers/multiphase/compressibleInterFoam/UEqn.H

@@ -14,6 +14,8 @@
     //- fvc::div(muEff*(fvc::interpolate(dev(fvc::grad(U))) & mesh.Sf()))
     );
 
+    UEqn.relax();
+
     if (momentumPredictor)
     {
         solve
@@ -23,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/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(oneField(), 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

@@ -29,8 +29,9 @@ 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
+    momentum equation is solved.
-    selectable incompressible LES model.
+
+    Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
 
 \*---------------------------------------------------------------------------*/
 
@@ -39,7 +40,7 @@ Description
 #include "subCycle.H"
 #include "interfaceProperties.H"
 #include "twoPhaseMixture.H"
-#include "incompressible/LESModel/LESModel.H"
+#include "turbulenceModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -69,8 +70,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 +88,8 @@ int main(int argc, char *argv[])
 
         rho = alpha1*rho1 + alpha2*rho2;
 
+        turbulence->correct();
+
         runTime.write();
 
         Info<< "ExecutionTime = "

applications/solvers/multiphase/compressibleInterFoam/createFields.H

@@ -0,0 +1,152 @@
+    Info<< "Reading field pd\n" << endl;
+    volScalarField pd
+    (
+        IOobject
+        (
+            "pd",
+            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<< "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);
+
+    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 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;
+
+    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/compressibleLesInterFoam/Make/files

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

applications/solvers/multiphase/interDyMFoam/Make/options

@@ -1,11 +1,9 @@
 EXE_INC = \
-    -I../rasInterFoam \
     -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)/turbulenceModels/RAS/incompressible/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
@@ -16,6 +14,7 @@ EXE_LIBS = \
     -linterfaceProperties \
     -lincompressibleTransportModels \
     -lincompressibleRASModels \
+    -lincompressibleLESModels \
     -lfiniteVolume \
     -ldynamicMesh \
     -lmeshTools \

applications/solvers/multiphase/interDyMFoam/createFields.H

@@ -12,12 +12,12 @@
         mesh
     );
 
-    Info<< "Reading field gamma\n" << endl;
+    Info<< "Reading field alpha1\n" << endl;
-    volScalarField gamma
+    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);
 
     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,
+        alpha1*rho1 + (scalar(1) - alpha1)*rho2,
-        gamma.boundaryField().types()
+        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,13 +83,13 @@
     );
 
 
-    // Construct interface from gamma distribution
+    // Construct interface from alpha1 distribution
-    interfaceProperties interface(gamma, U, twoPhaseProperties);
+    interfaceProperties interface(alpha1, U, twoPhaseProperties);
 
-    // Construct incompressible RAS model
+    // Construct incompressible turbulence model
-    autoPtr<incompressible::RASModel> turbulence
+    autoPtr<incompressible::turbulenceModel> turbulence
     (
-        incompressible::RASModel::New(U, phi, twoPhaseProperties)
+        incompressible::turbulenceModel::New(U, phi, twoPhaseProperties)
     );
 
     wordList pcorrTypes(pd.boundaryField().types());

applications/solvers/multiphase/interDyMFoam/interDyMFoam.C

@@ -39,7 +39,7 @@ Description
 #include "subCycle.H"
 #include "interfaceProperties.H"
 #include "twoPhaseMixture.H"
-#include "incompressible/RASModel/RASModel.H"
+#include "turbulenceModel.H"
 #include "probes.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -58,8 +58,7 @@ int main(int argc, char *argv[])
     #include "CourantNo.H"
     #include "setInitialDeltaT.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
     Info<< "\nStarting time loop\n" << endl;
 
     while (runTime.run())
@@ -106,7 +105,7 @@ int main(int argc, char *argv[])
 
         twoPhaseProperties.correct();
 
-        #include "gammaEqnSubCycle.H"
+        #include "alphaEqnSubCycle.H"
 
         #include "UEqn.H"
 

applications/solvers/multiphase/interDyMFoam/pEqn.H

@@ -7,7 +7,7 @@
 
     phi = phiU +
     (
-        fvc::interpolate(interface.sigmaK())*fvc::snGrad(gamma)
+        fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
       - ghf*fvc::snGrad(rho)
     )*rAUf*mesh.magSf();
 

applications/solvers/multiphase/interFoam/Make/options

@@ -2,9 +2,12 @@ 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

applications/solvers/multiphase/interFoam/UBlendingFactor.H

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

applications/solvers/multiphase/interFoam/UEqn.H

@@ -1,14 +1,21 @@
-    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)
+      - fvm::laplacian(muEff, U)
-      - (fvc::grad(U) & fvc::grad(muf))
+      - (fvc::grad(U) & fvc::grad(muEff))
-    //- fvc::div(muf*(fvc::interpolate(dev(fvc::grad(U))) & mesh.Sf()))
+    //- fvc::div(muEff*(fvc::interpolate(dev(fvc::grad(U))) & mesh.Sf()))
     );
 
+    UEqn.relax();
+
     if (momentumPredictor)
     {
         solve
@@ -18,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 gCorr=0; gCorr<nAlphaCorr; gCorr++)
+    {
+        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/createFields.H

@@ -12,12 +12,12 @@
         mesh
     );
 
-    Info<< "Reading field gamma\n" << endl;
+    Info<< "Reading field alpha1\n" << endl;
-    volScalarField gamma
+    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);
 
     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,
+        alpha1*rho1 + (scalar(1) - alpha1)*rho2,
-        gamma.boundaryField().types()
+        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
@@ -107,5 +107,12 @@
     setRefCell(pd, mesh.solutionDict().subDict("PISO"), pdRefCell, pdRefValue);
 
 
-    // Construct interface from gamma distribution
+    // Construct interface from alpha1 distribution
-    interfaceProperties interface(gamma, U, twoPhaseProperties);
+    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

@@ -31,6 +31,8 @@ Description
     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 +42,7 @@ Description
 #include "subCycle.H"
 #include "interfaceProperties.H"
 #include "twoPhaseMixture.H"
+#include "turbulenceModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -57,7 +60,7 @@ int main(int argc, char *argv[])
     #include "CourantNo.H"
     #include "setInitialDeltaT.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 
@@ -74,7 +77,7 @@ int main(int argc, char *argv[])
 
         twoPhaseProperties.correct();
 
-        #include "gammaEqnSubCycle.H"
+        #include "alphaEqnSubCycle.H"
 
         #include "UEqn.H"
 
@@ -88,6 +91,8 @@ int main(int argc, char *argv[])
 
         p = pd + rho*gh;
 
+        turbulence->correct();
+
         runTime.write();
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"

applications/solvers/multiphase/interFoam/pEqn.H

@@ -12,7 +12,7 @@
 
     phi = phiU +
         (
-            fvc::interpolate(interface.sigmaK())*fvc::snGrad(gamma)
+            fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
           - ghf*fvc::snGrad(rho)
         )*rUAf*mesh.magSf();
 

applications/solvers/multiphase/interPhaseChangeFoam/Make/options

@@ -2,13 +2,13 @@ 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/incompressible/turbulenceModel \
-    -I$(LIB_SRC)/turbulenceModels/LES/LESdeltas/lnInclude \
     -IphaseChangeTwoPhaseMixtures/phaseChangeTwoPhaseMixture \
     -I$(LIB_SRC)/finiteVolume/lnInclude
 
 EXE_LIBS = \
     -linterfaceProperties \
     -lincompressibleTransportModels \
+    -lincompressibleRASModels \
     -lincompressibleLESModels \
     -lfiniteVolume

applications/solvers/multiphase/interPhaseChangeFoam/UEqn.H

@@ -1,6 +1,6 @@
     surfaceScalarField muf =
         twoPhaseProperties->muf()
-      + fvc::interpolate(rho*turbulence->nuSgs());
+      + fvc::interpolate(rho*turbulence->nut());
 
     fvVectorMatrix UEqn
     (
@@ -23,7 +23,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 alphaScheme("div(phi,alpha)");
-    word gammarScheme("div(phirb,gamma)");
+    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,
+                alpha1,
-                gammaScheme
+                alphaScheme
             )
           + fvc::flux
             (
-                -fvc::flux(-phir, scalar(1) - gamma, gammarScheme),
+                -fvc::flux(-phir, scalar(1) - alpha1, alpharScheme),
-                gamma,
+                alpha1,
-                gammarScheme
+                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(alpha1, phi, phiAlpha, 1, 0);
-        //MULES::explicitSolve(oneField(), gamma, phi, phiGamma, Sp, Su, 1, 0);
+        //MULES::explicitSolve(oneField(), alpha1, phi, phiAlpha, Sp, Su, 1, 0);
-        MULES::implicitSolve(oneField(), gamma, phi, phiGamma, Sp, Su, 1, 0);
+        MULES::implicitSolve(oneField(), 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()
+        << alpha1.weightedAverage(mesh.V()).value()
-        << "  Min(gamma) = " << min(gamma).value()
+        << "  Min(alpha1) = " << min(alpha1).value()
-        << "  Max(gamma) = " << max(gamma).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);
+            subCycle<volScalarField> alphaSubCycle(alpha1, nAlphaSubCycles);
-            !(++gammaSubCycle).end();
+            !(++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/createFields.H

@@ -12,12 +12,12 @@
         mesh
     );
 
-    Info<< "Reading field gamma\n" << endl;
+    Info<< "Reading field alpha1\n" << endl;
-    volScalarField gamma
+    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);
 
     const dimensionedScalar& rho1 = twoPhaseProperties->rho1();
     const dimensionedScalar& rho2 = twoPhaseProperties->rho2();
@@ -60,8 +60,8 @@
             mesh,
             IOobject::READ_IF_PRESENT
         ),
-        gamma*rho1 + (scalar(1) - gamma)*rho2,
+        alpha1*rho1 + (scalar(1) - alpha1)*rho2,
-        gamma.boundaryField().types()
+        alpha1.boundaryField().types()
     );
     rho.oldTime();
 
@@ -88,11 +88,11 @@
     );
 
 
-    // Construct interface from gamma distribution
+    // Construct interface from alpha1 distribution
-    interfaceProperties interface(gamma, U, twoPhaseProperties());
+    interfaceProperties interface(alpha1, U, twoPhaseProperties());
 
-    // Construct LES model
+    // Construct incompressible turbulence model
-    autoPtr<incompressible::LESModel> turbulence
+    autoPtr<incompressible::turbulenceModel> turbulence
     (
-        incompressible::LESModel::New(U, phi, twoPhaseProperties())
+        incompressible::turbulenceModel::New(U, phi, twoPhaseProperties())
     );

applications/solvers/multiphase/interPhaseChangeFoam/interPhaseChangeFoam.C

@@ -35,6 +35,8 @@ Description
     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 +44,7 @@ Description
 #include "subCycle.H"
 #include "interfaceProperties.H"
 #include "phaseChangeTwoPhaseMixture.H"
-#include "incompressible/LESModel/LESModel.H"
+#include "turbulenceModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -60,7 +62,7 @@ int main(int argc, char *argv[])
     #include "CourantNo.H"
     #include "setInitialDeltaT.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 
@@ -75,9 +77,7 @@ int main(int argc, char *argv[])
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        twoPhaseProperties->correct();
+        #include "alphaEqnSubCycle.H"
-
-        #include "gammaEqnSubCycle.H"
 
         turbulence->correct();
 
@@ -95,6 +95,8 @@ int main(int argc, char *argv[])
             #include "continuityErrs.H"
         }
 
+        twoPhaseProperties->correct();
+
         runTime.write();
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"

applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H

@@ -13,7 +13,7 @@
 
     phi = phiU +
         (
-            fvc::interpolate(interface.sigmaK())*fvc::snGrad(gamma)
+            fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
           - ghf*fvc::snGrad(rho)
         )*rUAf*mesh.magSf();
 

applications/solvers/multiphase/interPhaseChangeFoam/phaseChangeTwoPhaseMixtures/phaseChangeTwoPhaseMixture/phaseChangeTwoPhaseMixture.H

@@ -106,7 +106,7 @@ public:
         (
             const volVectorField& U,
             const surfaceScalarField& phi,
-            const word& alpha1Name
+            const word& alpha1Name = "alpha1"
         );