New solver: rhoPorousMRFPimpleFoam

Added compressibility support to MRFZone.
Required generalisation of the oneField concept.
Also clean-up the use of one and zero for consistency.
Simple mixerVessel2D tutorial case supplied.

applications/solvers/compressible/rhoPorousMRFPimpleFoam/Make/files

+rhoPorousMRFPimpleFoam.C
+
+EXE = $(FOAM_APPBIN)/rhoPorousMRFPimpleFoam
+

applications/solvers/compressible/rhoPorousMRFPimpleFoam/Make/options

+EXE_INC = \
+    -I../rhoPimpleFoam \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \
+    -I$(LIB_SRC)/finiteVolume/cfdTools \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude
+
+EXE_LIBS = \
+    -lbasicThermophysicalModels \
+    -lspecie \
+    -lcompressibleRASModels \
+    -lcompressibleLESModels \
+    -lfiniteVolume \
+    -lmeshTools

applications/solvers/compressible/rhoPorousMRFPimpleFoam/UEqn.H

+// Solve the Momentum equation
+
+tmp<fvVectorMatrix> UEqn
+(
+    pZones.ddt(rho, U)
+  + fvm::div(phi, U)
+  + turbulence->divDevRhoReff(U)
+);
+
+if (oCorr == nOuterCorr-1)
+{
+    UEqn().relax(1);
+}
+else
+{
+    UEqn().relax();
+}
+
+mrfZones.addCoriolis(rho, UEqn());
+pZones.addResistance(UEqn());
+
+volScalarField rUA = 1.0/UEqn().A();
+
+if (momentumPredictor)
+{
+    if (oCorr == nOuterCorr-1)
+    {
+        solve(UEqn() == -fvc::grad(p), mesh.solver("UFinal"));
+    }
+    else
+    {
+        solve(UEqn() == -fvc::grad(p));
+    }
+}
+else
+{
+    U = rUA*(UEqn().H() - fvc::grad(p));
+    U.correctBoundaryConditions();
+}

applications/solvers/compressible/rhoPorousMRFPimpleFoam/createFields.H

+    Info<< "Reading thermophysical properties\n" << endl;
+
+    autoPtr<basicPsiThermo> pThermo
+    (
+        basicPsiThermo::New(mesh)
+    );
+    basicPsiThermo& thermo = pThermo();
+
+    volScalarField& p = thermo.p();
+    volScalarField& h = thermo.h();
+    const volScalarField& psi = thermo.psi();
+
+    volScalarField rho
+    (
+        IOobject
+        (
+            "rho",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        thermo.rho()
+    );
+
+    Info<< "Reading field U\n" << endl;
+    volVectorField U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    #include "compressibleCreatePhi.H"
+
+    dimensionedScalar pMin
+    (
+        mesh.solutionDict().subDict("PIMPLE").lookup("pMin")
+    );
+
+    Info<< "Creating turbulence model\n" << endl;
+    autoPtr<compressible::turbulenceModel> turbulence
+    (
+        compressible::turbulenceModel::New
+        (
+            rho,
+            U,
+            phi,
+            thermo
+        )
+    );
+
+    //dimensionedScalar initialMass = fvc::domainIntegrate(rho);
+
+
+    Info<< "Creating field DpDt\n" << endl;
+    volScalarField DpDt =
+        fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
+
+    MRFZones mrfZones(mesh);
+    mrfZones.correctBoundaryVelocity(U);
+
+    porousZones pZones(mesh);
+    Switch pressureImplicitPorosity(false);

applications/solvers/compressible/rhoPorousMRFPimpleFoam/pEqn.H

+rho = thermo.rho();
+
+volScalarField rUA = 1.0/UEqn().A();
+U = rUA*UEqn().H();
+
+if (nCorr <= 1)
+{
+    UEqn.clear();
+}
+
+if (transonic)
+{
+    surfaceScalarField phid
+    (
+        "phid",
+        fvc::interpolate(psi)
+       *(
+            (fvc::interpolate(U) & mesh.Sf())
+          + fvc::ddtPhiCorr(rUA, rho, U, phi)
+        )
+    );
+    mrfZones.relativeFlux(fvc::interpolate(psi), phid);
+
+    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    {
+        fvScalarMatrix pEqn
+        (
+            fvm::ddt(psi, p)
+          + fvm::div(phid, p)
+          - fvm::laplacian(rho*rUA, p)
+        );
+
+        if
+        (
+            oCorr == nOuterCorr-1
+            && corr == nCorr-1
+            && nonOrth == nNonOrthCorr
+        )
+        {
+            pEqn.solve(mesh.solver("pFinal"));
+        }
+        else
+        {
+            pEqn.solve();
+        }
+
+        if (nonOrth == nNonOrthCorr)
+        {
+            phi == pEqn.flux();
+        }
+    }
+}
+else
+{
+    phi =
+        fvc::interpolate(rho)*
+        (
+            (fvc::interpolate(U) & mesh.Sf())
+        //+ fvc::ddtPhiCorr(rUA, rho, U, phi)
+        );
+    mrfZones.relativeFlux(fvc::interpolate(rho), phi);
+
+    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+    {
+        // Pressure corrector
+        fvScalarMatrix pEqn
+        (
+            fvm::ddt(psi, p)
+          + fvc::div(phi)
+          - fvm::laplacian(rho*rUA, p)
+        );
+
+        if
+        (
+            oCorr == nOuterCorr-1
+         && corr == nCorr-1
+         && nonOrth == nNonOrthCorr
+        )
+        {
+            pEqn.solve(mesh.solver("pFinal"));
+        }
+        else
+        {
+            pEqn.solve();
+        }
+
+        if (nonOrth == nNonOrthCorr)
+        {
+            phi += pEqn.flux();
+        }
+    }
+}
+
+#include "rhoEqn.H"
+#include "compressibleContinuityErrs.H"
+
+//if (oCorr != nOuterCorr-1)
+{
+    // Explicitly relax pressure for momentum corrector
+    p.relax();
+
+    rho = thermo.rho();
+    rho.relax();
+    Info<< "rho max/min : " << max(rho).value()
+        << " " << min(rho).value() << endl;
+}
+
+U -= rUA*fvc::grad(p);
+U.correctBoundaryConditions();
+
+DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p);
+
+bound(p, pMin);
+
+// For closed-volume cases adjust the pressure and density levels
+// to obey overall mass continuity
+/*
+if (closedVolume)
+{
+    p += (initialMass - fvc::domainIntegrate(psi*p))
+        /fvc::domainIntegrate(psi);
+}
+*/

