ShihQuadraticKE: Rewritten and tested on the boundaryWallFunctions and pitzDaily cases

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

@@ -47,25 +47,29 @@ addToRunTimeSelectionTable(RASModel, ShihQuadraticKE, dictionary);
 
 void ShihQuadraticKE::correctNut()
 {
-    nut_ = Cmu_*sqr(k_)/epsilon_;
+    correctNonlinearStress(fvc::grad(U_));
-    #include "wallNonlinearViscosityI.H"
 }
 
 
 void ShihQuadraticKE::correctNonlinearStress(const volTensorField& gradU)
 {
-    nonlinearStress_ = symm
+    volSymmTensorField S(symm(gradU));
-    (
+    volTensorField W(skew(gradU));
-        pow(k_, 3.0)/sqr(epsilon_)
+
+    volScalarField sBar((k_/epsilon_)*sqrt(2.0)*mag(symm(gradU)));
+    volScalarField wBar((k_/epsilon_)*sqrt(2.0)*mag(skew(gradU)));
+
+    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_))
        *(
-            Ctau1_/fEta_
+           beta1_*dev(symm(S&S))
-           *(
+         + beta2_*symm((W&S) - (S&W))
-                (gradU & gradU)
+         + beta3_*dev(symm(W&W))
-              + (gradU & gradU)().T()
-            )
-          + Ctau2_/fEta_*(gradU & T(gradU))
-          + Ctau3_/fEta_*(T(gradU) & gradU)
-        )
         );
 }
 
@@ -132,7 +136,7 @@ ShihQuadraticKE::ShihQuadraticKE
             1.3
         )
     ),
-    A1_
+    Cmu1_
     (
         dimensioned<scalar>::lookupOrAddToDict
         (
@@ -141,43 +145,7 @@ ShihQuadraticKE::ShihQuadraticKE
             1.25
         )
     ),
-    A2_
+    Cmu2_
-    (
-        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
         (
@@ -186,23 +154,40 @@ ShihQuadraticKE::ShihQuadraticKE
             0.9
         )
     ),
-
+    A1_
-    kappa_
     (
         dimensioned<scalar>::lookupOrAddToDict
         (
-            "kappa_",
+            "A2",
             coeffDict_,
-            0.41
+            1000.0
         )
     ),
-    E_
+    beta1_
     (
         dimensioned<scalar>::lookupOrAddToDict
         (
-            "E",
+            "beta1",
             coeffDict_,
-            9.8
+            3.0
+        )
+    ),
+    beta2_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "beta2",
+            coeffDict_,
+            15.0
+        )
+    ),
+    beta3_
+    (
+        dimensioned<scalar>::lookupOrAddToDict
+        (
+            "beta3",
+            coeffDict_,
+            -19.0
         )
     ),
 
@@ -230,28 +215,14 @@ ShihQuadraticKE::ShihQuadraticKE
             IOobject::AUTO_WRITE
         ),
         mesh_
-    ),
+    )
-
-    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_ + pow(eta_, 3.0))
 {
     bound(k_, kMin_);
-    // already bounded: bound(epsilon_, epsilonMin_);
+    bound(epsilon_, epsilonMin_);
 
     if (type == typeName)
     {
         correctNut();
-        correctNonlinearStress(fvc::grad(U));
         printCoeffs(type);
     }
 }
@@ -267,15 +238,12 @@ bool ShihQuadraticKE::read()
         C2_.readIfPresent(coeffDict());
         sigmak_.readIfPresent(coeffDict());
         sigmaEps_.readIfPresent(coeffDict());
+        Cmu1_.readIfPresent(coeffDict());
+        Cmu2_.readIfPresent(coeffDict());
         A1_.readIfPresent(coeffDict());
-        A2_.readIfPresent(coeffDict());
+        beta1_.readIfPresent(coeffDict());
-        Ctau1_.readIfPresent(coeffDict());
+        beta2_.readIfPresent(coeffDict());
-        Ctau2_.readIfPresent(coeffDict());
+        beta3_.readIfPresent(coeffDict());
-        Ctau3_.readIfPresent(coeffDict());
-        alphaKsi_.readIfPresent(coeffDict());
-
-        kappa_.readIfPresent(coeffDict());
-        E_.readIfPresent(coeffDict());
 
         return true;
     }
@@ -298,16 +266,14 @@ void ShihQuadraticKE::correct()
     tmp<volTensorField> tgradU = fvc::grad(U_);
     const volTensorField& gradU = tgradU();
 
