adding run-time selectable thermal wall function capability to compressible RAS models

src/turbulenceModels/RAS/compressible/LRR/LRR.C

@@ -204,11 +204,20 @@ LRR::LRR
             IOobject::AUTO_WRITE
         ),
         autoCreateMut("mut", mesh_)
+    ),
+    alphat_
+    (
+        IOobject
+        (
+            "alphat",
+            runTime_.timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        autoCreateAlphat("alphat", mesh_)
     )
 {
-    mut_ == Cmu_*rho_*sqr(k_)/(epsilon_ + epsilonSmall_);
-    mut_.correctBoundaryConditions();
-
     if (couplingFactor_.value() < 0.0 || couplingFactor_.value() > 1.0)
     {
         FatalErrorIn
@@ -221,6 +230,12 @@ LRR::LRR
             << exit(FatalError);
     }
 
+    mut_ == Cmu_*rho_*sqr(k_)/(epsilon_ + epsilonSmall_);
+    mut_.correctBoundaryConditions();
+
+    alphat_ == mut_/Prt_;
+    alphat_.correctBoundaryConditions();
+
     printCoeffs();
 }
 
@@ -312,6 +327,11 @@ void LRR::correct()
         // Re-calculate viscosity
         mut_ == rho_*Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_);
         mut_.correctBoundaryConditions();
+
+        // Re-calculate thermal diffusivity
+        alphat_ = mut_/Prt_;
+        alphat_.correctBoundaryConditions();
+
         return;
     }
 
@@ -402,6 +422,10 @@ void LRR::correct()
     mut_ == rho_*Cmu_*sqr(k_)/epsilon_;
     mut_.correctBoundaryConditions();
 
+    // Re-calculate thermal diffusivity
+    alphat_ = mut_/Prt_;
+    alphat_.correctBoundaryConditions();
+
 
     // Correct wall shear stresses
 

src/turbulenceModels/RAS/compressible/LRR/LRR.H

@@ -97,6 +97,7 @@ class LRR
         volScalarField k_;
         volScalarField epsilon_;
         volScalarField mut_;
+        volScalarField alphat_;
 
 
 public:
@@ -152,7 +153,7 @@ public:
         {
             return tmp<volScalarField>
             (
-                new volScalarField("alphaEff", alphah_*mut_ + alpha())
+                new volScalarField("alphaEff", alphah_*alphat_ + alpha())
             );
         }
 

src/turbulenceModels/RAS/compressible/LaunderGibsonRSTM/LaunderGibsonRSTM.C

@@ -226,11 +226,20 @@ LaunderGibsonRSTM::LaunderGibsonRSTM
             IOobject::AUTO_WRITE
         ),
         autoCreateMut("mut", mesh_)
+    ),
+    alphat_
+    (
+        IOobject
+        (
+            "alphat",
+            runTime_.timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        autoCreateAlphat("alphat", mesh_)
     )
 {
-    mut_ == Cmu_*rho_*sqr(k_)/(epsilon_ + epsilonSmall_);
-    mut_.correctBoundaryConditions();
-
     if (couplingFactor_.value() < 0.0 || couplingFactor_.value() > 1.0)
     {
         FatalErrorIn
@@ -243,6 +252,12 @@ LaunderGibsonRSTM::LaunderGibsonRSTM
             << exit(FatalError);
     }
 
+    mut_ == Cmu_*rho_*sqr(k_)/(epsilon_ + epsilonSmall_);
+    mut_.correctBoundaryConditions();
+
+    alphat_ == mut_/Prt_;
+    alphat_.correctBoundaryConditions();
+
     printCoeffs();
 }
 
@@ -337,6 +352,11 @@ void LaunderGibsonRSTM::correct()
         // Re-calculate viscosity
         mut_ == rho_*Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_);
         mut_.correctBoundaryConditions();
+
+        // Re-calculate thermal diffusivity
+        alphat_ = mut_/Prt_;
+        alphat_.correctBoundaryConditions();
+
         return;
     }
 
@@ -441,6 +461,9 @@ void LaunderGibsonRSTM::correct()
     mut_ == Cmu_*rho_*sqr(k_)/epsilon_;
     mut_.correctBoundaryConditions();
 
