Added generic turbulenceModel base class to incompressible turbulence models.

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/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/channelOodles/channelOodles.C

@@ -1,155 +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 for flow in a channel.
-
-\*---------------------------------------------------------------------------*/
-
-#include "fvCFD.H"
-#include "incompressible/singlePhaseTransportModel/singlePhaseTransportModel.H"
-#include "incompressible/LESModel/LESModel.H"
-#include "IFstream.H"
-#include "OFstream.H"
-#include "Random.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-int main(int argc, char *argv[])
-{
-    #include "setRootCase.H"
-    #include "createTime.H"
-    #include "createMesh.H"
-    #include "readTransportProperties.H"
-    #include "createFields.H"
-    #include "initContinuityErrs.H"
-    #include "createGradP.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)
-         ==
-            flowDirection*gradP
-        );
-
-        if (momentumPredictor)
-        {
-            solve(UEqn == -fvc::grad(p));
-        }
-
-
-        // --- PISO loop
-
-        volScalarField rUA = 1.0/UEqn.A();
-
-        for (int corr=0; corr<nCorr; corr++)
-        {
-            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();
-        }
-
-
-        // Correct driving force for a constant mass flow rate
-
-        // Extract the velocity in the flow direction
-        dimensionedScalar magUbarStar =
-            (flowDirection & U)().weightedAverage(mesh.V());
-
-        // Calculate the pressure gradient increment needed to
-        // adjust the average flow-rate to the correct value
-        dimensionedScalar gragPplus =
-            (magUbar - magUbarStar)/rUA.weightedAverage(mesh.V());
-
-        U += flowDirection*rUA*gragPplus;
-
-        gradP += gragPplus;
-
-        Info<< "Uncorrected Ubar = " << magUbarStar.value() << tab
-            << "pressure gradient = " << gradP.value() << endl;
-
-        runTime.write();
-
-        #include "writeGradP.H"
-
-        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
-            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
-            << nl << endl;
-    }
-
-    Info<< "End\n" << endl;
-
-    return(0);
-}
-
-
-// ************************************************************************* //

applications/solvers/incompressible/channelOodles/createFields.H

@@ -1,43 +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::LESModel> sgsModel
-    (
-        incompressible::LESModel::New(U, phi, laminarTransport)
-    );

applications/solvers/incompressible/channelOodles/createGradP.H

@@ -1,24 +0,0 @@
-    dimensionedScalar gradP
-    (
-        "gradP",
-        dimensionSet(0, 1, -2, 0, 0),
-        0.0
-    );
-
-
-    IFstream gradPFile
-    (
-        runTime.path()/runTime.timeName()/"uniform"/"gradP.raw"
-    );
-
-    if(gradPFile.good())
-    {
-        gradPFile >> gradP;
-        Info<< "Reading average pressure gradient" <<endl
-            << endl;
-    }
-    else
-    {
-        Info<< "Initializing with 0 pressure gradient" <<endl
-            << endl;
-    };

applications/solvers/incompressible/channelOodles/readTransportProperties.H

