Improved linearisation of coupling. Robert Keser

src/turbulenceModels/incompressible/RAS/coupledKOmegaSST/coupledKOmegaSST.C

@@ -86,6 +86,31 @@ tmp<volScalarField> coupledKOmegaSST::F2() const
 }
 
 
+tmp<volScalarField> coupledKOmegaSST::F3() const
+{
+    tmp<volScalarField> arg3 = min
+    (
+        150*nu()/(omega_*sqr(y_)),
+        scalar(10)
+    );
+
+    return 1 - tanh(pow4(arg3));
+}
+
+
+tmp<volScalarField> coupledKOmegaSST::F23() const
+{
+    tmp<volScalarField> f23(F2());
+
+    if (F3_)
+    {
+        f23() *= F3();
+    }
+
+    return f23;
+}
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 coupledKOmegaSST::coupledKOmegaSST
@@ -99,16 +124,16 @@ coupledKOmegaSST::coupledKOmegaSST
 
     alphaK1_
     (
-        dimensioned<scalar>::lookupOrAddToDict
+        dimensionedScalar::lookupOrAddToDict
         (
             "alphaK1",
             coeffDict_,
-            0.85034
+            0.85
         )
     ),
     alphaK2_
     (
-        dimensioned<scalar>::lookupOrAddToDict
+        dimensionedScalar::lookupOrAddToDict
         (
             "alphaK2",
             coeffDict_,
@@ -117,7 +142,7 @@ coupledKOmegaSST::coupledKOmegaSST
     ),
     alphaOmega1_
     (
-        dimensioned<scalar>::lookupOrAddToDict
+        dimensionedScalar::lookupOrAddToDict
         (
             "alphaOmega1",
             coeffDict_,
@@ -126,34 +151,34 @@ coupledKOmegaSST::coupledKOmegaSST
     ),
     alphaOmega2_
     (
-        dimensioned<scalar>::lookupOrAddToDict
+        dimensionedScalar::lookupOrAddToDict
         (
             "alphaOmega2",
             coeffDict_,
-            0.85616
+            0.856
         )
     ),
     gamma1_
     (
-        dimensioned<scalar>::lookupOrAddToDict
+        dimensionedScalar::lookupOrAddToDict
         (
             "gamma1",
             coeffDict_,
-            0.5532
+            5.0/9.0
         )
     ),
     gamma2_
     (
-        dimensioned<scalar>::lookupOrAddToDict
+        dimensionedScalar::lookupOrAddToDict
         (
             "gamma2",
             coeffDict_,
-            0.4403
+            0.44
         )
     ),
     beta1_
     (
-        dimensioned<scalar>::lookupOrAddToDict
+        dimensionedScalar::lookupOrAddToDict
         (
             "beta1",
             coeffDict_,
@@ -162,7 +187,7 @@ coupledKOmegaSST::coupledKOmegaSST
     ),
     beta2_
     (
-        dimensioned<scalar>::lookupOrAddToDict
+        dimensionedScalar::lookupOrAddToDict
         (
             "beta2",
             coeffDict_,
@@ -171,7 +196,7 @@ coupledKOmegaSST::coupledKOmegaSST
     ),
     betaStar_
     (
-        dimensioned<scalar>::lookupOrAddToDict
+        dimensionedScalar::lookupOrAddToDict
         (
             "betaStar",
             coeffDict_,
@@ -180,22 +205,40 @@ coupledKOmegaSST::coupledKOmegaSST
     ),
     a1_
     (
-        dimensioned<scalar>::lookupOrAddToDict
+        dimensionedScalar::lookupOrAddToDict
         (
             "a1",
             coeffDict_,
             0.31
         )
     ),
+    b1_
+    (
+        dimensionedScalar::lookupOrAddToDict
+        (
+            "b1",
+            coeffDict_,
+            1.0
+        )
+    ),
     c1_
     (
-        dimensioned<scalar>::lookupOrAddToDict
+        dimensionedScalar::lookupOrAddToDict
         (
             "c1",
             coeffDict_,
             10.0
         )
     ),
+    F3_
+    (
+        Switch::lookupOrAddToDict
+        (
+            "F3",
+            coeffDict_,
+            false
+        )
+    ),
 
     y_(mesh_),
 
@@ -249,9 +292,18 @@ coupledKOmegaSST::coupledKOmegaSST
         dimensionedVector2("zero", dimless, vector2::zero)
     )
 {
+    bound(k_, k0_);
     bound(omega_, omega0_);
 
-    nut_ = a1_*k_/max(a1_*omega_, F2()*sqrt(2.0)*mag(symm(fvc::grad(U_))));
+    nut_ =
+    (
+        a1_*k_/
+        max
+        (
+            a1_*omega_,
+            b1_*F23()*sqrt(2.0)*mag(symm(fvc::grad(U_)))
+        )
+    );
     nut_.correctBoundaryConditions();
 
     printCoeffs();
@@ -325,7 +377,9 @@ bool coupledKOmegaSST::read()
         beta2_.readIfPresent(coeffDict());
         betaStar_.readIfPresent(coeffDict());
         a1_.readIfPresent(coeffDict());
+        b1_.readIfPresent(coeffDict());
         c1_.readIfPresent(coeffDict());
+        F3_.readIfPresent("F3", coeffDict());
 
         return true;
     }
@@ -358,8 +412,8 @@ void coupledKOmegaSST::correct()
         y_.correct();
     }
 
-    volScalarField S2 = magSqr(symm(fvc::grad(U_)));
+    const volScalarField S2(2*magSqr(symm(fvc::grad(U_))));
-    volScalarField G("RASModel::G", nut_*2*S2);
+    volScalarField G("RASModel::G", nut_*S2);
 
     // Make coupled matrix
     fvBlockMatrix<vector2> kOmegaEqn(kOmega_);
@@ -367,12 +421,17 @@ void coupledKOmegaSST::correct()
     // Update omega and G at the wall
     omega_.boundaryField().updateCoeffs();
 
-    volScalarField CDkOmega =
+    const volScalarField CDkOmega
-        (2*alphaOmega2_)*(fvc::grad(k_) & fvc::grad(omega_))/omega_;
+    (
+        (2*alphaOmega2_)*(fvc::grad(k_) & fvc::grad(omega_))/omega_
+    );
 
-    volScalarField F1 = this->F1(CDkOmega);
+    const volScalarField F1(this->F1(CDkOmega));
 
-    Info<< "Coupled k-omega" << endl;
+    // Record coupling equations when omega is set.  Same equations
+    // will be eliminated from k coupling as well due to wall functions
+    // HJ, 13/Nov/2015
+    labelList eliminated;
 
     // Turbulent frequency equation
     {
@@ -382,14 +441,20 @@ void coupledKOmegaSST::correct()
           + fvm::div(phi_, omega_)
           + fvm::SuSp(-fvc::div(phi_), omega_)
           - fvm::laplacian(DomegaEff(F1), omega_)
-          + fvm::SuSp
+         ==
+            gamma(F1)
+           *min
             (
-                beta(F1)*omega_
+                S2,
-              + (F1 - scalar(1))*CDkOmega/omega_,
+                (c1_/a1_)*betaStar_*omega_*max(a1_*omega_, b1_*F23()*sqrt(S2))
+            )
+          - fvm::Sp(2*beta(F1)*omega_, omega_)
+          + beta(F1)*sqr(omega_)
+          - fvm::SuSp
+            (
+                (F1 - scalar(1))*CDkOmega/omega_,
                 omega_
             )
-         ==
-            gamma(F1)*2*S2
         );
 
         omegaEqn.relax();
@@ -398,23 +463,25 @@ void coupledKOmegaSST::correct()
         kOmegaEqn.insertEquation(1, omegaEqn);
 
         // Add coupling term
-        volScalarField coupling
+        volScalarField couplingOmega
         (
-            "coupling",
+            "couplingOmega",
-            -gamma(F1)*2*S2/k_
+            -beta(F1)*sqr(omega_)/k_
         );
-        scalarField& couplingIn = coupling.internalField();
+        scalarField& couplingOmegaIn = couplingOmega.internalField();
 
-        // Eliminate coupling in wall function cells
+        // Collect cell indices where wall functions are active from the
-        labelList eliminated = omegaEqn.eliminatedEqns().toc();
+        // omega wall functions
+        eliminated = omegaEqn.eliminatedEqns().toc();
 
+        // Eliminate coupling in wall function cells for omega
         forAll (eliminated, cellI)
         {
-            couplingIn[eliminated[cellI]] *= 0;
+            couplingOmegaIn[eliminated[cellI]] *= 0;
         }
 
         // Insert coupling
-        kOmegaEqn.insertEquationCoupling(1, 0, coupling);
+        kOmegaEqn.insertEquationCoupling(1, 0, couplingOmega);
    }
 
     // Turbulent kinetic energy equation
@@ -425,17 +492,32 @@ void coupledKOmegaSST::correct()
           + fvm::div(phi_, k_)
           + fvm::SuSp(-fvc::div(phi_), k_)
           - fvm::laplacian(DkEff(F1), k_)
-          + fvm::SuSp
+         ==
-            (
+            min(G, c1_*betaStar_*k_*omega_)
-                betaStar_*omega_
+          - fvm::Sp(betaStar_*omega_, k_)
-              - min(G/k_, c1_*betaStar_*omega_),
-                k_
-            )
         );
 
         kEqn.relax();
+        kEqn.completeAssembly();
 
         kOmegaEqn.insertEquation(0, kEqn);
+
+        // Add coupling term
+        volScalarField couplingK
+        (
+            "couplingK",
+           betaStar_*k_
+        );
+        scalarField& couplingKIn = couplingK.internalField();
+
+        // Eliminate coupling in wall function cells for k
+        forAll (eliminated, cellI)
+        {
+            couplingKIn[eliminated[cellI]] *= 0;
+        }
+
+        // Insert coupling
+//         kOmegaEqn.insertEquationCoupling(0, 1, couplingK);
     }
 
     // Update source coupling: coupling terms eliminated from source
