TurbulenceModels: LienCubicKE, ShihQuadraticKE and LienLeschziner models rewritten

The implementation now correspond to the definitions in the readily
available reference:

http://personalpages.manchester.ac.uk/staff/david.d.apsley/specturb.pdf

in which a large number of linear and non-linear models are presented in
a clean and consistent manner.  This has made re-implementation and
checking far easier than working from the original references which anyway
are no longer available for the LienCubicKE and ShihQuadraticKE models.

applications/test/tensor/Test-tensor.C

@@ -57,19 +57,18 @@ int main()
     Info<< "Check symmetric transformation "
         << transform(t1, st1) << endl;
 
-    Info<< "Check for dot product of symmetric tensors "
+    Info<< "Check dot product of symmetric tensors "
         << (st1 & st2) << endl;
 
-    Info<< "Check for inner sqr of a symmetric tensor "
+    Info<< "Check inner sqr of a symmetric tensor "
         << innerSqr(st1) << " " << innerSqr(st1) - (st1 & st1) << endl;
 
-    Info<< "Check for symmetric part of dot product of symmetric tensors "
+    Info<< "Check symmetric part of dot product of symmetric tensors "
         << twoSymm(st1&st2) - ((st1&st2) + (st2&st1)) << endl;
 
     tensor sk1 = skew(t6);
-    tensor sk2 = skew(t7);
+    Info<< "Check dot product of symmetric and skew tensors "
-    Info<< "Check for symmetric part of dot product of skew tensors "
+        << twoSymm(st1&sk1) - ((st1&sk1) - (sk1&st1)) << endl;
-        << twoSymm(sk1&sk2) - ((sk1&sk2) - (sk2&sk1)) << endl;
 
     vector v1(1, 2, 3);
 

src/TurbulenceModels/incompressible/Make/files

@@ -4,13 +4,11 @@ turbulentTransportModels/turbulentTransportModels.C
 turbulentTransportModels/RAS/qZeta/qZeta.C
 turbulentTransportModels/RAS/kkLOmega/kkLOmega.C
 turbulentTransportModels/RAS/LamBremhorstKE/LamBremhorstKE.C
-turbulentTransportModels/RAS/LienCubicKE/LienCubicKE.C
+turbulentTransportModels/RAS/LienLeschziner/LienLeschziner.C
-turbulentTransportModels/RAS/LienCubicKELowRe/LienCubicKELowRe.C
-turbulentTransportModels/RAS/LienLeschzinerLowRe/LienLeschzinerLowRe.C
 turbulentTransportModels/RAS/ShihQuadraticKE/ShihQuadraticKE.C
+turbulentTransportModels/RAS/LienCubicKE/LienCubicKE.C
 
 BCs = turbulentTransportModels/RAS/derivedFvPatchFields
-
 turbulentTransportModels/derivedFvPatchFields/wallFunctions/alphatWallFunctions/alphatJayatillekeWallFunction/alphatJayatillekeWallFunctionFvPatchScalarField.C
 
 LIB = $(FOAM_LIBBIN)/libincompressibleTurbulenceModels

src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/LamBremhorstKE/LamBremhorstKE.C

