ENH: InterfaceCompositionModels: add new thermo types for kineticGasEvaporation

More consistent use of combineReduce, simpler and/or reductions

See merge request Development/openfoam!566

applications/solvers/combustion/PDRFoam/setDeltaT.H

@@ -36,11 +36,11 @@ Description
 if (adjustTimeStep)
 {
     scalar maxDeltaTFact = maxCo/(CoNum + StCoNum + SMALL);
-    scalar deltaTFact = Foam::min(Foam::min(maxDeltaTFact, 1.0 + 0.1*maxDeltaTFact), 1.2);
+    scalar deltaTFact = min(min(maxDeltaTFact, 1.0 + 0.1*maxDeltaTFact), 1.2);
 
     runTime.setDeltaT
     (
-        Foam::min
+        min
         (
             deltaTFact*runTime.deltaTValue(),
             maxDeltaT

applications/solvers/combustion/chemFoam/setDeltaT.H

@@ -1,5 +1,5 @@
 if (adjustTimeStep)
 {
-    runTime.setDeltaT(Foam::min(dtChem, maxDeltaT));
+    runTime.setDeltaT(min(dtChem, maxDeltaT));
     Info<< "deltaT = " <<  runTime.deltaT().value() << endl;
 }

applications/solvers/combustion/fireFoam/setMultiRegionDeltaT.H

@@ -54,9 +54,9 @@ if (adjustTimeStep)
 
     runTime.setDeltaT
     (
-        Foam::min
+        min
         (
-            dt0*Foam::min(Foam::min(TFactorFluid, Foam::min(TFactorFilm, TFactorSolid)), 1.2),
+            dt0*min(min(TFactorFluid, min(TFactorFilm, TFactorSolid)), 1.2),
             maxDeltaT
         )
     );

applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionTwoPhaseEulerFoam/fluid/compressibleMultiRegionCourantNo.H

@@ -18,13 +18,13 @@
             const surfaceScalarField& phi2 =
                 phaseSystemFluid[regioni].phase2().phiRef();
 
-            sumPhi = Foam::max
+            sumPhi = max
             (
                 sumPhi,
                 fvc::surfaceSum(mag(phi1))().primitiveField()
             );
 
-            sumPhi = Foam::max
+            sumPhi = max
             (
                 sumPhi,
                 fvc::surfaceSum(mag(phi2))().primitiveField()
@@ -43,7 +43,7 @@
               / fluidRegions[regioni].V().field()
             )*runTime.deltaTValue(),
 
-            CoNum = Foam::max(UrCoNum, CoNum);
+            CoNum = max(UrCoNum, CoNum);
         }
     }
 

applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleMultiRegionCourantNo.H

@@ -2,7 +2,7 @@
 
     forAll(fluidRegions, regioni)
     {
-        CoNum = Foam::max
+        CoNum = max
         (
             compressibleCourantNo
             (

applications/solvers/heatTransfer/chtMultiRegionFoam/include/setInitialMultiRegionDeltaT.H

@@ -47,10 +47,10 @@ if (adjustTimeStep)
 
         runTime.setDeltaT
         (
-            Foam::min
+            min
             (
-                Foam::min(maxCo/CoNum, maxDi/DiNum)*runTime.deltaT().value(),
-                Foam::min(runTime.deltaTValue(), maxDeltaT)
+                min(maxCo/CoNum, maxDi/DiNum)*runTime.deltaT().value(),
+                min(runTime.deltaTValue(), maxDeltaT)
             )
         );
         Info<< "deltaT = " <<  runTime.deltaT().value() << endl;

applications/solvers/heatTransfer/chtMultiRegionFoam/include/setMultiRegionDeltaT.H

@@ -49,17 +49,17 @@ if (adjustTimeStep)
     scalar maxDeltaTSolid = maxDi/(DiNum + SMALL);
 
     scalar deltaTFluid =
-        Foam::min
+        min
         (
-            Foam::min(maxDeltaTFluid, 1.0 + 0.1*maxDeltaTFluid),
+            min(maxDeltaTFluid, 1.0 + 0.1*maxDeltaTFluid),
             1.2
         );
 
     runTime.setDeltaT
     (
-        Foam::min
+        min
         (
-            Foam::min(deltaTFluid, maxDeltaTSolid)*runTime.deltaT().value(),
+            min(deltaTFluid, maxDeltaTSolid)*runTime.deltaT().value(),
             maxDeltaT
         )
     );

applications/solvers/heatTransfer/chtMultiRegionFoam/solid/solidRegionDiffusionNo.H

@@ -22,7 +22,7 @@ forAll(solidRegions, i)
     tmp<volScalarField> trho = thermo.rho();
     const volScalarField& rho = trho();
 
-    DiNum = Foam::max
+    DiNum = max
     (
         solidRegionDiffNo
         (

applications/solvers/heatTransfer/solidFoam/solidDiffusionNo.H

@@ -17,7 +17,7 @@ scalar DiNum = -GREAT;
     tmp<volScalarField> trho = thermo.rho();
     const volScalarField& rho = trho();
 
-    DiNum = Foam::max
+    DiNum = max
     (
         solidRegionDiffNo
         (

applications/solvers/lagrangian/reactingParcelFoam/setMultiRegionDeltaT.H

@@ -36,13 +36,13 @@ Description
 if (adjustTimeStep)
 {
     const scalar maxDeltaTFact =
-        Foam::min(maxCo/(CoNum + SMALL), maxCo/(surfaceFilm.CourantNumber() + SMALL));
+        min(maxCo/(CoNum + SMALL), maxCo/(surfaceFilm.CourantNumber() + SMALL));
     const scalar deltaTFact =
-        Foam::min(Foam::min(maxDeltaTFact, 1.0 + 0.1*maxDeltaTFact), 1.2);
+        min(min(maxDeltaTFact, 1.0 + 0.1*maxDeltaTFact), 1.2);
 
     runTime.setDeltaT
     (
-        Foam::min
+        min
         (
             deltaTFact*runTime.deltaTValue(),
             maxDeltaT

applications/solvers/multiphase/MPPICInterFoam/alphaEqn.H

@@ -214,7 +214,7 @@
                 phiCN,
                 alphaPhi10,
                 Sp,
-                (Su + divU*Foam::min(alpha1(), scalar(1)))(),
+                (Su + divU*min(alpha1(), scalar(1)))(),
                 oneField(),
                 zeroField()
             );

applications/solvers/multiphase/VoF/alphaEqn.H

@@ -214,7 +214,7 @@
                 phiCN,
                 alphaPhi10,
                 Sp,
-                (Su + divU*Foam::min(alpha1(), scalar(1)))(),
+                (Su + divU*min(alpha1(), scalar(1)))(),
                 oneField(),
                 zeroField()
             );

applications/solvers/multiphase/VoF/setDeltaT.H

@@ -36,13 +36,13 @@ Description
 if (adjustTimeStep)
 {
     scalar maxDeltaTFact =
-        Foam::min(maxCo/(CoNum + SMALL), maxAlphaCo/(alphaCoNum + SMALL));
+        min(maxCo/(CoNum + SMALL), maxAlphaCo/(alphaCoNum + SMALL));
 
-    scalar deltaTFact = Foam::min(Foam::min(maxDeltaTFact, 1.0 + 0.1*maxDeltaTFact), 1.2);
+    scalar deltaTFact = min(min(maxDeltaTFact, 1.0 + 0.1*maxDeltaTFact), 1.2);
 
     runTime.setDeltaT
     (
-        Foam::min
+        min
         (
             deltaTFact*runTime.deltaTValue(),
             maxDeltaT

applications/solvers/multiphase/cavitatingFoam/setDeltaT.H

@@ -36,13 +36,13 @@ Description
 if (adjustTimeStep)
 {
     scalar maxDeltaTFact =
-        Foam::min(maxCo/(CoNum + SMALL), maxAcousticCo/(acousticCoNum + SMALL));
+        min(maxCo/(CoNum + SMALL), maxAcousticCo/(acousticCoNum + SMALL));
 
-    scalar deltaTFact = Foam::min(Foam::min(maxDeltaTFact, 1.0 + 0.1*maxDeltaTFact), 1.2);
+    scalar deltaTFact = min(min(maxDeltaTFact, 1.0 + 0.1*maxDeltaTFact), 1.2);
 
     runTime.setDeltaT
     (
-        Foam::min
+        min
         (
             deltaTFact*runTime.deltaTValue(),
             maxDeltaT

applications/solvers/multiphase/cavitatingFoam/setInitialDeltaT.H

@@ -37,11 +37,11 @@ if (adjustTimeStep)
     if (CoNum > SMALL)
     {
         scalar maxDeltaTFact =
-            Foam::min(maxCo/(CoNum + SMALL), maxAcousticCo/(acousticCoNum + SMALL));
+            min(maxCo/(CoNum + SMALL), maxAcousticCo/(acousticCoNum + SMALL));
 
         runTime.setDeltaT
         (
-            Foam::min
+            min
             (
                 maxDeltaTFact*runTime.deltaTValue(),
                 maxDeltaT

applications/solvers/multiphase/compressibleInterFoam/alphaSuSp.H

@@ -26,12 +26,12 @@ forAll(dgdt, celli)
 {
     if (dgdt[celli] > 0.0)
     {
-        Sp[celli] -= dgdt[celli]/Foam::max(1.0 - alpha1[celli], 1e-4);
-        Su[celli] += dgdt[celli]/Foam::max(1.0 - alpha1[celli], 1e-4);
+        Sp[celli] -= dgdt[celli]/max(1.0 - alpha1[celli], 1e-4);
+        Su[celli] += dgdt[celli]/max(1.0 - alpha1[celli], 1e-4);
     }
     else if (dgdt[celli] < 0.0)
     {
-        Sp[celli] += dgdt[celli]/Foam::max(alpha1[celli], 1e-4);
+        Sp[celli] += dgdt[celli]/max(alpha1[celli], 1e-4);
     }
 }
 

applications/solvers/multiphase/compressibleInterFoam/compressibleInterIsoFoam/alphaEqn.H

@@ -1,6 +1,6 @@
 // Update alpha1
 #include "alphaSuSp.H"
-advector.advect(Sp,(Su + divU*Foam::min(alpha1(), scalar(1)))());
+advector.advect(Sp,(Su + divU*min(alpha1(), scalar(1)))());
 
 // Update rhoPhi
 rhoPhi = advector.getRhoPhi(rho1, rho2);
@@ -10,6 +10,6 @@ alpha2 = 1.0 - alpha1;
 
 Info<< "Phase-1 volume fraction = "
     << alpha1.weightedAverage(mesh.Vsc()).value()
-    << "  Min(" << alpha1.name() << ") = " << Foam::min(alpha1).value()
-    << "  Max(" << alpha1.name() << ") - 1 = " << Foam::max(alpha1).value() - 1
+    << "  Min(" << alpha1.name() << ") = " << min(alpha1).value()
+    << "  Max(" << alpha1.name() << ") - 1 = " << max(alpha1).value() - 1
     << endl;

applications/solvers/multiphase/compressibleInterFoam/compressibleInterIsoFoam/alphaSuSp.H

@@ -26,12 +26,12 @@ forAll(dgdt, celli)
 {
     if (dgdt[celli] > 0.0)
     {
-        Sp[celli] -= dgdt[celli]/Foam::max(1.0 - alpha1[celli], 1e-4);
-        Su[celli] += dgdt[celli]/Foam::max(1.0 - alpha1[celli], 1e-4);
+        Sp[celli] -= dgdt[celli]/max(1.0 - alpha1[celli], 1e-4);
+        Su[celli] += dgdt[celli]/max(1.0 - alpha1[celli], 1e-4);
     }
     else if (dgdt[celli] < 0.0)
     {
-        Sp[celli] += dgdt[celli]/Foam::max(alpha1[celli], 1e-4);
+        Sp[celli] += dgdt[celli]/max(alpha1[celli], 1e-4);
     }
 }
 

applications/solvers/multiphase/compressibleInterFoam/overCompressibleInterDyMFoam/alphaSuSp.H

@@ -26,12 +26,12 @@ forAll(dgdt, celli)
 {
     if (dgdt[celli] > 0.0)
     {
-        Sp[celli] -= dgdt[celli]/Foam::max(1.0 - alpha1[celli], 1e-4);
-        Su[celli] += dgdt[celli]/Foam::max(1.0 - alpha1[celli], 1e-4);
+        Sp[celli] -= dgdt[celli]/max(1.0 - alpha1[celli], 1e-4);
+        Su[celli] += dgdt[celli]/max(1.0 - alpha1[celli], 1e-4);
     }
     else if (dgdt[celli] < 0.0)
     {
-        Sp[celli] += dgdt[celli]/Foam::max(alpha1[celli], 1e-4);
+        Sp[celli] += dgdt[celli]/max(alpha1[celli], 1e-4);
     }
 }
 

applications/solvers/multiphase/icoReactingMultiphaseInterFoam/setDeltaT.H

@@ -36,13 +36,13 @@ Description
 if (adjustTimeStep)
 {
     scalar maxDeltaTFact =
-        Foam::min
+        min
         (
             maxCo/(CoNum + SMALL),
-            Foam::min
+            min
             (
                 maxAlphaCo/(alphaCoNum + SMALL),
-                Foam::min
+                min
                 (
                     maxAlphaDdt/(ddtAlphaNum + SMALL),
                     maxDi/(DiNum + SMALL)
@@ -50,11 +50,11 @@ if (adjustTimeStep)
             )
         );
 
-    scalar deltaTFact = Foam::min(Foam::min(maxDeltaTFact, 1.0 + 0.1*maxDeltaTFact), 1.2);
+    scalar deltaTFact = min(min(maxDeltaTFact, 1.0 + 0.1*maxDeltaTFact), 1.2);
 
     runTime.setDeltaT
     (
-        Foam::min
+        min
         (
             deltaTFact*runTime.deltaTValue(),
             maxDeltaT

applications/solvers/multiphase/reactingTwoPhaseEulerFoam/CourantNos.H

@@ -8,5 +8,5 @@
 
     Info<< "Max Ur Courant Number = " << UrCoNum << endl;
 
-    CoNum = Foam::max(CoNum, UrCoNum);
+    CoNum = max(CoNum, UrCoNum);
 }

applications/solvers/multiphase/twoPhaseEulerFoam/CourantNos.H

@@ -8,5 +8,5 @@
 
     Info<< "Max Ur Courant Number = " << UrCoNum << endl;
 
-    CoNum = Foam::max(CoNum, UrCoNum);
+    CoNum = max(CoNum, UrCoNum);
 }

applications/utilities/mesh/advanced/snappyRefineMesh/snappyRefineMesh.C

@@ -142,26 +142,26 @@ scalar getEdgeStats(const primitiveMesh& mesh, const direction excludeCmpt)
 
         if (mag(eVec & x) > 1-edgeTol)
         {
-            minX = Foam::min(minX, eMag);
-            maxX = Foam::max(maxX, eMag);
+            minX = min(minX, eMag);
+            maxX = max(maxX, eMag);
             nX++;
         }
         else if (mag(eVec & y) > 1-edgeTol)
         {
-            minY = Foam::min(minY, eMag);
-            maxY = Foam::max(maxY, eMag);
+            minY = min(minY, eMag);
+            maxY = max(maxY, eMag);
             nY++;
         }
         else if (mag(eVec & z) > 1-edgeTol)
         {
-            minZ = Foam::min(minZ, eMag);
-            maxZ = Foam::max(maxZ, eMag);
+            minZ = min(minZ, eMag);
+            maxZ = max(maxZ, eMag);
             nZ++;
         }
         else
         {
-            minOther = Foam::min(minOther, eMag);
-            maxOther = Foam::max(maxOther, eMag);
+            minOther = min(minOther, eMag);
+            maxOther = max(maxOther, eMag);
         }
     }
 
@@ -179,19 +179,19 @@ scalar getEdgeStats(const primitiveMesh& mesh, const direction excludeCmpt)
 
     if (excludeCmpt == 0)
     {
-        return Foam::min(minY, Foam::min(minZ, minOther));
+        return min(minY, min(minZ, minOther));
     }
     else if (excludeCmpt == 1)
     {
-        return Foam::min(minX, Foam::min(minZ, minOther));
+        return min(minX, min(minZ, minOther));
     }
     else if (excludeCmpt == 2)
     {
-        return Foam::min(minX, Foam::min(minY, minOther));
+        return min(minX, min(minY, minOther));
     }
     else
     {
-        return Foam::min(minX, Foam::min(minY, Foam::min(minZ, minOther)));
+        return min(minX, min(minY, min(minZ, minOther)));
     }
 }
 
@@ -771,7 +771,7 @@ int main(int argc, char *argv[])
     {
         Info<< "Read existing refinement level from file "
             << refLevel.objectPath() << nl
-            << "   min level : " << Foam::min(refLevel) << nl
+            << "   min level : " << min(refLevel) << nl
             << "   max level : " << maxLevel << nl
             << endl;
     }

applications/utilities/mesh/conversion/cfx4ToFoam/cfx4ToFoam.C

@@ -271,7 +271,7 @@ int main(int argc, char *argv[])
                     if (blockPFacePointi != blockPFacePointi2)
                     {
                         sqrMergeTol =
-                            Foam::min
+                            min
                             (
                                 sqrMergeTol,
                                 magSqr
@@ -338,16 +338,16 @@ int main(int argc, char *argv[])
                         blockNFacePoints[blockNFacePointi]
                       + blockOffsets[blockNlabel];
 
-                    label minPN = Foam::min(PpointLabel, NpointLabel);
+                    label minPN = min(PpointLabel, NpointLabel);
 
                     if (pointMergeList[PpointLabel] != -1)
                     {
-                        minPN = Foam::min(minPN, pointMergeList[PpointLabel]);
+                        minPN = min(minPN, pointMergeList[PpointLabel]);
                     }
 
                     if (pointMergeList[NpointLabel] != -1)
                     {
-                        minPN = Foam::min(minPN, pointMergeList[NpointLabel]);
+                        minPN = min(minPN, pointMergeList[NpointLabel]);
                     }
 
                     pointMergeList[PpointLabel]
@@ -411,7 +411,7 @@ int main(int argc, char *argv[])
 
                         pointMergeList[PpointLabel]
                       = pointMergeList[NpointLabel]
-                      = Foam::min
+                      = min
                         (
                             pointMergeList[PpointLabel],
                             pointMergeList[NpointLabel]
@@ -757,7 +757,7 @@ int main(int argc, char *argv[])
     );
 
     // Set the precision of the points data to 10
-    IOstream::defaultPrecision(Foam::max(10u, IOstream::defaultPrecision()));
+    IOstream::defaultPrecision(max(10u, IOstream::defaultPrecision()));
 
     Info<< "Writing polyMesh" << endl;
     pShapeMesh.removeFiles();

applications/utilities/mesh/conversion/ensightToFoam/Make/files

@@ -0,0 +1,4 @@
+ensightMeshReader.C
+ensightToFoam.C
+
+EXE = $(FOAM_APPBIN)/ensightToFoam

applications/utilities/mesh/conversion/ensightToFoam/Make/options

@@ -0,0 +1,11 @@
+EXE_INC = \
+    -I$(LIB_SRC)/fileFormats/lnInclude \
+    -I$(LIB_SRC)/surfMesh/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/conversion/lnInclude
+
+EXE_LIBS = \
+    -lfileFormats \
+    -lsurfMesh \
+    -lmeshTools \
+    -lconversion

applications/utilities/mesh/conversion/ensightToFoam/ensightMeshReader.H

@@ -0,0 +1,194 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2022 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::fileFormats::ensightMeshReader
+
+Description
+
+Notes
+
+SourceFiles
+    ensightMeshReader.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_ensightMeshReader_H
+#define Foam_ensightMeshReader_H
+
+#include "meshReader.H"
+//#include "ensightReadFile.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// Forward Declarations
+class ensightReadFile;
+
+namespace fileFormats
+{
+
+/*---------------------------------------------------------------------------*\
+                Class fileFormats::ensightMeshReader Declaration
+\*---------------------------------------------------------------------------*/
+
+class ensightMeshReader
+:
+    public meshReader
+{
+    // Private Data
+
+        //- Merge distance
+        const scalar mergeTol_;
+
+        //- Check and correct handedness
+        const bool setHandedness_;
+
+
+protected:
+
+    // Protected Data
+
+        //- mesh point to original node_id
+        labelList nodeIds_;
+
+        //- mesh cell to original element_id
+        labelList elementIds_;
+
+
+    // Protected Member Functions
+
+        //- Rotate face so lowest vertex is first
+        const face& rotateFace
+        (
+            const face& f,
+            face& rotatedFace
+        ) const;
+
+        //- Read set of vertices. Optional mapping
+        void readVerts
+        (
+            ensightReadFile& is,
+            const label nVerts,
+            const Map<label>& nodeIdToPoints,
+            DynamicList<label>& verts
+        ) const;
+
+        //- Read set of element/node IDs
+        void readIDs
+        (
+            ensightReadFile& is,
+            const bool doRead,
+            const label nShapes,
+            labelList& foamToElem,
+            Map<label>& elemToFoam
+        ) const;
+
+        //- Swap handedness of hex if needed
+        void setHandedness
+        (
+            const cellModel& model,
+            DynamicList<label>& verts,
+            const pointField& points
+        ) const;
+
+        //- Read a single part until eof (return true) or until start of next
+        //  part (return false)
+        bool readGoldPart
+        (
+            ensightReadFile& is,
+            const bool read_node_ids,
+            const bool read_elem_ids,
+
+            pointField& points,
+            labelList& pointToNodeIds,
+            Map<label>& nodeIdToPoints,
+
+            // 3D-elems : cells (cell-to-faces)
+            faceListList& cellFaces,
+            labelList& cellToElemIds,
+            Map<label>& elemIdToCells,
+
+            // 2D-elems : faces
+            faceList& faces,
+            labelList& faceToElemIDs,
+            Map<label>& elemIdToFaces
+        ) const;
+
+        //- Read the mesh from the file(s)
+        virtual bool readGeometry(const scalar scaleFactor = 1.0);
+
+
+public:
+
+    //- Runtime type information
+    TypeName("ensightMeshReader");
+
+
+    // Constructors
+
+        //- Construct from case name
+        ensightMeshReader
+        (
+            const fileName& geomFile,
+            const objectRegistry& registry,
+            const scalar mergeTol = SMALL,
+            const scalar scaleFactor = 1.0,
+            const bool setHandedness = true
+        );
+
+
+    //- Destructor
+    virtual ~ensightMeshReader() = default;
+
+
+    // Access
+
+        //- Original node id (if supplied) or -1
+        const labelList& nodeIds() const noexcept
+        {
+            return nodeIds_;
+        }
+
+        //- Original element id (if supplied) or -1
+        const labelList& elementIds() const noexcept
+        {
+            return elementIds_;
+        }
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace fileFormats
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

applications/utilities/mesh/conversion/ensightToFoam/ensightToFoam.C

@@ -0,0 +1,118 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2022 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    ensightToFoam
+
+Group
+    grpMeshConversionUtilities
+
+Description
+    Convert an Ensight Gold mesh into OpenFOAM format.
+
+Usage
+    \b ensightToFoam [OPTION] \<ensightGeometryFile\>
+
+    Options:
+      - \par -mergeTol \<factor\>
+        Specify an alternative merging tolerance as a fraction of
+        the bounding box of the points.
+
+      - \par -scale \<factor\>
+        Specify an optional geometry scaling factor.
+
+      - \par -keepHandedness
+        Do not automatically flip negative volume cells
+
+See also
+    Foam::meshReader and Foam::fileFormats::STARCDMeshReader
+
+\*---------------------------------------------------------------------------*/
+
+#include "argList.H"
+#include "Time.H"
+#include "ensightMeshReader.H"
+
+using namespace Foam;
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    argList::addNote
+    (
+        "Convert Ensight mesh to OpenFOAM"
+    );
+
+    argList::noParallel();
+    argList::addArgument(".geo file", "The file containing the geometry");
+    argList::addOption
+    (
+        "mergeTol",
+        "factor",
+        "Merge tolerance as a fraction of bounding box - 0 to disable merging"
+    );
+    argList::addOption
+    (
+        "scale",
+        "factor",
+        "Geometry scaling factor - default is 1"
+    );
+    argList::addBoolOption
+    (
+        "keepHandedness",
+        "Do not automatically flip inverted cells"
+        " (default is to do a geometric test)"
+    );
+
+    argList args(argc, argv);
+    Time runTime(args.rootPath(), args.caseName());
+
+    // Increase the precision of the points data
+    IOstream::defaultPrecision(max(10u, IOstream::defaultPrecision()));
+
+    const fileName geomFile(args.get<fileName>(1));
+
+    {
+        fileFormats::ensightMeshReader reader
+        (
+            geomFile,
+            runTime,
+            args.getOrDefault<scalar>("mergeTol", 1e-10),
+            args.getOrDefault<scalar>("scale", 1.0),
+            args.found("keepHandedness")
+        );
+
+        autoPtr<polyMesh> mesh = reader.mesh(runTime);
+        mesh().setInstance(runTime.constant());
+        mesh().write();
+    }
+
+    Info<< "\nEnd\n" << endl;
+
+    return 0;
+}
+
+// ************************************************************************* //

applications/utilities/mesh/conversion/gambitToFoam/gambitToFoam.L

@@ -350,7 +350,7 @@ mtype                 {space}"MTYPE:"{space}
         // Find out how many labels are expected. If less or equal to
         // seven, read them all and finish with it. If there is more,
         // set read of the next line
-        label labelsToRead = Foam::min(8, nVertices);
+        label labelsToRead = min(8, nVertices);
         label labelI = 0;
         for (; labelI < labelsToRead; labelI++)
         {
@@ -387,7 +387,7 @@ mtype                 {space}"MTYPE:"{space}
 
         labelList& curLabels = cellLabels[curNumberOfCells];
 
-        label labelsToRead = Foam::min
+        label labelsToRead = min
         (
             (nCellContinuationLines + 1)*7,
             curLabels.size()
@@ -869,7 +869,7 @@ int main(int argc, char *argv[])
     );
 
     // Set the precision of the points data to 10
-    IOstream::defaultPrecision(Foam::max(10u, IOstream::defaultPrecision()));
+    IOstream::defaultPrecision(max(10u, IOstream::defaultPrecision()));
 
     Info<< "Writing polyMesh" << endl;
     pShapeMesh.removeFiles();

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

@@ -269,7 +269,7 @@ void readCells
     label maxUnvPoint = 0;
     forAll(unvPointID, pointi)
     {
-        maxUnvPoint = Foam::max(maxUnvPoint, unvPointID[pointi]);
+        maxUnvPoint = max(maxUnvPoint, unvPointID[pointi]);
     }
     labelList unvToFoam(invert(maxUnvPoint+1, unvPointID));
 
@@ -784,7 +784,7 @@ int main(int argc, char *argv[])
     label maxUnvPoint = 0;
     forAll(unvPointID, pointi)
     {
-        maxUnvPoint = Foam::max(maxUnvPoint, unvPointID[pointi]);
+        maxUnvPoint = max(maxUnvPoint, unvPointID[pointi]);
     }
     labelList unvToFoam(invert(maxUnvPoint+1, unvPointID));
 
@@ -1281,7 +1281,7 @@ int main(int argc, char *argv[])
     }
 
     // Set the precision of the points data to 10
-    IOstream::defaultPrecision(Foam::max(10u, IOstream::defaultPrecision()));
+    IOstream::defaultPrecision(max(10u, IOstream::defaultPrecision()));
 
     mesh.write();
 

applications/utilities/mesh/conversion/kivaToFoam/checkPatch.H

@@ -8,10 +8,10 @@
         {
             if (kivaVersion == kiva3v)
             {
-                regionIndex = Foam::max
+                regionIndex = max
                 (
-                    Foam::max(idface[quadFace[0]], idface[quadFace[1]]),
-                    Foam::max(idface[quadFace[2]], idface[quadFace[3]])
+                    max(idface[quadFace[0]], idface[quadFace[1]]),
+                    max(idface[quadFace[2]], idface[quadFace[3]])
                 );
 
                 if (regionIndex > 0)

applications/utilities/mesh/conversion/kivaToFoam/readKivaGrid.H

@@ -148,7 +148,7 @@ for (label i=0; i<nPoints; i++)
                     end = pointMap[end];
                 }
 
-                label minLabel = Foam::min(start, end);
+                label minLabel = min(start, end);
 
                 pointMap[start] = pointMap[end] = minLabel;
             }
@@ -331,7 +331,7 @@ if
     {
         forAll(pf, pfi)
         {
-            minz = Foam::min(minz, points[pf[pfi]].z());
+            minz = min(minz, points[pf[pfi]].z());
         }
     }
 
@@ -344,7 +344,7 @@ if
         scalar minfz = GREAT;
         forAll(pf, pfi)
         {
-            minfz = Foam::min(minfz, points[pf[pfi]].z());
+            minfz = min(minfz, points[pf[pfi]].z());
         }
 
         if (minfz > minz)
@@ -371,7 +371,7 @@ if
         scalar minfz = GREAT;
         forAll(pf, pfi)
         {
-            minfz = Foam::min(minfz, points[pf[pfi]].z());
+            minfz = min(minfz, points[pf[pfi]].z());
         }
 
         if (minfz < zHeadMin)
@@ -570,7 +570,7 @@ polyMesh pShapeMesh
 );
 
 // Set the precision of the points data to 10
-IOstream::defaultPrecision(Foam::max(10u, IOstream::defaultPrecision()));
+IOstream::defaultPrecision(max(10u, IOstream::defaultPrecision()));
 
 Info << "Writing polyMesh" << endl;
 pShapeMesh.removeFiles();

applications/utilities/mesh/conversion/netgenNeutralToFoam/netgenNeutralToFoam.C

@@ -188,7 +188,7 @@ int main(int argc, char *argv[])
         }
 
 
-        maxPatch = Foam::max(maxPatch, patchi);
+        maxPatch = max(maxPatch, patchi);
 
         triFace tri(readLabel(str)-1, readLabel(str)-1, readLabel(str)-1);
 
@@ -319,7 +319,7 @@ int main(int argc, char *argv[])
     );
 
     // Set the precision of the points data to 10
-    IOstream::defaultPrecision(Foam::max(10u, IOstream::defaultPrecision()));
+    IOstream::defaultPrecision(max(10u, IOstream::defaultPrecision()));
 
     Info<< "Writing mesh ..." << endl;
     mesh.removeFiles();

applications/utilities/mesh/generation/extrude/extrudeToRegionMesh/extrudeToRegionMesh.C

@@ -556,8 +556,8 @@ void calcEdgeMinMaxZone
             forAll(eFaces, i)
             {
                 label zoneI = mappedZoneID[eFaces[i]];
-                minZoneID[edgeI] = Foam::min(minZoneID[edgeI], zoneI);
-                maxZoneID[edgeI] = Foam::max(maxZoneID[edgeI], zoneI);
+                minZoneID[edgeI] = min(minZoneID[edgeI], zoneI);
+                maxZoneID[edgeI] = max(maxZoneID[edgeI], zoneI);
             }
         }
     }
@@ -813,8 +813,8 @@ void addCoupledPatches
             forAll(eFaces, i)
             {
                 label proci = procID[eFaces[i]];
-                minProcID[edgeI] = Foam::min(minProcID[edgeI], proci);
-                maxProcID[edgeI] = Foam::max(maxProcID[edgeI], proci);
+                minProcID[edgeI] = min(minProcID[edgeI], proci);
+                maxProcID[edgeI] = max(maxProcID[edgeI], proci);
             }
         }
     }
@@ -1291,7 +1291,7 @@ void extrudeGeometricProperties
             label celli = regionMesh.faceOwner()[facei];
             if (regionMesh.isInternalFace(facei))
             {
-                celli = Foam::max(celli, regionMesh.faceNeighbour()[facei]);
+                celli = max(celli, regionMesh.faceNeighbour()[facei]);
             }
 
             // Calculate layer from cell numbering (see createShellMesh)
@@ -2192,8 +2192,8 @@ int main(int argc, char *argv[])
 
             if (zone0 != zone1) // || (cos(angle) > blabla))
             {
-                label minZone = Foam::min(zone0,zone1);
-                label maxZone = Foam::max(zone0,zone1);
+                label minZone = min(zone0,zone1);
+                label maxZone = max(zone0,zone1);
                 label index = minZone*zoneNames.size()+maxZone;
 
                 ePatches.setSize(eFaces.size());

applications/utilities/mesh/generation/foamyMesh/foamyHexMeshBackgroundMesh/foamyHexMeshBackgroundMesh.C

@@ -727,14 +727,12 @@ int main(int argc, char *argv[])
 
     pointField mergedPoints;
     faceList mergedFaces;
-    labelList pointMergeMap;
     PatchTools::gatherAndMerge
     (
         tolDim,
         primitivePatch(SubList<face>(isoFaces), isoPoints),
         mergedPoints,
-        mergedFaces,
-        pointMergeMap
+        mergedFaces
     );
 
     if (Pstream::master())

applications/utilities/mesh/manipulation/checkMesh/writeFields.C

@@ -27,7 +27,7 @@ void maxFaceToCell
         const cell& cFaces = cells[cellI];
         forAll(cFaces, i)
         {
-            cellFld[cellI] = Foam::max(cellFld[cellI], faceData[cFaces[i]]);
+            cellFld[cellI] = max(cellFld[cellI], faceData[cFaces[i]]);
         }
     }
 
@@ -57,7 +57,7 @@ void minFaceToCell
         const cell& cFaces = cells[cellI];
         forAll(cFaces, i)
         {
-            cellFld[cellI] = Foam::min(cellFld[cellI], faceData[cFaces[i]]);
+            cellFld[cellI] = min(cellFld[cellI], faceData[cFaces[i]]);
         }
     }
 
@@ -88,8 +88,8 @@ void minFaceToCell
     // Internal faces
     forAll(own, facei)
     {
-        cellFld[own[facei]] = Foam::min(cellFld[own[facei]], faceData[facei]);
-        cellFld[nei[facei]] = Foam::min(cellFld[nei[facei]], faceData[facei]);
+        cellFld[own[facei]] = min(cellFld[own[facei]], faceData[facei]);
+        cellFld[nei[facei]] = min(cellFld[nei[facei]], faceData[facei]);
     }
 
     // Patch faces
@@ -100,7 +100,7 @@ void minFaceToCell
 
         forAll(fc, i)
         {
-            cellFld[fc[i]] = Foam::min(cellFld[fc[i]], fvp[i]);
+            cellFld[fc[i]] = min(cellFld[fc[i]], fvp[i]);
         }
     }
 
@@ -180,7 +180,7 @@ void Foam::writeFields
         (
             radToDeg
             (
-                Foam::acos(Foam::min(scalar(1), Foam::max(scalar(-1), faceOrthogonality)))
+                Foam::acos(min(scalar(1), max(scalar(-1), faceOrthogonality)))
             )
         );
 
@@ -540,7 +540,7 @@ void Foam::writeFields
                         ownCc,
                         fc
                     ).quality();
-                    ownVol = Foam::min(ownVol, tetQual);
+                    ownVol = min(ownVol, tetQual);
                 }
             }
             if (mesh.isInternalFace(facei))
@@ -556,7 +556,7 @@ void Foam::writeFields
                         fc,
                         neiCc
                     ).quality();
-                    neiVol = Foam::min(neiVol, tetQual);
+                    neiVol = min(neiVol, tetQual);
                 }
             }
         }
@@ -608,8 +608,8 @@ void Foam::writeFields
         // Internal faces
         forAll(own, facei)
         {
-            cellFld[own[facei]] = Foam::min(cellFld[own[facei]], ownPyrVol[facei]);
-            cellFld[nei[facei]] = Foam::min(cellFld[nei[facei]], neiPyrVol[facei]);
+            cellFld[own[facei]] = min(cellFld[own[facei]], ownPyrVol[facei]);
+            cellFld[nei[facei]] = min(cellFld[nei[facei]], neiPyrVol[facei]);
         }
 
         // Patch faces
@@ -620,7 +620,7 @@ void Foam::writeFields
             forAll(fc, i)
             {
                 const label meshFacei = fvp.patch().start();
-                cellFld[fc[i]] = Foam::min(cellFld[fc[i]], ownPyrVol[meshFacei]);
+                cellFld[fc[i]] = min(cellFld[fc[i]], ownPyrVol[meshFacei]);
             }
         }
 
@@ -631,7 +631,7 @@ void Foam::writeFields
         if (writeFaceFields)
         {
             scalarField minFacePyrVol(neiPyrVol);
-            minFacePyrVol = Foam::min
+            minFacePyrVol = min
             (
                 minFacePyrVol,
                 SubField<scalar>(ownPyrVol, mesh.nInternalFaces())

applications/utilities/mesh/manipulation/refineMesh/refineMesh.C

@@ -100,26 +100,26 @@ void printEdgeStats(const polyMesh& mesh)
 
             if (mag(eVec & x) > 1-edgeTol)
             {
-                minX = Foam::min(minX, eMag);
-                maxX = Foam::max(maxX, eMag);
+                minX = min(minX, eMag);
+                maxX = max(maxX, eMag);
                 nX++;
             }
             else if (mag(eVec & y) > 1-edgeTol)
             {
-                minY = Foam::min(minY, eMag);
-                maxY = Foam::max(maxY, eMag);
+                minY = min(minY, eMag);
+                maxY = max(maxY, eMag);
                 nY++;
             }
             else if (mag(eVec & z) > 1-edgeTol)
             {
-                minZ = Foam::min(minZ, eMag);
-                maxZ = Foam::max(maxZ, eMag);
+                minZ = min(minZ, eMag);
+                maxZ = max(maxZ, eMag);
                 nZ++;
             }
             else
             {
-                minOther = Foam::min(minOther, eMag);
-                maxOther = Foam::max(maxOther, eMag);
+                minOther = min(minOther, eMag);
+                maxOther = max(maxOther, eMag);
             }
         }
     }

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

@@ -145,10 +145,10 @@ void getBand
 
         // Note: mag not necessary for correct (upper-triangular) ordering.
         label diff = nei-own;
-        cellBandwidth[nei] = Foam::max(cellBandwidth[nei], diff);
+        cellBandwidth[nei] = max(cellBandwidth[nei], diff);
     }
 
-    bandwidth = Foam::max(cellBandwidth);
+    bandwidth = max(cellBandwidth);
 
     // Do not use field algebra because of conversion label to scalar
     profile = 0.0;
@@ -340,7 +340,7 @@ labelList getRegionFaceOrder
     }
 
     // Do region interfaces
-    label nRegions = Foam::max(cellToRegion)+1;
+    label nRegions = max(cellToRegion)+1;
     {
         // Sort in increasing region
         SortableList<label> sortKey(mesh.nFaces(), labelMax);
@@ -353,8 +353,8 @@ labelList getRegionFaceOrder
             if (ownRegion != neiRegion)
             {
                 sortKey[facei] =
-                    Foam::min(ownRegion, neiRegion)*nRegions
-                   +Foam::max(ownRegion, neiRegion);
+                    min(ownRegion, neiRegion)*nRegions
+                   +max(ownRegion, neiRegion);
             }
         }
         sortKey.sort();
@@ -566,7 +566,7 @@ labelList regionRenumber
 
     labelList cellOrder(cellToRegion.size());
 
-    label nRegions = Foam::max(cellToRegion)+1;
+    label nRegions = max(cellToRegion)+1;
 
     labelListList regionToCells(invertOneToMany(nRegions, cellToRegion));
 

applications/utilities/mesh/manipulation/setSet/setSet.C

@@ -320,10 +320,10 @@ bool doCommand
     const globalMeshData& parData = mesh.globalData();
 
     label typSize =
-        Foam::max
+        max
         (
             parData.nTotalCells(),
-            Foam::max
+            max
             (
                 parData.nTotalFaces(),
                 parData.nTotalPoints()
@@ -375,7 +375,7 @@ bool doCommand
 
                 topoSet& currentSet = currentSetPtr();
                 // Presize it according to current mesh data.
-                currentSet.resize(Foam::max(currentSet.size(), typSize));
+                currentSet.resize(max(currentSet.size(), typSize));
             }
         }
 

applications/utilities/mesh/manipulation/splitMeshRegions/splitMeshRegions.C

@@ -286,8 +286,8 @@ void addToInterface
 {
     edge interface
     (
-        Foam::min(ownRegion, neiRegion),
-        Foam::max(ownRegion, neiRegion)
+        min(ownRegion, neiRegion),
+        max(ownRegion, neiRegion)
     );
 
     auto iter = regionsToSize.find(interface);
@@ -544,8 +544,8 @@ void getInterfaceSizes
 
             edge interface
             (
-                Foam::min(ownRegion, neiRegion),
-                Foam::max(ownRegion, neiRegion)
+                min(ownRegion, neiRegion),
+                max(ownRegion, neiRegion)
             );
 
             faceToInterface[facei] = regionsToInterface[interface][zoneID];
@@ -567,8 +567,8 @@ void getInterfaceSizes
 
             edge interface
             (
-                Foam::min(ownRegion, neiRegion),
-                Foam::max(ownRegion, neiRegion)
+                min(ownRegion, neiRegion),
+                max(ownRegion, neiRegion)
             );
 
             faceToInterface[facei] = regionsToInterface[interface][zoneID];

applications/utilities/mesh/manipulation/subsetMesh/subsetMesh.C

@@ -345,7 +345,7 @@ void subsetTopoSets
 
         Info<< "Subsetting " << set.type() << " " << set.name() << endl;
 
-        labelHashSet subset(2*Foam::min(set.size(), map.size()));
+        labelHashSet subset(2*min(set.size(), map.size()));
 
         forAll(map, i)
         {

applications/utilities/parallelProcessing/reconstructParMesh/reconstructParMesh.C

@@ -977,7 +977,7 @@ int main(int argc, char *argv[])
                         for
                         (
                             label addedI=next;
-                            addedI<Foam::min(proci+step, nProcs);
+                            addedI<min(proci+step, nProcs);
                             addedI++
                         )
                         {

applications/utilities/parallelProcessing/redistributePar/redistributePar.C

@@ -334,11 +334,11 @@ void printMeshData(const polyMesh& mesh)
                 << nBndFaces-nProcFaces << endl;
         }
 
-        maxProcCells = Foam::max(maxProcCells, nLocalCells);
+        maxProcCells = max(maxProcCells, nLocalCells);
         totProcFaces += nProcFaces;
         totProcPatches += nei.size();
-        maxProcFaces = Foam::max(maxProcFaces, nProcFaces);
-        maxProcPatches = Foam::max(maxProcPatches, nei.size());
+        maxProcFaces = max(maxProcFaces, nProcFaces);
+        maxProcPatches = max(maxProcPatches, nei.size());
     }
 
     // Summary stats

applications/utilities/postProcessing/lagrangian/particleTracks/particleTracks.C

@@ -131,7 +131,7 @@ int main(int argc, char *argv[])
     args.readIfPresent("format", setFormat);
 
     args.readIfPresent("stride", sampleFrequency);
-    sampleFrequency = Foam::max(1, sampleFrequency);  // sanity
+    sampleFrequency = max(1, sampleFrequency);  // sanity
 
     // Setup the writer
     auto writerPtr =
@@ -179,7 +179,7 @@ int main(int argc, char *argv[])
                     maxIds.resize(origProc+1, -1);
                 }
 
-                maxIds[origProc] = Foam::max(maxIds[origProc], origId);
+                maxIds[origProc] = max(maxIds[origProc], origId);
             }
         }
 