@@ -447,14 +529,16 @@ void coupledKOmegaSST::correct()
     kOmegaEqn.retrieveSolution(0, k_.internalField());
     kOmegaEqn.retrieveSolution(1, omega_.internalField());
 
-    k_.correctBoundaryConditions();
-    omega_.correctBoundaryConditions();
-
     bound(k_, k0_);
     bound(omega_, omega0_);
 
+    k_.correctBoundaryConditions();
+    omega_.correctBoundaryConditions();
+
+
     // Re-calculate viscosity
-    nut_ = a1_*k_/max(a1_*omega_, F2()*sqrt(2*S2));
+    // Fixed sqrt(2) error.  HJ, 10/Jun/2015
+    nut_ = a1_*k_/max(a1_*omega_, b1_*F23()*sqrt(S2));
     nut_ = min(nut_, nuRatio()*nu());
     nut_.correctBoundaryConditions();
 }

src/turbulenceModels/incompressible/RAS/coupledKOmegaSST/coupledKOmegaSST.H

@@ -33,9 +33,32 @@ Description
         Menter, F., Esch, T.
         "Elements of Industrial Heat Transfer Prediction"
         16th Brazilian Congress of Mechanical Engineering (COBEM),
-        Nov. 2001
+        Nov. 2001.
     @endverbatim
 
