TurbulenceModels: Added DeardorffDiffStress LES model

Now using the Daly-Harlow generalized gradient diffusion model

src/TurbulenceModels/incompressible/turbulentTransportModels/turbulentTransportModels.C

@@ -82,5 +82,8 @@ makeLESModel(SpalartAllmarasDDES);
 #include "SpalartAllmarasIDDES.H"
 makeLESModel(SpalartAllmarasIDDES);
 
+#include "DeardorffDiffStress.H"
+makeLESModel(DeardorffDiffStress);
+
 
 // ************************************************************************* //

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

@@ -0,0 +1,215 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "DeardorffDiffStress.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace LESModels
+{
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+template<class BasicTurbulenceModel>
+void DeardorffDiffStress<BasicTurbulenceModel>::correctNut()
+{
+    this->nut_ = Ck_*sqrt(this->k())*this->delta();
+    this->nut_.correctBoundaryConditions();
+}
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+template<class BasicTurbulenceModel>
+DeardorffDiffStress<BasicTurbulenceModel>::DeardorffDiffStress
+(
+    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<LESModel<BasicTurbulenceModel> >
+    (
+        type,
+        alpha,
+        rho,
+        U,
+        alphaRhoPhi,
+        phi,
+        transport,
+        propertiesName
+    ),
+
+    Ck_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "Ck",
+            this->coeffDict_,
+            0.094
+        )
+    ),
+    Cm_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "Cm",
+            this->coeffDict_,
+            4.13
+        )
+    ),
+    Ce_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "Ce",
+            this->coeffDict_,
+            1.048
+        )
+    ),
+    Cs_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "Cs",
+            this->coeffDict_,
+            0.22
+        )
+    )
+{
+    if (type == typeName)
+    {
+        this->boundNormalStress(this->R_);
+        correctNut();
+        this->printCoeffs(type);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+template<class BasicTurbulenceModel>
+bool DeardorffDiffStress<BasicTurbulenceModel>::read()
+{
+    if (ReynoldsStress<LESModel<BasicTurbulenceModel> >::read())
+    {
+        Ck_.readIfPresent(this->coeffDict());
+        Cm_.readIfPresent(this->coeffDict());
+        Ce_.readIfPresent(this->coeffDict());
+        Cs_.readIfPresent(this->coeffDict());
+
+        return true;
+    }
+    else
+    {
+        return false;
+    }
+}
+
+
+template<class BasicTurbulenceModel>
+tmp<volScalarField> DeardorffDiffStress<BasicTurbulenceModel>::epsilon() const
+{
+    volScalarField k(this->k());
+
+    return tmp<volScalarField>
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                IOobject::groupName("epsilon", this->U_.group()),
+                this->runTime_.timeName(),
+                this->mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            this->Ce_*k*sqrt(k)/this->delta()
+        )
+    );
+}
+
+
+template<class BasicTurbulenceModel>
+void DeardorffDiffStress<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<LESModel<BasicTurbulenceModel> >::correct();
+
+    tmp<volTensorField> tgradU(fvc::grad(U));
+    const volTensorField& gradU = tgradU();
+
+    volSymmTensorField D(symm(gradU));
+
+    volSymmTensorField P(-twoSymm(R & gradU));
+
+    volScalarField k(this->k());
+
+    tmp<fvSymmTensorMatrix> BEqn
+    (
+        fvm::ddt(alpha, rho, R)
+      + fvm::div(alphaRhoPhi, R)
+      - fvm::laplacian(Cs_*(k/this->epsilon())*R, R)
+      + fvm::Sp(Cm_*alpha*rho*sqrt(k)/this->delta(), R)
+     ==
+        alpha*rho*P
+      + (4.0/5.0)*alpha*rho*k*D // Deardorff
+    //- 0.6*alpha*rho*dev(P) // LRR
+      - ((2.0/3.0)*(1.0 - Cm_/this->Ce_)*I)*(alpha*rho*this->epsilon())
+    );
+
+    BEqn().relax();
+    BEqn().solve();
+
+    this->boundNormalStress(this->R_);
+    correctNut();
+    this->correctWallShearStress(this->R_);
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace LESModels
+} // End namespace Foam
+
+// ************************************************************************* //