+    // Re-calculate thermal diffusivity
+    alphat_ = mut_/Prt_;
+    alphat_.correctBoundaryConditions();
 
     // Correct wall shear stresses
 

src/turbulenceModels/RAS/compressible/LaunderGibsonRSTM/LaunderGibsonRSTM.H

@@ -104,6 +104,7 @@ class LaunderGibsonRSTM
         volScalarField k_;
         volScalarField epsilon_;
         volScalarField mut_;
+        volScalarField alphat_;
 
 
 public:
@@ -161,7 +162,7 @@ public:
         {
             return tmp<volScalarField>
             (
-                new volScalarField("alphaEff", alphah_*mut_ + alpha())
+                new volScalarField("alphaEff", alphah_*alphat_ + alpha())
             );
         }
 

src/turbulenceModels/RAS/compressible/Make/files

@@ -14,6 +14,9 @@ kOmegaSST/kOmegaSST.C
 /* Wall functions */
 wallFunctions = derivedFvPatchFields/wallFunctions
 
+alphatWallFunctions = $(wallFunctions)/alphatWallFunctions
+$(alphatWallFunctions)/alphatWallFunction/alphatWallFunctionFvPatchScalarField.C
+
 mutWallFunctions = $(wallFunctions)/mutWallFunctions
 $(mutWallFunctions)/mutWallFunction/mutWallFunctionFvPatchScalarField.C
 $(mutWallFunctions)/mutRoughWallFunction/mutRoughWallFunctionFvPatchScalarField.C

src/turbulenceModels/RAS/compressible/RASModel/RASModel.C

@@ -47,7 +47,8 @@ void RASModel::printCoeffs()
 {
     if (printCoeffs_)
     {
-        Info<< type() << "Coeffs" << coeffDict_ << endl;
+        Info<< type() << "Coeffs" << coeffDict_ << nl
+            << "wallFunctionCoeffs" << wallFunctionDict_ << endl;
     }
 }
 
@@ -115,6 +116,15 @@ RASModel::RASModel
             0.09
         )
     ),
+    Prt_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "Prt",
+            wallFunctionDict_,
+            0.85
+        )
+    ),
 
     yPlusLam_(yPlusLam(kappa_.value(), E_.value())),
 
@@ -148,11 +158,9 @@ tmp<scalarField> RASModel::yPlus(const label patchNo) const
     tmp<scalarField> tYp(new scalarField(curPatch.size()));
     scalarField& Yp = tYp();
 
-    if (typeid(curPatch) == typeid(wallFvPatch))
+    if (isType<wallFvPatch>(curPatch))
     {
-        scalar Cmu(readScalar(coeffDict_.lookup("Cmu")));
-
-        Yp = pow(Cmu, 0.25)
+        Yp = pow(Cmu_.value(), 0.25)
             *y_[patchNo]
             *sqrt(k()().boundaryField()[patchNo].patchInternalField())
            /(
@@ -165,8 +173,8 @@ tmp<scalarField> RASModel::yPlus(const label patchNo) const
         WarningIn
         (
             "tmp<scalarField> RASModel::yPlus(const label patchNo) const"
-        )   << "Patch " << patchNo << " is not a wall.  Returning blank field"
-            << endl;
+        )   << "Patch " << patchNo << " is not a wall. Returning null field"
+            << nl << endl;
 
         Yp.setSize(0);
     }
@@ -191,8 +199,11 @@ bool RASModel::read()
         lookup("turbulence") >> turbulence_;
         coeffDict_ = subDict(type() + "Coeffs");
 
-        kappa_.readIfPresent(subDict("wallFunctionCoeffs"));
-        E_.readIfPresent(subDict("wallFunctionCoeffs"));
+        wallFunctionDict_ = subDict("wallFunctionCoeffs");
+        kappa_.readIfPresent(wallFunctionDict_);
+        E_.readIfPresent(wallFunctionDict_);
+        Cmu_.readIfPresent(wallFunctionDict_);
+        Prt_.readIfPresent(wallFunctionDict_);
 
         yPlusLam_ = yPlusLam(kappa_.value(), E_.value());
 

src/turbulenceModels/RAS/compressible/RASModel/RASModel.H

@@ -95,6 +95,7 @@ protected:
         dimensionedScalar kappa_;
         dimensionedScalar E_;
         dimensionedScalar Cmu_;
