Added compressibleLesInterFoam solver

applications/solvers/multiphase/compressibleLesInterFoam/Make/files

+compressibleLesInterFoam.C
+
+EXE = $(FOAM_USER_APPBIN)/compressibleLesInterFoam

applications/solvers/multiphase/compressibleLesInterFoam/Make/options

+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)/LESmodels \
+    -I$(LIB_SRC)/LESmodels/LESdeltas/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude
+
+EXE_LIBS = \
+    -linterfaceProperties \
+    -lincompressibleTransportModels \
+    -lincompressibleLESmodels \
+    -lfiniteVolume

applications/solvers/multiphase/compressibleLesInterFoam/UEqn.H

+    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

+{
+    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, 1.0)
+        );
+
+        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

+{
+    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

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2007 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

+    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, 0.0001);
+
+    // Construct interface from alpha1 distribution
+    interfaceProperties interface(alpha1, U, twoPhaseProperties);
+
+    // Construct LES model
+    autoPtr<LESmodel> turbulence
+    (
+        LESmodel::New(U, phi, twoPhaseProperties)
+    );

applications/solvers/multiphase/compressibleLesInterFoam/pEqn.H

+{
+    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, 0.0)*(psi1/rho1) + max(alpha2, 0.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

+   #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);
+    }