TUT: deleting flamePropagationWithObstacles

applications/solvers/combustion/PDRFoam/EaEqn.H

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

applications/solvers/combustion/PDRFoam/EauEqn.H

@@ -16,13 +16,7 @@ if (ign.ignited())
             )*rho/thermo.rhou()
           : -betav*dpdt*rho/thermo.rhou()
         )
-      - 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)
+      - fvm::laplacian(DbThermal, heau)
      ==
         betav*fvOptions(rho, heau)
     );

applications/solvers/combustion/PDRFoam/Make/files

@@ -1,33 +1,33 @@
 XiModels/XiModel/XiModel.C
 XiModels/XiModel/XiModelNew.C
-XiModels/fixed/fixed.C
-XiModels/algebraic/algebraic.C
-XiModels/transport/transport.C
+
+XiModels/transport/transportTwoEqs/transportTwoEqs.C
 
 XiModels/XiEqModels/XiEqModel/XiEqModel.C
 XiModels/XiEqModels/XiEqModel/XiEqModelNew.C
+
 XiModels/XiEqModels/Gulder/Gulder.C
-XiModels/XiEqModels/instabilityXiEq/instabilityXiEq.C
-XiModels/XiEqModels/SCOPEBlendXiEq/SCOPEBlendXiEq.C
-XiModels/XiEqModels/SCOPEXiEq/SCOPEXiEq.C
+XiModels/XiEqModels/instability2XiEq/instability2XiEq.C
+XiModels/XiEqModels/BLMgMaXiEq/BLMgMaXiEq.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/basicXiSubXiEq/basicXiSubXiEq.C
-PDRModels/XiGModels/basicXiSubG/basicXiSubG.C
+PDRModels/XiEqModels/normBasicXiSubXiEq/normBasicXiSubXiEq.C
+PDRModels/XiGModels/normBasicXiSubG/normBasicXiSubG.C
 
 laminarFlameSpeed/SCOPE/SCOPELaminarFlameSpeed.C
 
-/*PDRFoamAutoRefine.C*/
 PDRFoam.C
 
 EXE = $(FOAM_APPBIN)/PDRFoam

applications/solvers/combustion/PDRFoam/PDRFoam.C

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2017 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -81,8 +82,11 @@ 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"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -95,22 +99,24 @@ int main(int argc, char *argv[])
     );
 
     #include "postProcess.H"
-
     #include "addCheckCaseOptions.H"
     #include "setRootCaseLists.H"
     #include "createTime.H"
-    #include "createMesh.H"
-    #include "createControl.H"
+
+    #include "createDynamicFvMesh.H"
+    #include "createDyMControls.H"
+    #include "initContinuityErrs.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;
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -119,33 +125,62 @@ int main(int argc, char *argv[])
 
     while (runTime.run())
     {
-        #include "readTimeControls.H"
+        #include "readDyMControls.H"
+
         #include "compressibleCourantNo.H"
         #include "setDeltaT.H"
 
         ++runTime;
-        Info<< "\n\nTime = " << runTime.timeName() << endl;
 
-        #include "rhoEqn.H"
+        Info<< "\n\nTime = " << runTime.timeName() << endl;
 
         // --- 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"
             }
 
@@ -155,12 +190,14 @@ int main(int argc, char *argv[])
             }
         }
 
+        rho = thermo.rho();
+
         runTime.write();
 
         runTime.printExecutionTime(Info);
     }
 
-    Info<< "\n end\n";
+    Info<< "\nEnd\n" << endl;
 
     return 0;
 }

applications/solvers/combustion/PDRFoam/PDRFoamAutoRefine.C

@@ -1,216 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  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

@@ -1,150 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  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

@@ -1,147 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  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

@@ -0,0 +1,219 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/XiEqModels/normBasicXiSubXiEq/normBasicXiSubXiEq.H

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

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

@@ -1,115 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  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/normBasicXiSubG/normBasicXiSubG.C

@@ -0,0 +1,270 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/XiGModels/normBasicXiSubG/normBasicXiSubG.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2011-2013 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -24,37 +24,18 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::XiGModels::basicSubGrid
+    Foam::XiGModel::normBasicSubGrid
 
 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
-    basicSubGrid.C
+    normBasicSubGrid.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef basicSubGrid_H
-#define basicSubGrid_H
+#ifndef normBasicSubGrid_H
+#define normBasicSubGrid_H
 
 #include "XiGModel.H"
 
@@ -66,43 +47,56 @@ namespace XiGModels
 {
 
 /*---------------------------------------------------------------------------*\
-                        Class basicSubGrid Declaration
+                          Class normBasicSubGrid Declaration
 \*---------------------------------------------------------------------------*/
 
-class basicSubGrid
+class normBasicSubGrid
 :
     public XiGModel
 {
-    // Private data
+    // Private Data
 
         //- Sub-grid generation rate coefficient
-        scalar k1;
+        scalar k1_;
 
-        //- Xi generation rate model due to turbulence
-        autoPtr<XiGModel> XiGModel_;
+        //- 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_;
 
 
     // Private Member Functions
 
         //- No copy construct
-        basicSubGrid(const basicSubGrid&) = delete;
+        normBasicSubGrid(const normBasicSubGrid&) = delete;
 
         //- No copy assignment
-        void operator=(const basicSubGrid&) = delete;
+        void operator=(const normBasicSubGrid&) = delete;
 
 
 public:
 
     //- Runtime type information
-    TypeName("basicSubGridG");
+    TypeName("normBasicSubGridG");
 
 
     // Constructors
 
         //- Construct from components
-        basicSubGrid
+        normBasicSubGrid
         (
             const dictionary& XiGProperties,
+            const word& modelType,
             const psiuReactionThermo& thermo,
             const compressible::RASModel& turbulence,
             const volScalarField& Su
@@ -110,7 +104,7 @@ public:
 
 
     //- Destructor
-    virtual ~basicSubGrid();
+    virtual ~normBasicSubGrid() = default;
 
 
     // Member Functions

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

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

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

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

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

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

@@ -0,0 +1,296 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2011-2012 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -24,90 +24,114 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::XiModels::fixed
+    Foam::PDRDragModels::basicSch
 
 Description
-    Fixed value model for Xi. See \link XiModel.H \endlink for more details
-    on flame wrinkling modelling.
 
 SourceFiles
-    fixed.C
+    basicSch.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef fixed_H
-#define fixed_H
+#ifndef basicSch_H
+#define basicSch_H
 
-#include "XiModel.H"
+#include "PDRDragModel.H"
+#include "XiEqModel.H"
 
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
-namespace XiModels
+namespace PDRDragModels
 {
 
 /*---------------------------------------------------------------------------*\
-                          Class fixed Declaration
+                          Class basicSch Declaration
 \*---------------------------------------------------------------------------*/
 
-class fixed
+class basicSch
 :
-    public XiModel
+    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
-        fixed(const fixed&) = delete;
+        basicSch(const basicSch&) = delete;
 
         //- No copy assignment
-        void operator=(const fixed&) = delete;
+        void operator=(const basicSch&) = delete;
 
 
 public:
 
     //- Runtime type information
-    TypeName("fixed");
+    TypeName("basicSch");
 
 
     // Constructors
 
         //- Construct from components
-        fixed
+        basicSch
         (
-            const dictionary& XiProperties,
-            const psiuReactionThermo& thermo,
+            const dictionary& PDRProperties,
             const compressible::RASModel& turbulence,
-            const volScalarField& Su,
             const volScalarField& rho,
-            const volScalarField& b,
+            const volVectorField& U,
             const surfaceScalarField& phi
         );
 
 
     //- Destructor
-    virtual ~fixed();
+    virtual ~basicSch();
 
 
     // Member Functions
 
-        //- Correct the flame-wrinkling Xi
-        virtual void correct()
-        {}
+        //- 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& XiProperties);
+        virtual bool read(const dictionary& PDRProperties);
 
-        //- Write fields of the XiEq model
-        virtual void writeFields()
-        {}
+        //- 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 XiModels
+} // End namespace PDRDragModels
 } // End namespace Foam
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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

@@ -5,8 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2019-2020 OpenCFD Ltd.
+    Copyright (C) 2011-2015 OpenFOAM Foundation
+    Copyright (C) 2019-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -29,6 +29,7 @@ License
 #include "PDRkEpsilon.H"
 #include "PDRDragModel.H"
 #include "addToRunTimeSelectionTable.H"
+#include "Time.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -70,16 +71,40 @@ PDRkEpsilon::PDRkEpsilon
         modelName
     ),
 
-    C4_
+    C5_(coeffDict_.get<scalar>("C5")),
+    C6_(coeffDict_.get<scalar>("C6")),
+    maxLOverCellW_(coeffDict_.get<scalar>("maxLOverCellW")),
+    lCoef_(coeffDict_.get<scalar>("lCoef")),
+    noTurbUntil_
     (
-        dimensioned<scalar>::getOrAddToDict
+        dimensioned<scalar>::lookupOrAddToDict
         (
-            "C4",
+            "noTurbUntil",
             coeffDict_,
-            0.1
+            0.0
         )
-    )
-{}
+    ),
+    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  * * * * * * * * * * * * * * //
@@ -94,7 +119,9 @@ bool PDRkEpsilon::read()
 {
     if (RASModel::read())
     {
-        C4_.readIfPresent(coeffDict_);
+        noTurbUntil_.readIfPresent(coeffDict_);
+        LOverLobs_.readIfPresent(coeffDict_);
+        LOverMobs_.readIfPresent(coeffDict_);
         return true;
     }
 
@@ -110,16 +137,12 @@ void PDRkEpsilon::correct()
         nut_ = Cmu_*sqr(k_)/epsilon_;
         nut_.correctBoundaryConditions();
 
-        // Re-calculate thermal diffusivity
-        //***HGWalphat_ = mut_/Prt_;
-        //alphat_.correctBoundaryConditions();
-
         return;
     }
 
     RASModel::correct();
 
-    volScalarField divU(fvc::div(phi_/fvc::interpolate(rho_)));
+    volScalarField divU("divU",fvc::div(phi_/fvc::interpolate(rho_)));
 
     if (mesh_.moving())
     {
@@ -127,7 +150,9 @@ 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
@@ -141,59 +166,102 @@ 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 LI
-        (C4_*(Lobs + dimensionedScalar("minLength", dimLength, VSMALL)));
-
-    // Dissipation equation
-    tmp<fvScalarMatrix> epsEqn
+    volScalarField LD
     (
-        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",
+        (LOverLobs_)*(Lobs + dimensionedScalar("minLength", dimLength, VSMALL))
     );
 
-    epsEqn.ref().relax();
+    LD.primitiveFieldRef() = min(LD.primitiveField(), maxLOverCellW_*sqrt(Cw));
 
-    epsEqn.ref().boundaryManipulate(epsilon_.boundaryFieldRef());
+    const volScalarField LI(lCoef_*pow(k_, 3.0/2.0)/epsilon_);
 
-    solve(epsEqn);
-    bound(epsilon_, epsilonMin_);
-
-
-    // Turbulent kinetic energy equation
-
-    tmp<fvScalarMatrix> kEqn
+    IOdictionary combustionProperties
     (
-        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_)
+        IOobject
+        (
+            "combustionProperties",
+            mesh_.time().constant(),
+            mesh_,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        )
     );
 
-    kEqn.ref().relax();
-    solve(kEqn);
-    bound(k_, kMin_);
+    ignition ign(combustionProperties, mesh_.time(), U.mesh());
 
+    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,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -30,30 +31,6 @@ 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
 
@@ -63,6 +40,7 @@ SourceFiles
 #define compressiblePDRkEpsilon_H
 
 #include "kEpsilon.H"