@@ -1,28 +0,0 @@
-    Info<< "\nReading transportProperties\n" << endl;
-    IOdictionary transportProperties
-    (
-        IOobject
-        (
-            "transportProperties",
-            runTime.constant(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE
-        )
-    );
-
-
-    dimensionedScalar nu
-    (
-        transportProperties.lookup("nu")
-    );
-
-
-    //  Read centerline velocity for channel simulations
-    dimensionedVector Ubar
-    (
-        transportProperties.lookup("Ubar")
-    );
-
-    dimensionedScalar magUbar = mag(Ubar);
-    vector flowDirection = (Ubar/magUbar).value();

applications/solvers/incompressible/channelOodles/writeGradP.H

@@ -1,19 +0,0 @@
-    if (runTime.outputTime())
-    {
-        OFstream gradPFile
-        (
-            runTime.path()/runTime.timeName()/"uniform"/"gradP.raw"
-        );
-
-        if(gradPFile.good())
-        {
-            gradPFile << gradP << endl;
-        }
-        else
-        {
-            FatalErrorIn(args.executable())
-                << "Cannot open file "
-                << runTime.path()/runTime.timeName()/"uniform"/"gradP.raw"
-                << exit(FatalError);
-        };
-    };

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/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/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/icoDyMFoam/icoDyMFoam.C

@@ -1,159 +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
-    icoDyMFoam
-
-Description
-    Transient solver for incompressible, laminar flow of Newtonian fluids
-    with moving mesh.
-
-\*---------------------------------------------------------------------------*/
-
-#include "fvCFD.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"
-
-        // Make the fluxes absolute
-        fvc::makeAbsolute(phi, U);
-
-#       include "setDeltaT.H"
-
-        runTime++;
-
-        Info<< "Time = " << runTime.timeName() << nl << endl;
-
-        mesh.update();
-
-        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"
-        }
-
-        // --- PIMPLE loop
-        for (int ocorr=0; ocorr<nOuterCorr; ocorr++)
-        {
-            p.storePrevIter();
-
-#           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());
-
-                if (p.needReference())
-                {
-                    fvc::makeRelative(phi, U);
-                    adjustPhi(phi, U, p);
-                    fvc::makeAbsolute(phi, U);
-                }
-
-                for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
-                {
-                    fvScalarMatrix pEqn
-                    (
-                        fvm::laplacian(rAU, 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"
-
-                // 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();
-            }
-        }
-
-        runTime.write();
-
-        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
-            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
-            << nl << endl;
-    }
-
-    Info<< "End\n" << endl;
-
-    return(0);
-}
-
-
-// ************************************************************************* //

applications/solvers/incompressible/icoDyMFoam/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/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/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)/meshTools/lnInclude \
-    -I$(LIB_SRC)/sampling/lnInclude
-
-EXE_LIBS = \
-    -lincompressibleLESModels \
-    -lincompressibleTransportModels \
-    -lfiniteVolume \
-    -lmeshTools

applications/solvers/incompressible/oodles/createFields.H

@@ -1,43 +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::LESModel> sgsModel
-    (
-        incompressible::LESModel::New(U, phi, laminarTransport)
-    );

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/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;
 
+        // --- Pressure-velocity PIMPLE corrector loop
+        for (int oCorr=0; oCorr<nOuterCorr; oCorr++)
+        {
             if (nOuterCorr != 1)
             {
                 p.storePrevIter();
             }
 
-        // --- Pressure-velocity PIMPLE corrector loop
-        for (int oCorr=0; oCorr<nOuterCorr; oCorr++)
-        {
             #include "UEqn.H"
 
             // --- PISO loop

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/Make/options

@@ -1,15 +0,0 @@
-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)/transportModels \
-    -I$(LIB_SRC)/finiteVolume/lnInclude
-
-EXE_LIBS = \
-    -ldynamicFvMesh \
-    -ldynamicMesh \
-    -lmeshTools \
-    -lincompressibleRASModels \
-    -lincompressibleTransportModels \
-    -lfiniteVolume

applications/solvers/incompressible/turbDyMFoam/UEqn.H

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

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/createFields.H

@@ -1,59 +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)
-    );
-
-    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/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/Make/options

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

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/incompressible/turbFoam/turbFoam.C

