TurbulenceModels: Added LRR model with Daly-Harlow generalized gradient diffusion

src/TurbulenceModels/compressible/turbulentFluidThermoModels/turbulentFluidThermoModels.C

@@ -69,6 +69,10 @@ makeRASModel(LaunderSharmaKE);
 #include "kOmegaSST.H"
 makeRASModel(kOmegaSST);
 
+#include "LRR.H"
+makeRASModel(LRR);
+
+
 #include "Smagorinsky.H"
 makeLESModel(Smagorinsky);
 
@@ -87,5 +91,8 @@ makeLESModel(SpalartAllmarasDDES);
 #include "SpalartAllmarasIDDES.H"
 makeLESModel(SpalartAllmarasIDDES);
 
+#include "DeardorffDiffStress.H"
+makeLESModel(DeardorffDiffStress);
+
 
 // ************************************************************************* //

src/TurbulenceModels/incompressible/turbulentTransportModels/turbulentTransportModels.C

@@ -64,6 +64,10 @@ makeRASModel(LaunderSharmaKE);
 #include "kOmegaSST.H"
 makeRASModel(kOmegaSST);
 
+#include "LRR.H"
+makeRASModel(LRR);
+
+
 #include "Smagorinsky.H"
 makeLESModel(Smagorinsky);
 

src/TurbulenceModels/turbulenceModels/LES/DeardorffDiffStress/DeardorffDiffStress.C

@@ -185,7 +185,7 @@ void DeardorffDiffStress<BasicTurbulenceModel>::correct()
 
     volScalarField k(this->k());
 
-    tmp<fvSymmTensorMatrix> BEqn
+    tmp<fvSymmTensorMatrix> REqn
     (
         fvm::ddt(alpha, rho, R)
       + fvm::div(alphaRhoPhi, R)
@@ -197,10 +197,10 @@ void DeardorffDiffStress<BasicTurbulenceModel>::correct()
       - ((2.0/3.0)*(1.0 - Cm_/this->Ce_)*I)*(alpha*rho*this->epsilon())
     );
 
-    BEqn().relax();
-    BEqn().solve();
+    REqn().relax();
+    REqn().solve();
 
-    this->boundNormalStress(this->R_);
+    this->boundNormalStress(R);
     correctNut();
 }
 

src/TurbulenceModels/turbulenceModels/RAS/LRR/LRR.C