applications/utilities/postProcessing/lagrangian/steadyParticleTracks/createControls.H

@@ -15,10 +15,5 @@ wordRes acceptFields(propsDict.get<wordRes>("fields"));
 wordRes excludeFields;
 propsDict.readIfPresent("exclude", excludeFields);
 
-const dictionary formatOptions
-(
-    propsDict.subOrEmptyDict("formatOptions", keyType::LITERAL)
-);
-
 
 // ************************************************************************* //

applications/utilities/preProcessing/PDR/pdrFields/PDRarraysCalc.C

@@ -1056,7 +1056,7 @@ void calc_drag_etc
     const scalar expon =
     (
         br > 0.0
-      ? Foam::min(Foam::max((surr_br / br - 0.25) * 4.0 / 3.0, scalar(0)), scalar(1))
+      ? min(max((surr_br / br - 0.25) * 4.0 / 3.0, scalar(0)), scalar(1))
       : 0.0
     );
 
@@ -1114,16 +1114,16 @@ void Foam::PDRarrays::blockageSummary() const
                 totVolBlock += v_block(ijk) * pdrBlock.V(ijk);
                 totArea += surf(ijk);
 
-                totCount += Foam::max(0, obs_count(ijk));
+                totCount += max(0, obs_count(ijk));
 
-                totDrag.x() += Foam::max(0, drag_s(ijk).xx());
-                totDrag.y() += Foam::max(0, drag_s(ijk).yy());
-                totDrag.z() += Foam::max(0, drag_s(ijk).zz());
+                totDrag.x() += max(0, drag_s(ijk).xx());
+                totDrag.y() += max(0, drag_s(ijk).yy());
+                totDrag.z() += max(0, drag_s(ijk).zz());
 
                 for (direction cmpt=0; cmpt < vector::nComponents; ++cmpt)
                 {
-                    totBlock[cmpt] += Foam::max(0, area_block_s(ijk)[cmpt]);
-                    totBlock[cmpt] += Foam::max(0, area_block_r(ijk)[cmpt]);
+                    totBlock[cmpt] += max(0, area_block_s(ijk)[cmpt]);
+                    totBlock[cmpt] += max(0, area_block_r(ijk)[cmpt]);
                 }
             }
         }

applications/utilities/preProcessing/PDR/pdrFields/PDRutilsOverlap.C

@@ -302,7 +302,7 @@ void Foam::PDRutils::circle_overlap
                     scalar da = ac - 0.5 * (a1 + a2);
                     scalar db = bc - 0.5 * (b1 + b2);
                     scalar dc = std::hypot(da, db);
-                    scalar rat1 = Foam::min(Foam::max((dc / sqrt(area) - 0.3) * 1.4, 0), 1);
+                    scalar rat1 = min(max((dc / sqrt(area) - 0.3) * 1.4, 0), 1);
                     scalar drg0 = c_drag(ia,ib).xx();
                     scalar drg1 = c_drag(ia,ib).yy();
                     scalar drg = std::hypot(drg0, drg1);
@@ -449,8 +449,8 @@ scalar block_overlap
             {
                 PDRobstacle over;
 
-                over.pt = Foam::max(blk1.pt, blk2.pt);
-                over.span = Foam::min(max1, max2) - over.pt;
+                over.pt = max(blk1.pt, blk2.pt);
+                over.span = min(max1, max2) - over.pt;
 
                 assert(cmptProduct(over.span) > 0.0);
 
@@ -603,11 +603,11 @@ scalar block_cylinder_overlap
 
                     over.x() = a_centre - 0.5 * a_lblk;
                     over.y() = b_centre - 0.5 * b_lblk;
-                    over.z() = Foam::max(blk1.z(), cyl2.z());
+                    over.z() = max(blk1.z(), cyl2.z());
 
                     over.span.x() = a_lblk;
                     over.span.y() = b_lblk;
-                    over.span.z() = Foam::min(max1.z(), cyl2.z() + cyl2.len()) - over.z();
+                    over.span.z() = min(max1.z(), cyl2.z() + cyl2.len()) - over.z();
                     assert(over.x() > -200.0);
                     assert(over.x() < 2000.0);
                 }
@@ -668,11 +668,11 @@ scalar block_cylinder_overlap
 
                     over.z() = a_centre - a_lblk * 0.5;
                     over.x() = b_centre - b_lblk * 0.5;
-                    over.y() = Foam::max(blk1.y(), cyl2.y());
+                    over.y() = max(blk1.y(), cyl2.y());
 
                     over.span.z() = a_lblk;
                     over.span.x() = b_lblk;
-                    over.span.y() = Foam::min(max1.y(), cyl2.y() + cyl2.len()) - over.y();
+                    over.span.y() = min(max1.y(), cyl2.y() + cyl2.len()) - over.y();
                 }
                 break;
 
@@ -734,11 +734,11 @@ scalar block_cylinder_overlap
 
                     over.y() = a_centre - a_lblk * 0.5;
                     over.z() = b_centre - b_lblk * 0.5;
-                    over.x() = Foam::max(blk1.x(), cyl2.x());
+                    over.x() = max(blk1.x(), cyl2.x());
 
                     over.span.y() = a_lblk;
                     over.span.z() = b_lblk;
-                    over.span.x() = Foam::min(max1.x(), cyl2.x() + cyl2.len()) - over.x();
+                    over.span.x() = min(max1.x(), cyl2.x() + cyl2.len()) - over.x();
                 }
                 break;
             }

applications/utilities/preProcessing/engineSwirl/engineSwirl.C

@@ -71,9 +71,9 @@ int main(int argc, char *argv[])
         {
             scalar b = j1(swirlProfile*r/cylinderRadius).value();
             scalar vEff = omega*b;
-            r = Foam::max(r, SMALL);
+            r = max(r, SMALL);
             U[celli] = ((vEff/r)*(c & yT))*xT + (-(vEff/r)*(c & xT))*yT;
-            Umax = Foam::max(Umax, mag(U[celli]));
+            Umax = max(Umax, mag(U[celli]));
         }
     }
 

applications/utilities/surface/surfaceCheck/surfaceCheck.C

@@ -889,7 +889,7 @@ int main(int argc, char *argv[])
             writeParts
             (
                 surf,
-                Foam::min(outputThreshold, numZones),
+                min(outputThreshold, numZones),
                 faceZone,
                 surfFilePath,
                 surfFileStem
@@ -955,7 +955,7 @@ int main(int argc, char *argv[])
             writeParts
             (
                 surf,
-                Foam::min(outputThreshold, numNormalZones),
+                min(outputThreshold, numNormalZones),
                 normalZone,
                 surfFilePath,
                 surfFileStem + "_normal"

applications/utilities/surface/surfaceHookUp/surfaceHookUp.C

@@ -290,7 +290,7 @@ int main(int argc, char *argv[])
     const IOdictionary dict(dictIO);
 
     const scalar dist(args.get<scalar>(1));
-    const scalar matchTolerance(Foam::max(1e-6*dist, SMALL));
+    const scalar matchTolerance(max(1e-6*dist, SMALL));
     const label maxIters = 100;
 
     Info<< "Hooking distance = " << dist << endl;

applications/utilities/surface/surfaceInflate/surfaceInflate.C

@@ -275,7 +275,7 @@ label detectIntersectionPoints
 
     // 1. Extrusion offset vectors intersecting new surface location
     {
-        scalar tol = Foam::max(tolerance, 10*s.tolerance());
+        scalar tol = max(tolerance, 10*s.tolerance());
 
         // Collect all the edge vectors. Slightly shorten the edges to prevent
         // finding lots of intersections. The fast triangle intersection routine
@@ -295,7 +295,7 @@ label detectIntersectionPoints
             &&  !localFaces[hits[pointI].index()].found(pointI)
             )
             {
-                scale[pointI] = Foam::max(0.0, scale[pointI]-0.2);
+                scale[pointI] = max(0.0, scale[pointI]-0.2);
 
                 isPointOnHitEdge.set(pointI);
                 nHits++;
@@ -328,7 +328,7 @@ label detectIntersectionPoints
                         << pt
                         << endl;
 
-                    scale[e[0]] = Foam::max(0.0, scale[e[0]]-0.2);
+                    scale[e[0]] = max(0.0, scale[e[0]]-0.2);
                     nHits++;
                 }
                 if (isPointOnHitEdge.set(e[1]))
@@ -340,7 +340,7 @@ label detectIntersectionPoints
                         << pt
                         << endl;
 
-                    scale[e[1]] = Foam::max(0.0, scale[e[1]]-0.2);
+                    scale[e[1]] = max(0.0, scale[e[1]]-0.2);
                     nHits++;
                 }
             }
@@ -416,7 +416,7 @@ void minSmooth
         const edge& e = edges[edgeI];
         scalar w = mag(points[mp[e[0]]]-points[mp[e[1]]]);
 
-        edgeWeights[edgeI] = 1.0/(Foam::max(w, SMALL));
+        edgeWeights[edgeI] = 1.0/(max(w, SMALL));
     }
 
     tmp<scalarField> tavgFld = avg(s, fld, edgeWeights);