src/TurbulenceModels/turbulenceModels/LES/DeardorffDiffStress/DeardorffDiffStress.H

@@ -0,0 +1,166 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     |
+    \\  /    A nd           | Copyright (C) 2011-2012 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::LESModels::DeardorffDiffStress
+
+Group
+    grpLESTurbulence
+
+Description
+    Differential SGS Stress Equation Model for incompressible and
+    compressible flows
+
+    Reference:
+    \verbatim
+        Deardorff, J. W. (1973).
+        The use of subgrid transport equations in a three-dimensional model
+        of atmospheric turbulence.
+        Journal of Fluids Engineering, 95(3), 429-438.
+    \endverbatim
+
+    This SGS model uses a full balance equation for the SGS stress tensor to
+    simulate the behaviour of B.
+
+SourceFiles
+    DeardorffDiffStress.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef DeardorffDiffStress_H
+#define DeardorffDiffStress_H
+
+#include "LESModel.H"
+#include "ReynoldsStress.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+namespace LESModels
+{
+
+/*---------------------------------------------------------------------------*\
+                    Class DeardorffDiffStress Declaration
+\*---------------------------------------------------------------------------*/
+
+template<class BasicTurbulenceModel>
+class DeardorffDiffStress
+:
+    public ReynoldsStress<LESModel<BasicTurbulenceModel> >
+{
+    // Private Member Functions
+
+        // Disallow default bitwise copy construct and assignment
+        DeardorffDiffStress(const DeardorffDiffStress&);
+        DeardorffDiffStress& operator=(const DeardorffDiffStress&);
+
+
+protected:
+
+    // Protected data
+
+        // Model constants
+
+            dimensionedScalar Ck_;
+            dimensionedScalar Cm_;
+            dimensionedScalar Ce_;
+            dimensionedScalar Cs_;
+
+
+    // 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("DeardorffDiffStress");
+
+
+    // Constructors
+
+        //- Constructor from components
+        DeardorffDiffStress
+        (
+            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 ~DeardorffDiffStress()
+    {}
+
+
+    // Member Functions
+
+        //- Read model coefficients if they have changed
+        virtual bool read();
+
+        //- Return the turbulence kinetic energy dissipation rate
+        virtual tmp<volScalarField> epsilon() const;
+
+        //- Return the effective diffusivity for B
+        tmp<volScalarField> DBEff() const
+        {
+            return tmp<volScalarField>
+            (
+                new volScalarField("DBEff", this->nut_ + this->nu())
+            );
+        }
+
+        //- Correct sub-grid stress, eddy-Viscosity and related properties
+        virtual void correct();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace LESModels
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#ifdef NoRepository
+#   include "DeardorffDiffStress.C"
+#endif
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/TurbulenceModels/turbulenceModels/LES/WALE/WALE.H

@@ -100,7 +100,7 @@ protected:
         //  calculated from the given velocity gradient
         tmp<volScalarField> k(const volTensorField& gradU) const;
 
-        //- Update the SGS eddy viscosity
+        //- Update the SGS eddy-viscosity
         virtual void correctNut();
 
 

src/TurbulenceModels/turbulenceModels/LES/kEqn/kEqn.H

@@ -50,7 +50,7 @@ Description
 
     The default model coefficients correspond to the following:
     \verbatim
-        NicenoKEqnCoeffs
+        kEqnCoeffs
         {
             Ck                  0.094;
             Ce                  1.048;
@@ -165,7 +165,7 @@ public:
             );
         }
 
-        //- Correct Eddy-Viscosity and related properties
+        //- Correct eddy-Viscosity and related properties
         virtual void correct();
 };
 

src/TurbulenceModels/turbulenceModels/ReynoldsStress/ReynoldsStress.C

@@ -26,6 +26,82 @@ License
 #include "ReynoldsStress.H"
 #include "fvc.H"
 #include "fvm.H"
+#include "wallFvPatch.H"
+
+// * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
+
+template<class BasicTurbulenceModel>
+void Foam::ReynoldsStress<BasicTurbulenceModel>::boundNormalStress
+(
+    volSymmTensorField& R
+) const
+{
+    scalar kMin = this->kMin_.value();
+
+    R.max
+    (
+        dimensionedSymmTensor
+        (
+            "zero",
+            R.dimensions(),
+            symmTensor
+            (
+                kMin, -GREAT, -GREAT,
+                kMin, -GREAT,
+                kMin
+            )
+        )
+    );
+}
+
+
+template<class BasicTurbulenceModel>
+void Foam::ReynoldsStress<BasicTurbulenceModel>::correctWallShearStress
+(
+    volSymmTensorField& R
+) const
+{
+    const fvPatchList& patches = this->mesh_.boundary();
+
+    forAll(patches, patchi)
+    {
+        const fvPatch& curPatch = patches[patchi];
+
+        if (isA<wallFvPatch>(curPatch))
+        {
+            symmTensorField& Rw = R.boundaryField()[patchi];
+
+            const scalarField& nutw = this->nut_.boundaryField()[patchi];
+
+            const vectorField snGradU
+            (
+                this->U_.boundaryField()[patchi].snGrad()
+            );
+
+            const vectorField& faceAreas
+                = this->mesh_.Sf().boundaryField()[patchi];
+
+            const scalarField& magFaceAreas
+                = this->mesh_.magSf().boundaryField()[patchi];
+
+            forAll(curPatch, facei)
+            {
+                // Calculate near-wall velocity gradient
+                tensor gradUw
+                    = (faceAreas[facei]/magFaceAreas[facei])*snGradU[facei];
+
+                // Calculate near-wall shear-stress tensor
+                tensor tauw = -nutw[facei]*2*dev(symm(gradUw));
+
+                // Reset the shear components of the stress tensor
+                Rw[facei].xy() = tauw.xy();
+                Rw[facei].xz() = tauw.xz();
+                Rw[facei].yz() = tauw.yz();
+            }
+        }
+    }
+}
+
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -113,7 +189,7 @@ template<class BasicTurbulenceModel>
 Foam::tmp<Foam::volScalarField>
 Foam::ReynoldsStress<BasicTurbulenceModel>::k() const
 {
-    tmp<Foam::volScalarField> tk(tr(R_));
+    tmp<Foam::volScalarField> tk(0.5*tr(R_));
     tk().rename("k");
     return tk;
 }
@@ -168,7 +244,7 @@ Foam::ReynoldsStress<BasicTurbulenceModel>::divDevRhoReff
                 "laplacian(nuEff,U)"
             )
           - fvm::laplacian(this->alpha_*this->rho_*this->nuEff(), U)
-          - fvc::div(this->alpha_*this->rho_*nu()*dev2(T(fvc::grad(U))))
+          - fvc::div(this->alpha_*this->rho_*this->nu()*dev2(T(fvc::grad(U))))
         );
     }
     else
@@ -183,7 +259,7 @@ Foam::ReynoldsStress<BasicTurbulenceModel>::divDevRhoReff
                 "laplacian(nuEff,U)"
             )
           - fvm::laplacian(this->alpha_*this->rho_*this->nuEff(), U)
-          - fvc::div(this->alpha_*this->rho_*nu()*dev2(T(fvc::grad(U))))
+          - fvc::div(this->alpha_*this->rho_*this->nu()*dev2(T(fvc::grad(U))))
         );
     }
 

src/TurbulenceModels/turbulenceModels/ReynoldsStress/ReynoldsStress.H

@@ -69,6 +69,10 @@ protected:
 
     // Protected Member Functions
 
+        void boundNormalStress(volSymmTensorField& R) const;
+        void correctWallShearStress(volSymmTensorField& R) const;
+
+        //- Update the eddy-viscosity
         virtual void correctNut() = 0;