ENH: new SRFrhoSimpleFoam solver

META-INFO/api-info

@@ -1,2 +1,2 @@
 api=2012
-patch=210210
+patch=0

applications/solvers/combustion/PDRFoam/EaEqn.H

@@ -15,7 +15,7 @@
             )
           : -betav*dpdt
         )
-      - fvm::laplacian(DbThermal, hea)
+      - fvm::laplacian(Db, hea)
       + betav*fvOptions(rho, hea)
     );
 

applications/solvers/combustion/PDRFoam/EauEqn.H

@@ -16,7 +16,13 @@ if (ign.ignited())
             )*rho/thermo.rhou()
           : -betav*dpdt*rho/thermo.rhou()
         )
-      - fvm::laplacian(DbThermal, heau)
+      - fvm::laplacian(Db, heau)
+
+        // These terms cannot be used in partially-premixed combustion due to
+        // the resultant inconsistency between ft and heau transport.
+        // A possible solution would be to solve for ftu as well as ft.
+        //- fvm::div(muEff*fvc::grad(b)/(b + 0.001), heau)
+        //+ fvm::Sp(fvc::div(muEff*fvc::grad(b)/(b + 0.001)), heau)
      ==
         betav*fvOptions(rho, heau)
     );

applications/solvers/combustion/PDRFoam/Make/files

@@ -1,33 +1,33 @@
 XiModels/XiModel/XiModel.C
 XiModels/XiModel/XiModelNew.C
-
-XiModels/transport/transportTwoEqs/transportTwoEqs.C
+XiModels/fixed/fixed.C
+XiModels/algebraic/algebraic.C
+XiModels/transport/transport.C
 
 XiModels/XiEqModels/XiEqModel/XiEqModel.C
 XiModels/XiEqModels/XiEqModel/XiEqModelNew.C
-
 XiModels/XiEqModels/Gulder/Gulder.C
-XiModels/XiEqModels/instability2XiEq/instability2XiEq.C
-XiModels/XiEqModels/BLMgMaXiEq/BLMgMaXiEq.C
+XiModels/XiEqModels/instabilityXiEq/instabilityXiEq.C
+XiModels/XiEqModels/SCOPEBlendXiEq/SCOPEBlendXiEq.C
+XiModels/XiEqModels/SCOPEXiEq/SCOPEXiEq.C
 
 XiModels/XiGModels/XiGModel/XiGModel.C
 XiModels/XiGModels/XiGModel/XiGModelNew.C
 XiModels/XiGModels/KTS/KTS.C
 XiModels/XiGModels/instabilityG/instabilityG.C
-XiModels/XiGModels/instability2G/instability2G.C
 
 PDRModels/turbulence/PDRkEpsilon/PDRkEpsilon.C
 
 PDRModels/dragModels/PDRDragModel/PDRDragModel.C
 PDRModels/dragModels/PDRDragModel/PDRDragModelNew.C
 PDRModels/dragModels/basic/basic.C
-PDRModels/dragModels/basicSch/basicSch.C
 
-PDRModels/XiEqModels/normBasicXiSubXiEq/normBasicXiSubXiEq.C
-PDRModels/XiGModels/normBasicXiSubG/normBasicXiSubG.C
+PDRModels/XiEqModels/basicXiSubXiEq/basicXiSubXiEq.C
+PDRModels/XiGModels/basicXiSubG/basicXiSubG.C
 
 laminarFlameSpeed/SCOPE/SCOPELaminarFlameSpeed.C
 
+/*PDRFoamAutoRefine.C*/
 PDRFoam.C
 
 EXE = $(FOAM_APPBIN)/PDRFoam

applications/solvers/combustion/PDRFoam/PDRFoam.C

@@ -6,7 +6,6 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -82,11 +81,8 @@ Description
 #include "PDRDragModel.H"
 #include "ignition.H"
 #include "bound.H"
-#include "dynamicFvMesh.H"
-#include "dynamicRefineFvMesh.H"
 #include "pimpleControl.H"
 #include "fvOptions.H"
-#include "fvcSmooth.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -99,24 +95,22 @@ int main(int argc, char *argv[])
     );
 
     #include "postProcess.H"
+
     #include "addCheckCaseOptions.H"
     #include "setRootCaseLists.H"
     #include "createTime.H"
-
-    #include "createDynamicFvMesh.H"
-    #include "createDyMControls.H"
-    #include "initContinuityErrs.H"
-
+    #include "createMesh.H"
+    #include "createControl.H"
     #include "readCombustionProperties.H"
     #include "readGravitationalAcceleration.H"
     #include "createFields.H"
     #include "createFieldRefs.H"
-
+    #include "initContinuityErrs.H"
+    #include "createTimeControls.H"
     #include "compressibleCourantNo.H"
     #include "setInitialDeltaT.H"
 
     turbulence->validate();
-
     scalar StCoNum = 0.0;
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -125,62 +119,33 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
-        #include "readDyMControls.H"
-
+        #include "readTimeControls.H"
         #include "compressibleCourantNo.H"
         #include "setDeltaT.H"
 
         ++runTime;
-
         Info<< "\n\nTime = " << runTime.timeName() << endl;
 
+        #include "rhoEqn.H"
+
         // --- Pressure-velocity PIMPLE corrector loop
         while (pimple.loop())
         {
-            if (pimple.firstIter() || moveMeshOuterCorrectors)
-            {
-                fvc::makeAbsolute(phi, rho, U);
-
-                // Flux estimate for introduced faces.
-                volVectorField rhoU("rhoU", rho*U);
-
-                bool meshChanged = mesh.controlledUpdate();
-
-                if (runTime.write() && meshChanged)
-                {
-                    betav.write();
-                    betai.write();
-                    Blong.write();
-                    Bv.write();
-                    Lobs.write();
-                    CT.write();
-                    drag->writeFields();
-                    flameWrinkling->writeFields();
-                }
-
-                // Make the fluxes relative to the mesh motion
-                fvc::makeRelative(phi, rho, U);
-            }
-
-            if (pimple.firstIter())
-            {
-                #include "rhoEqn.H"
-            }
-
             #include "UEqn.H"
-            #include "bEqn.H"
-            #include "ftEqn.H"
-            #include "EauEqn.H"
-            #include "EaEqn.H"
-
-            if (!ign.ignited())
-            {
-                thermo.heu() == thermo.he();
-            }
 
             // --- Pressure corrector loop
             while (pimple.correct())
             {
+                #include "bEqn.H"
+                #include "ftEqn.H"
+                #include "EauEqn.H"
+                #include "EaEqn.H"
+
+                if (!ign.ignited())
+                {
+                    thermo.heu() == thermo.he();
+                }
+
                 #include "pEqn.H"
             }
 
@@ -190,14 +155,12 @@ int main(int argc, char *argv[])
             }
         }
 
-        rho = thermo.rho();
-
         runTime.write();
 
         runTime.printExecutionTime(Info);
     }
 
-    Info<< "\nEnd\n" << endl;
+    Info<< "\n end\n";
 
     return 0;
 }

applications/solvers/combustion/PDRFoam/PDRFoamAutoRefine.C

@@ -0,0 +1,216 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2016 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+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 3 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, see <http://www.gnu.org/licenses/>.
+
+Application
+    PDRFoam
+
+Description
+    Solver for compressible premixed/partially-premixed combustion with
+    turbulence modelling.
+
+    Combusting RANS code using the b-Xi two-equation model.
+    Xi may be obtained by either the solution of the Xi transport
+    equation or from an algebraic expression.  Both approaches are
+    based on Gulder's flame speed correlation which has been shown
+    to be appropriate by comparison with the results from the
+    spectral model.
+
+    Strain effects are incorporated directly into the Xi equation
+    but not in the algebraic approximation.  Further work need to be
+    done on this issue, particularly regarding the enhanced removal rate
+    caused by flame compression.  Analysis using results of the spectral
+    model will be required.
+
+    For cases involving very lean Propane flames or other flames which are
+    very strain-sensitive, a transport equation for the laminar flame
+    speed is present.  This equation is derived using heuristic arguments
+    involving the strain time scale and the strain-rate at extinction.
+    the transport velocity is the same as that for the Xi equation.
+
+    For large flames e.g. explosions additional modelling for the flame
+    wrinkling due to surface instabilities may be applied.
+
+    PDR (porosity/distributed resistance) modelling is included to handle
+    regions containing blockages which cannot be resolved by the mesh.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "dynamicFvMesh.H"
+#include "psiuReactionThermo.H"
+#include "turbulentFluidThermoModel.H"
+#include "laminarFlameSpeed.H"
+#include "XiModel.H"
+#include "PDRDragModel.H"
+#include "ignition.H"
+#include "bound.H"
+#include "dynamicRefineFvMesh.H"
+#include "pimpleControl.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    argList::addNote
+    (
+        "Solver for compressible premixed/partially-premixed combustion with"
+        " turbulence modelling."
+    );
+
+    #include "setRootCaseLists.H"
+    #include "createTime.H"
+    #include "createDynamicFvMesh.H"
+
+    pimpleControl pimple(mesh);
+
+    #include "readCombustionProperties.H"
+    #include "readGravitationalAcceleration.H"
+    #include "createFields.H"
+    #include "initContinuityErrs.H"
+    #include "createTimeControls.H"
+    #include "compressibleCourantNo.H"
+    #include "setInitialDeltaT.H"
+
+    turbulence->validate();
+    scalar StCoNum = 0.0;
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    bool hasChanged = false;
+
+    while (runTime.run())
+    {
+        #include "readTimeControls.H"
+        #include "compressibleCourantNo.H"
+        #include "setDeltaT.H"
+
+        // Indicators for refinement.
+        // Note: before ++runTime only for post-processing reasons.
+        tmp<volScalarField> tmagGradP = mag(fvc::grad(p));
+        volScalarField normalisedGradP
+        (
+            "normalisedGradP",
+            tmagGradP()/max(tmagGradP())
+        );
+        normalisedGradP.writeOpt() = IOobject::AUTO_WRITE;
+        tmagGradP.clear();
+
+        ++runTime;
+
+        Info<< "\n\nTime = " << runTime.timeName() << endl;
+
+        {
+            // Make the fluxes absolute
+            fvc::makeAbsolute(phi, rho, U);
+
+            // Test : disable refinement for some cells
+            bitSet& protectedCell =
+                refCast<dynamicRefineFvMesh>(mesh).protectedCell();
+
+            if (protectedCell.empty())
+            {
+                protectedCell.setSize(mesh.nCells());
+                protectedCell = false;
+            }
+
+            forAll(betav, celli)
+            {
+                if (betav[celli] < 0.99)
+                {
+                    protectedCell.set(celli);
+                }
+            }
+
+            // Flux estimate for introduced faces.
+            volVectorField rhoU("rhoU", rho*U);
+
+            // Do any mesh changes
+            bool meshChanged = mesh.update();
+
+
+            if (meshChanged)
+            {
+                hasChanged = true;
+            }
+
+            if (runTime.write() && hasChanged)
+            {
+                betav.write();
+                Lobs.write();
+                CT.write();
+                drag->writeFields();
+                flameWrinkling->writeFields();
+                hasChanged = false;
+            }
+
+            // Make the fluxes relative to the mesh motion
+            fvc::makeRelative(phi, rho, U);
+        }
+
+
+        #include "rhoEqn.H"
+
+        // --- Pressure-velocity PIMPLE corrector loop
+        while (pimple.loop())
+        {
+            #include "UEqn.H"
+
+
+            // --- Pressure corrector loop
+            while (pimple.correct())
+            {
+                #include "bEqn.H"
+                #include "ftEqn.H"
+                #include "huEqn.H"
+                #include "hEqn.H"
+
+                if (!ign.ignited())
+                {
+                    hu == h;
+                }
+
+                #include "pEqn.H"
+            }
+
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
+            }
+        }
+
+        runTime.write();
+
+        runTime.printExecutionTime(Info);
+    }
+
+    Info<< "\n end\n";
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/combustion/PDRFoam/PDRModels/XiEqModels/basicXiSubXiEq/basicXiSubXiEq.C