applications/solvers/compressible/rhoPorousMRFPimpleFoam/rhoPorousMRFPimpleFoam.C

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2010 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
+    rhoPorousMRFPimpleFoam
+
+Description
+    Transient solver for laminar or turbulent flow of compressible fluids
+    with support for porous media and MRF for HVAC and similar applications.
+
+    Uses the flexible PIMPLE (PISO-SIMPLE) solution for time-resolved and
+    pseudo-transient simulations.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "basicPsiThermo.H"
+#include "turbulenceModel.H"
+#include "bound.H"
+#include "MRFZones.H"
+#include "porousZones.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;
+
+    while (runTime.run())
+    {
+        #include "readTimeControls.H"
+        #include "readPIMPLEControls.H"
+        #include "compressibleCourantNo.H"
+        #include "setDeltaT.H"
+
+        runTime++;
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        if (nOuterCorr != 1)
+        {
+            p.storePrevIter();
+            rho.storePrevIter();
+        }
+
+        #include "rhoEqn.H"
+
+        // --- Pressure-velocity PIMPLE corrector loop
+        for (int oCorr=0; oCorr<nOuterCorr; oCorr++)
+        {
+            #include "UEqn.H"
+            #include "hEqn.H"
+
+            // --- PISO loop
+            for (int corr=0; corr<nCorr; corr++)
+            {
+                #include "pEqn.H"
+            }
+
+            turbulence->correct();
+        }
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/multiphase/compressibleInterDyMFoam/alphaEqns.H

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

applications/solvers/multiphase/compressibleInterFoam/alphaEqns.H

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

applications/solvers/multiphase/interMixingFoam/alphaEqns.H

@@ -70,12 +70,12 @@
         MULES::limiter
         (
             allLambda,
-            oneField(),
+            geometricOneField(),
             alpha1,
             phiAlpha1BD,
             phiAlpha1,
-            zeroField(),
-            zeroField(),
+            zero(),
+            zero(),
             1,
             0,
             3
@@ -107,12 +107,12 @@
         MULES::limiter
         (
             allLambda,
-            oneField(),
+            geometricOneField(),
             alpha2,
             phiAlpha2BD,
             phiAlpha2,
-            zeroField(),
-            zeroField(),
+            zero(),
+            zero(),
             1,
             0,
             3

applications/solvers/multiphase/interPhaseChangeFoam/alphaEqn.H

@@ -51,8 +51,8 @@
         );
 
         //MULES::explicitSolve(alpha1, phi, phiAlpha, 1, 0);
-        //MULES::explicitSolve(oneField(), alpha1, phi, phiAlpha, Sp, Su, 1, 0);
-        MULES::implicitSolve(oneField(), alpha1, phi, phiAlpha, Sp, Su, 1, 0);
+        //MULES::explicitSolve(geometricOneField(), alpha1, phi, phiAlpha, Sp, Su, 1, 0);
+        MULES::implicitSolve(geometricOneField(), alpha1, phi, phiAlpha, Sp, Su, 1, 0);
 
         rhoPhi +=
             (runTime.deltaT()/totalDeltaT)

src/OpenFOAM/fields/FieldFields/oneFieldField/oneFieldField.H

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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
+    Foam::oneFieldField
+
+Description
+    A class representing the concept of a field of oneFields used to
+    avoid unnecessary manipulations for objects which are known to be one at
+    compile-time.
+
+    Used for example as the density argument to a function written for
+    compressible to be used for incompressible flow.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef oneFieldField_H
+#define oneFieldField_H
+
+#include "oneField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class oneField Declaration
+\*---------------------------------------------------------------------------*/
+
+class oneFieldField
+:
+    public one
+{
+
+public:
+
+    // Constructors
+
+        //- Construct null
+        oneFieldField()
+        {}
+
+
+    // Member Operators
+
+        inline oneField operator[](const label) const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#   include "oneFieldFieldI.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/OpenFOAM/fields/FieldFields/oneFieldField/oneFieldFieldI.H

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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 "oneFieldField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+inline Foam::oneField Foam::oneFieldField::operator[](const label) const
+{
+    return oneField();
+}
+
+
+// ************************************************************************* //

src/OpenFOAM/fields/FieldFields/zeroFieldField/zeroFieldField.H

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2009 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
+    Foam::zeroFieldField
+
+Description
+    A class representing the concept of a field of zeroFields used to
+    avoid unnecessary manipulations for objects which are known to be zero at
+    compile-time.
+
+    Used for example as the density argument to a function written for
+    compressible to be used for incompressible flow.
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef zeroFieldField_H