Merge branch 'feature-one-equation-turbulence-model' into 'develop'

ENH: new RANS model for geophysical applications

See merge request Development/openfoam!493

src/atmosphericModels/Make/files

@@ -1,5 +1,6 @@
 /* Models */
-atmosphericTurbulentTransportModels.C
+turbulenceModels/atmosphericTurbulentTransportModels.C
+turbulenceModels/atmosphericTurbulentFluidThermoModels.C
 porosityModels/powerLawLopesdaCosta/powerLawLopesdaCosta.C
 
 

src/atmosphericModels/turbulenceModels/RAS/kL/kL.C

@@ -0,0 +1,494 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "kL.H"
+#include "fvOptions.H"
+#include "bound.H"
+#include "gravityMeshObject.H"
+#include "wallDist.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace RASModels
+{
+
+// * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
+
+template<class BasicTurbulenceModel>
+tmp<volScalarField> kL<BasicTurbulenceModel>::Cmu() const
+{
+    // (A:Eq. 31)
+    return (0.556 + 0.108*Rt_)/(1.0 + 0.308*Rt_ + 0.00837*sqr(Rt_));
+}
+
+
+template<class BasicTurbulenceModel>
+tmp<volScalarField> kL<BasicTurbulenceModel>::CmuPrime() const
+{
+    // (A:Eq. 32)
+    return 0.556/(1.0 + 0.277*Rt_);
+}
+
+
+template<class BasicTurbulenceModel>
+tmp<volScalarField> kL<BasicTurbulenceModel>::nutPrime() const
+{
+    // (A:Eq. 12)
+    return CmuPrime()*sqrt(k_)*L_;
+}
+
+
+template<class BasicTurbulenceModel>
+tmp<volScalarField> kL<BasicTurbulenceModel>::epsilonCanopy() const
+{
+    const auto* CdPtr =
+        this->mesh_.template findObject<volScalarField>("plantCd");
+    const auto* LADPtr =
+        this->mesh_.template findObject<volScalarField>("leafAreaDensity");
+    const volVectorField& U = this->U_;
+
+    if (CdPtr && LADPtr)
+    {
+        const auto& Cd = *CdPtr;
+        const auto& LAD = *LADPtr;
+
+        // (W:Eq. 13)
+        return Cd*LAD*mag(U)*k_;
+    }
+
+    return tmp<volScalarField>::New
+    (
+        IOobject
+        (
+            IOobject::groupName("epsilonCanopy", this->alphaRhoPhi_.group()),
+            this->runTime_.timeName(),
+            this->mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false
+        ),
+        this->mesh_,
+        dimensionedScalar(sqr(dimLength)/pow3(dimTime), Zero)
+    );
+}
+
+
+template<class BasicTurbulenceModel>
+tmp<volScalarField> kL<BasicTurbulenceModel>::epsilon() const
+{
+    // (W:Eq. 13)
+    tmp<volScalarField> tepsilonCanopy = epsilonCanopy();
+
+    // (A:Eq. 19)
+    tmp<volScalarField> tepsilonPlain = pow3(Cmu0_)*pow(k_, 1.5)/L_;
+
+    // (W:Eq. 13)
+    tmp<volScalarField> tepsilon = max(tepsilonPlain, tepsilonCanopy);
+    volScalarField& epsilon = tepsilon.ref();
+    bound(epsilon, this->epsilonMin_);
+
+    return tepsilon;
+}
+
+
+template<class BasicTurbulenceModel>
+tmp<volScalarField> kL<BasicTurbulenceModel>::canopyHeight() const
+{
+    const auto* canopyHeightPtr =
+        this->mesh_.template findObject<volScalarField>("canopyHeight");
+
+    if (canopyHeightPtr)
+    {
+        const auto& canopyHeight = *canopyHeightPtr;
+        return canopyHeight;
+    }
+
+    return tmp<volScalarField>::New
+    (
+        IOobject
+        (
+            IOobject::groupName("canopyHeight", this->alphaRhoPhi_.group()),
+            this->runTime_.timeName(),
+            this->mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE,
+            false
+        ),
+        this->mesh_,
+        dimensionedScalar(dimLength, Zero)
+    );
+}
+
+
+template<class BasicTurbulenceModel>
+tmp<volScalarField> kL<BasicTurbulenceModel>::L() const
+{
+    // (A:Eq. 22)
+    const volScalarField Lplain(kappa_*y_);
+
+    // Increase roughness for canopy (forest, vegetation etc)
+    tmp<volScalarField> tLcanopy = kappa_*canopyHeight();
+    const volScalarField& Lcanopy = tLcanopy;
+
+    // (W:Eq. 16)
+    return max(Lcanopy, Lplain);
+}
+
+
+template<class BasicTurbulenceModel>
+void kL<BasicTurbulenceModel>::stratification(const volScalarField& fVB)
+{
+    tmp<volScalarField> tLg = L();
+    const volScalarField& Lg = tLg.cref();
+
+    tmp<volScalarField> tcanopyHeight = canopyHeight();
+    const volScalarField& canopyHeight = tcanopyHeight;
+
+    tmp<volScalarField> tLcanopy = kappa_*canopyHeight;
+    const volScalarField& Lcanopy = tLcanopy;
+
+    const scalar Cmu0 = Cmu0_.value();
+    const scalar CbStable = CbStable_.value();
+    const scalar CbUnstable = CbUnstable_.value();
+
+    forAll(L_, i)
+    {
+        if (y_[i] > canopyHeight[i])
+        {
+            if (fVB[i] > 0)
+            {
+                // (A:Eq. 23)
+                const scalar Lb = CbStable*sqrt(k_[i])/sqrt(fVB[i]);
+
+                // (A:Eq. 26)
+                L_[i] = sqrt(sqr(Lg[i]*Lb)/(sqr(Lg[i]) + sqr(Lb)));
+            }
+            else
+            {
+                // For unstable/neutral boundary layer (A:p. 80)
+                // Smoothing function for turbulent Richardson
+                // number to ensure gentle transition into
+                // the regime of strong convection
+                Rt_[i] =
+                    min
+                    (
+                        max(Rt_[i], -1.0),
+                        Rt_[i] - sqr(Rt_[i] + 1.0)/(Rt_[i] - 1.0)
+                    );
+
+                // (A:Eq. 28)
+                L_[i] =
+                    Lg[i]
+                   *sqrt(1.0 - pow(Cmu0, 6.0)*pow(CbUnstable, -2.0)*Rt_[i]);
+            }
+        }
+        else
+        {
+            L_[i] = Lcanopy[i];
+        }
+    }
+
+    // Limit characteristic length scale
+    L_ = min(L_, Lmax_);
+}
+
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+template<class BasicTurbulenceModel>
+void kL<BasicTurbulenceModel>::correctNut()
+{
+    this->nut_ = Cmu()*sqrt(k_)*L_;
+    this->nut_.correctBoundaryConditions();
+    fv::options::New(this->mesh_).correct(this->nut_);
+
+    BasicTurbulenceModel::correctNut();
+}
+
+
+template<class BasicTurbulenceModel>
+tmp<fvScalarMatrix> kL<BasicTurbulenceModel>::kSource() const
+{
+    return tmp<fvScalarMatrix>::New
+    (
+        k_,
+        dimVolume*this->rho_.dimensions()*k_.dimensions()/dimTime
+    );
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class BasicTurbulenceModel>
+kL<BasicTurbulenceModel>::kL
+(
+    const alphaField& alpha,
+    const rhoField& rho,
+    const volVectorField& U,
+    const surfaceScalarField& alphaRhoPhi,
+    const surfaceScalarField& phi,
+    const transportModel& transport,
+    const word& propertiesName,
+    const word& type
+)
+:
+    eddyViscosity<RASModel<BasicTurbulenceModel>>
+    (
+        type,
+        alpha,
+        rho,
+        U,
+        alphaRhoPhi,
+        phi,
+        transport,
+        propertiesName
+    ),
+
+    kappa_
+    (
+        dimensioned<scalar>::getOrAddToDict
+        (
+            "kappa",
+            this->coeffDict_,
+            0.41
+        )
+    ),
+    sigmak_
+    (
+        dimensioned<scalar>::getOrAddToDict
+        (
+            "sigmak",
+            this->coeffDict_,
+            1.0
+        )
+    ),
+    beta_
+    (
+        dimensioned<scalar>::getOrAddToDict
+        (
+            "beta",
+            this->coeffDict_,
+            dimless/dimTemperature,
+            3.3e-03
+        )
+    ),
+    Cmu0_
+    (
+        dimensioned<scalar>::getOrAddToDict
+        (
+            "Cmu0",
+            this->coeffDict_,
+            0.556
+        )
+    ),
+    Lmax_
+    (
+        dimensioned<scalar>::getOrAddToDict
+        (
+            "Lmax",
+            this->coeffDict_,
+            dimLength,
+            GREAT
+        )
+    ),
+    CbStable_
+    (
+        dimensioned<scalar>::getOrAddToDict
+        (
+            "CbStable",
+            this->coeffDict_,
+            0.25
+        )
+    ),
+    CbUnstable_
+    (
+        dimensioned<scalar>::getOrAddToDict
+        (
+            "CbUnstable",
+            this->coeffDict_,
+            0.35
+        )
+    ),
+
+    k_
+    (
+        IOobject
+        (
+            IOobject::groupName("k", alphaRhoPhi.group()),
+            this->runTime_.timeName(),
+            this->mesh_,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        this->mesh_
+    ),
+    L_
+    (
+        IOobject
+        (
+            IOobject::groupName("L", alphaRhoPhi.group()),
+            this->runTime_.timeName(),
+            this->mesh_,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        this->mesh_,
+        dimensionedScalar(dimLength, scalar(1))
+    ),
+    Rt_
+    (
+        IOobject
+        (
+            IOobject::groupName("Rt", alphaRhoPhi.group()),
+            this->runTime_.timeName(),
+            this->mesh_,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        this->mesh_,
+        dimensionedScalar(dimless, Zero)
+    ),
+    g_(meshObjects::gravity::New(this->mesh_.time())),
+    y_(wallDist::New(this->mesh_).y())
+{
+    bound(k_, this->kMin_);
+
+    if (type == typeName)
+    {
+        this->printCoeffs(type);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class BasicTurbulenceModel>
+bool kL<BasicTurbulenceModel>::read()
+{
+    if (eddyViscosity<RASModel<BasicTurbulenceModel>>::read())
+    {
+        kappa_.readIfPresent(this->coeffDict());
+        sigmak_.readIfPresent(this->coeffDict());
+        beta_.readIfPresent(this->coeffDict());
+        Cmu0_.readIfPresent(this->coeffDict());
+        Lmax_.readIfPresent(this->coeffDict());
+        CbStable_.readIfPresent(this->coeffDict());
+        CbUnstable_.readIfPresent(this->coeffDict());
+
+        return true;
+    }
+
+    return false;
+}
+
+
+template<class BasicTurbulenceModel>
+void kL<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_;
+    const volVectorField& U = this->U_;
+    const volScalarField& nut = this->nut_;
+
+    fv::options& fvOptions(fv::options::New(this->mesh_));
+
+    eddyViscosity<RASModel<BasicTurbulenceModel>>::correct();
+
+    // Turbulent kinetic energy production rate
+    tmp<volTensorField> tgradU = fvc::grad(U);
+    const volScalarField::Internal G
+    (
+        this->GName(),
+        nut.v()*2*magSqr(dev(symm(tgradU.cref().v())))
+    );
+    tgradU.clear();
+
+    // Square of Brunt-Vaisala (buoyancy) frequency
+    const auto& T = U.mesh().lookupObject<volScalarField>("T");
+    tmp<volScalarField> tfBV = -beta_*(fvc::grad(T) & g_);
+    const volScalarField& fBV = tfBV.cref();
+
+    // Sink or source of TKE depending on stratification type (A:Eq. 15)
+    tmp<volScalarField> tPb = -fBV*nutPrime();
+    const volScalarField& Pb = tPb.cref();
+
+    // Turbulent kinetic energy dissipation rate due to plains and canopy
+    tmp<volScalarField> tepsilon = epsilon();
+    const volScalarField& epsilon = tepsilon.cref();
+
+    // Divergence of velocity
+    tmp<volScalarField> tdivU = fvc::div(fvc::absolute(this->phi(), U));
+    const volScalarField::Internal& divU = tdivU.cref().v();
+
+    // Turbulent kinetic energy equation
+    tmp<fvScalarMatrix> kEqn
+    (
+        fvm::ddt(alpha, rho, k_)
+      + fvm::div(alphaRhoPhi, k_)
+      - fvm::laplacian(alpha*rho*DkEff(), k_)
+     ==
+        alpha()*rho()*G
+      + fvm::SuSp((Pb - epsilon)/k_, k_)
+      - fvm::SuSp((2.0/3.0)*alpha()*rho()*divU, k_)
+      + kSource()
+      + fvOptions(alpha, rho, k_)
+    );
+
+    tdivU.clear();
+    tPb.clear();
+
+    kEqn.ref().relax();
+    fvOptions.constrain(kEqn.ref());
+    solve(kEqn);
+    fvOptions.correct(k_);
+    bound(k_, this->kMin_);
+
+    // Turbulent Richardson number (A:Eq. 29)
+    Rt_ = fBV*sqr(k_/tepsilon);
+
+    stratification(fBV);
+    tfBV.clear();
+
+    correctNut();
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace RASModels
+} // End namespace Foam
+
+// ************************************************************************* //

src/atmosphericModels/turbulenceModels/RAS/kL/kL.H

@@ -0,0 +1,313 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Class
+    Foam::RASModels::kL
+
+Group
+    grpRASTurbulence
+
+Description
+    A one-equation (turbulent kinetic energy \c k) turbulence closure
+    model for incompressible and compressible geophysical applications.
+
+    Turbulent kinetic energy (\c k) is computed with a transport equation
+    and the turbulent length scale (\c L) is computed with an algebraic
+    expression which depends on the local stratification.
+
+    References:
+    \verbatim
+        Standard model (tag:A):
+            Axell, L. B., & Liungman, O. (2001).
+            A one-equation turbulence model for geophysical applications:
+            comparison with data and the k−ε model.
+            Environmental Fluid Mechanics, 1(1), 71-106.
+            DOI:10.1023/A:1011560202388
+
+        Canopy-related models (tag:W):
+            Wilson, J. D., & Flesch, T. K. (1999).
+            Wind and remnant tree sway in forest cutblocks.
+            III. A windflow model to diagnose spatial variation.
+            Agricultural and Forest Meteorology, 93(4), 259-282.
+            DOI:10.1016/S0168-1923(98)00121-X
+    \endverbatim
+
+Usage
+    Example by using \c constant/turbulenceProperties:
+    \verbatim
+    RAS
+    {
+        // Mandatory entries
+        RASModel        kL;
+
+        // Optional entries
+        kLCoeffs
+        {
+            kappa       <scalar>;
+            sigmak      <scalar>;
+            beta        <scalar>;
+            Cmu0        <scalar>;
+            Lmax        <scalar>;
+            CbStable    <scalar>;
+            CbUnstable  <scalar>;
+        }
+
+        // Inherited entries
+        ...
+    }
+    \endverbatim
+
+    where the entries mean:
+    \table
+      Property   | Description                     | Type   | Reqd | Deflt
+      RASModel   | Type name: kL                   | word   | yes  | -
+      kappa      | von Karman constant             | scalar | no   | 0.41
+      sigmak     | Empirical model coefficient     | scalar | no   | 1.0
+      beta | Thermal expansion coefficient [1/K]   | scalar | no   | 3.3e-3
+      Cmu0       | Empirical model coefficient     | scalar | no   | 0.556
+      Lmax   | Maximum mixing-length scale [m]     | scalar | no   | GREAT
+      CbStable   | Stable stratification constant  | scalar | no   | 0.25
+      CbUnstable | Unstable stratification constant | scalar | no  | 0.35
+    \endtable
+
+    The inherited entries are elaborated in:
+      - \link eddyViscosity.H \endlink
+
+    \heading Input fields (mandatory)
+    \plaintable
+        k         | Turbulent kinetic energy                [m2/s2]
+        T         | Potential temperature                   [K]
+    \endplaintable
+
+    \heading Input fields (optional)
+    \plaintable
+        canopyHeight    | Canopy height                     [m]
+        plantCd         | Plant canopy drag coefficient     [-]
+        leafAreaDensity | Leaf area density                 [1/m]
+        Rt              | Turbulent Richardson number       [-]
+        L               | Characteristic length scale       [m]
+    \endplaintable
+
+Note
+  - Optional input fields can/should be input
+    by using \c readFields function object.
+
+SourceFiles
+    kL.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef kL_H
+#define kL_H
+
+#include "RASModel.H"
+#include "eddyViscosity.H"
+#include "uniformDimensionedFields.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace RASModels
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class kL Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class BasicTurbulenceModel>
+class kL
+:
+    public eddyViscosity<RASModel<BasicTurbulenceModel>>
+{
+    // Private Member Functions
+
+        //- Return modified stability function of Launder,
+        //- calibrated with the results of Hogstrom (A:Eq. 31)
+        tmp<volScalarField> Cmu() const;
+
+        //- Return modified stability function of Launder,
+        //- calibrated with the results of Hogstrom (A:Eq. 32)
+        tmp<volScalarField> CmuPrime() const;
+
+        //- Return eddy diffusivity (A:Eq. 12)
+        tmp<volScalarField> nutPrime() const;
+
+        //- Return turbulent kinetic energy dissipation rate due to canopy
+        tmp<volScalarField> epsilonCanopy() const;
+
+        //- Return turbulent kinetic energy
+        //- dissipation rate due to plains and canopy
+        tmp<volScalarField> epsilon() const;
+
+        //- Return canopy height
+        tmp<volScalarField> canopyHeight() const;
+
+        //- Return characteristic length scale
+        tmp<volScalarField> L() const;
+
+        //- Modify characteristic length scale
+        //- according to the specified stratification
+        void stratification(const volScalarField& fVB);
+
+        // Generated Methods
+
+            //- No copy construct
+            kL(const kL&) = delete;
+
+            //- No copy assignment
+            void operator=(const kL&) = delete;
+
+
+protected:
+
+    // Protected Data
+
+        // Model coefficients
+
+            //- von Karman constant
+            dimensionedScalar kappa_;
+
+            //- Empirical model coefficient
+            dimensionedScalar sigmak_;
+
+            //- Thermal expansion coefficient [1/K]
+            dimensionedScalar beta_;
+
+            //- Empirical model coefficient
+            dimensionedScalar Cmu0_;
+
+            //- Maximum mixing-length scalar [m]
+            dimensionedScalar Lmax_;
+
+            //- Stable stratification constant
+            dimensionedScalar CbStable_;
+
+            //- Unstable stratification constant
+            dimensionedScalar CbUnstable_;
+
+
+        // Fields
+
+            //- Turbulent kinetic energy    [m2/s2]
+            volScalarField k_;
+
+            //- Characteristic length scale [m]
+            volScalarField L_;
+
+            //- Turbulent Richardson number [-]
+            volScalarField Rt_;
+
+            //- Gravitational acceleration  [m2/s2]
+            const uniformDimensionedVectorField& g_;
+
+            //- Wall distance
+            //  Note: different to wall distance in parent RASModel
+            //  which is for near-wall cells only
+            const volScalarField& y_;
+
+
+    // Protected Member Functions
+
+        //- Correct the turbulence viscosity
+        virtual void correctNut();
+
+        //- Add explicit source for turbulent kinetic energy
+        virtual tmp<fvScalarMatrix> kSource() const;
+
+
+public:
+
+    typedef typename BasicTurbulenceModel::alphaField alphaField;
+    typedef typename BasicTurbulenceModel::rhoField rhoField;
+    typedef typename BasicTurbulenceModel::transportModel transportModel;
+
+
+    //- Runtime type information
+    TypeName("kL");
+
+
+    // Constructors
+
+        //- Construct from components
+        kL
+        (
+            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 ~kL() = default;
+
+
+    // Member Functions
+
+        //- Re-read model coefficients if they have changed
+        virtual bool read();
+
+        //- Return the effective diffusivity for k
+        tmp<volScalarField> DkEff() const
+        {
+            return tmp<volScalarField>::New
+            (
+                "DkEff",
+                (this->nut_/sigmak_ + this->nu())
+            );
+        }
+
+        //- Return the turbulent kinetic energy field
+        virtual tmp<volScalarField> k() const
+        {
+            return k_;
+        }
+
+        //- Solve the turbulence equations and correct the turbulence viscosity
+        virtual void correct();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace RASModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+    #include "kL.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/atmosphericModels/turbulenceModels/atmosphericTurbulentFluidThermoModels.C

@@ -0,0 +1,40 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | www.openfoam.com
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+    Copyright (C) 2021 OpenCFD Ltd.
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "turbulentFluidThermoModels.H"
+
+// -------------------------------------------------------------------------- //
+// RAS models
+// -------------------------------------------------------------------------- //
+
+#include "kEpsilonLopesdaCosta.H"
+makeRASModel(kEpsilonLopesdaCosta);
+
+#include "kL.H"
+makeRASModel(kL);
+
+// ************************************************************************* //

src/atmosphericModels/turbulenceModels/atmosphericTurbulentTransportModels.C

@@ -6,6 +6,7 @@
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
     Copyright (C) 2018 OpenFOAM Foundation
+    Copyright (C) 2021 OpenCFD Ltd.
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
@@ -34,4 +35,7 @@ License
 #include "kEpsilonLopesdaCosta.H"
 makeRASModel(kEpsilonLopesdaCosta);
 
+#include "kL.H"
+makeRASModel(kL);
+
 // ************************************************************************* //

tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/Allrun

@@ -10,6 +10,7 @@ cd "${0%/*}" || exit                                # Run from this directory
     setups="
     kEpsilon
     kOmegaSST
+    kL
     "
 
     # flag to enable computations in parallel mode

tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/plot

@@ -7,8 +7,9 @@ cd "${0%/*}" || exit                                # Run from this directory
 
     # operand setups
     setups="
-    kEpsilon-neutral-stability
+    kEpsilon
-    kOmegaSST-neutral-stability
+    kOmegaSST
+    kL
     "
 
 

tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/setups.orig/kL/0.orig/T

@@ -0,0 +1 @@
+../../common/0.orig/T

tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/setups.orig/kL/0.orig/U

@@ -0,0 +1 @@
+../../common/0.orig/U

tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/setups.orig/kL/0.orig/alphat

@@ -0,0 +1 @@
+../../common/0.orig/alphat

tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/setups.orig/kL/0.orig/k

@@ -0,0 +1 @@
+../../common/0.orig/k

tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/setups.orig/kL/0.orig/leafAreaDensity

@@ -0,0 +1 @@
+../../common/0.orig/leafAreaDensity

tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/setups.orig/kL/0.orig/nut

@@ -0,0 +1 @@
+../../common/0.orig/nut

tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/setups.orig/kL/0.orig/p_rgh

@@ -0,0 +1 @@
+../../common/0.orig/p_rgh

tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/setups.orig/kL/0.orig/plantCd

@@ -0,0 +1 @@
+../../common/0.orig/plantCd

tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/setups.orig/kL/0.orig/qPlant

@@ -0,0 +1,36 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2106                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      qPlant;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -3 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    "(bottom|top)"
+    {
+        type            fixedValue;
+        value           uniform 0;
+    }
+
+    "(inlet|outlet|left|right)"
+    {
+        type            cyclic;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/setups.orig/kL/constant/fvOptions

@@ -0,0 +1,57 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2106                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      fvOptions;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+pressureGradient
+{
+    type            vectorSemiImplicitSource;
+    volumeMode      specific;
+    selectionMode   all;
+    injectionRateSuSp
+    {
+        U               ( ( 0 0.00197805 0 ) 0 );
+    }
+}
+
+atmCoriolisUSource1
+{
+    type            atmCoriolisUSource;
+    atmCoriolisUSourceCoeffs
+    {
+        selectionMode   all;
+        Omega           ( 0 0 5.65156e-05 );
+    }
+}
+
+atmPlantCanopyUSource1
+{
+    type            atmPlantCanopyUSource;
+    atmPlantCanopyUSourceCoeffs
+    {
+        selectionMode   all;
+    }
+}
+
+atmPlantCanopyTSource1
+{
+    type            atmPlantCanopyTSource;
+    atmPlantCanopyTSourceCoeffs
+    {
+        selectionMode   all;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/setups.orig/kL/constant/g

@@ -0,0 +1 @@
+../../common/constant/g

tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/setups.orig/kL/constant/transportProperties

@@ -0,0 +1 @@
+../../common/constant/transportProperties

tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/setups.orig/kL/constant/turbulenceProperties

@@ -0,0 +1,36 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2106                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      turbulenceProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+simulationType      RAS;
+
+RAS
+{
+    RASModel        kL;
+
+    turbulence      on;
+
+    printCoeffs     on;
+
+    kLCoeffs
+    {
+        CbStable    0.25;
+        CbUnstable  0.35;
+        Lmax        41.0;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/precursor/setups.orig/kL/system

@@ -0,0 +1 @@
+../common/system

tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/Allrun

@@ -10,6 +10,7 @@ cd "${0%/*}" || exit                                # Run from this directory
     setups="
     kEpsilon
     kOmegaSST
+    kL
     "
 
     # flag to enable computations in parallel mode

tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/setups.orig/kL/0.orig/T

@@ -0,0 +1 @@
+../../common/0.orig/T

tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/setups.orig/kL/0.orig/U

@@ -0,0 +1 @@
+../../common/0.orig/U

tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/setups.orig/kL/0.orig/alphat

@@ -0,0 +1 @@
+../../common/0.orig/alphat

tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/setups.orig/kL/0.orig/k

@@ -0,0 +1 @@
+../../common/0.orig/k

tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/setups.orig/kL/0.orig/nut

@@ -0,0 +1 @@
+../../common/0.orig/nut

tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/setups.orig/kL/0.orig/p_rgh

@@ -0,0 +1 @@
+../../common/0.orig/p_rgh

tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/setups.orig/kL/constant/boundaryData

@@ -0,0 +1 @@
+../../common/constant/boundaryData

tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/setups.orig/kL/constant/fvOptions

@@ -0,0 +1,39 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2106                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      fvOptions;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+pressureGradient
+{
+    type                vectorSemiImplicitSource;
+    selectionMode       all;
+    volumeMode          specific;
+    injectionRateSuSp
+    {
+        U ((0 1.978046e-03 0) 0);
+    }
+}
+
+atmCoriolisUSource1
+{
+    type               atmCoriolisUSource;
+    atmCoriolisUSourceCoeffs
+    {
+        selectionMode   all;
+        Omega           (0 0 5.65156e-5);
+    }
+}
+
+
+// ************************************************************************* //

tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/setups.orig/kL/constant/g

@@ -0,0 +1 @@
+../../common/constant/g

tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/setups.orig/kL/constant/transportProperties

@@ -0,0 +1 @@
+../../common/constant/transportProperties

tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/setups.orig/kL/constant/turbulenceProperties

@@ -0,0 +1,29 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  v2106                                 |
+|   \\  /    A nd           | Website:  www.openfoam.com                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      turbulenceProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+simulationType        RAS;
+
+RAS
+{
+    RASModel          kL;
+
+    turbulence        on;
+
+    printCoeffs       on;
+}
+
+
+// ************************************************************************* //

tutorials/verificationAndValidation/atmosphericModels/atmFlatTerrain/successor/setups.orig/kL/system

@@ -0,0 +1 @@
+../common/system/