ENH: Turbulence models - ensure unlimited nut used in production terms

src/TurbulenceModels/turbulenceModels/RAS/RNGkEpsilon/RNGkEpsilon.C

@@ -251,23 +251,30 @@ void RNGkEpsilon<BasicTurbulenceModel>::correct()
     const rhoField& rho = this->rho_;
     const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_;
     const volVectorField& U = this->U_;
-    volScalarField& nut = this->nut_;
+    const volScalarField::Internal unlimitedNut(Cmu_*sqr(k_())/epsilon_());
+    const volScalarField& nut = this->nut_;
     fv::options& fvOptions(fv::options::New(this->mesh_));
 
     eddyViscosity<RASModel<BasicTurbulenceModel>>::correct();
 
-    volScalarField divU(fvc::div(fvc::absolute(this->phi(), U)));
+    const volScalarField::Internal divU
+    (
+        fvc::div(fvc::absolute(this->phi(), U))().v()
+    );
 
     tmp<volTensorField> tgradU = fvc::grad(U);
-    volScalarField S2((tgradU() && dev(twoSymm(tgradU()))));
+    const volScalarField::Internal GbyNu
+    (
+        tgradU().v() && dev(twoSymm(tgradU().v()))
+    );
     tgradU.clear();
 
-    volScalarField G(this->GName(), nut*S2);
+    const volScalarField::Internal G(this->GName(), unlimitedNut*GbyNu);
 
-    volScalarField eta(sqrt(mag(S2))*k_/epsilon_);
-    volScalarField eta3(eta*sqr(eta));
+    const volScalarField::Internal eta(sqrt(mag(GbyNu))*k_/epsilon_);
+    const volScalarField::Internal eta3(eta*sqr(eta));
 
-    volScalarField R
+    const volScalarField::Internal R
     (
         ((eta*(-eta/eta0_ + scalar(1)))/(beta_*eta3 + scalar(1)))
     );
@@ -282,9 +289,9 @@ void RNGkEpsilon<BasicTurbulenceModel>::correct()
       + fvm::div(alphaRhoPhi, epsilon_)
       - fvm::laplacian(alpha*rho*DepsilonEff(), epsilon_)
      ==
-        (C1_ - R)*alpha*rho*G*epsilon_/k_
-      - fvm::SuSp(((2.0/3.0)*C1_ - C3_)*alpha*rho*divU, epsilon_)
-      - fvm::Sp(C2_*alpha*rho*epsilon_/k_, epsilon_)
+        (C1_ - R)*alpha()*rho()*G*epsilon_()/k_()
+      - fvm::SuSp(((2.0/3.0)*C1_ - C3_)*alpha()*rho()*divU, epsilon_)
+      - fvm::Sp(C2_*alpha()*rho()*epsilon_()/k_(), epsilon_)
       + epsilonSource()
       + fvOptions(alpha, rho, epsilon_)
     );
@@ -305,9 +312,9 @@ void RNGkEpsilon<BasicTurbulenceModel>::correct()
       + fvm::div(alphaRhoPhi, k_)
       - fvm::laplacian(alpha*rho*DkEff(), k_)
      ==
-        alpha*rho*G
-      - fvm::SuSp((2.0/3.0)*alpha*rho*divU, k_)
-      - fvm::Sp(alpha*rho*epsilon_/k_, k_)
+        alpha()*rho()*GbyNu*nut()
+      - fvm::SuSp((2.0/3.0)*alpha()*rho()*divU, k_)
+      - fvm::Sp(alpha()*rho()*epsilon_()/k_(), k_)
       + kSource()
       + fvOptions(alpha, rho, k_)
     );

src/TurbulenceModels/turbulenceModels/RAS/kEpsilon/kEpsilon.C

@@ -231,22 +231,24 @@ void kEpsilon<BasicTurbulenceModel>::correct()
     const rhoField& rho = this->rho_;
     const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_;
     const volVectorField& U = this->U_;
-    volScalarField& nut = this->nut_;
+    const volScalarField::Internal unlimitedNut(Cmu_*sqr(k_())/epsilon_());
+    const volScalarField& nut = this->nut_;
+
     fv::options& fvOptions(fv::options::New(this->mesh_));
 
     eddyViscosity<RASModel<BasicTurbulenceModel>>::correct();
 
