From b7aa89e1777e4eb9220a505ae43f7355a73c4a0a Mon Sep 17 00:00:00 2001 From: Andrew Heather <> Date: Tue, 17 Dec 2019 12:19:13 +0000 Subject: [PATCH] ENH: Turbulence models - ensure unlimited nut used in production terms --- .../RAS/RNGkEpsilon/RNGkEpsilon.C | 33 +++++++++++-------- .../turbulenceModels/RAS/kEpsilon/kEpsilon.C | 14 ++++---- .../turbulenceModels/RAS/kOmega/kOmega.C | 26 ++++++++------- 3 files changed, 43 insertions(+), 30 deletions(-) diff --git a/src/TurbulenceModels/turbulenceModels/RAS/RNGkEpsilon/RNGkEpsilon.C b/src/TurbulenceModels/turbulenceModels/RAS/RNGkEpsilon/RNGkEpsilon.C index ac07349339..7efc1b848b 100644 --- a/src/TurbulenceModels/turbulenceModels/RAS/RNGkEpsilon/RNGkEpsilon.C +++ b/src/TurbulenceModels/turbulenceModels/RAS/RNGkEpsilon/RNGkEpsilon.C @@ -251,23 +251,30 @@ void RNGkEpsilon::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>::correct(); - volScalarField divU(fvc::div(fvc::absolute(this->phi(), U))); + const volScalarField::Internal divU + ( + fvc::div(fvc::absolute(this->phi(), U))().v() + ); tmp 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::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::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_) ); diff --git a/src/TurbulenceModels/turbulenceModels/RAS/kEpsilon/kEpsilon.C b/src/TurbulenceModels/turbulenceModels/RAS/kEpsilon/kEpsilon.C index ca5d0392df..33704af01f 100644 --- a/src/TurbulenceModels/turbulenceModels/RAS/kEpsilon/kEpsilon.C +++ b/src/TurbulenceModels/turbulenceModels/RAS/kEpsilon/kEpsilon.C @@ -231,22 +231,24 @@ void kEpsilon::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>::correct(); - volScalarField::Internal divU + const volScalarField::Internal divU ( fvc::div(fvc::absolute(this->phi(), U))().v() ); tmp 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::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() diff --git a/src/TurbulenceModels/turbulenceModels/RAS/kOmega/kOmega.C b/src/TurbulenceModels/turbulenceModels/RAS/kOmega/kOmega.C index e73958b9ce..89ceedb5bb 100644 --- a/src/TurbulenceModels/turbulenceModels/RAS/kOmega/kOmega.C +++ b/src/TurbulenceModels/turbulenceModels/RAS/kOmega/kOmega.C @@ -191,19 +191,23 @@ void kOmega::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>::correct(); - volScalarField divU(fvc::div(fvc::absolute(this->phi(), U))); + const volScalarField::Internal divU + ( + fvc::div(fvc::absolute(this->phi(), U))().v() + ); tmp 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::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::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_) );