@@ -427,7 +427,7 @@ void minSmooth
     {
         if (isAffectedPoint.test(pointI))
         {
-            newFld[pointI] = Foam::min
+            newFld[pointI] = min
             (
                 fld[pointI],
                 0.5*fld[pointI] + 0.5*avgFld[pointI]

applications/utilities/thermophysical/adiabaticFlameT/adiabaticFlameT.C

@@ -180,8 +180,8 @@ int main(int argc, char *argv[])
 
         scalar o2 = (1.0/equiv)*stoicO2;
         scalar n2 = (0.79/0.21)*o2;
-        scalar fres = Foam::max(1.0 - 1.0/equiv, 0.0);
-        scalar ores = Foam::max(1.0/equiv - 1.0, 0.0);
+        scalar fres = max(1.0 - 1.0/equiv, 0.0);
+        scalar ores = max(1.0/equiv - 1.0, 0.0);
         scalar fburnt = 1.0 - fres;
 
         thermo reactants

applications/utilities/thermophysical/equilibriumFlameT/equilibriumFlameT.C

@@ -196,12 +196,12 @@ int main(int argc, char *argv[])
         // Number of moles of species for one mole of fuel
         scalar o2 = (1.0/equiv)*stoicO2;
         scalar n2 = (0.79/0.21)*o2;
-        scalar fres = Foam::max(1.0 - 1.0/equiv, 0.0);
+        scalar fres = max(1.0 - 1.0/equiv, 0.0);
         scalar fburnt = 1.0 - fres;
 
         // Initial guess for number of moles of product species
         // ignoring product dissociation
-        scalar oresInit = Foam::max(1.0/equiv - 1.0, 0.0)*stoicO2;
+        scalar oresInit = max(1.0/equiv - 1.0, 0.0)*stoicO2;
         scalar co2Init = fburnt*stoicCO2;
         scalar h2oInit = fburnt*stoicH2O;
 
@@ -231,18 +231,18 @@ int main(int argc, char *argv[])
             if (j > 0)
             {
                 co = co2*
-                    Foam::min
+                    min
                     (
                         CO2BreakUp.Kn(P, equilibriumFlameTemperature, N)
-                       /::sqrt(Foam::max(ores, 0.001)),
+                       /::sqrt(max(ores, 0.001)),
                         1.0
                     );
 
                 h2 = h2o*
-                    Foam::min
+                    min
                     (
                         H2OBreakUp.Kn(P, equilibriumFlameTemperature, N)
-                       /::sqrt(Foam::max(ores, 0.001)),
+                       /::sqrt(max(ores, 0.001)),
                         1.0
                     );
 

src/OpenFOAM/containers/Lists/List/List.C

@@ -49,7 +49,7 @@ void Foam::List<T>::doResize(const label len)
         // With sign-check to avoid spurious -Walloc-size-larger-than
         T* nv = new T[len];
 
-        const label overlap = Foam::min(this->size_, len);
+        const label overlap = min(this->size_, len);
 
         if (overlap)
         {

src/OpenFOAM/db/IOstreams/Pstreams/UOPstream.H

@@ -327,7 +327,7 @@ public:
         //  \return True on success
         static bool write
         (
-            const commsTypes commsType,
+            const UPstream::commsTypes commsType,
             const int toProcNo,
             const char* buf,
             const std::streamsize bufSize,

src/OpenFOAM/db/IOstreams/Pstreams/UPstream.C

@@ -434,7 +434,7 @@ Foam::UPstream::commsTypes Foam::UPstream::defaultCommsType
 
 namespace Foam
 {
-    // Register re-reader
+    //- Registered reader for UPstream::defaultCommsType
     class addcommsTypeToOpt
     :
         public ::Foam::simpleRegIOobject

src/OpenFOAM/db/IOstreams/Pstreams/UPstream.H

@@ -66,12 +66,12 @@ public:
     //- Types of communications
     enum class commsTypes : char
     {
-        blocking,       //!< "blocking"
-        scheduled,      //!< "scheduled"
-        nonBlocking     //!< "nonBlocking"
+        blocking,       //!< "blocking"    : (MPI_Bsend, MPI_Recv)
+        scheduled,      //!< "scheduled"   : (MPI_Send, MPI_Recv)
+        nonBlocking     //!< "nonBlocking" : (MPI_Isend, MPI_Irecv)
     };
 
-    //- Names of the communication types
+    //- Enumerated names for the communication types
     static const Enum<commsTypes> commsTypeNames;
 
 

src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchTools.H

@@ -221,9 +221,48 @@ public:
     //
     //  \param[in] mergeDist Geometric merge tolerance for Foam::mergePoints
     //  \param[in] pp The patch to merge
-    //  \param[out] mergedPoints
-    //  \param[out] mergedFaces
-    //  \param[out] pointMergeMap
+    //  \param[out] mergedPoints merged points (master only, empty elsewhere)
+    //  \param[out] mergedFaces merged faces (master only, empty elsewhere)
+    //  \param[out] pointAddr Points globalIndex gather addressing
+    //      (master only, empty elsewhere)
+    //  \param[out] faceAddr Faces globalIndex gather addressing
+    //      (master only, empty elsewhere)
+    //  \param[out] pointMergeMap An old-to-new mapping from original
+    //      point index to the index into merged points.
+    //  \param[in]  useLocal  gather/merge patch localFaces/localPoints
+    //      instead of faces/points
+    //
+    //  \note
+    //  - OpenFOAM-v2112 and earlier: geometric merge on all patch points.
+    //  - OpenFOAM-v2206 and later: geometric merge on patch boundary points.
+    template<class FaceList, class PointField>
+    static void gatherAndMerge
+    (
+        const scalar mergeDist,
+        const PrimitivePatch<FaceList, PointField>& pp,
+        Field
+        <
+            typename PrimitivePatch<FaceList, PointField>::point_type
+        >& mergedPoints,
+        List
+        <
+            typename PrimitivePatch<FaceList, PointField>::face_type
+        >& mergedFaces,
+        globalIndex& pointAddr,
+        globalIndex& faceAddr,
+        labelList& pointMergeMap = const_cast<labelList&>(labelList::null()),
+        const bool useLocal = false
+    );
+
+    //- Gather points and faces onto master and merge into single patch.
+    //  Note: Normally uses faces/points (not localFaces/localPoints)
+    //
+    //  \param[in] mergeDist Geometric merge tolerance for Foam::mergePoints
+    //  \param[in] pp The patch to merge
+    //  \param[out] mergedPoints merged points (master only, empty elsewhere)
+    //  \param[out] mergedFaces merged faces (master only, empty elsewhere)
+    //  \param[out] pointMergeMap An old-to-new mapping from original
+    //      point index to the index into merged points.
     //  \param[in] useLocal gather/merge patch localFaces/localPoints
     //      instead of faces/points
     //
@@ -243,19 +282,20 @@ public:
         <
             typename PrimitivePatch<FaceList, PointField>::face_type
         >& mergedFaces,
-        labelList& pointMergeMap,
+        labelList& pointMergeMap = const_cast<labelList&>(labelList::null()),
         const bool useLocal = false
     );
 
     //- Gather (mesh!) points and faces onto master and merge collocated
     //  points into a single patch. Uses coupled point mesh
     //  structure so does not need tolerances.
-    //  On master and slave returns:
+    //  On master and sub-ranks returns:
     //  - pointToGlobal : for every local point index the global point index
     //  - uniqueMeshPointLabels : my local mesh points
     //  - globalPoints : global numbering for the global points
     //  - globalFaces : global numbering for the faces
-    //  On master only:
+    //  .
+    //  On master only returns:
     //  - mergedFaces : the merged faces
     //  - mergedPoints : the merged points
     template<class FaceList>

src/OpenFOAM/meshes/primitiveMesh/PatchTools/PatchToolsGatherAndMerge.C

@@ -46,17 +46,19 @@ void Foam::PatchTools::gatherAndMerge
     <
         typename PrimitivePatch<FaceList, PointField>::face_type
     >& mergedFaces,
+    globalIndex& pointAddr,
+    globalIndex& faceAddr,
     labelList& pointMergeMap,
     const bool useLocal
 )
 {
-    typedef typename PrimitivePatch<FaceList,PointField>::face_type FaceType;
+    typedef typename PrimitivePatch<FaceList, PointField>::face_type FaceType;
 
     // Faces from all ranks
-    const globalIndex faceAddr(pp.size(), globalIndex::gatherOnly{});
+    faceAddr = globalIndex(pp.size(), globalIndex::gatherOnly{});
 
     // Points from all ranks
-    const globalIndex pointAddr
+    pointAddr = globalIndex
     (
         (useLocal ? pp.localPoints().size() : pp.points().size()),
         globalIndex::gatherOnly{}
@@ -152,6 +154,40 @@ void Foam::PatchTools::gatherAndMerge
 }
 
 
+template<class FaceList, class PointField>
+void Foam::PatchTools::gatherAndMerge
+(
+    const scalar mergeDist,
+    const PrimitivePatch<FaceList, PointField>& pp,
+    Field
+    <
+        typename PrimitivePatch<FaceList, PointField>::point_type
+    >& mergedPoints,
+    List
+    <
+        typename PrimitivePatch<FaceList, PointField>::face_type
+    >& mergedFaces,
+    labelList& pointMergeMap,
+    const bool useLocal
+)
+{
+    globalIndex pointAddr;
+    globalIndex faceAddr;
+
+    PatchTools::gatherAndMerge<FaceList, PointField>
+    (
+        mergeDist,
+        pp,
+        mergedPoints,
+        mergedFaces,
+        pointAddr,
+        faceAddr,
+        pointMergeMap,
+        useLocal
+    );
+}
+
+
 template<class FaceList>
 void Foam::PatchTools::gatherAndMerge
 (

src/OpenFOAM/parallel/globalIndex/globalIndex.H

@@ -108,8 +108,8 @@ public:
 
     // Constructors
 
-        //- Default construct
-        globalIndex() = default;
+        //- Default construct (empty)
+        globalIndex() noexcept = default;
 
         //- Copy construct from a list of offsets.
         //- No communication required
@@ -184,9 +184,15 @@ public:
         //- Global max of localSizes
         inline label maxSize() const;
 
+
+    // Access
+
         //- Const-access to the offsets
         inline const labelList& offsets() const noexcept;
 
+        //- Write-access to the offsets, for changing after construction
+        inline labelList& offsets() noexcept;
+
 
     // Dimensions
 
@@ -202,8 +208,8 @@ public:
 
     // Edit
 
-        //- Write-access to the offsets, for changing after construction
-        inline labelList& offsets() noexcept;
+        //- Reset to be empty (no offsets)
+        inline void clear();
 
         //- Reset from local size, using gather/broadcast
         //- with default/specified communicator if parallel.

src/OpenFOAM/parallel/globalIndex/globalIndexI.H

@@ -176,6 +176,12 @@ inline Foam::labelList& Foam::globalIndex::offsets() noexcept
 }
 
 
+inline void Foam::globalIndex::clear()
+{
+    offsets_.clear();
+}
+
+
 inline const Foam::labelUList Foam::globalIndex::localStarts() const
 {
     const label len = (offsets_.size() - 1);

src/OpenFOAM/primitives/Scalar/scalar/scalar.C

@@ -40,6 +40,24 @@ Foam::scalar Foam::readScalar(Istream& is)
 }
 
 
+Foam::scalar Foam::readScalarOrDefault(Istream& is, const scalar defaultValue)
+{
+    if (is.good())
+    {
+        token tok(is);
+
+        if (tok.isNumber())
+        {
+            return tok.scalarToken();
+        }
+
+        is.putBack(tok);
+    }
+
+    return defaultValue;
+}
+
+
 Foam::scalar Foam::readRawScalar(Istream& is)
 {
     scalar val(0);

src/OpenFOAM/primitives/Scalar/scalar/scalar.H

@@ -103,6 +103,9 @@ namespace Foam
     //- Read scalar from stream.
     scalar readScalar(Istream& is);
 
+    //- Read scalar from stream if present or return default value
+    scalar readScalarOrDefault(Istream& is, const scalar defaultValue);
+
     //- Read raw scalar from binary stream.
     //  \note No internal check for binary vs ascii,
     //        the caller knows what they are doing
@@ -172,6 +175,9 @@ namespace Foam
     //- Read scalar from stream.
     scalar readScalar(Istream& is);
 
+    //- Read scalar from stream if present or return default value
+    scalar readScalarOrDefault(Istream& is, const scalar defaultValue);
+
     //- Read raw scalar from binary stream.
     //  \note No internal check for binary vs ascii,
     //        the caller knows what they are doing

src/Pstream/dummy/UOPstreamWrite.C

@@ -41,7 +41,7 @@ bool Foam::UOPstream::bufferIPCsend()
 
 bool Foam::UOPstream::write
 (
-    const commsTypes commsType,
+    const UPstream::commsTypes commsType,
     const int toProcNo,
     const char* buf,
     const std::streamsize bufSize,

src/Pstream/mpi/UOPstreamWrite.C

@@ -52,7 +52,7 @@ bool Foam::UOPstream::bufferIPCsend()
 
 bool Foam::UOPstream::write
 (
-    const commsTypes commsType,
+    const UPstream::commsTypes commsType,
     const int toProcNo,
     const char* buf,
     const std::streamsize bufSize,
@@ -65,7 +65,7 @@ bool Foam::UOPstream::write
         Pout<< "UOPstream::write : starting write to:" << toProcNo
             << " tag:" << tag
             << " comm:" << communicator << " size:" << label(bufSize)
-            << " commsType:" << UPstream::commsTypeNames[commsType]
+            << " commType:" << UPstream::commsTypeNames[commsType]
             << Foam::endl;
     }
     if (UPstream::warnComm != -1 && communicator != UPstream::warnComm)
@@ -73,7 +73,7 @@ bool Foam::UOPstream::write
         Pout<< "UOPstream::write : starting write to:" << toProcNo
             << " tag:" << tag
             << " comm:" << communicator << " size:" << label(bufSize)
-            << " commsType:" << UPstream::commsTypeNames[commsType]
+            << " commType:" << UPstream::commsTypeNames[commsType]
             << " warnComm:" << UPstream::warnComm
             << Foam::endl;
         error::printStack(Pout);
@@ -154,7 +154,7 @@ bool Foam::UOPstream::write
         {
             Pout<< "UOPstream::write : started write to:" << toProcNo
                 << " tag:" << tag << " size:" << label(bufSize)
-                << " commsType:" << UPstream::commsTypeNames[commsType]
+                << " commType:" << UPstream::commsTypeNames[commsType]
                 << " request:" << PstreamGlobals::outstandingRequests_.size()
                 << Foam::endl;
         }

src/TurbulenceModels/incompressible/turbulentTransportModels/turbulentTransportModels.C

@@ -126,6 +126,9 @@ makeLESModel(dynamicKEqn);
 #include "dynamicLagrangian.H"
 makeLESModel(dynamicLagrangian);
 
+#include "sigma.H"
+makeLESModel(sigma);
+
 #include "SpalartAllmarasDES.H"
 makeLESModel(SpalartAllmarasDES);
 

src/TurbulenceModels/schemes/DEShybrid/DEShybrid.H

@@ -7,6 +7,7 @@
 -------------------------------------------------------------------------------
     Copyright (C) 2015 OpenFOAM Foundation
     Copyright (C) 2015-2022 OpenCFD Ltd.
+    Copyright (C) 2022 Upstream CFD GmbH
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -28,7 +29,8 @@ Class
     Foam::DEShybrid
 
 Description
-    Hybrid convection scheme of Travin et al. for hybrid RAS/LES calculations.
+    Improved hybrid convection scheme of Travin et al. for hybrid RAS/LES
+    calculations with enhanced Grey Area Mitigation (GAM) behaviour.
 
     The scheme provides a blend between two convection schemes, based on local
     properties including the wall distance, velocity gradient and eddy
@@ -53,25 +55,26 @@ Description
 
     First published in:
     \verbatim
-        A. Travin, M. Shur, M. Strelets, P. Spalart (2000).
-        Physical and numerical upgrades in the detached-eddy simulation of
-        complex turbulent flows.
-        In Proceedings of the 412th Euromech Colloquium on LES and Complex
-        Transition and Turbulent Flows, Munich, Germany
+        Travin, A., Shur, M., Strelets, M., & Spalart, P. R. (2000).
+        Physical and numerical upgrades in the detached-eddy
+        simulation of complex turbulent flows.
+        In LES of Complex Transitional and Turbulent Flows.
+        Proceedings of the Euromech Colloquium 412. Munich, Germany
     \endverbatim
 
-    Original publication contained a typo for C_H3 constant. Corrected version
-    with minor changes for 2 lower limiters published in:
+    Original publication contained a typo for \c C_H3 constant.
+    Corrected version with minor changes for 2 lower limiters published in:
     \verbatim
-        P. Spalart, M. Shur, M. Strelets, A. Travin (2012).
-        Sensitivity of Landing-Gear Noise Predictions by Large-Eddy
-        Simulation to Numerics and Resolution.
-        AIAA Paper 2012-1174, 50th AIAA Aerospace Sciences Meeting,
-        Nashville / TN, Jan. 2012
+        Spalart, P., Shur, M., Strelets, M., & Travin, A. (2012).
+        Sensitivity of landing-gear noise predictions by large-eddy
+        simulation to numerics and resolution.
+        In 50th AIAA Aerospace Sciences Meeting Including the
+        New Horizons Forum and Aerospace Exposition. Nashville, US.
+        DOI:10.2514/6.2012-1174
     \endverbatim
 
-    Example of the DEShybrid scheme specification using linear within the LES
-    region and linearUpwind within the RAS region:
+    Example of the \c DEShybrid scheme specification using \c linear
+    within the LES region and \c linearUpwind within the RAS region:
     \verbatim
     divSchemes
     {
@@ -80,13 +83,14 @@ Description
         div(phi,U)      Gauss DEShybrid
             linear                    // scheme 1
             linearUpwind grad(U)      // scheme 2
-            hmax                      // LES delta name, e.g. 'delta', 'hmax'
-            0.65                      // DES coefficient, typically = 0.65
+            delta                     // LES delta name, e.g. 'delta', 'hmax'
+            0.65                      // CDES coefficient
             30                        // Reference velocity scale
             2                         // Reference length scale
             0                         // Minimum sigma limit (0-1)
             1                         // Maximum sigma limit (0-1)
-            1.0e-03;                  // Limiter of B function, typically 1e-03
+            1.0e-03                   // Limiter of B function, typically 1e-03
+            1.0;                      // nut limiter (if > 1, GAM extension is active)
         .
         .
     }
@@ -97,12 +101,12 @@ Notes
     be used in the detached/vortex shedding regions.
   - Scheme 2 should be an upwind/deferred correction/TVD scheme which will
     be used in the free-stream/Euler/boundary layer regions.
-  - the scheme is compiled into a separate library, and not available to
+  - The scheme is compiled into a separate library, and not available to
     solvers by default.  In order to use the scheme, add the library as a
     run-time loaded library in the \$FOAM\_CASE/system/controlDict
     dictionary, e.g.:
     \verbatim
-        libs ("libturbulenceModelSchemes.so");
+        libs (turbulenceModelSchemes);
     \endverbatim
 
 SourceFiles
@@ -110,15 +114,14 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef DEShybrid_H
-#define DEShybrid_H
+#ifndef Foam_DEShybrid_H
+#define Foam_DEShybrid_H
 
 #include "surfaceInterpolationScheme.H"
 #include "surfaceInterpolate.H"
 #include "fvcGrad.H"
 #include "blendedSchemeBase.H"
-#include "turbulentTransportModel.H"
-#include "turbulentFluidThermoModel.H"
+#include "turbulenceModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -135,6 +138,9 @@ class DEShybrid
     public surfaceInterpolationScheme<Type>,
     public blendedSchemeBase<Type>
 {
+    typedef GeometricField<Type, fvPatchField, volMesh> VolFieldType;
+    typedef GeometricField<Type, fvsPatchField, surfaceMesh> SurfaceFieldType;
+
     // Private Data
 
         //- Scheme 1
@@ -146,7 +152,7 @@ class DEShybrid
         //- Name of the LES delta field
         word deltaName_;
 
-        //- DES Coefficient
+        //- DES coefficient
         scalar CDES_;
 
         //- Reference velocity scale [m/s]
@@ -164,10 +170,14 @@ class DEShybrid
         //- Limiter of B function
         scalar OmegaLim_;
 
+        //- Limiter for modified GAM behaviour
+        scalar nutLim_;
+
         //- Scheme constants
         scalar CH1_;
         scalar CH2_;
         scalar CH3_;
+        scalar Cs_;
 
         //- No copy construct
         DEShybrid(const DEShybrid&) = delete;
@@ -178,51 +188,108 @@ class DEShybrid
 
     // Private Member Functions
 
+        //- Check the scheme coefficients
+        void checkValues()
+        {
+            if (U0_.value() <= 0)
+            {
+                FatalErrorInFunction
+                    << "U0 coefficient must be > 0. "
+                    << "Current value: " << U0_ << exit(FatalError);
+            }
+            if (L0_.value() <= 0)
+            {
+                FatalErrorInFunction
+                    << "L0 coefficient must be > 0. "
+                    << "Current value: " << L0_ << exit(FatalError);
+            }
+            if (sigmaMin_ < 0)
+            {
+                FatalErrorInFunction
+                    << "sigmaMin coefficient must be >= 0. "
+                    << "Current value: " << sigmaMin_ << exit(FatalError);
+            }
+            if (sigmaMax_ < 0)
+            {
+                FatalErrorInFunction
+                    << "sigmaMax coefficient must be >= 0. "
+                    << "Current value: " << sigmaMax_ << exit(FatalError);
+            }
+            if (sigmaMin_ > 1)
+            {
+                FatalErrorInFunction
+                    << "sigmaMin coefficient must be <= 1. "
+                    << "Current value: " << sigmaMin_ << exit(FatalError);
+            }
+            if (sigmaMax_ > 1)
+            {
+                FatalErrorInFunction
+                    << "sigmaMax coefficient must be <= 1. "
+                    << "Current value: " << sigmaMax_ << exit(FatalError);
+            }
+
+            if (debug)
+            {
+                Info<< type() << "coefficients:" << nl
+                    << "    delta : " << deltaName_ << nl
+                    << "    CDES : " << CDES_ << nl
+                    << "    U0 : " << U0_.value() << nl
+                    << "    L0 : " << L0_.value() << nl
+                    << "    sigmaMin : " << sigmaMin_ << nl
+                    << "    sigmaMax : " << sigmaMax_ << nl
+                    << "    OmegaLim : " << OmegaLim_ << nl
+                    << "    nutLim : " << nutLim_ << nl
+                    << "    CH1 : " << CH1_ << nl
+                    << "    CH2 : " << CH2_ << nl
+                    << "    CH3 : " << CH3_ << nl
+                    << "    Cs : " << Cs_ << nl
+                    << endl;
+            }
+        }
+
+
         //- Calculate the blending factor
         tmp<surfaceScalarField> calcBlendingFactor
         (
-            const GeometricField<Type, fvPatchField, volMesh>& vf,
-            const volScalarField& nuEff,
+            const VolFieldType& vf,
+            const volScalarField& nut,
+            const volScalarField& nu,
             const volVectorField& U,
             const volScalarField& delta
         ) const
         {
-            tmp<volTensorField> gradU(fvc::grad(U));
-            const volScalarField S(sqrt(2.0)*mag(symm(gradU())));
-            const volScalarField Omega(sqrt(2.0)*mag(skew(gradU())));
+            tmp<volTensorField> tgradU = fvc::grad(U);
+            const volTensorField& gradU = tgradU.cref();
+            const volScalarField S(sqrt(2.0)*mag(symm(gradU)));
+            const volScalarField Omega(sqrt(2.0)*mag(skew(tgradU)));
             const dimensionedScalar tau0_ = L0_/U0_;
 
-            const volScalarField B
-            (
+            tmp<volScalarField> tB =
                 CH3_*Omega*max(S, Omega)
-               /max(0.5*(sqr(S) + sqr(Omega)), sqr(OmegaLim_/tau0_))
-            );
+               /max(0.5*(sqr(S) + sqr(Omega)), sqr(OmegaLim_/tau0_));
 
-            const volScalarField K
-            (
-                max(Foam::sqrt(0.5*(sqr(S) + sqr(Omega))), 0.1/tau0_)
-            );
+            tmp<volScalarField> tg = tanh(pow4(tB));
 
-            const volScalarField lTurb
-            (
+            tmp<volScalarField> tK =
+                max(Foam::sqrt(0.5*(sqr(S) + sqr(Omega))), 0.1/tau0_);
+
+            tmp<volScalarField> tlTurb =
                 Foam::sqrt
                 (
                     max
                     (
-                        nuEff/(pow(0.09, 1.5)*K),
-                        dimensionedScalar("l0", sqr(dimLength), 0)
+                        (max(nut, min(sqr(Cs_*delta)*S, nutLim_*nut)) + nu)
+                       /(pow(0.09, 1.5)*tK),
+                        dimensionedScalar(sqr(dimLength), Zero)
                     )
-                )
-            );
-
-            const volScalarField g(tanh(pow4(B)));
+                );
 
             const volScalarField A
             (
                 CH2_*max
                 (
                     scalar(0),
-                    CDES_*delta/max(lTurb*g, SMALL*L0_) - 0.5
+                    CDES_*delta/max(tlTurb*tg, SMALL*L0_) - 0.5
                 )
             );
 
@@ -297,14 +364,8 @@ public:
         DEShybrid(const fvMesh& mesh, Istream& is)
         :
             surfaceInterpolationScheme<Type>(mesh),
-            tScheme1_
-            (
-                surfaceInterpolationScheme<Type>::New(mesh, is)
-            ),
-            tScheme2_
-            (
-                surfaceInterpolationScheme<Type>::New(mesh, is)
-            ),
+            tScheme1_(surfaceInterpolationScheme<Type>::New(mesh, is)),
+            tScheme2_(surfaceInterpolationScheme<Type>::New(mesh, is)),
             deltaName_(is),
             CDES_(readScalar(is)),
             U0_("U0", dimLength/dimTime, readScalar(is)),
@@ -312,46 +373,13 @@ public:
             sigmaMin_(readScalar(is)),
             sigmaMax_(readScalar(is)),
             OmegaLim_(readScalar(is)),
+            nutLim_(readScalarOrDefault(is, scalar(1))),
             CH1_(3.0),
             CH2_(1.0),
-            CH3_(2.0)
+            CH3_(2.0),
+            Cs_(0.18)
         {
-            if (U0_.value() <= 0)
-            {
-                FatalErrorInFunction
-                    << "U0 coefficient must be > 0. "
-                    << "Current value: " << U0_ << exit(FatalError);
-            }
-            if (L0_.value() <= 0)
-            {
-                FatalErrorInFunction
-                    << "L0 coefficient must be > 0. "
-                    << "Current value: " << L0_ << exit(FatalError);
-            }
-            if (sigmaMin_ < 0)
-            {
-                FatalErrorInFunction
-                    << "sigmaMin coefficient must be >= 0. "
-                    << "Current value: " << sigmaMin_ << exit(FatalError);
-            }
-            if (sigmaMax_ < 0)
-            {
-                FatalErrorInFunction
-                    << "sigmaMax coefficient must be >= 0. "
-                    << "Current value: " << sigmaMax_ << exit(FatalError);
-            }
-            if (sigmaMin_ > 1)
-            {
-                FatalErrorInFunction
-                    << "sigmaMin coefficient must be <= 1. "
-                    << "Current value: " << sigmaMin_ << exit(FatalError);
-            }
-            if (sigmaMax_ > 1)
-            {
-                FatalErrorInFunction
-                    << "sigmaMax coefficient must be <= 1. "
-                    << "Current value: " << sigmaMax_ << exit(FatalError);
-            }
+            checkValues();
         }
 
         //- Construct from mesh, faceFlux and Istream
@@ -378,46 +406,13 @@ public:
             sigmaMin_(readScalar(is)),
             sigmaMax_(readScalar(is)),
             OmegaLim_(readScalar(is)),
+            nutLim_(readScalarOrDefault(is, scalar(1))),
             CH1_(3.0),
             CH2_(1.0),
-            CH3_(2.0)
+            CH3_(2.0),
+            Cs_(0.18)
         {
-            if (U0_.value() <= 0)
-            {
-                FatalErrorInFunction
-                    << "U0 coefficient must be > 0. "
-                    << "Current value: " << U0_ << exit(FatalError);
-            }
-            if (L0_.value() <= 0)
-            {
-                FatalErrorInFunction
-                    << "L0 coefficient must be > 0. "
-                    << "Current value: " << U0_ << exit(FatalError);
-            }
-            if (sigmaMin_ < 0)
-            {
-                FatalErrorInFunction
-                    << "sigmaMin coefficient must be >= 0. "
-                    << "Current value: " << sigmaMin_ << exit(FatalError);
-            }
-            if (sigmaMax_ < 0)
-            {
-                FatalErrorInFunction
-                    << "sigmaMax coefficient must be >= 0. "
-                    << "Current value: " << sigmaMax_ << exit(FatalError);
-            }
-            if (sigmaMin_ > 1)
-            {
-                FatalErrorInFunction
-                    << "sigmaMin coefficient must be <= 1. "
-                    << "Current value: " << sigmaMin_ << exit(FatalError);
-            }
-            if (sigmaMax_ > 1)
-            {
-                FatalErrorInFunction
-                    << "sigmaMax coefficient must be <= 1. "
-                    << "Current value: " << sigmaMax_ << exit(FatalError);
-            }
+            checkValues();
         }
 
 
@@ -431,33 +426,24 @@ public:
         {
             const fvMesh& mesh = this->mesh();
 
-            typedef compressible::turbulenceModel cmpModel;
-            typedef incompressible::turbulenceModel icoModel;
+            // Retrieve LES delta from the mesh database
+            const auto& delta =
+                mesh.lookupObject<const volScalarField>(deltaName_);
 
-            // Lookup the LES delta from the mesh database
-            const volScalarField& delta = this->mesh().template
-                 lookupObject<const volScalarField>(deltaName_);
+            // Retrieve turbulence model from the mesh database
+            const auto* modelPtr =
+                mesh.cfindObject<turbulenceModel>
+                (
+                    turbulenceModel::propertiesName
+                );
 
-            // Could avoid the compressible/incompressible case by looking
-            // up all fields from the database - but retrieving from model
-            // ensures consistent fields are being employed e.g. for multiphase
-            // where group name is used
-
-            if (mesh.foundObject<icoModel>(icoModel::propertiesName))
+            if (modelPtr)
             {
-                const icoModel& model =
-                    mesh.lookupObject<icoModel>(icoModel::propertiesName);
-
-                return calcBlendingFactor(vf, model.nuEff(), model.U(), delta);
-            }
-            else if (mesh.foundObject<cmpModel>(cmpModel::propertiesName))
-            {
-                const cmpModel& model =
-                    mesh.lookupObject<cmpModel>(cmpModel::propertiesName);
+                const auto& model = *modelPtr;
 
                 return calcBlendingFactor
                 (
-                    vf, model.muEff()/model.rho(), model.U(), delta
+                    vf, model.nut(), model.nu(), model.U(), delta
                 );
             }
 
@@ -471,12 +457,9 @@ public:
 
 
         //- Return the interpolation weighting factors
-        tmp<surfaceScalarField> weights
-        (
-            const GeometricField<Type, fvPatchField, volMesh>& vf
-        ) const
+        tmp<surfaceScalarField> weights(const VolFieldType& vf) const
         {
-            surfaceScalarField bf(blendingFactor(vf));
+            const surfaceScalarField bf(blendingFactor(vf));
 
             return
                 (scalar(1) - bf)*tScheme1_().weights(vf)
@@ -485,14 +468,10 @@ public:
 
 
         //- Return the face-interpolate of the given cell field
-        //  with explicit correction
-        tmp<GeometricField<Type, fvsPatchField, surfaceMesh>>
-        interpolate
-        (
-            const GeometricField<Type, fvPatchField, volMesh>& vf
-        ) const
+        //- with explicit correction
+        tmp<SurfaceFieldType> interpolate(const VolFieldType& vf) const
         {
-            surfaceScalarField bf(blendingFactor(vf));
+            const surfaceScalarField bf(blendingFactor(vf));
 
             return
                 (scalar(1) - bf)*tScheme1_().interpolate(vf)
@@ -508,14 +487,10 @@ public:
 
 
         //- Return the explicit correction to the face-interpolate
-        //  for the given field
-        virtual tmp<GeometricField<Type, fvsPatchField, surfaceMesh>>
-        correction
-        (
-            const GeometricField<Type, fvPatchField, volMesh>& vf
-        ) const
+        //- for the given field
+        virtual tmp<SurfaceFieldType> correction(const VolFieldType& vf) const
         {
-            surfaceScalarField bf(blendingFactor(vf));
+            const surfaceScalarField bf(blendingFactor(vf));
 
             if (tScheme1_().corrected())
             {

src/TurbulenceModels/schemes/Make/options

@@ -1,19 +1,9 @@
 EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/transportModels \
-    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
-    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
-    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
-    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
-    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude
 
 LIB_LIBS = \
     -lfiniteVolume \
     -lmeshTools \
-    -lcompressibleTransportModels \
-    -lturbulenceModels \
-    -lincompressibleTurbulenceModels \
-    -lincompressibleTransportModels \
-    -lcompressibleTurbulenceModels \
-    -lfluidThermophysicalModels
+    -lturbulenceModels

src/TurbulenceModels/turbulenceModels/Base/SpalartAllmaras/SpalartAllmarasBase.C

@@ -0,0 +1,499 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2016-2022 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 "SpalartAllmarasBase.H"
+#include "wallDist.H"
+#include "bound.H"
+#include "fvOptions.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+template<class BasicEddyViscosityModel>
+tmp<volScalarField> SpalartAllmarasBase<BasicEddyViscosityModel>::chi() const
+{
+    return nuTilda_/this->nu();
+}
+
+
+template<class BasicEddyViscosityModel>
+tmp<volScalarField> SpalartAllmarasBase<BasicEddyViscosityModel>::fv1
+(
+    const volScalarField& chi
+) const
+{
+    const volScalarField chi3("chi3", pow3(chi));
+    return chi3/(chi3 + pow3(Cv1_));
+}
+
+
+template<class BasicEddyViscosityModel>
+tmp<volScalarField> SpalartAllmarasBase<BasicEddyViscosityModel>::fv2
+(
+    const volScalarField& chi,
+    const volScalarField& fv1
+) const
+{
+    return scalar(1) - chi/(scalar(1) + chi*fv1);
+}
+
+
+template<class BasicEddyViscosityModel>
+tmp<volScalarField> SpalartAllmarasBase<BasicEddyViscosityModel>::ft2
+(
+    const volScalarField& chi
+) const
+{
+    if (ft2_)
+    {
+        return Ct3_*exp(-Ct4_*sqr(chi));
+    }
+
+    return tmp<volScalarField>::New
+    (
+        IOobject
+        (
+            "ft2",
+            this->runTime_.timeName(),
+            this->mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        this->mesh_,
+        dimensionedScalar(dimless, Zero)
+    );
+}
+
+
+template<class BasicEddyViscosityModel>
+tmp<volScalarField> SpalartAllmarasBase<BasicEddyViscosityModel>::Omega
+(
+    const volTensorField& gradU
+) const
+{
+    return sqrt(2.0)*mag(skew(gradU));
+}
+
+
+template<class BasicEddyViscosityModel>
+tmp<volScalarField> SpalartAllmarasBase<BasicEddyViscosityModel>::r
+(
+    const volScalarField& nur,
+    const volScalarField& Stilda,
+    const volScalarField& dTilda
+) const
+{
+    const dimensionedScalar eps(Stilda.dimensions(), SMALL);
+
+    tmp<volScalarField> tr =
+        min(nur/(max(Stilda, eps)*sqr(kappa_*dTilda)), scalar(10));
+
+    tr.ref().boundaryFieldRef() == 0;
+
+    return tr;
+}
+
+
+template<class BasicEddyViscosityModel>
+tmp<volScalarField::Internal> SpalartAllmarasBase<BasicEddyViscosityModel>::fw
+(
+    const volScalarField& Stilda,
+    const volScalarField& dTilda
+) const
+{
+    const volScalarField::Internal r(this->r(nuTilda_, Stilda, dTilda)()());
+    const volScalarField::Internal g(r + Cw2_*(pow6(r) - r));
+
+    return g*pow((1 + pow6(Cw3_))/(pow6(g) + pow6(Cw3_)), 1.0/6.0);
+}
+
+
+template<class BasicEddyViscosityModel>
+tmp<volScalarField> SpalartAllmarasBase<BasicEddyViscosityModel>::Stilda
+(
+    const volScalarField& chi,
+    const volScalarField& fv1,
+    const volTensorField& gradU,
+    const volScalarField& dTilda
+) const
+{
+    const volScalarField Omega(this->Omega(gradU));
+
+    return
+        max
+        (
+            Omega + fv2(chi, fv1)*nuTilda_/sqr(kappa_*dTilda),
+            Cs_*Omega
+        );
+}
+
+
+template<class BasicEddyViscosityModel>
+void SpalartAllmarasBase<BasicEddyViscosityModel>::correctNut
+(
+    const volScalarField& fv1
+)
+{
+    this->nut_ = nuTilda_*fv1;
+    this->nut_.correctBoundaryConditions();
+    fv::options::New(this->mesh_).correct(this->nut_);
+}
+
+
+template<class BasicEddyViscosityModel>
+void SpalartAllmarasBase<BasicEddyViscosityModel>::correctNut()
+{
+    correctNut(fv1(this->chi()));
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class BasicEddyViscosityModel>
+SpalartAllmarasBase<BasicEddyViscosityModel>::SpalartAllmarasBase
+(
+    const word& type,
+    const alphaField& alpha,
+    const rhoField& rho,
+    const volVectorField& U,
+    const surfaceScalarField& alphaRhoPhi,
+    const surfaceScalarField& phi,
+    const transportModel& transport,
+    const word& propertiesName
+)
+:
+    BasicEddyViscosityModel
+    (
+        type,
+        alpha,
+        rho,
+        U,
+        alphaRhoPhi,
+        phi,
+        transport,
+        propertiesName
+    ),
+
+    sigmaNut_
+    (
+        dimensioned<scalar>::getOrAddToDict
+        (
+            "sigmaNut",
+            this->coeffDict_,
+            0.66666
+        )
+    ),
+    kappa_
+    (
+        dimensioned<scalar>::getOrAddToDict
+        (
+            "kappa",
+            this->coeffDict_,
+            0.41
+        )
+    ),
+    Cb1_
+    (
+        dimensioned<scalar>::getOrAddToDict
+        (
+            "Cb1",
+            this->coeffDict_,
+            0.1355
+        )
+    ),
+    Cb2_
+    (
+        dimensioned<scalar>::getOrAddToDict
+        (
+            "Cb2",
+            this->coeffDict_,
+            0.622
+        )
+    ),
+    Cw1_(Cb1_/sqr(kappa_) + (1.0 + Cb2_)/sigmaNut_),
+    Cw2_
+    (
+        dimensioned<scalar>::getOrAddToDict
+        (
+            "Cw2",
+            this->coeffDict_,
+            0.3
+        )
+    ),
+    Cw3_
+    (
+        dimensioned<scalar>::getOrAddToDict
+        (
+            "Cw3",
+            this->coeffDict_,
+            2.0
+        )
+    ),
+    Cv1_
+    (
+        dimensioned<scalar>::getOrAddToDict
+        (
+            "Cv1",
+            this->coeffDict_,
+            7.1
+        )
+    ),
+    Cs_
+    (
+        dimensioned<scalar>::getOrAddToDict
+        (
+            "Cs",
+            this->coeffDict_,
+            0.3
+        )
+    ),
+    ck_
+    (
+        dimensioned<scalar>::getOrAddToDict
+        (
+            "ck",
+            this->coeffDict_,
+            0.07
+        )
+    ),
+    ft2_
+    (
+        Switch::getOrAddToDict
+        (
+            "ft2",
+            this->coeffDict_,
+            false
+        )
+    ),
+    Ct3_
+    (
+        dimensioned<scalar>::getOrAddToDict
+        (
+            "Ct3",
+            this->coeffDict_,
+            1.2
+        )
+    ),
+    Ct4_
+    (
+        dimensioned<scalar>::getOrAddToDict
+        (
+            "Ct4",
+            this->coeffDict_,
+            0.5
+        )
+    ),
+
+    nuTilda_
+    (
+        IOobject
+        (
+            "nuTilda",
+            this->runTime_.timeName(),
+            this->mesh_,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        this->mesh_
+    ),
+
+    y_(wallDist::New(this->mesh_).y())
+{
+    if (ft2_)
+    {
+        Info<< "ft2 term: active" << nl;
+    }
+    else
+    {
+        Info<< "ft2 term: inactive" << nl;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class BasicEddyViscosityModel>
+bool SpalartAllmarasBase<BasicEddyViscosityModel>::read()
+{
+    if (BasicEddyViscosityModel::read())
+    {
+        sigmaNut_.readIfPresent(this->coeffDict());
+        kappa_.readIfPresent(this->coeffDict());
+
+        Cb1_.readIfPresent(this->coeffDict());
+        Cb2_.readIfPresent(this->coeffDict());
+        Cw1_ = Cb1_/sqr(kappa_) + (1.0 + Cb2_)/sigmaNut_;
+        Cw2_.readIfPresent(this->coeffDict());
+        Cw3_.readIfPresent(this->coeffDict());
+        Cv1_.readIfPresent(this->coeffDict());
+        Cs_.readIfPresent(this->coeffDict());
+        ck_.readIfPresent(this->coeffDict());
+
+        ft2_.readIfPresent("ft2", this->coeffDict());
+        Ct3_.readIfPresent(this->coeffDict());
+        Ct4_.readIfPresent(this->coeffDict());
+
+        if (ft2_)
+        {
+            Info<< "    ft2 term: active" << nl;
+        }
+        else
+        {
+            Info<< "    ft2 term: inactive" << nl;
+        }
+
+        return true;
+    }
+
+    return false;
+}
+
+
+template<class BasicEddyViscosityModel>
+tmp<volScalarField>
+SpalartAllmarasBase<BasicEddyViscosityModel>::DnuTildaEff() const
+{
+    return tmp<volScalarField>::New
+    (
+        IOobject::groupName("DnuTildaEff", this->alphaRhoPhi_.group()),
+        (nuTilda_ + this->nu())/sigmaNut_
+    );
+}
+
+
+template<class BasicEddyViscosityModel>
+tmp<volScalarField> SpalartAllmarasBase<BasicEddyViscosityModel>::k() const
+{
+    // (B:Eq. 4.50)
+    const scalar Cmu = 0.09;
+    const auto fv1 = this->fv1(chi());
+
+    return tmp<volScalarField>::New
+    (
+        IOobject::groupName("k", this->alphaRhoPhi_.group()),
+        cbrt(fv1)*nuTilda_*::sqrt(scalar(2)/Cmu)*mag(symm(fvc::grad(this->U_)))
+    );
+}
+
+template<class BasicEddyViscosityModel>
+tmp<volScalarField>
+SpalartAllmarasBase<BasicEddyViscosityModel>::epsilon() const
+{
+    // (B:Eq. 4.50)
+    const scalar Cmu = 0.09;
+    const auto fv1 = this->fv1(chi());
+    const dimensionedScalar nutSMALL(nuTilda_.dimensions(), SMALL);
+
+    return tmp<volScalarField>::New
+    (
+        IOobject::groupName("epsilon", this->alphaRhoPhi_.group()),
+        sqrt(fv1)*sqr(::sqrt(Cmu)*this->k())/(nuTilda_ + this->nut_ + nutSMALL)
+    );
+}
+
+
+template<class BasicEddyViscosityModel>
+tmp<volScalarField> SpalartAllmarasBase<BasicEddyViscosityModel>::omega() const
+{
+    // (P:p. 384)
+    const scalar betaStar = 0.09;
+    const dimensionedScalar k0(sqr(dimLength/dimTime), SMALL);
+
+    return tmp<volScalarField>::New
+    (
+        IOobject::groupName("omega", this->alphaRhoPhi_.group()),
+        this->epsilon()/(betaStar*(this->k() + k0))
+    );
+}
+
+
+template<class BasicEddyViscosityModel>
+void SpalartAllmarasBase<BasicEddyViscosityModel>::correct()
+{
+    if (!this->turbulence_)
+    {
+        return;
+    }
+
+    {
+        // Local references
+        const alphaField& alpha = this->alpha_;
+        const rhoField& rho = this->rho_;
+        const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_;
+        const volVectorField& U = this->U_;
+        fv::options& fvOptions(fv::options::New(this->mesh_));
+
+        BasicEddyViscosityModel::correct();
+
+        const volScalarField chi(this->chi());
+        const volScalarField fv1(this->fv1(chi));
+        const volScalarField ft2(this->ft2(chi));
+
+        tmp<volTensorField> tgradU = fvc::grad(U);
+        volScalarField dTilda(this->dTilda(chi, fv1, tgradU()));
+        volScalarField Stilda(this->Stilda(chi, fv1, tgradU(), dTilda));
+        tgradU.clear();
+
+        tmp<fvScalarMatrix> nuTildaEqn
+        (
+            fvm::ddt(alpha, rho, nuTilda_)
+          + fvm::div(alphaRhoPhi, nuTilda_)
+          - fvm::laplacian(alpha*rho*DnuTildaEff(), nuTilda_)
+          - Cb2_/sigmaNut_*alpha()*rho()*magSqr(fvc::grad(nuTilda_)()())
+         ==
+            Cb1_*alpha()*rho()*Stilda()*nuTilda_()*(scalar(1) - ft2())
+          - fvm::Sp
+            (
+                (Cw1_*fw(Stilda, dTilda) - Cb1_/sqr(kappa_)*ft2())
+               *alpha()*rho()*nuTilda_()/sqr(dTilda()),
+                nuTilda_
+            )
+          + fvOptions(alpha, rho, nuTilda_)
+        );
+
+        nuTildaEqn.ref().relax();
+        fvOptions.constrain(nuTildaEqn.ref());
+        solve(nuTildaEqn);
+        fvOptions.correct(nuTilda_);
+        bound(nuTilda_, dimensionedScalar(nuTilda_.dimensions(), Zero));
+        nuTilda_.correctBoundaryConditions();
+    }
+
+    correctNut();
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/TurbulenceModels/turbulenceModels/Base/SpalartAllmaras/SpalartAllmarasBase.H

@@ -0,0 +1,246 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2011-2016 OpenFOAM Foundation
+    Copyright (C) 2015-2022 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::LESModels::SpalartAllmarasBase
+
+Group
+    grpDESTurbulence
+
+Description
+    Base class to handle various characteristics for \c SpalartAllmaras based
+    LES/DES turbulence models for incompressible and compressible flows.
+
+    References:
+    \verbatim
+        Standard model:
+            Spalart, P.R., & Allmaras, S.R. (1994).
+            A one-equation turbulence model for aerodynamic flows.
+            La Recherche Aerospatiale, 1, 5-21.
+
+        Standard model:
+            Spalart, P. R., Jou, W. H., Strelets, M., & Allmaras, S. R. (1997).
+            Comments on the feasibility of LES for wings, and on a hybrid
+            RANS/LES approach.
+            Advances in DNS/LES, 1, 4-8.
+
+        Estimation expression for k and epsilon (tag:B), Eq. 4.50:
+            Bourgoin, A. (2019).
+            Bathymetry induced turbulence modelling the
+            Alderney Race site: regional approach with TELEMAC-LES.
+            Normandie Université.
+
+        Estimation expressions for omega (tag:P):
+            Pope, S. B. (2000).
+            Turbulent flows.
+            Cambridge, UK: Cambridge Univ. Press
+            DOI:10.1017/CBO9780511840531
+    \endverbatim
+
+SourceFiles
+    SpalartAllmarasBase.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_SpalartAllmarasBase_H
+#define Foam_SpalartAllmarasBase_H
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                    Class SpalartAllmarasBase Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class BasicEddyViscosityModel>
+class SpalartAllmarasBase
+:
+    public BasicEddyViscosityModel
+{
+    // Private Member Functions
+
+        //- No copy construct
+        SpalartAllmarasBase(const SpalartAllmarasBase&) = delete;
+
+        //- No copy assignment
+        void operator=(const SpalartAllmarasBase&) = delete;
+
+
+protected:
+
+    // Protected Data
+
+        // Model constants
+
+            dimensionedScalar sigmaNut_;
+            dimensionedScalar kappa_;
+
+            dimensionedScalar Cb1_;
+            dimensionedScalar Cb2_;
+            dimensionedScalar Cw1_;
+            dimensionedScalar Cw2_;
+            dimensionedScalar Cw3_;
+            dimensionedScalar Cv1_;
+            dimensionedScalar Cs_;
+            dimensionedScalar ck_;
+
+            Switch ft2_;
+            dimensionedScalar Ct3_;
+            dimensionedScalar Ct4_;
+
+
+        // Fields
+
+            //- Modified kinematic viscosity [m^2/s]
+            volScalarField nuTilda_;
+
+            //- Wall distance
+            //  Note: different to wall distance in parent RASModel
+            //  which is for near-wall cells only
+            const volScalarField& y_;
+
+
+    // Protected Member Functions
+
+        tmp<volScalarField> chi() const;
+
+        tmp<volScalarField> fv1(const volScalarField& chi) const;
+
+        tmp<volScalarField> fv2
+        (
+            const volScalarField& chi,
+            const volScalarField& fv1
+        ) const;
+
+        tmp<volScalarField> ft2(const volScalarField& chi) const;
+
+        tmp<volScalarField> Omega(const volTensorField& gradU) const;
+
+        tmp<volScalarField> r
+        (
+            const volScalarField& nur,
+            const volScalarField& Stilda,
+            const volScalarField& dTilda
+        ) const;
+
+        tmp<volScalarField::Internal> fw
+        (
+            const volScalarField& Stilda,
+            const volScalarField& dTilda
+        ) const;
+
+        virtual tmp<volScalarField> Stilda
+        (
+            const volScalarField& chi,
+            const volScalarField& fv1,
+            const volTensorField& gradU,
+            const volScalarField& dTilda
+        ) const;
+
+        //- Length scale
+        virtual tmp<volScalarField> dTilda
+        (
+            const volScalarField& chi,
+            const volScalarField& fv1,
+            const volTensorField& gradU
+        ) const = 0;
+
+        void correctNut(const volScalarField& fv1);
+        virtual void correctNut();
+
+
+public:
+
+    typedef typename BasicEddyViscosityModel::alphaField alphaField;
+    typedef typename BasicEddyViscosityModel::rhoField rhoField;
+    typedef typename BasicEddyViscosityModel::transportModel transportModel;
+
+
+    // Constructors
+
+        //- Construct from components
+        SpalartAllmarasBase
+        (
+            const word& type,
+            const alphaField& alpha,
+            const rhoField& rho,
+            const volVectorField& U,
+            const surfaceScalarField& alphaRhoPhi,
+            const surfaceScalarField& phi,
+            const transportModel& transport,
+            const word& propertiesName = turbulenceModel::propertiesName
+        );
+
+
+    //- Destructor
+    virtual ~SpalartAllmarasBase() = default;
+
+
+    // Member Functions
+
+        //- Re-read model coefficients if they have changed
+        virtual bool read();
+
+        //- Return the effective diffusivity for nuTilda
+        tmp<volScalarField> DnuTildaEff() const;
+
+        //- Return the (estimated) turbulent kinetic energy
+        virtual tmp<volScalarField> k() const;
+
+        //- Return the (estimated) turbulent kinetic energy dissipation rate
+        virtual tmp<volScalarField> epsilon() const;
+
+        //- Return the (estimated) specific dissipation rate
+        virtual tmp<volScalarField> omega() const;
+
+        //- Return the modified kinematic viscosity
+        tmp<volScalarField> nuTilda() const
+        {
+            return nuTilda_;
+        }
+
+        //- Correct nuTilda and related properties
+        virtual void correct();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "SpalartAllmarasBase.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/TurbulenceModels/turbulenceModels/Base/kOmegaSST/kOmegaSSTBase.C

@@ -6,7 +6,8 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2016-2020 OpenCFD Ltd.
+    Copyright (C) 2022 Upstream CFD GmbH
+    Copyright (C) 2016-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -47,7 +48,7 @@ tmp<volScalarField> kOmegaSSTBase<BasicEddyViscosityModel>::F1
     tmp<volScalarField> CDkOmegaPlus = max
     (
         CDkOmega,
-        dimensionedScalar("1.0e-10", dimless/sqr(dimTime), 1.0e-10)
+        dimensionedScalar(dimless/sqr(dimTime), 1.0e-10)
     );
 
     tmp<volScalarField> arg1 = min
@@ -132,6 +133,17 @@ void kOmegaSSTBase<BasicEddyViscosityModel>::correctNut()
 }
 
 
+template<class BasicEddyViscosityModel>
+Foam::tmp<Foam::volScalarField> kOmegaSSTBase<BasicEddyViscosityModel>::S2
+(
+    const volScalarField& F1,
+    const volTensorField& gradU
+) const
+{
+    return 2*magSqr(symm(gradU));
+}
+
+
 template<class BasicEddyViscosityModel>
 tmp<volScalarField::Internal> kOmegaSSTBase<BasicEddyViscosityModel>::Pk
 (
@@ -143,17 +155,32 @@ tmp<volScalarField::Internal> kOmegaSSTBase<BasicEddyViscosityModel>::Pk
 
 
 template<class BasicEddyViscosityModel>
-tmp<volScalarField::Internal>
-kOmegaSSTBase<BasicEddyViscosityModel>::epsilonByk
+tmp<volScalarField::Internal> kOmegaSSTBase<BasicEddyViscosityModel>::epsilonByk
 (
-    const volScalarField& F1,
-    const volTensorField& gradU
+    const volScalarField& /* F1 not used */,
+    const volTensorField& /* gradU not used */
 ) const
 {
     return betaStar_*omega_();
 }
 
 
+template<class BasicEddyViscosityModel>
+tmp<volScalarField::Internal> kOmegaSSTBase<BasicEddyViscosityModel>::GbyNu0
+(
+    const volTensorField& gradU,
+    const volScalarField& F1,
+    const volScalarField& S2
+) const
+{
+    return tmp<volScalarField::Internal>::New
+    (
+        IOobject::scopedName(this->type(), "GbyNu"),
+        gradU() && dev(twoSymm(gradU()))
+    );
+}
+
+
 template<class BasicEddyViscosityModel>
 tmp<volScalarField::Internal> kOmegaSSTBase<BasicEddyViscosityModel>::GbyNu
 (
@@ -165,8 +192,7 @@ tmp<volScalarField::Internal> kOmegaSSTBase<BasicEddyViscosityModel>::GbyNu
     return min
     (
         GbyNu0,
-        (c1_/a1_)*betaStar_*omega_()
-       *max(a1_*omega_(), b1_*F2*sqrt(S2))
+        (c1_/a1_)*betaStar_*omega_()*max(a1_*omega_(), b1_*F2*sqrt(S2))
     );
 }
 
@@ -174,13 +200,10 @@ tmp<volScalarField::Internal> kOmegaSSTBase<BasicEddyViscosityModel>::GbyNu
 template<class BasicEddyViscosityModel>
 tmp<fvScalarMatrix> kOmegaSSTBase<BasicEddyViscosityModel>::kSource() const
 {
-    return tmp<fvScalarMatrix>
+    return tmp<fvScalarMatrix>::New
     (
-        new fvScalarMatrix
-        (
-            k_,
-            dimVolume*this->rho_.dimensions()*k_.dimensions()/dimTime
-        )
+        k_,
+        dimVolume*this->rho_.dimensions()*k_.dimensions()/dimTime
     );
 }
 
@@ -188,13 +211,10 @@ tmp<fvScalarMatrix> kOmegaSSTBase<BasicEddyViscosityModel>::kSource() const
 template<class BasicEddyViscosityModel>
 tmp<fvScalarMatrix> kOmegaSSTBase<BasicEddyViscosityModel>::omegaSource() const
 {
-    return tmp<fvScalarMatrix>
+    return tmp<fvScalarMatrix>::New
     (
-        new fvScalarMatrix
-        (
-            omega_,
-            dimVolume*this->rho_.dimensions()*omega_.dimensions()/dimTime
-        )
+        omega_,
+        dimVolume*this->rho_.dimensions()*omega_.dimensions()/dimTime
     );
 }
 
@@ -207,13 +227,10 @@ tmp<fvScalarMatrix> kOmegaSSTBase<BasicEddyViscosityModel>::Qsas
     const volScalarField::Internal& beta
 ) const
 {
-    return tmp<fvScalarMatrix>
+    return tmp<fvScalarMatrix>::New
     (
-        new fvScalarMatrix
-        (
-            omega_,
-            dimVolume*this->rho_.dimensions()*omega_.dimensions()/dimTime
-        )
+        omega_,
+        dimVolume*this->rho_.dimensions()*omega_.dimensions()/dimTime
     );
 }
 
@@ -389,6 +406,7 @@ kOmegaSSTBase<BasicEddyViscosityModel>::kOmegaSSTBase
         ),
         this->mesh_
     ),
+
     decayControl_
     (
         Switch::getOrAddToDict
@@ -498,31 +516,30 @@ void kOmegaSSTBase<BasicEddyViscosityModel>::correct()
 
     BasicEddyViscosityModel::correct();
 
-    volScalarField::Internal divU(fvc::div(fvc::absolute(this->phi(), U)));
+    const volScalarField::Internal divU
+    (
+        fvc::div(fvc::absolute(this->phi(), U))
+    );
+
+    const volScalarField CDkOmega
+    (
+        (2*alphaOmega2_)*(fvc::grad(k_) & fvc::grad(omega_))/omega_
+    );
+
+    const volScalarField F1(this->F1(CDkOmega));
+    const volScalarField F23(this->F23());
 
     tmp<volTensorField> tgradU = fvc::grad(U);
-    volScalarField S2(2*magSqr(symm(tgradU())));
-    volScalarField::Internal GbyNu0
-    (
-        this->type() + ":GbyNu",
-        (tgradU() && dev(twoSymm(tgradU())))
-    );
+    const volScalarField S2(this->S2(F1, tgradU()));
+    volScalarField::Internal GbyNu0(this->GbyNu0(tgradU(), F1, S2));
     volScalarField::Internal G(this->GName(), nut*GbyNu0);
 
     // Update omega and G at the wall
     omega_.boundaryFieldRef().updateCoeffs();
 
-    volScalarField CDkOmega
-    (
-        (2*alphaOmega2_)*(fvc::grad(k_) & fvc::grad(omega_))/omega_
-    );
-
-    volScalarField F1(this->F1(CDkOmega));
-    volScalarField F23(this->F23());
-
     {
-        volScalarField::Internal gamma(this->gamma(F1));
-        volScalarField::Internal beta(this->beta(F1));
+        const volScalarField::Internal gamma(this->gamma(F1));
+        const volScalarField::Internal beta(this->beta(F1));
 
         GbyNu0 = GbyNu(GbyNu0, F23(), S2());
 
@@ -555,28 +572,30 @@ void kOmegaSSTBase<BasicEddyViscosityModel>::correct()
         bound(omega_, this->omegaMin_);
     }
 
-    // Turbulent kinetic energy equation
-    tmp<fvScalarMatrix> kEqn
-    (
-        fvm::ddt(alpha, rho, k_)
-      + fvm::div(alphaRhoPhi, k_)
-      - fvm::laplacian(alpha*rho*DkEff(F1), k_)
-     ==
-        alpha()*rho()*Pk(G)
-      - fvm::SuSp((2.0/3.0)*alpha()*rho()*divU, k_)
-      - fvm::Sp(alpha()*rho()*epsilonByk(F1, tgradU()), k_)
-      + alpha()*rho()*betaStar_*omegaInf_*kInf_
-      + kSource()
-      + fvOptions(alpha, rho, k_)
-    );
+    {
+        // Turbulent kinetic energy equation
+        tmp<fvScalarMatrix> kEqn
+        (
+            fvm::ddt(alpha, rho, k_)
+          + fvm::div(alphaRhoPhi, k_)
+          - fvm::laplacian(alpha*rho*DkEff(F1), k_)
+         ==
+            alpha()*rho()*Pk(G)
+          - fvm::SuSp((2.0/3.0)*alpha()*rho()*divU, k_)
+          - fvm::Sp(alpha()*rho()*epsilonByk(F1, tgradU()), k_)
+          + alpha()*rho()*betaStar_*omegaInf_*kInf_
+          + kSource()
+          + fvOptions(alpha, rho, k_)
+        );
 
-    tgradU.clear();
+        tgradU.clear();
 
-    kEqn.ref().relax();
-    fvOptions.constrain(kEqn.ref());
-    solve(kEqn);
-    fvOptions.correct(k_);
-    bound(k_, this->kMin_);
+        kEqn.ref().relax();
+        fvOptions.constrain(kEqn.ref());
+        solve(kEqn);
+        fvOptions.correct(k_);
+        bound(k_, this->kMin_);
+    }
 
     correctNut(S2);
 }

src/TurbulenceModels/turbulenceModels/Base/kOmegaSST/kOmegaSSTBase.H

@@ -6,7 +6,8 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2017-2021 OpenCFD Ltd.
+    Copyright (C) 2022 Upstream CFD GmbH
+    Copyright (C) 2017-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -141,7 +142,7 @@ class kOmegaSSTBase
 
 protected:
 
-    // Protected data
+    // Protected Data
 
         // Model coefficients
 
@@ -174,7 +175,10 @@ protected:
             //  which is for near-wall cells only
             const volScalarField& y_;
 
+            //- Turbulent kinetic energy field [m^2/s^2]
             volScalarField k_;
+
+            //- Specific dissipation rate field [1/s]
             volScalarField omega_;
 
 
@@ -188,6 +192,7 @@ protected:
 
     // Protected Member Functions
 
+        //- Set decay control with kInf and omegaInf
         void setDecayControl(const dictionary& dict);
 
         virtual tmp<volScalarField> F1(const volScalarField& CDkOmega) const;
@@ -195,6 +200,7 @@ protected:
         virtual tmp<volScalarField> F3() const;
         virtual tmp<volScalarField> F23() const;
 
+        //- Return the blended field
         tmp<volScalarField> blend
         (
             const volScalarField& F1,
@@ -205,6 +211,7 @@ protected:
             return F1*(psi1 - psi2) + psi2;
         }
 
+        //- Return the internal blended field
         tmp<volScalarField::Internal> blend
         (
             const volScalarField::Internal& F1,
@@ -253,6 +260,13 @@ protected:
 
         virtual void correctNut();
 
+        //- Return square of strain rate
+        virtual tmp<volScalarField> S2
+        (
+            const volScalarField& F1,
+            const volTensorField& gradU
+        ) const;
+
         //- Return k production rate
         virtual tmp<volScalarField::Internal> Pk
         (
@@ -266,6 +280,14 @@ protected:
             const volTensorField& gradU
         ) const;
 
+        //- Return (G/nu)_0
+        virtual tmp<volScalarField::Internal> GbyNu0
+        (
+            const volTensorField& gradU,
+            const volScalarField& F1,
+            const volScalarField& S2
+        ) const;
+
         //- Return G/nu
         virtual tmp<volScalarField::Internal> GbyNu
         (
@@ -321,22 +343,20 @@ public:
         //- Return the effective diffusivity for k
         tmp<volScalarField> DkEff(const volScalarField& F1) const
         {
-            return tmp<volScalarField>
+            return tmp<volScalarField>::New
             (
-                new volScalarField("DkEff", alphaK(F1)*this->nut_ + this->nu())
+                "DkEff",
+                alphaK(F1)*this->nut_ + this->nu()
             );
         }
 
         //- Return the effective diffusivity for omega
         tmp<volScalarField> DomegaEff(const volScalarField& F1) const
         {
-            return tmp<volScalarField>
+            return tmp<volScalarField>::New
             (
-                new volScalarField
-                (
-                    "DomegaEff",
-                    alphaOmega(F1)*this->nut_ + this->nu()
-                )
+                "DomegaEff",
+                alphaOmega(F1)*this->nut_ + this->nu()
             );
         }
 

src/TurbulenceModels/turbulenceModels/DES/DESModel/DESModel.C

@@ -6,6 +6,8 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2012-2015 OpenFOAM Foundation
+    Copyright (C) 2022 Upstream CFD GmbH
+    Copyright (C) 2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -59,15 +61,127 @@ DESModel<BasicTurbulenceModel>::DESModel
         phi,
         transport,
         propertiesName
-    )
-{}
+    ),
+    Ctrans_(dimless, Zero)
+{
+    // Note: Ctrans is model-specific and initialised in derived classes
+}
 
 
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 template<class BasicTurbulenceModel>
-DESModel<BasicTurbulenceModel>::~DESModel()
-{}
+bool DESModel<BasicTurbulenceModel>::read()
+{
+    if (LESeddyViscosity<BasicTurbulenceModel>::read())
+    {
+        Ctrans_.readIfPresent(this->coeffDict());
+
+        return true;
+    }
+
+    return false;
+}
+
+
+template<class BasicTurbulenceModel>
+tmp<volScalarField> DESModel<BasicTurbulenceModel>::Ssigma
+(
+    const volTensorField& gradU,
+    const dimensionedScalar& coeff
+)
+{
+    // Limiter
+    const dimensionedScalar eps0(dimless, SMALL);
+    const dimensionedScalar eps2(dimless/sqr(dimTime), SMALL);
+    const dimensionedScalar eps4(dimless/pow4(dimTime), SMALL);
+    const dimensionedScalar max2(dimless/sqr(dimTime), GREAT);
+    const dimensionedTensor maxTen2
+    (
+        dimless/sqr(dimTime),
+        tensor::max
+    );
+    const dimensionedTensor minTen2
+    (
+        dimless/sqr(dimTime),
+        tensor::min
+    );
+
+    const volTensorField G(max(min(gradU.T() & gradU, maxTen2), minTen2));
+
+    // Tensor invariants
+    const volScalarField I1(tr(G));
+    const volScalarField I2(0.5*(sqr(I1) - tr(G & G)));
+    tmp<volScalarField> tI3 = det(G);
+
+    const volScalarField alpha1(max(sqr(I1)/9.0 - I2/3.0, eps4));
+
+    tmp<volScalarField> talpha2 =
+        pow3(min(I1, max2))/27.0 - I1*I2/6.0 + 0.5*tI3;
+
+    const volScalarField alpha3
+    (
+        1.0/3.0
+       *acos
+        (
+            max
+            (
+                scalar(-1) + eps0,
+                min(scalar(1) - eps0, talpha2/pow(alpha1, 3.0/2.0))
+            )
+        )
+    );
+
+    const scalar piBy3 = constant::mathematical::pi/3.0;
+    const volScalarField sigma1
+    (
+        sqrt(max(I1/3.0 + 2.0*sqrt(alpha1)*cos(alpha3), eps2))
+    );
+    const volScalarField sigma2
+    (
+        sqrt(max(I1/3.0 - 2.0*sqrt(alpha1)*cos(piBy3 + alpha3), eps2))
+    );
+    const volScalarField sigma3
+    (
+        sqrt(max(I1/3.0 - 2.0*sqrt(alpha1)*cos(piBy3 - alpha3), eps2))
+    );
+
+    return
+        coeff
+       *sigma3
+       *(sigma1 - sigma2)
+       *(sigma2 - sigma3)
+       /max(sqr(sigma1), eps2);
+}
+
+
+template<class BasicTurbulenceModel>
+tmp<volScalarField> DESModel<BasicTurbulenceModel>::Ssigma
+(
+    const volTensorField& gradU
+) const
+{
+    return Ssigma(gradU, Ctrans_);
+}
+
+
+template<class BasicTurbulenceModel>
+tmp<volScalarField> DESModel<BasicTurbulenceModel>::fd() const
+{
+    return tmp<volScalarField>::New
+    (
+        IOobject
+        (
+            "fd",
+            this->mesh_.time().timeName(),
+            this->mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        this->mesh_,
+        dimensionedScalar(dimless, Zero)
+    );
+}
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/TurbulenceModels/turbulenceModels/DES/DESModel/DESModel.H

@@ -6,7 +6,8 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2012-2015 OpenFOAM Foundation
-    Copyright (C) 2016 OpenCFD Ltd.
+    Copyright (C) 2022 Upstream CFD GmbH
+    Copyright (C) 2016, 2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -57,9 +58,6 @@ class DESModel
     public DESModelBase,
     public LESeddyViscosity<BasicTurbulenceModel>
 {
-
-private:
-
     // Private Member Functions
 
         //- No copy construct
@@ -69,6 +67,14 @@ private:
         void operator=(const DESModel&) = delete;
 
 
+protected:
+
+    // Protected Data
+
+        //- Model-specific transition constant
+        dimensionedScalar Ctrans_;
+
+
 public:
 
     typedef typename BasicTurbulenceModel::alphaField alphaField;
@@ -92,13 +98,30 @@ public:
 
 
     //- Destructor
-    virtual ~DESModel();
+    virtual ~DESModel() = default;
 
 
     // Public Member Functions
 
+        //- Re-read model coefficients if they have changed
+        virtual bool read();
+
         //- Return the LES field indicator
         virtual tmp<volScalarField> LESRegion() const = 0;
+
+        //- Return modified strain rate
+        static tmp<volScalarField> Ssigma
+        (
+            const volTensorField& gradU,
+            const dimensionedScalar& coeff
+        );
+
+        //- Return modified strain rate
+        //  Note: uses Ctrans_ coefficient
+        virtual tmp<volScalarField> Ssigma(const volTensorField& gradU) const;
+
+        //- Return the shielding function
+        virtual tmp<volScalarField> fd() const;
 };
 
 

src/TurbulenceModels/turbulenceModels/DES/DESModel/DESModelBase.C

@@ -23,16 +23,6 @@ License
     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::DESModelBase
-
-Description
-    Base class for DES models providing an interfaces to the LESRegion
-    function.
-
-SourceFiles
-    DESModelBase.C
-
 \*---------------------------------------------------------------------------*/
 
 #include "DESModelBase.H"

src/TurbulenceModels/turbulenceModels/DES/DESModel/DESModelBase.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2015 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -28,8 +28,7 @@ Class
     Foam::DESModelBase
 
 Description
-    Base class for DES models providing an interfaces to the LESRegion
-    function.
+    Base class for DES models providing an interfaces to DES fields.
 
 SourceFiles
     DESModelBase.C
@@ -70,6 +69,9 @@ public:
 
         //- Return the LES field indicator
         virtual tmp<volScalarField> LESRegion() const = 0;
+
+        //- Return the shielding function
+        virtual tmp<volScalarField> fd() const = 0;
 };
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/TurbulenceModels/turbulenceModels/DES/SpalartAllmarasDDES/SpalartAllmarasDDES.C

@@ -6,7 +6,8 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2015-2020 OpenCFD Ltd.
+    Copyright (C) 2022 Upstream CFD GmbH
+    Copyright (C) 2015-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -37,46 +38,60 @@ namespace LESModels
 
 // * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
 
-template<class BasicTurbulenceModel>
-tmp<volScalarField> SpalartAllmarasDDES<BasicTurbulenceModel>::rd
-(
-    const volScalarField& magGradU
-) const
-{
-    tmp<volScalarField> tr
-    (
-        min
-        (
-            this->nuEff()
-           /(
-                max
-                (
-                    magGradU,
-                    dimensionedScalar("SMALL", magGradU.dimensions(), SMALL)
-                )
-               *sqr(this->kappa_*this->y_)
-            ),
-            scalar(10)
-        )
-    );
-    tr.ref().boundaryFieldRef() == 0.0;
-
-    return tr;
-}
-
-
 template<class BasicTurbulenceModel>
 tmp<volScalarField> SpalartAllmarasDDES<BasicTurbulenceModel>::fd
 (
     const volScalarField& magGradU
 ) const
 {
-    return 1 - tanh(pow(Cd1_*rd(magGradU), Cd2_));
+    return
+        1
+      - tanh(pow(this->Cd1_*this->r(this->nuEff(), magGradU, this->y_), Cd2_));
 }
 
 
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
+template<class BasicTurbulenceModel>
+tmp<volScalarField> SpalartAllmarasDDES<BasicTurbulenceModel>::Stilda
+(
+    const volScalarField& chi,
+    const volScalarField& fv1,
+    const volTensorField& gradU,
+    const volScalarField& dTilda
+) const
+{
+    if (this->useSigma_)
+    {
+        const volScalarField& lRAS(this->y_);
+        const volScalarField fv2(this->fv2(chi, fv1));
+        const volScalarField lLES(this->lengthScaleLES(chi, fv1));
+        const volScalarField Omega(this->Omega(gradU));
+        const volScalarField Ssigma(this->Ssigma(gradU));
+        const volScalarField SsigmaDES
+        (
+            Omega - fd(mag(gradU))*pos(lRAS - lLES)*(Omega - Ssigma)
+        );
+
+        return
+            max
+            (
+                SsigmaDES + fv2*this->nuTilda_/sqr(this->kappa_*dTilda),
+                this->Cs_*SsigmaDES
+            );
+    }
+
+    return
+        SpalartAllmarasBase<DESModel<BasicTurbulenceModel>>::Stilda
+        (
+            chi,
+            fv1,
+            gradU,
+            dTilda
+        );
+}
+
+
 template<class BasicTurbulenceModel>
 tmp<volScalarField> SpalartAllmarasDDES<BasicTurbulenceModel>::dTilda
 (
@@ -86,17 +101,12 @@ tmp<volScalarField> SpalartAllmarasDDES<BasicTurbulenceModel>::dTilda
 ) const
 {
     const volScalarField& lRAS(this->y_);
-    const volScalarField lLES(this->psi(chi, fv1)*this->CDES_*this->delta());
+    const volScalarField lLES(this->lengthScaleLES(chi, fv1));
+    const dimensionedScalar l0(dimLength, Zero);
 
     return max
     (
-        lRAS
-      - fd(mag(gradU))
-       *max
-        (
-            lRAS - lLES,
-            dimensionedScalar(dimLength, Zero)
-        ),
+        lRAS - fd(mag(gradU))*max(lRAS - lLES, l0),
         dimensionedScalar("small", dimLength, SMALL)
     );
 }
@@ -131,12 +141,19 @@ SpalartAllmarasDDES<BasicTurbulenceModel>::SpalartAllmarasDDES
 
     Cd1_
     (
-        dimensioned<scalar>::getOrAddToDict
-        (
-            "Cd1",
-            this->coeffDict_,
-            8
-        )
+        this->useSigma_
+          ? dimensioned<scalar>::getOrAddToDict
+            (
+                "Cd1Sigma",
+                this->coeffDict_,
+                10
+            )
+          : dimensioned<scalar>::getOrAddToDict
+            (
+                "Cd1",
+                this->coeffDict_,
+                8
+            )
     ),
     Cd2_
     (
@@ -172,6 +189,13 @@ bool SpalartAllmarasDDES<BasicTurbulenceModel>::read()
 }
 
 
+template<class BasicTurbulenceModel>
+tmp<volScalarField> SpalartAllmarasDDES<BasicTurbulenceModel>::fd() const
+{
+    return fd(mag(fvc::grad(this->U_)));
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace LESModels

src/TurbulenceModels/turbulenceModels/DES/SpalartAllmarasDDES/SpalartAllmarasDDES.H

@@ -6,7 +6,8 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2022 Upstream CFD GmbH
+    Copyright (C) 2019-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -32,15 +33,16 @@ Group
 
 Description
     SpalartAllmaras DDES turbulence model for incompressible and compressible
-    flows
+    flows.
 
     Reference:
     \verbatim
-        Spalart, P. R., Deck, S., Shur, M. L., Squires, K. D., Strelets, M. K.,
-        & Travin, A. (2006).
-        A new version of detached-eddy simulation, resistant to ambiguous grid
-        densities.
+        Spalart, P. R., Deck, S., Shur, M. L., Squires,
+        K. D., Strelets, M. K., & Travin, A. (2006).
+        A new version of detached-eddy simulation,
+        resistant to ambiguous grid densities.
         Theoretical and computational fluid dynamics, 20(3), 181-195.
+        DOI:10.1007/s00162-006-0015-0
     \endverbatim
 
 SourceFiles
@@ -48,10 +50,10 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef SpalartAllmarasDDES_H
-#define SpalartAllmarasDDES_H
+#ifndef Foam_SpalartAllmarasDDES_H
+#define Foam_SpalartAllmarasDDES_H
 
-#include "SpalartAllmarasDES.H"
+#include "SpalartAllmarasBase.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -71,10 +73,9 @@ class SpalartAllmarasDDES
 {
     // Private Member Functions
 
+        //- Return the shielding function
         tmp<volScalarField> fd(const volScalarField& magGradU) const;
 
-        tmp<volScalarField> rd(const volScalarField& magGradU) const;
-
         //- No copy construct
         SpalartAllmarasDDES(const SpalartAllmarasDDES&) = delete;
 
@@ -84,7 +85,7 @@ class SpalartAllmarasDDES
 
 protected:
 
-    // Protected data
+    // Protected Data
 
         // Model coefficients
 
@@ -94,7 +95,16 @@ protected:
 
     // Protected Member Functions
 
-        //- Length scale
+        //- Return the production term
+        virtual tmp<volScalarField> Stilda
+        (
+            const volScalarField& chi,
+            const volScalarField& fv1,
+            const volTensorField& gradU,
+            const volScalarField& dTilda
+        ) const;
+
+        //- Return the length scale
         virtual tmp<volScalarField> dTilda
         (
             const volScalarField& chi,
@@ -138,6 +148,9 @@ public:
 
         //- Read from dictionary
         virtual bool read();
+
+        //- Return the shielding function
+        virtual tmp<volScalarField> fd() const;
 };
 
 

src/TurbulenceModels/turbulenceModels/DES/SpalartAllmarasDES/SpalartAllmarasDES.C

@@ -6,7 +6,8 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2016-2020 OpenCFD Ltd.
+    Copyright (C) 2022 Upstream CFD GmbH
+    Copyright (C) 2016-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -38,120 +39,6 @@ namespace LESModels
 
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
-template<class BasicTurbulenceModel>
-tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::chi() const
-{
-    return nuTilda_/this->nu();
-}
-
-
-template<class BasicTurbulenceModel>
-tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::fv1
-(
-    const volScalarField& chi
-) const
-{
-    const volScalarField chi3("chi3", pow3(chi));
-    return chi3/(chi3 + pow3(Cv1_));
-}
-
-
-template<class BasicTurbulenceModel>
-tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::fv2
-(
-    const volScalarField& chi,
-    const volScalarField& fv1
-) const
-{
-    return 1.0 - chi/(1.0 + chi*fv1);
-}
-
-
-template<class BasicTurbulenceModel>
-tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::ft2
-(
-    const volScalarField& chi
-) const
-{
-    return Ct3_*exp(-Ct4_*sqr(chi));
-}
-
-
-template<class BasicTurbulenceModel>
-tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::Omega
-(
-    const volTensorField& gradU
-) const
-{
-    return sqrt(2.0)*mag(skew(gradU));
-}
-
-
-template<class BasicTurbulenceModel>
-tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::Stilda
-(
-    const volScalarField& chi,
-    const volScalarField& fv1,
-    const volScalarField& Omega,
-    const volScalarField& dTilda
-) const
-{
-    return
-    (
-        max
-        (
-            Omega
-          + fv2(chi, fv1)*nuTilda_/sqr(kappa_*dTilda),
-            Cs_*Omega
-        )
-    );
-}
-
-
-template<class BasicTurbulenceModel>
-tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::r
-(
-    const volScalarField& nur,
-    const volScalarField& Omega,
-    const volScalarField& dTilda
-) const
-{
-    tmp<volScalarField> tr
-    (
-        min
-        (
-            nur
-           /(
-                max
-                (
-                    Omega,
-                    dimensionedScalar("SMALL", Omega.dimensions(), SMALL)
-                )
-                *sqr(kappa_*dTilda)
-            ),
-            scalar(10)
-        )
-    );
-    tr.ref().boundaryFieldRef() == 0.0;
-
-    return tr;
-}
-
-
-template<class BasicTurbulenceModel>
-tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::fw
-(
-    const volScalarField& Omega,
-    const volScalarField& dTilda
-) const
-{
-    const volScalarField r(this->r(nuTilda_, Omega, dTilda));
-    const volScalarField g(r + Cw2_*(pow6(r) - r));
-
-    return g*pow((1 + pow6(Cw3_))/(pow6(g) + pow6(Cw3_)), 1.0/6.0);
-}
-
-
 template<class BasicTurbulenceModel>
 tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::psi
 (
@@ -159,26 +46,23 @@ tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::psi
     const volScalarField& fv1
 ) const
 {
-    tmp<volScalarField> tpsi
+    auto tpsi = tmp<volScalarField>::New
     (
-        new volScalarField
+        IOobject
         (
-            IOobject
-            (
-                type() + ":psi",
-                this->time().timeName(),
-                this->mesh(),
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
+            IOobject::scopedName(type(), "psi"),
+            this->time().timeName(),
             this->mesh(),
-            dimensionedScalar("one", dimless, 1)
-        )
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        this->mesh(),
+        dimensionedScalar(dimless, Zero)
     );
 
     if (lowReCorrection_)
     {
-        volScalarField& psi = tpsi.ref();
+        auto& psi = tpsi.ref();
 
         const volScalarField fv2(this->fv2(chi, fv1));
         const volScalarField ft2(this->ft2(chi));
@@ -189,7 +73,8 @@ tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::psi
                 min
                 (
                     scalar(100),
-                    (1 - Cb1_/(Cw1_*sqr(kappa_)*fwStar_)*(ft2 + (1 - ft2)*fv2))
+                    (1 - this->Cb1_/(this->Cw1_*sqr(this->kappa_)*fwStar_)
+                   *(ft2 + (1 - ft2)*fv2))
                    /max(SMALL, (fv1*max(scalar(1e-10), 1 - ft2)))
                 )
             );
@@ -199,6 +84,57 @@ tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::psi
 }
 
 
+template<class BasicTurbulenceModel>
+tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::lengthScaleLES
+(
+    const volScalarField& chi,
+    const volScalarField& fv1
+) const
+{
+    return psi(chi, fv1)*CDES_*this->delta();
+}
+
+
+template<class BasicTurbulenceModel>
+tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::Stilda
+(
+    const volScalarField& chi,
+    const volScalarField& fv1,
+    const volTensorField& gradU,
+    const volScalarField& dTilda
+) const
+{
+    if (useSigma_)
+    {
+        const volScalarField& lRAS = this->y_;
+        const volScalarField fv2(this->fv2(chi, fv1));
+        const volScalarField lLES(this->lengthScaleLES(chi, fv1));
+        const volScalarField Omega(this->Omega(gradU));
+        const volScalarField Ssigma(this->Ssigma(gradU));
+        const volScalarField SsigmaDES
+        (
+            pos(dTilda - lRAS)*Omega + (scalar(1) - pos(dTilda - lRAS))*Ssigma
+        );
+
+        return
+            max
+            (
+                SsigmaDES + fv2*this->nuTilda_/sqr(this->kappa_*dTilda),
+                this->Cs_*SsigmaDES
+            );
+    }
+
+    return
+        SpalartAllmarasBase<DESModel<BasicTurbulenceModel>>::Stilda
+        (
+            chi,
+            fv1,
+            gradU,
+            dTilda
+        );
+}
+
+
 template<class BasicTurbulenceModel>
 tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::dTilda
 (
@@ -207,33 +143,15 @@ tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::dTilda
     const volTensorField& gradU
 ) const
 {
-    tmp<volScalarField> tdTilda(psi(chi, fv1)*CDES_*this->delta());
-    min(tdTilda.ref().ref(), tdTilda(), y_);
+    // Initialise with LES length scale
+    tmp<volScalarField> tdTilda = lengthScaleLES(chi, fv1);
+
+    // Take min vs wall distance
+    min(tdTilda.ref(), tdTilda(), this->y_);
     return tdTilda;
 }
 
 
-template<class BasicTurbulenceModel>
-void SpalartAllmarasDES<BasicTurbulenceModel>::correctNut
-(
-    const volScalarField& fv1
-)
-{
-    this->nut_ = nuTilda_*fv1;
-    this->nut_.correctBoundaryConditions();
-    fv::options::New(this->mesh_).correct(this->nut_);
-
-    BasicTurbulenceModel::correctNut();
-}
-
-
-template<class BasicTurbulenceModel>
-void SpalartAllmarasDES<BasicTurbulenceModel>::correctNut()
-{
-    correctNut(fv1(this->chi()));
-}
-
-
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 template<class BasicTurbulenceModel>
@@ -249,7 +167,7 @@ SpalartAllmarasDES<BasicTurbulenceModel>::SpalartAllmarasDES
     const word& type
 )
 :
-    DESModel<BasicTurbulenceModel>
+    SpalartAllmarasBase<DESModel<BasicTurbulenceModel>>
     (
         type,
         alpha,
@@ -261,77 +179,13 @@ SpalartAllmarasDES<BasicTurbulenceModel>::SpalartAllmarasDES
         propertiesName
     ),
 
-    sigmaNut_
+    useSigma_
     (
-        dimensioned<scalar>::getOrAddToDict
+        Switch::getOrAddToDict
         (
-            "sigmaNut",
+            "useSigma",
             this->coeffDict_,
-            0.66666
-        )
-    ),
-    kappa_
-    (
-        dimensioned<scalar>::getOrAddToDict
-        (
-            "kappa",
-            this->coeffDict_,
-            0.41
-        )
-    ),
-    Cb1_
-    (
-        dimensioned<scalar>::getOrAddToDict
-        (
-            "Cb1",
-            this->coeffDict_,
-            0.1355
-        )
-    ),
-    Cb2_
-    (
-        dimensioned<scalar>::getOrAddToDict
-        (
-            "Cb2",
-            this->coeffDict_,
-            0.622
-        )
-    ),
-    Cw1_(Cb1_/sqr(kappa_) + (1.0 + Cb2_)/sigmaNut_),
-    Cw2_
-    (
-        dimensioned<scalar>::getOrAddToDict
-        (
-            "Cw2",
-            this->coeffDict_,
-            0.3
-        )
-    ),
-    Cw3_
-    (
-        dimensioned<scalar>::getOrAddToDict
-        (
-            "Cw3",
-            this->coeffDict_,
-            2.0
-        )
-    ),
-    Cv1_
-    (
-        dimensioned<scalar>::getOrAddToDict
-        (
-            "Cv1",
-            this->coeffDict_,
-            7.1
-        )
-    ),
-    Cs_
-    (
-        dimensioned<scalar>::getOrAddToDict
-        (
-            "Cs",
-            this->coeffDict_,
-            0.3
+            false
         )
     ),
     CDES_
@@ -343,15 +197,6 @@ SpalartAllmarasDES<BasicTurbulenceModel>::SpalartAllmarasDES
             0.65
         )
     ),
-    ck_
-    (
-        dimensioned<scalar>::getOrAddToDict
-        (
-            "ck",
-            this->coeffDict_,
-            0.07
-        )
-    ),
     lowReCorrection_
     (
         Switch::getOrAddToDict
@@ -361,24 +206,6 @@ SpalartAllmarasDES<BasicTurbulenceModel>::SpalartAllmarasDES
             true
         )
     ),
-    Ct3_
-    (
-        dimensioned<scalar>::getOrAddToDict
-        (
-            "Ct3",
-            this->coeffDict_,
-            1.2
-        )
-    ),
-    Ct4_
-    (
-        dimensioned<scalar>::getOrAddToDict
-        (
-            "Ct4",
-            this->coeffDict_,
-            0.5
-        )
-    ),
     fwStar_
     (
         dimensioned<scalar>::getOrAddToDict
@@ -387,25 +214,19 @@ SpalartAllmarasDES<BasicTurbulenceModel>::SpalartAllmarasDES
             this->coeffDict_,
             0.424
         )
-    ),
-
-    nuTilda_
-    (
-        IOobject
-        (
-            "nuTilda",
-            this->runTime_.timeName(),
-            this->mesh_,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        this->mesh_
-    ),
-
-    y_(wallDist::New(this->mesh_).y())
+    )
 {
+    // Note: Ctrans coeff is model specific; for S-A = 67.7
+    this->Ctrans_ =
+        dimensioned<scalar>::getOrAddToDict("Ctrans", this->coeffDict_, 67.7);
+
     if (type == typeName)
     {
+        WarningInFunction
+            << "This model is not recommended and will be deprecated in future "
+            << "releases. Please consider using the DDES/IDDES versions instead"
+            << endl;
+
         this->printCoeffs(type);
     }
 }
@@ -416,25 +237,11 @@ SpalartAllmarasDES<BasicTurbulenceModel>::SpalartAllmarasDES
 template<class BasicTurbulenceModel>
 bool SpalartAllmarasDES<BasicTurbulenceModel>::read()
 {
-    if (DESModel<BasicTurbulenceModel>::read())
+    if (SpalartAllmarasBase<DESModel<BasicTurbulenceModel>>::read())
     {
-        sigmaNut_.readIfPresent(this->coeffDict());
-        kappa_.readIfPresent(*this);
-
-        Cb1_.readIfPresent(this->coeffDict());
-        Cb2_.readIfPresent(this->coeffDict());
-        Cw1_ = Cb1_/sqr(kappa_) + (1.0 + Cb2_)/sigmaNut_;
-        Cw2_.readIfPresent(this->coeffDict());
-        Cw3_.readIfPresent(this->coeffDict());
-        Cv1_.readIfPresent(this->coeffDict());
-        Cs_.readIfPresent(this->coeffDict());
-
+        useSigma_.readIfPresent("useSigma", this->coeffDict());
         CDES_.readIfPresent(this->coeffDict());
-        ck_.readIfPresent(this->coeffDict());
-
         lowReCorrection_.readIfPresent("lowReCorrection", this->coeffDict());
-        Ct3_.readIfPresent(this->coeffDict());
-        Ct4_.readIfPresent(this->coeffDict());
         fwStar_.readIfPresent(this->coeffDict());
 
         return true;
@@ -445,124 +252,24 @@ bool SpalartAllmarasDES<BasicTurbulenceModel>::read()
 
 
 template<class BasicTurbulenceModel>
-tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::
-DnuTildaEff() const
-{
-    return tmp<volScalarField>
-    (
-        new volScalarField("DnuTildaEff", (nuTilda_ + this->nu())/sigmaNut_)
-    );
-}
-
-
-template<class BasicTurbulenceModel>
-tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::k() const
+tmp<volScalarField>
+SpalartAllmarasDES<BasicTurbulenceModel>::LESRegion() const
 {
     const volScalarField chi(this->chi());
     const volScalarField fv1(this->fv1(chi));
 
-    tmp<volScalarField> tdTilda
+    return tmp<volScalarField>::New
     (
-        new volScalarField
+        IOobject
         (
-            IOobject
-            (
-                "dTilda",
-                this->mesh_.time().timeName(),
-                this->mesh_,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
+            IOobject::scopedName("DES", "LESRegion"),
+            this->mesh_.time().timeName(),
             this->mesh_,
-            dimensionedScalar(dimLength, Zero),
-            zeroGradientFvPatchField<vector>::typeName
-        )
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        neg(dTilda(chi, fv1, fvc::grad(this->U_)) - this->y_)
     );
-    volScalarField& dTildaF = tdTilda.ref();
-    dTildaF = dTilda(chi, fv1, fvc::grad(this->U_));
-    dTildaF.correctBoundaryConditions();
-
-    return sqr(this->nut()/ck_/dTildaF);
-}
-
-
-template<class BasicTurbulenceModel>
-tmp<volScalarField> SpalartAllmarasDES<BasicTurbulenceModel>::LESRegion() const
-{
-    const volScalarField chi(this->chi());
-    const volScalarField fv1(this->fv1(chi));
-
-    tmp<volScalarField> tLESRegion
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                "DES::LESRegion",
-                this->mesh_.time().timeName(),
-                this->mesh_,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            neg(dTilda(chi, fv1, fvc::grad(this->U_)) - y_)
-        )
-    );
-
-    return tLESRegion;
-}
-
-
-template<class BasicTurbulenceModel>
-void SpalartAllmarasDES<BasicTurbulenceModel>::correct()
-{
-    if (!this->turbulence_)
-    {
-        return;
-    }
-
-    // Local references
-    const alphaField& alpha = this->alpha_;
-    const rhoField& rho = this->rho_;
-    const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_;
-    const volVectorField& U = this->U_;
-    fv::options& fvOptions(fv::options::New(this->mesh_));
-
-    DESModel<BasicTurbulenceModel>::correct();
-
-    const volScalarField chi(this->chi());
-    const volScalarField fv1(this->fv1(chi));
-    const volScalarField ft2(this->ft2(chi));
-
-    tmp<volTensorField> tgradU = fvc::grad(U);
-    const volScalarField Omega(this->Omega(tgradU()));
-    const volScalarField dTilda(this->dTilda(chi, fv1, tgradU()));
-    const volScalarField Stilda(this->Stilda(chi, fv1, Omega, dTilda));
-
-    tmp<fvScalarMatrix> nuTildaEqn
-    (
-        fvm::ddt(alpha, rho, nuTilda_)
-      + fvm::div(alphaRhoPhi, nuTilda_)
-      - fvm::laplacian(alpha*rho*DnuTildaEff(), nuTilda_)
-      - Cb2_/sigmaNut_*alpha*rho*magSqr(fvc::grad(nuTilda_))
-     ==
-        Cb1_*alpha*rho*Stilda*nuTilda_*(scalar(1) - ft2)
-      - fvm::Sp
-        (
-            (Cw1_*fw(Stilda, dTilda) - Cb1_/sqr(kappa_)*ft2)
-           *alpha*rho*nuTilda_/sqr(dTilda),
-            nuTilda_
-        )
-      + fvOptions(alpha, rho, nuTilda_)
-    );
-
-    nuTildaEqn.ref().relax();
-    fvOptions.constrain(nuTildaEqn.ref());
-    solve(nuTildaEqn);
-    fvOptions.correct(nuTilda_);
-    bound(nuTilda_, dimensionedScalar(nuTilda_.dimensions(), Zero));
-    nuTilda_.correctBoundaryConditions();
-
-    correctNut();
 }
 
 

src/TurbulenceModels/turbulenceModels/DES/SpalartAllmarasDES/SpalartAllmarasDES.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2015-2019 OpenCFD Ltd.
+    Copyright (C) 2015-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -32,7 +32,7 @@ Group
 
 Description
     SpalartAllmarasDES DES turbulence model for incompressible and
-    compressible flows
+    compressible flows.
 
     Reference:
     \verbatim
@@ -44,11 +44,12 @@ Description
 
     Including the low Reynolds number correction documented in
     \verbatim
-        Spalart, P. R., Deck, S., Shur, M.L., Squires, K.D., Strelets, M.Kh,
-        Travin, A. (2006).
-        A new version of detached-eddy simulation, resistant to ambiguous grid
-        densities.
-        Theor. Comput. Fluid Dyn., 20, 181-195.
+        Spalart, P. R., Deck, S., Shur, M. L., Squires,
+        K. D., Strelets, M. K., & Travin, A. (2006).
+        A new version of detached-eddy simulation,
+        resistant to ambiguous grid densities.
+        Theoretical and computational fluid dynamics, 20(3), 181-195.
+        DOI:10.1007/s00162-006-0015-0
     \endverbatim
 
 Note
@@ -61,9 +62,10 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef SpalartAllmarasDES_H
-#define SpalartAllmarasDES_H
+#ifndef Foam_SpalartAllmarasDES_H
+#define Foam_SpalartAllmarasDES_H
 
+#include "SpalartAllmarasBase.H"
 #include "DESModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -80,7 +82,7 @@ namespace LESModels
 template<class BasicTurbulenceModel>
 class SpalartAllmarasDES
 :
-    public DESModel<BasicTurbulenceModel>
+    public SpalartAllmarasBase<DESModel<BasicTurbulenceModel>>
 {
     // Private Member Functions
 
@@ -93,86 +95,47 @@ class SpalartAllmarasDES
 
 protected:
 
-    // Protected data
+    // Protected Data
+
+        //- Switch to activate grey-area enhanced sigma-(D)DES
+        Switch useSigma_;
 
         // Model constants
 
-            dimensionedScalar sigmaNut_;
-            dimensionedScalar kappa_;
-
-            dimensionedScalar Cb1_;
-            dimensionedScalar Cb2_;
-            dimensionedScalar Cw1_;
-            dimensionedScalar Cw2_;
-            dimensionedScalar Cw3_;
-            dimensionedScalar Cv1_;
-            dimensionedScalar Cs_;
+            //- DES coefficient
             dimensionedScalar CDES_;
-            dimensionedScalar ck_;
 
-
-            // Low Reynolds number correction
-
-                Switch lowReCorrection_;
-                dimensionedScalar Ct3_;
-                dimensionedScalar Ct4_;
-                dimensionedScalar fwStar_;
-
-
-        // Fields
-
-            volScalarField nuTilda_;
-
-            //- Wall distance
-            //  Note: different to wall distance in parent RASModel
-            //  which is for near-wall cells only
-            const volScalarField& y_;
+            //- Flag for low Reynolds number correction
+            Switch lowReCorrection_;
+            dimensionedScalar fwStar_;
 
 
     // Protected Member Functions
 
-        tmp<volScalarField> chi() const;
-
-        tmp<volScalarField> fv1(const volScalarField& chi) const;
-
-        tmp<volScalarField> fv2
+        //- Return the low Reynolds number correction function
+        virtual tmp<volScalarField> psi
         (
             const volScalarField& chi,
             const volScalarField& fv1
         ) const;
 
-        tmp<volScalarField> ft2(const volScalarField& chi) const;
+        //- Return the LES length scale
+        virtual tmp<volScalarField> lengthScaleLES
+        (
+            const volScalarField& chi,
+            const volScalarField& fv1
+        ) const;
 
-        tmp<volScalarField> Omega(const volTensorField& gradU) const;
-
-        tmp<volScalarField> Stilda
+        //- Return the production term
+        virtual tmp<volScalarField> Stilda
         (
             const volScalarField& chi,
             const volScalarField& fv1,
-            const volScalarField& Omega,
+            const volTensorField& gradU,
             const volScalarField& dTilda
         ) const;
 
-        tmp<volScalarField> r
-        (
-            const volScalarField& nur,
-            const volScalarField& Omega,
-            const volScalarField& dTilda
-        ) const;
-
-        tmp<volScalarField> fw
-        (
-            const volScalarField& Omega,
-            const volScalarField& dTilda
-        ) const;
-
-        tmp<volScalarField> psi
-        (
-            const volScalarField& chi,
-            const volScalarField& fv1
-        ) const;
-
-        //- Length scale
+        //- Return the length scale
         virtual tmp<volScalarField> dTilda
         (
             const volScalarField& chi,
@@ -180,9 +143,6 @@ protected:
             const volTensorField& gradU
         ) const;
 
-        void correctNut(const volScalarField& fv1);
-        virtual void correctNut();
-
 
 public:
 
@@ -220,22 +180,8 @@ public:
         //- Re-read model coefficients if they have changed
         virtual bool read();
 
-        //- Return the effective diffusivity for nuTilda
-        tmp<volScalarField> DnuTildaEff() const;
-
-        //- Return SGS kinetic energy
-        virtual tmp<volScalarField> k() const;
-
-        tmp<volScalarField> nuTilda() const
-        {
-            return nuTilda_;
-        }
-
         //- Return the LES field indicator
         virtual tmp<volScalarField> LESRegion() const;
-
-        //- Correct nuTilda and related properties
-        virtual void correct();
 };
 
 

src/TurbulenceModels/turbulenceModels/DES/SpalartAllmarasIDDES/SpalartAllmarasIDDES.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2015 OpenFOAM Foundation
-    Copyright (C) 2015-2020 OpenCFD Ltd.
+    Copyright (C) 2015-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -64,7 +64,7 @@ tmp<volScalarField> SpalartAllmarasIDDES<BasicTurbulenceModel>::ft
     const volScalarField& magGradU
 ) const
 {
-    return tanh(pow3(sqr(Ct_)*rd(this->nut_, magGradU)));
+    return tanh(pow3(sqr(Ct_)*this->r(this->nut_, magGradU, this->y_)));
 }
 
 
@@ -74,36 +74,7 @@ tmp<volScalarField> SpalartAllmarasIDDES<BasicTurbulenceModel>::fl
     const volScalarField& magGradU
 ) const
 {
-    return tanh(pow(sqr(Cl_)*rd(this->nu(), magGradU), 10));
-}
-
-
-template<class BasicTurbulenceModel>
-tmp<volScalarField> SpalartAllmarasIDDES<BasicTurbulenceModel>::rd
-(
-    const volScalarField& nur,
-    const volScalarField& magGradU
-) const
-{
-    tmp<volScalarField> tr
-    (
-        min
-        (
-            nur
-           /(
-                max
-                (
-                    magGradU,
-                    dimensionedScalar("SMALL", magGradU.dimensions(), SMALL)
-                )
-               *sqr(this->kappa_*this->y_)
-            ),
-            scalar(10)
-        )
-    );
-    tr.ref().boundaryFieldRef() == 0.0;
-
-    return tr;
+    return tanh(pow(sqr(Cl_)*this->r(this->nu(), magGradU, this->y_), 10));
 }
 
 
@@ -113,7 +84,7 @@ tmp<volScalarField> SpalartAllmarasIDDES<BasicTurbulenceModel>::fdt
     const volScalarField& magGradU
 ) const
 {
-    return 1 - tanh(pow(Cdt1_*rd(this->nut_, magGradU), Cdt2_));
+    return 1 - tanh(pow(Cdt1_*this->r(this->nut_, magGradU, this->y_), Cdt2_));
 }
 
 
@@ -130,31 +101,35 @@ tmp<volScalarField> SpalartAllmarasIDDES<BasicTurbulenceModel>::dTilda
     const volScalarField magGradU(mag(gradU));
     const volScalarField psi(this->psi(chi, fv1));
 
-    const volScalarField& lRAS(this->y_);
+    const volScalarField& lRAS = this->y_;
     const volScalarField lLES(psi*this->CDES_*this->delta());
 
     const volScalarField alpha(this->alpha());
     const volScalarField expTerm(exp(sqr(alpha)));
 
     tmp<volScalarField> fB = min(2*pow(expTerm, -9.0), scalar(1));
-    tmp<volScalarField> fe1 =
-        2*(pos0(alpha)*pow(expTerm, -11.09) + neg(alpha)*pow(expTerm, -9.0));
-    tmp<volScalarField> fe2 = 1 - max(ft(magGradU), fl(magGradU));
-    tmp<volScalarField> fe = max(fe1 - 1, scalar(0))*psi*fe2;
-
     const volScalarField fdTilda(max(1 - fdt(magGradU), fB));
 
-    // Simplified formulation from Gritskevich et al. paper (2011) where fe = 0
-    // return max
-    // (
-    //     fdTilda*lRAS + (1 - fdTilda)*lLES,
-    //     dimensionedScalar("SMALL", dimLength, SMALL)
-    // );
+    if (fe_)
+    {
+        tmp<volScalarField> fe1 =
+            2*(pos0(alpha)*pow(expTerm, -11.09) + neg(alpha)*pow(expTerm, -9.));
+        tmp<volScalarField> fe2 = 1 - max(ft(magGradU), fl(magGradU));
+        tmp<volScalarField> fe = max(fe1 - 1, scalar(0))*psi*fe2;
 
-    // Original formulation from Shur et al. paper (2008)
+        // Original formulation from Shur et al. paper (2008)
+        return max
+        (
+            fdTilda*(1 + fe)*lRAS + (1 - fdTilda)*lLES,
+            dimensionedScalar("SMALL", dimLength, SMALL)
+        );
+    }
+
+
+    // Simplified formulation from Gritskevich et al. paper (2011) where fe = 0
     return max
     (
-        fdTilda*(1 + fe)*lRAS + (1 - fdTilda)*lLES,
+        fdTilda*lRAS + (1 - fdTilda)*lLES,
         dimensionedScalar("SMALL", dimLength, SMALL)
     );
 }
@@ -223,6 +198,16 @@ SpalartAllmarasIDDES<BasicTurbulenceModel>::SpalartAllmarasIDDES
             1.63
         )
     ),
+    fe_
+    (
+        Switch::getOrAddToDict
+        (
+            "fe",
+            this->coeffDict_,
+            true
+        )
+    ),
+
     IDDESDelta_(setDelta())
 {
     if (type == typeName)
@@ -251,6 +236,17 @@ bool SpalartAllmarasIDDES<BasicTurbulenceModel>::read()
 }
 
 
+template<class BasicTurbulenceModel>
+tmp<volScalarField> SpalartAllmarasIDDES<BasicTurbulenceModel>::fd() const
+{
+    const volScalarField alpha(this->alpha());
+    const volScalarField expTerm(exp(sqr(alpha)));
+
+    tmp<volScalarField> fB = min(2*pow(expTerm, -9.0), scalar(1));
+    return max(1 - fdt(mag(fvc::grad(this->U_))), fB);
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace LESModels

src/TurbulenceModels/turbulenceModels/DES/SpalartAllmarasIDDES/SpalartAllmarasIDDES.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2011-2016 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -31,15 +31,16 @@ Group
     grpDESTurbulence
 
 Description
-    SpalartAllmaras IDDES turbulence model for incompressible and compressible
-    flows
+    SpalartAllmaras IDDES turbulence model
+    for incompressible and compressible flows.
 
     Reference:
     \verbatim
         Shur, M. L., Spalart, P. R., Strelets, M. K., & Travin, A. K. (2008).
-        A hybrid RANS-LES approach with delayed-DES and wall-modelled LES
-        capabilities.
-        International Journal of Heat and Fluid Flow, 29(6), 1638-1649.
+        A hybrid RANS-LES approach with delayed-DES
+        and wall-modelled LES capabilities.
+        International journal of heat and fluid flow, 29(6), 1638-1649.
+        DOI:10.1016/j.ijheatfluidflow.2008.07.001
     \endverbatim
 
 SourceFiles
@@ -47,8 +48,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef SpalartAllmarasIDDES_H
-#define SpalartAllmarasIDDES_H
+#ifndef Foam_SpalartAllmarasIDDES_H
+#define Foam_SpalartAllmarasIDDES_H
 
 #include "SpalartAllmarasDES.H"
 #include "IDDESDelta.H"
@@ -61,7 +62,7 @@ namespace LESModels
 {
 
 /*---------------------------------------------------------------------------*\
-                        Class SpalartAllmarasIDDES Declaration
+                    Class SpalartAllmarasIDDES Declaration
 \*---------------------------------------------------------------------------*/
 
 template<class BasicTurbulenceModel>
@@ -75,14 +76,10 @@ class SpalartAllmarasIDDES
         const IDDESDelta& setDelta() const;
 
         tmp<volScalarField> alpha() const;
-        tmp<volScalarField> ft(const volScalarField& magGradU) const;
-        tmp<volScalarField> fl(const volScalarField& magGradU) const;
 
-        tmp<volScalarField> rd
-        (
-            const volScalarField& nur,
-            const volScalarField& magGradU
-        ) const;
+        tmp<volScalarField> ft(const volScalarField& magGradU) const;
+
+        tmp<volScalarField> fl(const volScalarField& magGradU) const;
 
         //- Delay function
         tmp<volScalarField> fdt(const volScalarField& magGradU) const;
@@ -96,7 +93,7 @@ class SpalartAllmarasIDDES
 
 protected:
 
-    // Protected data
+    // Protected Data
 
         // Model coefficients
 
@@ -104,14 +101,16 @@ protected:
             dimensionedScalar Cdt2_;
             dimensionedScalar Cl_;
             dimensionedScalar Ct_;
+            Switch fe_;
 
-        // Fields
 
-            const IDDESDelta& IDDESDelta_;
+        //- IDDES delta
+        const IDDESDelta& IDDESDelta_;
+
 
     // Protected Member Functions
 
-        //- Length scale
+        //- Return the length scale
         virtual tmp<volScalarField> dTilda
         (
             const volScalarField& chi,
@@ -155,6 +154,9 @@ public:
 
         //- Re-read model coefficients if they have changed
         virtual bool read();
+
+        //- Return the shielding function
+        virtual tmp<volScalarField> fd() const;
 };
 
 

src/TurbulenceModels/turbulenceModels/DES/kOmegaSSTDDES/kOmegaSSTDDES.C

@@ -6,7 +6,8 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2015 OpenFOAM Foundation
-    Copyright (C) 2016-2020 OpenCFD Ltd.
+    Copyright (C) 2022 Upstream CFD GmbH
+    Copyright (C) 2016-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -37,46 +38,44 @@ namespace LESModels
 
 // * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
 
-template<class BasicTurbulenceModel>
-tmp<volScalarField> kOmegaSSTDDES<BasicTurbulenceModel>::rd
-(
-    const volScalarField& magGradU
-) const
-{
-    tmp<volScalarField> tr
-    (
-        min
-        (
-            this->nuEff()
-           /(
-                max
-                (
-                    magGradU,
-                    dimensionedScalar("SMALL", magGradU.dimensions(), SMALL)
-                )
-                *sqr(this->kappa_*this->y_)
-            ),
-            scalar(10)
-        )
-    );
-    tr.ref().boundaryFieldRef() == 0.0;
-
-    return tr;
-}
-
-
 template<class BasicTurbulenceModel>
 tmp<volScalarField> kOmegaSSTDDES<BasicTurbulenceModel>::fd
 (
     const volScalarField& magGradU
 ) const
 {
-    return 1 - tanh(pow(Cd1_*rd(magGradU), Cd2_));
+    return 1 - tanh(pow(Cd1_*this->r(this->nuEff(), magGradU), Cd2_));
 }
 
 
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
+template<class BasicTurbulenceModel>
+tmp<volScalarField> kOmegaSSTDDES<BasicTurbulenceModel>::S2
+(
+    const volScalarField& F1,
+    const volTensorField& gradU
+) const
+{
+    tmp<volScalarField> tS2 =
+        kOmegaSSTBase<DESModel<BasicTurbulenceModel>>::S2(F1, gradU);
+
+    if (this->useSigma_)
+    {
+        volScalarField& S2 = tS2.ref();
+        const volScalarField CDES(this->CDES(F1));
+        const volScalarField Ssigma(this->Ssigma(gradU));
+
+        S2 -=
+            fd(mag(gradU))
+           *pos(this->lengthScaleRAS() - this->lengthScaleLES(CDES))
+           *(S2 - sqr(Ssigma));
+    }
+
+    return tS2;
+}
+
+
 template<class BasicTurbulenceModel>
 tmp<volScalarField> kOmegaSSTDDES<BasicTurbulenceModel>::dTilda
 (
@@ -84,25 +83,46 @@ tmp<volScalarField> kOmegaSSTDDES<BasicTurbulenceModel>::dTilda
     const volScalarField& CDES
 ) const
 {
-    const volScalarField& k = this->k_;
-    const volScalarField& omega = this->omega_;
-
-    const volScalarField lRAS(sqrt(k)/(this->betaStar_*omega));
-    const volScalarField lLES(CDES*this->delta());
+    const volScalarField lRAS(this->lengthScaleRAS());
+    const volScalarField lLES(this->lengthScaleLES(CDES));
+    const dimensionedScalar l0(dimLength, Zero);
 
     return max
     (
-        lRAS
-      - fd(magGradU)
-       *max
-        (
-            lRAS - lLES,
-            dimensionedScalar(dimLength, Zero)
-        ),
+        lRAS - fd(magGradU)*max(lRAS - lLES, l0),
         dimensionedScalar("small", dimLength, SMALL)
     );
 }
 
+template<class BasicTurbulenceModel>
+tmp<volScalarField::Internal> kOmegaSSTDDES<BasicTurbulenceModel>::GbyNu0
+(
+    const volTensorField& gradU,
+    const volScalarField& F1,
+    const volScalarField& S2
+) const
+{
+    if (this->useSigma_)
+    {
+        return S2();
+    }
+
+    return
+        kOmegaSSTBase<DESModel<BasicTurbulenceModel>>::GbyNu0(gradU, F1, S2);
+}
+
+
+template<class BasicTurbulenceModel>
+tmp<volScalarField::Internal> kOmegaSSTDDES<BasicTurbulenceModel>::GbyNu
+(
+    const volScalarField::Internal& GbyNu0,
+    const volScalarField::Internal& F2,
+    const volScalarField::Internal& S2
+) const
+{
+    return GbyNu0; // Unlimited
+}
+
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -133,12 +153,19 @@ kOmegaSSTDDES<BasicTurbulenceModel>::kOmegaSSTDDES
 
     Cd1_
     (
-        dimensioned<scalar>::getOrAddToDict
-        (
-            "Cd1",
-            this->coeffDict_,
-            20
-        )
+        this->useSigma_
+          ? dimensioned<scalar>::getOrAddToDict
+            (
+                "Cd1Sigma",
+                this->coeffDict_,
+                22
+            )
+          : dimensioned<scalar>::getOrAddToDict
+            (
+                "Cd1",
+                this->coeffDict_,
+                20
+            )
     ),
     Cd2_
     (
@@ -174,6 +201,13 @@ bool kOmegaSSTDDES<BasicTurbulenceModel>::read()
 }
 
 
+template<class BasicTurbulenceModel>
+tmp<volScalarField> kOmegaSSTDDES<BasicTurbulenceModel>::fd() const
+{
+    return fd(mag(fvc::grad(this->U_)));
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace LESModels

src/TurbulenceModels/turbulenceModels/DES/kOmegaSSTDDES/kOmegaSSTDDES.H

@@ -6,7 +6,8 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2015 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2022 Upstream CFD GmbH
+    Copyright (C) 2019, 2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -31,14 +32,16 @@ Group
     grpDESTurbulence
 
 Description
-    k-omega-SST DDES turbulence model for incompressible and compressible flows
+    k-omega-SST DDES turbulence model for incompressible and compressible flows.
 
     Reference:
     \verbatim
-        Gritskevich, M.S., Garbaruk, A.V., Schuetze, J., Menter, F.R. (2011)
-        Development of DDES and IDDES Formulations for the k-omega
-        Shear Stress Transport Model, Flow, Turbulence and Combustion,
-        pp. 1-19
+        Gritskevich, M. S., Garbaruk, A. V.,
+        Schütze, J., & Menter, F. R. (2012).
+        Development of DDES and IDDES formulations for
+        the k-ω shear stress transport model.
+        Flow, turbulence and combustion, 88(3), 431-449.
+        DOI:10.1007/s10494-011-9378-4
     \endverbatim
 
 SourceFiles
@@ -46,8 +49,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef kOmegaSSTDDES_H
-#define kOmegaSSTDDES_H
+#ifndef Foam_kOmegaSSTDDES_H
+#define Foam_kOmegaSSTDDES_H
 
 #include "kOmegaSSTDES.H"
 
@@ -69,10 +72,9 @@ class kOmegaSSTDDES
 {
     // Private Member Functions
 
+        //- Return the shielding function field
         tmp<volScalarField> fd(const volScalarField& magGradU) const;
 
-        tmp<volScalarField> rd(const volScalarField& magGradU) const;
-
         //- No copy construct
         kOmegaSSTDDES(const kOmegaSSTDDES&) = delete;
 
@@ -82,7 +84,7 @@ class kOmegaSSTDDES
 
 protected:
 
-    // Protected data
+    // Protected Data
 
         // Model coefficients
 
@@ -92,6 +94,13 @@ protected:
 
     // Protected Member Functions
 
+        //- Return square of strain rate
+        virtual tmp<volScalarField> S2
+        (
+            const volScalarField& F1,
+            const volTensorField& gradU
+        ) const;
+
         //- Length scale
         virtual tmp<volScalarField> dTilda
         (
@@ -99,6 +108,21 @@ protected:
             const volScalarField& CDES
         ) const;
 
+        //- Return (G/nu)_0
+        virtual tmp<volScalarField::Internal> GbyNu0
+        (
+            const volTensorField& gradU,
+            const volScalarField& F1,
+            const volScalarField& S2
+        ) const;
+
+        //- Return G/nu
+        virtual tmp<volScalarField::Internal> GbyNu
+        (
+            const volScalarField::Internal& GbyNu0,
+            const volScalarField::Internal& F2,
+            const volScalarField::Internal& S2
+        ) const;
 
 public:
 
@@ -135,6 +159,9 @@ public:
 
         //- Re-read model coefficients if they have changed
         virtual bool read();
+
+        //- Return the shielding function
+        virtual tmp<volScalarField> fd() const;
 };
 
 

src/TurbulenceModels/turbulenceModels/DES/kOmegaSSTDES/kOmegaSSTDES.C

@@ -6,7 +6,8 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2015 OpenFOAM Foundation
-    Copyright (C) 2016-2020 OpenCFD Ltd.
+    Copyright (C) 2022 Upstream CFD GmbH
+    Copyright (C) 2016-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -55,6 +56,51 @@ void kOmegaSSTDES<BasicTurbulenceModel>::correctNut()
 }
 
 
+template<class BasicTurbulenceModel>
+tmp<volScalarField> kOmegaSSTDES<BasicTurbulenceModel>::r
+(
+    const volScalarField& nur,
+    const volScalarField& magGradU
+) const
+{
+    const dimensionedScalar eps(magGradU.dimensions(), SMALL);
+
+    tmp<volScalarField> tr =
+        min(nur/(max(magGradU, eps)*sqr(this->kappa_*this->y_)), scalar(10));
+
+    tr.ref().boundaryFieldRef() == 0;
+
+    return tr;
+}
+
+
+template<class BasicTurbulenceModel>
+tmp<volScalarField> kOmegaSSTDES<BasicTurbulenceModel>::S2
+(
+    const volScalarField& F1,
+    const volTensorField& gradU
+) const
+{
+    tmp<volScalarField> tS2 =
+        kOmegaSSTBase<DESModel<BasicTurbulenceModel>>::S2(F1, gradU);
+
+    if (this->useSigma_)
+    {
+        volScalarField& S2 = tS2.ref();
+        const volScalarField CDES(this->CDES(F1));
+        const volScalarField dTilda(this->dTilda(mag(gradU), CDES));
+        const volScalarField lengthScaleRAS(this->lengthScaleRAS());
+        const volScalarField Ssigma(this->Ssigma(gradU));
+
+        S2 =
+            pos(dTilda - lengthScaleRAS)*S2
+          + (scalar(1) - pos(dTilda - lengthScaleRAS))*sqr(Ssigma);
+    }
+
+    return tS2;
+}
+
+
 template<class BasicTurbulenceModel>
 tmp<volScalarField> kOmegaSSTDES<BasicTurbulenceModel>::dTilda
 (
@@ -62,10 +108,7 @@ tmp<volScalarField> kOmegaSSTDES<BasicTurbulenceModel>::dTilda
     const volScalarField& CDES
 ) const
 {
-    const volScalarField& k = this->k_;
-    const volScalarField& omega = this->omega_;
-
-    return min(CDES*this->delta(), sqrt(k)/(this->betaStar_*omega));
+    return min(lengthScaleLES(CDES), lengthScaleRAS());
 }
 
 
@@ -81,6 +124,24 @@ tmp<volScalarField::Internal> kOmegaSSTDES<BasicTurbulenceModel>::epsilonByk
 }
 
 
+template<class BasicTurbulenceModel>
+tmp<volScalarField::Internal> kOmegaSSTDES<BasicTurbulenceModel>::GbyNu0
+(
+    const volTensorField& gradU,
+    const volScalarField& F1,
+    const volScalarField& S2
+) const
+{
+    if (this->useSigma_)
+    {
+        return S2();
+    }
+
+    return
+        kOmegaSSTBase<DESModel<BasicTurbulenceModel>>::GbyNu0(gradU, F1, S2);
+}
+
+
 template<class BasicTurbulenceModel>
 tmp<volScalarField::Internal> kOmegaSSTDES<BasicTurbulenceModel>::GbyNu
 (
@@ -120,6 +181,15 @@ kOmegaSSTDES<BasicTurbulenceModel>::kOmegaSSTDES
         propertiesName
     ),
 
+    useSigma_
+    (
+        Switch::getOrAddToDict
+        (
+            "useSigma",
+            this->coeffDict_,
+            false
+        )
+    ),
     kappa_
     (
         dimensioned<scalar>::getOrAddToDict
@@ -148,8 +218,17 @@ kOmegaSSTDES<BasicTurbulenceModel>::kOmegaSSTDES
         )
     )
 {
+    // Note: Ctrans coeff is model specific; for k-w = 60
+    this->Ctrans_ =
+        dimensioned<scalar>::getOrAddToDict("Ctrans", this->coeffDict_, 60.0);
+
     if (type == typeName)
     {
+        WarningInFunction
+            << "This model is not recommended and will be deprecated in future "
+            << "releases. Please consider using the DDES/IDDES versions instead"
+            << endl;
+
         this->printCoeffs(type);
     }
 }
@@ -162,6 +241,7 @@ bool kOmegaSSTDES<BasicTurbulenceModel>::read()
 {
     if (kOmegaSSTBase<DESModel<BasicTurbulenceModel>>::read())
     {
+        useSigma_.readIfPresent("useSigma", this->coeffDict());
         kappa_.readIfPresent(this->coeffDict());
         CDESkom_.readIfPresent(this->coeffDict());
         CDESkeps_.readIfPresent(this->coeffDict());
@@ -173,6 +253,28 @@ bool kOmegaSSTDES<BasicTurbulenceModel>::read()
 }
 
 
+template<class BasicTurbulenceModel>
+Foam::tmp<Foam::volScalarField>
+kOmegaSSTDES<BasicTurbulenceModel>::lengthScaleRAS() const
+{
+    const volScalarField& k = this->k_;
+    const volScalarField& omega = this->omega_;
+
+    return sqrt(k)/(this->betaStar_*omega);
+}
+
+
+template<class BasicTurbulenceModel>
+Foam::tmp<Foam::volScalarField>
+kOmegaSSTDES<BasicTurbulenceModel>::lengthScaleLES
+(
+    const volScalarField& CDES
+) const
+{
+    return CDES*this->delta();
+}
+
+
 template<class BasicTurbulenceModel>
 tmp<volScalarField> kOmegaSSTDES<BasicTurbulenceModel>::LESRegion() const
 {
@@ -187,31 +289,11 @@ tmp<volScalarField> kOmegaSSTDES<BasicTurbulenceModel>::LESRegion() const
 
     const volScalarField F1(this->F1(CDkOmega));
 
-    tmp<volScalarField> tLESRegion
+    return tmp<volScalarField>::New
     (
-        new volScalarField
-        (
-            IOobject
-            (
-                "DES::LESRegion",
-                this->mesh_.time().timeName(),
-                this->mesh_,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            neg
-            (
-                dTilda
-                (
-                    mag(fvc::grad(U)),
-                    F1*CDESkom_ + (1 - F1)*CDESkeps_
-                )
-              - sqrt(k)/(this->betaStar_*omega)
-            )
-        )
+        IOobject::scopedName("DES", "LESRegion"),
+        neg(dTilda(mag(fvc::grad(U)), CDES(F1)) - lengthScaleRAS())
     );
-
-    return tLESRegion;
 }
 
 

src/TurbulenceModels/turbulenceModels/DES/kOmegaSSTDES/kOmegaSSTDES.H

@@ -31,13 +31,13 @@ Group
     grpDESTurbulence
 
 Description
-    k-omega-SST DES turbulence model for incompressible and compressible flows
+    k-omega-SST DES turbulence model for incompressible and compressible flows.
 
     Reference:
     \verbatim
-        Strelets, M. (2001)
-        Detached Eddy Simulation of Massively Separated Flows,
-        39th AIAA Aerospace Sciences Meeting and Exhibit, Reno, NV
+        Strelets, M. (2001).
+        Detached Eddy Simulation of Massively Separated Flows.
+        39th AIAA Aerospace Sciences Meeting and Exhibit, Reno, NV.
     \endverbatim
 
 Note
@@ -51,8 +51,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef kOmegaSSTDES_H
-#define kOmegaSSTDES_H
+#ifndef Foam_kOmegaSSTDES_H
+#define Foam_kOmegaSSTDES_H
 
 #include "DESModel.H"
 #include "kOmegaSSTBase.H"
@@ -65,7 +65,7 @@ namespace LESModels
 {
 
 /*---------------------------------------------------------------------------*\
-                        class kOmegaSSTDES Declaration
+                        Class kOmegaSSTDES Declaration
 \*---------------------------------------------------------------------------*/
 
 template<class BasicTurbulenceModel>
@@ -84,11 +84,16 @@ class kOmegaSSTDES
 
 protected:
 
-    // Protected data
+    // Protected Data
+
+        //- Switch to activate grey-area enhanced sigma-(D)DES
+        Switch useSigma_;
 
         // Model coefficients
 
             dimensionedScalar kappa_;
+
+            //- DES coefficients
             dimensionedScalar CDESkom_;
             dimensionedScalar CDESkeps_;
 
@@ -104,7 +109,20 @@ protected:
         virtual void correctNut(const volScalarField& S2);
         virtual void correctNut();
 
-        //- Length scale
+        tmp<volScalarField> r
+        (
+            const volScalarField& nur,
+            const volScalarField& magGradU
+        ) const;
+
+        //- Return square of strain rate
+        virtual tmp<volScalarField> S2
+        (
+            const volScalarField& F1,
+            const volTensorField& gradU
+        ) const;
+
+        //- Return length scale
         virtual tmp<volScalarField> dTilda
         (
             const volScalarField& magGradU,
@@ -118,6 +136,14 @@ protected:
             const volTensorField& gradU
         ) const;
 
+        //- Return (G/nu)_0
+        virtual tmp<volScalarField::Internal> GbyNu0
+        (
+            const volTensorField& gradU,
+            const volScalarField& F1,
+            const volScalarField& S2
+        ) const;
+
         //- Return G/nu
         virtual tmp<volScalarField::Internal> GbyNu
         (
@@ -163,6 +189,15 @@ public:
         //- Re-read model coefficients if they have changed
         virtual bool read();
 
+        //- RAS length scale
+        virtual tmp<volScalarField> lengthScaleRAS() const;
+
+        //- LES length scale
+        virtual tmp<volScalarField> lengthScaleLES
+        (
+            const volScalarField& CDES
+        ) const;
+
         //- Return the LES field indicator
         virtual tmp<volScalarField> LESRegion() const;
 };

src/TurbulenceModels/turbulenceModels/DES/kOmegaSSTIDDES/kOmegaSSTIDDES.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2015 OpenFOAM Foundation
-    Copyright (C) 2015-2020 OpenCFD Ltd.
+    Copyright (C) 2015-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -64,7 +64,7 @@ tmp<volScalarField> kOmegaSSTIDDES<BasicTurbulenceModel>::ft
     const volScalarField& magGradU
 ) const
 {
-    return tanh(pow3(sqr(Ct_)*rd(this->nut_, magGradU)));
+    return tanh(pow3(sqr(Ct_)*this->r(this->nut_, magGradU)));
 }
 
 
@@ -74,51 +74,22 @@ tmp<volScalarField> kOmegaSSTIDDES<BasicTurbulenceModel>::fl
     const volScalarField& magGradU
 ) const
 {
-    return tanh(pow(sqr(Cl_)*rd(this->nu(), magGradU), 10));
+    return tanh(pow(sqr(Cl_)*this->r(this->nu(), magGradU), 10));
 }
 
 
-template<class BasicTurbulenceModel>
-tmp<volScalarField> kOmegaSSTIDDES<BasicTurbulenceModel>::rd
-(
-    const volScalarField& nur,
-    const volScalarField& magGradU
-) const
-{
-    tmp<volScalarField> tr
-    (
-        min
-        (
-            nur
-           /(
-                max
-                (
-                    magGradU,
-                    dimensionedScalar("SMALL", magGradU.dimensions(), SMALL)
-                )
-                *sqr(this->kappa_*this->y_)
-            ),
-            scalar(10)
-        )
-    );
-    tr.ref().boundaryFieldRef() == 0.0;
-
-    return tr;
-}
-
-
-// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
-
 template<class BasicTurbulenceModel>
 tmp<volScalarField> kOmegaSSTIDDES<BasicTurbulenceModel>::fdt
 (
     const volScalarField& magGradU
 ) const
 {
-    return 1 - tanh(pow(Cdt1_*rd(this->nut_, magGradU), Cdt2_));
+    return 1 - tanh(pow(Cdt1_*this->r(this->nut_, magGradU), Cdt2_));
 }
 
 
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
 template<class BasicTurbulenceModel>
 tmp<volScalarField> kOmegaSSTIDDES<BasicTurbulenceModel>::dTilda
 (
@@ -126,35 +97,35 @@ tmp<volScalarField> kOmegaSSTIDDES<BasicTurbulenceModel>::dTilda
     const volScalarField& CDES
 ) const
 {
-    const volScalarField& k = this->k_;
-    const volScalarField& omega = this->omega_;
-
-    const volScalarField lRAS(sqrt(k)/(this->betaStar_*omega));
-    const volScalarField lLES(CDES*this->delta());
+    const volScalarField lRAS(this->lengthScaleRAS());
+    const volScalarField lLES(this->lengthScaleLES(CDES));
 
     const volScalarField alpha(this->alpha());
     const volScalarField expTerm(exp(sqr(alpha)));
 
     tmp<volScalarField> fB = min(2*pow(expTerm, -9.0), scalar(1));
-    tmp<volScalarField> fe1 =
-        2*(pos0(alpha)*pow(expTerm, -11.09) + neg(alpha)*pow(expTerm, -9.0));
-    tmp<volScalarField> fe2 = 1 - max(ft(magGradU), fl(magGradU));
-    tmp<volScalarField> fe = max(fe1 - 1, scalar(0))*fe2;
-
     const volScalarField fdTilda(max(1 - fdt(magGradU), fB));
 
-    // Simplified formulation from Gritskevich et al. paper (2011) where fe = 0
-    // return max
-    // (
-    //     fdTilda*lRAS + (1 - fdTilda)*lLES,
-    //     dimensionedScalar("SMALL", dimLength, SMALL)
-    // );
+    if (fe_)
+    {
+        tmp<volScalarField> fe1 =
+            2*(pos0(alpha)*pow(expTerm, -11.09) + neg(alpha)*pow(expTerm, -9.));
+        tmp<volScalarField> fe2 = 1 - max(ft(magGradU), fl(magGradU));
+        tmp<volScalarField> fe = max(fe1 - 1, scalar(0))*fe2;
 
-    // Original formulation from Shur et al. paper (2008)
+        // Original formulation from Shur et al. paper (2008)
+        return max
+        (
+            fdTilda*(1 + fe)*lRAS + (1 - fdTilda)*lLES,
+            dimensionedScalar(dimLength, SMALL)
+        );
+    }
+
+    // Simplified formulation from Gritskevich et al. paper (2011) where fe = 0
     return max
     (
-        fdTilda*(1 + fe)*lRAS + (1 - fdTilda)*lLES,
-        dimensionedScalar("SMALL", dimLength, SMALL)
+        fdTilda*lRAS + (1 - fdTilda)*lLES,
+        dimensionedScalar(dimLength, SMALL)
     );
 }
 
@@ -222,6 +193,16 @@ kOmegaSSTIDDES<BasicTurbulenceModel>::kOmegaSSTIDDES
             1.87
         )
     ),
+    fe_
+    (
+        Switch::getOrAddToDict
+        (
+            "fe",
+            this->coeffDict_,
+            true
+        )
+    ),
+
     IDDESDelta_(setDelta())
 {
     if (type == typeName)
@@ -250,6 +231,17 @@ bool kOmegaSSTIDDES<BasicTurbulenceModel>::read()
 }
 
 
+template<class BasicTurbulenceModel>
+tmp<volScalarField> kOmegaSSTIDDES<BasicTurbulenceModel>::fd() const
+{
+    const volScalarField alpha(this->alpha());
+    const volScalarField expTerm(exp(sqr(alpha)));
+
+    tmp<volScalarField> fB = min(2*pow(expTerm, -9.0), scalar(1));
+    return max(1 - fdt(mag(fvc::grad(this->U_))), fB);
+}
+
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace LESModels

src/TurbulenceModels/turbulenceModels/DES/kOmegaSSTIDDES/kOmegaSSTIDDES.H

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2015 OpenFOAM Foundation
-    Copyright (C) 2019 OpenCFD Ltd.
+    Copyright (C) 2019-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -31,15 +31,17 @@ Group
     grpDESTurbulence
 
 Description
-    k-omega-SST IDDES turbulence model for incompressible and compressible
-    flows
+    k-omega-SST IDDES turbulence model for
+    incompressible and compressible flows.
 
     Reference:
     \verbatim
-        Gritskevich, M.S., Garbaruk, A.V., Schuetze, J., Menter, F.R. (2011)
-        Development of DDES and IDDES Formulations for the k-omega
-        Shear Stress Transport Model, Flow, Turbulence and Combustion,
-        pp. 1-19
+        Gritskevich, M. S., Garbaruk, A. V.,
+        Schütze, J., & Menter, F. R. (2012).
+        Development of DDES and IDDES formulations for
+        the k-ω shear stress transport model.
+        Flow, turbulence and combustion, 88(3), 431-449.
+        DOI:10.1007/s10494-011-9378-4
     \endverbatim
 
 SourceFiles
@@ -47,8 +49,8 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef kOmegaSSTIDDES_H
-#define kOmegaSSTIDDES_H
+#ifndef Foam_kOmegaSSTIDDES_H
+#define Foam_kOmegaSSTIDDES_H
 
 #include "kOmegaSSTDES.H"
 
@@ -77,12 +79,6 @@ class kOmegaSSTIDDES
         tmp<volScalarField> ft(const volScalarField& magGradU) const;
         tmp<volScalarField> fl(const volScalarField& magGradU) const;
 
-        tmp<volScalarField> rd
-        (
-            const volScalarField& nur,
-            const volScalarField& magGradU
-        ) const;
-
         //- Delay function
         tmp<volScalarField> fdt(const volScalarField& magGradU) const;
 
@@ -95,7 +91,7 @@ class kOmegaSSTIDDES
 
 protected:
 
-    // Protected data
+    // Protected Data
 
         // Model coefficients
 
@@ -103,14 +99,16 @@ protected:
             dimensionedScalar Cdt2_;
             dimensionedScalar Cl_;
             dimensionedScalar Ct_;
+            Switch fe_;
 
-        // Fields
 
-            const IDDESDelta& IDDESDelta_;
+        //- IDDES delta
+        const IDDESDelta& IDDESDelta_;
+
 
     // Protected Member Functions
 
-        //- Length scale
+        //- Return the length scale
         virtual tmp<volScalarField> dTilda
         (
             const volScalarField& magGradU,
@@ -153,6 +151,9 @@ public:
 
         //- Re-read model coefficients if they have changed
         virtual bool read();
+
+        //- Return the shielding function
+        virtual tmp<volScalarField> fd() const;
 };
 
 

src/TurbulenceModels/turbulenceModels/LES/LESdeltas/DeltaOmegaTildeDelta/DeltaOmegaTildeDelta.C

@@ -0,0 +1,196 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2022 Upstream CFD GmbH
+    Copyright (C) 2022 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 "DeltaOmegaTildeDelta.H"
+#include "fvcCurl.H"
+#include "maxDeltaxyz.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace LESModels
+{
+    defineTypeNameAndDebug(DeltaOmegaTildeDelta, 0);
+    addToRunTimeSelectionTable(LESdelta, DeltaOmegaTildeDelta, dictionary);
+}
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::LESModels::DeltaOmegaTildeDelta::calcDelta()
+{
+    const fvMesh& mesh = turbulenceModel_.mesh();
+
+    const volVectorField& U0 = turbulenceModel_.U();
+    tmp<volVectorField> tvorticity = fvc::curl(U0);
+    const volVectorField& vorticity = tvorticity.cref();
+    const volVectorField nvecvort
+    (
+        vorticity
+      / (
+            max
+            (
+                mag(vorticity),
+                dimensionedScalar(dimless/dimTime, SMALL)
+            )
+        )
+    );
+    tvorticity.clear();
+
+    const cellList& cells = mesh.cells();
+    const vectorField& cellCentres = mesh.cellCentres();
+    const vectorField& faceCentres = mesh.faceCentres();
+
+    forAll(cells, celli)
+    {
+        const labelList& cFaces = cells[celli];
+        const point& cc = cellCentres[celli];
+        const vector& nv = nvecvort[celli];
+
+        scalar deltaMaxTmp = 0;
+
+        for (const label facei : cFaces)
+        {
+            const point& fc = faceCentres[facei];
+            const scalar tmp = 2.0*mag(nv ^ (fc - cc));
+
+            if (tmp > deltaMaxTmp)
+            {
+                deltaMaxTmp = tmp;
+            }
+        }
+
+        delta_[celli] = deltaCoeff_*deltaMaxTmp;
+    }
+
+    const label nD = mesh.nGeometricD();
+
+    if (nD == 2)
+    {
+        WarningInFunction
+            << "Case is 2D, LES is not strictly applicable" << nl
+            << endl;
+    }
+    else if (nD != 3)
+    {
+        FatalErrorInFunction
+            << "Case must be either 2D or 3D" << exit(FatalError);
+    }
+
+    // Handle coupled boundaries
+    delta_.correctBoundaryConditions();
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::LESModels::DeltaOmegaTildeDelta::DeltaOmegaTildeDelta
+(
+    const word& name,
+    const turbulenceModel& turbulence,
+    const dictionary& dict
+)
+:
+    LESdelta(name, turbulence),
+    hmaxPtr_(nullptr),
+    deltaCoeff_
+    (
+        dict.optionalSubDict(type() + "Coeffs").getOrDefault<scalar>
+        (
+            "deltaCoeff",
+            1.035
+        )
+    ),
+    requireUpdate_
+    (
+        dict.optionalSubDict(type() + "Coeffs").getOrDefault<bool>
+        (
+            "requireUpdate",
+            true
+        )
+    )
+{
+    if (dict.optionalSubDict(type() + "Coeffs").found("hmax"))
+    {
+        // User-defined hmax
+        hmaxPtr_ =
+            LESdelta::New
+            (
+                IOobject::groupName("hmax", turbulence.U().group()),
+                turbulence,
+                dict.optionalSubDict("hmaxCoeffs"),
+                "hmax"
+            );
+    }
+    else
+    {
+        Info<< "Employing " << maxDeltaxyz::typeName << " for hmax" << endl;
+
+        hmaxPtr_.reset
+        (
+            new maxDeltaxyz
+            (
+                IOobject::groupName("hmax", turbulence.U().group()),
+                turbulence,
+                dict.optionalSubDict("hmaxCoeffs")
+            )
+        );
+    }
+
+    calcDelta();
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::LESModels::DeltaOmegaTildeDelta::read(const dictionary& dict)
+{
+    const dictionary& coeffsDict = dict.optionalSubDict(type() + "Coeffs");
+
+    coeffsDict.readIfPresent<scalar>("deltaCoeff", deltaCoeff_);
+    coeffsDict.readIfPresent<bool>("requireUpdate", requireUpdate_);
+
+    calcDelta();
+}
+
+
+void Foam::LESModels::DeltaOmegaTildeDelta::correct()
+{
+    if (turbulenceModel_.mesh().changing() && requireUpdate_)
+    {
+        hmaxPtr_->correct();
+    }
+
+    calcDelta();
+}
+
+
+// ************************************************************************* //

src/TurbulenceModels/turbulenceModels/LES/LESdeltas/DeltaOmegaTildeDelta/DeltaOmegaTildeDelta.H

@@ -0,0 +1,142 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2022 Upstream CFD GmbH
+    Copyright (C) 2022 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::LESModels::DeltaOmegaTildeDelta
+
+Description
+    Delta formulation that accounts for the orientation of the vorticity
+    vector. In "2D-regions" (i.e. xy-plane), delta is of the order
+    max(delta_x,delta_y), so that the influence of delta_z is discarded.
+    This can help to accelerate the transition from RANS to LES in hybrid
+    RANS/LES simulations.
+
+    Reference:
+    \verbatim
+        Shur, M. L., Spalart, P. R., Strelets, M. K., & Travin, A. K. (2015).
+        An enhanced version of DES with rapid transition
+        from RANS to LES in separated flows.
+        Flow, turbulence and combustion, 95(4), 709-737.
+        DOI:10.1007/s10494-015-9618-0
+    \endverbatim
+
+SourceFiles
+    DeltaOmegaTildeDelta.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef LESModels_DeltaOmegaTildeDelta_H
+#define LESModels_DeltaOmegaTildeDelta_H
+
+#include "LESdelta.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace LESModels
+{
+
+/*---------------------------------------------------------------------------*\
+                    Class DeltaOmegaTildeDelta Declaration
+\*---------------------------------------------------------------------------*/
+
+class DeltaOmegaTildeDelta
+:
+    public LESdelta
+{
+    // Private Data
+
+        //- Run-time selectable delta for hmax
+        //  Defaults to the maxDeltaxyz model if not supplied
+        autoPtr<LESdelta> hmaxPtr_;
+
+        //- Model coefficient
+        scalar deltaCoeff_;
+
+        //- Flag to indicate whether hmax requires updating
+        bool requireUpdate_;
+
+
+    // Private Member Functions
+
+        //- Calculate the delta values
+        void calcDelta();
+
+        //- No copy construct
+        DeltaOmegaTildeDelta(const DeltaOmegaTildeDelta&) = delete;
+
+        //- No copy assignment
+        void operator=(const DeltaOmegaTildeDelta&) = delete;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("DeltaOmegaTilde");
+
+
+    // Constructors
+
+        //- Construct from name, turbulenceModel and dictionary
+        DeltaOmegaTildeDelta
+        (
+            const word& name,
+            const turbulenceModel& turbulence,
+            const dictionary&
+        );
+
+
+    //- Destructor
+    virtual ~DeltaOmegaTildeDelta() = default;
+
+
+    // Member Functions
+
+        //- Read the LESdelta dictionary
+        virtual void read(const dictionary&);
+
+        //- Return the hmax delta field
+        const volScalarField& hmax() const
+        {
+            return hmaxPtr_();
+        }
+
+        // Correct values
+        void correct();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace LESModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/TurbulenceModels/turbulenceModels/LES/LESdeltas/SLADelta/SLADelta.C

@@ -0,0 +1,412 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2022 Upstream CFD GmbH
+    Copyright (C) 2022 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 "SLADelta.H"
+#include "wallDist.H"
+#include "maxDeltaxyz.H"
+#include "fvcGrad.H"
+#include "fvcCurl.H"
+#include "addToRunTimeSelectionTable.H"
+
+#include "oneField.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+template<class Type, class WeightType = oneField>
+tmp<scalarField> sumNeighbours
+(
+    const GeometricField<Type, fvPatchField, volMesh>& field,
+    GeometricField<Type, fvPatchField, volMesh>& result,
+    const WeightType& weights = oneField()
+)
+{
+    const fvMesh& mesh = field.mesh();
+
+    result == dimensioned<Type>(field.dimensions(), Zero);
+
+    auto tscaling = tmp<scalarField>::New(mesh.nCells(), Zero);
+    auto& scaling = tscaling.ref();
+
+    const auto& faceNeighbour = mesh.faceNeighbour();
+    const auto& faceOwner = mesh.faceOwner();
+
+    forAll(faceNeighbour, i)
+    {
+        const label own = faceOwner[i];
+        const label nbr = faceNeighbour[i];
+
+        result[own] += field[nbr];
+        result[nbr] += field[own];
+
+        scaling[own] = scaling[own] + weights[nbr];
+        scaling[nbr] = scaling[nbr] + weights[own];
+    }
+
+    const auto& pbm = mesh.boundaryMesh();
+
+    for (const polyPatch& pp : pbm)
+    {
+        const auto& pf = field.boundaryField()[pp.index()];
+        const labelUList& faceCells = pp.faceCells();
+
+        if (pf.coupled())
+        {
+            const scalarField nbrField(pf.patchNeighbourField());
+
+            forAll(faceCells, facei)
+            {
+                const label celli = faceCells[facei];
+                result[celli] += nbrField[facei];
+                scaling[celli] = scaling[celli] + weights[celli];
+            }
+        }
+    }
+
+    return tscaling;
+}
+
+} // End namespace Foam
+
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace LESModels
+{
+    defineTypeNameAndDebug(SLADelta, 0);
+    addToRunTimeSelectionTable(LESdelta, SLADelta, dictionary);
+}
+}
+
+
+// * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
+
+void Foam::LESModels::SLADelta::calcDelta()
+{
+    const fvMesh& mesh = turbulenceModel_.mesh();
+
+    tmp<volScalarField> tnut = turbulenceModel_.nut();
+    const volScalarField& nut = tnut.cref();
+
+    tmp<volScalarField> tnu = turbulenceModel_.nu();
+    const volScalarField& nu = tnu.cref();
+
+    // Calculate vortex tilting measure, VTM
+    const volVectorField& U0 = turbulenceModel_.U();
+
+    tmp<volVectorField> tvorticity = fvc::curl(U0);
+    const volVectorField& vorticity = tvorticity.cref();
+
+    tmp<volSymmTensorField> tS = symm(fvc::grad(U0));
+    const volSymmTensorField& S = tS.cref();
+
+    const dimensionedScalar nuMin(nu.dimensions(), SMALL);
+    const dimensionedScalar eps(dimless/pow3(dimTime), SMALL);
+
+    volScalarField vtm
+    (
+        max(scalar(1), 0.2*nu/max(nut, nuMin))
+       *sqrt(6.0)*mag((S & vorticity)^vorticity)
+       /max(magSqr(vorticity)*sqrt(3*magSqr(S) - sqr(tr(S))), eps)
+    );
+    vtm.correctBoundaryConditions();
+    tS.clear();
+
+    const dimensionedScalar vortMin(dimless/dimTime, SMALL);
+    const volVectorField nVecVort(vorticity/(max(mag(vorticity), vortMin)));
+    tvorticity.clear();
+
+
+    // Calculate averaged VTM
+    volScalarField vtmAve("vtmAve", vtm);
+    tmp<scalarField> tweights = sumNeighbours(vtm, vtmAve);
+
+    // Add cell centre values
+    vtmAve += vtm;
+
+    // Weights normalisation (add 1 for cell centres)
+    vtmAve.primitiveFieldRef() /= tweights + 1;
+
+
+    // Calculate DDES shielding function, fd
+    const volScalarField& y = wallDist::New(mesh).y();
+
+    const dimensionedScalar magGradUeps(dimless/dimTime, SMALL);
+    const dimensionedScalar nuEps(nu.dimensions(), ROOTSMALL);
+
+    tmp<volScalarField> tmagGradU = mag(fvc::grad(U0));
+
+    tmp<volScalarField> trd =
+        min
+        (
+            (nut + nu)/(max(tmagGradU, magGradUeps)*sqr(kappa_*y) + nuEps),
+            scalar(10)
+        );
+    tnut.clear();
+    tnu.clear();
+
+    const volScalarField fd(1.0 - tanh(pow(Cd1_*trd, Cd2_)));
+
+
+    // Assemble delta
+    const cellList& cells = mesh.cells();
+    const vectorField& cellCentres = mesh.cellCentres();
+    const vectorField& faceCentres = mesh.faceCentres();
+
+    forAll(cells, celli)
+    {
+        const labelList& cFaces = cells[celli];
+        const point& cc = cellCentres[celli];
+        const vector& nv = nVecVort[celli];
+
+        scalar deltaMaxTmp = 0;
+
+        for (const label facei : cFaces)
+        {
+            const point& fc = faceCentres[facei];
+            const scalar tmp = 2.0*mag(nv ^ (fc - cc));
+
+            if (tmp > deltaMaxTmp)
+            {
+                deltaMaxTmp = tmp;
+            }
+        }
+
+        scalar FKH =
+            max
+            (
+                FKHmin_,
+                min
+                (
+                    FKHmax_,
+                    FKHmin_
+                  + (FKHmax_ - FKHmin_)*(vtmAve[celli] - a1_)/(a2_ - a1_)
+                )
+            );
+
+        if ((shielding_ == true) && (fd[celli] < (1.0 - epsilon_)))
+        {
+            FKH = 1;
+        }
+
+        delta_[celli] = deltaCoeff_*deltaMaxTmp*FKH;
+    }
+
+    const label nD = mesh.nGeometricD();
+
+    if (nD == 2)
+    {
+        WarningInFunction
+            << "Case is 2D, LES is not strictly applicable" << nl
+            << endl;
+    }
+    else if (nD != 3)
+    {
+        FatalErrorInFunction
+            << "Case must be either 2D or 3D" << exit(FatalError);
+    }
+
+    // Handle coupled boundaries
+    delta_.correctBoundaryConditions();
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+Foam::LESModels::SLADelta::SLADelta
+(
+    const word& name,
+    const turbulenceModel& turbulence,
+    const dictionary& dict
+)
+:
+    LESdelta(name, turbulence),
+    hmaxPtr_(nullptr),
+    deltaCoeff_
+    (
+        dict.optionalSubDict(type() + "Coeffs").getOrDefault<scalar>
+        (
+            "deltaCoeff",
+            1.035
+        )
+    ),
+    requireUpdate_
+    (
+        dict.optionalSubDict(type() + "Coeffs").getOrDefault<bool>
+        (
+            "requireUpdate",
+            true
+        )
+    ),
+    shielding_
+    (
+        dict.optionalSubDict(type() + "Coeffs").getOrDefault<bool>
+        (
+            "shielding",
+            true
+        )
+    ),
+    FKHmin_
+    (
+        dict.optionalSubDict(type() + "Coeffs").getOrDefault<scalar>
+        (
+            "FKHmin",
+            0.1
+        )
+    ),
+    FKHmax_
+    (
+        dict.optionalSubDict(type() + "Coeffs").getOrDefault<scalar>
+        (
+            "FKHmax",
+            1.0
+        )
+    ),
+    a1_
+    (
+        dict.optionalSubDict(type() + "Coeffs").getOrDefault<scalar>
+        (
+            "a1",
+            0.15
+        )
+    ),
+    a2_
+    (
+        dict.optionalSubDict(type() + "Coeffs").getOrDefault<scalar>
+        (
+            "a2",
+            0.30
+        )
+    ),
+    epsilon_
+    (
+        dict.optionalSubDict(type() + "Coeffs").getOrDefault<scalar>
+        (
+            "epsilon",
+            0.01
+        )
+    ),
+    kappa_
+    (
+        dict.optionalSubDict(type() + "Coeffs").getOrDefault<scalar>
+        (
+            "kappa",
+            0.41
+        )
+    ),
+    Cd1_
+    (
+        dict.optionalSubDict(type() + "Coeffs").getOrDefault<scalar>
+        (
+            "Cd1",
+            20
+        )
+    ),
+    Cd2_
+    (
+        dict.optionalSubDict(type() + "Coeffs").getOrDefault<scalar>
+        (
+            "Cd2",
+            3
+        )
+    )
+{
+    if (dict.optionalSubDict(type() + "Coeffs").found("hmax"))
+    {
+        // User-defined hmax
+        hmaxPtr_ =
+            LESdelta::New
+            (
+                IOobject::groupName("hmax", turbulence.U().group()),
+                turbulence,
+                dict.optionalSubDict("hmaxCoeffs"),
+                "hmax"
+            );
+    }
+    else
+    {
+        Info<< "Employing " << maxDeltaxyz::typeName << " for hmax" << endl;
+
+        hmaxPtr_.reset
+        (
+            new maxDeltaxyz
+            (
+                IOobject::groupName("hmax", turbulence.U().group()),
+                turbulence,
+                dict.optionalSubDict("hmaxCoeffs")
+            )
+        );
+    }
+
+    if (mag(a2_ - a1_) < SMALL)
+    {
+        FatalIOErrorInFunction(dict)
+            << "Model coefficients a1 = " << a1_
+            << ", and a2 = " << a2_ << " cannot be equal."
+            << abort(FatalIOError);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void Foam::LESModels::SLADelta::read(const dictionary& dict)
+{
+    const dictionary& coeffsDict(dict.optionalSubDict(type() + "Coeffs"));
+
+    coeffsDict.readIfPresent<scalar>("deltaCoeff", deltaCoeff_);
+    coeffsDict.readIfPresent<bool>("requireUpdate", requireUpdate_);
+    coeffsDict.readIfPresent<scalar>("FKHmin", FKHmin_);
+    coeffsDict.readIfPresent<scalar>("FKHmax", FKHmax_);
+    coeffsDict.readIfPresent<scalar>("a1", a1_);
+    coeffsDict.readIfPresent<scalar>("a2", a2_);
+    coeffsDict.readIfPresent<scalar>("epsilon", epsilon_);
+    coeffsDict.readIfPresent<scalar>("kappa", kappa_);
+    coeffsDict.readIfPresent<scalar>("Cd1", Cd1_);
+    coeffsDict.readIfPresent<scalar>("Cd2", Cd2_);
+
+    calcDelta();
+}
+
+
+void Foam::LESModels::SLADelta::correct()
+{
+    if (turbulenceModel_.mesh().changing() && requireUpdate_)
+    {
+        hmaxPtr_->correct();
+    }
+
+    calcDelta();
+}
+
+
+// ************************************************************************* //

src/TurbulenceModels/turbulenceModels/LES/LESdeltas/SLADelta/SLADelta.H

@@ -0,0 +1,149 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2022 Upstream CFD GmbH
+    Copyright (C) 2022 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::LESModels::SLADelta
+
+Description
+    Delta formulation that accounts for the orientation of the vorticity vector
+    and a flow-sensitised function to detect 2D flow regions. Delta value is
+    strongly reduced in these regions to accelerate the transition from RANS
+    to LES in hybrid RANS/LES simulations.
+
+    Reference:
+    \verbatim
+        Shur, M. L., Spalart, P. R., Strelets, M. K., & Travin, A. K. (2015).
+        An enhanced version of DES with rapid transition
+        from RANS to LES in separated flows.
+        Flow, turbulence and combustion, 95(4), 709-737.
+        DOI:10.1007/s10494-015-9618-0
+    \endverbatim
+
+SourceFiles
+    SLADelta.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef LESModels_SLADelta_H
+#define LESModels_SLADelta_H
+
+#include "LESdelta.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace LESModels
+{
+
+/*---------------------------------------------------------------------------*\
+                          Class SLADelta Declaration
+\*---------------------------------------------------------------------------*/
+
+class SLADelta
+:
+    public LESdelta
+{
+    // Private Data
+
+        //- Run-time selectable delta for hmax
+        //  Defaults to the maxDeltaxyz model if not supplied
+        autoPtr<LESdelta> hmaxPtr_;
+
+        // Model coefficients
+
+            scalar deltaCoeff_;
+            bool requireUpdate_;
+            bool shielding_;
+            scalar FKHmin_;
+            scalar FKHmax_;
+            scalar a1_;
+            scalar a2_;
+            scalar epsilon_;
+            scalar kappa_;
+            scalar Cd1_;
+            scalar Cd2_;
+
+
+    // Private Member Functions
+
+        // Calculate the delta values
+        void calcDelta();
+
+        //- No copy construct
+        SLADelta(const SLADelta&) = delete;
+
+        //- No copy assignment
+        void operator=(const SLADelta&) = delete;
+
+
+public:
+
+    //- Runtime type information
+    TypeName("SLADelta");
+
+
+    // Constructors
+
+        //- Construct from name, turbulenceModel and IOdictionary
+        SLADelta
+        (
+            const word& name,
+            const turbulenceModel& turbulence,
+            const dictionary&
+        );
+
+
+    //- Destructor
+    virtual ~SLADelta() = default;
+
+
+    // Member Functions
+
+        //- Read the LESdelta dictionary
+        void read(const dictionary&);
+
+        //- Return the hmax delta field
+        const volScalarField& hmax() const
+        {
+            return hmaxPtr_();
+        }
+
+        // Correct values
+        void correct();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace LESModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/TurbulenceModels/turbulenceModels/LES/sigma/sigma.C

@@ -0,0 +1,173 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2022 Upstream CFD GmbH
+    Copyright (C) 2022 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 "sigma.H"
+#include "DESModel.H"
+#include "fvOptions.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace LESModels
+{
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+template<class BasicTurbulenceModel>
+void sigma<BasicTurbulenceModel>::correctNut()
+{
+    this->nut_ =
+        sqr(this->delta())
+       *DESModel<BasicTurbulenceModel>::Ssigma(fvc::grad(this->U_), Csigma_);
+    this->nut_.correctBoundaryConditions();
+    fv::options::New(this->mesh_).correct(this->nut_);
+
+    BasicTurbulenceModel::correctNut();
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class BasicTurbulenceModel> sigma<BasicTurbulenceModel>::sigma
+(
+    const alphaField& alpha,
+    const rhoField& rho,
+    const volVectorField& U,
+    const surfaceScalarField& alphaRhoPhi,
+    const surfaceScalarField& phi,
+    const transportModel& transport,
+    const word& propertiesName,
+    const word& type
+)
+:
+    LESeddyViscosity<BasicTurbulenceModel>
+    (
+        type,
+        alpha,
+        rho,
+        U,
+        alphaRhoPhi,
+        phi,
+        transport,
+        propertiesName
+    ),
+
+    Ck_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "Ck",
+            this->coeffDict_,
+            0.094
+        )
+    ),
+
+    Cw_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "Cw",
+            this->coeffDict_,
+            0.325
+        )
+    ),
+
+    Csigma_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "Csigma",
+            this->coeffDict_,
+            1.68
+        )
+    )
+{
+    if (type == typeName)
+    {
+        this->printCoeffs(type);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class BasicTurbulenceModel>
+bool sigma<BasicTurbulenceModel>::read()
+{
+    if (LESeddyViscosity<BasicTurbulenceModel>::read())
+    {
+        Ck_.readIfPresent(this->coeffDict());
+        Cw_.readIfPresent(this->coeffDict());
+        Csigma_.readIfPresent(this->coeffDict());
+
+        return true;
+    }
+
+    return false;
+}
+
+
+template<class BasicTurbulenceModel>
+tmp<volScalarField> sigma<BasicTurbulenceModel>::k() const
+{
+    return tmp<volScalarField>::New
+    (
+        IOobject::groupName("k", this->U_.group()),
+        (2.0*Ck_/this->Ce_)
+       *sqr(this->delta())
+       *magSqr(dev(symm(fvc::grad(this->U_))))
+    );
+}
+
+
+template<class BasicTurbulenceModel>
+tmp<volScalarField> sigma<BasicTurbulenceModel>::epsilon() const
+{
+    return tmp<volScalarField>::New
+    (
+        IOobject::groupName("epsilon", this->U_.group()),
+        this->Ce_*k()*sqrt(k())/this->delta()
+    );
+}
+
+
+template<class BasicTurbulenceModel>
+void sigma<BasicTurbulenceModel>::correct()
+{
+    LESeddyViscosity<BasicTurbulenceModel>::correct();
+    correctNut();
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace LESModels
+} // End namespace Foam
+
+// ************************************************************************* //

src/TurbulenceModels/turbulenceModels/LES/sigma/sigma.H

@@ -0,0 +1,173 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2022 Upstream CFD GmbH
+    Copyright (C) 2022 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::LESModels::sigma
+
+Group
+    grpLESTurbulence
+
+Description
+    The sigma SGS model.
+
+    Reference:
+    \verbatim
+        Nicoud, F., Toda, H. B., Cabrit, O., Bose, S., & Lee, J. (2011).
+        Using singular values to build a subgrid-scale
+        model for large eddy simulations.
+        Physics of fluids, 23(8), 085106.
+        DOI:10.1063/1.3623274
+    \endverbatim
+
+    The default model coefficients correspond to the following:
+    \verbatim
+        sigmaCoeffs
+        {
+            ck          0.094;
+            Csigma      1.68;
+        }
+    \endverbatim
+
+Note
+    The default value of the Csigma constant implemented was calibrated for
+    OpenFOAM using decaying isotropic turbulence and is slightly higher than
+    the value suggested in the reference publication.
+
+SourceFiles
+    sigma.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef Foam_LESModels_sigma_H
+#define Foam_LESModels_sigma_H
+
+#include "LESModel.H"
+#include "LESeddyViscosity.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace LESModels
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class sigma Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class BasicTurbulenceModel>
+class sigma
+:
+    public LESeddyViscosity<BasicTurbulenceModel>
+{
+    // Private Member Functions
+
+        //- No copy construct
+        sigma(const sigma&) = delete;
+
+        //- No copy assignment
+        void operator=(const sigma&) = delete;
+
+
+protected:
+
+    // Protected Data
+
+        dimensionedScalar Ck_;
+        dimensionedScalar Cw_;
+        dimensionedScalar Csigma_;
+
+
+    // Protected Member Functions
+
+        //- Update the SGS eddy-viscosity
+        virtual void correctNut();
+
+
+public:
+
+    typedef typename BasicTurbulenceModel::alphaField alphaField;
+    typedef typename BasicTurbulenceModel::rhoField rhoField;
+    typedef typename BasicTurbulenceModel::transportModel transportModel;
+
+
+    //- Runtime type information
+    TypeName("sigma");
+
+
+    // Constructors
+
+        //- Construct from components
+        sigma
+        (
+            const alphaField& alpha,
+            const rhoField& rho,
+            const volVectorField& U,
+            const surfaceScalarField& alphaRhoPhi,
+            const surfaceScalarField& phi,
+            const transportModel& transport,
+            const word& propertiesName = turbulenceModel::propertiesName,
+            const word& type = typeName
+        );
+
+
+    //- Destructor
+    virtual ~sigma()
+    {}
+
+
+    // Member Functions
+
+        //- Read model coefficients if they have changed
+        virtual bool read();
+
+        //- Return SGS kinetic energy
+        virtual tmp<volScalarField> k() const;
+
+        //- Return SGS disipation rate
+        virtual tmp<volScalarField> epsilon() const;
+
+        //- Correct Eddy-Viscosity and related properties
+        virtual void correct();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace LESModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "sigma.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/TurbulenceModels/turbulenceModels/Make/files

@@ -10,7 +10,8 @@ $(LESdelta)/smoothDelta/smoothDelta.C
 $(LESdelta)/maxDeltaxyz/maxDeltaxyz.C
 $(LESdelta)/IDDESDelta/IDDESDelta.C
 $(LESdelta)/maxDeltaxyzCubeRootLESDelta/maxDeltaxyzCubeRootLESDelta.C
-
+$(LESdelta)/DeltaOmegaTildeDelta/DeltaOmegaTildeDelta.C
+$(LESdelta)/SLADelta/SLADelta.C
 
 LESfilters = LES/LESfilters
 

src/TurbulenceModels/turbulenceModels/RAS/SpalartAllmaras/SpalartAllmaras.C

@@ -27,9 +27,6 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "SpalartAllmaras.H"
-#include "fvOptions.H"
-#include "bound.H"
-#include "wallDist.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -41,101 +38,14 @@ namespace RASModels
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
 template<class BasicTurbulenceModel>
-tmp<volScalarField> SpalartAllmaras<BasicTurbulenceModel>::chi() const
-{
-    return nuTilda_/this->nu();
-}
-
-
-template<class BasicTurbulenceModel>
-tmp<volScalarField> SpalartAllmaras<BasicTurbulenceModel>::fv1
+Foam::tmp<Foam::volScalarField> SpalartAllmaras<BasicTurbulenceModel>::dTilda
 (
-    const volScalarField& chi
+    const volScalarField& chi,
+    const volScalarField& fv1,
+    const volTensorField& gradU
 ) const
 {
-    const volScalarField chi3(pow3(chi));
-
-    return chi3/(chi3 + pow3(Cv1_));
-}
-
-
-template<class BasicTurbulenceModel>
-tmp<volScalarField::Internal> SpalartAllmaras<BasicTurbulenceModel>::fv2
-(
-    const volScalarField::Internal& chi,
-    const volScalarField::Internal& fv1
-) const
-{
-    return scalar(1) - chi/(scalar(1) + chi*fv1);
-}
-
-
-template<class BasicTurbulenceModel>
-tmp<volScalarField::Internal> SpalartAllmaras<BasicTurbulenceModel>::Stilda()
-const
-{
-    const volScalarField chi(this->chi());
-
-    const volScalarField fv1(this->fv1(chi));
-
-    const volScalarField::Internal Omega
-    (
-        ::sqrt(scalar(2))*mag(skew(fvc::grad(this->U_)().v()))
-    );
-
-    return
-    (
-        max
-        (
-            Omega + fv2(chi(), fv1())*nuTilda_()/sqr(kappa_*y_),
-            Cs_*Omega
-        )
-    );
-}
-
-
-template<class BasicTurbulenceModel>
-tmp<volScalarField::Internal> SpalartAllmaras<BasicTurbulenceModel>::fw
-(
-    const volScalarField::Internal& Stilda
-) const
-{
-    const volScalarField::Internal r
-    (
-        min
-        (
-            nuTilda_
-           /(
-               max
-               (
-                   Stilda,
-                   dimensionedScalar(Stilda.dimensions(), SMALL)
-               )
-              *sqr(kappa_*y_)
-            ),
-            scalar(10)
-        )
-    );
-
-    const volScalarField::Internal g(r + Cw2_*(pow6(r) - r));
-
-    return
-        g*pow
-        (
-            (scalar(1) + pow6(Cw3_))/(pow6(g) + pow6(Cw3_)),
-            scalar(1)/scalar(6)
-        );
-}
-
-
-template<class BasicTurbulenceModel>
-void SpalartAllmaras<BasicTurbulenceModel>::correctNut()
-{
-    this->nut_ = nuTilda_*this->fv1(this->chi());
-    this->nut_.correctBoundaryConditions();
-    fv::options::New(this->mesh_).correct(this->nut_);
-
-    BasicTurbulenceModel::correctNut();
+    return this->y_;
 }
 
 
@@ -154,7 +64,7 @@ SpalartAllmaras<BasicTurbulenceModel>::SpalartAllmaras
     const word& type
 )
 :
-    eddyViscosity<RASModel<BasicTurbulenceModel>>
+    SpalartAllmarasBase<eddyViscosity<RASModel<BasicTurbulenceModel>>>
     (
         type,
         alpha,
@@ -164,97 +74,7 @@ SpalartAllmaras<BasicTurbulenceModel>::SpalartAllmaras
         phi,
         transport,
         propertiesName
-    ),
-
-    sigmaNut_
-    (
-        dimensioned<scalar>::getOrAddToDict
-        (
-            "sigmaNut",
-            this->coeffDict_,
-            scalar(2)/scalar(3)
-        )
-    ),
-    kappa_
-    (
-        dimensioned<scalar>::getOrAddToDict
-        (
-            "kappa",
-            this->coeffDict_,
-            0.41
-        )
-    ),
-
-    Cb1_
-    (
-        dimensioned<scalar>::getOrAddToDict
-        (
-            "Cb1",
-            this->coeffDict_,
-            0.1355
-        )
-    ),
-    Cb2_
-    (
-        dimensioned<scalar>::getOrAddToDict
-        (
-            "Cb2",
-            this->coeffDict_,
-            0.622
-        )
-    ),
-    Cw1_(Cb1_/sqr(kappa_) + (scalar(1) + Cb2_)/sigmaNut_),
-    Cw2_
-    (
-        dimensioned<scalar>::getOrAddToDict
-        (
-            "Cw2",
-            this->coeffDict_,
-            0.3
-        )
-    ),
-    Cw3_
-    (
-        dimensioned<scalar>::getOrAddToDict
-        (
-            "Cw3",
-            this->coeffDict_,
-            2.0
-        )
-    ),
-    Cv1_
-    (
-        dimensioned<scalar>::getOrAddToDict
-        (
-            "Cv1",
-            this->coeffDict_,
-            7.1
-        )
-    ),
-    Cs_
-    (
-        dimensioned<scalar>::getOrAddToDict
-        (
-            "Cs",
-            this->coeffDict_,
-            0.3
-        )
-    ),
-
-    nuTilda_
-    (
-        IOobject
-        (
-            "nuTilda",
-            this->runTime_.timeName(),
-            this->mesh_,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        this->mesh_
-    ),
-
-    y_(wallDist::New(this->mesh_).y())
+    )
 {
     if (type == typeName)
     {
@@ -263,153 +83,6 @@ SpalartAllmaras<BasicTurbulenceModel>::SpalartAllmaras
 }
 
 
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-template<class BasicTurbulenceModel>
-bool SpalartAllmaras<BasicTurbulenceModel>::read()
-{
-    if (eddyViscosity<RASModel<BasicTurbulenceModel>>::read())
-    {
-        sigmaNut_.readIfPresent(this->coeffDict());
-        kappa_.readIfPresent(this->coeffDict());
-
-        Cb1_.readIfPresent(this->coeffDict());
-        Cb2_.readIfPresent(this->coeffDict());
-        Cw1_ = Cb1_/sqr(kappa_) + (scalar(1) + Cb2_)/sigmaNut_;
-        Cw2_.readIfPresent(this->coeffDict());
-        Cw3_.readIfPresent(this->coeffDict());
-        Cv1_.readIfPresent(this->coeffDict());
-        Cs_.readIfPresent(this->coeffDict());
-
-        return true;
-    }
-
-    return false;
-}
-
-
-template<class BasicTurbulenceModel>
-tmp<volScalarField> SpalartAllmaras<BasicTurbulenceModel>::DnuTildaEff() const
-{
-    return tmp<volScalarField>::New
-    (
-        "DnuTildaEff",
-        (nuTilda_ + this->nu())/sigmaNut_
-    );
-}
-
-
-template<class BasicTurbulenceModel>
-tmp<volScalarField> SpalartAllmaras<BasicTurbulenceModel>::k() const
-{
-    // (B:Eq. 4.50)
-    const scalar Cmu = 0.09;
-
-    return tmp<volScalarField>::New
-    (
-        IOobject
-        (
-            IOobject::groupName("k", this->alphaRhoPhi_.group()),
-            this->runTime_.timeName(),
-            this->mesh_
-        ),
-        cbrt(this->fv1(this->chi()))
-        *nuTilda_
-        *::sqrt(scalar(2)/Cmu)
-        *mag(symm(fvc::grad(this->U_))),
-        this->nut_.boundaryField().types()
-    );
-}
-
-
-template<class BasicTurbulenceModel>
-tmp<volScalarField> SpalartAllmaras<BasicTurbulenceModel>::epsilon() const
-{
-    // (B:Eq. 4.50)
-    const scalar Cmu = 0.09;
-    const dimensionedScalar nutSMALL(sqr(dimLength)/dimTime, SMALL);
-
-    return tmp<volScalarField>::New
-    (
-        IOobject
-        (
-            IOobject::groupName("epsilon", this->alphaRhoPhi_.group()),
-            this->runTime_.timeName(),
-            this->mesh_
-        ),
-        pow(this->fv1(this->chi()), 0.5)
-        *pow(::sqrt(Cmu)*this->k(), 2)
-        /(nuTilda_ + this->nut_ + nutSMALL),
-        this->nut_.boundaryField().types()
-    );
-}
-
-
-template<class BasicTurbulenceModel>
-tmp<volScalarField> SpalartAllmaras<BasicTurbulenceModel>::omega() const
-{
-    // (P:p. 384)
-    const scalar betaStar = 0.09;
-    const dimensionedScalar k0(sqr(dimLength/dimTime), SMALL);
-
-    return tmp<volScalarField>::New
-    (
-        IOobject
-        (
-            IOobject::groupName("omega", this->alphaRhoPhi_.group()),
-            this->runTime_.timeName(),
-            this->mesh_
-        ),
-        this->epsilon()/(betaStar*(this->k() + k0)),
-        this->nut_.boundaryField().types()
-    );
-}
-
-
-template<class BasicTurbulenceModel>
-void SpalartAllmaras<BasicTurbulenceModel>::correct()
-{
-    if (!this->turbulence_)
-    {
-        return;
-    }
-
-    {
-        // Construct local convenience references
-        const alphaField& alpha = this->alpha_;
-        const rhoField& rho = this->rho_;
-        const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_;
-        fv::options& fvOptions(fv::options::New(this->mesh_));
-
-        eddyViscosity<RASModel<BasicTurbulenceModel>>::correct();
-
-        const volScalarField::Internal Stilda(this->Stilda());
-
-        tmp<fvScalarMatrix> nuTildaEqn
-        (
-            fvm::ddt(alpha, rho, nuTilda_)
-          + fvm::div(alphaRhoPhi, nuTilda_)
-          - fvm::laplacian(alpha*rho*DnuTildaEff(), nuTilda_)
-          - Cb2_/sigmaNut_*alpha*rho*magSqr(fvc::grad(nuTilda_))
-         ==
-            Cb1_*alpha()*rho()*Stilda*nuTilda_()
-          - fvm::Sp(Cw1_*alpha()*rho()*fw(Stilda)*nuTilda_()/sqr(y_), nuTilda_)
-          + fvOptions(alpha, rho, nuTilda_)
-        );
-
-        nuTildaEqn.ref().relax();
-        fvOptions.constrain(nuTildaEqn.ref());
-        solve(nuTildaEqn);
-        fvOptions.correct(nuTilda_);
-        bound(nuTilda_, dimensionedScalar(nuTilda_.dimensions(), Zero));
-        nuTilda_.correctBoundaryConditions();
-    }
-
-    // Update nut with latest available nuTilda
-    correctNut();
-}
-
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace RASModels

src/TurbulenceModels/turbulenceModels/RAS/SpalartAllmaras/SpalartAllmaras.H

@@ -52,18 +52,6 @@ Description
             Spalart-Allmaras One-Equation Model without ft2 Term (SA-noft2).
             https://turbmodels.larc.nasa.gov/spalart.html#sanoft2
             (Retrieved:12-01-2021).
-
-        Estimation expression for k and epsilon (tag:B), Eq. 4.50:
-            Bourgoin, A. (2019).
-            Bathymetry induced turbulence modelling the
-            Alderney Race site: regional approach with TELEMAC-LES.
-            Normandie Université.
-
-        Estimation expressions for omega (tag:P):
-            Pope, S. B. (2000).
-            Turbulent flows.
-            Cambridge, UK: Cambridge Univ. Press
-            DOI:10.1017/CBO9780511840531
     \endverbatim
 
 Usage
@@ -96,8 +84,6 @@ Note
   - The model is implemented without the trip-term since the model has almost
     always been used in fully turbulent applications rather than those where
     laminar-turbulent transition occurs.
-  - It has been argued that the \c ft2 term is not needed in the absence of the
-    trip-term, hence \c ft2 term is also not implementated.
   - The \c Stilda generation term should never be allowed to be zero or negative
     to avoid potential numerical issues and unphysical results for complex
     flows. To this end, a limiter proposed by Spalart (R:Note-1(b)) is applied
@@ -112,11 +98,12 @@ SourceFiles
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef SpalartAllmaras_H
-#define SpalartAllmaras_H
+#ifndef Foam_SpalartAllmaras_H
+#define Foam_SpalartAllmaras_H
 
 #include "RASModel.H"
 #include "eddyViscosity.H"
+#include "SpalartAllmarasBase.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -132,7 +119,7 @@ namespace RASModels
 template<class BasicTurbulenceModel>
 class SpalartAllmaras
 :
-    public eddyViscosity<RASModel<BasicTurbulenceModel>>
+    public SpalartAllmarasBase<eddyViscosity<RASModel<BasicTurbulenceModel>>>
 {
     // Private Member Functions
 
@@ -145,55 +132,16 @@ class SpalartAllmaras
 
 protected:
 
-    // Protected Data
-
-        // Model coefficients
-
-            dimensionedScalar sigmaNut_;
-            dimensionedScalar kappa_;
-            dimensionedScalar Cb1_;
-            dimensionedScalar Cb2_;
-            dimensionedScalar Cw1_;
-            dimensionedScalar Cw2_;
-            dimensionedScalar Cw3_;
-            dimensionedScalar Cv1_;
-            dimensionedScalar Cs_;
-
-
-        // Fields
-
-            //- Modified kinematic viscosity [m2/s]
-            volScalarField nuTilda_;
-
-            //- Wall distance
-            //  Note: different to wall distance in parent RASModel
-            //  which is for near-wall cells only
-            const volScalarField::Internal& y_;
-
-
     // Protected Member Functions
 
-        tmp<volScalarField> chi() const;
-
-        tmp<volScalarField> fv1(const volScalarField& chi) const;
-
-        tmp<volScalarField::Internal> fv2
+        //- Return the length scale
+        virtual tmp<volScalarField> dTilda
         (
-            const volScalarField::Internal& chi,
-            const volScalarField::Internal& fv1
+            const volScalarField& chi,
+            const volScalarField& fv1,
+            const volTensorField& gradU
         ) const;
 
-        tmp<volScalarField::Internal> Stilda() const;
-
-        tmp<volScalarField::Internal> fw
-        (
-            const volScalarField::Internal& Stilda
-        ) const;
-
-
-        //- Update nut with the latest available nuTilda
-        virtual void correctNut();
-
 
 public:
 
@@ -224,27 +172,6 @@ public:
 
     //- Destructor
     virtual ~SpalartAllmaras() = default;
-
-
-    // Member Functions
-
-        //- Re-read model coefficients if they have changed
-        virtual bool read();
-
-        //- Return the effective diffusivity for nuTilda
-        tmp<volScalarField> DnuTildaEff() const;
-
-        //- Return the (estimated) turbulent kinetic energy
-        virtual tmp<volScalarField> k() const;
-
-        //- Return the (estimated) turbulent kinetic energy dissipation rate
-        virtual tmp<volScalarField> epsilon() const;
-
-        //- Return the (estimated) specific dissipation rate
-        virtual tmp<volScalarField> omega() const;
-
-        //- Solve the turbulence equations and correct the turbulent viscosity
-        virtual void correct();
 };
 
 

src/dynamicMesh/motionSolvers/displacement/points0/points0MotionSolver.C

@@ -6,7 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2012-2017 OpenFOAM Foundation
-    Copyright (C) 2015-2021 OpenCFD Ltd.
+    Copyright (C) 2015-2022 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -84,7 +84,7 @@ Foam::points0MotionSolver::points0MotionSolver
 )
 :
     motionSolver(mesh, dict, type),
-    zoneMotion(dict, mesh),
+    zoneMotion(coeffDict(), mesh),
     points0_(points0IO(mesh))
 {
     if
@@ -128,7 +128,7 @@ Foam::points0MotionSolver::points0MotionSolver
 )
 :
     motionSolver(mesh, dict, type),
-    zoneMotion(dict, mesh),
+    zoneMotion(coeffDict(), mesh),
     points0_(points0)
 {
     if (points0_.size() != mesh.nPoints())

src/engine/include/StCorr.H

@@ -80,7 +80,7 @@
             );
             dimensionedScalar AkEst = gSum(mgb*mesh.V().field());
 
-            StCorr.value() = Foam::max(Foam::min((Ak/AkEst).value(), 10.0), 1.0);
+            StCorr.value() = max(min((Ak/AkEst).value(), 10.0), 1.0);
 
             Info<< "StCorr = " << StCorr.value() << endl;
         }