From 195dfbe16832cfe9ce967b83164bdd39c6ce8851 Mon Sep 17 00:00:00 2001 From: Henry Weller Date: Tue, 19 Dec 2023 22:50:10 +0000 Subject: [PATCH] MomentumTransportModels::k-omega model: Improved omega bounding using nutMaxCoeff Omega lower limit bounding is now based on a maximum turbulence viscosity nut rather than a minimum omega value which improves stability and robustness of the k-omega models in case of numerical boundedness problems. The maximum nut value is calculated by multiplying the laminar viscosity by nutMaxCoeff which defaults to 1e5 but can be set by the user in the momentumTransport dictionary. --- .../applyBoundaryLayer/applyBoundaryLayer.C | 1 - .../incompressible/RAS/kkLOmega/kkLOmega.C | 15 ++++++++++- .../incompressible/RAS/kkLOmega/kkLOmega.H | 3 +++ .../Base/kOmegaSST/kOmegaSSTBase.C | 26 ++++++++++-------- .../Base/kOmegaSST/kOmegaSSTBase.H | 3 +++ .../LES/LESModel/LESModel.C | 9 ++++--- .../LES/LESModel/LESModel.H | 16 ++--------- .../RAS/RASModel/RASModel.C | 13 +-------- .../RAS/RASModel/RASModel.H | 27 ------------------- .../RAS/kOmega/kOmega.C | 12 +++++++-- .../RAS/kOmega/kOmega.H | 8 ++++++ .../RAS/kOmega2006/kOmega2006.C | 14 ++++++++-- .../RAS/kOmega2006/kOmega2006.H | 8 ++++++ 13 files changed, 81 insertions(+), 74 deletions(-) diff --git a/applications/utilities/preProcessing/applyBoundaryLayer/applyBoundaryLayer.C b/applications/utilities/preProcessing/applyBoundaryLayer/applyBoundaryLayer.C index dfabdd3554..90ff99608a 100644 --- a/applications/utilities/preProcessing/applyBoundaryLayer/applyBoundaryLayer.C +++ b/applications/utilities/preProcessing/applyBoundaryLayer/applyBoundaryLayer.C @@ -195,7 +195,6 @@ int main(int argc, char *argv[]) refCast(turbulence()); omega = (1 - mask)*omega + mask*ce0*sqrt(k)/(Cmu*min(y, ybl)); - bound(omega, rasModel.omegaMin()); // Do not correct BC - G set by the wall-functions is not available // omega.correctBoundaryConditions(); diff --git a/src/MomentumTransportModels/incompressible/RAS/kkLOmega/kkLOmega.C b/src/MomentumTransportModels/incompressible/RAS/kkLOmega/kkLOmega.C index 701a355cce..f32d185932 100644 --- a/src/MomentumTransportModels/incompressible/RAS/kkLOmega/kkLOmega.C +++ b/src/MomentumTransportModels/incompressible/RAS/kkLOmega/kkLOmega.C @@ -492,6 +492,16 @@ kkLOmega::kkLOmega 1.17 ) ), + omegaMin_ + ( + dimensioned::lookupOrAddToDict + ( + "omegaMin", + coeffDict_, + dimless/dimTime, + small + ) + ), kt_ ( IOobject @@ -536,12 +546,13 @@ kkLOmega::kkLOmega runTime_.name(), mesh_ ), - kt_*max(omega_, omegaMin_) + D(kl_) + D(kt_) + kt_*omega_ ) { bound(kt_, kMin_); bound(kl_, kMin_); bound(omega_, omegaMin_); + epsilon_ = kt_*omega_ + D(kl_) + D(kt_); if (type == typeName) { @@ -751,9 +762,11 @@ void kkLOmega::correct() solve(ktEqn); bound(kt_, kMin_); + // Update total fluctuation kinetic energy dissipation rate epsilon_ = kt_*omega_ + Dl + Dt; + // Re-calculate turbulent viscosity nut_ = nuts + nutl; nut_.correctBoundaryConditions(); diff --git a/src/MomentumTransportModels/incompressible/RAS/kkLOmega/kkLOmega.H b/src/MomentumTransportModels/incompressible/RAS/kkLOmega/kkLOmega.H index 344a381177..0f8d9fef7f 100644 --- a/src/MomentumTransportModels/incompressible/RAS/kkLOmega/kkLOmega.H +++ b/src/MomentumTransportModels/incompressible/RAS/kkLOmega/kkLOmega.H @@ -192,6 +192,9 @@ protected: dimensionedScalar Sigmak_; dimensionedScalar Sigmaw_; + //- Lower limit for omega + dimensionedScalar omegaMin_; + // Fields diff --git a/src/MomentumTransportModels/momentumTransportModels/Base/kOmegaSST/kOmegaSSTBase.C b/src/MomentumTransportModels/momentumTransportModels/Base/kOmegaSST/kOmegaSSTBase.C index e711cd6492..1fe452591a 100644 --- a/src/MomentumTransportModels/momentumTransportModels/Base/kOmegaSST/kOmegaSSTBase.C +++ b/src/MomentumTransportModels/momentumTransportModels/Base/kOmegaSST/kOmegaSSTBase.C @@ -36,9 +36,17 @@ namespace Foam // * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * * // +template +void +kOmegaSST::boundOmega() +{ + omega_ = max(omega_, k_/(this->nutMaxCoeff_*this->nu())); +} + + template tmp -kOmegaSST::kOmegaSST::F1 +kOmegaSST::F1 ( const volScalarField& CDkOmega ) const @@ -68,8 +76,7 @@ kOmegaSST::kOmegaSST::F1 template tmp -kOmegaSST::kOmegaSST:: -F2() const +kOmegaSST::F2() const { tmp arg2 = min ( @@ -86,8 +93,7 @@ F2() const template tmp -kOmegaSST::kOmegaSST:: -F3() const +kOmegaSST::F3() const { tmp arg3 = min ( @@ -100,8 +106,7 @@ F3() const template tmp -kOmegaSST::kOmegaSST:: -F23() const +kOmegaSST::F23() const { tmp f23(F2()); @@ -127,8 +132,6 @@ void kOmegaSST::correctNut } -// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * // - template void kOmegaSST:: correctNut() @@ -380,7 +383,7 @@ kOmegaSST::kOmegaSST ) { bound(k_, this->kMin_); - bound(omega_, this->omegaMin_); + boundOmega(); } @@ -493,7 +496,7 @@ void kOmegaSST::correct() omegaEqn.ref().boundaryManipulate(omega_.boundaryFieldRef()); solve(omegaEqn); fvConstraints.constrain(omega_); - bound(omega_, this->omegaMin_); + boundOmega(); } // Turbulent kinetic energy equation @@ -515,6 +518,7 @@ void kOmegaSST::correct() solve(kEqn); fvConstraints.constrain(k_); bound(k_, this->kMin_); + boundOmega(); correctNut(S2, F23); } diff --git a/src/MomentumTransportModels/momentumTransportModels/Base/kOmegaSST/kOmegaSSTBase.H b/src/MomentumTransportModels/momentumTransportModels/Base/kOmegaSST/kOmegaSSTBase.H index 75976ab577..04ce0e440e 100644 --- a/src/MomentumTransportModels/momentumTransportModels/Base/kOmegaSST/kOmegaSSTBase.H +++ b/src/MomentumTransportModels/momentumTransportModels/Base/kOmegaSST/kOmegaSSTBase.H @@ -200,6 +200,9 @@ protected: return blend(F1, gamma1_, gamma2_); } + //- Bound omega + void boundOmega(); + virtual void correctNut ( const volScalarField& S2, diff --git a/src/MomentumTransportModels/momentumTransportModels/LES/LESModel/LESModel.C b/src/MomentumTransportModels/momentumTransportModels/LES/LESModel/LESModel.C index 33f332aa36..50d5b469a5 100644 --- a/src/MomentumTransportModels/momentumTransportModels/LES/LESModel/LESModel.C +++ b/src/MomentumTransportModels/momentumTransportModels/LES/LESModel/LESModel.C @@ -79,14 +79,14 @@ Foam::LESModel::LESModel ) ), - omegaMin_ + nutMaxCoeff_ ( dimensioned::lookupOrAddToDict ( - "omegaMin", + "nutMaxCoeff", LESDict_, - dimless/dimTime, - small + dimless, + 1e5 ) ), @@ -192,6 +192,7 @@ bool Foam::LESModel::read() delta_().read(LESDict_); kMin_.readIfPresent(LESDict_); + nutMaxCoeff_.readIfPresent(LESDict_); return true; } diff --git a/src/MomentumTransportModels/momentumTransportModels/LES/LESModel/LESModel.H b/src/MomentumTransportModels/momentumTransportModels/LES/LESModel/LESModel.H index c056801041..4fdab12756 100644 --- a/src/MomentumTransportModels/momentumTransportModels/LES/LESModel/LESModel.H +++ b/src/MomentumTransportModels/momentumTransportModels/LES/LESModel/LESModel.H @@ -82,8 +82,8 @@ protected: //- Lower limit of k dimensionedScalar kMin_; - //- Lower limit for omega - dimensionedScalar omegaMin_; + //- Upper limit coefficient for nut + dimensionedScalar nutMaxCoeff_; //- Run-time selectable generalised Newtonian viscosity model autoPtr @@ -176,18 +176,6 @@ public: return coeffDict_; } - //- Return the lower allowable limit for k (default: small) - const dimensionedScalar& kMin() const - { - return kMin_; - } - - //- Allow kMin to be changed - dimensionedScalar& kMin() - { - return kMin_; - } - //- Access function to filter width inline const volScalarField& delta() const { diff --git a/src/MomentumTransportModels/momentumTransportModels/RAS/RASModel/RASModel.C b/src/MomentumTransportModels/momentumTransportModels/RAS/RASModel/RASModel.C index d6916dd170..d9fe2b2433 100644 --- a/src/MomentumTransportModels/momentumTransportModels/RAS/RASModel/RASModel.C +++ b/src/MomentumTransportModels/momentumTransportModels/RAS/RASModel/RASModel.C @@ -79,17 +79,6 @@ Foam::RASModel::RASModel ) ), - omegaMin_ - ( - dimensioned::lookupOrAddToDict - ( - "omegaMin", - RASDict_, - dimless/dimTime, - small - ) - ), - nutMaxCoeff_ ( dimensioned::lookupOrAddToDict @@ -191,7 +180,7 @@ bool Foam::RASModel::read() coeffDict_ <<= RASDict_.optionalSubDict(type() + "Coeffs"); kMin_.readIfPresent(RASDict_); - omegaMin_.readIfPresent(RASDict_); + nutMaxCoeff_.readIfPresent(RASDict_); return true; } diff --git a/src/MomentumTransportModels/momentumTransportModels/RAS/RASModel/RASModel.H b/src/MomentumTransportModels/momentumTransportModels/RAS/RASModel/RASModel.H index 6a63883c7c..9cdb20a411 100644 --- a/src/MomentumTransportModels/momentumTransportModels/RAS/RASModel/RASModel.H +++ b/src/MomentumTransportModels/momentumTransportModels/RAS/RASModel/RASModel.H @@ -75,9 +75,6 @@ protected: //- Lower limit of k dimensionedScalar kMin_; - //- Lower limit for omega - dimensionedScalar omegaMin_; - //- Upper limit coefficient for nut dimensionedScalar nutMaxCoeff_; @@ -163,30 +160,6 @@ public: //- Read model coefficients if they have changed virtual bool read(); - //- Return the lower allowable limit for k (default: small) - const dimensionedScalar& kMin() const - { - return kMin_; - } - - //- Return the lower allowable limit for omega (default: small) - const dimensionedScalar& omegaMin() const - { - return omegaMin_; - } - - //- Allow kMin to be changed - dimensionedScalar& kMin() - { - return kMin_; - } - - //- Allow omegaMin to be changed - dimensionedScalar& omegaMin() - { - return omegaMin_; - } - //- Const access to the coefficients dictionary virtual const dictionary& coeffDict() const { diff --git a/src/MomentumTransportModels/momentumTransportModels/RAS/kOmega/kOmega.C b/src/MomentumTransportModels/momentumTransportModels/RAS/kOmega/kOmega.C index fb80ee2c9a..ec8fc1443e 100644 --- a/src/MomentumTransportModels/momentumTransportModels/RAS/kOmega/kOmega.C +++ b/src/MomentumTransportModels/momentumTransportModels/RAS/kOmega/kOmega.C @@ -37,6 +37,13 @@ namespace RASModels // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * // +template +void kOmega::boundOmega() +{ + omega_ = max(omega_, k_/(this->nutMaxCoeff_*this->nu())); +} + + template void kOmega::correctNut() { @@ -172,7 +179,7 @@ kOmega::kOmega ) { bound(k_, this->kMin_); - bound(omega_, this->omegaMin_); + boundOmega(); if (type == typeName) { @@ -260,7 +267,7 @@ void kOmega::correct() omegaEqn.ref().boundaryManipulate(omega_.boundaryFieldRef()); solve(omegaEqn); fvConstraints.constrain(omega_); - bound(omega_, this->omegaMin_); + boundOmega(); // Turbulent kinetic energy equation @@ -282,6 +289,7 @@ void kOmega::correct() solve(kEqn); fvConstraints.constrain(k_); bound(k_, this->kMin_); + boundOmega(); correctNut(); } diff --git a/src/MomentumTransportModels/momentumTransportModels/RAS/kOmega/kOmega.H b/src/MomentumTransportModels/momentumTransportModels/RAS/kOmega/kOmega.H index a7ed0a1052..200f603cac 100644 --- a/src/MomentumTransportModels/momentumTransportModels/RAS/kOmega/kOmega.H +++ b/src/MomentumTransportModels/momentumTransportModels/RAS/kOmega/kOmega.H @@ -96,8 +96,16 @@ protected: // Protected Member Functions + //- Bound omega + void boundOmega(); + + //- Correct the eddy-viscosity nut virtual void correctNut(); + + //- Source term for the k equation virtual tmp kSource() const; + + //- Source term for the omega equation virtual tmp omegaSource() const; diff --git a/src/MomentumTransportModels/momentumTransportModels/RAS/kOmega2006/kOmega2006.C b/src/MomentumTransportModels/momentumTransportModels/RAS/kOmega2006/kOmega2006.C index 37b159c3a3..c495f67cba 100644 --- a/src/MomentumTransportModels/momentumTransportModels/RAS/kOmega2006/kOmega2006.C +++ b/src/MomentumTransportModels/momentumTransportModels/RAS/kOmega2006/kOmega2006.C @@ -35,6 +35,15 @@ namespace Foam namespace RASModels { +// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * // + +template +void kOmega2006::boundOmega() +{ + omega_ = max(omega_, k_/(this->nutMaxCoeff_*this->nu())); +} + + template void kOmega2006::correctNut ( @@ -238,7 +247,7 @@ kOmega2006::kOmega2006 ) { bound(k_, this->kMin_); - bound(omega_, this->omegaMin_); + boundOmega(); if (type == typeName) { @@ -328,7 +337,7 @@ void kOmega2006::correct() omegaEqn.ref().boundaryManipulate(omega_.boundaryFieldRef()); solve(omegaEqn); fvConstraints.constrain(omega_); - bound(omega_, this->omegaMin_); + boundOmega(); // Turbulent kinetic energy equation @@ -350,6 +359,7 @@ void kOmega2006::correct() solve(kEqn); fvConstraints.constrain(k_); bound(k_, this->kMin_); + boundOmega(); correctNut(gradU); } diff --git a/src/MomentumTransportModels/momentumTransportModels/RAS/kOmega2006/kOmega2006.H b/src/MomentumTransportModels/momentumTransportModels/RAS/kOmega2006/kOmega2006.H index ef352b34ed..59f29e2d33 100644 --- a/src/MomentumTransportModels/momentumTransportModels/RAS/kOmega2006/kOmega2006.H +++ b/src/MomentumTransportModels/momentumTransportModels/RAS/kOmega2006/kOmega2006.H @@ -111,8 +111,16 @@ protected: // Protected Member Functions + //- Bound omega + void boundOmega(); + + //- Correct the eddy-viscosity nut virtual void correctNut(); + + //- Source term for the k equation virtual tmp kSource() const; + + //- Source term for the omega equation virtual tmp omegaSource() const;