@@ -183,13 +183,13 @@ bool LamBremhorstKE::read()
 
 void LamBremhorstKE::correct()
 {
-    eddyViscosity<incompressible::RASModel>::correct();
-
     if (!turbulence_)
     {
         return;
     }
 
+    eddyViscosity<incompressible::RASModel>::correct();
+
     volScalarField G(GName(), nut_*2*magSqr(symm(fvc::grad(U_))));
 
 
@@ -205,7 +205,6 @@ void LamBremhorstKE::correct()
 
 
     // Dissipation equation
-
     tmp<fvScalarMatrix> epsEqn
     (
         fvm::ddt(epsilon_)
@@ -222,7 +221,6 @@ void LamBremhorstKE::correct()
 
 
     // Turbulent kinetic energy equation
-
     tmp<fvScalarMatrix> kEqn
     (
         fvm::ddt(k_)

src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/LamBremhorstKE/LamBremhorstKE.H

@@ -85,7 +85,6 @@ protected:
         const volScalarField& y_;
 
         volScalarField Rt_;
-
         volScalarField fMu_;
 
 

src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/LienCubicKE/LienCubicKE.C

@@ -1,4 +1,4 @@
-/*---------------------------------------------------------------------------*\
+/*---------------------------------------------------------------------------* \
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
@@ -24,6 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "LienCubicKE.H"
+#include "wallDist.H"
 #include "bound.H"
 #include "addToRunTimeSelectionTable.H"
 
@@ -43,52 +44,76 @@ addToRunTimeSelectionTable(RASModel, LienCubicKE, dictionary);
 
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
+tmp<volScalarField> LienCubicKE::fMu() const
+{
+    const volScalarField yStar(sqrt(k_)*y_/nu());
+
+    return
+        (scalar(1) - exp(-Anu_*yStar))
+       *(scalar(1) + (2*kappa_/(pow(Cmu_, 0.75))/(yStar + SMALL)));
+}
+
+
+tmp<volScalarField> LienCubicKE::f2() const
+{
+    tmp<volScalarField> Rt = sqr(k_)/(nu()*epsilon_);
+
+    return scalar(1) - 0.3*exp(-sqr(Rt));
+}
+
+
+tmp<volScalarField> LienCubicKE::E(const volScalarField& f2) const
+{
+    const volScalarField yStar(sqrt(k_)*y_/nu());
+    const volScalarField le
+    (
+        kappa_*y_/(scalar(1) + (2*kappa_/(pow(Cmu_, 0.75))/(yStar + SMALL)))
+    );
+
+    return
+        (Ceps2_*pow(Cmu_, 0.75))
+       *(f2*sqrt(k_)*epsilon_/le)*exp(-AE_*sqr(yStar));
+}
+
+
 void LienCubicKE::correctNut()
 {
-    nut_ =
+    correctNonlinearStress(fvc::grad(U_));
-        Cmu_*sqr(k_)/epsilon_
-        // C5 term, implicit
-      + max
-        (
-            C5viscosity_,
-            dimensionedScalar("0", C5viscosity_.dimensions(), 0.0)
-        );
-
-    nut_.correctBoundaryConditions();
 }
 
 
 void LienCubicKE::correctNonlinearStress(const volTensorField& gradU)
 {
-    nonlinearStress_ = symm
+    volSymmTensorField S(symm(gradU));
-    (
+    volTensorField W(skew(gradU));
-        // quadratic terms
+
-        pow3(k_)/sqr(epsilon_)
+    volScalarField sBar((k_/epsilon_)*sqrt(2.0)*mag(S));
+    volScalarField wBar((k_/epsilon_)*sqrt(2.0)*mag(W));
+
+    volScalarField Cmu((2.0/3.0)/(Cmu1_ + sBar + Cmu2_*wBar));
+    volScalarField fMu(this->fMu());
+
+    nut_ = Cmu*fMu*sqr(k_)/epsilon_;
+    nut_.correctBoundaryConditions();
+
+    nonlinearStress_ =
+        fMu*k_
        *(
-            Ctau1_/fEta_
+            // Quadratic terms
+            sqr(k_/epsilon_)/(Cbeta_ + pow3(sBar))
            *(
-                (gradU & gradU)
+                Cbeta1_*dev(innerSqr(S))
-              + (gradU & gradU)().T()
+              + Cbeta2_*twoSymm(S&W)
+              + Cbeta3_*dev(symm(W&W))
             )
-          + Ctau2_/fEta_*(gradU & gradU.T())
+
-          + Ctau3_/fEta_*(gradU.T() & gradU)
+            // Cubic terms
-        )
+          - pow3(Cmu*k_/epsilon_)
-        // cubic term C4
+           *(
-      - 20.0*pow4(k_)/pow3(epsilon_)
+                (Cgamma1_*magSqr(S) - Cgamma2_*magSqr(W))*S
-       *pow3(Cmu_)
+              + Cgamma4_*twoSymm((innerSqr(S)&W))
-       *(
+            )
-            ((gradU & gradU) & gradU.T())
+        );
-          + ((gradU & gradU.T()) & gradU.T())
-          - ((gradU.T() & gradU) & gradU)
-          - ((gradU.T() & gradU.T()) & gradU)
-        )
-        // cubic term C5, explicit part
-      + min
-        (
-            C5viscosity_,
-            dimensionedScalar("0", C5viscosity_.dimensions(), 0.0)
-        )*gradU
-    );
 }
 
 
@@ -118,20 +143,20 @@ LienCubicKE::LienCubicKE
         propertiesName
     ),
 
-    C1_
+    Ceps1_
     (
         dimensioned<scalar>::lookupOrAddToDict
         (
-            "C1",
+            "Ceps1",
             coeffDict_,
             1.44
         )
     ),
-    C2_
+    Ceps2_
     (
         dimensioned<scalar>::lookupOrAddToDict
         (
-            "C2",
+            "Ceps2",
             coeffDict_,
             1.92
         )
@@ -154,58 +179,121 @@ LienCubicKE::LienCubicKE
             1.3
         )
     ),
-    A1_
+    Cmu1_
     (
         dimensioned<scalar>::lookupOrAddToDict
         (
-            "A1",
+            "Cmu1",
             coeffDict_,
             1.25
         )
     ),
-    A2_
+    Cmu2_
     (
         dimensioned<scalar>::lookupOrAddToDict
         (
-            "A2",
+            "Cmu2",
+            coeffDict_,
+            0.9
+        )
+    ),
+    Cbeta_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "Cbeta",
             coeffDict_,
             1000.0
         )
     ),
-    Ctau1_
+    Cbeta1_
     (
         dimensioned<scalar>::lookupOrAddToDict
         (
-            "Ctau1",
+            "Cbeta1",
             coeffDict_,
-            -4.0
+            3.0
         )
     ),
-    Ctau2_
+    Cbeta2_
     (
         dimensioned<scalar>::lookupOrAddToDict
         (
-            "Ctau2",
+            "Cbeta2",
             coeffDict_,
-            13.0
+            15.0
         )
     ),
-    Ctau3_
+    Cbeta3_
     (
         dimensioned<scalar>::lookupOrAddToDict
         (
-            "Ctau3",
+            "Cbeta3",
             coeffDict_,
-            -2.0
+            -19.0
         )
     ),
-    alphaKsi_
+    Cgamma1_
     (
         dimensioned<scalar>::lookupOrAddToDict
         (
-            "alphaKsi",
+            "Cgamma1",
             coeffDict_,
-            0.9
+            16.0
+        )
+    ),
+    Cgamma2_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "Cgamma2",
+            coeffDict_,
+            16.0
+        )
+    ),
+    Cgamma4_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "Cgamma4",
+            coeffDict_,
+            -80.0
+        )
+    ),
+    Cmu_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "Cmu",
+            coeffDict_,
+            0.09
+        )
+    ),
+    kappa_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "kappa",
+            coeffDict_,
+            0.41
+        )
+    ),
+    Anu_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "Anu",
+            coeffDict_,
+            0.0198
+        )
+    ),
+    AE_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "AE",
+            coeffDict_,
+            0.00375
         )
     ),
 
@@ -221,6 +309,7 @@ LienCubicKE::LienCubicKE
         ),
         mesh_
     ),
+
     epsilon_
     (
         IOobject
@@ -234,35 +323,14 @@ LienCubicKE::LienCubicKE
         mesh_
     ),
 
-    eta_
+    y_(wallDist::New(mesh_).y())
-    (
-        k_/bound(epsilon_, epsilonMin_)
-       *sqrt(2.0*magSqr(0.5*(fvc::grad(U) + T(fvc::grad(U)))))
-    ),
-    ksi_
-    (
-        k_/epsilon_
-       *sqrt(2.0*magSqr(0.5*(fvc::grad(U) - T(fvc::grad(U)))))
-    ),
-    Cmu_(2.0/(3.0*(A1_ + eta_ + alphaKsi_*ksi_))),
-    fEta_(A2_ + pow3(eta_)),
-
-    C5viscosity_
-    (
-       -2.0*pow3(Cmu_)*pow4(k_)/pow3(epsilon_)
-       *(
-            magSqr(fvc::grad(U) + T(fvc::grad(U)))
-          - magSqr(fvc::grad(U) - T(fvc::grad(U)))
-        )
-    )
 {
     bound(k_, kMin_);
-    // already bounded: bound(epsilon_, epsilonMin_);
+    bound(epsilon_, epsilonMin_);
 
     if (type == typeName)
     {
         correctNut();
-        correctNonlinearStress(fvc::grad(U));
         printCoeffs(type);
     }
 }
@@ -274,16 +342,23 @@ bool LienCubicKE::read()
 {
     if (nonlinearEddyViscosity<incompressible::RASModel>::read())
     {
-        C1_.readIfPresent(coeffDict());
+        Ceps1_.readIfPresent(coeffDict());
-        C2_.readIfPresent(coeffDict());
+        Ceps2_.readIfPresent(coeffDict());
         sigmak_.readIfPresent(coeffDict());
         sigmaEps_.readIfPresent(coeffDict());
-        A1_.readIfPresent(coeffDict());
+        Cmu1_.readIfPresent(coeffDict());
-        A2_.readIfPresent(coeffDict());
+        Cmu2_.readIfPresent(coeffDict());
-        Ctau1_.readIfPresent(coeffDict());
+        Cbeta_.readIfPresent(coeffDict());
-        Ctau2_.readIfPresent(coeffDict());
+        Cbeta1_.readIfPresent(coeffDict());
-        Ctau3_.readIfPresent(coeffDict());
+        Cbeta2_.readIfPresent(coeffDict());
-        alphaKsi_.readIfPresent(coeffDict());
+        Cbeta3_.readIfPresent(coeffDict());
+        Cgamma1_.readIfPresent(coeffDict());
+        Cgamma2_.readIfPresent(coeffDict());
+        Cgamma4_.readIfPresent(coeffDict());
+        Cmu_.readIfPresent(coeffDict());
+        kappa_.readIfPresent(coeffDict());
+        Anu_.readIfPresent(coeffDict());
+        AE_.readIfPresent(coeffDict());
 
         return true;
     }
@@ -296,13 +371,13 @@ bool LienCubicKE::read()
 
 void LienCubicKE::correct()
 {
-    nonlinearEddyViscosity<incompressible::RASModel>::correct();
-
     if (!turbulence_)
     {
         return;
     }
 
+    nonlinearEddyViscosity<incompressible::RASModel>::correct();
+
     tmp<volTensorField> tgradU = fvc::grad(U_);
     const volTensorField& gradU = tgradU();
 
@@ -316,6 +391,8 @@ void LienCubicKE::correct()
     // Update epsilon and G at the wall
     epsilon_.boundaryField().updateCoeffs();
 
+    const volScalarField f2(this->f2());
+
     // Dissipation equation
     tmp<fvScalarMatrix> epsEqn
     (
@@ -323,8 +400,9 @@ void LienCubicKE::correct()
       + fvm::div(phi_, epsilon_)
       - fvm::laplacian(DepsilonEff(), epsilon_)
       ==
-        C1_*G*epsilon_/k_
+        Ceps1_*G*epsilon_/k_
-      - fvm::Sp(C2_*epsilon_/k_, epsilon_)
+      - fvm::Sp(Ceps2_*f2*epsilon_/k_, epsilon_)
+      + E(f2)
     );
 
     epsEqn().relax();
@@ -349,18 +427,7 @@ void LienCubicKE::correct()
     bound(k_, kMin_);
 
 
-    // Re-calculate viscosity
+    // Re-calculate viscosity and non-linear stress
-
-    eta_ = k_/epsilon_*sqrt(2.0*magSqr(0.5*(gradU + gradU.T())));
-    ksi_ = k_/epsilon_*sqrt(2.0*magSqr(0.5*(gradU - gradU.T())));
-    Cmu_ = 2.0/(3.0*(A1_ + eta_ + alphaKsi_*ksi_));
-    fEta_ = A2_ + pow3(eta_);
-
-    C5viscosity_ =
-       -2.0*pow3(Cmu_)*pow4(k_)/pow3(epsilon_)
-       *(magSqr(gradU + gradU.T()) - magSqr(gradU - gradU.T()));
-
-    correctNut();
     correctNonlinearStress(gradU);
 }
 

src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/LienCubicKE/LienCubicKE.H

@@ -28,7 +28,8 @@ Group
     grpIcoRASTurbulence
 
 Description
-    Lien cubic non-linear k-epsilon turbulence model for incompressible flows.
+    Lien cubic non-linear low-Reynolds k-epsilon turbulence models for
+    incompressible flows.
 
     This turbulence model is described in:
     \verbatim
@@ -38,8 +39,17 @@ Description
         Engineering Turbulence Modelling and Experiments 3, 91-100.
     \endverbatim
 
+    Implemented according to the specification in:
+    <a href=
+    "http://personalpages.manchester.ac.uk/staff/david.d.apsley/specturb.pdf"
+    >Apsley: Turbulence Models 2002</a>
+
+    In addition to the low-Reynolds number damping functions support for
+    wall-functions is also included to allow for low- and high-Reynolds number
+    operation.
+
 See Also
-    Foam::incompressible::RASModels::LienCubicKELowRe
+    Foam::incompressible::RASModels::ShihQuadraticKE
 
 SourceFiles
     LienCubicKE.C
@@ -76,16 +86,25 @@ protected:
 
         // Model coefficients
 
-            dimensionedScalar C1_;
+            dimensionedScalar Ceps1_;
-            dimensionedScalar C2_;
+            dimensionedScalar Ceps2_;
             dimensionedScalar sigmak_;
             dimensionedScalar sigmaEps_;
-            dimensionedScalar A1_;
+            dimensionedScalar Cmu1_;
-            dimensionedScalar A2_;
+            dimensionedScalar Cmu2_;
-            dimensionedScalar Ctau1_;
+            dimensionedScalar Cbeta_;
-            dimensionedScalar Ctau2_;
+            dimensionedScalar Cbeta1_;
-            dimensionedScalar Ctau3_;
+            dimensionedScalar Cbeta2_;
-            dimensionedScalar alphaKsi_;
+            dimensionedScalar Cbeta3_;
+            dimensionedScalar Cgamma1_;
+            dimensionedScalar Cgamma2_;
+            dimensionedScalar Cgamma4_;
+
+            dimensionedScalar Cmu_;
+            dimensionedScalar kappa_;
+
+            dimensionedScalar Anu_;
+            dimensionedScalar AE_;
 
 
         // Fields
@@ -93,15 +112,18 @@ protected:
             volScalarField k_;
             volScalarField epsilon_;
 
-            volScalarField eta_;
+            //- Wall distance
-            volScalarField ksi_;
+            //  Note: different to wall distance in parent RASModel
-            volScalarField Cmu_;
+            //  which is for near-wall cells only
-            volScalarField fEta_;
+            const volScalarField& y_;
-            volScalarField C5viscosity_;
 
 
     // Protected Member Functions
 
+        tmp<volScalarField> fMu() const;
+        tmp<volScalarField> f2() const;
+        tmp<volScalarField> E(const volScalarField& f2) const;
+
         virtual void correctNut();
         virtual void correctNonlinearStress(const volTensorField& gradU);
 
@@ -175,7 +197,7 @@ public:
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace RASModels
-} // End namespace incompressible
+} // Edn namespace incompressible
 } // End namespace Foam
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/LienCubicKELowRe/LienCubicKELowRe.C

@@ -1,446 +0,0 @@
-/*---------------------------------------------------------------------------* \
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
-     \\/     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 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 "LienCubicKELowRe.H"
-#include "wallDist.H"
-#include "wallFvPatch.H"
-#include "bound.H"
-#include "addToRunTimeSelectionTable.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace incompressible
-{
-namespace RASModels
-{
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-defineTypeNameAndDebug(LienCubicKELowRe, 0);
-addToRunTimeSelectionTable(RASModel, LienCubicKELowRe, dictionary);
-
-// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
-
-tmp<volScalarField> LienCubicKELowRe::fMu()
-{
-    return
-        (scalar(1) - exp(-Am_*yStar_))
-       /(scalar(1) - exp(-Aepsilon_*yStar_) + SMALL);
-}
-
-
-void LienCubicKELowRe::correctNut()
-{
-    nut_ =
-        Cmu_*fMu()*sqr(k_)/epsilon_
-        // C5 term, implicit
-      + max
-        (
-            C5viscosity_,
-            dimensionedScalar("0", C5viscosity_.dimensions(), 0.0)
-        );
-
-    nut_.correctBoundaryConditions();
-}
-
-
-void LienCubicKELowRe::correctNonlinearStress(const volTensorField& gradU)
-{
-    nonlinearStress_ = symm
-    (
-        // quadratic terms
-        pow3(k_)/sqr(epsilon_)
-       *(
-            Ctau1_/fEta_
-           *(
-                (gradU & gradU)
-              + (gradU & gradU)().T()
-            )
-          + Ctau2_/fEta_*(gradU & gradU.T())
-          + Ctau3_/fEta_*(gradU.T() & gradU)
-        )
-        // cubic term C4
-      - 20.0*pow4(k_)/pow3(epsilon_)
-       *pow3(Cmu_)
-       *(
-            ((gradU & gradU) & gradU.T())
-          + ((gradU & gradU.T()) & gradU.T())
-          - ((gradU.T() & gradU) & gradU)
-          - ((gradU.T() & gradU.T()) & gradU)
-        )
-        // cubic term C5, explicit part
-      + min
-        (
-            C5viscosity_,
-            dimensionedScalar("0", C5viscosity_.dimensions(), 0.0)
-        )*gradU
-    );
-}
-
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-LienCubicKELowRe::LienCubicKELowRe
-(
-    const geometricOneField& alpha,
-    const geometricOneField& rho,
-    const volVectorField& U,
-    const surfaceScalarField& alphaRhoPhi,
-    const surfaceScalarField& phi,
-    const transportModel& transport,
-    const word& propertiesName,
-    const word& type
-)
-:
-    nonlinearEddyViscosity<incompressible::RASModel>
-    (
-        type,
-        alpha,
-        rho,
-        U,
-        alphaRhoPhi,
-        phi,
-        transport,
-        propertiesName
-    ),
-
-    C1_
-    (
-        dimensioned<scalar>::lookupOrAddToDict
-        (
-            "C1",
-            coeffDict_,
-            1.44
-        )
-    ),
-    C2_
-    (
-        dimensioned<scalar>::lookupOrAddToDict
-        (
-            "C2",
-            coeffDict_,
-            1.92
-        )
-    ),
-    sigmak_
-    (
-        dimensioned<scalar>::lookupOrAddToDict
-        (
-            "sigmak",
-            coeffDict_,
-            1.0
-        )
-    ),
-    sigmaEps_
-    (
-        dimensioned<scalar>::lookupOrAddToDict
-        (
-            "sigmaEps",
-            coeffDict_,
-            1.3
-        )
-    ),
-    A1_
-    (
-        dimensioned<scalar>::lookupOrAddToDict
-        (
-            "A1",
-            coeffDict_,
-            1.25
-        )
-    ),
-    A2_
-    (
-        dimensioned<scalar>::lookupOrAddToDict
-        (
-            "A2",
-            coeffDict_,
-            1000.0
-        )
-    ),
-    Ctau1_
-    (
-        dimensioned<scalar>::lookupOrAddToDict
-        (
-            "Ctau1",
-            coeffDict_,
-            -4.0
-        )
-    ),
-    Ctau2_
-    (
-        dimensioned<scalar>::lookupOrAddToDict
-        (
-            "Ctau2",
-            coeffDict_,
-            13.0
-        )
-    ),
-    Ctau3_
-    (
-        dimensioned<scalar>::lookupOrAddToDict
-        (
-            "Ctau3",
-            coeffDict_,
-            -2.0
-        )
-    ),
-    alphaKsi_
-    (
-        dimensioned<scalar>::lookupOrAddToDict
-        (
-            "alphaKsi",
-            coeffDict_,
-            0.9
-        )
-    ),
-    CmuWall_
-    (
-        dimensioned<scalar>::lookupOrAddToDict
-        (
-            "Cmu",
-            coeffDict_,
-            0.09
-        )
-    ),
-    kappa_
-    (
-        dimensioned<scalar>::lookupOrAddToDict
-        (
-            "kappa",
-            coeffDict_,
-            0.41
-        )
-    ),
-    Am_
-    (
-        dimensioned<scalar>::lookupOrAddToDict
-        (
-            "Am",
-            coeffDict_,
-            0.016
-        )
-    ),
-    Aepsilon_
-    (
-        dimensioned<scalar>::lookupOrAddToDict
-        (
-            "Aepsilon",
-            coeffDict_,
-            0.263
-        )
-    ),
-    Amu_
-    (
-        dimensioned<scalar>::lookupOrAddToDict
-        (
-            "Amu",
-            coeffDict_,
-            0.00222
-        )
-    ),
-
-    k_
-    (
-        IOobject
-        (
-            IOobject::groupName("k", U.group()),
-            runTime_.timeName(),
-            mesh_,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh_
-    ),
-
-    epsilon_
-    (
-        IOobject
-        (
-            IOobject::groupName("epsilon", U.group()),
-            runTime_.timeName(),
-            mesh_,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh_
-    ),
-
-    y_(wallDist::New(mesh_).y()),
-
-    eta_
-    (
-        k_/bound(epsilon_, epsilonMin_)
-       *sqrt(2.0*magSqr(0.5*(fvc::grad(U) + T(fvc::grad(U)))))
-    ),
-    ksi_
-    (
-        k_/epsilon_
-       *sqrt(2.0*magSqr(0.5*(fvc::grad(U) - T(fvc::grad(U)))))
-    ),
-    Cmu_(2.0/(3.0*(A1_ + eta_ + alphaKsi_*ksi_))),
-    fEta_(A2_ + pow3(eta_)),
-
-    C5viscosity_
-    (
-       -2.0*pow3(Cmu_)*pow4(k_)/pow3(epsilon_)
-       *(
-            magSqr(fvc::grad(U) + T(fvc::grad(U)))
-          - magSqr(fvc::grad(U) - T(fvc::grad(U)))
-        )
-    ),
-
-    yStar_(sqrt(k_)*y_/nu() + SMALL)
-{
-    bound(k_, kMin_);
-    // already bounded: bound(epsilon_, epsilonMin_);
-
-    if (type == typeName)
-    {
-        correctNut();
-        correctNonlinearStress(fvc::grad(U));
-        printCoeffs(type);
-    }
-}
-
-
-// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-
-bool LienCubicKELowRe::read()
-{
-    if (nonlinearEddyViscosity<incompressible::RASModel>::read())
-    {
-        C1_.readIfPresent(coeffDict());
-        C2_.readIfPresent(coeffDict());
-        sigmak_.readIfPresent(coeffDict());
-        sigmaEps_.readIfPresent(coeffDict());
-        A1_.readIfPresent(coeffDict());
-        A2_.readIfPresent(coeffDict());
-        Ctau1_.readIfPresent(coeffDict());
-        Ctau2_.readIfPresent(coeffDict());
-        Ctau3_.readIfPresent(coeffDict());
-        alphaKsi_.readIfPresent(coeffDict());
-        CmuWall_.readIfPresent(coeffDict());
-        kappa_.readIfPresent(coeffDict());
-        Am_.readIfPresent(coeffDict());
-        Aepsilon_.readIfPresent(coeffDict());
-        Amu_.readIfPresent(coeffDict());
-
-        return true;
-    }
-    else
-    {
-        return false;
-    }
-}
-
-
-void LienCubicKELowRe::correct()
-{
-    nonlinearEddyViscosity<incompressible::RASModel>::correct();
-
-    if (!turbulence_)
-    {
-        return;
-    }
-
-    tmp<volTensorField> tgradU = fvc::grad(U_);
-    const volTensorField& gradU = tgradU();
-
-    yStar_ = sqrt(k_)*y_/nu() + SMALL;
-    tmp<volScalarField> Rt = sqr(k_)/(nu()*epsilon_);
-
-    const volScalarField f2
-    (
-        scalar(1) - 0.3*exp(-sqr(Rt))
-    );
-
-    volScalarField G
-    (
-        GName(),
-        (nut_*twoSymm(gradU) - nonlinearStress_) && gradU
-    );
-
-    // Dissipation equation
-    tmp<fvScalarMatrix> epsEqn
-    (
-        fvm::ddt(epsilon_)
-      + fvm::div(phi_, epsilon_)
-      - fvm::laplacian(DepsilonEff(), epsilon_)
-      ==
-        C1_*G*epsilon_/k_
-        // E-term
-      + C2_*f2*pow(Cmu_, 0.75)*pow(k_, scalar(0.5))
-       /(kappa_*y_*(scalar(1) - exp(-Aepsilon_*yStar_)))
-       *exp(-Amu_*sqr(yStar_))*epsilon_
-      - fvm::Sp(C2_*f2*epsilon_/k_, epsilon_)
-    );
-
-    epsEqn().relax();
-    solve(epsEqn);
-    bound(epsilon_, epsilonMin_);
-
-
-    // Turbulent kinetic energy equation
-    tmp<fvScalarMatrix> kEqn
-    (
-        fvm::ddt(k_)
-      + fvm::div(phi_, k_)
-      - fvm::laplacian(DkEff(), k_)
-      ==
-        G
-      - fvm::Sp(epsilon_/k_, k_)
-    );
-
-    kEqn().relax();
-    solve(kEqn);
-    bound(k_, kMin_);
-
-
-    // Re-calculate viscosity
-
-    eta_ = k_/epsilon_*sqrt(2.0*magSqr(0.5*(gradU + gradU.T())));
-    ksi_ = k_/epsilon_*sqrt(2.0*magSqr(0.5*(gradU - gradU.T())));
-    Cmu_ = 2.0/(3.0*(A1_ + eta_ + alphaKsi_*ksi_));
-    fEta_ = A2_ + pow3(eta_);
-
-    C5viscosity_ =
-       -2.0*pow3(Cmu_)*pow4(k_)/pow3(epsilon_)
-       *(magSqr(gradU + gradU.T()) - magSqr(gradU - gradU.T()));
-
-    correctNut();
-    correctNonlinearStress(gradU);
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace RASModels
-} // End namespace incompressible
-} // End namespace Foam
-
-// ************************************************************************* //

src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/LienCubicKELowRe/LienCubicKELowRe.H

@@ -1,202 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
-     \\/     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 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::incompressible::RASModels::LienCubicKELowRe
-
-Group
-    grpIcoRASTurbulence
-
-Description
-    Lien cubic non-linear low-Reynolds k-epsilon turbulence models for
-    incompressible flows.
-
-    This turbulence model is described in:
-    \verbatim
-        Lien, F.S., Chen, W.L. & Leschziner, M.A. (1996).
-        Low-Reynolds-number eddy-viscosity modeling based on non-linear
-        stress-strain/vorticity relations.
-        Engineering Turbulence Modelling and Experiments 3, 91-100.
-    \endverbatim
-
-See Also
-    Foam::incompressible::RASModels::LienCubicKE
-
-SourceFiles
-    LienCubicKELowRe.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef LienCubicKELowRe_H
-#define LienCubicKELowRe_H
-
-#include "turbulentTransportModel.H"
-#include "nonlinearEddyViscosity.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-namespace incompressible
-{
-namespace RASModels
-{
-
-/*---------------------------------------------------------------------------*\
-                           Class LienCubicKELowRe Declaration
-\*---------------------------------------------------------------------------*/
-
-class LienCubicKELowRe
-:
-    public nonlinearEddyViscosity<incompressible::RASModel>
-{
-
-protected:
-
-    // Protected data
-
-        // Model coefficients
-
-            dimensionedScalar C1_;
-            dimensionedScalar C2_;
-            dimensionedScalar sigmak_;
-            dimensionedScalar sigmaEps_;
-            dimensionedScalar A1_;
-            dimensionedScalar A2_;
-            dimensionedScalar Ctau1_;
-            dimensionedScalar Ctau2_;
-            dimensionedScalar Ctau3_;
-            dimensionedScalar alphaKsi_;
-
-            dimensionedScalar CmuWall_;
-            dimensionedScalar kappa_;
-
-            dimensionedScalar Am_;
-            dimensionedScalar Aepsilon_;
-            dimensionedScalar Amu_;
-
-
-        // Fields
-
-            volScalarField k_;
-            volScalarField epsilon_;
-
-            //- Wall distance
-            //  Note: different to wall distance in parent RASModel
-            //  which is for near-wall cells only
-            const volScalarField& y_;
-
-            volScalarField eta_;
-            volScalarField ksi_;
-            volScalarField Cmu_;
-            volScalarField fEta_;
-            volScalarField C5viscosity_;
-
-            volScalarField yStar_;
-
-
-    // Protected Member Functions
-
-        tmp<volScalarField> fMu();
-
-        virtual void correctNut();
-        virtual void correctNonlinearStress(const volTensorField& gradU);
-
-
-public:
-
-    //- Runtime type information
-    TypeName("LienCubicKELowRe");
-
-    // Constructors
-
-        //- Construct from components
-        LienCubicKELowRe
-        (
-            const geometricOneField& alpha,
-            const geometricOneField& rho,
-            const volVectorField& U,
-            const surfaceScalarField& alphaRhoPhi,
-            const surfaceScalarField& phi,
-            const transportModel& transport,
-            const word& propertiesName = turbulenceModel::propertiesName,
-            const word& type = typeName
-        );
-
-
-    //- Destructor
-    virtual ~LienCubicKELowRe()
-    {}
-
-
-    // Member Functions
-
-        //- Read RASProperties dictionary
-        virtual bool read();
-
-        //- Return the effective diffusivity for k
-        tmp<volScalarField> DkEff() const
-        {
-            return tmp<volScalarField>
-            (
-                new volScalarField("DkEff", nut_/sigmak_ + nu())
-            );
-        }
-
-        //- Return the effective diffusivity for epsilon
-        tmp<volScalarField> DepsilonEff() const
-        {
-            return tmp<volScalarField>
-            (
-                new volScalarField("DepsilonEff", nut_/sigmaEps_ + nu())
-            );
-        }
-
-        //- Return the turbulence kinetic energy
-        virtual tmp<volScalarField> k() const
-        {
-            return k_;
-        }
-
-        //- Return the turbulence kinetic energy dissipation rate
-        virtual tmp<volScalarField> epsilon() const
-        {
-            return epsilon_;
-        }
-
-        //- Solve the turbulence equations and correct the turbulence viscosity
-        virtual void correct();
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace RASModels
-} // Edn namespace incompressible
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/LienLeschziner/LienLeschziner.C

@@ -23,9 +23,8 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "LienLeschzinerLowRe.H"
+#include "LienLeschziner.H"
 #include "wallDist.H"
-#include "wallFvPatch.H"
 #include "bound.H"
 #include "addToRunTimeSelectionTable.H"
 
@@ -40,20 +39,44 @@ namespace RASModels
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
-defineTypeNameAndDebug(LienLeschzinerLowRe, 0);
+defineTypeNameAndDebug(LienLeschziner, 0);
-addToRunTimeSelectionTable(RASModel, LienLeschzinerLowRe, dictionary);
+addToRunTimeSelectionTable(RASModel, LienLeschziner, dictionary);
 
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
-tmp<volScalarField> LienLeschzinerLowRe::fMu()
+tmp<volScalarField> LienLeschziner::fMu() const
 {
+    const volScalarField yStar(sqrt(k_)*y_/nu());
+
     return
-        (scalar(1) - exp(-Am_*yStar_))
+        (scalar(1) - exp(-Anu_*yStar))
-       /(scalar(1) - exp(-Aepsilon_*yStar_) + SMALL);
+       /((scalar(1) + SMALL) - exp(-Aeps_*yStar));
 }
 
 
-void LienLeschzinerLowRe::correctNut()
+tmp<volScalarField> LienLeschziner::f2() const
+{
+    tmp<volScalarField> Rt = sqr(k_)/(nu()*epsilon_);
+
+    return scalar(1) - 0.3*exp(-sqr(Rt));
+}
+
+
+tmp<volScalarField> LienLeschziner::E(const volScalarField& f2) const
+{
+    const volScalarField yStar(sqrt(k_)*y_/nu());
+    const volScalarField le
+    (
+        kappa_*y_*((scalar(1) + SMALL) - exp(-Aeps_*yStar))
+    );
+
+    return
+        (Ceps2_*pow(Cmu_, 0.75))
+       *(f2*sqrt(k_)*epsilon_/le)*exp(-AE_*sqr(yStar));
+}
+
+
+void LienLeschziner::correctNut()
 {
     nut_ = Cmu_*fMu()*sqr(k_)/epsilon_;
     nut_.correctBoundaryConditions();
@@ -62,7 +85,7 @@ void LienLeschzinerLowRe::correctNut()
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-LienLeschzinerLowRe::LienLeschzinerLowRe
+LienLeschziner::LienLeschziner
 (
     const geometricOneField& alpha,
     const geometricOneField& rho,
@@ -86,20 +109,20 @@ LienLeschzinerLowRe::LienLeschzinerLowRe
         propertiesName
     ),
 
-    C1_
+    Ceps1_
     (
         dimensioned<scalar>::lookupOrAddToDict
         (
-            "C1",
+            "Ceps1",
             coeffDict_,
             1.44
         )
     ),
-    C2_
+    Ceps2_
     (
         dimensioned<scalar>::lookupOrAddToDict
         (
-            "C2",
+            "Ceps2",
             coeffDict_,
             1.92
         )
@@ -140,29 +163,29 @@ LienLeschzinerLowRe::LienLeschzinerLowRe
             0.41
         )
     ),
-    Am_
+    Anu_
     (
         dimensioned<scalar>::lookupOrAddToDict
         (
-            "Am",
+            "Anu",
             coeffDict_,
             0.016
         )
     ),
-    Aepsilon_
+    Aeps_
     (
         dimensioned<scalar>::lookupOrAddToDict
         (
-            "Aepsilon",
+            "Aeps",
             coeffDict_,
             0.263
         )
     ),
-    Amu_
+    AE_
     (
         dimensioned<scalar>::lookupOrAddToDict
         (
-            "Amu",
+            "AE",
             coeffDict_,
             0.00222
         )
@@ -194,9 +217,7 @@ LienLeschzinerLowRe::LienLeschzinerLowRe
         mesh_
     ),
 
-    y_(wallDist::New(mesh_).y()),
+    y_(wallDist::New(mesh_).y())
-
-    yStar_(sqrt(k_)*y_/nu() + SMALL)
 {
     bound(k_, kMin_);
     bound(epsilon_, epsilonMin_);
@@ -211,19 +232,19 @@ LienLeschzinerLowRe::LienLeschzinerLowRe
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-bool LienLeschzinerLowRe::read()
+bool LienLeschziner::read()
 {
     if (eddyViscosity<incompressible::RASModel>::read())
     {
-        C1_.readIfPresent(coeffDict());
+        Ceps1_.readIfPresent(coeffDict());
-        C2_.readIfPresent(coeffDict());
+        Ceps2_.readIfPresent(coeffDict());
         sigmak_.readIfPresent(coeffDict());
         sigmaEps_.readIfPresent(coeffDict());
         Cmu_.readIfPresent(coeffDict());
         kappa_.readIfPresent(coeffDict());
-        Am_.readIfPresent(coeffDict());
+        Anu_.readIfPresent(coeffDict());
-        Aepsilon_.readIfPresent(coeffDict());
+        Aeps_.readIfPresent(coeffDict());
-        Amu_.readIfPresent(coeffDict());
+        AE_.readIfPresent(coeffDict());
 
         return true;
     }
@@ -234,23 +255,27 @@ bool LienLeschzinerLowRe::read()
 }
 
 
-void LienLeschzinerLowRe::correct()
+void LienLeschziner::correct()
 {
-    eddyViscosity<incompressible::RASModel>::correct();
-
     if (!turbulence_)
     {
         return;
     }
 
+    eddyViscosity<incompressible::RASModel>::correct();
+
     tmp<volTensorField> tgradU = fvc::grad(U_);
-    volScalarField G(GName(), nut_*(tgradU() && twoSymm(tgradU())));
+    volScalarField G
+    (
+        GName(),
+        nut_*(tgradU() && twoSymm(tgradU()))
+    );
     tgradU.clear();
 
-    scalar Cmu75 = pow(Cmu_.value(), 0.75);
+    // Update epsilon and G at the wall
-    yStar_ = sqrt(k_)*y_/nu() + SMALL;
+    epsilon_.boundaryField().updateCoeffs();
-    tmp<volScalarField> Rt = sqr(k_)/(nu()*epsilon_);
+
-    const volScalarField f2(scalar(1) - 0.3*exp(-sqr(Rt)));
+    const volScalarField f2(this->f2());
 
     // Dissipation equation
     tmp<fvScalarMatrix> epsEqn
@@ -258,18 +283,14 @@ void LienLeschzinerLowRe::correct()
         fvm::ddt(epsilon_)
       + fvm::div(phi_, epsilon_)
       - fvm::laplacian(DepsilonEff(), epsilon_)
-      ==
+     ==
-        C1_*G*epsilon_/k_
+        Ceps1_*G*epsilon_/k_
-
+      - fvm::Sp(Ceps2_*f2*epsilon_/k_, epsilon_)
-        // E-term
+      + E(f2)
-        + C2_*f2*Cmu75*sqrt(k_)
-         /(kappa_*y_*(scalar(1) - exp(-Aepsilon_*yStar_)))
-         *exp(-Amu_*sqr(yStar_))*epsilon_
-
-      - fvm::Sp(C2_*f2*epsilon_/k_, epsilon_)
     );
 
     epsEqn().relax();
+    epsEqn().boundaryManipulate(epsilon_.boundaryField());
     solve(epsEqn);
     bound(epsilon_, epsilonMin_);
 
@@ -280,7 +301,7 @@ void LienLeschzinerLowRe::correct()
         fvm::ddt(k_)
       + fvm::div(phi_, k_)
       - fvm::laplacian(DkEff(), k_)
-      ==
+     ==
         G
       - fvm::Sp(epsilon_/k_, k_)
     );

src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/LienLeschziner/LienLeschziner.H

@@ -22,13 +22,13 @@ License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 Class
-    Foam::incompressible::RASModels::LienLeschzinerLowRe
+    Foam::incompressible::RASModels::LienLeschziner
 
 Group
     grpIcoRASTurbulence
 
 Description
-    Lien and Leschziner low-Reynolds k-epsilon turbulence model for
+    Lien and Leschziner low-Reynolds number k-epsilon turbulence model for
     incompressible flows.
 
     This turbulence model is described in:
@@ -40,13 +40,22 @@ Description
         Journal of fluids engineering, 115(4), 717-725.
     \endverbatim
 
+    Implemented according to the specification in:
+    <a href=
+    "http://personalpages.manchester.ac.uk/staff/david.d.apsley/specturb.pdf"
+    >Apsley: Turbulence Models 2002</a>
+
+    In addition to the low-Reynolds number damping functions support for
+    wall-functions is also included to allow for low- and high-Reynolds number
+    operation.
+
 SourceFiles
-    LienLeschzinerLowRe.C
+    LienLeschziner.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef LienLeschzinerLowRe_H
+#ifndef LienLeschziner_H
-#define LienLeschzinerLowRe_H
+#define LienLeschziner_H
 
 #include "turbulentTransportModel.H"
 #include "eddyViscosity.H"
@@ -61,10 +70,10 @@ namespace RASModels
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class LienLeschzinerLowRe Declaration
+                           Class LienLeschziner Declaration
 \*---------------------------------------------------------------------------*/
 
-class LienLeschzinerLowRe
+class LienLeschziner
 :
     public eddyViscosity<incompressible::RASModel>
 {
@@ -75,16 +84,17 @@ protected:
 
         // Model coefficients
 
-            dimensionedScalar C1_;
+            dimensionedScalar Ceps1_;
-            dimensionedScalar C2_;
+            dimensionedScalar Ceps2_;
             dimensionedScalar sigmak_;
             dimensionedScalar sigmaEps_;
+
             dimensionedScalar Cmu_;
             dimensionedScalar kappa_;
 
-            dimensionedScalar Am_;
+            dimensionedScalar Anu_;
-            dimensionedScalar Aepsilon_;
+            dimensionedScalar Aeps_;
-            dimensionedScalar Amu_;
+            dimensionedScalar AE_;
 
 
         // Fields
@@ -97,24 +107,24 @@ protected:
             //  which is for near-wall cells only
             const volScalarField& y_;
 
-            volScalarField yStar_;
-
 
     // Protected Member Functions
 
-        tmp<volScalarField> fMu();
+        tmp<volScalarField> fMu() const;
+        tmp<volScalarField> f2() const;
+        tmp<volScalarField> E(const volScalarField& f2) const;
 
         virtual void correctNut();
 
 
 public:
 
-    TypeName("LienLeschzinerLowRe");
+    TypeName("LienLeschziner");
 
     // Constructors
 
         //- Construct from components
-        LienLeschzinerLowRe
+        LienLeschziner
         (
             const geometricOneField& alpha,
             const geometricOneField& rho,
@@ -128,7 +138,7 @@ public:
 
 
     //- Destructor
-    virtual ~LienLeschzinerLowRe()
+    virtual ~LienLeschziner()
     {}
 
 

src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/ShihQuadraticKE/ShihQuadraticKE.C

@@ -56,20 +56,20 @@ void ShihQuadraticKE::correctNonlinearStress(const volTensorField& gradU)
     volSymmTensorField S(symm(gradU));
     volTensorField W(skew(gradU));
 
-    volScalarField sBar((k_/epsilon_)*sqrt(2.0)*mag(symm(gradU)));
+    volScalarField sBar((k_/epsilon_)*sqrt(2.0)*mag(S));
-    volScalarField wBar((k_/epsilon_)*sqrt(2.0)*mag(skew(gradU)));
+    volScalarField wBar((k_/epsilon_)*sqrt(2.0)*mag(W));
 
-    volScalarField Cmu(2.0/(3.0*(Cmu1_ + sBar + Cmu2_*wBar)));
+    volScalarField Cmu((2.0/3.0)/(Cmu1_ + sBar + Cmu2_*wBar));
 
     nut_ = Cmu*sqr(k_)/epsilon_;
     nut_.correctBoundaryConditions();
 
     nonlinearStress_ =
-        pow3(k_)/((A1_ + pow3(sBar))*sqr(epsilon_))
+        k_*sqr(k_/epsilon_)/(Cbeta_ + pow3(sBar))
        *(
-           beta1_*dev(innerSqr(S))
+           Cbeta1_*dev(innerSqr(S))
-         + beta2_*twoSymm(S&W)
+         + Cbeta2_*twoSymm(S&W)
-         + beta3_*dev(symm(W&W))
+         + Cbeta3_*dev(symm(W&W))
         );
 }
 
@@ -100,20 +100,20 @@ ShihQuadraticKE::ShihQuadraticKE
         propertiesName
     ),
 
-    C1_
+    Ceps1_
     (
         dimensioned<scalar>::lookupOrAddToDict
         (
-            "C1",
+            "Ceps1",
             coeffDict_,
             1.44
         )
     ),
-    C2_
+    Ceps2_
     (
         dimensioned<scalar>::lookupOrAddToDict
         (
-            "C2",
+            "Ceps2",
             coeffDict_,
             1.92
         )
@@ -140,7 +140,7 @@ ShihQuadraticKE::ShihQuadraticKE
     (
         dimensioned<scalar>::lookupOrAddToDict
         (
-            "A1",
+            "Cmu1",
             coeffDict_,
             1.25
         )
@@ -149,43 +149,43 @@ ShihQuadraticKE::ShihQuadraticKE
     (
         dimensioned<scalar>::lookupOrAddToDict
         (
-            "alphaKsi",
+            "Cmu2",
             coeffDict_,
             0.9
         )
     ),
-    A1_
+    Cbeta_
     (
         dimensioned<scalar>::lookupOrAddToDict
         (
-            "A2",
+            "Cbeta",
             coeffDict_,
             1000.0
         )
     ),
-    beta1_
+    Cbeta1_
     (
         dimensioned<scalar>::lookupOrAddToDict
         (
-            "beta1",
+            "Cbeta1",
             coeffDict_,
             3.0
         )
     ),
-    beta2_
+    Cbeta2_
     (
         dimensioned<scalar>::lookupOrAddToDict
         (
-            "beta2",
+            "Cbeta2",
             coeffDict_,
             15.0
         )
     ),
-    beta3_
+    Cbeta3_
     (
         dimensioned<scalar>::lookupOrAddToDict
         (
-            "beta3",
+            "Cbeta3",
             coeffDict_,
             -19.0
         )
@@ -234,16 +234,16 @@ bool ShihQuadraticKE::read()
 {
     if (nonlinearEddyViscosity<incompressible::RASModel>::read())
     {
-        C1_.readIfPresent(coeffDict());
+        Ceps1_.readIfPresent(coeffDict());
-        C2_.readIfPresent(coeffDict());
+        Ceps2_.readIfPresent(coeffDict());
         sigmak_.readIfPresent(coeffDict());
         sigmaEps_.readIfPresent(coeffDict());
         Cmu1_.readIfPresent(coeffDict());
         Cmu2_.readIfPresent(coeffDict());
-        A1_.readIfPresent(coeffDict());
+        Cbeta_.readIfPresent(coeffDict());
-        beta1_.readIfPresent(coeffDict());
+        Cbeta1_.readIfPresent(coeffDict());
-        beta2_.readIfPresent(coeffDict());
+        Cbeta2_.readIfPresent(coeffDict());
-        beta3_.readIfPresent(coeffDict());
+        Cbeta3_.readIfPresent(coeffDict());
 
         return true;
     }
@@ -256,13 +256,13 @@ bool ShihQuadraticKE::read()
 
 void ShihQuadraticKE::correct()
 {
-    nonlinearEddyViscosity<incompressible::RASModel>::correct();
-
     if (!turbulence_)
     {
         return;
     }
 
+    nonlinearEddyViscosity<incompressible::RASModel>::correct();
+
     tmp<volTensorField> tgradU = fvc::grad(U_);
     const volTensorField& gradU = tgradU();
 
@@ -283,8 +283,8 @@ void ShihQuadraticKE::correct()
       + fvm::div(phi_, epsilon_)
       - fvm::laplacian(DepsilonEff(), epsilon_)
       ==
-        C1_*G*epsilon_/k_
+        Ceps1_*G*epsilon_/k_
-      - fvm::Sp(C2_*epsilon_/k_, epsilon_)
+      - fvm::Sp(Ceps2_*epsilon_/k_, epsilon_)
     );
 
     epsEqn().relax();

src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/ShihQuadraticKE/ShihQuadraticKE.H

@@ -38,6 +38,11 @@ Description
         NASA technical memorandum 105993.
     \endverbatim
 
+    Implemented according to the specification in:
+    <a href=
+    "http://personalpages.manchester.ac.uk/staff/david.d.apsley/specturb.pdf"
+    >Apsley: Turbulence Models 2002</a>
+
 SourceFiles
     ShihQuadraticKE.C
 
@@ -73,16 +78,16 @@ protected:
 
         // Model coefficients
 
-            dimensionedScalar C1_;
+            dimensionedScalar Ceps1_;
-            dimensionedScalar C2_;
+            dimensionedScalar Ceps2_;
             dimensionedScalar sigmak_;
             dimensionedScalar sigmaEps_;
             dimensionedScalar Cmu1_;
             dimensionedScalar Cmu2_;
-            dimensionedScalar A1_;
+            dimensionedScalar Cbeta_;
-            dimensionedScalar beta1_;
+            dimensionedScalar Cbeta1_;
-            dimensionedScalar beta2_;
+            dimensionedScalar Cbeta2_;
-            dimensionedScalar beta3_;
+            dimensionedScalar Cbeta3_;
 
 
         // Fields

src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/qZeta/qZeta.C

@@ -239,19 +239,16 @@ bool qZeta::read()
 
 void qZeta::correct()
 {
-    eddyViscosity<incompressible::RASModel>::correct();
-
     if (!turbulence_)
     {
         return;
     }
 
-    tmp<volScalarField> S2 = 2*magSqr(symm(fvc::grad(U_)));
+    eddyViscosity<incompressible::RASModel>::correct();
 
-    volScalarField G(GName(), nut_/(2.0*q_)*S2);
+    volScalarField G(GName(), nut_/(2.0*q_)*2*magSqr(symm(fvc::grad(U_))));
     const volScalarField E(nu()*nut_/q_*fvc::magSqrGradGrad(U_));
 
-
     // Zeta equation
     tmp<fvScalarMatrix> zetaEqn
     (