+        dimensionedScalar Prt_;
 
         scalar yPlusLam_;
 
@@ -244,6 +245,12 @@ public:
                 return Cmu_;
             }
 
+            //- Return turbulent Prandtl number for use in wall-functions
+            dimensionedScalar Prt() const
+            {
+                return Prt_;
+            }
+
             //- Return the near wall distances
             const nearWallDist& y() const
             {

src/turbulenceModels/RAS/compressible/RNGkEpsilon/RNGkEpsilon.C

@@ -174,11 +174,26 @@ RNGkEpsilon::RNGkEpsilon
             IOobject::AUTO_WRITE
         ),
         autoCreateMut("mut", mesh_)
+    ),
+    alphat_
+    (
+        IOobject
+        (
+            "alphat",
+            runTime_.timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        autoCreateAlphat("alphat", mesh_)
     )
 {
     mut_ == Cmu_*rho_*sqr(k_)/(epsilon_ + epsilonSmall_);
     mut_.correctBoundaryConditions();
 
+    alphat_ == mut_/Prt_;
+    alphat_.correctBoundaryConditions();
+
     printCoeffs();
 }
 
@@ -265,6 +280,11 @@ void RNGkEpsilon::correct()
         // Re-calculate viscosity
         mut_ == rho_*Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_);
         mut_.correctBoundaryConditions();
+
+        // Re-calculate thermal diffusivity
+        alphat_ = mut_/Prt_;
+        alphat_.correctBoundaryConditions();
+
         return;
     }
 
@@ -332,6 +352,10 @@ void RNGkEpsilon::correct()
     // Re-calculate viscosity
     mut_ == rho_*Cmu_*sqr(k_)/epsilon_;
     mut_.correctBoundaryConditions();
+
+    // Re-calculate thermal diffusivity
+    alphat_ = mut_/Prt_;
+    alphat_.correctBoundaryConditions();
 }
 
 

src/turbulenceModels/RAS/compressible/RNGkEpsilon/RNGkEpsilon.H

@@ -87,6 +87,7 @@ class RNGkEpsilon
         volScalarField k_;
         volScalarField epsilon_;
         volScalarField mut_;
+        volScalarField alphat_;
 
 
 public:
@@ -142,7 +143,7 @@ public:
         {
             return tmp<volScalarField>
             (
-                new volScalarField("alphaEff", alphah_*mut_ + alpha())
+                new volScalarField("alphaEff", alphah_*alphat_ + alpha())
             );
         }
 

src/turbulenceModels/RAS/compressible/backwardsCompatibilityWallFunctions/backwardsCompatibilityWallFunctions.C

@@ -27,6 +27,7 @@ License
 #include "backwardsCompatibilityWallFunctions.H"
 
 #include "calculatedFvPatchField.H"
+#include "alphatWallFunctionFvPatchScalarField.H"
 #include "mutWallFunctionFvPatchScalarField.H"
 #include "epsilonWallFunctionFvPatchScalarField.H"
 #include "kQRWallFunctionFvPatchField.H"
@@ -41,6 +42,76 @@ namespace compressible
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+tmp<volScalarField> autoCreateAlphat
+(
+    const word& fieldName,
+    const fvMesh& mesh
+)
+{
+    IOobject alphatHeader
+    (
+        fieldName,
+        mesh.time().timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE,
+        false
+    );
+
+    if (alphatHeader.headerOk())
+    {
+        return tmp<volScalarField>(new volScalarField(alphatHeader, mesh));
+    }
+    else
+    {
+        Info<< "--> Upgrading " << fieldName << " to employ run-time "
+            << "selectable wall functions" << endl;
+
+        const fvBoundaryMesh& bm = mesh.boundary();
+
+        wordList alphatBoundaryTypes(bm.size());
+
+        forAll(bm, patchI)
+        {
+            if (isType<wallFvPatch>(bm[patchI]))
+            {
+                alphatBoundaryTypes[patchI] =
+                    RASModels::alphatWallFunctionFvPatchScalarField::typeName;
+            }
+            else
+            {
+                alphatBoundaryTypes[patchI] =
+                    calculatedFvPatchField<scalar>::typeName;
+            }
+        }
+
+        tmp<volScalarField> alphat
+        (
+            new volScalarField
+            (
+                IOobject
+                (
+                    fieldName,
+                    mesh.time().timeName(),
+                    mesh,
+                    IOobject::NO_READ,
+                    IOobject::NO_WRITE,
+                    false
+                ),
+                mesh,
+                dimensionedScalar("zero", dimDensity*dimArea/dimTime, 0.0),
+                alphatBoundaryTypes
+            )
+        );
+
+        Info<< "    Writing updated " << fieldName << endl;
+        alphat().write();
+
+        return alphat;
+    }
+}
+
+
 tmp<volScalarField> autoCreateMut
 (
     const word& fieldName,

src/turbulenceModels/RAS/compressible/backwardsCompatibilityWallFunctions/backwardsCompatibilityWallFunctions.H

@@ -53,6 +53,13 @@ namespace compressible
         const fvMesh& mesh
     );
 