-    // generation term
-    tmp<volScalarField> S2 = symm(gradU) && gradU;
-
     volScalarField G
     (
         GName(),
-        Cmu_*sqr(k_)/epsilon_*S2 - (nonlinearStress_ && gradU)
+        (nut_*twoSymm(gradU) - nonlinearStress_) && gradU
     );
 
-    #include "nonLinearWallFunctionsI.H"
+    // Update epsilon and G at the wall
+    epsilon_.boundaryField().updateCoeffs();
 
     // Dissipation equation
     tmp<fvScalarMatrix> epsEqn
@@ -322,14 +288,13 @@ void ShihQuadraticKE::correct()
 
     epsEqn().relax();
 
-    #include "wallDissipationI.H"
+    epsEqn().boundaryManipulate(epsilon_.boundaryField());
 
     solve(epsEqn);
     bound(epsilon_, epsilonMin_);
 
 
     // Turbulent kinetic energy equation
-
     tmp<fvScalarMatrix> kEqn
     (
         fvm::ddt(k_)
@@ -345,14 +310,7 @@ void ShihQuadraticKE::correct()
     bound(k_, kMin_);
 
 
-    // Re-calculate viscosity
+    // Re-calculate viscosity and non-linear stress
-
-    eta_ = k_/epsilon_*sqrt(2.0*magSqr(0.5*(gradU + T(gradU))));
-    ksi_ = k_/epsilon_*sqrt(2.0*magSqr(0.5*(gradU - T(gradU))));
-    Cmu_ = 2.0/(3.0*(A1_ + eta_ + alphaKsi_*ksi_));
-    fEta_ = A2_ + pow(eta_, 3.0);
-
-    correctNut();
     correctNonlinearStress(gradU);
 }
 

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

@@ -77,15 +77,12 @@ protected:
             dimensionedScalar C2_;
             dimensionedScalar sigmak_;
             dimensionedScalar sigmaEps_;
+            dimensionedScalar Cmu1_;
+            dimensionedScalar Cmu2_;
             dimensionedScalar A1_;
-            dimensionedScalar A2_;
+            dimensionedScalar beta1_;
-            dimensionedScalar Ctau1_;
+            dimensionedScalar beta2_;
-            dimensionedScalar Ctau2_;
+            dimensionedScalar beta3_;
-            dimensionedScalar Ctau3_;
-            dimensionedScalar alphaKsi_;
-
-            dimensionedScalar kappa_;
-            dimensionedScalar E_;
 
 
         // Fields
@@ -93,11 +90,6 @@ protected:
             volScalarField k_;
             volScalarField epsilon_;
 
-            volScalarField eta_;
-            volScalarField ksi_;
-            volScalarField Cmu_;
-            volScalarField fEta_;
-
 
     // Protected Member Functions
 

src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/include/nonLinearWallFunctionsI.H