@@ -0,0 +1,150 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2016 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+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 3 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, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "basicXiSubXiEq.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace XiEqModels
+{
+    defineTypeNameAndDebug(basicSubGrid, 0);
+    addToRunTimeSelectionTable(XiEqModel, basicSubGrid, dictionary);
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::XiEqModels::basicSubGrid::basicSubGrid
+(
+    const dictionary& XiEqProperties,
+    const psiuReactionThermo& thermo,
+    const compressible::RASModel& turbulence,
+    const volScalarField& Su
+)
+:
+    XiEqModel(XiEqProperties, thermo, turbulence, Su),
+
+    B_
+    (
+        IOobject
+        (
+            "B",
+            Su.mesh().facesInstance(),
+            Su.mesh(),
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        ),
+        Su.mesh()
+    ),
+
+    XiEqModel_(XiEqModel::New(XiEqModelCoeffs_, thermo, turbulence, Su))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::XiEqModels::basicSubGrid::~basicSubGrid()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField> Foam::XiEqModels::basicSubGrid::XiEq() const
+{
+    const fvMesh& mesh = Su_.mesh();
+    const volVectorField& U = mesh.lookupObject<volVectorField>("U");
+
+    const volScalarField& Nv = mesh.lookupObject<volScalarField>("Nv");
+    const volSymmTensorField& nsv =
+        mesh.lookupObject<volSymmTensorField>("nsv");
+
+    volScalarField magU(mag(U));
+    volVectorField Uhat
+    (
+        U/(mag(U) + dimensionedScalar("Usmall", U.dimensions(), 1e-4))
+    );
+
+    const scalarField Cw = pow(mesh.V(), 2.0/3.0);
+
+    volScalarField N
+    (
+        IOobject
+        (
+            "N",
+            mesh.time().constant(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedScalar(Nv.dimensions(), Zero)
+    );
+    N.primitiveFieldRef() = Nv.primitiveField()*Cw;
+
+    volSymmTensorField ns
+    (
+        IOobject
+        (
+            "ns",
+            U.mesh().time().timeName(),
+            U.mesh(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        U.mesh(),
+        dimensionedSymmTensor(nsv.dimensions(), Zero)
+    );
+    ns.primitiveFieldRef() = nsv.primitiveField()*Cw;
+
+    volScalarField n(max(N - (Uhat & ns & Uhat), scalar(1e-4)));
+    volScalarField b((Uhat & B_ & Uhat)/sqrt(n));
+    volScalarField up(sqrt((2.0/3.0)*turbulence_.k()));
+
+    volScalarField XiSubEq
+    (
+        scalar(1)
+      + max(2.2*sqrt(b), min(0.34*magU/up*sqrt(b), scalar(1.6)))
+      * min(n, scalar(1))
+    );
+
+    return (XiSubEq*XiEqModel_->XiEq());
+}
+
+
+bool Foam::XiEqModels::basicSubGrid::read(const dictionary& XiEqProperties)
+{
+    XiEqModel::read(XiEqProperties);
+
+    return XiEqModel_->read(XiEqModelCoeffs_);
+}
+
+
+// ************************************************************************* //

applications/solvers/combustion/PDRFoam/PDRModels/XiEqModels/basicXiSubXiEq/basicXiSubXiEq.H

@@ -0,0 +1,147 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2012 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+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 3 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, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::XiEqModels::basicSubGrid
+
+Description
+    Basic sub-grid obstacle flame-wrinkling enhancement factor model.
+    Details supplied by J Puttock 2/7/06.
+
+    <b> Sub-grid flame area generation </b>
+
+    \f$ n = N - \hat{\dwea{\vec{U}}}.n_{s}.\hat{\dwea{\vec{U}}} \f$
+    \f$ n_{r} = \sqrt{n} \f$
+
+    where:
+
+        \f$ \hat{\dwea{\vec{U}}} = \dwea{\vec{U}} / \vert \dwea{\vec{U}}
+        \vert \f$
+
+        \f$ b = \hat{\dwea{\vec{U}}}.B.\hat{\dwea{\vec{U}}} / n_{r} \f$
+
+    where:
+
+        \f$ B \f$ is the file "B".
+
+        \f$ N \f$ is the file "N".
+
+        \f$  n_{s} \f$ is the file "ns".
+
+    The flame area enhancement factor \f$ \Xi_{sub} \f$ is expected to
+    approach:
+
+    \f[
+        \Xi_{{sub}_{eq}} =
+            1 + max(2.2 \sqrt{b}, min(0.34 \frac{\vert \dwea{\vec{U}}
+            \vert}{{\vec{U}}^{'}}, 1.6)) \times min(\frac{n}{4}, 1)
+    \f]
+
+
+SourceFiles
+    basicSubGrid.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef basicSubGrid_H
+#define basicSubGrid_H
+
+#include "XiEqModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace XiEqModels
+{
+
+/*---------------------------------------------------------------------------*\
+                          Class basicSubGrid Declaration
+\*---------------------------------------------------------------------------*/
+
+class basicSubGrid
+:
+    public XiEqModel
+{
+    // Private data
+
+        //- tblock
+        volSymmTensorField B_;
+
+        //- Equilibrium Xi model due to turbulence
+        autoPtr<XiEqModel> XiEqModel_;
+
+
+    // Private Member Functions
+
+        //- No copy construct
+        basicSubGrid(const basicSubGrid&) = delete;
+
+        //- No copy assignment
+        void operator=(const basicSubGrid&) = delete;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("basicSubGrid");
+
+
+    // Constructors
+
+        //- Construct from components
+        basicSubGrid
+        (
+            const dictionary& XiEqProperties,
+            const psiuReactionThermo& thermo,
+            const compressible::RASModel& turbulence,
+            const volScalarField& Su
+        );
+
+
+    //- Destructor
+    virtual ~basicSubGrid();
+
+
+    // Member Functions
+
+        //- Return the flame-wrinkling XiEq
+        virtual tmp<volScalarField> XiEq() const;
+
+        //- Update properties from given dictionary
+        virtual bool read(const dictionary& XiEqProperties);
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace XiEqModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/combustion/PDRFoam/PDRModels/XiEqModels/normBasicXiSubXiEq/normBasicXiSubXiEq.C

@@ -1,219 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2021 OpenCFD Ltd.
--------------------------------------------------------------------------------
-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 3 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, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "normBasicXiSubXiEq.H"
-#include "addToRunTimeSelectionTable.H"
-#include "ignition.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace XiEqModels
-{
-    defineTypeNameAndDebug(normBasicSubGrid, 0);
-    addToRunTimeSelectionTable(XiEqModel, normBasicSubGrid, dictionary);
-}
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::XiEqModels::normBasicSubGrid::normBasicSubGrid
-(
-    const dictionary& XiEqProperties,
-    const word& modelType,
-    const psiuReactionThermo& thermo,
-    const compressible::RASModel& turbulence,
-    const volScalarField& Su
-)
-:
-    XiEqModel(XiEqProperties, modelType,thermo, turbulence, Su),
-    Cxpe1_(XiEqModelCoeffs_.get<scalar>("Cxpe1")),
-    Cxpe2_(XiEqModelCoeffs_.get<scalar>("Cxpe2")),
-    Cxpe3_(XiEqModelCoeffs_.get<scalar>("Cxpe3")),
-    Cxpe4_(XiEqModelCoeffs_.get<scalar>("Cxpe4"))
-{}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::XiEqModels::normBasicSubGrid::~normBasicSubGrid()
-{}
-
-
-// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
-
-Foam::tmp<Foam::volScalarField> Foam::XiEqModels::normBasicSubGrid::XiEq() const
-{
-    const fvMesh& mesh = Su_.mesh();
-    const volVectorField& U = mesh.lookupObject<volVectorField>("U");
-    const volScalarField& b = mesh.lookupObject<volScalarField>("b");
-    const volScalarField& Nv = mesh.lookupObject<volScalarField>("Nv");
-    const volSymmTensorField& nsv =
-        mesh.lookupObject<volSymmTensorField>("nsv");
-    const volSymmTensorField& Bv =
-        mesh.lookupObject<volSymmTensorField>("Bv");
-
-    volScalarField magU(mag(U));
-
-    const scalarField Cw(pow(mesh.V(), 2.0/3.0));
-
-    auto tN = tmp<volScalarField>::New
-    (
-        IOobject
-        (
-            "tN",
-            mesh.time().constant(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh,
-        dimensionedScalar(Nv.dimensions(), Zero)
-    );
-    auto& N = tN.ref();
-
-    N.primitiveFieldRef() = Nv.primitiveField()*Cw;
-
-    auto tns = tmp<volSymmTensorField>::New
-    (
-        IOobject
-        (
-            "tns",
-            mesh.time().timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh,
-        dimensionedSymmTensor(nsv.dimensions(), Zero)
-    );
-    auto& ns = tns.ref();
-
-    auto tB = tmp<volSymmTensorField>::New
-    (
-        IOobject
-        (
-            "tB",
-            mesh.time().timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh,
-        dimensionedSymmTensor(Bv.dimensions(), Zero)
-    );
-    auto& B = tB.ref();
-
-    //calculating flame normal
-
-    volVectorField flNormal
-    (
-        "flNormal",
-        fvc::reconstruct(fvc::snGrad(b)*mesh.magSf())
-    );
-
-    volScalarField mgb("mgb", mag(flNormal));
-
-    dimensionedScalar dMgb("dMgb", mgb.dimensions(), SMALL);
-
-    const volScalarField bc(b*(1.0-b));
-
-    dMgb += 1.0e-8*
-        (bc*mgb)().weightedAverage(mesh.V())
-       /(bc.weightedAverage(mesh.V()) + SMALL);
-
-    mgb += dMgb;
-    flNormal /= mgb;
-
-    B.primitiveFieldRef() = Bv.primitiveField()*sqrt(Cw);
-    volScalarField Ntemp("Ntemp", N);
-    volScalarField Np("Np", max(N - (flNormal & ns & flNormal), scalar(1)));
-
-    // B_ is Bv*sqrt(Cw)
-    volScalarField bl("bl",(flNormal & B & flNormal)/sqrt(Np));
-    bl.min(1.0);
-
-    volScalarField up(sqrt((2.0/3.0)*turbulence_.k()));
-
-    IOdictionary combustionProperties
-    (
-        IOobject
-        (
-            "combustionProperties",
-            mesh.time().constant(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE
-        )
-    );
-
-
-    ignition ign(combustionProperties, mesh.time(), mesh);
-
-    dimensionedVector ignLoc("ignLoc", dimLength, ign.sites()[0].location());
-
-    dimensionedScalar filtRad2
-    (
-        "filtRad2",
-        dimLength,
-        6.0*ign.sites()[0].diameter()
-    );
-
-    const volScalarField filDist(mag(mesh.C() - ignLoc));
-
-    const volScalarField filterMult
-    (
-        pos(filDist - filtRad2)*neg(bl - 0.99)*pos(N - 1e-3)
-    );
-
-    tmp<volScalarField> XiSubEq
-    (
-        scalar(1)
-      + min( min(Cxpe1_, Cxpe2_*magU/up)*sqrt(bl), Cxpe3_)
-      * filterMult
-    );
-
-    return XiSubEq;
-}
-
-
-bool Foam::XiEqModels::normBasicSubGrid::read(const dictionary& XiEqProperties)
-{
-    XiEqModel::read(XiEqProperties);
-
-    XiEqModelCoeffs_.readEntry("Cxpe1", Cxpe1_);
-    XiEqModelCoeffs_.readEntry("Cxpe2", Cxpe2_);
-    XiEqModelCoeffs_.readEntry("Cxpe3", Cxpe3_);
-    XiEqModelCoeffs_.readEntry("Cxpe4", Cxpe4_);
-
-    return true;
-}
-
-
-// ************************************************************************* //

applications/solvers/combustion/PDRFoam/PDRModels/XiGModels/basicXiSubG/basicXiSubG.C

@@ -0,0 +1,115 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2016 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+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 3 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, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "basicXiSubG.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace XiGModels
+{
+    defineTypeNameAndDebug(basicSubGrid, 0);
+    addToRunTimeSelectionTable(XiGModel, basicSubGrid, dictionary);
+};
+};
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::XiGModels::basicSubGrid::basicSubGrid
+(
+    const dictionary& XiGProperties,
+    const psiuReactionThermo& thermo,
+    const compressible::RASModel& turbulence,
+    const volScalarField& Su
+)
+:
+    XiGModel(XiGProperties, thermo, turbulence, Su),
+    k1(XiGModelCoeffs_.get<scalar>("k1")),
+    XiGModel_(XiGModel::New(XiGModelCoeffs_, thermo, turbulence, Su))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::XiGModels::basicSubGrid::~basicSubGrid()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField> Foam::XiGModels::basicSubGrid::G() const
+{
+    const objectRegistry& db = Su_.db();
+    const volVectorField& U = db.lookupObject<volVectorField>("U");
+    const volScalarField& Nv = db.lookupObject<volScalarField>("Nv");
+    const volScalarField& Lobs = db.lookupObject<volScalarField>("Lobs");
+
+    tmp<volScalarField> tGtot = XiGModel_->G();
+    volScalarField& Gtot = tGtot.ref();
+
+    const scalarField Cw = pow(Su_.mesh().V(), 2.0/3.0);
+    scalarField N(Nv.primitiveField()*Cw);
+
+    forAll(N, celli)
+    {
+        if (N[celli] > 1e-3)
+        {
+            Gtot[celli] += k1*mag(U[celli])/Lobs[celli];
+        }
+    }
+
+    return tGtot;
+}
+
+
+Foam::tmp<Foam::volScalarField> Foam::XiGModels::basicSubGrid::Db() const
+{
+    const objectRegistry& db = Su_.db();
+    const volScalarField& Xi = db.lookupObject<volScalarField>("Xi");
+    const volScalarField& rho = db.lookupObject<volScalarField>("rho");
+    const volScalarField& mgb = db.lookupObject<volScalarField>("mgb");
+    const volScalarField& Lobs = db.lookupObject<volScalarField>("Lobs");
+
+    return XiGModel_->Db()
+        + rho*Su_*(Xi - 1.0)*mgb*(0.5*Lobs)*Lobs/(mgb*Lobs + 1.0);
+}
+
+
+bool Foam::XiGModels::basicSubGrid::read(const dictionary& XiGProperties)
+{
+    XiGModel::read(XiGProperties);
+
+    XiGModelCoeffs_.readEntry("k1", k1);
+
+    return true;
+}
+
+
+// ************************************************************************* //

applications/solvers/combustion/PDRFoam/PDRModels/XiGModels/basicXiSubG/basicXiSubG.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2021 OpenCFD Ltd.
+    Copyright (C) 2011-2013 OpenFOAM Foundation
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -24,18 +24,37 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::XiGModel::normBasicSubGrid
+    Foam::XiGModels::basicSubGrid
 
 Description
 
+    Basic sub-grid obstacle flame-wrinkling generation rate coefficient model.
+    Details supplied by J Puttock 2/7/06.
+
+    \f$ G_{sub} \f$ denotes the generation coefficient and it is given by
+
+    \f[
+        G_{sub} = k_{1} /frac{\vert \dwea{\vec{U}} \vert}{L_{obs}}
+                 \frac{/Xi_{{sub}_{eq}}-1}{/Xi_{sub}}
+    \f]
+
+    and the removal:
+
+    \f[
+        - k_{1} /frac{\vert \dwea{\vec{U}} \vert}{L_{sub}}
+        \frac{\Xi_{sub}-1}{\Xi_{sub}}
+    \f]
+
+    Finally, \f$ G_{sub} \f$ is added to generation rate \f$ G_{in} \f$
+    due to the turbulence.
 
 SourceFiles
-    normBasicSubGrid.C
+    basicSubGrid.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef normBasicSubGrid_H
-#define normBasicSubGrid_H
+#ifndef basicSubGrid_H
+#define basicSubGrid_H
 
 #include "XiGModel.H"
 
@@ -47,56 +66,43 @@ namespace XiGModels
 {
 
 /*---------------------------------------------------------------------------*\
-                          Class normBasicSubGrid Declaration
+                        Class basicSubGrid Declaration
 \*---------------------------------------------------------------------------*/
 
-class normBasicSubGrid
+class basicSubGrid
 :
     public XiGModel
 {
-    // Private Data
+    // Private data
 
         //- Sub-grid generation rate coefficient
-        scalar k1_;
+        scalar k1;
 
-        //- Sub-grid generation rate coefficient - * sqrt(b)
-        scalar kb1_;
-
-        //- Sub-grid generation rate coefficient - * b
-        scalar kbe_;
-
-        //- Sub-grid generation rate upper limit coefficient - * b
-        scalar kbx_;
-
-        //- Sub-grid generation rate coefficient
-        scalar k2_;
-
-        //- Maximum Lobs/CellWidth
-        scalar LOverCw_;
+        //- Xi generation rate model due to turbulence
+        autoPtr<XiGModel> XiGModel_;
 
 
     // Private Member Functions
 
         //- No copy construct
-        normBasicSubGrid(const normBasicSubGrid&) = delete;
+        basicSubGrid(const basicSubGrid&) = delete;
 
         //- No copy assignment
-        void operator=(const normBasicSubGrid&) = delete;
+        void operator=(const basicSubGrid&) = delete;
 
 
 public:
 
     //- Runtime type information
-    TypeName("normBasicSubGridG");
+    TypeName("basicSubGridG");
 
 
     // Constructors
 
         //- Construct from components
-        normBasicSubGrid
+        basicSubGrid
         (
             const dictionary& XiGProperties,
-            const word& modelType,
             const psiuReactionThermo& thermo,
             const compressible::RASModel& turbulence,
             const volScalarField& Su
@@ -104,7 +110,7 @@ public:
 
 
     //- Destructor
-    virtual ~normBasicSubGrid() = default;
+    virtual ~basicSubGrid();
 
 
     // Member Functions

applications/solvers/combustion/PDRFoam/PDRModels/XiGModels/normBasicXiSubG/normBasicXiSubG.C

@@ -1,270 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2021 OpenCFD Ltd.
--------------------------------------------------------------------------------
-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 3 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, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "normBasicXiSubG.H"
-#include "zeroGradientFvPatchField.H"
-#include "addToRunTimeSelectionTable.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace XiGModels
-{
-    defineTypeNameAndDebug(normBasicSubGrid, 0);
-    addToRunTimeSelectionTable(XiGModel, normBasicSubGrid, dictionary);
-}
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::XiGModels::normBasicSubGrid::normBasicSubGrid
-(
-    const dictionary& XiGProperties,
-    const word& modelType,
-    const psiuReactionThermo& thermo,
-    const compressible::RASModel& turbulence,
-    const volScalarField& Su
-)
-:
-    XiGModel(XiGProperties, modelType, thermo, turbulence, Su),
-//     Bv_
-//     (
-//         IOobject
-//         (
-//             "Bv",
-//             Su.mesh().facesInstance(),
-//             Su.mesh(),
-//             IOobject::MUST_READ,
-//             IOobject::NO_WRITE
-//         ),
-//         Su.mesh()
-//     ),
-    k1_(XiGModelCoeffs_.get<scalar>("k1")),
-    kb1_(XiGModelCoeffs_.get<scalar>("kb1")),
-    kbe_(XiGModelCoeffs_.get<scalar>("kbe")),
-    kbx_(XiGModelCoeffs_.get<scalar>("kbx")),
-    k2_(XiGModelCoeffs_.get<scalar>("k2")),
-    LOverCw_(XiGModelCoeffs_.get<scalar>("LOverCw"))
-{}
-
-
-// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
-
-Foam::tmp<Foam::volScalarField> Foam::XiGModels::normBasicSubGrid::G() const
-{
-    const objectRegistry& db = Su_.db();
-    const fvMesh& mesh = Su_.mesh();
-    const volVectorField& U = db.lookupObject<volVectorField>("U");
-    const volScalarField& b = db.lookupObject<volScalarField>("b");
-    const volScalarField& Nv = db.lookupObject<volScalarField>("Nv");
-    const volScalarField& St = db.lookupObject<volScalarField>("St");
-    const volSymmTensorField& nsv = db.lookupObject<volSymmTensorField>("nsv");
-    const volScalarField& Lobs = db.lookupObject<volScalarField>("Lobs");
-    const volSymmTensorField& Bv = db.lookupObject<volSymmTensorField>("Bv");
-
-    const scalarField Cw(pow(Su_.mesh().V(), 2.0/3.0));
-    volScalarField CwVol
-    (
-        IOobject
-        (
-            "CwVol",
-            mesh.time().timeName(),
-            mesh
-        ),
-        mesh,
-        dimensionSet(dimLength),
-        Cw,
-        zeroGradientFvPatchField<scalar>::typeName
-    );
-    CwVol.correctBoundaryConditions();
-
-    if (!db.foundObject<volScalarField>("Ep"))
-    {
-        FatalErrorInFunction
-            << "Looking for Ep in db that does not exist" << nl
-            << Foam::abort(FatalError);
-    }
-
-    const volScalarField& Ep = db.lookupObject<volScalarField>("Ep");
-    const volScalarField& Xp = db.lookupObject<volScalarField>("Xp");
-    const volScalarField& Xi = db.lookupObject<volScalarField>("Xi");
-
-    //tmp<volScalarField> tGtot = XiGModel_->G();
-    auto tGtot = tmp<volScalarField>::New
-    (
-        IOobject
-        (
-            "tGtot",
-            Su_.mesh().time().timeName(),
-            Su_.mesh(),
-            IOobject::NO_READ,
-            IOobject::NO_WRITE,
-            false
-        ),
-        Su_.mesh(),
-        dimensionedScalar(inv(dimTime), Zero)
-    );
-    auto& Gtot = tGtot.ref();
-
-    //calculating flame normal
-
-    volVectorField flNormal(fvc::reconstruct(fvc::snGrad(b)*mesh.magSf()));
-    volScalarField mgb("mgb", mag(flNormal));
-    dimensionedScalar dMgb("dMgb", mgb.dimensions(), SMALL);
-    {
-        volScalarField bc(b*(1.0-b));
-
-        dMgb += 1.0e-8*
-            (bc*mgb)().weightedAverage(mesh.V())
-           /(bc.weightedAverage(mesh.V()) + SMALL);
-    }
-    mgb += dMgb;
-    flNormal /= mgb;
-
-
-    auto tN = tmp<volScalarField>::New
-    (
-        IOobject
-        (
-            "tN",
-            Su_.mesh().time().timeName(),
-            Su_.mesh(),
-            IOobject::NO_READ,
-            IOobject::NO_WRITE,
-            false
-        ),
-        Su_.mesh(),
-        dimensionedScalar(Nv.dimensions(), Zero)
-    );
-    auto& N = tN.ref();
-
-    N.primitiveFieldRef() = Nv.primitiveField()*Cw;
-
-    auto tns = tmp<volSymmTensorField>::New
-    (
-        IOobject
-        (
-            "tns",
-            Su_.mesh().time().timeName(),
-            Su_.mesh(),
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        Su_.mesh(),
-        dimensionedSymmTensor(nsv.dimensions(), Zero)
-    );
-    auto& ns = tns.ref();
-
-    ns.primitiveFieldRef() = nsv.primitiveField()*Cw;
-
-    auto tB = tmp<volSymmTensorField>::New
-    (
-        IOobject
-        (
-            "tB",
-            Su_.mesh().time().timeName(),
-            Su_.mesh(),
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        Su_.mesh(),
-        dimensionedSymmTensor(Bv.dimensions(), Zero)
-    );
-    auto& B = tB.ref();
-
-    B.primitiveFieldRef() = Bv.primitiveField()*sqrt(Cw);
-
-    volScalarField Np(max(N - (flNormal & ns & flNormal), scalar(1)));
-
-    // B_ is Bv*sqrt(Cw)
-    volScalarField bl("bl",(flNormal & B & flNormal)/sqrt(Np));
-    bl.min(1.0);
-
-    volScalarField flSpeed("flSpeed", ((U & flNormal) + St)*b/(b+SMALL)) ;
-
-    volScalarField up("up", sqrt((2.0/3.0)*turbulence_.k()));
-
-    const volScalarField Gtot1
-    (
-        "Gtot1",
-        (
-            k1_ + kb1_*min(pow(bl, kbe_), kbx_)
-        )*mag(flSpeed)/(max(Lobs, LOverCw_*CwVol))
-    );
-
-    const volScalarField Gtot2("Gtot2", k2_*Ep*Su_*Xi/(Xp - 0.999));
-
-    const volScalarField value(pos(N - 1.e-3)*neg(bl - 0.99));
-
-    Gtot = value*Gtot1+(1.0 - value)*Gtot2;
-
-
-    /// if (Xi.mesh().time().outputTime())
-    /// {
-    ///     Gtot.write();
-    ///     bl.write();
-    ///     Lobs.write();
-    ///     flSpeed.write();
-    ///     N.write();
-    /// }
-
-    return tGtot;
-}
-
-
-Foam::tmp<Foam::volScalarField> Foam::XiGModels::normBasicSubGrid::Db() const
-{
-    // Not used //
-    const objectRegistry& db = Su_.db();
-    const volScalarField& Xi = db.lookupObject<volScalarField>("Xi");
-    const volScalarField& rho = db.lookupObject<volScalarField>("rho");
-    const volScalarField& mgb = db.lookupObject<volScalarField>("mgb");
-    const volScalarField& Lobs = db.lookupObject<volScalarField>("Lobs");
-    const volScalarField& Db = db.lookupObject<volScalarField>("Db");
-
-    //return  turbulence_.muEff()
-    return Db + rho*Su_*(Xi - 1.0)*mgb*(0.5*Lobs)*Lobs/(mgb*Lobs + 1.0);
-}
-
-
-bool Foam::XiGModels::normBasicSubGrid::read(const dictionary& XiGProperties)
-{
-    XiGModel::read(XiGProperties);
-
-    XiGModelCoeffs_.readEntry("k1", k1_);
-    XiGModelCoeffs_.readEntry("kb1", kb1_);
-    XiGModelCoeffs_.readEntry("kbe", kbe_);
-    XiGModelCoeffs_.readEntry("kbx", kbx_);
-    XiGModelCoeffs_.readEntry("k2", k2_);
-
-    return true;
-}
-
-
-// ************************************************************************* //

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

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2011-2015 OpenFOAM Foundation
+    Copyright (C) 2011-2017 OpenFOAM Foundation
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -81,7 +81,7 @@ Foam::PDRDragModel::~PDRDragModel()
 
 bool Foam::PDRDragModel::read(const dictionary& PDRProperties)
 {
-    PDRDragModelCoeffs_ = PDRProperties.subDict(type() + "Coeffs");
+    PDRDragModelCoeffs_ = PDRProperties.optionalSubDict(type() + "Coeffs");
 
     PDRDragModelCoeffs_.readEntry("drag", on_);
 
@@ -89,10 +89,4 @@ bool Foam::PDRDragModel::read(const dictionary& PDRProperties)
 }
 
 
-void Foam::PDRDragModel::writeFields() const
-{
-    NotImplemented;
-}
-
-
 // ************************************************************************* //

applications/solvers/combustion/PDRFoam/PDRModels/dragModels/PDRDragModel/PDRDragModel.H

@@ -57,9 +57,10 @@ class PDRDragModel
 :
     public regIOobject
 {
+
 protected:
 
-    // Protected Data
+    // Protected data
 
         dictionary PDRDragModelCoeffs_;
 
@@ -71,7 +72,9 @@ protected:
         bool on_;
 
 
-    // Protected Member Functions
+private:
+
+    // Private Member Functions
 
         //- No copy construct
         PDRDragModel(const PDRDragModel&) = delete;
@@ -143,7 +146,7 @@ public:
     // Member Functions
 
         //- Return true if the drag model is switched on
-        bool on() const noexcept
+        bool on() const
         {
             return on_;
         }
@@ -158,16 +161,17 @@ public:
         using regIOobject::read;
 
         //- Update properties from given dictionary
-        virtual bool read(const dictionary& PDRProperties);
+        virtual bool read(const dictionary& PDRProperties) = 0;
 
-        //- Write Data
         virtual bool writeData(Ostream&) const
         {
             return true;
         }
 
-        //- Write fields
-        virtual void writeFields() const;
+        virtual void writeFields() const
+        {
+            NotImplemented;
+        }
 };
 
 

applications/solvers/combustion/PDRFoam/PDRModels/dragModels/PDRDragModel/PDRDragModelNew.C

@@ -56,7 +56,8 @@ Foam::autoPtr<Foam::PDRDragModel> Foam::PDRDragModel::New
         ) << exit(FatalIOError);
     }
 
-    return autoPtr<PDRDragModel>(cstrIter()(dict, turbulence, rho, U, phi));
+    return autoPtr<PDRDragModel>
+        (cstrIter()(dict, turbulence, rho, U, phi));
 }
 
 

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

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2011 OpenFOAM Foundation
+    Copyright (C) 2011-2016 OpenFOAM Foundation
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -93,20 +93,24 @@ Foam::PDRDragModels::basic::~basic()
 
 Foam::tmp<Foam::volSymmTensorField> Foam::PDRDragModels::basic::Dcu() const
 {
-    auto tDragDcu = tmp<volSymmTensorField>::New
+    tmp<volSymmTensorField> tDragDcu
     (
-        IOobject
+        new volSymmTensorField
         (
-            "tDragDcu",
-            U_.mesh().time().constant(),
+            IOobject
+            (
+                "tDragDcu",
+                U_.mesh().time().constant(),
+                U_.mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
             U_.mesh(),
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        U_.mesh(),
-        dimensionedSymmTensor(dimMass/dimTime/dimVolume, Zero)
+            dimensionedSymmTensor(dimMass/dimTime/dimVolume, Zero)
+        )
     );
-    auto& DragDcu = tDragDcu.ref();
+
+    volSymmTensorField& DragDcu = tDragDcu.ref();
 
     if (on_)
     {
@@ -123,20 +127,24 @@ Foam::tmp<Foam::volSymmTensorField> Foam::PDRDragModels::basic::Dcu() const
 
 Foam::tmp<Foam::volScalarField> Foam::PDRDragModels::basic::Gk() const
 {
-    auto tGk = tmp<volScalarField>::New
+    tmp<volScalarField> tGk
     (
-        IOobject
+        new volScalarField
         (
-            "tGk",
-            U_.mesh().time().constant(),
+            IOobject
+            (
+                "tGk",
+                U_.mesh().time().constant(),
+                U_.mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
             U_.mesh(),
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        U_.mesh(),
-        dimensionedScalar(dimMass/dimLength/pow3(dimTime), Zero)
+            dimensionedScalar(dimMass/dimLength/pow3(dimTime), Zero)
+        )
     );
-    auto& Gk = tGk.ref();
+
+    volScalarField& Gk = tGk.ref();
 
     if (on_)
     {
@@ -151,7 +159,6 @@ Foam::tmp<Foam::volScalarField> Foam::PDRDragModels::basic::Gk() const
           + Csk*betav*turbulence_.muEff()*sqr(Aw_)*magSqr(U_);
     }
 
-
     return tGk;
 }
 
@@ -173,5 +180,4 @@ void Foam::PDRDragModels::basic::writeFields() const
     CR_.write();
 }
 
-
 // ************************************************************************* //

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

@@ -99,7 +99,7 @@ class basic
 :
     public PDRDragModel
 {
-    // Private Data
+    // Private data
 
         dimensionedScalar Csu;
         dimensionedScalar Csk;

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

@@ -1,296 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2021 OpenCFD Ltd.
--------------------------------------------------------------------------------
-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 3 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, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "basicSch.H"
-#include "addToRunTimeSelectionTable.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace PDRDragModels
-{
-    defineTypeNameAndDebug(basicSch, 0);
-    addToRunTimeSelectionTable(PDRDragModel, basicSch, dictionary);
-}
-}
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::PDRDragModels::basicSch::basicSch
-(
-    const dictionary& PDRProperties,
-    const compressible::RASModel& turbulence,
-    const volScalarField& rho,
-    const volVectorField& U,
-    const surfaceScalarField& phi
-)
-:
-    PDRDragModel(PDRProperties,turbulence, rho, U, phi),
-    Csu("Csu", dimless, PDRDragModelCoeffs_),
-    Csk("Csk", dimless, PDRDragModelCoeffs_),
-    Aw_
-    (
-        IOobject
-        (
-            "Aw",
-            U_.mesh().facesInstance(),
-            U_.mesh(),
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE
-        ),
-        U_.mesh()
-    ),
-
-    CR_
-    (
-        IOobject
-        (
-            "CR",
-            U_.mesh().facesInstance(),
-            U_.mesh(),
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE
-        ),
-        U_.mesh()
-    ),
-    nrCoef_(PDRDragModelCoeffs_.get<scalar>("nrCoef")),
-    nrExp2_(PDRDragModelCoeffs_.get<scalar>("nrExp2")),
-    lCoef_(PDRDragModelCoeffs_.get<scalar>("lCoef")),
-    maxSchFac_(PDRDragModelCoeffs_.get<scalar>("maxSchFac")),
-    subGridSchelkin_(PDRDragModelCoeffs_.get<bool>("subGridSchelkin"))
-{}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::PDRDragModels::basicSch::~basicSch()
-{}
-
-
-// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
-
-Foam::tmp<Foam::volSymmTensorField> Foam::PDRDragModels::basicSch::Dcu() const
-{
-    auto tDragDcu = tmp<volSymmTensorField>::New
-    (
-        IOobject
-        (
-            "tDragDcu",
-            U_.mesh().time().constant(),
-            U_.mesh(),
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        U_.mesh(),
-        dimensionedSymmTensor(dimMass/dimTime/dimVolume, Zero)
-    );
-    auto& DragDcu = tDragDcu.ref();
-
-    if (on_)
-    {
-        const volScalarField& betav =
-            U_.db().lookupObject<volScalarField>("betav");
-
-        DragDcu =
-            (0.5*rho_)*CR_*mag(U_) + (Csu*I)*betav*turbulence_.muEff()*sqr(Aw_);
-    }
-
-    return tDragDcu;
-}
-
-
-Foam::tmp<Foam::volScalarField> Foam::PDRDragModels::basicSch::Gk() const
-{
-    auto tGk = tmp<volScalarField>::New
-    (
-        IOobject
-        (
-            "tGk",
-            U_.mesh().time().constant(),
-            U_.mesh(),
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        U_.mesh(),
-        dimensionedScalar(dimMass/dimLength/pow3(dimTime), Zero)
-    );
-    auto& Gk = tGk.ref();
-
-    if (on_)
-    {
-        const volScalarField& betav =
-            U_.db().lookupObject<volScalarField>("betav");
-
-        const volSymmTensorField& CT =
-            U_.db().lookupObject<volSymmTensorField>("CT");
-
-        Gk =
-            (0.5*rho_)*mag(U_)*(U_ & CT & U_)
-          + Csk*betav*turbulence_.muEff()*sqr(Aw_)*magSqr(U_);
-
-        if (subGridSchelkin_)
-        {
-            Gk *= this->SchelkinFactor();
-        }
-    }
-
-    return tGk;
-}
-
-
-Foam::tmp<Foam::volScalarField>
-Foam::PDRDragModels::basicSch::SchelkinFactor() const
-{
-    const volScalarField& Su_ = U_.db().lookupObject<volScalarField>("Su");
-    const volScalarField& rhou_ = U_.db().lookupObject<volScalarField>("rhou");
-    const volScalarField& muu_ = U_.db().lookupObject<volScalarField>("muu");
-
-    auto tfac = tmp<volScalarField>::New
-    (
-        IOobject
-        (
-            "tfac",
-            U_.mesh().time().constant(),
-            U_.mesh(),
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        U_.mesh(),
-        dimensionedScalar(dimless, Zero)
-    );
-    auto& schFac = tfac.ref();
-
-    const volScalarField& k = turbulence_.k();
-    const volScalarField& epsilon = turbulence_.epsilon();
-
-    const volScalarField up(sqrt((2.0/3.0)*k));
-    const volScalarField l(lCoef_*sqrt(3.0/2.0)*up*k/epsilon);
-
-    volScalarField Rs(Su_*l*rhou_/muu_);
-
-    if (subGridSchelkin_)
-    {
-        schFac = max
-        (
-            1.0,
-            min
-            (
-                maxSchFac_,
-                pow(Rs, 2.0 * SchelkinExponent(nrCoef_, nrExp2_, Su_))
-            )
-        );
-    }
-
-    return tfac;
-}
-
-
-Foam::tmp<Foam::volScalarField> Foam::PDRDragModels::basicSch::SchelkinExponent
-(
-    const scalar nrCoef,
-    const scalar nrExp,
-    const volScalarField& Su
-) const
-{
-    const fvMesh& mesh = Su.mesh();
-
-    const volVectorField& U = mesh.lookupObject<volVectorField>("U");
-
-    const volScalarField& Nv = mesh.lookupObject<volScalarField>("Nv");
-    const volSymmTensorField& nsv =
-        mesh.lookupObject<volSymmTensorField>("nsv");
-
-    auto tN = tmp<volScalarField>::New
-    (
-        IOobject
-        (
-            "tN",
-            mesh.time().timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE,
-            false
-        ),
-        mesh,
-        dimensionedScalar(Nv.dimensions(), Zero)
-    );
-    auto& N = tN.ref();
-
-    N.primitiveFieldRef() = Nv.primitiveField()*pow(mesh.V(), 2.0/3.0);
-
-    auto tns = tmp<volSymmTensorField>::New
-    (
-        IOobject
-        (
-            "tns",
-            mesh.time().timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh,
-        dimensionedSymmTensor(nsv.dimensions(), Zero)
-    );
-    auto& ns = tns.ref();
-
-    ns.primitiveFieldRef() = nsv.primitiveField()*pow(mesh.V(), 2.0/3.0);
-
-    const volVectorField Uhat
-    (
-        U/(mag(U) + dimensionedScalar("Usmall", U.dimensions(), 1e-4))
-    );
-
-    const volScalarField nr(sqrt(max(N - (Uhat & ns & Uhat), scalar(1.0))));
-
-    //Re use tN
-    N.primitiveFieldRef() =
-        nrCoef*((scalar(1.0) - pow(nrExp, nr))/(1.0 - nrExp) - scalar(1.0));
-
-    return tN;
-
-}
-
-
-bool Foam::PDRDragModels::basicSch::read(const dictionary& PDRProperties)
-{
-    PDRDragModel::read(PDRProperties);
-
-    PDRDragModelCoeffs_.readEntry("Csu", Csu.value());
-    PDRDragModelCoeffs_.readEntry("Csk", Csk.value());
-
-    return true;
-}
-
-
-void Foam::PDRDragModels::basicSch::writeFields() const
-{
-    Aw_.write();
-    CR_.write();
-}
-
-// ************************************************************************* //

applications/solvers/combustion/PDRFoam/PDRModels/dragModels/basicSch/basicSch.H

@@ -1,141 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2021 OpenCFD Ltd.
--------------------------------------------------------------------------------
-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 3 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, see <http://www.gnu.org/licenses/>.
-
-Class
-    Foam::PDRDragModels::basicSch
-
-Description
-
-SourceFiles
-    basicSch.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef basicSch_H
-#define basicSch_H
-
-#include "PDRDragModel.H"
-#include "XiEqModel.H"
-
-
-namespace Foam
-{
-namespace PDRDragModels
-{
-
-/*---------------------------------------------------------------------------*\
-                          Class basicSch Declaration
-\*---------------------------------------------------------------------------*/
-
-class basicSch
-:
-    public PDRDragModel
-{
-    // Private Data
-
-        dimensionedScalar Csu;
-        dimensionedScalar Csk;
-        volScalarField Aw_;
-        volSymmTensorField CR_;
-
-        //- Schelkin effect Model constants
-        const scalar nrCoef_;
-        const scalar nrExp2_;
-        const scalar lCoef_;
-        const scalar maxSchFac_;
-
-        //- Use sub-grid Schelkin effect
-        bool subGridSchelkin_;
-
-
-    // Private Member Functions
-
-        //- No copy construct
-        basicSch(const basicSch&) = delete;
-
-        //- No copy assignment
-        void operator=(const basicSch&) = delete;
-
-
-public:
-
-    //- Runtime type information
-    TypeName("basicSch");
-
-
-    // Constructors
-
-        //- Construct from components
-        basicSch
-        (
-            const dictionary& PDRProperties,
-            const compressible::RASModel& turbulence,
-            const volScalarField& rho,
-            const volVectorField& U,
-            const surfaceScalarField& phi
-        );
-
-
-    //- Destructor
-    virtual ~basicSch();
-
-
-    // Member Functions
-
-        //- Return the momentum drag coefficient
-        virtual tmp<volSymmTensorField> Dcu() const;
-
-        //- Return the momentum drag turbulence generation rate
-        virtual tmp<volScalarField> Gk() const;
-
-        //- Update properties from given dictionary
-        virtual bool read(const dictionary& PDRProperties);
-
-        //- Write fields
-        void writeFields() const;
-
-        //- Return the Schelkin factor for drag turbulence generation rate
-        tmp<volScalarField> SchelkinFactor() const;
-
-        //- Return the sub-grid Schelkin effect exponent
-        tmp<volScalarField> SchelkinExponent
-        (
-            const scalar,
-            const scalar,
-            const volScalarField&
-        ) const;
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace PDRDragModels
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

applications/solvers/combustion/PDRFoam/PDRModels/turbulence/PDRkEpsilon/PDRkEpsilon.C

@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2019-2021 OpenCFD Ltd.
+    Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -29,7 +29,6 @@ License
 #include "PDRkEpsilon.H"
 #include "PDRDragModel.H"
 #include "addToRunTimeSelectionTable.H"
-#include "Time.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -71,40 +70,16 @@ PDRkEpsilon::PDRkEpsilon
         modelName
     ),
 
-    C5_(coeffDict_.get<scalar>("C5")),
-    C6_(coeffDict_.get<scalar>("C6")),
-    maxLOverCellW_(coeffDict_.get<scalar>("maxLOverCellW")),
-    lCoef_(coeffDict_.get<scalar>("lCoef")),
-    noTurbUntil_
+    C4_
     (
-        dimensioned<scalar>::lookupOrAddToDict
+        dimensioned<scalar>::getOrAddToDict
         (
-            "noTurbUntil",
+            "C4",
             coeffDict_,
-            0.0
+            0.1
         )
-    ),
-    LOverLobs_
-    (
-        dimensioned<scalar>::lookupOrAddToDict
-        (
-            "LOverLobs",
-            coeffDict_,
-            0.2
-        )
-    ),
-    LOverMobs_
-    (
-        dimensioned<scalar>::lookupOrAddToDict
-        (
-            "LOverMobs",
-            coeffDict_,
-            0.0
-        )
-    ),
-    bMin_(coeffDict_.get<scalar>("bMin")),
-    flameFilteredG_(coeffDict_.get<bool>("flameFilteredG"))
- {}
+    )
+{}
 
 
 // * * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * //
@@ -119,9 +94,7 @@ bool PDRkEpsilon::read()
 {
     if (RASModel::read())
     {
-        noTurbUntil_.readIfPresent(coeffDict_);
-        LOverLobs_.readIfPresent(coeffDict_);
-        LOverMobs_.readIfPresent(coeffDict_);
+        C4_.readIfPresent(coeffDict_);
         return true;
     }
 
@@ -137,12 +110,16 @@ void PDRkEpsilon::correct()
         nut_ = Cmu_*sqr(k_)/epsilon_;
         nut_.correctBoundaryConditions();
 
+        // Re-calculate thermal diffusivity
+        //***HGWalphat_ = mut_/Prt_;
+        //alphat_.correctBoundaryConditions();
+
         return;
     }
 
     RASModel::correct();
 
-    volScalarField divU("divU",fvc::div(phi_/fvc::interpolate(rho_)));
+    volScalarField divU(fvc::div(phi_/fvc::interpolate(rho_)));
 
     if (mesh_.moving())
     {
@@ -150,9 +127,7 @@ void PDRkEpsilon::correct()
     }
 
     tmp<volTensorField> tgradU = fvc::grad(U_);
-
     volScalarField G(GName(), rho_*nut_*(tgradU() && dev(twoSymm(tgradU()))));
-
     tgradU.clear();
 
     // Update epsilon and G at the wall
@@ -166,102 +141,59 @@ void PDRkEpsilon::correct()
     const volScalarField& Lobs =
         U_.db().lookupObject<volScalarField>("Lobs");
 
-    const volVectorField& U = mesh_.lookupObject<volVectorField>("U");
-
-    const volScalarField& b = mesh_.lookupObject<volScalarField>("b");
-
-    const scalarField Cw(pow(mesh_.V(), 2.0/3.0));
-
     const PDRDragModel& drag =
         U_.db().lookupObject<PDRDragModel>("PDRDragModel");
 
     volScalarField GR(drag.Gk());
 
-    volScalarField LD
+    volScalarField LI
+        (C4_*(Lobs + dimensionedScalar("minLength", dimLength, VSMALL)));
+
+    // Dissipation equation
+    tmp<fvScalarMatrix> epsEqn
     (
-        "LD",
-        (LOverLobs_)*(Lobs + dimensionedScalar("minLength", dimLength, VSMALL))
+        betav*fvm::ddt(rho_, epsilon_)
+      + fvm::div(phi_, epsilon_)
+      - fvm::laplacian(rho_*DepsilonEff(), epsilon_)
+     ==
+        C1_*betav*G*epsilon_/k_
+      + 1.5*pow(Cmu_, 3.0/4.0)*GR*sqrt(k_)/LI
+      - fvm::SuSp(((2.0/3.0)*C1_)*betav*rho_*divU, epsilon_)
+      - fvm::Sp(C2_*betav*rho_*epsilon_/k_, epsilon_)
     );
 
-    LD.primitiveFieldRef() = min(LD.primitiveField(), maxLOverCellW_*sqrt(Cw));
+    epsEqn.ref().relax();
 
-    const volScalarField LI(lCoef_*pow(k_, 3.0/2.0)/epsilon_);
+    epsEqn.ref().boundaryManipulate(epsilon_.boundaryFieldRef());
 
-    IOdictionary combustionProperties
+    solve(epsEqn);
+    bound(epsilon_, epsilonMin_);
+
+
+    // Turbulent kinetic energy equation
+
+    tmp<fvScalarMatrix> kEqn
     (
-        IOobject
-        (
-            "combustionProperties",
-            mesh_.time().constant(),
-            mesh_,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE
-        )
+        betav*fvm::ddt(rho_, k_)
+      + fvm::div(phi_, k_)
+      - fvm::laplacian(rho_*DkEff(), k_)
+     ==
+        betav*G + GR
+      - fvm::SuSp((2.0/3.0)*betav*rho_*divU, k_)
+      - fvm::Sp(betav*rho_*epsilon_/k_, k_)
     );
 
-    ignition ign(combustionProperties, mesh_.time(), U.mesh());
+    kEqn.ref().relax();
+    solve(kEqn);
+    bound(k_, kMin_);
 
-    dimensionedVector ignLoc("ignLoc", dimLength, ign.sites()[0].location());
-
-    dimensionedScalar filtRad2
-    (
-        "filtRad2",
-        dimLength,
-        6.0*ign.sites()[0].diameter()
-    );
-
-    if (flameFilteredG_)
-    {
-        const volScalarField filDist(mag(mesh_.C() - ignLoc));
-        const volScalarField filterG(pos(filDist - filtRad2));
-        const volScalarField filterdivU(pos(filDist - filtRad2));
-        const volScalarField filterGR(pos(filDist - filtRad2)*pos(b - bMin_));
-
-        G *= filterG;
-        GR *= filterGR;
-        divU *= filterdivU;
-    }
-
-    volScalarField Cl(C5_ + (C6_*((LI - LD)/LI)));
-    Cl.max(0.0);
-
-    if (mesh_.time() > noTurbUntil_)
-    {
-        tmp<fvScalarMatrix> epsEqn
-        (
-            betav*fvm::ddt(rho_, epsilon_)
-          + fvm::div(phi_, epsilon_)
-          - fvm::laplacian(rho_*DepsilonEff(), epsilon_)
-        ==
-            C1_*betav*G*epsilon_/k_
-          + Cl*(epsilon_/k_)*GR
-          - fvm::SuSp(((2.0/3.0)*C1_+C3_)*betav*rho_*divU, epsilon_)
-          - fvm::Sp(C2_*betav*rho_*epsilon_/k_, epsilon_)
-        );
-
-        epsEqn.ref().relax();
-        epsEqn.ref().boundaryManipulate(epsilon_.boundaryFieldRef());
-        solve(epsEqn);
-        bound(epsilon_, epsilonMin_);
-
-        tmp<fvScalarMatrix> kEqn
-        (
-            betav*fvm::ddt(rho_, k_)
-        + fvm::div(phi_, k_)
-        - fvm::laplacian(rho_*DkEff(), k_)
-        ==
-          (betav*G + GR)
-        - fvm::SuSp((2.0/3.0)*betav*rho_*divU, k_)
-        - fvm::Sp(betav*rho_*epsilon_/k_, k_)
-        );
-
-        kEqn.ref().relax();
-        solve(kEqn);
-        bound(k_, kMin_);
-    }
     // Re-calculate viscosity
     nut_ = Cmu_*sqr(k_)/epsilon_;
     nut_.correctBoundaryConditions();
+
+    // Re-calculate thermal diffusivity
+    //***HGWalphat_ = mut_/Prt_;
+    //alphat_.correctBoundaryConditions();
 }
 
 

applications/solvers/combustion/PDRFoam/PDRModels/turbulence/PDRkEpsilon/PDRkEpsilon.H

@@ -6,7 +6,6 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -31,6 +30,30 @@ Description
     Standard k-epsilon turbulence model with additional source terms
     corresponding to PDR basic drag model (\link basic.H \endlink)
 
+    The default model coefficients correspond to the following:
+    \verbatim
+        PDRkEpsilonCoeffs
+        {
+            Cmu         0.09;
+            C1          1.44;
+            C2          1.92;
+            C3          -0.33;  // only for compressible
+            C4          0.1;
+            sigmak      1.0;    // only for compressible
+            sigmaEps    1.3;
+            Prt         1.0;    // only for compressible
+        }
+    \endverbatim
+
+    The turbulence source term \f$ G_{R} \f$ appears in the
+    \f$ \kappa-\epsilon \f$ equation for the generation of turbulence due to
+    interaction with unresolved obstacles.
+
+    In the \f$ \epsilon  \f$ equation \f$ C_{1} G_{R} \f$ is added as a source
+    term.
+
+    In the \f$ \kappa \f$ equation \f$ G_{R} \f$ is added as a source term.
+
 SourceFiles
     PDRkEpsilon.C
 
@@ -40,7 +63,6 @@ SourceFiles
 #define compressiblePDRkEpsilon_H
 
 #include "kEpsilon.H"
-#include "ignition.H"
 #include "turbulentFluidThermoModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -66,18 +88,10 @@ class PDRkEpsilon
         >
     >
 {
-    // Private Data
+    // Private data
 
         // Model coefficients
-        scalar C5_;
-        scalar C6_;
-        scalar maxLOverCellW_;
-        scalar lCoef_;
-        dimensionedScalar noTurbUntil_;
-        dimensionedScalar LOverLobs_;
-        dimensionedScalar LOverMobs_;
-        scalar bMin_;
-        bool flameFilteredG_;
+        dimensionedScalar C4_;
 
 
 public:

applications/solvers/combustion/PDRFoam/StCourantNo.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2011 OpenFOAM Foundation
+    Copyright (C) 2011-2016 OpenFOAM Foundation
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -36,9 +36,10 @@ Description
 
     if (mesh.nInternalFaces())
     {
-        scalarField sumPhi
+       scalarField sumPhi
         (
-            fvc::surfaceSum(mag(phiSt))().primitiveField()/rho.primitiveField()
+            fvc::surfaceSum(mag(phiSt))().primitiveField()
+          / rho.primitiveField()
         );
 
         StCoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue();

applications/solvers/combustion/PDRFoam/UEqn.H

@@ -1,6 +1,9 @@
+    MRF.correctBoundaryVelocity(U);
+
     fvVectorMatrix UEqn
     (
         betav*fvm::ddt(rho, U) + fvm::div(phi, U)
+      + MRF.DDt(rho, U)
       + turbulence->divDevRhoReff(U)
      ==
         betav*rho*g

applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/BLMgMaXiEq/BLMgMaXiEq.C

@@ -1,258 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2021 OpenCFD Ltd.
--------------------------------------------------------------------------------
-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 3 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, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "BLMgMaXiEq.H"
-#include "addToRunTimeSelectionTable.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace XiEqModels
-{
-    defineTypeNameAndDebug(BLMgMaXiEq, 0);
-    addToRunTimeSelectionTable(XiEqModel, BLMgMaXiEq, dictionary);
-}
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::XiEqModels::BLMgMaXiEq::BLMgMaXiEq
-(
-    const dictionary& XiEqProperties,
-    const word& modelType,
-    const psiuReactionThermo& thermo,
-    const compressible::RASModel& turbulence,
-    const volScalarField& Su
-)
-:
-    XiEqModel(XiEqProperties, modelType, thermo, turbulence, Su),
-    kaCoef_(XiEqModelCoeffs_.get<scalar>("kaCoef")),
-    lowK0_(XiEqModelCoeffs_.get<scalar>("lowK0")),
-    lowKg_(XiEqModelCoeffs_.get<scalar>("lowKg")),
-    XiEqCoef_(XiEqModelCoeffs_.get<scalar>("XiEqCoef")),
-    alphaCoefP_(XiEqModelCoeffs_.get<scalar>("alphaCoefP")),
-    betaCoefP_(XiEqModelCoeffs_.get<scalar>("betaCoefP")),
-    alphaCoefN_(XiEqModelCoeffs_.get<scalar>("alphaCoefN")),
-    betaCoefN_(XiEqModelCoeffs_.get<scalar>("betaCoefN")),
-    maLim_(XiEqModelCoeffs_.get<scalar>("maLim")),
-    maLim1_(XiEqModelCoeffs_.get<scalar>("maLim1")),
-    quenchCoef_(XiEqModelCoeffs_.get<scalar>("quenchCoef")),
-    quenchExp_(XiEqModelCoeffs_.get<scalar>("quenchExp")),
-    quenchM_(XiEqModelCoeffs_.get<scalar>("quenchM")),
-    quenchRate1_(XiEqModelCoeffs_.get<scalar>("quenchRate1")),
-    quenchRate2_(XiEqModelCoeffs_.get<scalar>("quenchRate2")),
-    lCoef_(XiEqModelCoeffs_.get<scalar>("lCoef")),
-    SuMin_(0.01*Su.average()),
-    uPrimeCoef_(XiEqModelCoeffs_.get<scalar>("uPrimeCoef")),
-    nrExp_(XiEqModelCoeffs_.get<scalar>("nrExp")),
-    subGridSchelkin_(XiEqModelCoeffs_.get<bool>("subGridSchelkin")),
-    MaModel
-    (
-        IOdictionary
-        (
-            IOobject
-            (
-                "combustionProperties",
-                Su.mesh().time().constant(),
-                Su.mesh(),
-                IOobject::MUST_READ
-            )
-        ),
-        thermo
-    )
-{}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::XiEqModels::BLMgMaXiEq::~BLMgMaXiEq()
-{}
-
-
-// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
-
-Foam::tmp<Foam::volScalarField> Foam::XiEqModels::BLMgMaXiEq::XiEq() const
-{
-    const volScalarField& k = turbulence_.k();
-    const volScalarField& epsilon = turbulence_.epsilon();
-    volScalarField up("up", sqrt((2.0/3.0)*k));
-    if (subGridSchelkin_)
-    {
-        up.primitiveFieldRef() +=
-            calculateSchelkinEffect(uPrimeCoef_, nrExp_);
-    }
-
-    volScalarField l(lCoef_*sqrt(3.0/2.0)*up*k/epsilon);
-    volScalarField Rl(up*l*thermo_.rhou()/thermo_.muu());
-
-    volScalarField upBySu("upBySu", up/(Su_ + SuMin_));
-    volScalarField K("K", kaCoef_*upBySu*upBySu/sqrt(Rl));
-    volScalarField Ma("Ma", MaModel.Ma());
-
-    volScalarField regime("regime", MaModel.Ma()*scalar(0.0));
-
-    auto tXiEq = tmp<volScalarField>::New
-    (
-        IOobject
-        (
-            "XiEq",
-            epsilon.time().timeName(),
-            epsilon.db(),
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        epsilon.mesh(),
-        dimensionedScalar(dimless, Zero)
-    );
-
-    const objectRegistry& db = Su_.db();
-
-    const volScalarField& b = db.lookupObject<volScalarField>("b");
-
-    const volScalarField multiP1(0.0*pos(b - 0.99) + 1.0*neg(b - 0.99));
-    const volScalarField multiP2(1.0*pos(b - 0.01) + 0.0*neg(b - 0.01));
-
-    volScalarField& xieq = tXiEq.ref();
-
-    forAll(xieq, celli)
-    {
-        scalar alpha;
-        scalar beta;
-        scalar gulderMa;
-
-        if (Ma[celli]>= 0)
-        {
-            gulderMa =
-                1.0
-              + (0.1402 - 0.007*Ma[celli])
-              * upBySu[celli]*sqrt(upBySu[celli]/K[celli]);
-
-            regime[celli] = multiP1[celli]*multiP2[celli];
-
-        }
-        else
-        {
-            gulderMa =
-                1.0
-              + (0.005*Ma[celli]*Ma[celli]+0.01*Ma[celli] + 0.125)
-              * upBySu[celli]*sqrt(upBySu[celli]/K[celli]);
-
-             regime[celli] = 2*multiP1[celli]*multiP2[celli];
-        }
-
-
-        if (K[celli] < (lowK0_ + lowKg_*Ma[celli]) )
-        {
-            xieq[celli] = gulderMa;
-        }
-        else
-        {
-            if (Ma[celli] >= 0.0)
-            {
-                alpha = alphaCoefP_*(maLim_ - Ma[celli]);
-                beta = betaCoefP_*(maLim_ - Ma[celli]);
-                regime[celli] = 3*multiP1[celli]*multiP2[celli];
-            }
-            else
-            {
-                alpha = alphaCoefN_*(maLim1_ - Ma[celli]) ;
-                beta = betaCoefN_*(maLim_ + Ma[celli]);
-                regime[celli] = 4*multiP1[celli]*multiP2[celli];
-            }
-            xieq[celli] = XiEqCoef_*alpha*pow(K[celli], beta)*upBySu[celli];
-        }
-
-        if (Ma[celli] > -3.0 && Ma[celli] < 11.0)
-        {
-            scalar K0p8 = quenchCoef_*pow( Ma[celli] - quenchM_, quenchExp_);
-            scalar quenchRate = quenchRate1_ + quenchRate2_*Ma[celli];
-            if (K[celli] > (K0p8 - 0.223/quenchRate))
-            {
-                xieq[celli] *= 0.8*exp(-quenchRate*(K[celli] - K0p8));
-                regime[celli] = 5*multiP1[celli]*multiP2[celli];
-            }
-        }
-    }
-
-    forAll(xieq.boundaryField(), patchi)
-    {
-        scalarField& xieqp = xieq.boundaryFieldRef()[patchi];
-        const scalarField& Kp = K.boundaryField()[patchi];
-        const scalarField& Map = Ma.boundaryField()[patchi];
-        const scalarField& upBySup = upBySu.boundaryField()[patchi];
-
-        forAll(xieqp, facei)
-        {
-            scalar alpha;
-            scalar beta;
-
-            if (Map[facei] > 0.0)
-            {
-                alpha = alphaCoefP_*(maLim_ - Map[facei]);
-                beta = betaCoefP_*(maLim_ - Map[facei]);
-            }
-            else
-            {
-                alpha = alphaCoefN_*(maLim_ - Map[facei]);
-                beta = betaCoefN_*(maLim_ + Map[facei]);
-            }
-            xieqp[facei] =
-                XiEqCoef_*alpha*pow(Kp[facei], beta)*upBySup[facei];
-
-            if (Map[facei] > -3.0 && Map[facei] < 11.0)
-            {
-                scalar K0p8 = quenchCoef_*pow(Map[facei] - quenchM_, quenchExp_);
-                scalar quenchRate = quenchRate1_ + quenchRate2_*Ma[facei];
-
-                if (Kp[facei] > (K0p8 - 0.223/quenchRate))
-                {
-                    xieqp[facei] *= 0.8*exp(-quenchRate*(Kp[facei] - K0p8));
-                }
-            }
-            else
-            {
-                Info<<
-                    "Markstein Number out of range for Quench Formulation" << endl;
-            }
-        }
-    }
-
-    return tXiEq;
-}
-
-
-bool Foam::XiEqModels::BLMgMaXiEq::read(const dictionary& XiEqProperties)
-{
-    XiEqModel::read(XiEqProperties);
-
-    return true;
-}
-
-
-// ************************************************************************* //

applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/Gulder/Gulder.C

@@ -5,8 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2011-2012 OpenFOAM Foundation
-    Copyright (C) 2021 OpenCFD Ltd.
+    Copyright (C) 2011-2016 OpenFOAM Foundation
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -46,21 +45,22 @@ namespace XiEqModels
 Foam::XiEqModels::Gulder::Gulder
 (
     const dictionary& XiEqProperties,
-    const word& modelType,
     const psiuReactionThermo& thermo,
     const compressible::RASModel& turbulence,
     const volScalarField& Su
 )
 :
-    XiEqModel(XiEqProperties, modelType, thermo, turbulence, Su),
+    XiEqModel(XiEqProperties, thermo, turbulence, Su),
     XiEqCoef_(XiEqModelCoeffs_.get<scalar>("XiEqCoef")),
     SuMin_(0.01*Su.average()),
     uPrimeCoef_(XiEqModelCoeffs_.get<scalar>("uPrimeCoef")),
-    nrExp_(XiEqModelCoeffs_.get<scalar>("nrExp")),
-    subGridSchelkin_
-    (
-        XiEqModelCoeffs_.get<bool>("subGridSchelkin")
-    )
+    subGridSchelkin_(XiEqModelCoeffs_.get<bool>("subGridSchelkin"))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::XiEqModels::Gulder::~Gulder()
 {}
 
 
@@ -73,7 +73,7 @@ Foam::tmp<Foam::volScalarField> Foam::XiEqModels::Gulder::XiEq() const
 
     if (subGridSchelkin_)
     {
-        up.primitiveFieldRef() += calculateSchelkinEffect(uPrimeCoef_, nrExp_);
+        up.primitiveFieldRef() += calculateSchelkinEffect(uPrimeCoef_);
     }
 
     volScalarField tauEta(sqrt(mag(thermo_.muu()/(thermo_.rhou()*epsilon))));
@@ -97,7 +97,6 @@ bool Foam::XiEqModels::Gulder::read(const dictionary& XiEqProperties)
 
     XiEqModelCoeffs_.readEntry("XiEqCoef", XiEqCoef_);
     XiEqModelCoeffs_.readEntry("uPrimeCoef", uPrimeCoef_);
-    XiEqModelCoeffs_.readEntry("nrExp", nrExp_);
     XiEqModelCoeffs_.readEntry("subGridSchelkin", subGridSchelkin_);
 
     return true;

applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/Gulder/Gulder.H

@@ -55,7 +55,7 @@ class Gulder
 :
     public XiEqModel
 {
-    // Private Data
+    // Private data
 
         //- Model constant
         scalar XiEqCoef_;
@@ -66,8 +66,6 @@ class Gulder
         //- Schelkin effect Model constant
         scalar uPrimeCoef_;
 
-        scalar nrExp_;
-
         //- Use sub-grid Schelkin effect
         bool subGridSchelkin_;
 
@@ -93,7 +91,6 @@ public:
         Gulder
         (
             const dictionary& XiEqProperties,
-            const word& modelType,
             const psiuReactionThermo& thermo,
             const compressible::RASModel& turbulence,
             const volScalarField& Su
@@ -101,7 +98,7 @@ public:
 
 
     //- Destructor
-    virtual ~Gulder() = default;
+    virtual ~Gulder();
 
 
     // Member Functions
@@ -111,6 +108,7 @@ public:
 
         //- Update properties from given dictionary
         virtual bool read(const dictionary& XiEqProperties);
+
 };
 
 

applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/SCOPEBlendXiEq/SCOPEBlendXiEq.C

@@ -0,0 +1,95 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2015 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+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 3 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, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "SCOPEBlendXiEq.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace XiEqModels
+{
+    defineTypeNameAndDebug(SCOPEBlend, 0);
+    addToRunTimeSelectionTable(XiEqModel, SCOPEBlend, dictionary);
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::XiEqModels::SCOPEBlend::SCOPEBlend
+(
+    const dictionary& XiEqProperties,
+    const psiuReactionThermo& thermo,
+    const compressible::RASModel& turbulence,
+    const volScalarField& Su
+)
+:
+    XiEqModel(XiEqProperties, thermo, turbulence, Su),
+    XiEqModelL_
+    (
+        XiEqModel::New
+        (
+            XiEqModelCoeffs_.subDict("XiEqModelL"),
+            thermo,
+            turbulence,
+            Su
+        )
+    ),
+    XiEqModelH_
+    (
+        XiEqModel::New
+        (
+            XiEqModelCoeffs_.subDict("XiEqModelH"),
+            thermo,
+            turbulence,
+            Su
+        )
+    )
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::XiEqModels::SCOPEBlend::~SCOPEBlend()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField> Foam::XiEqModels::SCOPEBlend::XiEq() const
+{
+    return pow
+    (
+        pow4(1.0/XiEqModelL_->XiEq()) + pow4(1.0/XiEqModelH_->XiEq()),
+        -0.25
+    );
+}
+
+
+// ************************************************************************* //

applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/SCOPEBlendXiEq/SCOPEBlendXiEq.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2021 OpenCFD Ltd.
+    Copyright (C) 2011-2012 OpenFOAM Foundation
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -24,25 +24,22 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::XiEqModels::instability2XiEq
+    Foam::XiEqModels::SCOPEBlend
 
 Description
-
+    Simple SCOPEBlendXiEq model for XiEq based on SCOPEXiEqs correlation
+    with a linear correction function to give a plausible profile for XiEq.
+    See \link SCOPELaminarFlameSpeed.H \endlink for details on the SCOPE
+    laminar flame speed model.
 
 SourceFiles
-    instability2XiEq.C
+    SCOPEBlend.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef instability2XiEq_H
-#define instability2XiEq_H
+#ifndef SCOPEBlend_H
+#define SCOPEBlend_H
 
-#include "laminarFlameSpeed.H"
-#include "SCOPELaminarFlameSpeed.H"
-#include "ignitionSite.H"
-#include "ignition.H"
-#include "Time.H"
-#include "fvMesh.H"
 #include "XiEqModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -53,53 +50,43 @@ namespace XiEqModels
 {
 
 /*---------------------------------------------------------------------------*\
-                          Class instability2 Declaration
+                          Class SCOPEBlend Declaration
 \*---------------------------------------------------------------------------*/
 
-class instability2XiEq
+class SCOPEBlend
 :
     public XiEqModel
 {
     // Private data
 
-        //- Laminar burning speed
-        laminarFlameSpeedModels::SCOPE saModel_;
+        //- Low turbulence intensity equilibrium Xi model
+        autoPtr<XiEqModel> XiEqModelL_;
 
-        //- GIn (initial instability G)divided by Su^2. Read from fuel file
-        scalar CIn_;
-
-        //- Default CIn if not in fuel file
-        scalar defaultCIn_;
-
-        //- Determines how fast XiEqIn fades out as turbulence comes in
-        scalar XiEqInFade_;
-
-        //- Equilibrium Xi model due to all other effects
-        autoPtr<XiEqModel> XiEqModel_;
+        //- High turbulence intensity equilibrium Xi model
+        autoPtr<XiEqModel> XiEqModelH_;
 
 
     // Private Member Functions
 
         //- No copy construct
-        instability2XiEq(const instability2XiEq&) = delete;
+        SCOPEBlend(const SCOPEBlend&) = delete;
 
         //- No copy assignment
-        void operator=(const instability2XiEq&) = delete;
+        void operator=(const SCOPEBlend&) = delete;
 
 
 public:
 
     //- Runtime type information
-    TypeName("instability2XiEq");
+    TypeName("SCOPEBlend");
 
 
     // Constructors
 
         //- Construct from components
-        instability2XiEq
+        SCOPEBlend
         (
             const dictionary& XiEqProperties,
-            const word& modelType,
             const psiuReactionThermo& thermo,
             const compressible::RASModel& turbulence,
             const volScalarField& Su
@@ -107,7 +94,7 @@ public:
 
 
     //- Destructor
-    virtual ~instability2XiEq() = default;
+    virtual ~SCOPEBlend();
 
 
     // Member Functions
@@ -116,7 +103,11 @@ public:
         virtual tmp<volScalarField> XiEq() const;
 
         //- Update properties from given dictionary
-        virtual bool read(const dictionary& XiEqProperties);
+        virtual bool read(const dictionary& XiEqProperties)
+        {
+            return true;
+        }
+
 };
 
 

applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/SCOPEXiEq/SCOPEXiEq.C

@@ -0,0 +1,159 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2016 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+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 3 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, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "SCOPEXiEq.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace XiEqModels
+{
+    defineTypeNameAndDebug(SCOPEXiEq, 0);
+    addToRunTimeSelectionTable(XiEqModel, SCOPEXiEq, dictionary);
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::XiEqModels::SCOPEXiEq::SCOPEXiEq
+(
+    const dictionary& XiEqProperties,
+    const psiuReactionThermo& thermo,
+    const compressible::RASModel& turbulence,
+    const volScalarField& Su
+)
+:
+    XiEqModel(XiEqProperties, thermo, turbulence, Su),
+    XiEqCoef_(XiEqModelCoeffs_.get<scalar>("XiEqCoef")),
+    XiEqExp_(XiEqModelCoeffs_.get<scalar>("XiEqExp")),
+    lCoef_(XiEqModelCoeffs_.get<scalar>("lCoef")),
+    SuMin_(0.01*Su.average()),
+    uPrimeCoef_(XiEqModelCoeffs_.get<scalar>("uPrimeCoef")),
+    subGridSchelkin_(XiEqModelCoeffs_.get<bool>("subGridSchelkin")),
+    MaModel
+    (
+        Su.mesh().lookupObject<IOdictionary>("combustionProperties"),
+        thermo
+    )
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::XiEqModels::SCOPEXiEq::~SCOPEXiEq()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField> Foam::XiEqModels::SCOPEXiEq::XiEq() const
+{
+    const volScalarField& k = turbulence_.k();
+    const volScalarField& epsilon = turbulence_.epsilon();
+
+    volScalarField up(sqrt((2.0/3.0)*k));
+    if (subGridSchelkin_)
+    {
+        up.primitiveFieldRef() += calculateSchelkinEffect(uPrimeCoef_);
+    }
+
+    volScalarField l(lCoef_*sqrt(3.0/2.0)*up*k/epsilon);
+    volScalarField Rl(up*l*thermo_.rhou()/thermo_.muu());
+
+    volScalarField upBySu(up/(Su_ + SuMin_));
+    volScalarField K(0.157*upBySu/sqrt(Rl));
+    volScalarField Ma(MaModel.Ma());
+
+    tmp<volScalarField> tXiEq
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "XiEq",
+                epsilon.time().timeName(),
+                epsilon.db(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            epsilon.mesh(),
+            dimensionedScalar(dimless, Zero)
+        )
+    );
+    volScalarField& xieq = tXiEq.ref();
+
+    forAll(xieq, celli)
+    {
+        if (Ma[celli] > 0.01)
+        {
+            xieq[celli] =
+                XiEqCoef_*pow(K[celli]*Ma[celli], -XiEqExp_)*upBySu[celli];
+        }
+    }
+
+    volScalarField::Boundary& xieqBf = xieq.boundaryFieldRef();
+
+    forAll(xieq.boundaryField(), patchi)
+    {
+        scalarField& xieqp = xieqBf[patchi];
+        const scalarField& Kp = K.boundaryField()[patchi];
+        const scalarField& Map = Ma.boundaryField()[patchi];
+        const scalarField& upBySup = upBySu.boundaryField()[patchi];
+
+        forAll(xieqp, facei)
+        {
+            if (Ma[facei] > 0.01)
+            {
+                xieqp[facei] =
+                    XiEqCoef_*pow(Kp[facei]*Map[facei], -XiEqExp_)
+                   *upBySup[facei];
+            }
+        }
+    }
+
+    return tXiEq;
+}
+
+
+bool Foam::XiEqModels::SCOPEXiEq::read(const dictionary& XiEqProperties)
+{
+    XiEqModel::read(XiEqProperties);
+
+    XiEqModelCoeffs_.readEntry("XiEqCoef", XiEqCoef_);
+    XiEqModelCoeffs_.readEntry("XiEqExp", XiEqExp_);
+    XiEqModelCoeffs_.readEntry("lCoef", lCoef_);
+    XiEqModelCoeffs_.readEntry("uPrimeCoef", uPrimeCoef_);
+    XiEqModelCoeffs_.readEntry("subGridSchelkin", subGridSchelkin_);
+
+    return true;
+}
+
+
+// ************************************************************************* //

applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/SCOPEXiEq/SCOPEXiEq.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2021 OpenCFD Ltd.
+    Copyright (C) 2011-2013 OpenFOAM Foundation
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -24,22 +24,21 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::XiEqModel::BLMgMaXiEq
+    Foam::XiEqModels::SCOPEXiEq
 
 Description
-    Model for XiEq based on Bradley, Lawes and Mansour (2011)
-    Cobustion and Falme, 158, 123 correlation
+    Simple SCOPEXiEq model for XiEq based on SCOPEXiEqs correlation
     with a linear correction function to give a plausible profile for XiEq.
     See \link SCOPELaminarFlameSpeed.H \endlink for details on the SCOPE laminar
     flame speed model.
 
 SourceFiles
-    BLMgMaXiEq.C
+    SCOPEXiEq.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef BLMgMaXiEq_H
-#define BLMgMaXiEq_H
+#ifndef SCOPEXiEq_H
+#define SCOPEXiEq_H
 
 #include "XiEqModel.H"
 #include "SCOPELaminarFlameSpeed.H"
@@ -52,37 +51,29 @@ namespace XiEqModels
 {
 
 /*---------------------------------------------------------------------------*\
-                          Class BLMgMaXiEq Declaration
+                          Class SCOPEXiEq Declaration
 \*---------------------------------------------------------------------------*/
 
-class BLMgMaXiEq
+class SCOPEXiEq
 :
     public XiEqModel
 {
     // Private data
 
-        // Model constants
-
-        scalar kaCoef_;
-        scalar lowK0_;
-        scalar lowKg_;
+        // Model constant
         scalar XiEqCoef_;
-        scalar alphaCoefP_;
-        scalar betaCoefP_;
-        scalar alphaCoefN_;
-        scalar betaCoefN_;
-        scalar maLim_;
-        scalar maLim1_;
-        scalar quenchCoef_, quenchExp_, quenchM_;
-        scalar quenchRate1_, quenchRate2_;
+
+        // Model constant
+        scalar XiEqExp_;
+
+        // Model constant
         scalar lCoef_;
 
         //- Minimum Su
         dimensionedScalar SuMin_;
 
-        //- Schelkin effect Model constants
+        //- Schelkin effect Model constant
         scalar uPrimeCoef_;
-        scalar nrExp_;
 
         //- Use sub-grid Schelkin effect
         bool subGridSchelkin_;
@@ -96,25 +87,24 @@ class BLMgMaXiEq
     // Private Member Functions
 
         //- No copy construct
-        BLMgMaXiEq(const BLMgMaXiEq&) = delete;
+        SCOPEXiEq(const SCOPEXiEq&) = delete;
 
         //- No copy assignment
-        void operator=(const BLMgMaXiEq&) = delete;
+        void operator=(const SCOPEXiEq&) = delete;
 
 
 public:
 
     //- Runtime type information
-    TypeName("BLMgMaXiEq");
+    TypeName("SCOPEXiEq");
 
 
     // Constructors
 
         //- Construct from components
-        BLMgMaXiEq
+        SCOPEXiEq
         (
             const dictionary& XiEqProperties,
-            const word& modelType,
             const psiuReactionThermo& thermo,
             const compressible::RASModel& turbulence,
             const volScalarField& Su
@@ -122,7 +112,7 @@ public:
 
 
     //- Destructor
-    virtual ~BLMgMaXiEq();
+    virtual ~SCOPEXiEq();
 
 
     // Member Functions

applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/XiEqModel/XiEqModel.C

@@ -5,8 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2021 OpenCFD Ltd.
+    Copyright (C) 2011-2017 OpenFOAM Foundation
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -42,7 +41,6 @@ namespace Foam
 Foam::XiEqModel::XiEqModel
 (
     const dictionary& XiEqProperties,
-    const word& modelType,
     const psiuReactionThermo& thermo,
     const compressible::RASModel& turbulence,
     const volScalarField& Su
@@ -52,7 +50,7 @@ Foam::XiEqModel::XiEqModel
     (
         XiEqProperties.subDict
         (
-            XiEqProperties.get<word>(modelType) + "Coeffs"
+            XiEqProperties.get<word>("XiEqModel") + "Coeffs"
         )
     ),
     thermo_(thermo),
@@ -71,51 +69,56 @@ Foam::XiEqModel::~XiEqModel()
 
 bool Foam::XiEqModel::read(const dictionary& XiEqProperties)
 {
-
-    XiEqModelCoeffs_ = XiEqProperties.subDict(type() + "Coeffs");
+    XiEqModelCoeffs_ = XiEqProperties.optionalSubDict(type() + "Coeffs");
 
     return true;
 }
 
 
 void Foam::XiEqModel::writeFields() const
-{}
+{
+    //***HGW It is not clear why B is written here
+    if (Su_.mesh().foundObject<volSymmTensorField>("B"))
+    {
+        const volSymmTensorField& B =
+            Su_.mesh().lookupObject<volSymmTensorField>("B");
+        B.write();
+    }
+}
 
 
-Foam::tmp<Foam::volScalarField>Foam::XiEqModel::calculateSchelkinEffect
-(
-    const scalar uPrimeCoef,
-    const scalar nrExp
-) const
+Foam::tmp<Foam::volScalarField>
+Foam::XiEqModel::calculateSchelkinEffect(const scalar uPrimeCoef) const
 {
     const fvMesh& mesh = Su_.mesh();
 
     const volVectorField& U = mesh.lookupObject<volVectorField>("U");
-
     const volSymmTensorField& CT = mesh.lookupObject<volSymmTensorField>("CT");
     const volScalarField& Nv = mesh.lookupObject<volScalarField>("Nv");
     const volSymmTensorField& nsv =
         mesh.lookupObject<volSymmTensorField>("nsv");
 
-    auto tN = tmp<volScalarField>::New
+    tmp<volScalarField> tN
     (
-        IOobject
+        new volScalarField
         (
-            "tN",
-            mesh.time().timeName(),
+            IOobject
+            (
+                "tN",
+                mesh.time().timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
             mesh,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE,
-            false
-        ),
-        mesh,
-        dimensionedScalar(Nv.dimensions(), Zero)
+            dimensionedScalar(Nv.dimensions(), Zero)
+        )
     );
-    auto& N = tN.ref();
-
+    volScalarField& N = tN.ref();
     N.primitiveFieldRef() = Nv.primitiveField()*pow(mesh.V(), 2.0/3.0);
 
-    auto tns = tmp<volSymmTensorField>::New
+    volSymmTensorField ns
     (
         IOobject
         (
@@ -128,8 +131,6 @@ Foam::tmp<Foam::volScalarField>Foam::XiEqModel::calculateSchelkinEffect
         mesh,
         dimensionedSymmTensor(nsv.dimensions(), Zero)
     );
-    auto& ns = tns.ref();
-
     ns.primitiveFieldRef() = nsv.primitiveField()*pow(mesh.V(), 2.0/3.0);
 
     const volVectorField Uhat
@@ -143,8 +144,10 @@ Foam::tmp<Foam::volScalarField>Foam::XiEqModel::calculateSchelkinEffect
 
     const scalarField upLocal(uPrimeCoef*sqrt((U & CT & U)*cellWidth));
 
-    //Re use tN
-    N.primitiveFieldRef() = upLocal*(max(scalar(1.0), pow(nr, nrExp)) - 1.0);
+    const scalarField deltaUp(upLocal*(max(scalar(1), pow(nr, 0.5)) - 1.0));
+
+    // Re use tN
+    N.primitiveFieldRef() = upLocal*(max(scalar(1), pow(nr, 0.5)) - 1.0);
 
     return tN;
 }

applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/XiEqModel/XiEqModel.H

@@ -6,7 +6,6 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -29,6 +28,12 @@ Class
 
 Description
     Base-class for all XiEq models used by the b-XiEq combustion model.
+    The available models are :
+        \link basicXiSubXiEq.H \endlink
+        \link Gulder.H \endlink
+        \link instabilityXiEq.H \endlink
+        \link SCOPEBlendXiEq.H \endlink
+        \link SCOPEXiEq.H \endlink
 
 SourceFiles
     XiEqModel.C
@@ -54,9 +59,10 @@ namespace Foam
 
 class XiEqModel
 {
+
 protected:
 
-    // Protected Data
+    // Protected data
 
         //- Dictionary
         dictionary XiEqModelCoeffs_;
@@ -71,7 +77,9 @@ protected:
         const volScalarField& Su_;
 
 
-    // Protected Member Functions
+private:
+
+    // Private Member Functions
 
         //- No copy construct
         XiEqModel(const XiEqModel&) = delete;
@@ -95,14 +103,12 @@ public:
             dictionary,
             (
                 const dictionary& XiEqProperties,
-                const word& modelType,
                 const psiuReactionThermo& thermo,
                 const compressible::RASModel& turbulence,
                 const volScalarField& Su
             ),
             (
                 XiEqProperties,
-                modelType,
                 thermo,
                 turbulence,
                 Su
@@ -116,7 +122,6 @@ public:
         static autoPtr<XiEqModel> New
         (
             const dictionary& XiEqProperties,
-            const word& modelType,
             const psiuReactionThermo& thermo,
             const compressible::RASModel& turbulence,
             const volScalarField& Su
@@ -129,7 +134,6 @@ public:
         XiEqModel
         (
             const dictionary& XiEqProperties,
-            const word& modelType,
             const psiuReactionThermo& thermo,
             const compressible::RASModel& turbulence,
             const volScalarField& Su
@@ -143,14 +147,13 @@ public:
     // Member Functions
 
         //- Return the flame-wrinkling XiEq
-        virtual tmp<volScalarField> XiEq() const = 0;
+        virtual tmp<volScalarField> XiEq() const
+        {
+            return turbulence_.muEff();
+        }
 
-        // Only used by sperseded modles
-        tmp<volScalarField> calculateSchelkinEffect
-        (
-            const scalar uPrimeCoef,
-            const scalar nrExp
-        ) const;
+        //- Return the sub-grid Schelkin effect
+        tmp<volScalarField> calculateSchelkinEffect(const scalar) const;
 
         //- Update properties from given dictionary
         virtual bool read(const dictionary& XiEqProperties) = 0;

applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/XiEqModel/XiEqModelNew.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2019-2021 OpenCFD Ltd.
+    Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -33,15 +33,14 @@ License
 Foam::autoPtr<Foam::XiEqModel> Foam::XiEqModel::New
 (
     const dictionary& dict,
-    const word& modelName,
     const psiuReactionThermo& thermo,
     const compressible::RASModel& turbulence,
     const volScalarField& Su
 )
 {
-    const word modelType(dict.get<word>(modelName));
+    const word modelType(dict.get<word>("XiEqModel"));
 
-    Info<< "Selecting flame-wrinkling model Eq " << modelType << endl;
+    Info<< "Selecting flame-wrinkling model " << modelType << endl;
 
     auto cstrIter = dictionaryConstructorTablePtr_->cfind(modelType);
 
@@ -50,16 +49,13 @@ Foam::autoPtr<Foam::XiEqModel> Foam::XiEqModel::New
         FatalIOErrorInLookup
         (
             dict,
-            modelName,
+            "XiEqModel",
             modelType,
             *dictionaryConstructorTablePtr_
         ) << exit(FatalIOError);
     }
 
-    return autoPtr<XiEqModel>(cstrIter()
-    (
-        dict, modelName, thermo, turbulence, Su)
-    );
+    return autoPtr<XiEqModel>(cstrIter()(dict, thermo, turbulence, Su));
 }
 
 

applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/instability2XiEq/instability2XiEq.C

@@ -1,133 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2021 OpenCFD Ltd.
--------------------------------------------------------------------------------
-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 3 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, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "instability2XiEq.H"
-#include "addToRunTimeSelectionTable.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace XiEqModels
-{
-    defineTypeNameAndDebug(instability2XiEq, 0);
-    addToRunTimeSelectionTable(XiEqModel, instability2XiEq, dictionary);
-}
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::XiEqModels::instability2XiEq::instability2XiEq
-(
-    const dictionary& XiEqProperties,
-    const word& modelType,
-    const psiuReactionThermo& thermo,
-    const compressible::RASModel& turbulence,
-    const volScalarField& Su
-)
-:
-    XiEqModel(XiEqProperties, modelType, thermo, turbulence, Su),
-    saModel_
-    (
-        IOdictionary
-        (
-            IOobject
-            (
-                "combustionProperties",
-                Su.mesh().time().constant(),
-                Su.mesh(),
-                IOobject::MUST_READ
-            )
-        ),
-        thermo
-    ),
-    CIn_(saModel_.CIn()),
-    defaultCIn_(XiEqModelCoeffs_.get<scalar>("defaultCIn")),
-    XiEqInFade_(XiEqModelCoeffs_.get<scalar>("XiEqInFade")),
-    XiEqModel_
-    (
-        XiEqModel::New(XiEqModelCoeffs_, modelType, thermo, turbulence, Su)
-    )
-{
-    if (CIn_ <= 0.0)
-    {
-        CIn_ = defaultCIn_;
-    }
-}
-
-
-// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
-
-Foam::tmp<Foam::volScalarField> Foam::XiEqModels::instability2XiEq::XiEq() const
-{
-    IOdictionary combustionProperties
-    (
-        IOobject
-        (
-            "combustionProperties",
-            Su_.mesh().time().constant(),
-            Su_.mesh(),
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE
-        )
-    );
-
-    ignition ign(combustionProperties, Su_.mesh().time(), Su_.mesh());
-
-    //const scalar ignTim = ign.sites()[0].tmIgn();
-    const scalar curTime = Su_.mesh().time().value();
-    const scalar deltaT = Su_.mesh().time().deltaTValue();
-    const scalar ignTim = curTime - deltaT - ign.sites()[0].time();
-
-    volScalarField turbXiEq(XiEqModel_->XiEq());
-
-    volScalarField XiEqIn1("XiEqIn1", 0.0*turbXiEq);
-
-    dimensionedScalar CIn("CIn", dimensionSet(0, -2, 1, 0, 0, 0, 0), CIn_);
-    dimensionedScalar ignTm("ignTm", dimTime, ignTim);
-    XiEqIn1 =  exp(CIn*Su_*Su_*ignTm) - 1.0;
-
-    return
-    (
-        1.0 + sqrt(XiEqInFade_*sqr(XiEqIn1) + sqr(turbXiEq - 1.0))
-    );
-}
-
-
-bool Foam::XiEqModels::instability2XiEq::read(const dictionary& XiEqProperties)
-{
-    XiEqModel::read(XiEqProperties);
-
-    XiEqModelCoeffs_.readEntry("defaultCIn", defaultCIn_);
-    XiEqModelCoeffs_.readEntry("XiEqInFade", XiEqInFade_);
-
-    return XiEqModel_->read(XiEqModelCoeffs_);
-}
-
-
-// ************************************************************************* //

applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/instabilityXiEq/instabilityXiEq.C

@@ -0,0 +1,84 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2015 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+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 3 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, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "instabilityXiEq.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace XiEqModels
+{
+    defineTypeNameAndDebug(instability, 0);
+    addToRunTimeSelectionTable(XiEqModel, instability, dictionary);
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::XiEqModels::instability::instability
+(
+    const dictionary& XiEqProperties,
+    const psiuReactionThermo& thermo,
+    const compressible::RASModel& turbulence,
+    const volScalarField& Su
+)
+:
+    XiEqModel(XiEqProperties, thermo, turbulence, Su),
+    XiEqIn(XiEqModelCoeffs_.get<scalar>("XiEqIn")),
+    XiEqModel_(XiEqModel::New(XiEqModelCoeffs_, thermo, turbulence, Su))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::XiEqModels::instability::~instability()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField> Foam::XiEqModels::instability::XiEq() const
+{
+    volScalarField turbXiEq(XiEqModel_->XiEq());
+    return XiEqIn/turbXiEq + turbXiEq;
+}
+
+
+bool Foam::XiEqModels::instability::read(const dictionary& XiEqProperties)
+{
+    XiEqModel::read(XiEqProperties);
+
+    XiEqModelCoeffs_.readEntry("XiEqIn", XiEqIn);
+
+    return XiEqModel_->read(XiEqModelCoeffs_);
+}
+
+
+// ************************************************************************* //

applications/solvers/combustion/PDRFoam/XiModels/XiEqModels/instabilityXiEq/instabilityXiEq.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2021 OpenCFD Ltd.
+    Copyright (C) 2011-2012 OpenFOAM Foundation
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -24,18 +24,20 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::XiEqModels::normBasicSubGrid
+    Foam::XiEqModels::instability
 
 Description
-
+    This is the equilibrium level of the flame wrinkling generated by
+    instability. It is a constant (default 2.5). It is used in
+    \link XiModel.H \endlink.
 
 SourceFiles
-    normBasicSubGrid.C
+    instability.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef normBasicSubGrid_H
-#define normBasicSubGrid_H
+#ifndef instability_H
+#define instability_H
 
 #include "XiEqModel.H"
 
@@ -47,44 +49,43 @@ namespace XiEqModels
 {
 
 /*---------------------------------------------------------------------------*\
-                          Class normBasicSubGrid Declaration
+                          Class instability Declaration
 \*---------------------------------------------------------------------------*/
 
-class normBasicSubGrid
+class instability
 :
     public XiEqModel
 {
-    // Private Data
+    // Private data
 
-        // Constants in the equilibrium Xp equation
-        scalar Cxpe1_;
-        scalar Cxpe2_;
-        scalar Cxpe3_;
-        scalar Cxpe4_;
+        //- Equilibrium Xi due to instability only
+        scalar XiEqIn;
+
+        //- Equilibrium Xi model due to all other effects
+        autoPtr<XiEqModel> XiEqModel_;
 
 
     // Private Member Functions
 
         //- No copy construct
-        normBasicSubGrid(const normBasicSubGrid&) = delete;
+        instability(const instability&) = delete;
 
         //- No copy assignment
-        void operator=(const normBasicSubGrid&) = delete;
+        void operator=(const instability&) = delete;
 
 
 public:
 
     //- Runtime type information
-    TypeName("normBasicSubGrid");
+    TypeName("instability");
 
 
     // Constructors
 
         //- Construct from components
-        normBasicSubGrid
+        instability
         (
             const dictionary& XiEqProperties,
-            const word& modelType,
             const psiuReactionThermo& thermo,
             const compressible::RASModel& turbulence,
             const volScalarField& Su
@@ -92,7 +93,7 @@ public:
 
 
     //- Destructor
-    virtual ~normBasicSubGrid();
+    virtual ~instability();
 
 
     // Member Functions
@@ -102,6 +103,7 @@ public:
 
         //- Update properties from given dictionary
         virtual bool read(const dictionary& XiEqProperties);
+
 };
 
 

applications/solvers/combustion/PDRFoam/XiModels/XiGModels/KTS/KTS.C

@@ -6,7 +6,6 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -46,17 +45,22 @@ namespace XiGModels
 Foam::XiGModels::KTS::KTS
 (
     const dictionary& XiGProperties,
-    const word& modelType,
     const psiuReactionThermo& thermo,
     const compressible::RASModel& turbulence,
     const volScalarField& Su
 )
 :
-    XiGModel(XiGProperties, modelType, thermo, turbulence, Su),
+    XiGModel(XiGProperties, thermo, turbulence, Su),
     GEtaCoef_(XiGModelCoeffs_.get<scalar>("GEtaCoef"))
 {}
 
 
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::XiGModels::KTS::~KTS()
+{}
+
+
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
 Foam::tmp<Foam::volScalarField> Foam::XiGModels::KTS::G() const
@@ -70,15 +74,6 @@ Foam::tmp<Foam::volScalarField> Foam::XiGModels::KTS::G() const
 }
 
 
-Foam::tmp<Foam::volScalarField> Foam::XiGModels::KTS::Db() const
-{
-    const objectRegistry& db = Su_.db();
-    const volScalarField& Db1 = db.lookupObject<volScalarField>("Db");
-    //return turbulence_.muEff();
-    return Db1;
-}
-
-
 bool Foam::XiGModels::KTS::read(const dictionary& XiGProperties)
 {
     XiGModel::read(XiGProperties);

applications/solvers/combustion/PDRFoam/XiModels/XiGModels/KTS/KTS.H

@@ -6,7 +6,6 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2012 OpenFOAM Foundation
-    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -28,8 +27,8 @@ Class
     Foam::XiGModels::KTS
 
 Description
-    Simple Kolmogorov time-scale (KTS) model for the flame-wrinkling
-    generation rate.
+    Simple Kolmogorov time-scale (KTS) model for the flame-wrinling generation
+    rate.
 
 SourceFiles
     KTS.C
@@ -56,7 +55,7 @@ class KTS
 :
     public XiGModel
 {
-    // Private Data
+    // Private data
 
         scalar GEtaCoef_;
 
@@ -82,7 +81,6 @@ public:
         KTS
         (
             const dictionary& XiGProperties,
-            const word& modelType,
             const psiuReactionThermo& thermo,
             const compressible::RASModel& turbulence,
             const volScalarField& Su
@@ -90,7 +88,7 @@ public:
 
 
     //- Destructor
-    virtual ~KTS() = default;
+    virtual ~KTS();
 
 
     // Member Functions
@@ -98,9 +96,6 @@ public:
         //- Return the flame-wrinkling generation rate
         virtual tmp<volScalarField> G() const;
 
-        //- Return the flame diffusivity
-        virtual tmp<volScalarField> Db() const;
-
         //- Update properties from given dictionary
         virtual bool read(const dictionary& XiGProperties);
 };

applications/solvers/combustion/PDRFoam/XiModels/XiGModels/XiGModel/XiGModel.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2011-2015 OpenFOAM Foundation
+    Copyright (C) 2011-2017 OpenFOAM Foundation
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -41,7 +41,6 @@ namespace Foam
 Foam::XiGModel::XiGModel
 (
     const dictionary& XiGProperties,
-    const word& modelType,
     const psiuReactionThermo& thermo,
     const compressible::RASModel& turbulence,
     const volScalarField& Su
@@ -51,7 +50,7 @@ Foam::XiGModel::XiGModel
     (
         XiGProperties.subDict
         (
-            XiGProperties.get<word>(modelType) + "Coeffs"
+            XiGProperties.get<word>("XiGModel") + "Coeffs"
         )
     ),
     thermo_(thermo),
@@ -70,7 +69,7 @@ Foam::XiGModel::~XiGModel()
 
 bool Foam::XiGModel::read(const dictionary& XiGProperties)
 {
-    XiGModelCoeffs_ = XiGProperties.subDict(type() + "Coeffs");
+    XiGModelCoeffs_ = XiGProperties.optionalSubDict(type() + "Coeffs");
 
     return true;
 }

applications/solvers/combustion/PDRFoam/XiModels/XiGModels/XiGModel/XiGModel.H

@@ -28,6 +28,9 @@ Class
 
 Description
     Base-class for all Xi generation models used by the b-Xi combustion model.
+    See Technical Report SH/RE/01R for details on the PDR modelling. For details
+    on the use of XiGModel see \link XiModel.H \endlink. The model available is
+    \link instabilityG.H \endlink
 
 SourceFiles
     XiGModel.C
@@ -53,9 +56,10 @@ namespace Foam
 
 class XiGModel
 {
+
 protected:
 
-    // Protected Data
+    // Protected data
 
         dictionary XiGModelCoeffs_;
 
@@ -64,7 +68,9 @@ protected:
         const volScalarField& Su_;
 
 
-    // Protected Member Functions
+private:
+
+    // Private Member Functions
 
         //- No copy construct
         XiGModel(const XiGModel&) = delete;
@@ -88,14 +94,12 @@ public:
             dictionary,
             (
                 const dictionary& XiGProperties,
-                const word& modelType,
                 const psiuReactionThermo& thermo,
                 const compressible::RASModel& turbulence,
                 const volScalarField& Su
             ),
             (
                 XiGProperties,
-                modelType,
                 thermo,
                 turbulence,
                 Su
@@ -109,7 +113,6 @@ public:
         static autoPtr<XiGModel> New
         (
             const dictionary& XiGProperties,
-            const word& modelType,
             const psiuReactionThermo& thermo,
             const compressible::RASModel& turbulence,
             const volScalarField& Su
@@ -122,7 +125,6 @@ public:
         XiGModel
         (
             const dictionary& XiGProperties,
-            const word& modelType,
             const psiuReactionThermo& thermo,
             const compressible::RASModel& turbulence,
             const volScalarField& Su
@@ -135,11 +137,14 @@ public:
 
     // Member Functions
 
-        //- Return the flame-wrinkling genration rate
+        //- Return the flame-wrinkling generation rate
         virtual tmp<volScalarField> G() const = 0;
 
         //- Return the flame diffusivity
-        virtual tmp<volScalarField> Db() const = 0;
+        virtual tmp<volScalarField> Db() const
+        {
+            return turbulence_.muEff();
+        }
 
         //- Update properties from given dictionary
         virtual bool read(const dictionary& XiGProperties) = 0;

applications/solvers/combustion/PDRFoam/XiModels/XiGModels/XiGModel/XiGModelNew.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2019-2021 OpenCFD Ltd.
+    Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -33,15 +33,14 @@ License
 Foam::autoPtr<Foam::XiGModel> Foam::XiGModel::New
 (
     const dictionary& dict,
-    const word& modelName,
     const psiuReactionThermo& thermo,
     const compressible::RASModel& turbulence,
     const volScalarField& Su
 )
 {
-    const word modelType(dict.get<word>(modelName));
+    const word modelType(dict.get<word>("XiGModel"));
 
-    Info<< "Selecting flame-wrinkling model for G " << modelType << endl;
+    Info<< "Selecting flame-wrinkling model " << modelType << endl;
 
     auto cstrIter = dictionaryConstructorTablePtr_->cfind(modelType);
 
@@ -50,16 +49,13 @@ Foam::autoPtr<Foam::XiGModel> Foam::XiGModel::New
         FatalIOErrorInLookup
         (
             dict,
-            modelName,
+            "XiGModel",
             modelType,
             *dictionaryConstructorTablePtr_
         ) << exit(FatalIOError);
     }
 
-    return autoPtr<XiGModel>(cstrIter()
-    (
-        dict, modelName, thermo, turbulence, Su)
-    );
+    return autoPtr<XiGModel>(cstrIter()(dict, thermo, turbulence, Su));
 }
 
 

applications/solvers/combustion/PDRFoam/XiModels/XiGModels/instability2G/instability2G.C

@@ -1,165 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2021 OpenCFD Ltd.
--------------------------------------------------------------------------------
-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 3 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, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "IFstream.H"
-#include "instability2G.H"
-#include "addToRunTimeSelectionTable.H"
-#include "fvCFD.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace XiGModels
-{
-    defineTypeNameAndDebug(instability2G, 0);
-    addToRunTimeSelectionTable(XiGModel, instability2G, dictionary);
-}
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::XiGModels::instability2G::instability2G
-(
-    const dictionary& XiGProperties,
-    const word& modelType,
-    const psiuReactionThermo& thermo,
-    const compressible::RASModel& turbulence,
-    const volScalarField& Su
-)
-:
-    XiGModel(XiGProperties, modelType, thermo, turbulence, Su),
-    saModel_
-    (
-        IOdictionary
-        (
-            IOobject
-            (
-                "combustionProperties",
-                Su.mesh().time().constant(),
-                Su.mesh(),
-                IOobject::MUST_READ
-            )
-        ),
-        thermo
-    ),
-    CIn_(saModel_.CIn()),
-    defaultCIn_(XiGModelCoeffs_.get<scalar>("defaultCIn")),
-    GInFade_(XiGModelCoeffs_.get<scalar>("GInFade")),
-    GInMult_(XiGModelCoeffs_.get<scalar>("GInMult")),
-    lambdaIn_("lambdaIn", XiGModelCoeffs_),
-    XiGModel_
-    (
-        XiGModel::New(XiGModelCoeffs_, modelType, thermo, turbulence, Su)
-    )
-{
-    if (CIn_ <= 0.0)
-    {
-        CIn_ = defaultCIn_;
-    }
-}
-
-
-// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
-
-Foam::tmp<Foam::volScalarField> Foam::XiGModels::instability2G::G() const
-{
-    IOdictionary combustionProperties
-    (
-        IOobject
-        (
-            "combustionProperties",
-            Su_.mesh().time().constant(),
-            Su_.mesh(),
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE
-        )
-    );
-
-    ignition ign(combustionProperties, Su_.mesh().time(), Su_.mesh());
-
-    scalar curTime = Su_.mesh().time().value();
-    scalar deltaT = Su_.mesh().time().deltaTValue();
-    const scalar ignTim = curTime - deltaT - ign.sites()[0].time();
-
-    volScalarField turbXiG(XiGModel_->G());
-
-    volScalarField GIn("GIn", 0.0*turbXiG);
-
-    forAll (GIn, i)
-    {
-        GIn[i] = CIn_*Su_[i]*Su_[i]*exp(CIn_*Su_[i]*Su_[i]*ignTim)*GInMult_;
-    }
-
-    dimensionedScalar CIn("CIn", dimensionSet(0, -2, 1, 0, 0, 0, 0), CIn_);
-    dimensionedScalar ignTm("ignTm", dimTime, ignTim);
-
-    GIn = CIn*Su_*Su_*exp(CIn*Su_*Su_*ignTm)*GInMult_;
-
-    GIn *=
-    (
-        GIn /
-        (
-            GIn
-          + GInFade_*turbXiG
-          + dimensionedScalar("GSmall", inv(dimTime), SMALL)
-        )
-    );
-
-    return (GIn + turbXiG);
-}
-
-
-Foam::tmp<Foam::volScalarField> Foam::XiGModels::instability2G::Db() const
-{
-    const objectRegistry& db = Su_.db();
-    const volScalarField& Xi = db.lookupObject<volScalarField>("Xi");
-    const volScalarField& rho = db.lookupObject<volScalarField>("rho");
-    const volScalarField& mgb = db.lookupObject<volScalarField>("mgb");
-    const volScalarField& Db1 = db.lookupObject<volScalarField>("Db");
-
-    //return  turbulence_.muEff()
-    return Db1
-        + rho*Su_*(Xi - 1.0)*mgb*(0.5*lambdaIn_)/(mgb + 1.0/lambdaIn_);
-}
-
-
-bool Foam::XiGModels::instability2G::read(const dictionary& XiGProperties)
-{
-    XiGModel::read(XiGProperties);
-
-    XiGModelCoeffs_.readEntry("defaultCIn", defaultCIn_);
-    XiGModelCoeffs_.readEntry("GInFade", GInFade_);
-    XiGModelCoeffs_.readEntry("GInMult", GInMult_);
-    XiGModelCoeffs_.readEntry("lambdaIn", lambdaIn_);
-
-    return true;
-}
-
-
-// ************************************************************************* //

applications/solvers/combustion/PDRFoam/XiModels/XiGModels/instability2G/instability2G.H

@@ -1,142 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2021 OpenCFD Ltd.
--------------------------------------------------------------------------------
-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 3 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, see <http://www.gnu.org/licenses/>.
-
-Class
-    Foam::XiGModels::instability2G
-
-Description
-    Flame-surface instabilityG flame-wrinkling generation rate coefficient model
-    used in \link XiModel.H \endlink.
-
-SourceFiles
-    instability2G.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef instability2G_H
-#define instability2G_H
-
-#include "laminarFlameSpeed.H"
-#include "SCOPELaminarFlameSpeed.H"
-#include "XiGModel.H"
-#include "ignitionSite.H"
-#include "ignition.H"
-#include "Time.H"
-#include "fvMesh.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace XiGModels
-{
-
-/*---------------------------------------------------------------------------*\
-                          Class instability2G Declaration
-\*---------------------------------------------------------------------------*/
-
-class instability2G
-:
-    public XiGModel
-{
-    // Private Data
-
-        laminarFlameSpeedModels::SCOPE saModel_;
-
-        // GIn (inituial instability G)divided by Su^2. Read from fuel file
-        scalar CIn_;
-
-        //- Default CIn if not in fuel file
-        scalar defaultCIn_;
-
-        // Determine how fast GIn fades out as turbulence starts
-        scalar GInFade_;
-
-        // Set GIn large so that XiEq determines Xi value.
-        // Son increase byfactor:
-        scalar GInMult_;
-
-        //- The instability2G length-scale
-        dimensionedScalar lambdaIn_;
-
-        //- Xi generation rate model due to all other processes
-        autoPtr<XiGModel> XiGModel_;
-
-
-    // Private Member Functions
-
-        //- No copy construct
-        instability2G(const instability2G&) = delete;
-
-        //- No copy assignment
-        void operator=(const instability2G&) = delete;
-
-
-public:
-
-    //- Runtime type information
-    TypeName("instability2G");
-
-
-    // Constructors
-
-        //- Construct from components
-        instability2G
-        (
-            const dictionary& XiGProperties,
-            const word& modelType,
-            const psiuReactionThermo& thermo,
-            const compressible::RASModel& turbulence,
-            const volScalarField& Su
-        );
-
-
-    //- Destructor
-    virtual ~instability2G() = default;
-
-
-    // Member Functions
-
-        //- Return the flame-wrinkling generation rate
-        virtual tmp<volScalarField> G() const;
-
-        //- Return the flame diffusivity
-        virtual tmp<volScalarField> Db() const;
-
-        //- Update properties from given dictionary
-        virtual bool read(const dictionary& XiGProperties);
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace XiGModels
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

applications/solvers/combustion/PDRFoam/XiModels/XiGModels/instabilityG/instabilityG.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2021 OpenCFD Ltd.
+    Copyright (C) 2011-2015 OpenFOAM Foundation
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -45,16 +45,21 @@ namespace XiGModels
 Foam::XiGModels::instabilityG::instabilityG
 (
     const dictionary& XiGProperties,
-    const word& modelType,
     const psiuReactionThermo& thermo,
     const compressible::RASModel& turbulence,
     const volScalarField& Su
 )
 :
-    XiGModel(XiGProperties, modelType, thermo, turbulence, Su),
-    GIn_("GIn", XiGModelCoeffs_),
-    lambdaIn_("lambdaIn", XiGModelCoeffs_),
-    XiGModel_(XiGModel::New(XiGModelCoeffs_,modelType,thermo, turbulence, Su))
+    XiGModel(XiGProperties, thermo, turbulence, Su),
+    GIn_("GIn", dimless/dimTime, XiGModelCoeffs_),
+    lambdaIn_("lambdaIn", dimLength, XiGModelCoeffs_),
+    XiGModel_(XiGModel::New(XiGModelCoeffs_, thermo, turbulence, Su))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::XiGModels::instabilityG::~instabilityG()
 {}
 
 
@@ -71,14 +76,11 @@ Foam::tmp<Foam::volScalarField> Foam::XiGModels::instabilityG::Db() const
 {
     const objectRegistry& db = Su_.db();
     const volScalarField& Xi = db.lookupObject<volScalarField>("Xi");
-    const volScalarField& Xp = db.lookupObject<volScalarField>("Xp");
     const volScalarField& rho = db.lookupObject<volScalarField>("rho");
     const volScalarField& mgb = db.lookupObject<volScalarField>("mgb");
-    const volScalarField& Db1 = db.lookupObject<volScalarField>("Db");
 
-    //// OLD return turbulence_.muEff()
-    return Db1
-        + rho*Su_*(Xp*Xi - 1.0)*mgb*(0.5*lambdaIn_)/(mgb + 1.0/lambdaIn_);
+    return XiGModel_->Db()
+        + rho*Su_*(Xi - 1.0)*mgb*(0.5*lambdaIn_)/(mgb + 1.0/lambdaIn_);
 }
 
 

applications/solvers/combustion/PDRFoam/XiModels/XiGModels/instabilityG/instabilityG.H

@@ -6,7 +6,6 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2012 OpenFOAM Foundation
-    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -31,6 +30,8 @@ Description
     Flame-surface instabilityG flame-wrinkling generation rate coefficient model
     used in \link XiModel.H \endlink.
 
+    See Technical Report SH/RE/01R for details on the PDR modelling.
+
 SourceFiles
     instabilityG.C
 
@@ -56,9 +57,9 @@ class instabilityG
 :
     public XiGModel
 {
-    // Private Data
+    // Private data
 
-        //- Flame instabilityG wrinling generation rate coefficient
+        //- Flame instabilityG wrinkling generation rate coefficient
         dimensionedScalar GIn_;
 
         //- InstabilityG length-scale
@@ -89,7 +90,6 @@ public:
         instabilityG
         (
             const dictionary& XiGProperties,
-            const word& modelType,
             const psiuReactionThermo& thermo,
             const compressible::RASModel& turbulence,
             const volScalarField& Su
@@ -97,7 +97,7 @@ public:
 
 
     //- Destructor
-    virtual ~instabilityG() = default;
+    virtual ~instabilityG();
 
 
     // Member Functions

applications/solvers/combustion/PDRFoam/XiModels/XiModel/XiModel.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2011-2015 OpenFOAM Foundation
+    Copyright (C) 2011-2017 OpenFOAM Foundation
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -73,19 +73,7 @@ Foam::XiModel::XiModel
             IOobject::AUTO_WRITE
         ),
         b.mesh()
-    ),
-    Xp_
-    (
-        IOobject
-        (
-            "Xp",
-            b.time().timeName(),
-            b.db(),
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        b.mesh()
-     )
+    )
 {}
 
 

applications/solvers/combustion/PDRFoam/XiModels/XiModel/XiModel.H

@@ -95,8 +95,8 @@ SourceFiles
 #include "IOdictionary.H"
 #include "psiuReactionThermo.H"
 #include "turbulentFluidThermoModel.H"
-#include "fvcDiv.H"
 #include "multivariateSurfaceInterpolationScheme.H"
+#include "fvcDiv.H"
 #include "runTimeSelectionTables.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -110,9 +110,10 @@ namespace Foam
 
 class XiModel
 {
+
 protected:
 
-    // Protected Data
+    // Protected data
 
         dictionary XiModelCoeffs_;
 
@@ -126,11 +127,10 @@ protected:
         //- Flame wrinkling field
         volScalarField Xi_;
 
-        //- Flame wrinkling field for the subgrid obstacles
-        volScalarField Xp_;
 
+private:
 
-    // Protected Member Functions
+    // Private Member Functions
 
         //- No copy construct
         XiModel(const XiModel&) = delete;
@@ -215,31 +215,19 @@ public:
             return Xi_;
         }
 
-        //- Return the obstacle-flame-wrinkling Xp
-        virtual const volScalarField& Xp() const
-        {
-             return Xp_;
-        }
-
         //- Return the flame diffusivity
         virtual tmp<volScalarField> Db() const
         {
-            const objectRegistry& db = Su_.db();
-            const volScalarField& Db1 = db.lookupObject<volScalarField>("Db");
-            //return turbulence_.muEff();
-            return Db1;
+            return turbulence_.muEff();
         }
 
         //- Add Xi to the multivariateSurfaceInterpolationScheme table
         //  if required
-        void addXi
+        virtual void addXi
         (
-            multivariateSurfaceInterpolationScheme<scalar>::fieldTable& fields
+            multivariateSurfaceInterpolationScheme<scalar>::fieldTable&
         )
-        {
-            fields.add(Xi_);
-            fields.add(Xp_);
-        }
+        {}
 
         //- Correct the flame-wrinkling Xi
         virtual void correct() = 0;

applications/solvers/combustion/PDRFoam/XiModels/algebraic/algebraic.C

@@ -0,0 +1,100 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2015 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+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 3 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, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "algebraic.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace XiModels
+{
+    defineTypeNameAndDebug(algebraic, 0);
+    addToRunTimeSelectionTable(XiModel, algebraic, dictionary);
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::XiModels::algebraic::algebraic
+(
+    const dictionary& XiProperties,
+    const psiuReactionThermo& thermo,
+    const compressible::RASModel& turbulence,
+    const volScalarField& Su,
+    const volScalarField& rho,
+    const volScalarField& b,
+    const surfaceScalarField& phi
+)
+:
+    XiModel(XiProperties, thermo, turbulence, Su, rho, b, phi),
+    XiShapeCoef(XiModelCoeffs_.get<scalar>("XiShapeCoef")),
+    XiEqModel_(XiEqModel::New(XiProperties, thermo, turbulence, Su)),
+    XiGModel_(XiGModel::New(XiProperties, thermo, turbulence, Su))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::XiModels::algebraic::~algebraic()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField> Foam::XiModels::algebraic::Db() const
+{
+    return XiGModel_->Db();
+}
+
+
+void Foam::XiModels::algebraic::correct()
+{
+    volScalarField XiEqEta(XiEqModel_->XiEq());
+    volScalarField GEta(XiGModel_->G());
+
+    volScalarField R(GEta*XiEqEta/(XiEqEta - 0.999));
+
+    volScalarField XiEqStar(R/(R - GEta));
+
+    Xi_ == 1.0 + (1.0 + (2*XiShapeCoef)*(0.5 - b_))*(XiEqStar - 1.0);
+}
+
+
+bool Foam::XiModels::algebraic::read(const dictionary& XiProperties)
+{
+    XiModel::read(XiProperties);
+
+    XiModelCoeffs_.readEntry("XiShapeCoef", XiShapeCoef);
+
+    return true;
+}
+
+
+// ************************************************************************* //

applications/solvers/combustion/PDRFoam/XiModels/algebraic/algebraic.H

@@ -0,0 +1,133 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2012 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+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 3 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, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::XiModels::algebraic
+
+Description
+    Simple algebraic model for Xi based on Gulders correlation
+    with a linear correction function to give a plausible profile for Xi.
+    See report TR/HGW/10 for details on the Weller two equations model.
+    See \link XiModel.H \endlink for more details on flame wrinkling modelling.
+
+SourceFiles
+    algebraic.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef algebraic_H
+#define algebraic_H
+
+#include "XiModel.H"
+#include "XiEqModel.H"
+#include "XiGModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace XiModels
+{
+
+/*---------------------------------------------------------------------------*\
+                          Class algebraic Declaration
+\*---------------------------------------------------------------------------*/
+
+class algebraic
+:
+    public XiModel
+{
+    // Private data
+
+        scalar XiShapeCoef;
+
+        autoPtr<XiEqModel> XiEqModel_;
+        autoPtr<XiGModel> XiGModel_;
+
+
+    // Private Member Functions
+
+        //- No copy construct
+        algebraic(const algebraic&) = delete;
+
+        //- No copy assignment
+        void operator=(const algebraic&) = delete;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("algebraic");
+
+
+    // Constructors
+
+        //- Construct from components
+        algebraic
+        (
+            const dictionary& XiProperties,
+            const psiuReactionThermo& thermo,
+            const compressible::RASModel& turbulence,
+            const volScalarField& Su,
+            const volScalarField& rho,
+            const volScalarField& b,
+            const surfaceScalarField& phi
+        );
+
+
+    //- Destructor
+    virtual ~algebraic();
+
+
+    // Member Functions
+
+        //- Return the flame diffusivity
+        virtual tmp<volScalarField> Db() const;
+
+        //- Correct the flame-wrinkling Xi
+        virtual void correct();
+
+        //- Update properties from given dictionary
+        virtual bool read(const dictionary& XiProperties);
+
+        //- Write fields of the XiEq model
+        virtual void writeFields()
+        {
+            XiEqModel_().writeFields();
+        }
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace XiModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/combustion/PDRFoam/XiModels/fixed/fixed.C

@@ -6,7 +6,6 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -24,26 +23,52 @@ License
     You should have received a copy of the GNU General Public License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
-Global
-    readTimeControls
-
-Description
-    Read the control parameters used by setDeltaT.
-    Modified for PDRFoam to include maxDeltaTRatio.
-
 \*---------------------------------------------------------------------------*/
 
-const bool adjustTimeStep =
-    runTime.controlDict().getOrDefault("adjustTimeStep", false);
+#include "fixed.H"
+#include "addToRunTimeSelectionTable.H"
 
-scalar maxCo =
-    runTime.controlDict().getOrDefault<scalar>("maxCo", 1);
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
-scalar maxDeltaT =
-    runTime.controlDict().getOrDefault<scalar>("maxDeltaT", GREAT);
+namespace Foam
+{
+namespace XiModels
+{
+    defineTypeNameAndDebug(fixed, 0);
+    addToRunTimeSelectionTable(XiModel, fixed, dictionary);
+}
+}
 
-scalar maxDeltaTRatio =
-    runTime.controlDict().getOrDefault<scalar>("maxDeltaTRatio", 1.2);
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::XiModels::fixed::fixed
+(
+    const dictionary& XiProperties,
+    const psiuReactionThermo& thermo,
+    const compressible::RASModel& turbulence,
+    const volScalarField& Su,
+    const volScalarField& rho,
+    const volScalarField& b,
+    const surfaceScalarField& phi
+)
+:
+    XiModel(XiProperties, thermo, turbulence, Su, rho, b, phi)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::XiModels::fixed::~fixed()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+bool Foam::XiModels::fixed::read(const dictionary& XiProperties)
+{
+    return XiModel::read(XiProperties);
+}
 
 
 // ************************************************************************* //

applications/solvers/combustion/PDRFoam/XiModels/fixed/fixed.H

@@ -0,0 +1,117 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2012 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+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 3 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, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::XiModels::fixed
+
+Description
+    Fixed value model for Xi. See \link XiModel.H \endlink for more details
+    on flame wrinkling modelling.
+
+SourceFiles
+    fixed.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef fixed_H
+#define fixed_H
+
+#include "XiModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace XiModels
+{
+
+/*---------------------------------------------------------------------------*\
+                          Class fixed Declaration
+\*---------------------------------------------------------------------------*/
+
+class fixed
+:
+    public XiModel
+{
+    // Private Member Functions
+
+        //- No copy construct
+        fixed(const fixed&) = delete;
+
+        //- No copy assignment
+        void operator=(const fixed&) = delete;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("fixed");
+
+
+    // Constructors
+
+        //- Construct from components
+        fixed
+        (
+            const dictionary& XiProperties,
+            const psiuReactionThermo& thermo,
+            const compressible::RASModel& turbulence,
+            const volScalarField& Su,
+            const volScalarField& rho,
+            const volScalarField& b,
+            const surfaceScalarField& phi
+        );
+
+
+    //- Destructor
+    virtual ~fixed();
+
+
+    // Member Functions
+
+        //- Correct the flame-wrinkling Xi
+        virtual void correct()
+        {}
+
+        //- Update properties from given dictionary
+        virtual bool read(const dictionary& XiProperties);
+
+        //- Write fields of the XiEq model
+        virtual void writeFields()
+        {}
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace XiModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/solvers/combustion/PDRFoam/XiModels/transport/transport.C

@@ -0,0 +1,147 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2015 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+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 3 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, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "transport.H"
+#include "surfaceInterpolate.H"
+#include "fvmDdt.H"
+#include "fvcLaplacian.H"
+#include "fvmDiv.H"
+#include "fvmSup.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace XiModels
+{
+    defineTypeNameAndDebug(transport, 0);
+    addToRunTimeSelectionTable(XiModel, transport, dictionary);
+}
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::XiModels::transport::transport
+(
+    const dictionary& XiProperties,
+    const psiuReactionThermo& thermo,
+    const compressible::RASModel& turbulence,
+    const volScalarField& Su,
+    const volScalarField& rho,
+    const volScalarField& b,
+    const surfaceScalarField& phi
+)
+:
+    XiModel(XiProperties, thermo, turbulence, Su, rho, b, phi),
+    XiShapeCoef(XiModelCoeffs_.get<scalar>("XiShapeCoef")),
+    XiEqModel_(XiEqModel::New(XiProperties, thermo, turbulence, Su)),
+    XiGModel_(XiGModel::New(XiProperties, thermo, turbulence, Su))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+Foam::XiModels::transport::~transport()
+{}
+
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volScalarField> Foam::XiModels::transport::Db() const
+{
+    return XiGModel_->Db();
+}
+
+
+void Foam::XiModels::transport::correct
+(
+    const fv::convectionScheme<scalar>& mvConvection
+)
+{
+    volScalarField XiEqEta(XiEqModel_->XiEq());
+    volScalarField GEta(XiGModel_->G());
+
+    volScalarField R(GEta*XiEqEta/(XiEqEta - 0.999));
+
+    volScalarField XiEqStar(R/(R - GEta));
+
+    volScalarField XiEq
+    (
+        1.0 + (1.0 + (2*XiShapeCoef)*(0.5 - b_))*(XiEqStar - 1.0)
+    );
+
+    volScalarField G(R*(XiEq - 1.0)/XiEq);
+
+    const objectRegistry& db = b_.db();
+    const volScalarField& betav = db.lookupObject<volScalarField>("betav");
+    const volScalarField& mgb = db.lookupObject<volScalarField>("mgb");
+    const surfaceScalarField& phiSt =
+        db.lookupObject<surfaceScalarField>("phiSt");
+    const volScalarField& Db = db.lookupObject<volScalarField>("Db");
+    const surfaceScalarField& nf = db.lookupObject<surfaceScalarField>("nf");
+
+    surfaceScalarField phiXi
+    (
+        "phiXi",
+        phiSt
+      + (
+          - fvc::interpolate(fvc::laplacian(Db, b_)/mgb)*nf
+          + fvc::interpolate(rho_)*fvc::interpolate(Su_*(1.0/Xi_ - Xi_))*nf
+        )
+    );
+
+    solve
+    (
+        betav*fvm::ddt(rho_, Xi_)
+      + mvConvection.fvmDiv(phi_, Xi_)
+      + fvm::div(phiXi, Xi_)
+      - fvm::Sp(fvc::div(phiXi), Xi_)
+     ==
+        betav*rho_*R
+      - fvm::Sp(betav*rho_*(R - G), Xi_)
+    );
+
+    // Correct boundedness of Xi
+    // ~~~~~~~~~~~~~~~~~~~~~~~~~
+    Xi_.max(1.0);
+    Xi_ = min(Xi_, 2.0*XiEq);
+}
+
+
+bool Foam::XiModels::transport::read(const dictionary& XiProperties)
+{
+    XiModel::read(XiProperties);
+
+    XiModelCoeffs_.readEntry("XiShapeCoef", XiShapeCoef);
+
+    return true;
+}
+
+
+// ************************************************************************* //

applications/solvers/combustion/PDRFoam/XiModels/transport/transport.H

@@ -6,7 +6,6 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -25,17 +24,21 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::XiModels::transportTwoEqs
+    Foam::XiModels::transport
 
 Description
+    Simple transport model for Xi based on Gulders correlation
+    with a linear correction function to give a plausible profile for Xi.
+    See report TR/HGW/10 for details on the Weller two equations model.
+    See \link XiModel.H \endlink for more details on flame wrinkling modelling.
 
 SourceFiles
-    transportTwoEqs.C
+    transport.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef transportTwoEqs_H
-#define transportTwoEqs_H
+#ifndef transport_H
+#define transport_H
 
 #include "XiModel.H"
 #include "XiEqModel.H"
@@ -49,58 +52,40 @@ namespace XiModels
 {
 
 /*---------------------------------------------------------------------------*\
-                          Class transportTwoEqs Declaration
+                          Class transport Declaration
 \*---------------------------------------------------------------------------*/
 
-class transportTwoEqs
+class transport
 :
     public XiModel
 {
     // Private data
 
-        scalar XiShapeCoef_;
-        scalar CpfiDot_;
-        scalar CpfiCross_;
-        scalar GEtaExp_;
+        scalar XiShapeCoef;
 
-        //- Maximum Lobs/CellWidth
-        scalar LOverCw_;
-
-        //- Equilibrium for Xi (turbulence)
         autoPtr<XiEqModel> XiEqModel_;
-
-        //- Generation for Xi (turbulence)
         autoPtr<XiGModel> XiGModel_;
 
-        //- Equilibrium for Xp (obstacles)
-        autoPtr<XiEqModel> XpEqModel_;
-
-        //- Generation for Xp (obstacles)
-        autoPtr<XiGModel> XpGModel_;
-
-        //- Dissipation length scale for subgrid obstacles
-        volScalarField Ep_;
-
 
     // Private Member Functions
 
         //- No copy construct
-        transportTwoEqs(const transportTwoEqs&) = delete;
+        transport(const transport&) = delete;
 
         //- No copy assignment
-        void operator=(const transportTwoEqs&) = delete;
+        void operator=(const transport&) = delete;
 
 
 public:
 
     //- Runtime type information
-    TypeName("transportTwoEqs");
+    TypeName("transport");
 
 
     // Constructors
 
         //- Construct from components
-        transportTwoEqs
+        transport
         (
             const dictionary& XiProperties,
             const psiuReactionThermo& thermo,
@@ -113,7 +98,7 @@ public:
 
 
     //- Destructor
-    virtual ~transportTwoEqs();
+    virtual ~transport();
 
 
     // Member Functions
@@ -121,6 +106,15 @@ public:
         //- Return the flame diffusivity
         virtual tmp<volScalarField> Db() const;
 
+        //- Add Xi to the multivariateSurfaceInterpolationScheme table
+        virtual void addXi
+        (
+            multivariateSurfaceInterpolationScheme<scalar>::fieldTable& fields
+        )
+        {
+            fields.add(Xi_);
+        }
+
         //- Correct the flame-wrinkling Xi
         virtual void correct()
         {
@@ -137,7 +131,6 @@ public:
         virtual void writeFields()
         {
             XiEqModel_().writeFields();
-            XpEqModel_().writeFields();
         }
 };
 

applications/solvers/combustion/PDRFoam/XiModels/transport/transportTwoEqs/transportTwoEqs.C

@@ -1,258 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2021 OpenCFD Ltd.
--------------------------------------------------------------------------------
-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 3 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, see <http://www.gnu.org/licenses/>.
-
-\*---------------------------------------------------------------------------*/
-
-#include "transportTwoEqs.H"
-#include "addToRunTimeSelectionTable.H"
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace XiModels
-{
-    defineTypeNameAndDebug(transportTwoEqs, 0);
-    addToRunTimeSelectionTable(XiModel, transportTwoEqs, dictionary);
-}
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-Foam::XiModels::transportTwoEqs::transportTwoEqs
-(
-    const dictionary& XiProperties,
-    const psiuReactionThermo& thermo,
-    const compressible::RASModel& turbulence,
-    const volScalarField& Su,
-    const volScalarField& rho,
-    const volScalarField& b,
-    const surfaceScalarField& phi
-)
-:
-    XiModel(XiProperties, thermo, turbulence, Su, rho, b, phi),
-    XiShapeCoef_(XiModelCoeffs_.get<scalar>("XiShapeCoef")),
-    CpfiDot_(XiModelCoeffs_.get<scalar>("CpfiDot")),
-    CpfiCross_(XiModelCoeffs_.get<scalar>("CpfiCross")),
-    GEtaExp_(XiModelCoeffs_.get<scalar>("GEtaExp")),
-    LOverCw_(XiModelCoeffs_.get<scalar>("LOverCw")),
-    XiEqModel_
-    (
-        XiEqModel::New(XiProperties, "XiEqModel", thermo, turbulence, Su)
-    ),
-    XiGModel_(XiGModel::New(XiProperties, "XiGModel", thermo, turbulence, Su)),
-    XpEqModel_
-    (
-        XiEqModel::New(XiProperties, "XpEqModel", thermo, turbulence, Su)
-    ),
-    XpGModel_
-    (
-        XiGModel::New(XiProperties, "XpGModel", thermo, turbulence, Su)
-    ),
-    Ep_
-    (
-        IOobject
-        (
-            "Ep",
-            b.time().timeName(),
-            b.db(),
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        b.mesh()
-    )
-{}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::XiModels::transportTwoEqs::~transportTwoEqs()
-{}
-
-
-// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
-
-
-Foam::tmp<Foam::volScalarField> Foam::XiModels::transportTwoEqs::Db() const
-{
-    return XiGModel_->Db();
-}
-
-
-void Foam::XiModels::transportTwoEqs::correct
-(
-    const fv::convectionScheme<scalar>& mvConvection
-)
-{
-    const volScalarField XiEqEta(XiEqModel_->XiEq());
-    volScalarField GEta(XiGModel_->G());
-
-    GEta *= max( 1.0, exp( GEtaExp_*(1.0 - (Xi_ - 1.0)/(XiEqEta - 0.999)))) ;
-
-    const volScalarField R(GEta*XiEqEta/(XiEqEta - 0.999));
-
-    const volScalarField XiEqStar(R/(R - GEta));
-
-    const volScalarField XiEq
-    (
-        1.0 + (1.0 + (2*XiShapeCoef_)*(0.5 - b_))*(XiEqStar - 1.0)
-    );
-
-    const volScalarField G(R*(XiEq - 1.0)/XiEq);
-
-
-    const objectRegistry& db = b_.db();
-    const volScalarField& betav = db.lookupObject<volScalarField>("betav");
-    const volScalarField& p = db.lookupObject<volScalarField>("p");
-    const volScalarField& mgb = db.lookupObject<volScalarField>("mgb");
-    const surfaceScalarField& phiSt =
-        db.lookupObject<surfaceScalarField>("phiSt");
-    const volScalarField& Db = db.lookupObject<volScalarField>("Db");
-    const surfaceScalarField& nf = db.lookupObject<surfaceScalarField>("nf");
-
-    surfaceScalarField phiXi
-    (
-        "phiXi",
-        phiSt
-      + (
-          - fvc::interpolate(fvc::laplacian(Db, b_)/mgb)*nf
-          + fvc::interpolate(rho_)
-          * fvc::interpolate(Su_*(1.0/(Xi_*Xp_) - (Xi_*Xp_)))*nf
-        )
-    );
-
-
-    dimensionedScalar zero
-    (
-        "zero",
-        dimensionSet(2, -6, -2, 0, 0, 0, 0),
-        scalar(0.0)
-    );
-
-    const volScalarField Gpfi
-    (
-        CpfiDot_
-      * sqrt(max(fvc::grad(rho_)&fvc::grad(p), zero))
-      / rho_*b_*(1.0-b_)
-
-      + CpfiCross_
-      * sqrt(mag(fvc::grad(rho_)^fvc::grad(p)))
-      / rho_*b_*(1.0-b_)
-    );
-
-    fvScalarMatrix  XiEqn_
-    (
-        betav*fvm::ddt(rho_, Xi_)
-      + mvConvection.fvmDiv(phi_, Xi_)
-      + fvm::div(phiXi, Xi_)
-      - fvm::Sp(fvc::div(phiXi), Xi_)
-     ==
-        betav*rho_*(R + Gpfi )
-      - fvm::Sp(betav*rho_*(R - G), Xi_)
-    );
-
-    XiEqn_.relax();
-    XiEqn_.solve();
-
-    // Correct boundedness of Xi
-    // ~~~~~~~~~~~~~~~~~~~~~~~~~
-    Xi_.max(1.0);
-    Xi_ = min(Xi_, 2.0*XiEq);
-
-    // Calculation of Xp generated by obstacles
-    volScalarField XpEqEta("XpEqEta",XpEqModel_->XiEq());
-
-    const volScalarField GpEta("GpEta", XpGModel_->G());
-
-    const volScalarField Rp("Rp", GpEta*XpEqEta/(XpEqEta - 0.999));
-
-    const volScalarField XpEq
-    (
-        "XpEq",
-        1.0 + (1.0 + (2*XiShapeCoef_)*(0.5 - b_))*(XpEqEta - 1.0)
-    );
-
-    const volScalarField Gpp("Gpp", Rp*(XpEq - 1.0)/XpEq);
-
-
-    fvScalarMatrix XpEqn_
-    (
-        betav*fvm::ddt(rho_, Xp_)
-      + mvConvection.fvmDiv(phi_, Xp_)
-      + fvm::div(phiXi, Xp_)
-      - fvm::Sp(fvc::div(phiXi), Xp_)
-     ==
-        betav*rho_*Rp
-      - fvm::Sp(betav*rho_*(Rp - Gpp), Xp_)
-    );
-
-    XpEqn_.relax();
-    XpEqn_.solve();
-
-    Xp_.max(1.0);
-    Xp_ = min(Xp_, 20.0*XpEq);
-
-    // Calculate Ep
-    const volScalarField& Lobs = db.lookupObject<volScalarField>("Lobs");
-    const scalarField Cw(pow(Su_.mesh().V(), 2.0/3.0));
-    volScalarField LI(Lobs);
-
-    LI.primitiveFieldRef() = max(LI.primitiveField(),LOverCw_*sqrt(Cw));
-
-    fvScalarMatrix EpEqn_
-    (
-        betav*fvm::ddt(rho_, Ep_)
-      + mvConvection.fvmDiv(phi_, Ep_)
-      + fvm::div(phiXi, Ep_)
-      - fvm::Sp(fvc::div(phiXi), Ep_)
-     ==
-        betav*rho_*Gpp*Xp_/LI
-      - fvm::Sp(betav*rho_*Rp, Ep_)
-   );
-
-    EpEqn_.relax();
-    EpEqn_.solve();
-
-    Ep_.max(0.0);
-    Ep_.min(100000.0);
-}
-
-
-bool Foam::XiModels::transportTwoEqs::read(const dictionary& XiProperties)
-{
-    XiModel::read(XiProperties);
-
-    XiModelCoeffs_.readEntry("XiShapeCoef", XiShapeCoef_);
-    XiModelCoeffs_.readEntry("CpfiDot", CpfiDot_);
-    XiModelCoeffs_.readEntry("CpfiCross", CpfiCross_);
-    XiModelCoeffs_.readEntry("GEtaExp", GEtaExp_);
-
-    return true;
-}
-
-
-// ************************************************************************* //

applications/solvers/combustion/PDRFoam/bEqn.H

@@ -5,20 +5,11 @@ tmp<fv::convectionScheme<scalar>> mvConvection
         mesh,
         fields,
         phi,
-        mesh.divScheme("div(phi,ft_b_h_hu)")
+        mesh.divScheme("div(phi,ft_b_ha_hau)")
     )
 );
 
-rhou = thermo.rhou();
-muu = thermo.muu();
-
-volScalarField Db("Db", turbulence->mut());
-const volScalarField DbThermal
-(
-    "DbAlpha",
-    turbulence->alphaEff()
-    //thermo.alpha() + turbulence->mut()/Prt
-);
+volScalarField Db("Db", turbulence->muEff());
 
 if (ign.ignited())
 {
@@ -26,31 +17,30 @@ if (ign.ignited())
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     Su = unstrainedLaminarFlameSpeed()();
 
-    // Wrinkling due to turbulence
     const volScalarField& Xi = flameWrinkling->Xi();
 
-    // Wrinkling due to obstacles
-    const volScalarField& Xp = flameWrinkling->Xp();
-
-    // Progress variable
+    // progress variable
     // ~~~~~~~~~~~~~~~~~
-    c = (1.0 - b);
+    volScalarField c("c", 1.0 - b);
 
+    // Unburnt gas density
+    // ~~~~~~~~~~~~~~~~~~~
+    volScalarField rhou(thermo.rhou());
 
     // Calculate flame normal etc.
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
+    //volVectorField n(fvc::grad(b));
     volVectorField n(fvc::reconstruct(fvc::snGrad(b)*mesh.magSf()));
 
     volScalarField mgb("mgb", mag(n));
-    volScalarField ddtb("ddtb",fvc::ddt(b));
 
     dimensionedScalar dMgb("dMgb", mgb.dimensions(), SMALL);
 
     {
         volScalarField bc(b*c);
 
-        dMgb += 1.0e-8*
+        dMgb += 1.0e-3*
             (bc*mgb)().weightedAverage(mesh.V())
            /(bc.weightedAverage(mesh.V()) + SMALL);
     }
@@ -64,7 +54,6 @@ if (ign.ignited())
     surfaceScalarField nf("nf", mesh.Sf() & nfVec);
     n /= mgb;
 
-
     #include "StCorr.H"
 
     // Calculate turbulent flame speed flux
@@ -75,6 +64,8 @@ if (ign.ignited())
 
     Db = flameWrinkling->Db();
 
+    // Create b equation
+    // ~~~~~~~~~~~~~~~~~
     fvScalarMatrix bEqn
     (
         betav*fvm::ddt(rho, b)
@@ -84,7 +75,6 @@ if (ign.ignited())
       - fvm::laplacian(Db, b)
      ==
         betav*fvOptions(rho, b)
-        // + betav*rhou*St*StCorr*mag(fvc::grad(b))
     );
 
 
@@ -92,16 +82,17 @@ if (ign.ignited())
     // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     #include "ignite.H"
 
+    // Solve for b
+    // ~~~~~~~~~~~
     bEqn.relax();
+
     fvOptions.constrain(bEqn);
+
     bEqn.solve();
-    b.max(0.0);
-    b.min(1.0);
 
     fvOptions.correct(b);
 
-    Info<< "min(b) = " << min(b).value() << nl
-        << "max(b) = " << max(b).value() << endl;
+    Info<< "min(b) = " << min(b).value() << endl;
 
     if (composition.contains("ft"))
     {
@@ -118,14 +109,8 @@ if (ign.ignited())
             << endl;
     }
 
-
-    // Correct the flame-wrinkling, smoothing
+    // Correct the flame-wrinkling
     flameWrinkling->correct(mvConvection);
-    St = Xp*Xi*Su*combustFlag;
 
-    if (PDRProperties.getOrDefault("smoothSt", false))
-    {
-        const scalar coeff = PDRProperties.get<scalar>("smoothStCoeff");
-        fvc::smooth(St, coeff);
-    }
+    St = Xi*Su;
 }

applications/solvers/combustion/PDRFoam/createFieldRefs.H

@@ -1,2 +1 @@
 const volScalarField& psi = thermo.psi();
-scalar pMin = thermo.getOrDefault<scalar>("pMin", 0);

applications/solvers/combustion/PDRFoam/createFields.H

@@ -1,303 +1,218 @@
-    Info<< "Reading thermophysical properties\n" << endl;
+Info<< "Reading thermophysical properties\n" << endl;
 
-    autoPtr<psiuReactionThermo> pThermo
+autoPtr<psiuReactionThermo> pThermo
+(
+    psiuReactionThermo::New(mesh)
+);
+psiuReactionThermo& thermo = pThermo();
+thermo.validate(args.executable(), "ha", "ea");
+
+basicSpecieMixture& composition = thermo.composition();
+
+volScalarField rho
+(
+    IOobject
     (
-        psiuReactionThermo::New(mesh)
-    );
-    psiuReactionThermo& thermo = pThermo();
-    thermo.validate(args.executable(), "ha", "ea");
-    basicMultiComponentMixture& composition = thermo.composition();
+        "rho",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    thermo.rho()
+);
 
-    volScalarField rho
+volScalarField& p = thermo.p();
+
+volScalarField& b = composition.Y("b");
+Info<< "min(b) = " << min(b).value() << endl;
+
+Info<< "\nReading field U\n" << endl;
+volVectorField U
+(
+    IOobject
     (
-        IOobject
-        (
-            "rho",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        thermo.rho()
-    );
-    volScalarField rhou
+        "U",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+#include "compressibleCreatePhi.H"
+
+mesh.setFluxRequired(p.name());
+
+Info<< "Creating turbulence model\n" << endl;
+autoPtr<compressible::RASModel> turbulence
+(
+    compressible::New<compressible::RASModel>
     (
-        IOobject
-        (
-            "rhou",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        thermo.rhou()
-    );
-    volScalarField muu
-    (
-        IOobject
-        (
-            "muu",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        thermo.muu()
-    );
-
-    volScalarField& p = thermo.p();
-
-    volScalarField& b = composition.Y("b");
-
-    volScalarField combustFlag
-    (
-        IOobject
-        (
-            "combustFlag",
-            mesh.facesInstance(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh
-    );
-
-    Info<< "Reading fields betav, betai\n" << endl;
-    volScalarField betav
-    (
-        IOobject
-        (
-            "betav",
-            mesh.facesInstance(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh
-    );
-
-    volSymmTensorField betai
-    (
-        IOobject
-        (
-            "betai",
-            mesh.facesInstance(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh
-    );
-
-    Info<< "Reading field Bv\n" << endl;
-    volSymmTensorField Bv
-    (
-        IOobject
-        (
-            "Bv",
-            mesh.facesInstance(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh
-    );
-
-    Info<< "Reading fields Nv,nsv\n" << endl;
-    volScalarField Nv
-    (
-        IOobject
-        (
-            "Nv",
-            mesh.facesInstance(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh
-    );
-
-    volSymmTensorField nsv
-    (
-        IOobject
-        (
-            "nsv",
-            mesh.facesInstance(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh
-    );
-
-    Info<< "Reading fields Blong\n" << endl;
-    volSymmTensorField Blong
-    (
-        IOobject
-        (
-            "Blong",
-            mesh.facesInstance(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh
-    );
-
-    Info<< "\nCreating field U\n" << endl;
-    volVectorField U
-    (
-        IOobject
-        (
-            "U",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-    #include "compressibleCreatePhi.H"
-    mesh.setFluxRequired(p.name());
-
-    Info<< "Creating turbulence model\n" << endl;
-    autoPtr<compressible::RASModel> turbulence
-    (
-        compressible::New<compressible::RASModel>
-        (
-            rho,
-            U,
-            phi,
-            thermo
-        )
-    );
-
-    #include "createDpdt.H"
-
-    #include "createK.H"
-
-    Info<< "Creating the unstrained laminar flame speed\n" << endl;
-    autoPtr<laminarFlameSpeed> unstrainedLaminarFlameSpeed
-    (
-        laminarFlameSpeed::New(thermo)
-    );
-
-    Info<< "Reading strained laminar flame speed field Su\n" << endl;
-    volScalarField Su
-    (
-        IOobject
-        (
-            "Su",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-
-    Info<< "Reading field Lobs\n" << endl;
-    volScalarField Lobs
-    (
-        IOobject
-        (
-            "Lobs",
-            mesh.facesInstance(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh
-    );
-
-    Info<< "Reading field CT\n" << endl;
-    volSymmTensorField CT
-    (
-        IOobject
-        (
-            "CT",
-            mesh.facesInstance(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh
-    );
-
-    IOdictionary PDRProperties
-    (
-        IOobject
-        (
-            "PDRProperties",
-            runTime.constant(),
-            mesh,
-            IOobject::MUST_READ_IF_MODIFIED,
-            IOobject::NO_WRITE
-        )
-    );
-
-    //- Create the drag model
-    autoPtr<PDRDragModel> drag = PDRDragModel::New
-    (
-        PDRProperties,
-        turbulence,
         rho,
         U,
-        phi
-    );
+        phi,
+        thermo
+    )
+);
 
-    //- Create the flame-wrinkling model
-    autoPtr<XiModel> flameWrinkling = XiModel::New
+
+Info<< "Creating the unstrained laminar flame speed\n" << endl;
+autoPtr<laminarFlameSpeed> unstrainedLaminarFlameSpeed
+(
+    laminarFlameSpeed::New(thermo)
+);
+
+
+Info<< "Reading strained laminar flame speed field Su\n" << endl;
+volScalarField Su
+(
+    IOobject
     (
-        PDRProperties,
-        thermo,
-        turbulence,
-        Su,
-        rho,
-        b,
-        phi
-    );
+        "Su",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
 
-    Info<< "Calculating turbulent flame speed field St\n" << endl;
-    volScalarField St
+Info<< "Reading field betav\n" << endl;
+volScalarField betav
+(
+    IOobject
     (
-        IOobject
-        (
-            "St",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::AUTO_WRITE
-        ),
-        flameWrinkling->Xi()*Su
-    );
+        "betav",
+        mesh.facesInstance(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE
+    ),
+    mesh
+);
 
-    multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
-
-    if (composition.contains("ft"))
-    {
-        fields.add(composition.Y("ft"));
-    }
-
-    fields.add(b);
-    fields.add(thermo.he());
-    fields.add(thermo.heu());
-    flameWrinkling->addXi(fields);
-
-    #include "createMRF.H"
-    #include "createFvOptions.H"
-
-    volScalarField c
+Info<< "Reading field Lobs\n" << endl;
+volScalarField Lobs
+(
+    IOobject
     (
-        IOobject
-        (
-            "c",
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        scalar(1) - b
-    );
+        "Lobs",
+        mesh.facesInstance(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE
+    ),
+    mesh
+);
+
+Info<< "Reading field CT\n" << endl;
+volSymmTensorField CT
+(
+    IOobject
+    (
+        "CT",
+        mesh.facesInstance(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE
+    ),
+    mesh
+);
+
+Info<< "Reading field Nv\n" << endl;
+volScalarField Nv
+(
+    IOobject
+    (
+        "Nv",
+        mesh.facesInstance(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE
+    ),
+    mesh
+);
+
+Info<< "Reading field nsv\n" << endl;
+volSymmTensorField nsv
+(
+    IOobject
+    (
+        "nsv",
+        mesh.facesInstance(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE
+    ),
+    mesh
+);
+
+IOdictionary PDRProperties
+(
+    IOobject
+    (
+        "PDRProperties",
+        runTime.constant(),
+        mesh,
+        IOobject::MUST_READ_IF_MODIFIED,
+        IOobject::NO_WRITE
+    )
+);
+
+//- Create the drag model
+autoPtr<PDRDragModel> drag = PDRDragModel::New
+(
+    PDRProperties,
+    *turbulence,
+    rho,
+    U,
+    phi
+);
+
+//- Create the flame-wrinkling model
+autoPtr<XiModel> flameWrinkling = XiModel::New
+(
+    PDRProperties,
+    thermo,
+    *turbulence,
+    Su,
+    rho,
+    b,
+    phi
+);
+
+Info<< "Calculating turbulent flame speed field St\n" << endl;
+volScalarField St
+(
+    IOobject
+    (
+        "St",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    flameWrinkling->Xi()*Su
+);
+
+
+multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
+
+if (composition.contains("ft"))
+{
+    fields.add(composition.Y("ft"));
+}
+
+fields.add(b);
+fields.add(thermo.he());
+fields.add(thermo.heu());
+flameWrinkling->addXi(fields);
+
+#include "createDpdt.H"
+
+#include "createK.H"
+
+#include "createMRF.H"
+#include "createFvOptions.H"

applications/solvers/combustion/PDRFoam/laminarFlameSpeed/SCOPE/SCOPELaminarFlameSpeed.C

@@ -5,8 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2011-2012 OpenFOAM Foundation
-    Copyright (C) 2021 OpenCFD Ltd.
+    Copyright (C) 2011-2017 OpenFOAM Foundation
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -48,34 +47,6 @@ namespace laminarFlameSpeedModels
 }
 
 
-// * * * * * * * * * * * * * * * Local Functions * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-static void readPolynomialCoeffs
-(
-    Polynomial<7>& coeffs,
-    const word& keyword,
-    const dictionary& dict
-)
-{
-    {
-        coeffs[0] = dict.get<scalar>(keyword + "0");
-        coeffs[1] = dict.get<scalar>(keyword + "1");
-        coeffs[2] = dict.get<scalar>(keyword + "2");
-        coeffs[3] = dict.get<scalar>(keyword + "3");
-        coeffs[4] = dict.get<scalar>(keyword + "4");
-        coeffs[5] = dict.get<scalar>(keyword + "5");
-        coeffs[6] = dict.get<scalar>(keyword + "6");
-    }
-
-    // TBD: support direct reading of all coeffs?
-}
-
-} // End namespace Foam
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::laminarFlameSpeedModels::SCOPE::polynomial::polynomial
@@ -104,22 +75,35 @@ Foam::laminarFlameSpeedModels::SCOPE::SCOPE
     (
         dictionary
         (
-            IFstream(dict.get<fileName>("fuelFile"))()
-        ).subDict(typeName + "Coeffs")
+            IFstream
+            (
+                dict.get<fileName>("fuelFile")
+            )()
+        ).optionalSubDict(typeName + "Coeffs")
+    ),
+    LFL_
+    (
+        coeffsDict_.getCompat<scalar>
+        (
+            "lowerFlammabilityLimit",
+            {{"lowerFlamabilityLimit", 1712}}
+        )
+    ),
+    UFL_
+    (
+        coeffsDict_.getCompat<scalar>
+        (
+            "upperFlammabilityLimit",
+            {{"upperFlamabilityLimit", 1712}}
+        )
     ),
-    LFL_(coeffsDict_.get<scalar>("lowerFlamabilityLimit")),
-    UFL_(coeffsDict_.get<scalar>("upperFlamabilityLimit")),
     SuPolyL_(coeffsDict_.subDict("lowerSuPolynomial")),
     SuPolyU_(coeffsDict_.subDict("upperSuPolynomial")),
-    Texp_(),
-    pexp_(),
-    CIn_(coeffsDict_.getOrDefault<scalar>("CIn", 0)),
+    Texp_(coeffsDict_.get<scalar>("Texp")),
+    pexp_(coeffsDict_.get<scalar>("pexp")),
     MaPolyL_(coeffsDict_.subDict("lowerMaPolynomial")),
     MaPolyU_(coeffsDict_.subDict("upperMaPolynomial"))
 {
-    readPolynomialCoeffs(Texp_, "Texp", coeffsDict_);
-    readPolynomialCoeffs(pexp_, "pexp", coeffsDict_);
-
     SuPolyL_.ll = max(SuPolyL_.ll, LFL_) + SMALL;
     SuPolyU_.ul = min(SuPolyU_.ul, UFL_) - SMALL;
 
@@ -132,7 +116,7 @@ Foam::laminarFlameSpeedModels::SCOPE::SCOPE
     if (debug)
     {
         Info<< "phi     Su  (T = Tref, p = pref)" << endl;
-        const label n = 200;
+        label n = 200;
         for (int i=0; i<n; i++)
         {
             scalar phi = (2.0*i)/n;
@@ -156,7 +140,7 @@ inline Foam::scalar Foam::laminarFlameSpeedModels::SCOPE::polyPhi
     const polynomial& a
 )
 {
-    const scalar x = phi - 1.0;
+    scalar x = phi - 1.0;
 
     return
         a[0]
@@ -211,34 +195,28 @@ inline Foam::scalar Foam::laminarFlameSpeedModels::SCOPE::SuRef
     }
 }
 
+
 inline Foam::scalar Foam::laminarFlameSpeedModels::SCOPE::Ma
 (
     scalar phi
 ) const
 {
-    if (phi < LFL_ || phi > UFL_)
+    if (phi < MaPolyL_.ll)
     {
-        // Return 0 beyond the flamibility limits
-        return scalar(0);
-    }
-    else if (phi < MaPolyL_.ll)
-    {
-        // Use linear interpolation between the low end of the
-        // lower polynomial and the lower flammability limit
-        return MaPolyL_.llv*(phi - LFL_)/(MaPolyL_.ll - LFL_);
+        // Beyond the lower limit assume Ma is constant
+        return MaPolyL_.llv;
     }
     else if (phi > MaPolyU_.ul)
     {
-        // Use linear interpolation between the upper end of the
-        // upper polynomial and the upper flammability limit
-        return MaPolyU_.ulv*(UFL_ - phi)/(UFL_ - MaPolyU_.ul);
+        // Beyond the upper limit assume Ma is constant
+        return MaPolyU_.ulv;
     }
-    else if (phi < MaPolyL_.lu)
+    else if (phi < SuPolyL_.lu)
     {
         // Evaluate the lower polynomial
         return polyPhi(phi, MaPolyL_);
     }
-    else if (phi > MaPolyU_.lu)
+    else if (phi > SuPolyU_.lu)
     {
         // Evaluate the upper polynomial
         return polyPhi(phi, MaPolyU_);
@@ -263,13 +241,10 @@ inline Foam::scalar Foam::laminarFlameSpeedModels::SCOPE::Su0pTphi
     scalar phi
 ) const
 {
-    constexpr scalar Tref = 300.0;
-    constexpr scalar pRef = 1.013e5;
+    static const scalar Tref = 300.0;
+    static const scalar pRef = 1.013e5;
 
-    const scalar Texp = Texp_.value(phi-1.0);
-    const scalar pexp = pexp_.value(phi-1.0);
-
-    return SuRef(phi)*pow((Tu/Tref), Texp)*pow((p/pRef), pexp);
+    return SuRef(phi)*pow((Tu/Tref), Texp_)*pow((p/pRef), pexp_);
 }
 
 
@@ -280,29 +255,35 @@ Foam::tmp<Foam::volScalarField> Foam::laminarFlameSpeedModels::SCOPE::Su0pTphi
     scalar phi
 ) const
 {
-    auto tSu0 = tmp<volScalarField>::New
+    tmp<volScalarField> tSu0
     (
-        IOobject
+        new volScalarField
         (
-            "Su0",
-            p.time().timeName(),
-            p.db(),
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        p.mesh(),
-        dimensionedScalar(dimVelocity, Zero)
+            IOobject
+            (
+                "Su0",
+                p.time().timeName(),
+                p.db(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            p.mesh(),
+            dimensionedScalar(dimVelocity, Zero)
+        )
     );
-    auto& Su0 = tSu0.ref();
+
+    volScalarField& Su0 = tSu0.ref();
 
     forAll(Su0, celli)
     {
         Su0[celli] = Su0pTphi(p[celli], Tu[celli], phi);
     }
 
-    forAll(Su0.boundaryField(), patchi)
+    volScalarField::Boundary& Su0Bf = Su0.boundaryFieldRef();
+
+    forAll(Su0Bf, patchi)
     {
-        scalarField& Su0p = Su0.boundaryFieldRef()[patchi];
+        scalarField& Su0p = Su0Bf[patchi];
         const scalarField& pp = p.boundaryField()[patchi];
         const scalarField& Tup = Tu.boundaryField()[patchi];
 
@@ -323,29 +304,35 @@ Foam::tmp<Foam::volScalarField> Foam::laminarFlameSpeedModels::SCOPE::Su0pTphi
     const volScalarField& phi
 ) const
 {
-    auto tSu0 = tmp<volScalarField>::New
+    tmp<volScalarField> tSu0
     (
-        IOobject
+        new volScalarField
         (
-            "Su0",
-            p.time().timeName(),
-            p.db(),
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        p.mesh(),
-        dimensionedScalar(dimVelocity, Zero)
+            IOobject
+            (
+                "Su0",
+                p.time().timeName(),
+                p.db(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            p.mesh(),
+            dimensionedScalar(dimVelocity, Zero)
+        )
     );
-    auto& Su0 = tSu0.ref();
+
+    volScalarField& Su0 = tSu0.ref();
 
     forAll(Su0, celli)
     {
         Su0[celli] = Su0pTphi(p[celli], Tu[celli], phi[celli]);
     }
 
-    forAll(Su0.boundaryField(), patchi)
+    volScalarField::Boundary& Su0Bf = Su0.boundaryFieldRef();
+
+    forAll(Su0Bf, patchi)
     {
-        scalarField& Su0p = Su0.boundaryFieldRef()[patchi];
+        scalarField& Su0p = Su0Bf[patchi];
         const scalarField& pp = p.boundaryField()[patchi];
         const scalarField& Tup = Tu.boundaryField()[patchi];
         const scalarField& phip = phi.boundaryField()[patchi];
@@ -366,40 +353,40 @@ Foam::tmp<Foam::volScalarField> Foam::laminarFlameSpeedModels::SCOPE::Su0pTphi
 }
 
 
-Foam::scalar Foam::laminarFlameSpeedModels::SCOPE::CIn() const noexcept
-{
-    return CIn_ ;
-}
-
-
 Foam::tmp<Foam::volScalarField> Foam::laminarFlameSpeedModels::SCOPE::Ma
 (
     const volScalarField& phi
 ) const
 {
-    auto tMa = tmp<volScalarField>::New
+    tmp<volScalarField> tMa
     (
-        IOobject
+        new volScalarField
         (
-            "Ma",
-            phi.time().timeName(),
-            phi.db(),
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        phi.mesh(),
-        dimensionedScalar(dimless, Zero)
+            IOobject
+            (
+                "Ma",
+                phi.time().timeName(),
+                phi.db(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            phi.mesh(),
+            dimensionedScalar(dimless, Zero)
+        )
     );
-    auto& ma = tMa.ref();
+
+    volScalarField& ma = tMa.ref();
 
     forAll(ma, celli)
     {
         ma[celli] = Ma(phi[celli]);
     }
 
-    forAll(ma.boundaryField(), patchi)
+    volScalarField::Boundary& maBf = ma.boundaryFieldRef();
+
+    forAll(maBf, patchi)
     {
-        scalarField& map = ma.boundaryFieldRef()[patchi];
+        scalarField& map = maBf[patchi];
         const scalarField& phip = phi.boundaryField()[patchi];
 
         forAll(map, facei)
@@ -423,8 +410,7 @@ Foam::laminarFlameSpeedModels::SCOPE::Ma() const
         (
             dimensionedScalar
             (
-                "stoichiometricAirFuelMassRatio",
-                psiuReactionThermo_
+                "stoichiometricAirFuelMassRatio", dimless, psiuReactionThermo_
             )*ft/(scalar(1) - ft)
         );
     }
@@ -432,18 +418,21 @@ Foam::laminarFlameSpeedModels::SCOPE::Ma() const
     {
         const fvMesh& mesh = psiuReactionThermo_.p().mesh();
 
-        return tmp<volScalarField>::New
+        return tmp<volScalarField>
         (
-            IOobject
+            new volScalarField
             (
-                "Ma",
-                mesh.time().timeName(),
+                IOobject
+                (
+                    "Ma",
+                    mesh.time().timeName(),
+                    mesh,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE
+                ),
                 mesh,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            mesh,
-            dimensionedScalar("Ma", dimless, Ma(equivalenceRatio_))
+                dimensionedScalar("Ma", dimless, Ma(equivalenceRatio_))
+            )
         );
     }
 }
@@ -462,8 +451,7 @@ Foam::laminarFlameSpeedModels::SCOPE::operator()() const
             psiuReactionThermo_.Tu(),
             dimensionedScalar
             (
-                "stoichiometricAirFuelMassRatio",
-                psiuReactionThermo_
+                "stoichiometricAirFuelMassRatio", dimless, psiuReactionThermo_
             )*ft/(scalar(1) - ft)
         );
     }

applications/solvers/combustion/PDRFoam/laminarFlameSpeed/SCOPE/SCOPELaminarFlameSpeed.H

@@ -6,7 +6,6 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2012 OpenFOAM Foundation
-    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -66,7 +65,6 @@ SourceFiles
 #define SCOPE_H
 
 #include "laminarFlameSpeed.H"
-#include "Polynomial.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -107,7 +105,7 @@ class SCOPE
             scalar lu;
 
             //- Construct from dictionary
-            explicit polynomial(const dictionary& polyDict);
+            polynomial(const dictionary& polyDict);
         };
 
 
@@ -125,14 +123,11 @@ class SCOPE
         //- Upper Su polynomial
         polynomial SuPolyU_;
 
-        //- Temperature correction exponents
-        Polynomial<7> Texp_;
+        //- Temperature correction exponent
+        scalar Texp_;
 
-        //- Pressure correction exponents
-        Polynomial<7> pexp_;
-
-        //- Additional coefficient
-        scalar CIn_;
+        //- Pressure correction exponent
+        scalar pexp_;
 
         //- Lower Ma polynomial
         polynomial MaPolyL_;
@@ -141,25 +136,25 @@ class SCOPE
         polynomial MaPolyU_;
 
 
-    // Private Member Functions
+    // Private member functions
 
         //- Polynomial evaluated from the given equivalence ratio
-        //- and polynomial coefficients
+        //  and polynomial coefficients
         static inline scalar polyPhi(scalar phi, const polynomial& a);
 
         //- Laminar flame speed evaluated from the given equivalence ratio
-        //- at the reference temperature and pressure
+        //  at the reference temperature and pressure
         inline scalar SuRef(scalar phi) const;
 
         //- Markstein evaluated from the given equivalence ratio
         inline scalar Ma(scalar phi) const;
 
         //- Laminar flame speed evaluated from the given equivalence ratio
-        //- corrected for temperature and pressure dependence
+        //  corrected for temperature and pressure dependence
         inline scalar Su0pTphi(scalar p, scalar Tu, scalar phi) const;
 
         //- Laminar flame speed evaluated from the given uniform
-        //- equivalence ratio corrected for temperature and pressure dependence
+        //  equivalence ratio corrected for temperature and pressure dependence
         tmp<volScalarField> Su0pTphi
         (
             const volScalarField& p,
@@ -168,7 +163,7 @@ class SCOPE
         ) const;
 
         //- Laminar flame speed evaluated from the given equivalence ratio
-        //- distribution corrected for temperature and pressure dependence
+        //  distribution corrected for temperature and pressure dependence
         tmp<volScalarField> Su0pTphi
         (
             const volScalarField& p,
@@ -177,14 +172,13 @@ class SCOPE
         ) const;
 
         //- Return the Markstein number
-        //- evaluated from the given equivalence ratio
+        //  evaluated from the given equivalence ratio
         tmp<volScalarField> Ma(const volScalarField& phi) const;
 
-        //- No copy construct
-        SCOPE(const SCOPE&) = delete;
+        //- Construct as copy (not implemented)
+        SCOPE(const SCOPE&);
 
-        //- No copy assignment
-        void operator=(const SCOPE&) = delete;
+        void operator=(const SCOPE&);
 
 
 public:
@@ -197,8 +191,8 @@ public:
         //- Construct from dictionary and psiuReactionThermo
         SCOPE
         (
-            const dictionary& dict,
-            const psiuReactionThermo& ct
+            const dictionary&,
+            const psiuReactionThermo&
         );
 
 
@@ -208,9 +202,6 @@ public:
 
     // Member functions
 
-        //- Return CIn
-        scalar CIn() const noexcept;
-
         //- Return the Markstein number
         tmp<volScalarField> Ma() const;
 

applications/solvers/combustion/PDRFoam/pEqn.H

@@ -1,7 +1,6 @@
 rho = thermo.rho();
 
 volScalarField rAU(1.0/UEqn.A());
-
 volVectorField HbyA(constrainHbyA(invA & UEqn.H(), U, p));
 
 if (pimple.transonic())
@@ -9,7 +8,12 @@ if (pimple.transonic())
     surfaceScalarField phid
     (
         "phid",
-        fvc::interpolate(psi)*(fvc::flux(HbyA))
+        fvc::interpolate(psi)
+       *(
+            fvc::flux(HbyA)
+          + fvc::interpolate(rho*rAU)*fvc::ddtCorr(rho, U, phi)
+           /fvc::interpolate(rho)
+        )
     );
 
     while (pimple.correctNonOrthogonal())
@@ -36,7 +40,10 @@ else
     surfaceScalarField phiHbyA
     (
         "phiHbyA",
-        fvc::flux(rho*HbyA)
+        (
+            fvc::flux(rho*HbyA)
+          + fvc::interpolate(rho*rAU)*fvc::ddtCorr(rho, U, phi)
+        )
     );
 
     while (pimple.correctNonOrthogonal())
@@ -59,8 +66,6 @@ else
     }
 }
 
-p.max(pMin);
-
 #include "rhoEqn.H"
 #include "continuityErrs.H"
 

applications/solvers/combustion/PDRFoam/rhoEqn.H

@@ -1,3 +1,38 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+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 3 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, see <http://www.gnu.org/licenses/>.
+
+Global
+    rhoEqn
+
+Description
+    Solve the continuity for density.
+
+\*---------------------------------------------------------------------------*/
+
 {
     solve(betav*fvm::ddt(rho) + fvc::div(phi));
 }
+
+// ************************************************************************* //

applications/solvers/combustion/PDRFoam/setDeltaT.H

@@ -36,8 +36,7 @@ Description
 if (adjustTimeStep)
 {
     scalar maxDeltaTFact = maxCo/(CoNum + StCoNum + SMALL);
-    scalar deltaTFact =
-        min(min(maxDeltaTFact, 1.0 + 0.1*maxDeltaTFact), maxDeltaTRatio);
+    scalar deltaTFact = min(min(maxDeltaTFact, 1.0 + 0.1*maxDeltaTFact), 1.2);
 
     runTime.setDeltaT
     (

applications/solvers/combustion/chemFoam/Make/options

@@ -6,13 +6,11 @@ EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude
+    -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude
 
 EXE_LIBS = \
     -lfiniteVolume \
     -lmeshTools \
     -lfluidThermophysicalModels \
     -lreactionThermophysicalModels \
-    -lchemistryModel \
-    -lthermophysicalProperties
+    -lchemistryModel

applications/solvers/compressible/rhoSimpleFoam/SRFrhoSimpleFoam/EEqn.H

@@ -0,0 +1,26 @@
+{
+    volScalarField& he = thermo.he();
+
+    fvScalarMatrix EEqn
+    (
+        fvm::div(phi, he)
+      + (
+            he.name() == "e"
+          ? fvc::div(phi, volScalarField("Ekp", 0.5*magSqr(U) + p/rho))
+          : fvc::div(phi, volScalarField("K", 0.5*magSqr(U)))
+        )
+      - fvm::laplacian(turbulence->alphaEff(), he)
+     ==
+        fvOptions(rho, he)
+    );
+
+    EEqn.relax();
+
+    fvOptions.constrain(EEqn);
+
+    EEqn.solve();
+
+    fvOptions.correct(he);
+
+    thermo.correct();
+}

applications/solvers/compressible/rhoSimpleFoam/SRFrhoSimpleFoam/Make/files

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

applications/solvers/compressible/rhoSimpleFoam/SRFrhoSimpleFoam/Make/options

@@ -0,0 +1,21 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/cfdTools \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
+
+EXE_LIBS = \
+    -lfiniteVolume \
+    -lfvOptions \
+    -lmeshTools \
+    -lsampling \
+    -lcompressibleTransportModels \
+    -lfluidThermophysicalModels \
+    -lspecie \
+    -lturbulenceModels \
+    -lcompressibleTurbulenceModels \
+    -latmosphericModels

applications/solvers/compressible/rhoSimpleFoam/SRFrhoSimpleFoam/SRFrhoSimpleFoam.C

@@ -0,0 +1,99 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2017 OpenFOAM Foundation
+-------------------------------------------------------------------------------
+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 3 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, see <http://www.gnu.org/licenses/>.
+
+Application
+    SRFSimpleFoam
+
+Group
+    grpIncompressibleSolvers
+
+Description
+    Steady-state solver for incompressible, turbulent flow of non-Newtonian
+    fluids in a single rotating frame.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "fluidThermo.H"
+#include "turbulentFluidThermoModel.H"
+#include "simpleControl.H"
+#include "pressureControl.H"
+#include "fvOptions.H"
+#include "SRFModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    argList::addNote
+    (
+        "Steady-state solver for compressible, turbulent flow"
+        " of non-Newtonian fluids in a single rotating frame."
+    );
+
+    #include "postProcess.H"
+
+    #include "addCheckCaseOptions.H"
+    #include "setRootCaseLists.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createControl.H"
+    #include "createFields.H"
+    #include "createFieldRefs.H"
+    #include "initContinuityErrs.H"
+
+    turbulence->validate();
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (simple.loop())
+    {
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        // --- Pressure-velocity SIMPLE corrector
+        {
+            #include "UrelEqn.H"
+            #include "EEqn.H"
+            #include "pEqn.H"
+        }
+
+        U = Urel + SRF->U();
+
+        turbulence->correct();
+
+        runTime.write();
+
+        runTime.printExecutionTime(Info);
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/compressible/rhoSimpleFoam/SRFrhoSimpleFoam/UrelEqn.H

@@ -0,0 +1,22 @@
+    // Relative momentum predictor
+
+    tmp<fvVectorMatrix> tUrelEqn
+    (
+        fvm::div(phi, Urel)
+      + turbulence->divDevRhoReff(Urel)
+      + rho*SRF->Su()
+     ==
+        fvOptions(rho, Urel)
+    );
+    fvVectorMatrix& UrelEqn = tUrelEqn.ref();
+
+    UrelEqn.relax();
+
+    fvOptions.constrain(UrelEqn);
+
+    if (simple.momentumPredictor())
+    {
+        solve(UrelEqn == -fvc::grad(p));
+
+        fvOptions.correct(Urel);
+    }

applications/solvers/compressible/rhoSimpleFoam/SRFrhoSimpleFoam/createFieldRefs.H

@@ -0,0 +1 @@
+const volScalarField& psi = thermo.psi();

applications/solvers/compressible/rhoSimpleFoam/SRFrhoSimpleFoam/createFields.H

@@ -0,0 +1,89 @@
+Info<< "Reading thermophysical properties\n" << endl;
+
+autoPtr<fluidThermo> pThermo
+(
+    fluidThermo::New(mesh)
+);
+fluidThermo& thermo = pThermo();
+thermo.validate(args.executable(), "h", "e");
+
+volScalarField& p = thermo.p();
+
+volScalarField rho
+(
+    IOobject
+    (
+        "rho",
+        runTime.timeName(),
+        mesh,
+        IOobject::READ_IF_PRESENT,
+        IOobject::AUTO_WRITE
+    ),
+    thermo.rho()
+);
+
+Info<< "Reading field Urel\n" << endl;
+volVectorField Urel
+(
+    IOobject
+    (
+        "Urel",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+
+Info<< "Reading/calculating face flux field phi\n" << endl;
+surfaceScalarField phi
+(
+    IOobject
+    (
+        "phi",
+        runTime.timeName(),
+        mesh,
+        IOobject::READ_IF_PRESENT,
+        IOobject::AUTO_WRITE
+    ),
+    linearInterpolate(rho*Urel) & mesh.Sf()
+);
+
+pressureControl pressureControl(p, rho, simple.dict());
+
+mesh.setFluxRequired(p.name());
+
+Info<< "Creating SRF model\n" << endl;
+autoPtr<SRF::SRFModel> SRF(SRF::SRFModel::New(Urel));
+
+// Construct the absolute velocity
+volVectorField U
+(
+    IOobject
+    (
+        "U",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    Urel + SRF->U()
+);
+
+
+Info<< "Creating turbulence model\n" << endl;
+autoPtr<compressible::turbulenceModel> turbulence
+(
+    compressible::turbulenceModel::New
+    (
+        rho,
+        U,
+        phi,
+        thermo
+    )
+);
+
+dimensionedScalar initialMass = fvc::domainIntegrate(rho);
+
+#include "createFvOptions.H"

applications/solvers/compressible/rhoSimpleFoam/SRFrhoSimpleFoam/pEqn.H

@@ -0,0 +1,109 @@
+volScalarField rAUrel(1.0/UrelEqn.A());
+surfaceScalarField rhorAUrelf("rhorAUf", fvc::interpolate(rho*rAUrel));
+volVectorField HbyA(constrainHbyA(rAUrel*UrelEqn.H(), Urel, p));
+tUrelEqn.clear();
+
+bool closedVolume = false;
+
+surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(rho)*fvc::flux(HbyA));
+
+// Update the pressure BCs to ensure flux consistency
+constrainPressure(p, rho, Urel, phiHbyA, rhorAUrelf);
+
+if (simple.transonic())
+{
+    surfaceScalarField phid
+    (
+        "phid",
+        (fvc::interpolate(psi)/fvc::interpolate(rho))*phiHbyA
+    );
+
+    phiHbyA -= fvc::interpolate(psi*p)*phiHbyA/fvc::interpolate(rho);
+
+    while (simple.correctNonOrthogonal())
+    {
+        fvScalarMatrix pEqn
+        (
+            fvc::div(phiHbyA)
+          + fvm::div(phid, p)
+          - fvm::laplacian(rhorAUrelf, p)
+          ==
+            fvOptions(psi, p, rho.name())
+        );
+
+        // Relax the pressure equation to ensure diagonal-dominance
+        pEqn.relax();
+
+        pEqn.setReference
+        (
+            pressureControl.refCell(),
+            pressureControl.refValue()
+        );
+
+        pEqn.solve();
+
+        if (simple.finalNonOrthogonalIter())
+        {
+            phi = phiHbyA + pEqn.flux();
+        }
+    }
+}
+else
+{
+    closedVolume = adjustPhi(phiHbyA, Urel, p);
+
+    while (simple.correctNonOrthogonal())
+    {
+        fvScalarMatrix pEqn
+        (
+            fvc::div(phiHbyA)
+          - fvm::laplacian(rhorAUrelf, p)
+          ==
+            fvOptions(psi, p, rho.name())
+        );
+
+        pEqn.setReference
+        (
+            pressureControl.refCell(),
+            pressureControl.refValue()
+        );
+
+        pEqn.solve();
+
+        if (simple.finalNonOrthogonalIter())
+        {
+            phi = phiHbyA + pEqn.flux();
+        }
+    }
+}
+
+#include "incompressible/continuityErrs.H"
+
+// Explicitly relax pressure for momentum corrector
+p.relax();
+
+Urel = HbyA - rAUrel*fvc::grad(p);
+Urel.correctBoundaryConditions();
+fvOptions.correct(Urel);
+
+bool pLimited = pressureControl.limit(p);
+
+// 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);
+}
+
+if (pLimited || closedVolume)
+{
+    p.correctBoundaryConditions();
+}
+
+rho = thermo.rho();
+
+if (!simple.transonic())
+{
+    rho.relax();
+}

applications/solvers/multiphase/interPhaseChangeFoam/alphaEqnSubCycle.H

@@ -9,12 +9,7 @@
         phic += (interface.cAlpha()*icAlpha)*fvc::interpolate(mag(U));
     }
 
-    volScalarField divU
-    (
-        mesh.moving()
-      ? fvc::div(phi + mesh.phi())
-      : fvc::div(phi)
-    );
+    volScalarField divU(fvc::div(phi));
 
     if (nAlphaSubCycles > 1)
     {

applications/test/HashPtrTable/Test-HashPtrTable.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011 OpenFOAM Foundation
-    Copyright (C) 2017-2021 OpenCFD Ltd.
+    Copyright (C) 2017-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -125,7 +125,6 @@ int main()
     myTable.set("pi", new double(3.14159));
     myTable.set("natlog", new double(2.718282));
     myTable.insert("sqrt2", autoPtr<double>::New(1.414214));
-    myTable.insert("euler", autoPtr<double>::New(0.577216));
 
     HashTable<std::unique_ptr<double>, word, string::hash> myTable1;
 
@@ -140,24 +139,14 @@ int main()
     myTable2.set("natlog", autoPtr<double>(new double(2.718282)));
     myTable2.insert("sqrt2", autoPtr<double>::New(1.414214));
 
-    Info<< "Initial table" << nl;
+    // Info<< myTable << endl;
     printTable(myTable);
-
-    Info<< "print" << nl;
     Info<< myTable2 << nl;
 
-    {
-        auto iter2 = myTable2.find("pi");
-        Info<< nl "Got pi=";
-        if (iter2.good())
-        {
-            Info<< **iter2 << nl;
-        }
-        else
-        {
-            Info<< "not-found" << nl;
-        }
-    }
+    auto iter2 = myTable2.find("pi");
+
+    Info<<"PI: " << **iter2 << nl;
+
 
     HashPtrTable<double> copy(myTable);
 
@@ -176,9 +165,6 @@ int main()
     myTable.erase("abc");
     myTable.erase("unknownKey");
 
-    myTable.release("euler");
-
-    Info<< "After erasure" << nl;
     printTable(myTable);
 
     HashPtrTable<double> moved(std::move(copy));

applications/test/MathFunctions/Make/files

@@ -1,3 +0,0 @@
-Test-MathFunctions.C
-
-EXE = $(FOAM_USER_APPBIN)/Test-MathFunctions

applications/test/MathFunctions/Make/options

@@ -1 +0,0 @@
-EXE_INC = -I../TestTools

applications/test/MathFunctions/Test-MathFunctions.C

@@ -1,179 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | www.openfoam.com
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-    Copyright (C) 2021 OpenCFD Ltd.
--------------------------------------------------------------------------------
-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 3 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, see <http://www.gnu.org/licenses/>.
-
-Application
-    Test-MathFunctions
-
-Description
-    Tests for \c Math namespace member functions
-    using \c doubleScalar base type.
-
-\*---------------------------------------------------------------------------*/
-
-#include "MathFunctions.H"
-#include "mathematicalConstants.H"
-#include "IOmanip.H"
-#include "TestTools.H"
-
-using namespace Foam;
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-// Execute each member functions of Foam::Math::, and print output
-template<class Type>
-void test_member_funcs(Type)
-{
-    Random rndGen(1234);
-    const label numberOfTests = 10000;
-    const scalar tolerance = 1e-3;
-
-    for (label i = 0; i < numberOfTests; ++i)
-    {
-        Info<< "# Inverse error functions:" << endl;
-        const Type x1 = rndGen.sample01<Type>();
-        const Type x2 = rndGen.sample01<Type>();
-        Info<< "    # Operands:" << nl
-            << "      x1 = " << x1 << nl
-            << "      x2 = " << x2 << endl;
-
-        cmp
-        (
-            "       # erf(erfinv(x1)) = x1 = ",
-            x1,
-            Foam::erf(Foam::Math::erfInv(x1)),
-            tolerance
-        );
-
-        cmp
-        (
-            "       # erf(erfinv(-x2)) = -x2 = ",
-            x2,
-            Foam::erf(Foam::Math::erfInv(x2)),
-            tolerance
-        );
-
-        Info<< "# Incomplete gamma functions:" << endl;
-        const Type a = rndGen.position(SMALL, Type(100));
-        const Type x = rndGen.position(Type(0), Type(100));
-        Info<< "    # Operands:" << nl
-            << "      a = " << a << nl
-            << "      x = " << x << endl;
-
-        cmp
-        (
-            "       # incGammaRatio_Q(a,x) + incGammaRatio_P(a,x) = 1 = ",
-            Foam::Math::incGammaRatio_Q(a,x) + Foam::Math::incGammaRatio_P(a,x),
-            scalar(1),
-            tolerance
-        );
-
-        cmp
-        (
-            "       # incGamma_Q(1, x) = exp(-x) = ",
-            Foam::Math::incGamma_Q(1, x),
-            Foam::exp(-x),
-            tolerance
-        );
-
-        cmp
-        (
-            "       # incGamma_Q(0.5, x) = sqrt(pi)*erfc(sqrt(x)) = ",
-            Foam::Math::incGamma_Q(0.5, x),
-            Foam::sqrt(Foam::constant::mathematical::pi)
-           *Foam::erfc(Foam::sqrt(x)),
-            tolerance
-        );
-
-        cmp
-        (
-            "       # incGamma_P(1,x) = 1 - exp(x) = ",
-            Foam::Math::incGamma_P(1, x),
-            1 - Foam::exp(-x),
-            tolerance
-        );
-
-        cmp
-        (
-            "       # incGamma_P(0.5, x) = sqrt(pi)*erf(sqrt(x)) = ",
-            Foam::Math::incGamma_P(0.5, x),
-            Foam::sqrt(Foam::constant::mathematical::pi)
-           *Foam::erf(Foam::sqrt(x)),
-            tolerance
-        );
-    }
-}
-
-
-// Do compile-time recursion over the given types
-template<std::size_t I = 0, typename... Tp>
-inline typename std::enable_if<I == sizeof...(Tp), void>::type
-run_tests(const std::tuple<Tp...>& types, const List<word>& typeID){}
-
-
-template<std::size_t I = 0, typename... Tp>
-inline typename std::enable_if<I < sizeof...(Tp), void>::type
-run_tests(const std::tuple<Tp...>& types, const List<word>& typeID)
-{
-    Info<< nl << "    ## Test member functions: "<< typeID[I] <<" ##" << nl;
-    test_member_funcs(std::get<I>(types));
-
-    run_tests<I + 1, Tp...>(types, typeID);
-}
-
-
-// * * * * * * * * * * * * * * * Main Program  * * * * * * * * * * * * * * * //
-
-int main()
-{
-    Info<< setprecision(15);
-
-    const std::tuple<doubleScalar> types
-    (
-        std::make_tuple(Zero)
-    );
-
-    const List<word> typeID
-    ({
-        "doubleScalar"
-    });
-
-    run_tests(types, typeID);
-
-
-    if (nFail_)
-    {
-        Info<< nl << "        #### "
-            << "Failed in " << nFail_ << " tests "
-            << "out of total " << nTest_ << " tests "
-            << "####\n" << endl;
-        return 1;
-    }
-
-    Info<< nl << "        #### Passed all " << nTest_ <<" tests ####\n" << endl;
-    return 0;
-}
-
-
-// ************************************************************************* //

applications/test/TestTools/TestTools.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2020-2021 OpenCFD Ltd.
+    Copyright (C) 2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -30,10 +30,6 @@ Description
 
 using namespace Foam;
 
-#include "floatScalar.H"
-#include "doubleScalar.H"
-#include "complex.H"
-#include "Matrix.H"
 #include "Random.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

applications/utilities/mesh/conversion/fluent3DMeshToFoam/fluent3DMeshToFoam.L

@@ -6,7 +6,6 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -890,7 +889,7 @@ int main(int argc, char *argv[])
             {
                 faces[facei].flip();
             }
-            std::swap(owner[facei], neighbour[facei]);
+            Swap(owner[facei], neighbour[facei]);
         }
     }
 

applications/utilities/mesh/conversion/ideasUnvToFoam/ideasUnvToFoam.C

@@ -6,7 +6,6 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -890,8 +889,8 @@ int main(int argc, char *argv[])
             if (own[facei] == -1 && nei[facei] != -1)
             {
                 // Boundary face with incorrect orientation
-                boundaryFaces[facei].flip();
-                std::swap(own[facei], nei[facei]);
+                boundaryFaces[facei] = boundaryFaces[facei].reverseFace();
+                Swap(own[facei], nei[facei]);
                 nReverse++;
             }
         }

applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshCalcDualMesh.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2012-2016 OpenFOAM Foundation
-    Copyright (C) 2018-2021 OpenCFD Ltd.
+    Copyright (C) 2018-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -2348,7 +2348,7 @@ void Foam::conformalVoronoiMesh::createFacesOwnerNeighbourAndPatches
                 {
                     forAll(procPatchFaces, fI)
                     {
-                        procPatchFaces[fI].flip();
+                        procPatchFaces[fI] = procPatchFaces[fI].reverseFace();
                     }
                 }
 

applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/conformalVoronoiMeshZones.C

@@ -6,7 +6,6 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2013-2016 OpenFOAM Foundation
-    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -83,15 +82,25 @@ void Foam::conformalVoronoiMesh::selectSeparatedCoupledFaces
     boolList& selected
 ) const
 {
-    for (const polyPatch& pp : mesh.boundaryMesh())
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
+
+    forAll(patches, patchi)
     {
         // Check all coupled. Avoid using .coupled() so we also pick up AMI.
-
-        const auto* cpp = isA<coupledPolyPatch>(pp);
-
-        if (cpp && (cpp->separated() || !cpp->parallel()))
+        if (isA<coupledPolyPatch>(patches[patchi]))
         {
-            SubList<bool>(selected, pp.size(), pp.start()) = true;
+            const coupledPolyPatch& cpp = refCast<const coupledPolyPatch>
+            (
+                patches[patchi]
+            );
+
+            if (cpp.separated() || !cpp.parallel())
+            {
+                forAll(cpp, i)
+                {
+                    selected[cpp.start()+i] = true;
+                }
+            }
         }
     }
 }

applications/utilities/mesh/generation/foamyMesh/conformalVoronoiMesh/conformalVoronoiMesh/indexedCell/indexedCellI.H

@@ -6,7 +6,6 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2012-2016 OpenFOAM Foundation
-    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -326,7 +325,7 @@ inline Foam::tetCell CGAL::indexedCell<Gt, Cb>::vertexGlobalIndices
         {
             if (tVGI[j - 1] > tVGI[j])
             {
-                std::swap(tVGI[j - 1], tVGI[j]);
+                Foam::Swap(tVGI[j - 1], tVGI[j]);
             }
         }
     }
@@ -355,8 +354,8 @@ CGAL::indexedCell<Gt, Cb>::globallyOrderedCellVertices
         {
             if (tVGI[j - 1] > tVGI[j])
             {
-                std::swap(tVGI[j - 1], tVGI[j]);
-                std::swap(vertexMap[j - 1], vertexMap[j]);
+                Foam::Swap(tVGI[j - 1], tVGI[j]);
+                Foam::Swap(vertexMap[j - 1], vertexMap[j]);
             }
         }
     }

applications/utilities/mesh/manipulation/createPatch/createPatch.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
-    Copyright (C) 2016-2021 OpenCFD Ltd.
+    Copyright (C) 2016-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -205,14 +205,18 @@ void filterPatches(polyMesh& mesh, const wordHashSet& addedPatchNames)
 // Dump for all patches the current match
 void dumpCyclicMatch(const fileName& prefix, const polyMesh& mesh)
 {
-    for (const polyPatch& pp : mesh.boundaryMesh())
-    {
-        const cyclicPolyPatch* cpp = isA<cyclicPolyPatch>(pp);
+    const polyBoundaryMesh& patches = mesh.boundaryMesh();
 
-        if (cpp && cpp->owner())
+    forAll(patches, patchi)
+    {
+        if
+        (
+            isA<cyclicPolyPatch>(patches[patchi])
+         && refCast<const cyclicPolyPatch>(patches[patchi]).owner()
+        )
         {
-            const auto& cycPatch = *cpp;
-            const auto& nbrPatch = cycPatch.neighbPatch();
+            const cyclicPolyPatch& cycPatch =
+                refCast<const cyclicPolyPatch>(patches[patchi]);
 
             // Dump patches
             {
@@ -227,6 +231,7 @@ void dumpCyclicMatch(const fileName& prefix, const polyMesh& mesh)
                 );
             }
 
+            const cyclicPolyPatch& nbrPatch = cycPatch.neighbPatch();
             {
                 OFstream str(prefix+nbrPatch.name()+".obj");
                 Pout<< "Dumping " << nbrPatch.name()
@@ -320,16 +325,22 @@ void syncPoints
 
     if (Pstream::parRun())
     {
-        const labelList& procPatches = mesh.globalData().processorPatches();
-
         // Send
-        for (const label patchi : procPatches)
+
+        forAll(patches, patchi)
         {
             const polyPatch& pp = patches[patchi];
-            const auto& procPatch = refCast<const processorPolyPatch>(pp);
 
-            if (pp.nPoints() && procPatch.owner())
+            if
+            (
+                isA<processorPolyPatch>(pp)
+             && pp.nPoints() > 0
+             && refCast<const processorPolyPatch>(pp).owner()
+            )
             {
+                const processorPolyPatch& procPatch =
+                    refCast<const processorPolyPatch>(pp);
+
                 // Get data per patchPoint in neighbouring point numbers.
                 pointField patchInfo(procPatch.nPoints(), nullValue);
 
@@ -357,13 +368,20 @@ void syncPoints
 
         // Receive and set.
 
-        for (const label patchi : procPatches)
+        forAll(patches, patchi)
         {
             const polyPatch& pp = patches[patchi];
-            const auto& procPatch = refCast<const processorPolyPatch>(pp);
 
-            if (pp.nPoints() && !procPatch.owner())
+            if
+            (
+                isA<processorPolyPatch>(pp)
+             && pp.nPoints() > 0
+             && !refCast<const processorPolyPatch>(pp).owner()
+            )
             {
+                const processorPolyPatch& procPatch =
+                    refCast<const processorPolyPatch>(pp);
+
                 pointField nbrPatchInfo(procPatch.nPoints());
                 {
                     // We do not know the number of points on the other side
@@ -401,19 +419,22 @@ void syncPoints
     }
 
     // Do the cyclics.
-    for (const polyPatch& pp : patches)
+    forAll(patches, patchi)
     {
-        const cyclicPolyPatch* cpp = isA<cyclicPolyPatch>(pp);
+        const polyPatch& pp = patches[patchi];
 
-        if (cpp && cpp->owner())
+        if
+        (
+            isA<cyclicPolyPatch>(pp)
+         && refCast<const cyclicPolyPatch>(pp).owner()
+        )
         {
-            // Owner does all.
-
-            const auto& cycPatch = *cpp;
-            const auto& nbrPatch = cycPatch.neighbPatch();
+            const cyclicPolyPatch& cycPatch =
+                refCast<const cyclicPolyPatch>(pp);
 
             const edgeList& coupledPoints = cycPatch.coupledPoints();
             const labelList& meshPts = cycPatch.meshPoints();
+            const cyclicPolyPatch& nbrPatch = cycPatch.neighbPatch();
             const labelList& nbrMeshPts = nbrPatch.meshPoints();
 
             pointField half0Values(coupledPoints.size());

applications/utilities/mesh/manipulation/moveDynamicMesh/moveDynamicMesh.C

@@ -117,37 +117,39 @@ void writeWeights(const polyMesh& mesh)
 
     for (const polyPatch& pp : mesh.boundaryMesh())
     {
-        const auto* cpp = isA<cyclicAMIPolyPatch>(pp);
-
-        if (cpp && cpp->owner())
+        if (isA<cyclicAMIPolyPatch>(pp))
         {
-            const auto& cycPatch = *cpp;
-            const auto& nbrPatch = cycPatch.neighbPatch();
+            const cyclicAMIPolyPatch& cpp =
+                refCast<const cyclicAMIPolyPatch>(pp);
 
-            const AMIPatchToPatchInterpolation& ami = cycPatch.AMI();
+            if (cpp.owner())
+            {
+                Info<< "Calculating AMI weights between owner patch: "
+                    << cpp.name() << " and neighbour patch: "
+                    << cpp.neighbPatch().name() << endl;
 
-            Info<< "Calculating AMI weights between owner patch: "
-                << cycPatch.name() << " and neighbour patch: "
-                << nbrPatch.name() << endl;
+                const AMIPatchToPatchInterpolation& ami =
+                    cpp.AMI();
 
-            writeWeights
-            (
-                mesh,
-                ami.tgtWeightsSum(),
-                nbrPatch,
-                outputDir,
-                "patch" + Foam::name(pp.index()) + "-tgt",
-                mesh.time()
-            );
-            writeWeights
-            (
-                mesh,
-                ami.srcWeightsSum(),
-                cycPatch,
-                outputDir,
-                "patch" + Foam::name(pp.index()) + "-src",
-                mesh.time()
-            );
+                writeWeights
+                (
+                    mesh,
+                    ami.tgtWeightsSum(),
+                    cpp.neighbPatch(),
+                    outputDir,
+                    "patch" + Foam::name(pp.index()) + "-tgt",
+                    mesh.time()
+                );
+                writeWeights
+                (
+                    mesh,
+                    ami.srcWeightsSum(),
+                    cpp,
+                    outputDir,
+                    "patch" + Foam::name(pp.index()) + "-src",
+                    mesh.time()
+                );
+            }
         }
     }
 }

applications/utilities/mesh/manipulation/renumberMesh/renumberMesh.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2016-2021 OpenCFD Ltd.
+    Copyright (C) 2016-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -441,7 +441,7 @@ autoPtr<mapPolyMesh> reorderMesh
         if (nei < own)
         {
             newFaces[facei].flip();
-            std::swap(newOwner[facei], newNeighbour[facei]);
+            Swap(newOwner[facei], newNeighbour[facei]);
             flipFaceFlux.insert(facei);
         }
     }

applications/utilities/parallelProcessing/decomposePar/domainDecompositionTemplates.C

@@ -6,7 +6,6 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2014-2016 OpenFOAM Foundation
-    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -45,18 +44,22 @@ void Foam::domainDecomposition::processInterCyclics
     // Processor boundaries from split cyclics
     forAll(patches, patchi)
     {
-        const auto& pp = patches[patchi];
-        const auto* cpp = isA<cyclicPolyPatch>(pp);
-
-        if (cpp && cpp->owner() == owner)
+        if (isA<cyclicPolyPatch>(patches[patchi]))
         {
-            // cyclic: check opposite side on this processor
-            const auto& cycPatch = *cpp;
-            const auto& nbrPatch = cycPatch.neighbPatch();
+            const cyclicPolyPatch& pp = refCast<const cyclicPolyPatch>
+            (
+                patches[patchi]
+            );
+
+            if (pp.owner() != owner)
+            {
+                continue;
+            }
 
             // cyclic: check opposite side on this processor
             const labelUList& patchFaceCells = pp.faceCells();
-            const labelUList& nbrPatchFaceCells = nbrPatch.faceCells();
+            const labelUList& nbrPatchFaceCells =
+                pp.neighbPatch().faceCells();
 
             // Store old sizes. Used to detect which inter-proc patches
             // have been added to.

applications/utilities/thermophysical/chemkinToFoam/Make/options

@@ -1,11 +1,9 @@
 EXE_INC = \
     -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude
+    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude
 
 EXE_LIBS = \
     -lreactionThermophysicalModels \
     -lcompressibleTransportModels \
     -lfluidThermophysicalModels \
-    -lspecie \
-    -lthermophysicalProperties
+    -lspecie

bin/paraFoam

@@ -7,7 +7,7 @@
 #    \\/     M anipulation  |
 #------------------------------------------------------------------------------
 #     Copyright (C) 2011-2016 OpenFOAM Foundation
-#     Copyright (C) 2016-2021 OpenCFD Ltd.
+#     Copyright (C) 2016-2020 OpenCFD Ltd.
 #------------------------------------------------------------------------------
 # License
 #     This file is part of OpenFOAM, distributed under GPL-3.0-or-later.
@@ -140,8 +140,8 @@ do
             esac
         else
             echo "Ignore bad/invalid plugin-path: $pluginPath" 1>&2
+            unset pluginPath
         fi
-        unset pluginPath
         ;;
 
     -block*)  # Silently accepts -blockMesh
@@ -238,11 +238,8 @@ then
             ;;
         esac
 
-        if [ -z "$pluginError" ]
+        if [ -n "$pluginError" ]
         then
-            # Ensure plugin is also in the lib-path
-            LD_LIBRARY_PATH="${PV_PLUGIN_PATH}:$LD_LIBRARY_PATH"
-        else
             cat<< NO_PLUGIN 1>&2
 $pluginError
     See '${0##*/} -help-build' for more information

bin/tools/foamCreateModuleInclude

@@ -152,7 +152,7 @@ then
 
     # Check that it appears to be an OpenFOAM installation
     # could also check [ -d "$projectDir/META-INFO" ]
-    if [ -d "$projectDir" ] && [ -f "$projectDir/etc/bashrc" ]
+    if [ -d "$projectDir" ] && [ -f "etc/bashrc" ]
     then
         echo "Appears to be an OpenFOAM installation" 1>&2
     else
@@ -276,52 +276,17 @@ fi
 echo "Using openfoam: $WM_PROJECT_DIR" 1>&2
 echo "==> $moduleOutput" 1>&2
 
-# Check if directory is the OpenFOAM or ThirdParty tree
-dirInTree()
-{
-    local dir="$1"
-    [ "${dir#$WM_PROJECT_DIR}" != "$dir" ] || \
-    [ "${dir#$WM_THIRD_PARTY_DIR}" != "$dir" ]
-}
+# Remove user directories as being unreliable
 
-# User directories are entirely unreliable
-foamOldDirs="$FOAM_USER_APPBIN $FOAM_USER_LIBBIN"
-
-# Site locations must be within the OpenFOAM/ThirdParty tree
-for dir in "$WM_PROJECT_SITE" "$FOAM_SITE_APPBIN" "$FOAM_SITE_LIBBIN"
-do
-    if [ -d "$dir" ] && dirInTree "$dir"
-    then
-        continue
-    fi
-    # Discard non-directories or non-OpenFOAM/ThirdParty locations
-    foamOldDirs="$foamOldDirs $dir"
-done
-
-# Changes associated with the '[A-Z].*_ARCH_PATH' variables must be
-# within the OpenFOAM/ThirdParty tree
-# - avoids adding items that may have come from another module
-for envname in $(env | sed -n -e 's/^\(.*ARCH_PATH\)=.*$/\1/p')
-do
-    eval "dir=\$$envname" || continue
-    if [ -d "$dir" ] && dirInTree "$dir"
-    then
-        continue
-    fi
-    # Discard non-directories or non-OpenFOAM/ThirdParty locations
-    foamOldDirs="$foamOldDirs $dir"
-done
+foamOldDirs="$ADIOS2_ARCH_PATH $BOOST_ARCH_PATH $CGAL_ARCH_PATH $FFTW_ARCH_PATH \
+$MESA_ARCH_PATH $LLVM_ARCH_PATH \
+$MPI_ARCH_PATH $SCOTCH_ARCH_PATH \
+$FOAM_SITE_APPBIN $FOAM_SITE_LIBBIN $WM_PROJECT_SITE \
+$FOAM_USER_APPBIN $FOAM_USER_LIBBIN"
 
 if [ "${optParaview:-false}" != true ]
 then
-    for dir in "$ParaView_DIR"
-    do
-        if [ -d "$dir" ] && dirInTree "$dir"
-        then
-            continue
-        fi
-        foamOldDirs="$foamOldDirs $ParaView_DIR"
-    done
+    foamOldDirs="$foamOldDirs $ParaView_DIR"
 fi
 
 foamClean="$WM_PROJECT_DIR/bin/foamCleanPath"
@@ -371,14 +336,6 @@ unset FOAM_SITE_APPBIN FOAM_SITE_LIBBIN WM_PROJECT_SITE
 # Should work without (use default)
 unset WM_DIR
 
-# Remove any sourcing cruft
-# Remove: 'BASH_.*'
-
-for envname in $(env | sed -n -e 's/^\(BASH_.*\)=.*$/\1/p')
-do
-    eval "unset $envname"
-done
-
 # Third-party cruft - only used for orig compilation
 # - as per spack installation
 

bin/tools/foamPackRelease

@@ -30,12 +30,11 @@
 #
 #     foamPackRelease -with-api=1912 -pkg-modules origin/develop
 #
-# Debian-style, without OpenFOAM sub-directory
+# Debian-style without OpenFOAM sub-directory
 #
 #     foamPackRelease -debian origin/develop
-#
-#     == -debian=openfoam_{version}
-#     == -name=openfoam_{api}.{patch} -no-prefix
+#     foamPackRelease -debian=openfoam2002_200129 origin/develop
+#     foamPackRelease -name=openfoam2002_200129.orig -no-prefix origin/develop
 #
 #------------------------------------------------------------------------------
 Script="${0##*/}"
@@ -301,14 +300,18 @@ fi
 
 case "$tarName" in
 (debian)
-    tarName="openfoam_${packageApi}"
+    tarName="openfoam${packageApi}"
 
     if [ "$withPatchNum" = false ]
     then
         echo "Ignoring patch number for output name" 1>&2
-    elif [ "${patch:-0}" != 0 ]
-    then
-        tarName="${tarName}.${patch}"
+    else
+        # Start debian with patch=1, not patch=0
+        if [ "$patch" = 0 ]
+        then
+            patch=1
+        fi
+        tarName="${tarName}${versionSeparator}${patch}"
     fi
     tarName="${tarName}.orig"  # Append .orig
     ;;

etc/colourTables

@@ -44,7 +44,7 @@ coldAndHot
 fire
 {
     // ParaView: Black-Body Radiation
-    interpolate rgb;
+    interpolate rbg;
 
     table
     (
@@ -70,7 +70,7 @@ rainbow
 greyscale
 {
     // ParaView: grayscale
-    interpolate rgb;
+    interpolate rbg;
 
     table
     (
@@ -82,7 +82,7 @@ greyscale
 xray
 {
     // ParaView: "X ray"
-    interpolate rgb;
+    interpolate rbg;
 
     table
     (

etc/controlDict

@@ -36,13 +36,6 @@ InfoSwitches
     // The default ASCII write precision
     writePrecision  6;
 
-    // The output level (verbosity) of messages
-    //  - 0 : suppress all output
-    //  - 1 : normal output
-    //  - >= 2 : report source file and line number (default)
-    outputLevel     2;
-
-    // Report information as baseIOdictionary are read in
     writeDictionaries 0;
 
     // Report optional dictionary entries. For value > 1, treat as fatal.
@@ -678,6 +671,7 @@ DebugSwitches
     lduMesh             0;
     leastSquares        0;
     leastSquaresVectors 0;
+    level               2;
     limitWith           0;
     limited             0;
     limitedCubic        0;

src/OpenFOAM/Make/files

@@ -60,6 +60,8 @@ $(ints)/lists/labelListIOList.C
 primitives/Scalar/doubleScalar/doubleScalar.C
 primitives/Scalar/floatScalar/floatScalar.C
 primitives/Scalar/scalar/scalar.C
+primitives/Scalar/scalar/incGamma.C
+primitives/Scalar/scalar/invIncGamma.C
 primitives/Scalar/lists/scalarList.C
 primitives/Scalar/lists/scalarIOList.C
 primitives/Scalar/lists/scalarListIOList.C
@@ -115,12 +117,6 @@ primitives/functions/Function1/quarterCosineRamp/quarterCosineRamp.C
 primitives/functions/Function1/halfCosineRamp/halfCosineRamp.C
 primitives/functions/Polynomial/polynomialFunction.C
 
-/* Math functions */
-math = primitives/functions/Math
-$(math)/erfInv.C
-$(math)/incGamma.C
-$(math)/invIncGamma.C
-
 primitives/subModelBase/subModelBase.C
 
 strings = primitives/strings

src/OpenFOAM/containers/HashTables/HashPtrTable/HashPtrTable.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2017-2021 OpenCFD Ltd.
+    Copyright (C) 2017-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -26,9 +26,8 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "HashPtrTable.H"
-#include "autoPtr.H"
 #include "error.H"
+#include "HashPtrTable.H"
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -78,28 +77,6 @@ Foam::HashPtrTable<T, Key, Hash>::~HashPtrTable()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-template<class T, class Key, class Hash>
-Foam::autoPtr<T> Foam::HashPtrTable<T, Key, Hash>::release(iterator& iter)
-{
-    if (iter.good())
-    {
-        autoPtr<T> old(iter.val());
-        iter.val() = nullptr;
-        return old;
-    }
-
-    return nullptr;
-}
-
-
-template<class T, class Key, class Hash>
-Foam::autoPtr<T> Foam::HashPtrTable<T, Key, Hash>::release(const Key& key)
-{
-    iterator iter(this->find(key));
-    return this->release(iter);
-}
-
-
 template<class T, class Key, class Hash>
 Foam::autoPtr<T> Foam::HashPtrTable<T, Key, Hash>::remove(iterator& iter)
 {

src/OpenFOAM/containers/HashTables/HashPtrTable/HashPtrTable.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2017-2021 OpenCFD Ltd.
+    Copyright (C) 2017-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -93,7 +93,7 @@ public:
     // Constructors
 
         //- Default construct with default table capacity
-        HashPtrTable() = default;
+        inline HashPtrTable();
 
         //- Construct given initial table capacity
         inline explicit HashPtrTable(const label size);
@@ -130,26 +130,15 @@ public:
 
     // Edit
 
-        //- Release ownership of the pointer and replace with a nullptr
-        //  Includes a safeguard against the end-iterator.
-        //
-        //  \return the entry's old value as an encapsulated pointer.
-        autoPtr<T> release(iterator& iter);
-
-        //- Release ownership of the pointer and replace with a nullptr
-        //
-        //  \return the entry's old value as an encapsulated pointer.
-        autoPtr<T> release(const Key& key);
-
         //- Remove entry specified by given iterator.
         //  Includes a safeguard against the end-iterator.
         //
-        //  \return the entry's old value as an encapsulated pointer.
+        //  \return entry as an encapsulated pointer.
         autoPtr<T> remove(iterator& iter);
 
         //- Remove entry specified by given key.
         //
-        //  \return the entry's old value as an encapsulated pointer.
+        //  \return entry as an encapsulated pointer.
         autoPtr<T> remove(const Key& key);
 
         //- Erase entry specified by given iterator and delete the

src/OpenFOAM/containers/HashTables/HashPtrTable/HashPtrTableI.H

@@ -29,6 +29,13 @@ License
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
+template<class T, class Key, class Hash>
+inline Foam::HashPtrTable<T, Key, Hash>::HashPtrTable()
+:
+    parent_type()
+{}
+
+
 template<class T, class Key, class Hash>
 inline Foam::HashPtrTable<T, Key, Hash>::HashPtrTable(const label size)
 :

src/OpenFOAM/containers/HashTables/HashTable/HashTable.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2017-2021 OpenCFD Ltd.
+    Copyright (C) 2017-2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -661,9 +661,9 @@ void Foam::HashTable<T, Key, Hash>::swap(HashTable<T, Key, Hash>& rhs)
         return;  // Self-swap is a no-op
     }
 
-    std::swap(size_, rhs.size_);
-    std::swap(capacity_, rhs.capacity_);
-    std::swap(table_, rhs.table_);
+    Foam::Swap(size_, rhs.size_);
+    Foam::Swap(capacity_, rhs.capacity_);
+    Foam::Swap(table_, rhs.table_);
 }
 
 

src/OpenFOAM/containers/Lists/FixedList/FixedListI.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2017-2021 OpenCFD Ltd.
+    Copyright (C) 2017-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.

src/OpenFOAM/containers/PtrLists/PtrListOps/PtrListOps.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2019-2021 OpenCFD Ltd.
+    Copyright (C) 2019 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -136,59 +136,8 @@ struct greater
 };
 
 
-//- List of values generated by applying the access operation
-//- to each list item.
-//
-//  For example,
-//  \code
-//  PtrListOps::get(mesh.boundaryMesh(), nameOp<polyPatch>());
-//  \endcode
-template<class ReturnType, class T, class AccessOp>
-List<ReturnType> get
-(
-    const UPtrList<T>& list,
-    const AccessOp& aop
-);
-
-
-//- List of names generated by calling \c name() for each list item
-//- and filtered for matches
-//
-//  For example,
-//  \code
-//  wordRes matches(...);
-//  PtrListOps::names(mesh.boundaryMesh(), matches);
-//
-//  PtrListOps::names(mesh.boundaryMesh(), predicates::always());
-//  \endcode
-template<class T, class UnaryMatchPredicate>
-List<word> names
-(
-    const UPtrList<T>& list,
-    const UnaryMatchPredicate& matcher
-);
-
-
-//- Find first list item with 'name()' that matches, -1 on failure
-template<class T, class UnaryMatchPredicate>
-label firstMatching
-(
-    const UPtrList<T>& list,
-    const UnaryMatchPredicate& matcher
-);
-
-
-//- Extract list indices for all items with 'name()' that matches
-template<class T, class UnaryMatchPredicate>
-labelList findMatching
-(
-    const UPtrList<T>& list,
-    const UnaryMatchPredicate& matcher
-);
-
 } // End namespace ListOps
 
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/OpenFOAM/containers/PtrLists/PtrListOps/PtrListOpsTemplates.C

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2019-2021 OpenCFD Ltd.
+    Copyright (C) 2019-2020 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -103,118 +103,4 @@ void Foam::shuffle(UPtrList<T>& list)
 }
 
 
-
-// Templated implementation for types(), names(), etc - file-scope
-template<class ReturnType, class T, class AccessOp>
-Foam::List<ReturnType> Foam::PtrListOps::get
-(
-    const UPtrList<T>& list,
-    const AccessOp& aop
-)
-{
-    const label len = list.size();
-
-    List<ReturnType> output(len);
-
-    label count = 0;
-    for (label i = 0; i < len; ++i)
-    {
-        const T* ptr = list.get(i);
-
-        if (bool(ptr))
-        {
-            output[count++] = aop(*ptr);
-        }
-    }
-
-    output.resize(count);
-
-    return output;
-}
-
-
-template<class T, class UnaryMatchPredicate>
-Foam::List<Foam::word> Foam::PtrListOps::names
-(
-    const UPtrList<T>& list,
-    const UnaryMatchPredicate& matcher
-)
-{
-    // Possible: const auto aop = nameOp<T>();
-
-    const label len = list.size();
-
-    List<word> output(len);
-
-    label count = 0;
-    for (label i = 0; i < len; ++i)
-    {
-        const T* ptr = list.get(i);
-
-        if (bool(ptr))
-        {
-            if (matcher(ptr->name()))
-            {
-                output[count++] = (ptr->name());
-            }
-        }
-    }
-
-    output.resize(count);
-
-    return output;
-}
-
-
-template<class T, class UnaryMatchPredicate>
-Foam::label Foam::PtrListOps::firstMatching
-(
-    const UPtrList<T>& list,
-    const UnaryMatchPredicate& matcher
-)
-{
-    const label len = list.size();
-
-    for (label i = 0; i < len; ++i)
-    {
-        const T* ptr = list.get(i);
-
-        if (bool(ptr) && matcher(ptr->name()))
-        {
-            return i;
-        }
-    }
-
-    return -1;
-}
-
-
-template<class T, class UnaryMatchPredicate>
-Foam::labelList Foam::PtrListOps::findMatching
-(
-    const UPtrList<T>& list,
-    const UnaryMatchPredicate& matcher
-)
-{
-    const label len = list.size();
-
-    labelList output(len);
-
-    label count = 0;
-    for (label i = 0; i < len; ++i)
-    {
-        const T* ptr = list.get(i);
-
-        if (bool(ptr) && matcher(ptr->name()))
-        {
-            output[count++] = i;
-        }
-    }
-
-    output.resize(count);
-
-    return output;
-}
-
-
 // ************************************************************************* //