+    //- alphat
+    tmp<volScalarField> autoCreateAlphat
+    (
+        const word& fieldName,
+        const fvMesh& mesh
+    );
+
     //- epsilon
     tmp<volScalarField> autoCreateEpsilon
     (

src/turbulenceModels/RAS/compressible/derivedFvPatchFields/wallFunctions/alphatWallFunctions/alphatWallFunction/alphatWallFunctionFvPatchScalarField.C

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "alphatWallFunctionFvPatchScalarField.H"
+#include "RASModel.H"
+#include "fvPatchFieldMapper.H"
+#include "volFields.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace compressible
+{
+namespace RASModels
+{
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+alphatWallFunctionFvPatchScalarField::
+alphatWallFunctionFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedValueFvPatchScalarField(p, iF)
+{}
+
+
+alphatWallFunctionFvPatchScalarField::
+alphatWallFunctionFvPatchScalarField
+(
+    const alphatWallFunctionFvPatchScalarField& ptf,
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const fvPatchFieldMapper& mapper
+)
+:
+    fixedValueFvPatchScalarField(ptf, p, iF, mapper)
+{}
+
+
+alphatWallFunctionFvPatchScalarField::
+alphatWallFunctionFvPatchScalarField
+(
+    const fvPatch& p,
+    const DimensionedField<scalar, volMesh>& iF,
+    const dictionary& dict
+)
+:
+    fixedValueFvPatchScalarField(p, iF, dict)
+{}
+
+
+alphatWallFunctionFvPatchScalarField::
+alphatWallFunctionFvPatchScalarField
+(
+    const alphatWallFunctionFvPatchScalarField& awfpsf
+)
+:
+    fixedValueFvPatchScalarField(awfpsf)
+{}
+
+
+alphatWallFunctionFvPatchScalarField::
+alphatWallFunctionFvPatchScalarField
+(
+    const alphatWallFunctionFvPatchScalarField& awfpsf,
+    const DimensionedField<scalar, volMesh>& iF
+)
+:
+    fixedValueFvPatchScalarField(awfpsf, iF)
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void alphatWallFunctionFvPatchScalarField::updateCoeffs()
+{
+    const RASModel& ras = db().lookupObject<RASModel>("RASProperties");
+    const scalar Prt = ras.Prt().value();
+
+    const scalarField& mutw =
+        patch().lookupPatchField<volScalarField, scalar>("mut");
+
+    operator==(mutw/Prt);
+}
+
+
+void alphatWallFunctionFvPatchScalarField::write(Ostream& os) const
+{
+    fvPatchField<scalar>::write(os);
+    writeEntry("value", os);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+makePatchTypeField(fvPatchScalarField, alphatWallFunctionFvPatchScalarField);
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace RASModels
+} // End namespace compressible
+} // End namespace Foam
+
+// ************************************************************************* //

src/turbulenceModels/RAS/compressible/derivedFvPatchFields/wallFunctions/alphatWallFunctions/alphatWallFunction/alphatWallFunctionFvPatchScalarField.H

+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 1991-2008 OpenCFD Ltd.
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 2 of the License, or (at your
+    option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Class
+    Foam::compressible::RASModels::alphatWallFunctionFvPatchScalarField
+
+Description
+    Boundary condition for turbulent thermal diffusivity when using wall
+    functions
+    - replicates OpenFOAM v1.5 (and earlier) behaviour
+