@@ -1,129 +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
-    turbFoam
-
-Description
-    Transient solver for incompressible, turbulent flow.
-
-\*---------------------------------------------------------------------------*/
-
-#include "fvCFD.H"
-#include "incompressible/singlePhaseTransportModel/singlePhaseTransportModel.H"
-#include "incompressible/RASModel/RASModel.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-int main(int argc, char *argv[])
-{
-
-#   include "setRootCase.H"
-
-#   include "createTime.H"
-#   include "createMesh.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"
-
-        // Pressure-velocity PISO corrector
-        {
-            // Momentum predictor
-
-            fvVectorMatrix UEqn
-            (
-                fvm::ddt(U)
-              + fvm::div(phi, U)
-              + turbulence->divDevReff(U)
-            );
-
-            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);
-
-                // Non-orthogonal pressure corrector loop
-                for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
-                {
-                    // Pressure corrector
-
-                    fvScalarMatrix pEqn
-                    (
-                        fvm::laplacian(rUA, p) == fvc::div(phi)
-                    );
-
-                    pEqn.setReference(pRefCell, pRefValue);
-                    pEqn.solve();
-
-                    if (nonOrth == nNonOrthCorr)
-                    {
-                        phi -= pEqn.flux();
-                    }
-                }
-
-#               include "continuityErrs.H"
-
-                U -= rUA*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/multiphase/compressibleLesInterFoam/Make/files

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

applications/solvers/multiphase/compressibleLesInterFoam/Make/options

@@ -1,15 +0,0 @@
-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/LES \
-    -I$(LIB_SRC)/turbulenceModels/LES/LESdeltas/lnInclude \
-    -I$(LIB_SRC)/finiteVolume/lnInclude
-
-EXE_LIBS = \
-    -linterfaceProperties \
-    -lincompressibleTransportModels \
-    -lincompressibleLESModels \
-    -lfiniteVolume

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/compressibleLesInterFoam/alphaEqns.H

@@ -1,76 +0,0 @@
-{
-    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/compressibleLesInterFoam/alphaEqnsSubCycle.H

@@ -1,43 +0,0 @@
-{
-    label nAlphaCorr
-    (
-        readLabel(piso.lookup("nAlphaCorr"))
-    );
-
-    label nAlphaSubCycles
-    (
-        readLabel(piso.lookup("nAlphaSubCycles"))
-    );
-
-    surfaceScalarField phic = mag(phi/mesh.magSf());
-    phic = min(interface.cGamma()*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/compressibleLesInterFoam/compressibleLesInterFoam.C

@@ -1,105 +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
-    compressibleLesInterFoam
-
-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.
-
-\*---------------------------------------------------------------------------*/
-
-#include "fvCFD.H"
-#include "MULES.H"
-#include "subCycle.H"
-#include "interfaceProperties.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 "readEnvironmentalProperties.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"
-        #include "setDeltaT.H"
-
-        runTime++;
-
-        Info<< "Time = " << runTime.timeName() << nl << endl;
-
-        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/compressibleLesInterFoam/createFields.H

@@ -1,152 +0,0 @@
-    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 LES model
-    autoPtr<incompressible::LESModel> turbulence
-    (
-        incompressible::LESModel::New(U, phi, twoPhaseProperties)
-    );

applications/solvers/multiphase/compressibleLesInterFoam/pEqn.H

@@ -1,74 +0,0 @@
-{
-    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;
-}

applications/solvers/multiphase/compressibleLesInterFoam/readControls.H

@@ -1,20 +0,0 @@
-   #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);
-    }

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/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/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/gammaEqn.H

@@ -1,67 +0,0 @@
-{
-    word gammaScheme("div(phi,gamma)");
-    word gammarScheme("div(phirb,gamma)");
-
-    surfaceScalarField phir("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
-            );
-
-        Pair<tmp<volScalarField> > vDotAlphal =
-            twoPhaseProperties->vDotAlphal();
-        const volScalarField& vDotcAlphal = vDotAlphal[0]();
-        const volScalarField& vDotvAlphal = vDotAlphal[1]();
-
-        volScalarField Sp
-        (
-            IOobject
-            (
-                "Sp",
-                runTime.timeName(),
-                mesh
-            ),
-            vDotvAlphal - vDotcAlphal
-        );
-
-        volScalarField Su
-        (
-            IOobject
-            (
-                "Su",
-                runTime.timeName(),
-                mesh
-            ),
-            // Divergence term is handled explicitly to be
-            // consistent with the explicit transport solution
-            divU*gamma
-          + 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);
-
-        rhoPhi +=
-            (runTime.deltaT()/totalDeltaT)
-           *(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/interPhaseChangeFoam/gammaEqnSubCycle.H

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

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"
         );
 
 

applications/solvers/multiphase/lesCavitatingFoam/CourantNo.H

@@ -1,59 +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
-
-Global
-    CourantNo
-
-Description
-    Calculates and outputs the mean and maximum Courant Numbers.
-
-\*---------------------------------------------------------------------------*/
-
-scalar CoNum = 0.0;
-scalar meanCoNum = 0.0;
-scalar acousticCoNum = 0.0;
-
-if (mesh.nInternalFaces())
-{
-    surfaceScalarField SfUfbyDelta = 
-        mesh.surfaceInterpolation::deltaCoeffs()*mag(phiv);
-
-    CoNum = max(SfUfbyDelta/mesh.magSf())
-        .value()*runTime.deltaT().value();
-
-    meanCoNum = (sum(SfUfbyDelta)/sum(mesh.magSf()))
-        .value()*runTime.deltaT().value();
-
-    acousticCoNum = max
-    (
-        mesh.surfaceInterpolation::deltaCoeffs()/sqrt(fvc::interpolate(psi))
-    ).value()*runTime.deltaT().value();
-}
-
-Info<< "phiv Courant Number mean: " << meanCoNum
-    << " max: " << CoNum
-    << " acoustic max: " << acousticCoNum
-    << endl;
-
-// ************************************************************************* //

applications/solvers/multiphase/lesCavitatingFoam/Make/files

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

applications/solvers/multiphase/lesCavitatingFoam/Make/options

@@ -1,16 +0,0 @@
-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
-
-EXE_LIBS = \
-    -lincompressibleTransportModels \
-    -lincompressibleLESModels \
-    -lfiniteVolume \
-    -lbarotropicCompressibilityModel
-

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/continuityErrs.H

@@ -1,22 +0,0 @@
-{
-    volScalarField thermoRho = psi*p + (1.0 - gamma)*rhol0;
-
-    dimensionedScalar totalMass = fvc::domainIntegrate(rho);
-
-    scalar sumLocalContErr = 
-    (
-        fvc::domainIntegrate(mag(rho - thermoRho))/totalMass
-    ).value();
-
-    scalar globalContErr = 
-    (
-        fvc::domainIntegrate(rho - thermoRho)/totalMass
-    ).value();
-
-    cumulativeContErr += globalContErr;
-
-    Info<< "time step continuity errors : sum local = " << sumLocalContErr
-        << ", global = " << globalContErr
-        << ", cumulative = " << cumulativeContErr
-        << endl;
-}

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 (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
-
-\*---------------------------------------------------------------------------*/
-
-#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/devOneEqEddy/devOneEqEddy.H

@@ -1,148 +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
-
-Class
-    devOneEqEddy
-
-Description
-<pre>
-    One Equation Eddy Viscosity Model
-    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-    Eddy viscosity SGS model using a modeled balance equation to simulate the
-    behaviour of k, hence,
-
-        d/dt(k) + div(U*k) - div(nuEff*grad(k))
-        =
-        -B*L - ce*k^3/2/delta
-
-    and
-
-        B = 2/3*k*I - 2*nuEff*dev(D)
-
-    where
-
-        D = symm(grad(U));
-        nuSgs = ck*sqrt(k)*delta
-        nuEff = nuSgs + nu
-</pre>
-
-SourceFiles
-    devOneEqEddy.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef devOneEqEddy_H
-#define devOneEqEddy_H
-
-#include "GenEddyVisc.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace incompressible
-{
-namespace LESModels
-{
-
-/*---------------------------------------------------------------------------*\
-                           Class devOneEqEddy Declaration
-\*---------------------------------------------------------------------------*/
-
-class devOneEqEddy
-:
-    public GenEddyVisc
-{
-    // Private data
-
-        volScalarField k_;
-
-        dimensionedScalar ck_;
-
-
-    // Private Member Functions
-
-        // Disallow default bitwise copy construct and assignment
-        devOneEqEddy(const devOneEqEddy&);
-        devOneEqEddy& operator=(const devOneEqEddy&);
-
-
-public:
-
-    //- Runtime type information
-    TypeName("devOneEqEddy");
-
-    // Constructors
-
-        //- Constructor from components
-        devOneEqEddy
-        (
-            const volVectorField& U,
-            const surfaceScalarField& phi,
-            transportModel& transport
-        );
-
-
-    // Destructor
-
-        ~devOneEqEddy()
-        {}
-
-
-    // Member Functions
-
-        //- Return SGS kinetic energy
-        tmp<volScalarField> k() const
-        {
-            return k_;
-        }
-
-        //- Return the effective diffusivity for k
-        tmp<volScalarField> DkEff() const
-        {
-            return tmp<volScalarField>
-            (
-                new volScalarField("DkEff", nuSgs_ + nu())
-            );
-        }
-
-        //- Correct Eddy-Viscosity and related properties
-        void correct(const tmp<volTensorField>& gradU);
-
-        //- Read turbulenceProperties dictionary
-        bool read();
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace LESModelsModels
-} // End namespace incompressible
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

applications/solvers/multiphase/lesCavitatingFoam/gammaPsi.H

@@ -1,10 +0,0 @@
-{
-    gamma = max(min((rho - rholSat)/(rhovSat - rholSat), scalar(1)), scalar(0));
-
-    Info<< "max-min gamma: " << max(gamma).value()
-        << " " << min(gamma).value() << endl;
-
-    psiModel->correct();
-
-    //Info<< "min a: " << 1.0/sqrt(max(psi)).value() << endl;
-}

applications/solvers/multiphase/lesCavitatingFoam/lesCavitatingFoam.C

@@ -1,94 +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
-    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/lesCavitatingFoam/readControls.H

@@ -1,9 +0,0 @@
-#include "readTimeControls.H"
-
-scalar maxAcousticCo
-(
-    readScalar(runTime.controlDict().lookup("maxAcousticCo"))
-);
-
-
-#include "readPISOControls.H"

applications/solvers/multiphase/lesCavitatingFoam/readThermodynamicProperties.H

@@ -1,27 +0,0 @@
-    Info<< "Reading thermodynamicProperties\n" << endl;
-
-    IOdictionary thermodynamicProperties
-    (
-        IOobject
-        (
-            "thermodynamicProperties",
-            runTime.constant(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE
-        )
-    );
-
-    dimensionedScalar psil(thermodynamicProperties.lookup("psil"));
-
-    dimensionedScalar rholSat(thermodynamicProperties.lookup("rholSat"));
-
-    dimensionedScalar psiv(thermodynamicProperties.lookup("psiv"));
-
-    dimensionedScalar pSat(thermodynamicProperties.lookup("pSat"));
-
-    dimensionedScalar rhovSat("rhovSat", psiv*pSat);
-
-    dimensionedScalar rhol0("rhol0", rholSat - pSat*psil);
-
-    dimensionedScalar rhoMin(thermodynamicProperties.lookup("rhoMin"));

applications/solvers/multiphase/lesCavitatingFoam/resetPhiPatches.H

@@ -1,15 +0,0 @@
-fvsPatchScalarFieldField& phiPatches = phi.boundaryField();
-const fvPatchScalarFieldField& rhoPatches = rho.boundaryField();
-const fvPatchVectorFieldField& Upatches = U.boundaryField();
-const fvsPatchVectorFieldField& SfPatches = mesh.Sf().boundaryField();
-
-forAll(phiPatches, patchI)
-{
-    if (phi.boundaryField().types()[patchI] == "calculated")
-    {
-        calculatedFvsPatchScalarField& phiPatch =
-            refCast<calculatedFvsPatchScalarField>(phiPatches[patchI]);
-
-        phiPatch == ((rhoPatches[patchI]*Upatches[patchI]) & SfPatches[patchI]);
-    }
-}

applications/solvers/multiphase/lesCavitatingFoam/resetPhivPatches.H

@@ -1,14 +0,0 @@
-surfaceScalarField::GeometricBoundaryField& phivPatches = phiv.boundaryField();
-const volVectorField::GeometricBoundaryField& Upatches = U.boundaryField();
-const surfaceVectorField::GeometricBoundaryField& SfPatches = mesh.Sf().boundaryField();
-
-forAll(phivPatches, patchI)
-{
-    if (phiv.boundaryField().types()[patchI] == "calculated")
-    {
-        calculatedFvsPatchScalarField& phivPatch =
-            refCast<calculatedFvsPatchScalarField>(phivPatches[patchI]);
-
-        phivPatch == (Upatches[patchI] & SfPatches[patchI]);
-    }
-}

applications/solvers/multiphase/lesCavitatingFoam/rhoEqn.H

@@ -1,16 +0,0 @@
-{
-    fvScalarMatrix rhoEqn
-    (
-        fvm::ddt(rho)
-      + fvm::div(phiv, rho)
-    );
-
-    rhoEqn.solve();
-
-    phi = rhoEqn.flux();
-
-    Info<< "max-min rho: " << max(rho).value()
-        << " " << min(rho).value() << endl;
-
-    rho == max(rho, rhoMin);
-}

applications/solvers/multiphase/lesCavitatingFoam/setDeltaT.H

@@ -1,54 +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
-
-Global
-    setDeltaT
-
-Description
-    Reset the timestep to maintain a constant maximum courant Number.
-    Reduction of time-step is imediate but increase is damped to avoid
-    unstable oscillations.
-
-\*---------------------------------------------------------------------------*/
-
-if (adjustTimeStep)
-{
-    scalar maxDeltaTFact =
-        min(maxCo/(CoNum + SMALL), maxAcousticCo/(acousticCoNum + SMALL));
-
-    scalar deltaTFact = min(min(maxDeltaTFact, 1.0 + 0.1*maxDeltaTFact), 1.2);
-
-    runTime.setDeltaT
-    (
-        min
-        (
-            deltaTFact*runTime.deltaT().value(),
-            maxDeltaT
-        )
-    );
-    
-    Info<< "deltaT = " <<  runTime.deltaT().value() << endl;
-}
-
-// ************************************************************************* //

applications/solvers/multiphase/lesCavitatingFoam/setInitialDeltaT.H

@@ -1,54 +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
-
-Global
-    setInitialDeltaT
-
-Description
-    Set the initial timestep corresponding to the timestep adjustment
-    algorithm in setDeltaT
-
-\*---------------------------------------------------------------------------*/
-
-if (adjustTimeStep)
-{
-#   include "CourantNo.H"
-
-    if (CoNum > SMALL)
-    {
-        scalar maxDeltaTFact =
-            min(maxCo/(CoNum + SMALL), maxAcousticCo/(acousticCoNum + SMALL));
-
-        runTime.setDeltaT
-        (
-            min
-            (
-                maxDeltaTFact*runTime.deltaT().value(),
-                maxDeltaT
-            )
-        );
-    }
-}
-
-// ************************************************************************* //

applications/solvers/multiphase/lesInterFoam/Make/files

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

applications/solvers/multiphase/lesInterFoam/Make/options

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

applications/solvers/multiphase/lesInterFoam/UEqn.H

@@ -1,32 +0,0 @@
-    surfaceScalarField muEff
-    (
-        "muEff",
-        twoPhaseProperties.muf()
-      + fvc::interpolate(rho*turbulence->nuSgs())
-    );
-
-    fvVectorMatrix UEqn
-    (
-        fvm::ddt(rho, U)
-      + 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()))
-    );
-
-    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/lesInterFoam/createFields.H

@@ -1,115 +0,0 @@
-    Info<< "Reading field pd\n" << endl;
-    volScalarField pd
-    (
-        IOobject
-        (
-            "pd",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-    Info<< "Reading field gamma\n" << endl;
-    volScalarField gamma
-    (
-        IOobject
-        (
-            "gamma",
-            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"
-
-    Info<< "Reading transportProperties\n" << endl;
-    twoPhaseMixture twoPhaseProperties(U, phi, "gamma");
-
-    const dimensionedScalar& rho1 = twoPhaseProperties.rho1();
-    const dimensionedScalar& rho2 = twoPhaseProperties.rho2();
-
-
-    // Need to store rho for ddt(rho, U)
-    volScalarField rho
-    (
-        IOobject
-        (
-            "rho",
-            runTime.timeName(),
-            mesh,
-            IOobject::READ_IF_PRESENT
-        ),
-        gamma*rho1 + (scalar(1) - gamma)*rho2,
-        gamma.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.
-    surfaceScalarField rhoPhi
-    (
-        IOobject
-        (
-            "rho*phi",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        rho1*phi
-    );
-
-
-    label pdRefCell = 0;
-    scalar pdRefValue = 0.0;
-    setRefCell(pd, mesh.solutionDict().subDict("PISO"), pdRefCell, pdRefValue);
-
-
-    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
-    );
-
-
-    // Construct interface from gamma distribution
-    interfaceProperties interface(gamma, U, twoPhaseProperties);
-
-    // Construct LES model
-    autoPtr<incompressible::LESModel> turbulence
-    (
-        incompressible::LESModel::New(U, phi, twoPhaseProperties)
-    );

applications/solvers/multiphase/lesInterFoam/lesInterFoam.C

@@ -1,103 +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
-    lesInterFoam
-
-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.
-
-\*---------------------------------------------------------------------------*/
-
-#include "fvCFD.H"
-#include "MULES.H"
-#include "subCycle.H"
-#include "interfaceProperties.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 "readEnvironmentalProperties.H"
-    #include "readPISOControls.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 "readTimeControls.H"
-        #include "CourantNo.H"
-        #include "setDeltaT.H"
-
-        runTime++;
-        Info<< "Time = " << runTime.timeName() << nl << endl;
-
-        #include "gammaEqnSubCycle.H"
-
-        turbulence->correct();
-
-        #include "UEqn.H"
-
-        // --- PISO loop
-        for (int corr=0; corr < nCorr; corr++)
-        {
-            #include "pEqn.H"
-        }
-
-        #include "continuityErrs.H"
-
-        p = pd + rho*gh;
-
-        runTime.write();
-
-        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
-            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
-            << nl << endl;
-    }
-
-    Info<< "End\n" << endl;
-
-    return(0);
-}
-
-
-// ************************************************************************* //

applications/solvers/multiphase/multiphaseInterFoam/Make/options

@@ -6,9 +6,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 \
+    -lincompressibleTransportModels \
+    -lincompressibleRASModels \
+    -lincompressibleLESModels \
     -lfiniteVolume

applications/solvers/multiphase/multiphaseInterFoam/UEqn.H

@@ -1,10 +1,34 @@
-    surfaceScalarField muf = mixture.muf();
+    surfaceScalarField muEff
+    (
+        "muEff",
+        mixture.muf()
+      + fvc::interpolate(rho*turbulence->nut())
+    );
 
     fvVectorMatrix UEqn
     (
         fvm::ddt(rho, U)
       + fvm::div(mixture.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*(mesh.Sf() & fvc::interpolate(fvc::grad(U)().T())))
+    //- fvc::div(muEff*(fvc::interpolate(dev(fvc::grad(U))) & mesh.Sf()))
     );
+
+    UEqn.relax();
+
+    if (momentumPredictor)
+    {
+        solve
+        (
+            UEqn
+         ==
+            fvc::reconstruct
+            (
+                (
+                    mixture.surfaceTensionForce()
+                  - ghf*fvc::snGrad(rho)
+                  - fvc::snGrad(pd)
+                ) * mesh.magSf()
+            )
+        );
+    }

applications/solvers/multiphase/multiphaseInterFoam/createFields.H

@@ -46,9 +46,31 @@
 
 
     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);
+
+
+    // Construct incompressible turbulence model
+    autoPtr<incompressible::turbulenceModel> turbulence
+    (
+        incompressible::turbulenceModel::New(U, phi, mixture)
+    );

applications/solvers/multiphase/multiphaseInterFoam/multiphaseInterFoam.C

@@ -29,16 +29,18 @@ Description
     Solver for n incompressible fluids which captures the interfaces and
     includes surface-tension and contact-angle effects for each.
 
+    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"
@@ -51,7 +53,7 @@ int main(int argc, char *argv[])
     #include "CourantNo.H"
     #include "setInitialDeltaT.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
 
@@ -79,6 +81,10 @@ int main(int argc, char *argv[])
 
         #include "continuityErrs.H"
 
+        p = pd + rho*gh;
+
+        turbulence->correct();
+
         runTime.write();
 
         Info<< "ExecutionTime = "

applications/solvers/multiphase/rasCavitatingFoam/CourantNo.H

@@ -1,59 +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
-
-Global
-    CourantNo
-
-Description
-    Calculates and outputs the mean and maximum Courant Numbers.
-
-\*---------------------------------------------------------------------------*/
-
-scalar CoNum = 0.0;
-scalar meanCoNum = 0.0;
-scalar acousticCoNum = 0.0;
-
-if (mesh.nInternalFaces())
-{
-    surfaceScalarField SfUfbyDelta = 
-        mesh.surfaceInterpolation::deltaCoeffs()*mag(phiv);
-
-    CoNum = max(SfUfbyDelta/mesh.magSf())
-        .value()*runTime.deltaT().value();
-
-    meanCoNum = (sum(SfUfbyDelta)/sum(mesh.magSf()))
-        .value()*runTime.deltaT().value();
-
-    acousticCoNum = max
-    (
-        mesh.surfaceInterpolation::deltaCoeffs()/sqrt(fvc::interpolate(psi))
-    ).value()*runTime.deltaT().value();
-}
-
-Info<< "phiv Courant Number mean: " << meanCoNum
-    << " max: " << CoNum
-    << " acoustic max: " << acousticCoNum
-    << endl;
-
-// ************************************************************************* //

applications/solvers/multiphase/rasCavitatingFoam/Make/files

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

applications/solvers/multiphase/rasCavitatingFoam/Make/options

@@ -1,14 +0,0 @@
-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/RAS \
-    -I$(LIB_SRC)/thermophysicalModels/barotropicCompressibilityModel/lnInclude
-
-EXE_LIBS = \
-    -lincompressibleTransportModels \
-    -lincompressibleRASModels \
-    -lfiniteVolume \
-    -lbarotropicCompressibilityModel
-

applications/solvers/multiphase/rasCavitatingFoam/UEqn.H

@@ -1,20 +0,0 @@
-    surfaceScalarField muEff
-    (
-        "muEff",
-        twoPhaseProperties.muf()
-      + fvc::interpolate(rho*turbulence->nut())
-    );
-
-    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/rasCavitatingFoam/continuityErrs.H

@@ -1,22 +0,0 @@
-{
-    volScalarField thermoRho = psi*p + (1.0 - gamma)*rhol0;
-
-    dimensionedScalar totalMass = fvc::domainIntegrate(rho);
-
-    scalar sumLocalContErr = 
-    (
-        fvc::domainIntegrate(mag(rho - thermoRho))/totalMass
-    ).value();
-
-    scalar globalContErr = 
-    (
-        fvc::domainIntegrate(rho - thermoRho)/totalMass
-    ).value();
-
-    cumulativeContErr += globalContErr;
-
-    Info<< "time step continuity errors : sum local = " << sumLocalContErr
-        << ", global = " << globalContErr
-        << ", cumulative = " << cumulativeContErr
-        << endl;
-}