+#include "ignition.H"
 #include "turbulentFluidThermoModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -88,10 +66,18 @@ class PDRkEpsilon
         >
     >
 {
-    // Private data
+    // Private Data
 
         // Model coefficients
-        dimensionedScalar C4_;
+        scalar C5_;
+        scalar C6_;
+        scalar maxLOverCellW_;
+        scalar lCoef_;
+        dimensionedScalar noTurbUntil_;
+        dimensionedScalar LOverLobs_;
+        dimensionedScalar LOverMobs_;
+        scalar bMin_;
+        bool flameFilteredG_;
 
 
 public:

applications/solvers/combustion/PDRFoam/StCourantNo.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2011 OpenFOAM Foundation
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -36,10 +36,9 @@ 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,9 +1,6 @@
-    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

@@ -0,0 +1,258 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/BLMgMaXiEq/BLMgMaXiEq.H

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

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

@@ -5,7 +5,8 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2011-2012 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -45,22 +46,21 @@ namespace XiEqModels
 Foam::XiEqModels::Gulder::Gulder
 (
     const dictionary& XiEqProperties,
+    const word& modelType,
     const psiuReactionThermo& thermo,
     const compressible::RASModel& turbulence,
     const volScalarField& Su
 )
 :
-    XiEqModel(XiEqProperties, thermo, turbulence, Su),
+    XiEqModel(XiEqProperties, modelType, thermo, turbulence, Su),
     XiEqCoef_(XiEqModelCoeffs_.get<scalar>("XiEqCoef")),
     SuMin_(0.01*Su.average()),
     uPrimeCoef_(XiEqModelCoeffs_.get<scalar>("uPrimeCoef")),
-    subGridSchelkin_(XiEqModelCoeffs_.get<bool>("subGridSchelkin"))
-{}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::XiEqModels::Gulder::~Gulder()
+    nrExp_(XiEqModelCoeffs_.get<scalar>("nrExp")),
+    subGridSchelkin_
+    (
+        XiEqModelCoeffs_.get<bool>("subGridSchelkin")
+    )
 {}
 
 
@@ -73,7 +73,7 @@ Foam::tmp<Foam::volScalarField> Foam::XiEqModels::Gulder::XiEq() const
 
     if (subGridSchelkin_)
     {
-        up.primitiveFieldRef() += calculateSchelkinEffect(uPrimeCoef_);
+        up.primitiveFieldRef() += calculateSchelkinEffect(uPrimeCoef_, nrExp_);
     }
 
     volScalarField tauEta(sqrt(mag(thermo_.muu()/(thermo_.rhou()*epsilon))));
@@ -97,6 +97,7 @@ 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,6 +66,8 @@ class Gulder
         //- Schelkin effect Model constant
         scalar uPrimeCoef_;
 
+        scalar nrExp_;
+
         //- Use sub-grid Schelkin effect
         bool subGridSchelkin_;
 
@@ -91,6 +93,7 @@ public:
         Gulder
         (
             const dictionary& XiEqProperties,
+            const word& modelType,
             const psiuReactionThermo& thermo,
             const compressible::RASModel& turbulence,
             const volScalarField& Su
@@ -98,7 +101,7 @@ public:
 
 
     //- Destructor
-    virtual ~Gulder();
+    virtual ~Gulder() = default;
 
 
     // Member Functions
@@ -108,7 +111,6 @@ public:
 
         //- Update properties from given dictionary
         virtual bool read(const dictionary& XiEqProperties);
-
 };
 
 

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

@@ -1,95 +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
--------------------------------------------------------------------------------
-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/SCOPEXiEq/SCOPEXiEq.C

@@ -1,159 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  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/XiEqModel/XiEqModel.C

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

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

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -28,12 +29,6 @@ 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
@@ -59,10 +54,9 @@ namespace Foam
 
 class XiEqModel
 {
-
 protected:
 
-    // Protected data
+    // Protected Data
 
         //- Dictionary
         dictionary XiEqModelCoeffs_;
@@ -77,9 +71,7 @@ protected:
         const volScalarField& Su_;
 
 
-private:
-
-    // Private Member Functions
+    // Protected Member Functions
 
         //- No copy construct
         XiEqModel(const XiEqModel&) = delete;
@@ -103,12 +95,14 @@ public:
             dictionary,
             (
                 const dictionary& XiEqProperties,
+                const word& modelType,
                 const psiuReactionThermo& thermo,
                 const compressible::RASModel& turbulence,
                 const volScalarField& Su
             ),
             (
                 XiEqProperties,
+                modelType,
                 thermo,
                 turbulence,
                 Su
@@ -122,6 +116,7 @@ public:
         static autoPtr<XiEqModel> New
         (
             const dictionary& XiEqProperties,
+            const word& modelType,
             const psiuReactionThermo& thermo,
             const compressible::RASModel& turbulence,
             const volScalarField& Su
@@ -134,6 +129,7 @@ public:
         XiEqModel
         (
             const dictionary& XiEqProperties,
+            const word& modelType,
             const psiuReactionThermo& thermo,
             const compressible::RASModel& turbulence,
             const volScalarField& Su
@@ -147,13 +143,14 @@ public:
     // Member Functions
 
         //- Return the flame-wrinkling XiEq
-        virtual tmp<volScalarField> XiEq() const
-        {
-            return turbulence_.muEff();
-        }
+        virtual tmp<volScalarField> XiEq() const = 0;
 
-        //- Return the sub-grid Schelkin effect
-        tmp<volScalarField> calculateSchelkinEffect(const scalar) const;
+        // Only used by sperseded modles
+        tmp<volScalarField> calculateSchelkinEffect
+        (
+            const scalar uPrimeCoef,
+            const scalar nrExp
+        ) 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 OpenCFD Ltd.
+    Copyright (C) 2019-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -33,14 +33,15 @@ 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>("XiEqModel"));
+    const word modelType(dict.get<word>(modelName));
 
-    Info<< "Selecting flame-wrinkling model " << modelType << endl;
+    Info<< "Selecting flame-wrinkling model Eq " << modelType << endl;
 
     auto cstrIter = dictionaryConstructorTablePtr_->cfind(modelType);
 
@@ -49,13 +50,16 @@ Foam::autoPtr<Foam::XiEqModel> Foam::XiEqModel::New
         FatalIOErrorInLookup
         (
             dict,
-            "XiEqModel",
+            modelName,
             modelType,
             *dictionaryConstructorTablePtr_
         ) << exit(FatalIOError);
     }
 
-    return autoPtr<XiEqModel>(cstrIter()(dict, thermo, turbulence, Su));
+    return autoPtr<XiEqModel>(cstrIter()
+    (
+        dict, modelName, thermo, turbulence, Su)
+    );
 }
 
 

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

@@ -0,0 +1,133 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/instability2XiEq/instability2XiEq.H

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2011-2012 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -24,22 +24,25 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::XiEqModels::SCOPEBlend
+    Foam::XiEqModels::instability2XiEq
 
 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
-    SCOPEBlend.C
+    instability2XiEq.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef SCOPEBlend_H
-#define SCOPEBlend_H
+#ifndef instability2XiEq_H
+#define instability2XiEq_H
 
+#include "laminarFlameSpeed.H"
+#include "SCOPELaminarFlameSpeed.H"
+#include "ignitionSite.H"
+#include "ignition.H"
+#include "Time.H"
+#include "fvMesh.H"
 #include "XiEqModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -50,43 +53,53 @@ namespace XiEqModels
 {
 
 /*---------------------------------------------------------------------------*\
-                          Class SCOPEBlend Declaration
+                          Class instability2 Declaration
 \*---------------------------------------------------------------------------*/
 
-class SCOPEBlend
+class instability2XiEq
 :
     public XiEqModel
 {
     // Private data
 
-        //- Low turbulence intensity equilibrium Xi model
-        autoPtr<XiEqModel> XiEqModelL_;
+        //- Laminar burning speed
+        laminarFlameSpeedModels::SCOPE saModel_;
 
-        //- High turbulence intensity equilibrium Xi model
-        autoPtr<XiEqModel> XiEqModelH_;
+        //- 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_;
 
 
     // Private Member Functions
 
         //- No copy construct
-        SCOPEBlend(const SCOPEBlend&) = delete;
+        instability2XiEq(const instability2XiEq&) = delete;
 
         //- No copy assignment
-        void operator=(const SCOPEBlend&) = delete;
+        void operator=(const instability2XiEq&) = delete;
 
 
 public:
 
     //- Runtime type information
-    TypeName("SCOPEBlend");
+    TypeName("instability2XiEq");
 
 
     // Constructors
 
         //- Construct from components
-        SCOPEBlend
+        instability2XiEq
         (
             const dictionary& XiEqProperties,
+            const word& modelType,
             const psiuReactionThermo& thermo,
             const compressible::RASModel& turbulence,
             const volScalarField& Su
@@ -94,7 +107,7 @@ public:
 
 
     //- Destructor
-    virtual ~SCOPEBlend();
+    virtual ~instability2XiEq() = default;
 
 
     // Member Functions
@@ -103,11 +116,7 @@ public:
         virtual tmp<volScalarField> XiEq() const;
 
         //- Update properties from given dictionary
-        virtual bool read(const dictionary& XiEqProperties)
-        {
-            return true;
-        }
-
+        virtual bool read(const dictionary& XiEqProperties);
 };
 
 

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

@@ -1,84 +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
--------------------------------------------------------------------------------
-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/XiGModels/KTS/KTS.C

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -45,22 +46,17 @@ namespace XiGModels
 Foam::XiGModels::KTS::KTS
 (
     const dictionary& XiGProperties,
+    const word& modelType,
     const psiuReactionThermo& thermo,
     const compressible::RASModel& turbulence,
     const volScalarField& Su
 )
 :
-    XiGModel(XiGProperties, thermo, turbulence, Su),
+    XiGModel(XiGProperties, modelType, thermo, turbulence, Su),
     GEtaCoef_(XiGModelCoeffs_.get<scalar>("GEtaCoef"))
 {}
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-Foam::XiGModels::KTS::~KTS()
-{}
-
-
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
 Foam::tmp<Foam::volScalarField> Foam::XiGModels::KTS::G() const
@@ -74,6 +70,15 @@ 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,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2012 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -27,8 +28,8 @@ Class
     Foam::XiGModels::KTS
 
 Description
-    Simple Kolmogorov time-scale (KTS) model for the flame-wrinling generation
-    rate.
+    Simple Kolmogorov time-scale (KTS) model for the flame-wrinkling
+    generation rate.
 
 SourceFiles
     KTS.C
@@ -55,7 +56,7 @@ class KTS
 :
     public XiGModel
 {
-    // Private data
+    // Private Data
 
         scalar GEtaCoef_;
 
@@ -81,6 +82,7 @@ public:
         KTS
         (
             const dictionary& XiGProperties,
+            const word& modelType,
             const psiuReactionThermo& thermo,
             const compressible::RASModel& turbulence,
             const volScalarField& Su
@@ -88,7 +90,7 @@ public:
 
 
     //- Destructor
-    virtual ~KTS();
+    virtual ~KTS() = default;
 
 
     // Member Functions
@@ -96,6 +98,9 @@ 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-2017 OpenFOAM Foundation
+    Copyright (C) 2011-2015 OpenFOAM Foundation
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -41,6 +41,7 @@ namespace Foam
 Foam::XiGModel::XiGModel
 (
     const dictionary& XiGProperties,
+    const word& modelType,
     const psiuReactionThermo& thermo,
     const compressible::RASModel& turbulence,
     const volScalarField& Su
@@ -50,7 +51,7 @@ Foam::XiGModel::XiGModel
     (
         XiGProperties.subDict
         (
-            XiGProperties.get<word>("XiGModel") + "Coeffs"
+            XiGProperties.get<word>(modelType) + "Coeffs"
         )
     ),
     thermo_(thermo),
@@ -69,7 +70,7 @@ Foam::XiGModel::~XiGModel()
 
 bool Foam::XiGModel::read(const dictionary& XiGProperties)
 {
-    XiGModelCoeffs_ = XiGProperties.optionalSubDict(type() + "Coeffs");
+    XiGModelCoeffs_ = XiGProperties.subDict(type() + "Coeffs");
 
     return true;
 }

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

@@ -28,9 +28,6 @@ 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
@@ -56,10 +53,9 @@ namespace Foam
 
 class XiGModel
 {
-
 protected:
 
-    // Protected data
+    // Protected Data
 
         dictionary XiGModelCoeffs_;
 
@@ -68,9 +64,7 @@ protected:
         const volScalarField& Su_;
 
 
-private:
-
-    // Private Member Functions
+    // Protected Member Functions
 
         //- No copy construct
         XiGModel(const XiGModel&) = delete;
@@ -94,12 +88,14 @@ public:
             dictionary,
             (
                 const dictionary& XiGProperties,
+                const word& modelType,
                 const psiuReactionThermo& thermo,
                 const compressible::RASModel& turbulence,
                 const volScalarField& Su
             ),
             (
                 XiGProperties,
+                modelType,
                 thermo,
                 turbulence,
                 Su
@@ -113,6 +109,7 @@ public:
         static autoPtr<XiGModel> New
         (
             const dictionary& XiGProperties,
+            const word& modelType,
             const psiuReactionThermo& thermo,
             const compressible::RASModel& turbulence,
             const volScalarField& Su
@@ -125,6 +122,7 @@ public:
         XiGModel
         (
             const dictionary& XiGProperties,
+            const word& modelType,
             const psiuReactionThermo& thermo,
             const compressible::RASModel& turbulence,
             const volScalarField& Su
@@ -137,14 +135,11 @@ public:
 
     // Member Functions
 
-        //- Return the flame-wrinkling generation rate
+        //- Return the flame-wrinkling genration rate
         virtual tmp<volScalarField> G() const = 0;
 
         //- Return the flame diffusivity
-        virtual tmp<volScalarField> Db() const
-        {
-            return turbulence_.muEff();
-        }
+        virtual tmp<volScalarField> Db() const = 0;
 
         //- 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 OpenCFD Ltd.
+    Copyright (C) 2019-2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -33,14 +33,15 @@ 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>("XiGModel"));
+    const word modelType(dict.get<word>(modelName));
 
-    Info<< "Selecting flame-wrinkling model " << modelType << endl;
+    Info<< "Selecting flame-wrinkling model for G " << modelType << endl;
 
     auto cstrIter = dictionaryConstructorTablePtr_->cfind(modelType);
 
@@ -49,13 +50,16 @@ Foam::autoPtr<Foam::XiGModel> Foam::XiGModel::New
         FatalIOErrorInLookup
         (
             dict,
-            "XiGModel",
+            modelName,
             modelType,
             *dictionaryConstructorTablePtr_
         ) << exit(FatalIOError);
     }
 
-    return autoPtr<XiGModel>(cstrIter()(dict, thermo, turbulence, Su));
+    return autoPtr<XiGModel>(cstrIter()
+    (
+        dict, modelName, thermo, turbulence, Su)
+    );
 }
 
 

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

@@ -0,0 +1,165 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2011-2012 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -24,106 +24,115 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::XiModels::algebraic
+    Foam::XiGModels::instability2G
 
 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.
+    Flame-surface instabilityG flame-wrinkling generation rate coefficient model
+    used in \link XiModel.H \endlink.
 
 SourceFiles
-    algebraic.C
+    instability2G.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef algebraic_H
-#define algebraic_H
+#ifndef instability2G_H
+#define instability2G_H
 
-#include "XiModel.H"
-#include "XiEqModel.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 XiModels
+namespace XiGModels
 {
 
 /*---------------------------------------------------------------------------*\
-                          Class algebraic Declaration
+                          Class instability2G Declaration
 \*---------------------------------------------------------------------------*/
 
-class algebraic
+class instability2G
 :
-    public XiModel
+    public XiGModel
 {
-    // Private data
+    // Private Data
 
-        scalar XiShapeCoef;
+        laminarFlameSpeedModels::SCOPE saModel_;
 
-        autoPtr<XiEqModel> XiEqModel_;
+        // 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
-        algebraic(const algebraic&) = delete;
+        instability2G(const instability2G&) = delete;
 
         //- No copy assignment
-        void operator=(const algebraic&) = delete;
+        void operator=(const instability2G&) = delete;
 
 
 public:
 
     //- Runtime type information
-    TypeName("algebraic");
+    TypeName("instability2G");
 
 
     // Constructors
 
         //- Construct from components
-        algebraic
+        instability2G
         (
-            const dictionary& XiProperties,
+            const dictionary& XiGProperties,
+            const word& modelType,
             const psiuReactionThermo& thermo,
             const compressible::RASModel& turbulence,
-            const volScalarField& Su,
-            const volScalarField& rho,
-            const volScalarField& b,
-            const surfaceScalarField& phi
+            const volScalarField& Su
         );
 
 
     //- Destructor
-    virtual ~algebraic();
+    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;
 
-        //- 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();
-        }
-
+        virtual bool read(const dictionary& XiGProperties);
 };
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-} // End namespace XiModels
+} // End namespace XiGModels
 } // End namespace Foam
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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

@@ -5,7 +5,7 @@
     \\  /    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.
@@ -45,21 +45,16 @@ namespace XiGModels
 Foam::XiGModels::instabilityG::instabilityG
 (
     const dictionary& XiGProperties,
+    const word& modelType,
     const psiuReactionThermo& thermo,
     const compressible::RASModel& turbulence,
     const volScalarField& 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()
+    XiGModel(XiGProperties, modelType, thermo, turbulence, Su),
+    GIn_("GIn", XiGModelCoeffs_),
+    lambdaIn_("lambdaIn", XiGModelCoeffs_),
+    XiGModel_(XiGModel::New(XiGModelCoeffs_,modelType,thermo, turbulence, Su))
 {}
 
 
@@ -76,11 +71,14 @@ 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");
 
-    return XiGModel_->Db()
-        + rho*Su_*(Xi - 1.0)*mgb*(0.5*lambdaIn_)/(mgb + 1.0/lambdaIn_);
+    //// OLD return turbulence_.muEff()
+    return Db1
+        + rho*Su_*(Xp*Xi - 1.0)*mgb*(0.5*lambdaIn_)/(mgb + 1.0/lambdaIn_);
 }
 
 

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

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2012 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -30,8 +31,6 @@ 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
 
@@ -57,9 +56,9 @@ class instabilityG
 :
     public XiGModel
 {
-    // Private data
+    // Private Data
 
-        //- Flame instabilityG wrinkling generation rate coefficient
+        //- Flame instabilityG wrinling generation rate coefficient
         dimensionedScalar GIn_;
 
         //- InstabilityG length-scale
@@ -90,6 +89,7 @@ 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();
+    virtual ~instabilityG() = default;
 
 
     // Member Functions

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

@@ -5,7 +5,7 @@
     \\  /    A nd           | www.openfoam.com
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
-    Copyright (C) 2011-2017 OpenFOAM Foundation
+    Copyright (C) 2011-2015 OpenFOAM Foundation
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -73,7 +73,19 @@ 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 "multivariateSurfaceInterpolationScheme.H"
 #include "fvcDiv.H"
+#include "multivariateSurfaceInterpolationScheme.H"
 #include "runTimeSelectionTables.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -110,10 +110,9 @@ namespace Foam
 
 class XiModel
 {
-
 protected:
 
-    // Protected data
+    // Protected Data
 
         dictionary XiModelCoeffs_;
 
@@ -127,10 +126,11 @@ protected:
         //- Flame wrinkling field
         volScalarField Xi_;
 
+        //- Flame wrinkling field for the subgrid obstacles
+        volScalarField Xp_;
 
-private:
 
-    // Private Member Functions
+    // Protected Member Functions
 
         //- No copy construct
         XiModel(const XiModel&) = delete;
@@ -215,19 +215,31 @@ 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
         {
-            return turbulence_.muEff();
+            const objectRegistry& db = Su_.db();
+            const volScalarField& Db1 = db.lookupObject<volScalarField>("Db");
+            //return turbulence_.muEff();
+            return Db1;
         }
 
         //- Add Xi to the multivariateSurfaceInterpolationScheme table
         //  if required
-        virtual void addXi
+        void addXi
         (
-            multivariateSurfaceInterpolationScheme<scalar>::fieldTable&
+            multivariateSurfaceInterpolationScheme<scalar>::fieldTable& fields
         )
-        {}
+        {
+            fields.add(Xi_);
+            fields.add(Xp_);
+        }
 
         //- Correct the flame-wrinkling Xi
         virtual void correct() = 0;

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

@@ -1,100 +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
--------------------------------------------------------------------------------
-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/transport/transport.C

@@ -1,147 +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
--------------------------------------------------------------------------------
-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/transportTwoEqs/transportTwoEqs.C

@@ -0,0 +1,258 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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/XiModels/transport/transportTwoEqs/transportTwoEqs.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -24,21 +25,17 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::XiModels::transport
+    Foam::XiModels::transportTwoEqs
 
 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
-    transport.C
+    transportTwoEqs.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef transport_H
-#define transport_H
+#ifndef transportTwoEqs_H
+#define transportTwoEqs_H
 
 #include "XiModel.H"
 #include "XiEqModel.H"
@@ -52,40 +49,58 @@ namespace XiModels
 {
 
 /*---------------------------------------------------------------------------*\
-                          Class transport Declaration
+                          Class transportTwoEqs Declaration
 \*---------------------------------------------------------------------------*/
 
-class transport
+class transportTwoEqs
 :
     public XiModel
 {
     // Private data
 
-        scalar XiShapeCoef;
+        scalar XiShapeCoef_;
+        scalar CpfiDot_;
+        scalar CpfiCross_;
+        scalar GEtaExp_;
 
+        //- 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
-        transport(const transport&) = delete;
+        transportTwoEqs(const transportTwoEqs&) = delete;
 
         //- No copy assignment
-        void operator=(const transport&) = delete;
+        void operator=(const transportTwoEqs&) = delete;
 
 
 public:
 
     //- Runtime type information
-    TypeName("transport");
+    TypeName("transportTwoEqs");
 
 
     // Constructors
 
         //- Construct from components
-        transport
+        transportTwoEqs
         (
             const dictionary& XiProperties,
             const psiuReactionThermo& thermo,
@@ -98,7 +113,7 @@ public:
 
 
     //- Destructor
-    virtual ~transport();
+    virtual ~transportTwoEqs();
 
 
     // Member Functions
@@ -106,15 +121,6 @@ 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()
         {
@@ -131,6 +137,7 @@ public:
         virtual void writeFields()
         {
             XiEqModel_().writeFields();
+            XpEqModel_().writeFields();
         }
 };
 

applications/solvers/combustion/PDRFoam/bEqn.H

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

applications/solvers/combustion/PDRFoam/createFieldRefs.H

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

applications/solvers/combustion/PDRFoam/createFields.H

@@ -1,218 +1,303 @@
-Info<< "Reading thermophysical properties\n" << endl;
+    Info<< "Reading thermophysical properties\n" << endl;
 
-autoPtr<psiuReactionThermo> pThermo
-(
-    psiuReactionThermo::New(mesh)
-);
-psiuReactionThermo& thermo = pThermo();
-thermo.validate(args.executable(), "ha", "ea");
-
-basicSpecieMixture& composition = thermo.composition();
-
-volScalarField rho
-(
-    IOobject
+    autoPtr<psiuReactionThermo> pThermo
     (
-        "rho",
-        runTime.timeName(),
-        mesh,
-        IOobject::NO_READ,
-        IOobject::AUTO_WRITE
-    ),
-    thermo.rho()
-);
+        psiuReactionThermo::New(mesh)
+    );
+    psiuReactionThermo& thermo = pThermo();
+    thermo.validate(args.executable(), "ha", "ea");
+    basicMultiComponentMixture& composition = thermo.composition();
 
-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
+    volScalarField rho
     (
-        "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
+        (
+            "rho",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        thermo.rho()
+    );
+    volScalarField rhou
     (
+        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,
-        thermo
-    )
-);
+        phi
+    );
 
-
-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
+    //- Create the flame-wrinkling model
+    autoPtr<XiModel> flameWrinkling = XiModel::New
     (
-        "Su",
-        runTime.timeName(),
-        mesh,
-        IOobject::MUST_READ,
-        IOobject::AUTO_WRITE
-    ),
-    mesh
-);
+        PDRProperties,
+        thermo,
+        turbulence,
+        Su,
+        rho,
+        b,
+        phi
+    );
 
-Info<< "Reading field betav\n" << endl;
-volScalarField betav
-(
-    IOobject
+    Info<< "Calculating turbulent flame speed field St\n" << endl;
+    volScalarField St
     (
-        "betav",
-        mesh.facesInstance(),
-        mesh,
-        IOobject::MUST_READ,
-        IOobject::NO_WRITE
-    ),
-    mesh
-);
+        IOobject
+        (
+            "St",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        flameWrinkling->Xi()*Su
+    );
 
-Info<< "Reading field Lobs\n" << endl;
-volScalarField Lobs
-(
-    IOobject
+    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
     (
-        "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"
+        IOobject
+        (
+            "c",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        scalar(1) - b
+    );

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

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

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

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2012 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -65,6 +66,7 @@ SourceFiles
 #define SCOPE_H
 
 #include "laminarFlameSpeed.H"
+#include "Polynomial.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -105,7 +107,7 @@ class SCOPE
             scalar lu;
 
             //- Construct from dictionary
-            polynomial(const dictionary& polyDict);
+            explicit polynomial(const dictionary& polyDict);
         };
 
 
@@ -123,11 +125,14 @@ class SCOPE
         //- Upper Su polynomial
         polynomial SuPolyU_;
 
-        //- Temperature correction exponent
-        scalar Texp_;
+        //- Temperature correction exponents
+        Polynomial<7> Texp_;
 
-        //- Pressure correction exponent
-        scalar pexp_;
+        //- Pressure correction exponents
+        Polynomial<7> pexp_;
+
+        //- Additional coefficient
+        scalar CIn_;
 
         //- Lower Ma polynomial
         polynomial MaPolyL_;
@@ -136,25 +141,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,
@@ -163,7 +168,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,
@@ -172,13 +177,14 @@ 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;
 
-        //- Construct as copy (not implemented)
-        SCOPE(const SCOPE&);
+        //- No copy construct
+        SCOPE(const SCOPE&) = delete;
 
-        void operator=(const SCOPE&);
+        //- No copy assignment
+        void operator=(const SCOPE&) = delete;
 
 
 public:
@@ -191,8 +197,8 @@ public:
         //- Construct from dictionary and psiuReactionThermo
         SCOPE
         (
-            const dictionary&,
-            const psiuReactionThermo&
+            const dictionary& dict,
+            const psiuReactionThermo& ct
         );
 
 
@@ -202,6 +208,9 @@ public:
 
     // Member functions
 
+        //- Return CIn
+        scalar CIn() const noexcept;
+
         //- Return the Markstein number
         tmp<volScalarField> Ma() const;
 

applications/solvers/combustion/PDRFoam/pEqn.H

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

applications/solvers/combustion/PDRFoam/readTimeControls.H

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -23,52 +24,26 @@ 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.
+
 \*---------------------------------------------------------------------------*/
 
-#include "fixed.H"
-#include "addToRunTimeSelectionTable.H"
+const bool adjustTimeStep =
+    runTime.controlDict().getOrDefault("adjustTimeStep", false);
 
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+scalar maxCo =
+    runTime.controlDict().getOrDefault<scalar>("maxCo", 1);
 
-namespace Foam
-{
-namespace XiModels
-{
-    defineTypeNameAndDebug(fixed, 0);
-    addToRunTimeSelectionTable(XiModel, fixed, dictionary);
-}
-}
+scalar maxDeltaT =
+    runTime.controlDict().getOrDefault<scalar>("maxDeltaT", GREAT);
 
-
-// * * * * * * * * * * * * * * * * 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);
-}
+scalar maxDeltaTRatio =
+    runTime.controlDict().getOrDefault<scalar>("maxDeltaTRatio", 1.2);
 
 
 // ************************************************************************* //

applications/solvers/combustion/PDRFoam/rhoEqn.H

@@ -1,38 +1,3 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  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,7 +36,8 @@ Description
 if (adjustTimeStep)
 {
     scalar maxDeltaTFact = maxCo/(CoNum + StCoNum + SMALL);
-    scalar deltaTFact = min(min(maxDeltaTFact, 1.0 + 0.1*maxDeltaTFact), 1.2);
+    scalar deltaTFact =
+        min(min(maxDeltaTFact, 1.0 + 0.1*maxDeltaTFact), maxDeltaTRatio);
 
     runTime.setDeltaT
     (

tutorials/combustion/PDRFoam/flamePropagationWithObstacles/0.orig/Su

@@ -1,53 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  v2012                                 |
-|   \\  /    A nd           | Website:  www.openfoam.com                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       volScalarField;
-    location    "0";
-    object      Su;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-dimensions      [0 1 -1 0 0 0 0];
-
-internalField   uniform 0.5;
-
-boundaryField
-{
-    outer
-    {
-        type                uniformInletOutlet;
-        uniformInletValue   0.5;
-        value               uniform 0.5;
-    }
-    ground
-    {
-        type            zeroGradient;
-    }
-    blockedFaces
-    {
-        type            zeroGradient;
-    }
-    baffleWall
-    {
-        type            zeroGradient;
-    }
-    baffleCyclic_half0
-    {
-        type            cyclic;
-    }
-    baffleCyclic_half1
-    {
-        type            cyclic;
-    }
-}
-
-
-// ************************************************************************* //

tutorials/combustion/PDRFoam/flamePropagationWithObstacles/0.orig/T

@@ -1,55 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  v2012                                 |
-|   \\  /    A nd           | Website:  www.openfoam.com                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       volScalarField;
-    location    "0";
-    object      T;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-dimensions      [0 0 0 1 0 0 0];
-
-internalField   uniform 300;
-
-boundaryField
-{
-    outer
-    {
-        type                uniformInletOutlet;
-        uniformInletValue   300;
-        value               uniform 300;
-    }
-    ground
-    {
-        type            zeroGradient;
-    }
-    blockedFaces
-    {
-        type            zeroGradient;
-        value           uniform 300;
-    }
-    baffleWall
-    {
-        type            zeroGradient;
-        value           uniform 300;
-    }
-    baffleCyclic_half0
-    {
-        type            cyclic;
-    }
-    baffleCyclic_half1
-    {
-        type            cyclic;
-    }
-}
-
-
-// ************************************************************************* //

tutorials/combustion/PDRFoam/flamePropagationWithObstacles/0.orig/Tu

@@ -1,55 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  v2012                                 |
-|   \\  /    A nd           | Website:  www.openfoam.com                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       volScalarField;
-    location    "0";
-    object      Tu;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-dimensions      [0 0 0 1 0 0 0];
-
-internalField   uniform 300;
-
-boundaryField
-{
-    outer
-    {
-        type                uniformInletOutlet;
-        uniformInletValue   300;
-        value               uniform 300;
-    }
-    ground
-    {
-        type            zeroGradient;
-    }
-    blockedFaces
-    {
-        type            zeroGradient;
-        value           uniform 300;
-    }
-    baffleWall
-    {
-        type            zeroGradient;
-        value           uniform 300;
-    }
-    baffleCyclic_half0
-    {
-        type            cyclic;
-    }
-    baffleCyclic_half1
-    {
-        type            cyclic;
-    }
-}
-
-
-// ************************************************************************* //

tutorials/combustion/PDRFoam/flamePropagationWithObstacles/0.orig/U

@@ -1,61 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  v2012                                 |
-|   \\  /    A nd           | Website:  www.openfoam.com                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       volVectorField;
-    location    "0";
-    object      U;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-dimensions      [0 1 -1 0 0 0 0];
-
-internalField   uniform (0 0 0);
-
-boundaryField
-{
-    outer
-    {
-        type                uniformInletOutlet;
-        uniformInletValue   (0 0 0);
-        value               uniform (0 0 0);
-    }
-    ground
-    {
-        type            noSlip;
-    }
-    blockedFaces
-    {
-        type            noSlip;
-    }
-    baffleWall
-    {
-        type            activePressureForceBaffleVelocity;
-        value           uniform (0 0 0);
-        cyclicPatch     baffleCyclic_half0;
-        opening         1;
-        openingTime     0.01;
-        maxOpenFractionDelta 0.1;
-        openFraction    0;
-        minThresholdValue 8000;
-        forceBased      0;
-    }
-    baffleCyclic_half0
-    {
-        type            cyclic;
-    }
-    baffleCyclic_half1
-    {
-        type            cyclic;
-    }
-}
-
-
-// ************************************************************************* //

tutorials/combustion/PDRFoam/flamePropagationWithObstacles/0.orig/Xi

@@ -1,55 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  v2012                                 |
-|   \\  /    A nd           | Website:  www.openfoam.com                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       volScalarField;
-    location    "0";
-    object      Xi;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-dimensions      [0 0 0 0 0 0 0];
-
-internalField   uniform 1;
-
-boundaryField
-{
-    outer
-    {
-        type                uniformInletOutlet;
-        uniformInletValue   1;
-        value               uniform 1;
-    }
-    ground
-    {
-        type            zeroGradient;
-    }
-    blockedFaces
-    {
-        type            zeroGradient;
-        value           uniform 1;
-    }
-    baffleWall
-    {
-        type            zeroGradient;
-        value           uniform 1;
-    }
-    baffleCyclic_half0
-    {
-        type            cyclic;
-    }
-    baffleCyclic_half1
-    {
-        type            cyclic;
-    }
-}
-
-
-// ************************************************************************* //

tutorials/combustion/PDRFoam/flamePropagationWithObstacles/0.orig/alphat

@@ -1,58 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  v2012                                 |
-|   \\  /    A nd           | Website:  www.openfoam.com                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       volScalarField;
-    location    "0";
-    object      alphat;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-dimensions      [1 -1 -1 0 0 0 0];
-
-internalField   uniform 0;
-
-boundaryField
-{
-    outer
-    {
-        type            calculated;
-        value           uniform 0;
-    }
-    ground
-    {
-        type            compressible::alphatWallFunction;
-        Prt             0.85;
-        value           uniform 0;
-    }
-    blockedFaces
-    {
-        type            compressible::alphatWallFunction;
-        Prt             0.85;
-        value           nonuniform 0();
-    }
-    baffleWall
-    {
-        type            compressible::alphatWallFunction;
-        Prt             0.85;
-        value           nonuniform 0();
-    }
-    baffleCyclic_half0
-    {
-        type            cyclic;
-    }
-    baffleCyclic_half1
-    {
-        type            cyclic;
-    }
-}
-
-
-// ************************************************************************* //

tutorials/combustion/PDRFoam/flamePropagationWithObstacles/0.orig/b

@@ -1,55 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  v2012                                 |
-|   \\  /    A nd           | Website:  www.openfoam.com                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       volScalarField;
-    location    "0";
-    object      b;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-dimensions      [0 0 0 0 0 0 0];
-
-internalField   uniform 1;
-
-boundaryField
-{
-    outer
-    {
-        type                uniformInletOutlet;
-        uniformInletValue   1;
-        value               uniform 1;
-    }
-    ground
-    {
-        type            zeroGradient;
-    }
-    blockedFaces
-    {
-        type            zeroGradient;
-        value           uniform 1;
-    }
-    baffleWall
-    {
-        type            zeroGradient;
-        value           uniform 1;
-    }
-    baffleCyclic_half0
-    {
-        type            cyclic;
-    }
-    baffleCyclic_half1
-    {
-        type            cyclic;
-    }
-}
-
-
-// ************************************************************************* //

tutorials/combustion/PDRFoam/flamePropagationWithObstacles/0.orig/epsilon

@@ -1,65 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  v2012                                 |
-|   \\  /    A nd           | Website:  www.openfoam.com                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       volScalarField;
-    location    "0";
-    object      epsilon;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-dimensions      [0 2 -3 0 0 0 0];
-
-internalField   uniform 0.1;
-
-boundaryField
-{
-    outer
-    {
-        type                uniformInletOutlet;
-        uniformInletValue   0.1;
-        value               uniform 0.1;
-    }
-    ground
-    {
-        type            epsilonWallFunction;
-        Cmu             0.09;
-        kappa           0.41;
-        E               9.8;
-        value           uniform 0.1;
-    }
-    blockedFaces
-    {
-        type            epsilonWallFunction;
-        Cmu             0.09;
-        kappa           0.41;
-        E               9.8;
-        value           uniform 0.1;
-    }
-    baffleWall
-    {
-        type            epsilonWallFunction;
-        Cmu             0.09;
-        kappa           0.41;
-        E               9.8;
-        value           uniform 1e-05;
-    }
-    baffleCyclic_half0
-    {
-        type            cyclic;
-    }
-    baffleCyclic_half1
-    {
-        type            cyclic;
-    }
-}
-
-
-// ************************************************************************* //

tutorials/combustion/PDRFoam/flamePropagationWithObstacles/0.orig/ft

@@ -1,54 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  v2012                                 |
-|   \\  /    A nd           | Website:  www.openfoam.com                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       volScalarField;
-    location    "0";
-    object      ft;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-dimensions      [0 0 0 0 0 0 0];
-
-internalField   uniform 0.0623;
-
-boundaryField
-{
-    outer
-    {
-        type                uniformInletOutlet;
-        uniformInletValue   0.0623;
-        value               uniform 0.0623;
-    }
-    ground
-    {
-        type            zeroGradient;
-    }
-    blockedFaces
-    {
-        type            zeroGradient;
-    }
-    baffleWall
-    {
-        type            zeroGradient;
-        value           uniform 0.06024096;
-    }
-    baffleCyclic_half0
-    {
-        type            cyclic;
-    }
-    baffleCyclic_half1
-    {
-        type            cyclic;
-    }
-}
-
-
-// ************************************************************************* //

tutorials/combustion/PDRFoam/flamePropagationWithObstacles/0.orig/k

@@ -1,56 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  v2012                                 |
-|   \\  /    A nd           | Website:  www.openfoam.com                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       volScalarField;
-    location    "0";
-    object      k;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-dimensions      [0 2 -2 0 0 0 0];
-
-internalField   uniform 1.5;
-
-boundaryField
-{
-    outer
-    {
-        type                uniformInletOutlet;
-        uniformInletValue   1.5;
-        value               uniform 1.5;
-    }
-    ground
-    {
-        type            kqRWallFunction;
-        value           uniform 1.5;
-    }
-    blockedFaces
-    {
-        type            kqRWallFunction;
-        value           uniform 1.5;
-    }
-    baffleWall
-    {
-        type            kqRWallFunction;
-        value           uniform 1.5;
-    }
-    baffleCyclic_half0
-    {
-        type            cyclic;
-    }
-    baffleCyclic_half1
-    {
-        type            cyclic;
-    }
-}
-
-
-// ************************************************************************* //

tutorials/combustion/PDRFoam/flamePropagationWithObstacles/0.orig/k.old

@@ -1,56 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  v2012                                 |
-|   \\  /    A nd           | Website:  www.openfoam.com                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       volScalarField;
-    location    "0";
-    object      k;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-dimensions      [0 2 -2 0 0 0 0];
-
-internalField   uniform 1.5;
-
-boundaryField
-{
-    outer
-    {
-        type                uniformInletOutlet;
-        uniformInletValue   1.5;
-        value               uniform 1.5;
-    }
-    ground
-    {
-        type            kqRWallFunction;
-        value           uniform 1.5;
-    }
-    blockedFaces
-    {
-        type            kqRWallFunction;
-        value           uniform 1.5;
-    }
-    baffleWall
-    {
-        type            kqRWallFunction;
-        value           uniform 1.5;
-    }
-    baffleCyclic_half0
-    {
-        type            cyclic;
-    }
-    baffleCyclic_half1
-    {
-        type            cyclic;
-    }
-}
-
-
-// ************************************************************************* //

tutorials/combustion/PDRFoam/flamePropagationWithObstacles/0.orig/nut

@@ -1,64 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  v2012                                 |
-|   \\  /    A nd           | Website:  www.openfoam.com                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       volScalarField;
-    location    "0";
-    object      nut;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-dimensions      [0 2 -1 0 0 0 0];
-
-internalField   uniform 0;
-
-boundaryField
-{
-    outer
-    {
-        type            calculated;
-        value           uniform 0;
-    }
-    ground
-    {
-        type            nutkWallFunction;
-        Cmu             0.09;
-        kappa           0.41;
-        E               9.8;
-        value           uniform 0;
-    }
-    blockedFaces
-    {
-        type            nutkWallFunction;
-        Cmu             0.09;
-        kappa           0.41;
-        E               9.8;
-        value           nonuniform 0();
-    }
-    baffleWall
-    {
-        type            nutkWallFunction;
-        Cmu             0.09;
-        kappa           0.41;
-        E               9.8;
-        value           nonuniform 0();
-    }
-    baffleCyclic_half0
-    {
-        type            cyclic;
-    }
-    baffleCyclic_half1
-    {
-        type            cyclic;
-    }
-}
-
-
-// ************************************************************************* //

tutorials/combustion/PDRFoam/flamePropagationWithObstacles/0.orig/p

@@ -1,55 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  v2012                                 |
-|   \\  /    A nd           | Website:  www.openfoam.com                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       volScalarField;
-    location    "0";
-    object      p;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-dimensions      [1 -1 -2 0 0 0 0];
-
-internalField   uniform 100000;
-
-boundaryField
-{
-    outer
-    {
-        type            waveTransmissive;
-        gamma           1.3;
-        fieldInf        100000;
-        lInf            5;
-        value           uniform 100000;
-    }
-    ground
-    {
-        type            zeroGradient;
-    }
-    blockedFaces
-    {
-        type            zeroGradient;
-    }
-    baffleWall
-    {
-        type            zeroGradient;
-    }
-    baffleCyclic_half0
-    {
-        type            cyclic;
-    }
-    baffleCyclic_half1
-    {
-        type            cyclic;
-    }
-}
-
-
-// ************************************************************************* //

tutorials/combustion/PDRFoam/flamePropagationWithObstacles/Allclean

@@ -1,10 +0,0 @@
-#!/bin/sh
-cd "${0%/*}" || exit                                # Run from this directory
-. ${WM_PROJECT_DIR:?}/bin/tools/CleanFunctions      # Tutorial clean functions
-#------------------------------------------------------------------------------
-
-cleanCase0
-
-rm -rf VTK
-
-#------------------------------------------------------------------------------

tutorials/combustion/PDRFoam/flamePropagationWithObstacles/Allrun

@@ -1,15 +0,0 @@
-#!/bin/sh
-cd "${0%/*}" || exit                                # Run from this directory
-. ${WM_PROJECT_DIR:?}/bin/tools/RunFunctions        # Tutorial run functions
-#------------------------------------------------------------------------------
-
-restore0Dir
-runApplication blockMesh
-runApplication changeDictionary
-runApplication topoSet
-
-runApplication PDRMesh -overwrite
-
-runApplication $(getApplication)
-
-#------------------------------------------------------------------------------

tutorials/combustion/PDRFoam/flamePropagationWithObstacles/README

@@ -1,8 +0,0 @@
-PDR test case
-
-Step to introduce the PDR fields:
-
-1) Create zero-size patches for wall or/and coupled baffles in
-the boundary file.
-2) Specify the boundary contitions for these patches in the fields.
-3) Create the new PDR mesh using the PDRMesh utility.

tutorials/combustion/PDRFoam/flamePropagationWithObstacles/constant/PDRProperties

@@ -1,103 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  v2012                                 |
-|   \\  /    A nd           | Website:  www.openfoam.com                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       dictionary;
-    location    "constant";
-    object      PDRProperties;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-PDRDragModel basic;
-
-basicCoeffs
-{
-    drag     on;
-    Csu      0.5;
-    Csk      0.05;
-}
-
-
-XiModel transport;
-
-transportCoeffs
-{
-    XiShapeCoef 1;
-}
-
-
-XiEqModel instability;
-
-instabilityCoeffs
-{
-    XiEqIn 2.5;
-
-    XiEqModel basicSubGrid;
-
-    basicSubGridCoeffs
-    {
-        XiEqModel SCOPEBlend;
-
-        SCOPEBlendCoeffs
-        {
-            XiEqModelL
-            {
-                XiEqModel       Gulder;
-
-                GulderCoeffs
-                {
-                    XiEqCoef   0.62;
-                    uPrimeCoef      1.0;
-                    subGridSchelkin true;
-                }
-            }
-
-            XiEqModelH
-            {
-                XiEqModel       SCOPEXiEq;
-
-                SCOPEXiEqCoeffs
-                {
-                    XiEqCoef   1.6;
-                    XiEqExp    0.33333;
-                    lCoef      0.336;
-                    uPrimeCoef      1.0;
-                    subGridSchelkin true;
-                }
-            }
-        }
-    }
-}
-
-
-XiGModel instabilityG;
-
-instabilityGCoeffs
-{
-    lambdaIn        lambdaIn   [0 1 0 0 0 0 0] 0.6;
-    GIn             GIn        [0 0 -1 0 0 0 0] 1.917;
-
-    XiGModel basicSubGridG;
-
-    basicSubGridGCoeffs
-    {
-        k1 0.5;
-
-        XiGModel KTS;
-
-        KTSCoeffs
-        {
-             GEtaCoef   0.28;
-        }
-    }
-}
-
-
-// ************************************************************************* //

tutorials/combustion/PDRFoam/flamePropagationWithObstacles/constant/combustionProperties

@@ -1,46 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  v2012                                 |
-|   \\  /    A nd           | Website:  www.openfoam.com                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       dictionary;
-    location    "constant";
-    object      combustionProperties;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-laminarFlameSpeedCorrelation SCOPE;
-
-fuel        Propane;
-
-fuelFile    "fuels/propane.dat";
-
-ignite      yes;
-
-ignitionSites
-(
-
-    {
-        location        (1.5 1.5 0.5);
-        diameter        0.1;
-        start           1E-05;
-        duration        0.05;
-        strength        60.0;
-    }
-);
-
-ignitionSphereFraction  1;
-
-ignitionThickness   ignitionThickness [0 1 0 0 0 0 0] 0;
-
-ignitionCircleFraction  0;
-
-ignitionKernelArea ignitionKernelArea [0 2 0 0 0 0 0] 0;
-
-// ************************************************************************* //

tutorials/combustion/PDRFoam/flamePropagationWithObstacles/constant/dynamicMeshDict

@@ -1,53 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  v2012                                 |
-|   \\  /    A nd           | Website:  www.openfoam.com                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       dictionary;
-    object      dynamicMeshDict;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-dynamicFvMesh      dynamicRefineFvMesh;
-
-// Refine every refineInterval timesteps
-refineInterval 1;
-
-// Maximum refinement level (starts from 0)
-maxRefinement 2;
-
-// Maximum cell limit (approximate)
-maxCells 10000;
-
-// volScalarField to base refinement on
-field normalisedGradP;
-
-nBufferLayers 1;
-
-dumpLevel true;
-
-lowerRefineLevel 0.5;
-upperRefineLevel 1.5;
-
-unrefineLevel 0.5;
-
-nBufferLayers 1;
-// Newly introduced patch points optionally get projected onto a surface
-//projectSurfaces ("fixedWalls4.stl");
-//projectPatches (fixedWalls);
-// Maximum project distance
-//projectDistance 1;
-
-// Fluxes to adapt. For newly created faces or split faces the flux
-// gets estimated from an interpolated volVectorField ('velocity')
-// First is name of the flux to adapt, second is velocity that will
-// be interpolated and inner-producted with the face area vector.
-correctFluxes ((phi rhoU) (phi_0 none));
-
-// ************************************************************************* //

tutorials/combustion/PDRFoam/flamePropagationWithObstacles/constant/g

@@ -1,22 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  v2012                                 |
-|   \\  /    A nd           | Website:  www.openfoam.com                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       uniformDimensionedVectorField;
-    location    "constant";
-    object      g;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-dimensions      [0 1 -2 0 0 0 0];
-value           (0 0 -9.8);
-
-
-// ************************************************************************* //

tutorials/combustion/PDRFoam/flamePropagationWithObstacles/constant/thermophysicalProperties

@@ -1,95 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  v2012                                 |
-|   \\  /    A nd           | Website:  www.openfoam.com                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       dictionary;
-    location    "constant";
-    object      thermophysicalProperties;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-thermoType
-{
-    type            heheuPsiThermo;
-    mixture         inhomogeneousMixture;
-    transport       sutherland;
-    thermo          janaf;
-    equationOfState perfectGas;
-    specie          specie;
-    energy          absoluteEnthalpy;
-}
-
-stoichiometricAirFuelMassRatio  stoichiometricAirFuelMassRatio [0 0 0 0 0 0 0] 15.5776;
-
-fuel
-{
-    specie
-    {
-        molWeight       44.0962;
-    }
-    thermodynamics
-    {
-        Tlow            200;
-        Thigh           5000;
-        Tcommon         1000;
-        highCpCoeffs    ( 7.53414 0.0188722 -6.27185e-06 9.14756e-10 -4.78381e-14 -16467.5 -17.8923 );
-        lowCpCoeffs     ( 0.933554 0.0264246 6.10597e-06 -2.19775e-08 9.51493e-12 -13958.5 19.2017  );
-    }
-    transport
-    {
-        As              1.67212e-06;
-        Ts              170.672;
-    }
-}
-
-
-oxidant
-{
-    specie
-    {
-        molWeight       28.8504;
-    }
-    thermodynamics
-    {
-        Tlow            200;
-        Thigh           6000;
-        Tcommon         1000;
-        highCpCoeffs    ( 3.10131 0.00124137 -4.18816e-07 6.64158e-11 -3.91274e-15 -985.266 5.35597 );
-        lowCpCoeffs     ( 3.58378 -0.000727005 1.67057e-06 -1.09203e-10 -4.31765e-13 -1050.53 3.11239 );
-    }
-    transport
-    {
-        As              1.67212e-06;
-        Ts              170.672;
-    }
-}
-
-burntProducts
-{
-    specie
-    {
-        molWeight      28.3233;
-    }
-    thermodynamics
-    {
-        Tlow           200;
-        Thigh          6000;
-        Tcommon        1000;
-        highCpCoeffs   ( 3.10558 0.00179747 -5.94696e-07 9.05605e-11 -5.08443e-15 -11003.6 5.12104 );
-        lowCpCoeffs    ( 3.49796 0.000638555 -1.83881e-07 1.20989e-09 -7.68691e-13 -11080.5 3.18188 );
-    }
-    transport
-    {
-        As             1.67212e-06;
-        Ts             170.672;
-    }
-}
-
-// ************************************************************************* //

tutorials/combustion/PDRFoam/flamePropagationWithObstacles/constant/turbulenceProperties

@@ -1,42 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  v2012                                 |
-|   \\  /    A nd           | Website:  www.openfoam.com                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       dictionary;
-    location    "constant";
-    object      turbulenceProperties;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-simulationType  RAS;
-
-RAS
-{
-    RASModel PDRkEpsilon;
-
-    PDRkEpsilonCoeffs
-    {
-        Cmu              0.09;
-        C1               1.44;
-        C2               1.92;
-        C3               0;
-        C4               0.1;
-        alphah           1;
-        alphak           1;
-        alphaEps         0.76923;
-    }
-
-    turbulence      on;
-
-    printCoeffs     on;
-}
-
-
-// ************************************************************************* //

tutorials/combustion/PDRFoam/flamePropagationWithObstacles/fuels/propane.dat

@@ -1,72 +0,0 @@
-SCOPECoeffs
-{
-    lowerFlammabilityLimit      0.522;
-    upperFlammabilityLimit      2.362;
-
-    lowerSuPolynomial
-    {
-        lowerLimit 0.6;
-        upperLimit 1.2;
-        coefficients
-        (
-        0.45
-        0.713364389
-        -2.359004778
-        -2.629755677
-        -29.11933997
-        -55.34224225
-        0.0
-        );
-    }
-
-    upperSuPolynomial
-    {
-        lowerLimit 1.1;
-        upperLimit 1.64;
-        coefficients
-        (
-        0.419664668
-        2.280811555
-        -11.53872754
-        2.98656153
-        21.04870808
-        -16.09645303
-        0.0
-        );
-    }
-
-    Texp                    2;
-    pexp                    -0.5;
-
-    lowerMaPolynomial
-    {
-        lowerLimit 0.5;
-        upperLimit 0.75;
-        coefficients
-        (
-        3.739047
-        -5.12414
-        -38.1444
-        -134.813
-        -224.633
-        -144.27
-        0
-        );
-    }
-
-    upperMaPolynomial
-    {
-        lowerLimit 0.75;
-        upperLimit 1.0;
-        coefficients
-        (
-        3.739047
-        -5.12414
-        -38.1444
-        -134.813
-        -224.633
-        -144.27
-        0
-        );
-    }
-}

tutorials/combustion/PDRFoam/flamePropagationWithObstacles/system/PDRMeshDict

@@ -1,38 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  v2012                                 |
-|   \\  /    A nd           | Website:  www.openfoam.com                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       dictionary;
-    object      PDRMeshDict;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-//- Per faceSet the patch the faces should go into blocked baffles
-blockedFaces ((blockedFacesSet blockedFaces));
-
-//- Per faceSet the patch the faces should go into coupled baffles
-coupledFaces
-{
-
-    coupledFacesSet
-    {
-        wallPatch                   baffleWall;
-        cyclicMasterPatch           baffleCyclic_half0;
-    }
-}
-
-//- Name of cellSet that holds the cells to fully remove
-blockedCells blockedCellsSet;
-
-//- All exposed faces that are not specified in blockedFaces go into
-//  this patch
-defaultPatch outer;
-
-// ************************************************************************* //

tutorials/combustion/PDRFoam/flamePropagationWithObstacles/system/blockMeshDict

@@ -1,269 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  v2012                                 |
-|   \\  /    A nd           | Website:  www.openfoam.com                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       dictionary;
-    object      blockMeshDict;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-scale   1.0;
-
-vertices
-(
-    (0.0  0.0 0.0)
-    (1.0  0.0 0.0)
-    (2.0  0.0 0.0)
-    (3.0  0.0 0.0)
-    (0.0  1.0 0.0)
-    (1.0  1.0 0.0)
-    (2.0  1.0 0.0)
-    (3.0  1.0 0.0)
-    (0.0  2.0 0.0)
-    (1.0  2.0 0.0)
-    (2.0  2.0 0.0)
-    (3.0  2.0 0.0)
-    (0.0  3.0 0.0)
-    (1.0  3.0 0.0)
-    (2.0  3.0 0.0)
-    (3.0  3.0 0.0)
-    (0.0  0.0 1.0)
-    (1.0  0.0 1.0)
-    (2.0  0.0 1.0)
-    (3.0  0.0 1.0)
-    (0.0  1.0 1.0)
-    (1.0  1.0 1.0)
-    (2.0  1.0 1.0)
-    (3.0  1.0 1.0)
-    (0.0  2.0 1.0)
-    (1.0  2.0 1.0)
-    (2.0  2.0 1.0)
-    (3.0  2.0 1.0)
-    (0.0  3.0 1.0)
-    (1.0  3.0 1.0)
-    (2.0  3.0 1.0)
-    (3.0  3.0 1.0)
-    (0.0  0.0 2.0)
-    (1.0  0.0 2.0)
-    (2.0  0.0 2.0)
-    (3.0  0.0 2.0)
-    (0.0  1.0 2.0)
-    (1.0  1.0 2.0)
-    (2.0  1.0 2.0)
-    (3.0  1.0 2.0)
-    (0.0  2.0 2.0)
-    (1.0  2.0 2.0)
-    (2.0  2.0 2.0)
-    (3.0  2.0 2.0)
-    (0.0  3.0 2.0)
-    (1.0  3.0 2.0)
-    (2.0  3.0 2.0)
-    (3.0  3.0 2.0)
-    (-3.53553390593274  -3.53553390593274 0)
-    (-0.58113883008419  -4.74341649025257 0)
-    (3.58113883008419  -4.74341649025257 0)
-    (6.53553390593274  -3.53553390593274 0)
-    (-4.74341649025257  -0.58113883008419 0)
-    (-2.53553390593274  -2.53553390593274 0)
-    (5.53553390593274  -2.53553390593274 0)
-    (7.74341649025257  -0.58113883008419 0)
-    (-4.74341649025257  3.58113883008419 0)
-    (-2.53553390593274  5.53553390593274 0)
-    (5.53553390593274  5.53553390593274 0)
-    (7.74341649025257  3.58113883008419 0)
-    (-3.53553390593274  6.53553390593274 0)
-    (-0.58113883008419  7.74341649025257 0)
-    (3.58113883008419  7.74341649025257 0)
-    (6.53553390593274  6.53553390593274 0)
-    (-3.19801074533416  -3.19801074533416 3.1320071635561)
-    (-0.336306209562122  -4.00891862868637 3.67261241912424)
-    (3.33630620956212  -4.00891862868637 3.67261241912424)
-    (6.19801074533416  -3.19801074533416 3.1320071635561)
-    (-4.00891862868637  -0.336306209562122 3.67261241912424)
-    (7.00891862868637  -0.336306209562122 3.67261241912424)
-    (-4.00891862868637  3.33630620956212 3.67261241912424)
-    (7.00891862868637  3.33630620956212 3.67261241912424)
-    (-3.19801074533416  6.19801074533416 3.1320071635561)
-    (-0.336306209562122  7.00891862868637 3.67261241912424)
-    (3.33630620956212  7.00891862868637 3.67261241912424)
-    (6.19801074533416  6.19801074533416 3.1320071635561)
-    (-2.57247877713763  -2.57247877713763 5.42997170285018)
-    (0.0194193243090797  -2.94174202707276 5.92232270276368)
-    (2.98058067569092  -2.94174202707276 5.92232270276368)
-    (5.57247877713763  -2.57247877713763 5.42997170285018)
-    (-2.94174202707276  0.0194193243090797 5.92232270276368)
-    (-0.178511301977579  -0.178511301977579 6.71404520791032)
-    (3.17851130197758  -0.178511301977579 6.71404520791032)
-    (5.94174202707276  0.0194193243090797 5.92232270276368)
-    (-2.94174202707276  2.98058067569092 5.92232270276368)
-    (-0.178511301977579  3.17851130197758 6.71404520791032)
-    (3.17851130197758  3.17851130197758 6.71404520791032)
-    (5.94174202707276  2.98058067569092 5.92232270276368)
-    (-2.57247877713763  5.57247877713763 5.42997170285018)
-    (0.0194193243090797  5.94174202707276 5.92232270276368)
-    (2.98058067569092  5.94174202707276 5.92232270276368)
-    (5.57247877713763  5.57247877713763 5.42997170285018)
-);
-
-blocks
-(
-    hex  (0 1 5 4 16 17 21 20 ) (5 5 5) simpleGrading (1 1 1)
-    hex  (1 2 6 5 17 18 22 21 ) (5 5 5) simpleGrading (1 1 1)
-    hex  (2 3 7 6 18 19 23 22 ) (5 5 5) simpleGrading (1 1 1)
-    hex  (4 5 9 8 20 21 25 24 ) (5 5 5) simpleGrading (1 1 1)
-    hex  (5 6 10 9 21 22 26 25 ) (5 5 5) simpleGrading (1 1 1)
-    hex  (6 7 11 10 22 23 27 26 )(5 5 5) simpleGrading (1 1 1)
-    hex  (8 9 13 12 24 25 29 28 ) (5 5 5) simpleGrading (1 1 1)
-    hex  (9 10 14 13 25 26 30 29 ) (5 5 5) simpleGrading (1 1 1)
-    hex  (10 11 15 14 26 27 31 30 ) (5 5 5) simpleGrading (1 1 1)
-    hex  (16 17 21 20 32 33 37 36 ) (5 5 5) simpleGrading (1 1 1)
-    hex  (17 18 22 21 33 34 38 37 ) (5 5 5) simpleGrading (1 1 1)
-    hex  (18 19 23 22 34 35 39 38 ) (5 5 5) simpleGrading (1 1 1)
-    hex  (20 21 25 24 36 37 41 40 ) (5 5 5) simpleGrading (1 1 1)
-    hex  (21 22 26 25 37 38 42 41 ) (5 5 5) simpleGrading (1 1 1)
-    hex  (22 23 27 26 38 39 43 42 ) (5 5 5) simpleGrading (1 1 1)
-    hex  (24 25 29 28 40 41 45 44 ) (5 5 5) simpleGrading (1 1 1)
-    hex  (25 26 30 29 41 42 46 45 ) (5 5 5) simpleGrading (1 1 1)
-    hex  (26 27 31 30 42 43 47 46 ) (5 5 5) simpleGrading (1 1 1)
-    hex  (47 46 42 43 91 90 86 87) (5 5 7) simpleGrading (1 1 2.985984)
-    hex  (46 45 41 42 90 89 85 86) (5 5 7) simpleGrading (1 1 2.985984)
-    hex (45 44 40 41 89 88 84 85) (5 5 7) simpleGrading (1 1 2.985984)
-    hex (43 42 38 39 87 86 82 83) (5 5 7) simpleGrading (1 1 2.985984)
-    hex (42 41 37 38 86 85 81 82) (5 5 7) simpleGrading (1 1 2.985984)
-    hex (41 40 36 37 85 84 80 81) (5 5 7) simpleGrading (1 1 2.985984)
-    hex (39 38 34 35 83 82 78 79) (5 5 7) simpleGrading (1 1 2.985984)
-    hex (38 37 33 34 82 81 77 78) (5 5 7) simpleGrading (1 1 2.985984)
-    hex (37 36 32 33 81 80 76 77) (5 5 7) simpleGrading (1 1 2.985984)
-    hex (0 16 20 4 48 64 68 52) (5 5 7) simpleGrading (1 1 2.985984)
-    hex (16 32 36 20 64 76 80 68) (5 5 7) simpleGrading (1 1 2.985984)
-    hex (4 20 24 8 52 68 70 56) (5 5 7) simpleGrading (1 1 2.985984)
-    hex (20 36 40 24 68 80 84 70) (5 5 7) simpleGrading (1 1 2.985984)
-    hex (8 24 28 12 56 70 72 60) (5 5 7) simpleGrading (1 1 2.985984)
-    hex (24 40 44 28 70 84 88 72) (5 5 7) simpleGrading (1 1 2.985984)
-    hex (47 43 27 31 91 87 71 75) (5 5 7) simpleGrading (1 1 2.985984)
-    hex (43 39 23 27 87 83 69 71) (5 5 7) simpleGrading (1 1 2.985984)
-    hex (39 35 19 23 83 79 67 69) (5 5 7) simpleGrading (1 1 2.985984)
-    hex (31 27 11 15 75 71 59 63) (5 5 7) simpleGrading (1 1 2.985984)
-    hex (27 23 7 11 71 69 55 59) (5 5 7) simpleGrading (1 1 2.985984)
-    hex (23 19 3 7 69 67 51 55) (5 5 7) simpleGrading (1 1 2.985984)
-    hex (0 1 17 16 48 49 65 64) (5 5 7) simpleGrading (1 1 2.985984)
-    hex (1 2 18 17 49 50 66 65) (5 5 7) simpleGrading (1 1 2.985984)
-    hex (2 3 19 18 50 51 67 66) (5 5 7) simpleGrading (1 1 2.985984)
-    hex (16 17 33 32 64 65 77 76) (5 5 7) simpleGrading (1 1 2.985984)
-    hex (17 18 34 33 65 66 78 77) (5 5 7) simpleGrading (1 1 2.985984)
-    hex (18 19 35 34 66 67 79 78) (5 5 7) simpleGrading (1 1 2.985984)
-    hex (47 31 30 46 91 75 74 90) (5 5 7) simpleGrading (1 1 2.985984)
-    hex (31 15 14 30 75 63 62 74) (5 5 7) simpleGrading (1 1 2.985984)
-    hex (46 30 29 45 90 74 73 89) (5 5 7) simpleGrading (1 1 2.985984)
-    hex (30 14 13 29 74 62 61 73) (5 5 7) simpleGrading (1 1 2.985984)
-    hex (45 29 28 44 89 73 72 88) (5 5 7) simpleGrading (1 1 2.985984)
-    hex (29 13 12 28 73 61 60 72) (5 5 7) simpleGrading (1 1 2.985984)
-);
-
-boundary
-(
-    outer
-    {
-        type    patch;
-        faces
-        (
-            (91 90 86 87)
-            (90 89 85 86)
-            (89 88 84 85)
-            (87 86 82 83)
-            (86 85 81 82)
-            (85 84 80 81)
-            (83 82 78 79)
-            (82 81 77 78)
-            (81 80 76 77)
-            (48 64 68 52)
-            (64 76 80 68)
-            (52 68 70 56)
-            (68 80 84 70)
-            (56 70 72 60)
-            (70 84 88 72)
-            (91 87 71 75)
-            (87 83 69 71)
-            (83 79 67 69)
-            (75 71 59 63)
-            (71 69 55 59)
-            (69 67 51 55)
-            (48 49 65 64)
-            (49 50 66 65)
-            (50 51 67 66)
-            (64 65 77 76)
-            (65 66 78 77)
-            (66 67 79 78)
-            (91 75 74 90)
-            (75 63 62 74)
-            (90 74 73 89)
-            (74 62 61 73)
-            (89 73 72 88)
-            (73 61 60 72)
-        );
-    }
-
-    ground
-    {
-        type wall;
-        faces
-        (
-            (0 4 5 1)
-            (1 5 6 2)
-            (2 6 7 3)
-            (4 8 9 5)
-            (5 9 10 6)
-            (6 10 11 7)
-            (8 12 13 9)
-            (9 13 14 10)
-            (10 14 15 11)
-            (4 0 48 52)
-            (8 4 52 56)
-            (12 8 56 60)
-            (11 15 63 59)
-            (7 11 59 55)
-            (3 7 55 51)
-            (0 1 49 48)
-            (1 2 50 49)
-            (2 3 51 50)
-            (15 14 62 63)
-            (14 13 61 62)
-            (13 12 60 61)
-        );
-    }
-
-    blockedFaces
-    {
-        type wall;
-        faces ();
-    }
-
-    baffleWall
-    {
-        type wall;
-        faces ();
-    }
-
-    baffleCyclic_half0
-    {
-        type cyclic;
-        neighbourPatch  baffleCyclic_half1;
-        faces ();
-    }
-
-    baffleCyclic_half1
-    {
-        type cyclic;
-        neighbourPatch  baffleCyclic_half0;
-        faces ();
-    }
-);
-
-// ************************************************************************* //

tutorials/combustion/PDRFoam/flamePropagationWithObstacles/system/changeDictionaryDict

@@ -1,462 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  v2012                                 |
-|   \\  /    A nd           | Website:  www.openfoam.com                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       dictionary;
-    object      changeDictionaryDict;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-p
-{
-    boundaryField
-    {
-        blockedFaces
-        {
-            type        zeroGradient;
-        }
-        baffleWall
-        {
-            type        zeroGradient;
-        }
-        baffleCyclic_half0
-        {
-            type        cyclic;
-        }
-        baffleCyclic_half1
-        {
-            type        cyclic;
-        }
-    }
-}
-
-U
-{
-    boundaryField
-    {
-        blockedFaces
-        {
-            type            fixedValue;
-            value           uniform (0 0 0);
-        }
-        baffleWall
-        {
-            type                activePressureForceBaffleVelocity;
-            cyclicPatch         baffleCyclic_half0;
-            orientation         1;
-            openingTime         0.01;
-            maxOpenFractionDelta 0.1;
-            openFraction        0;
-            minThresholdValue   8000;
-            forceBased          0;
-            value               uniform (0 0 0);
-        }
-        baffleCyclic_half0
-        {
-            type        cyclic;
-        }
-        baffleCyclic_half1
-        {
-            type        cyclic;
-        }
-    }
-}
-
-epsilon
-{
-    boundaryField
-    {
-        blockedFaces
-        {
-            type            epsilonWallFunction;
-            value           uniform 0.1;
-        }
-        baffleWall
-        {
-            type            epsilonWallFunction;
-            value           uniform 1e-5;
-        }
-        baffleCyclic_half0
-        {
-            type        cyclic;
-        }
-        baffleCyclic_half1
-        {
-            type        cyclic;
-        }
-    }
-}
-
-k
-{
-    boundaryField
-    {
-        blockedFaces
-        {
-            type        kqRWallFunction;
-            value       uniform 1.5;
-     }
-        baffleWall
-        {
-            type        kqRWallFunction;
-            value       uniform 1.5;
-        }
-        baffleCyclic_half0
-        {
-            type        cyclic;
-        }
-        baffleCyclic_half1
-        {
-            type        cyclic;
-        }
-    }
-}
-
-Su
-{
-    boundaryField
-    {
-        blockedFaces
-        {
-            type        zeroGradient;
-        }
-        baffleWall
-        {
-            type        zeroGradient;
-        }
-        baffleCyclic_half0
-        {
-            type        cyclic;
-        }
-        baffleCyclic_half1
-        {
-            type        cyclic;
-        }
-    }
-}
-
-Xi
-{
-    boundaryField
-    {
-        blockedFaces
-        {
-            type        zeroGradient;
-            value       uniform 1;
-        }
-         baffleWall
-        {
-            type        zeroGradient;
-            value       uniform 1;
-        }
-        baffleCyclic_half0
-        {
-            type        cyclic;
-        }
-        baffleCyclic_half1
-        {
-            type        cyclic;
-        }
-    }
-}
-
-T
-{
-    boundaryField
-    {
-        blockedFaces
-        {
-            type        zeroGradient;
-            value       uniform 300;
-        }
-        baffleWall
-        {
-            type        zeroGradient;
-            value       uniform 300;
-        }
-        baffleCyclic_half0
-        {
-            type        cyclic;
-        }
-        baffleCyclic_half1
-        {
-            type        cyclic;
-        }
-    }
-}
-
-Tu
-{
-    boundaryField
-    {
-        blockedFaces
-        {
-            type        zeroGradient;
-            value       uniform 300;
-        }
-        baffleWall
-        {
-            type        zeroGradient;
-            value       uniform 300;
-        }
-        baffleCyclic_half0
-        {
-            type        cyclic;
-        }
-        baffleCyclic_half1
-        {
-            type        cyclic;
-        }
-    }
-}
-
-b
-{
-    boundaryField
-    {
-        blockedFaces
-        {
-            type        zeroGradient;
-            value       uniform 1;
-        }
-        baffleWall
-        {
-            type        zeroGradient;
-            value       uniform 1;
-        }
-        baffleCyclic_half0
-        {
-            type        cyclic;
-        }
-        baffleCyclic_half1
-        {
-            type        cyclic;
-        }
-    }
-}
-
-ft
-{
-    boundaryField
-    {
-        blockedFaces
-        {
-            type        zeroGradient;
-        }
-        baffleWall
-        {
-            type        zeroGradient;
-            value       uniform 0.06024096;
-        }
-        baffleCyclic_half0
-        {
-            type        cyclic;
-        }
-        baffleCyclic_half1
-        {
-            type        cyclic;
-        }
-    }
-}
-
-Aw
-{
-    boundaryField
-    {
-        blockedFaces
-        {
-            type        zeroGradient;
-            value       uniform 0.0;
-        }
-        baffleWall
-        {
-            type        zeroGradient;
-            value       uniform 0.0;
-        }
-        baffleCyclic_half0
-        {
-            type        cyclic;
-        }
-        baffleCyclic_half1
-        {
-            type        cyclic;
-        }
-    }
-}
-B
-{
-    boundaryField
-    {
-        blockedFaces
-        {
-            type        zeroGradient;
-            value       uniform (0 0 0 0 0 0);
-        }
-        baffleWall
-        {
-            type        zeroGradient;
-            value       uniform (0 0 0 0 0 0);
-        }
-        baffleCyclic_half0
-        {
-            type        cyclic;
-        }
-        baffleCyclic_half1
-        {
-            type        cyclic;
-        }
-    }
-}
-betav
-{
-    boundaryField
-    {
-        blockedFaces
-        {
-            type        zeroGradient;
-            value       uniform 0;
-        }
-        baffleWall
-        {
-            type        zeroGradient;
-            value       uniform 0;
-        }
-        baffleCyclic_half0
-        {
-            type        cyclic;
-        }
-        baffleCyclic_half1
-        {
-            type        cyclic;
-        }
-    }
-}
-CR
-{
-    boundaryField
-    {
-        blockedFaces
-        {
-            type        zeroGradient;
-            value       uniform (0 0 0 0 0 0);
-        }
-        baffleWall
-        {
-            type        zeroGradient;
-            value       uniform (0 0 0 0 0 0);
-        }
-        baffleCyclic_half0
-        {
-            type        cyclic;
-        }
-        baffleCyclic_half1
-        {
-            type        cyclic;
-        }
-    }
-}
-CT
-{
-    boundaryField
-    {
-        blockedFaces
-        {
-            type        zeroGradient;
-            value       uniform (0 0 0 0 0 0);
-        }
-        baffleWall
-        {
-            type        zeroGradient;
-            value       uniform (0 0 0 0 0 0);
-        }
-        baffleCyclic_half0
-        {
-            type        cyclic;
-        }
-        baffleCyclic_half1
-        {
-            type        cyclic;
-        }
-    }
-}
-Nv
-{
-    boundaryField
-    {
-        blockedFaces
-        {
-            type        zeroGradient;
-            value       uniform 0;
-        }
-        baffleWall
-        {
-            type        zeroGradient;
-            value       uniform 0;
-        }
-        baffleCyclic_half0
-        {
-            type        cyclic;
-        }
-        baffleCyclic_half1
-        {
-            type        cyclic;
-        }
-    }
-}
-nsv
-{
-    boundaryField
-    {
-        blockedFaces
-        {
-            type        zeroGradient;
-            value       uniform (0 0 0 0 0 0);
-        }
-        baffleWall
-        {
-            type        zeroGradient;
-            value       uniform (0 0 0 0 0 0);
-        }
-        baffleCyclic_half0
-        {
-            type        cyclic;
-        }
-        baffleCyclic_half1
-        {
-            type        cyclic;
-        }
-    }
-}
-Lobs
-{
-    boundaryField
-    {
-        blockedFaces
-        {
-            type        zeroGradient;
-            value       uniform 0;
-        }
-        baffleWall
-        {
-            type        zeroGradient;
-            value       uniform 0;
-        }
-        baffleCyclic_half0
-        {
-            type        cyclic;
-        }
-        baffleCyclic_half1
-        {
-            type        cyclic;
-        }
-    }
-}
-
-// ************************************************************************* //

tutorials/combustion/PDRFoam/flamePropagationWithObstacles/system/controlDict

@@ -1,54 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  v2012                                 |
-|   \\  /    A nd           | Website:  www.openfoam.com                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       dictionary;
-    location    "system";
-    object      controlDict;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-application     PDRFoam;
-
-startFrom       latestTime;
-
-startTime       0;
-
-stopAt          endTime;
-
-endTime         0.2;
-
-deltaT          5e-06;
-
-writeControl    adjustable;
-
-writeInterval   0.01;
-
-purgeWrite      0;
-
-writeFormat     ascii;
-
-writePrecision  8;
-
-writeCompression off;
-
-timeFormat      general;
-
-timePrecision   6;
-
-runTimeModifiable yes;
-
-adjustTimeStep  yes;
-
-maxCo           0.2;
-
-maxDeltaT       1;
-
-// ************************************************************************* //

tutorials/combustion/PDRFoam/flamePropagationWithObstacles/system/fvSchemes

@@ -1,72 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  v2012                                 |
-|   \\  /    A nd           | Website:  www.openfoam.com                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       dictionary;
-    location    "system";
-    object      fvSchemes;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-ddtSchemes
-{
-    default         Euler;
-}
-
-gradSchemes
-{
-    default         Gauss linear;
-    grad(p)         Gauss linear;
-}
-
-divSchemes
-{
-    default         none;
-
-    div(phi,U)      Gauss limitedLinearV 1;
-    div(phi,K)      Gauss limitedLinear 1;
-    div(phid,p)     Gauss limitedLinear 1;
-    div(phi,k)      Gauss limitedLinear 1;
-    div(phi,epsilon) Gauss limitedLinear 1;
-
-    div(phiXi,Xi)   Gauss limitedLinear 1;
-    div(phiSt,b)    Gauss limitedLinear01 1;
-
-    div(phi,ft_b_ha_hau) Gauss multivariateSelection
-    {
-        ft              limitedLinear01 1;
-        b               limitedLinear01 1;
-        Xi              limitedLinear 1;
-        ha              limitedLinear 1;
-        hau             limitedLinear 1;
-    };
-
-    div((Su*grad(b)))               Gauss linear;
-    div((U+((Su*Xi)*grad(b))))      Gauss linear;
-    div(((rho*nuEff)*dev2(T(grad(U)))))   Gauss linear;
-    div(U)                          Gauss linear;
-}
-
-laplacianSchemes
-{
-    default         Gauss linear limited corrected 0.33;
-}
-
-interpolationSchemes
-{
-    default         linear;
-}
-
-snGradSchemes
-{
-    default         limited corrected 0.33;
-}
-
-// ************************************************************************* //

tutorials/combustion/PDRFoam/flamePropagationWithObstacles/system/fvSolution

@@ -1,76 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  v2012                                 |
-|   \\  /    A nd           | Website:  www.openfoam.com                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       dictionary;
-    location    "system";
-    object      fvSolution;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-solvers
-{
-    rho
-    {
-        solver          PCG;
-        preconditioner  DIC;
-        tolerance       1e-05;
-        relTol          0.01;
-    };
-
-    rhoFinal
-    {
-        $rho;
-        tolerance       1e-06;
-        relTol          0;
-    };
-
-    p
-    {
-        solver          PCG;
-        preconditioner  DIC;
-        tolerance       1e-6;
-        relTol          0.01;
-    };
-
-    pFinal
-    {
-        $p;
-        tolerance       1e-6;
-        relTol          0;
-    };
-
-    "(b|Xi|ft|ha|hau|k|epsilon)"
-    {
-        solver          PBiCGStab;
-        preconditioner  DILU;
-        tolerance       1e-06;
-        relTol          0.1;
-    }
-
-    "(b|Xi|ft|ha|hau|k|epsilon)Final"
-    {
-        solver          PBiCGStab;
-        preconditioner  DILU;
-        tolerance       1e-06;
-        relTol          0;
-    }
-}
-
-PIMPLE
-{
-    nCorrectors     2;
-    nOuterCorrectors  1;
-    nNonOrthogonalCorrectors 0;
-    momentumPredictor yes;
-}
-
-
-// ************************************************************************* //

tutorials/combustion/PDRFoam/flamePropagationWithObstacles/system/topoSetDict

@@ -1,63 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  v2012                                 |
-|   \\  /    A nd           | Website:  www.openfoam.com                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       dictionary;
-    location    "system";
-    object      topoSetDict.1;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-actions
-(
-    {
-        name    blockedFacesSet;
-        type    faceSet;
-        action  new;
-        source  boxToFace;
-        box     (1 1 0.99) (2 2 1);
-    }
-
-    {
-        name    blockedCellsSet;
-        type    cellSet;
-        action  new;
-        source  boxToCell;
-        box     (2.4 2.4 0) (3 3 1);
-    }
-
-    {
-        name    blockedCellsSet;
-        type    cellSet;
-        action  add;
-        source  boxToCell;
-        box     (0 0 0) (0.6 0.6 1);
-    }
-
-    {
-        name    blockedFacesSet;
-        type    faceSet;
-        action  add;
-        source  cellToFace;
-        set     blockedCellsSet;
-        option  all;
-    }
-
-    {
-        name    coupledFacesSet;
-        type    faceSet;
-        action  new;
-        source  boxToFace;
-        box     (1.0 1.99 0) (2.0 2.09 0.6);
-    }
-
-);
-
-// ************************************************************************* //

tutorials/combustion/PDRFoam/pipeLattice/Allrun

@@ -12,6 +12,8 @@ runApplication PDRblockMesh
 
 runApplication PDRsetFields
 
+runApplication PDRMesh
+
 # Get location for ignition
 
 echo "Ignition point from constant/combustionProperties > system/ignitionPoint"

tutorials/combustion/PDRFoam/pipeLattice/constant/PDRProperties

@@ -15,19 +15,18 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+smoothSt        false;
 StSmoothCoef    1000.0;
 
-//No smooth if >100
-
 schelkin
 {
-    subGridSchelkin true;
-    uPrimeCoef  1.0;
-    lCoef       0.336;
-    maxSchFac   100.0;
-    nrCoef      0.1;
-    nrExp       0.0;
-    nrExp2      0.05;
+    subGridSchelkin     true;
+    uPrimeCoef          1.0;
+    lCoef               0.336;
+    maxSchFac           100.0;
+    nrCoef              0.1;
+    nrExp               0.0;
+    nrExp2              0.05;
 }
 
 PDRDragModel    basic;
@@ -47,17 +46,8 @@ basicSchCoeffs
     Csk         0.0;
 }
 
-XiModel         transport;
+XiModel    transportTwoEqs;
 
-algebraicCoeffs
-{
-    XiShapeCoef 1;
-    XpShapeCoef 1;
-    CpfiDot     0.0;
-    CpfiCross   0.0;
-    GEtaExp     0.0;
-    LOverCw     0.01;
-}
 
 transportTwoEqsCoeffs
 {
@@ -69,51 +59,11 @@ transportTwoEqsCoeffs
     LOverCw     0.01;
 }
 
-k3Coeffs
-{
-    k3Obs       10.0;
-    k3Open      1.0;
-}
-
-transportThreeEqsCoeffs
-{
-    $transportTwoEqsCoeffs;
-    $k3Coeffs;
-}
-
-transportFourEqsCoeffs
-{
-    $transportTwoEqsCoeffs;
-    $k3Coeffs;
-}
-
-transportOneEqObsCoeffs
-{
-    XiShapeCoef 1;
-    Cpfi        0.0;
-    GEtaExp     0.0;
-}
-
-transportXp
-{
-}
-
-transportCoeffs
-{
-    XiShapeCoef 1;
-    GEtaExp     0.28;
-}
-
-fixedCoeffs
-{
-}
-
-
 /*---------------------------------------------------------------------------*\
                           XiEqModel : Model for XiEq
 \*---------------------------------------------------------------------------*/
 
-XiEqModel       instability;
+XiEqModel       instability2XiEq;
 
 BLMcoeffs
 {
@@ -141,12 +91,9 @@ instability2XiEqCoeffs
     BLMgMaXiEqCoeffs
     {
         $schelkin;
-        gulderCoef  1.0; //this value is not usssed 1.0.
         kaCoef      0.25;
         lowK0       0.1;
         lowKg       0.0;
-        gMaCoef     0.032; //not used
-        gMaCoef1    0.0;  // not used
         $BLMcoeffs;
     }
 
@@ -161,45 +108,6 @@ instability2XiEqCoeffs
      }
 }
 
-instabilityCoeffs
-{
-
-    XiEqIn      2.5;
-    XiEqModel   Gulder;
-    GulderCoeffs
-    {
-        uPrimeCoef      1.0;
-        subGridSchelkin true;
-        XiEqCoef        0.62;
-    }
-
-    SCOPEBlendCoeffs
-    {
-        XiEqModelL
-        {
-            XiEqModel   Gulder;
-
-            GulderCoeffs
-            {
-                $schelkin;
-                XiEqCoef    0.62;
-            }
-        }
-
-        XiEqModelH
-        {
-            XiEqModel       SCOPEXiEq;
-
-            SCOPEXiEqCoeffs
-            {
-                $schelkin;
-                XiEqCoef    1.6;
-                XiEqExp     0.33333;
-            }
-        }
-    }
-}
-
 
 /*---------------------------------------------------------------------------*\
                      XiGModel : Model for generation of Xi
@@ -207,25 +115,10 @@ instabilityCoeffs
 
 XiGModel            instabilityG;
 
-instability2GCoeffs
-{
-    lambdaIn        0.0001;
-    defaultCIn      10.0;
-    GInMult         5.0;
-    GInFade         4.0;
-
-    XiGModel KTS;
-
-    KTSCoeffs
-    {
-        GEtaCoef    0.28;
-    }
-}
-
 instabilityGCoeffs
 {
-    lambdaIn        4.5e-3;
-    GIn             1.917;
+    lambdaIn        lambdaIn   [0 1 0 0 0 0 0] 4.5e-3;
+    GIn             GIn        [0 0 -1 0 0 0 0] 1.917;
 
     XiGModel        KTS;
     KTSCoeffs
@@ -249,11 +142,6 @@ normBasicSubGridCoeffs
     Cxpe4       1.0;
 }
 
-basicSubGridCoeffs
-{
-}
-
-
 /*---------------------------------------------------------------------------*\
                      XpGModel : Model for generation of Xp
 \*---------------------------------------------------------------------------*/
@@ -273,14 +161,4 @@ normBasicSubGridGCoeffs
     Cxpe3       400.0;
     Cxpe4       1.0;
 }
-
-basicSubGridGCoeffs
-{
-    XpGModel    KTS;
-    KTSCoeffs
-    {
-        GEtaCoef    0.28;
-    }
-}
-
 // ************************************************************************* //

tutorials/combustion/PDRFoam/pipeLattice/constant/thermophysicalProperties

@@ -15,15 +15,11 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-//thermoType  hhuMixtureThermo<inhomogeneousMixture<sutherlandTransport<specieThermo<janafThermo<perfectGas>>>>>;
 
 thermoType
 {
     type            heheuPsiThermo;
     mixture         inhomogeneousMixture;
-    //mixture         sprayHomogeneousMixtureDisp;
-    //mixture         sprayInhomogeneousMixture;
-    //mixture         inhomogeneousMixture;
     transport       sutherland;
     thermo          janaf;
     equationOfState perfectGas;
@@ -31,26 +27,6 @@ thermoType
     energy          absoluteEnthalpy;
 }
 
-inertSpecie     AIR;
-
-liquids
-{
-    //liquidComponents ( H2OPS );
-    H2OPS
-    {
-        defaultCoeffs   yes;
-    }
-    water
-    {
-        defaultCoeffs   yes;
-    }
-}
-
-solids
-{
-    //solidComponents ( );
-}
-
 stoichiometricAirFuelMassRatio 17.1667;
 
 fuel
@@ -119,27 +95,7 @@ burntProducts
     }
 }
 
-steam
-{
-    specie
-    {
-        nMoles          1;
-        molWeight       18.02;
-    }
-    thermodynamics
-    {
-        Tlow            200.00;
-        Thigh           5000.00;
-        Tcommon         1000.00;
-        highCpCoeffs    ( 2.67215e+00 3.05629e-03 -8.73026e-07 1.20100e-10 -6.39162e-15 -29899.199219 6.862820);
-        lowCpCoeffs     ( 3.38684e+00 3.47498e-03 -6.35470e-06 6.96858e-09 -2.50659e-12 -30208.099609 2.590230);
-    }
-    transport
-    {
-        As      1.6721e-06;
-        Ts      1.7067e+02;
-    }
-}
-
+dpdt        false;
+pMin        1e3;
 
 // ************************************************************************* //