@@ -1,137 +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/>.
-
-Global
-    nonLinearwallFunctions
-
-Description
-    Calculate wall generation and dissipation from wall-functions
-    for non-linear models.
-
-\*---------------------------------------------------------------------------*/
-
-{
-    labelList cellBoundaryFaceCount(epsilon_.size(), 0);
-
-    scalar yPlusLam = nutkWallFunctionFvPatchScalarField::yPlusLam
-    (
-        kappa_.value(),
-        E_.value()
-    );
-
-    const fvPatchList& patches = mesh_.boundary();
-
-    //- Initialise the near-wall G and epsilon fields to zero
-    forAll(patches, patchi)
-    {
-        const fvPatch& curPatch = patches[patchi];
-
-        if (isA<wallFvPatch>(curPatch))
-        {
-            forAll(curPatch, facei)
-            {
-                label faceCelli = curPatch.faceCells()[facei];
-
-                epsilon_[faceCelli] = 0.0;
-                G[faceCelli] = 0.0;
-            }
-        }
-    }
-
-    const volScalarField nuLam(this->nu());
-
-    //- Accumulate the wall face contributions to epsilon and G
-    //  Increment cellBoundaryFaceCount for each face for averaging
-    forAll(patches, patchi)
-    {
-        const fvPatch& curPatch = patches[patchi];
-
-        if (isA<wallFvPatch>(curPatch))
-        {
-            #include "checkPatchFieldTypes.H"
-
-            const scalarField& nuw = nuLam.boundaryField()[patchi];
-            const scalarField& nutw = nut_.boundaryField()[patchi];
-
-            const scalarField magFaceGradU
-            (
-                mag(U_.boundaryField()[patchi].snGrad())
-            );
-
-            forAll(curPatch, facei)
-            {
-                label faceCelli = curPatch.faceCells()[facei];
-
-                //- Using local Cmu !
-                scalar Cmu25 = pow025(Cmu_[faceCelli]);
-                scalar Cmu75 = pow(Cmu_[faceCelli], 0.75);
-
-                scalar yPlus =
-                    Cmu25*y_[patchi][facei]
-                    *sqrt(k_[faceCelli])
-                    /nuw[facei];
-
-                // For corner cells (with two boundary or more faces),
-                // epsilon and G in the near-wall cell are calculated
-                // as an average
-
-                cellBoundaryFaceCount[faceCelli]++;
-
-                epsilon_[faceCelli] +=
-                     Cmu75*pow(k_[faceCelli], 1.5)
-                    /(kappa_.value()*y_[patchi][facei]);
-
-                if (yPlus > yPlusLam)
-                {
-                    G[faceCelli] +=
-                        (nutw[facei] + nuw[facei])
-                        *magFaceGradU[facei]
-                        *Cmu25*sqrt(k_[faceCelli])
-                        /(kappa_.value()*y_[patchi][facei])
-                      - (nonlinearStress_[faceCelli] && gradU[faceCelli]);
-                }
-            }
-        }
-    }
-
-    // Perform the averaging
-
-    forAll(patches, patchi)
-    {
-        const fvPatch& curPatch = patches[patchi];
-
-        if (isA<wallFvPatch>(curPatch))
-        {
-            forAll(curPatch, facei)
-            {
-                label faceCelli = curPatch.faceCells()[facei];
-
-                epsilon_[faceCelli] /= cellBoundaryFaceCount[faceCelli];
-                G[faceCelli] /= cellBoundaryFaceCount[faceCelli];
-            }
-        }
-    }
-}
-
-
-// ************************************************************************* //

src/TurbulenceModels/incompressible/turbulentTransportModels/RAS/include/wallNonlinearViscosityI.H

@@ -1,78 +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/>.
-
-Global
-    wallNonlinearViscosity
-
-Description
-    Calculate wall viscosity for non-linear models
-
-\*---------------------------------------------------------------------------*/
-
-{
-    const fvPatchList& patches = mesh_.boundary();
-
-    const scalar yPlusLam = nutkWallFunctionFvPatchScalarField::yPlusLam
-    (
-        kappa_.value(),
-        E_.value()
-    );
-
-    const volScalarField nuLam(this->nu());
-
-    forAll(patches, patchi)
-    {
-        const fvPatch& curPatch = patches[patchi];
-
-        if (isA<wallFvPatch>(curPatch))
-        {
-            const scalarField& nuw = nuLam.boundaryField()[patchi];
-            scalarField& nutw = nut_.boundaryField()[patchi];
-
-            forAll(curPatch, facei)
-            {
-                label faceCelli = curPatch.faceCells()[facei];
-
-                //- Using local Cmu
-                scalar Cmu25 = pow025(Cmu_[faceCelli]);
-
-                scalar yPlus =
-                    Cmu25*y_[patchi][facei]*sqrt(k_[faceCelli])/nuw[facei];
-
-                if (yPlus > yPlusLam)
-                {
-                    nutw[facei] =
-                        nuw[facei]
-                       *(yPlus*kappa_.value()/log(E_.value()*yPlus) - 1);
-                }
-                else
-                {
-                    nutw[facei] = 0.0;
-                }
-            }
-        }
-    }
-}
-
-
-// ************************************************************************* //

tutorials/incompressible/simpleFoam/pitzDaily/constant/turbulenceProperties

@@ -19,7 +19,9 @@ simulationType RAS;
 
 RAS
 {
-    RASModel        v2f; //kEpsilon; //realizableKE; //kOmega; //kOmegaSST;
+    // Tested with ShihQuadraticKE, v2f, kEpsilon, realizableKE,
+    // kOmega, kOmegaSST;
+    RASModel        v2f;
 
     turbulence      on;
 

tutorials/incompressible/simpleFoam/pitzDaily/system/fvSchemes

@@ -34,6 +34,7 @@ divSchemes
     div(phi,omega)  bounded Gauss limitedLinear 1;
     div(phi,v2)     bounded Gauss limitedLinear 1;
     div((nuEff*dev2(T(grad(U))))) Gauss linear;
+    div(nonlinearStress) Gauss linear;
 }
 
 laplacianSchemes