-    volScalarField::Internal divU
+    const volScalarField::Internal divU
     (
         fvc::div(fvc::absolute(this->phi(), U))().v()
     );
 
     tmp<volTensorField> tgradU = fvc::grad(U);
-    volScalarField::Internal G
+    const volScalarField::Internal GbyNu
     (
-        this->GName(),
-        nut.v()*(dev(twoSymm(tgradU().v())) && tgradU().v())
+        tgradU().v() && dev(twoSymm(tgradU().v()))
     );
+    const volScalarField::Internal G(this->GName(), unlimitedNut*GbyNu);
     tgradU.clear();
 
     // Update epsilon and G at the wall
@@ -280,7 +282,7 @@ void kEpsilon<BasicTurbulenceModel>::correct()
       + fvm::div(alphaRhoPhi, k_)
       - fvm::laplacian(alpha*rho*DkEff(), k_)
      ==
-        alpha()*rho()*G
+        alpha()*rho()*GbyNu*nut()
       - fvm::SuSp((2.0/3.0)*alpha()*rho()*divU, k_)
       - fvm::Sp(alpha()*rho()*epsilon_()/k_(), k_)
       + kSource()

src/TurbulenceModels/turbulenceModels/RAS/kOmega/kOmega.C

@@ -191,19 +191,23 @@ void kOmega<BasicTurbulenceModel>::correct()
     const rhoField& rho = this->rho_;
     const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_;
     const volVectorField& U = this->U_;
-    volScalarField& nut = this->nut_;
+    const volScalarField::Internal unlimitedNut(k_()/omega_());
+    const volScalarField& nut = this->nut_;
     fv::options& fvOptions(fv::options::New(this->mesh_));
 
     eddyViscosity<RASModel<BasicTurbulenceModel>>::correct();
 
-    volScalarField divU(fvc::div(fvc::absolute(this->phi(), U)));
+    const volScalarField::Internal divU
+    (
+        fvc::div(fvc::absolute(this->phi(), U))().v()
+    );
 
     tmp<volTensorField> tgradU = fvc::grad(U);
-    volScalarField G
+    const volScalarField::Internal GbyNu
     (
-        this->GName(),
-        nut*(tgradU() && dev(twoSymm(tgradU())))
+        tgradU().v() && dev(twoSymm(tgradU().v()))
     );
+    const volScalarField::Internal G(this->GName(), unlimitedNut*GbyNu);
     tgradU.clear();
 
     // Update omega and G at the wall
@@ -216,9 +220,9 @@ void kOmega<BasicTurbulenceModel>::correct()
       + fvm::div(alphaRhoPhi, omega_)
       - fvm::laplacian(alpha*rho*DomegaEff(), omega_)
      ==
-        gamma_*alpha*rho*G*omega_/k_
-      - fvm::SuSp(((2.0/3.0)*gamma_)*alpha*rho*divU, omega_)
-      - fvm::Sp(beta_*alpha*rho*omega_, omega_)
+        gamma_*alpha()*rho()*G*omega_()/k_()
+      - fvm::SuSp(((2.0/3.0)*gamma_)*alpha()*rho()*divU, omega_)
+      - fvm::Sp(beta_*alpha()*rho()*omega_(), omega_)
       + fvOptions(alpha, rho, omega_)
     );
 
@@ -237,9 +241,9 @@ void kOmega<BasicTurbulenceModel>::correct()
       + fvm::div(alphaRhoPhi, k_)
       - fvm::laplacian(alpha*rho*DkEff(), k_)
      ==
-        alpha*rho*G
-      - fvm::SuSp((2.0/3.0)*alpha*rho*divU, k_)
-      - fvm::Sp(Cmu_*alpha*rho*omega_, k_)
+        alpha()*rho()*GbyNu*nut()
+      - fvm::SuSp((2.0/3.0)*alpha()*rho()*divU, k_)
+      - fvm::Sp(Cmu_*alpha()*rho()*omega_(), k_)
       + fvOptions(alpha, rho, k_)
     );