+    with updated coefficients from
+    @verbatim
+        Menter, F. R., Kuntz, M., and Langtry, R.,
+        "Ten Years of Industrial Experience with the SST Turbulence Model",
+        Turbulence, Heat and Mass Transfer 4, 2003,
+        pp. 625 - 632.
+    @endverbatim
+
+    but with the consistent production terms from the 2001 paper as form in the
+    2003 paper is a typo, see
+    @verbatim
+        http://turbmodels.larc.nasa.gov/sst.html
+    @endverbatim
+
+    and the addition of the optional F3 term for rough walls from
+    \verbatim
+        Hellsten, A.
+        "Some Improvements in Menter’s k-omega-SST turbulence model"
+        29th AIAA Fluid Dynamics Conference,
+        AIAA-98-2554,
+        June 1998.
+    \endverbatim
+
     Note that this implementation is written in terms of alpha diffusion
     coefficients rather than the more traditional sigma (alpha = 1/sigma) so
     that the blending can be applied to all coefficuients in a consistent
@@ -45,13 +68,6 @@ Description
     Also note that the error in the last term of equation (2) relating to
     sigma has been corrected.
 
-    Wall-functions are applied in this implementation by using equations (14)
-    to specify the near-wall omega as appropriate.
-
-    The blending functions (15) and (16) are not currently used because of the
-    uncertainty in their origin, range of applicability and that is y+ becomes
-    sufficiently small blending u_tau in this manner clearly becomes nonsense.
-
     The default model coefficients correspond to the following:
     @verbatim
         coupledKOmegaSSTCoeffs
@@ -63,14 +79,17 @@ Description
             beta1       0.075;
             beta2       0.0828;
             betaStar    0.09;
-            gamma1      0.5532;
+            gamma1      5/9;
-            gamma2      0.4403;
+            gamma2      0.44;
             a1          0.31;
+            b1          1.0;
             c1          10.0;
+            F3          no;
         }
     @endverbatim
 
 Author
+    Robert Keser, FMENA.  All rights reserved.
     Hrvoje Jasak, Wikki Ltd.  All rights reserved.
 
 SourceFiles
@@ -119,8 +138,11 @@ class coupledKOmegaSST
             dimensionedScalar betaStar_;
 
             dimensionedScalar a1_;
+            dimensionedScalar b1_;
             dimensionedScalar c1_;
 
+            Switch F3_;
+
 
         //- Wall distance field
         //  Note: different to wall distance in parent RASModel
@@ -141,6 +163,8 @@ class coupledKOmegaSST
 
         tmp<volScalarField> F1(const volScalarField& CDkOmega) const;
         tmp<volScalarField> F2() const;
+        tmp<volScalarField> F3() const;
+        tmp<volScalarField> F23() const;
 
         tmp<volScalarField> blend
         (
@@ -203,7 +227,6 @@ public:
 
 
     //- Destructor
-
     virtual ~coupledKOmegaSST()
     {}