@@ -0,0 +1,295 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "LRR.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace RASModels
+{
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+template<class BasicTurbulenceModel>
+void LRR<BasicTurbulenceModel>::correctNut()
+{
+    this->nut_ = this->Cmu_*sqr(k_)/epsilon_;
+    this->nut_.correctBoundaryConditions();
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class BasicTurbulenceModel>
+LRR<BasicTurbulenceModel>::LRR
+(
+    const alphaField& alpha,
+    const rhoField& rho,
+    const volVectorField& U,
+    const surfaceScalarField& alphaRhoPhi,
+    const surfaceScalarField& phi,
+    const transportModel& transport,
+    const word& propertiesName,
+    const word& type
+)
+:
+    ReynoldsStress<RASModel<BasicTurbulenceModel> >
+    (
+        type,
+        alpha,
+        rho,
+        U,
+        alphaRhoPhi,
+        phi,
+        transport,
+        propertiesName
+    ),
+
+    Cmu_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "Cmu",
+            this->coeffDict_,
+            0.09
+        )
+    ),
+    Clrr1_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "Clrr1",
+            this->coeffDict_,
+            1.8
+        )
+    ),
+    Clrr2_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "Clrr2",
+            this->coeffDict_,
+            0.6
+        )
+    ),
+    C1_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "C1",
+            this->coeffDict_,
+            1.44
+        )
+    ),
+    C2_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "C2",
+            this->coeffDict_,
+            1.92
+        )
+    ),
+    Cs_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "Cs",
+            this->coeffDict_,
+            0.25
+        )
+    ),
+    Ceps_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "Ceps",
+            this->coeffDict_,
+            0.15
+        )
+    ),
+
+    k_
+    (
+        IOobject
+        (
+            "k",
+            this->runTime_.timeName(),
+            this->mesh_,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        0.5*tr(this->R_)
+    ),
+    epsilon_
+    (
+        IOobject
+        (
+            "epsilon",
+            this->runTime_.timeName(),
+            this->mesh_,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        this->mesh_
+    )
+{
+    if (type == typeName)
+    {
+        this->boundNormalStress(this->R_);
+        bound(epsilon_, this->epsilonMin_);
+        k_ = 0.5*tr(this->R_);
+        correctNut();
+        this->printCoeffs(type);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class BasicTurbulenceModel>
+bool LRR<BasicTurbulenceModel>::read()
+{
+    if (ReynoldsStress<RASModel<BasicTurbulenceModel> >::read())
+    {
+        Cmu_.readIfPresent(this->coeffDict());
+        Clrr1_.readIfPresent(this->coeffDict());
+        Clrr2_.readIfPresent(this->coeffDict());
+        C1_.readIfPresent(this->coeffDict());
+        C2_.readIfPresent(this->coeffDict());
+        Cs_.readIfPresent(this->coeffDict());
+        Ceps_.readIfPresent(this->coeffDict());
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+template<class BasicTurbulenceModel>
+void LRR<BasicTurbulenceModel>::correct()
+{
+    if (!this->turbulence_)
+    {
+        return;
+    }
+
+    // Local references
+    const alphaField& alpha = this->alpha_;
+    const rhoField& rho = this->rho_;
+    const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_;
+    const volVectorField& U = this->U_;
+    volSymmTensorField& R = this->R_;
+
+    ReynoldsStress<RASModel<BasicTurbulenceModel> >::correct();
+
+    tmp<volTensorField> tgradU(fvc::grad(U));
+    const volTensorField& gradU = tgradU();
+
+    volSymmTensorField P(-twoSymm(R & gradU));
+    volScalarField G(this->GName(), 0.5*mag(tr(P)));
+
+    // Update epsilon and G at the wall
+    epsilon_.boundaryField().updateCoeffs();
+
+    // Dissipation equation
+    tmp<fvScalarMatrix> epsEqn
+    (
+        fvm::ddt(alpha, rho, epsilon_)
+      + fvm::div(alphaRhoPhi, epsilon_)
+      - fvm::laplacian(Ceps_*alpha*rho*(k_/epsilon_)*R, epsilon_)
+     ==
+        C1_*alpha*rho*G*epsilon_/k_
+      - fvm::Sp(C2_*alpha*rho*epsilon_/k_, epsilon_)
+    );
+
+    epsEqn().relax();
+
+    epsEqn().boundaryManipulate(epsilon_.boundaryField());
+
+    solve(epsEqn);
+    bound(epsilon_, this->epsilonMin_);
+
+
+    // Reynolds stress equation
+
+    const fvPatchList& patches = this->mesh_.boundary();
+
+    forAll(patches, patchi)
+    {
+        const fvPatch& curPatch = patches[patchi];
+
+        if (isA<wallFvPatch>(curPatch))
+        {
+            forAll(curPatch, facei)
+            {
+                label faceCelli = curPatch.faceCells()[facei];
+                P[faceCelli] *= min
+                (
+                    G[faceCelli]/(0.5*mag(tr(P[faceCelli])) + SMALL),
+                    1.0
+                );
+            }
+        }
+    }
+
+
+    tmp<fvSymmTensorMatrix> REqn
+    (
+        fvm::ddt(alpha, rho, R)
+      + fvm::div(alphaRhoPhi, R)
+      - fvm::laplacian(Cs_*alpha*rho*(k_/epsilon_)*R, R)
+      + fvm::Sp(Clrr1_*alpha*rho*epsilon_/k_, R)
+      ==
+        alpha*rho*P
+      - (2.0/3.0*(1 - Clrr1_)*I)*alpha*rho*epsilon_
+      - Clrr2_*alpha*rho*dev(P)
+    );
+
+    REqn().relax();
+    solve(REqn);
+
+    this->boundNormalStress(R);
+
+    k_ = 0.5*tr(R);
+
+    correctNut();
+
+    // Correct wall shear-stresses when applying wall-functions
+    this->correctWallShearStress(R);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace RASModels
+} // End namespace Foam
+
+// ************************************************************************* //

src/TurbulenceModels/turbulenceModels/RAS/LRR/LRR.H

@@ -0,0 +1,197 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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::RASModels::LRR
+
+Group
+    grpIcoRASTurbulence
+
+Description
+    Launder, Reece and Rodi Reynolds-stress turbulence model for
+    incompressible and compressible flows.
+
+    Reference:
+    \verbatim
+        Launder, B. E., Reece, G. J., & Rodi, W. (1975).
+        Progress in the development of a Reynolds-stress turbulence closure.
+        Journal of fluid mechanics, 68(03), 537-566.
+    \endverbatim
+
+    Including the recommended generalized gradient diffusion model of
+    Daly and Harlow:
+    \verbatim
+        Daly, B. J., & Harlow, F. H. (1970).
+        Transport equations in turbulence.
+        Physics of Fluids (1958-1988), 13(11), 2634-2649.
+    \endverbatim
+
+    The default model coefficients are:
+    \verbatim
+        LRRCoeffs
+        {
+            Cmu         0.09;
+            Clrr1       1.8;
+            Clrr2       0.6;
+            C1          1.44;
+            C2          1.92;
+            Cs          0.25;
+            Ceps        0.15;
+            couplingFactor  0.0;
+        }
+    \endverbatim
+
+SourceFiles
+    LRR.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef LRR_H
+#define LRR_H
+
+#include "RASModel.H"
+#include "ReynoldsStress.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace RASModels
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class LRR Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class BasicTurbulenceModel>
+class LRR
+:
+    public ReynoldsStress<RASModel<BasicTurbulenceModel> >
+{
+    // Private Member Functions
+
+        // Disallow default bitwise copy construct and assignment
+        LRR(const LRR&);
+        LRR& operator=(const LRR&);
+
+
+protected:
+
+    // Protected data
+
+        // Model coefficients
+
+            dimensionedScalar Cmu_;
+
+            dimensionedScalar Clrr1_;
+            dimensionedScalar Clrr2_;
+
+            dimensionedScalar C1_;
+            dimensionedScalar C2_;
+            dimensionedScalar Cs_;
+            dimensionedScalar Ceps_;
+
+
+        // Fields
+
+            volScalarField k_;
+            volScalarField epsilon_;
+
+
+    // Protected Member Functions
+
+        //- Update the eddy-viscosity
+        virtual void correctNut();
+
+
+public:
+
+    typedef typename BasicTurbulenceModel::alphaField alphaField;
+    typedef typename BasicTurbulenceModel::rhoField rhoField;
+    typedef typename BasicTurbulenceModel::transportModel transportModel;
+
+
+    //- Runtime type information
+    TypeName("LRR");
+
+
+    // Constructors
+
+        //- Construct from components
+        LRR
+        (
+            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 ~LRR()
+    {}
+
+
+    // Member Functions
+
+        //- Read model coefficients if they have changed
+        virtual bool read();
+
+        //- 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 eddy-Viscosity and
+        //  related properties
+        virtual void correct();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace RASModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "LRR.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

tutorials/incompressible/boundaryFoam/boundaryWallFunctions/system/controlDict

@@ -23,13 +23,13 @@ startTime       0;
 
 stopAt          endTime;
 
-endTime         1000;
+endTime         3000;
 
 deltaT          1;
 
 writeControl    timeStep;
 
-writeInterval   100;
+writeInterval   200;
 
 purgeWrite      0;