From 3e783accd293d33e94f13b527a8e463d2fa9dcbc Mon Sep 17 00:00:00 2001 From: Andrew Heather Date: Mon, 30 Nov 2015 16:17:21 +0000 Subject: [PATCH] ENH: Turbulence - refactored kOmegaSST models to enable easier addition of derived types --- .../RAS/kOmegaSST/kOmegaSST.C | 405 +-------------- .../RAS/kOmegaSST/kOmegaSST.H | 165 +----- .../RAS/kOmegaSST/kOmegaSSTBase.C | 480 ++++++++++++++++++ .../RAS/kOmegaSST/kOmegaSSTBase.H | 320 ++++++++++++ .../{ => doc}/turbulenceModelDoc.H | 0 5 files changed, 819 insertions(+), 551 deletions(-) create mode 100644 src/TurbulenceModels/turbulenceModels/RAS/kOmegaSST/kOmegaSSTBase.C create mode 100644 src/TurbulenceModels/turbulenceModels/RAS/kOmegaSST/kOmegaSSTBase.H rename src/TurbulenceModels/turbulenceModels/{ => doc}/turbulenceModelDoc.H (100%) diff --git a/src/TurbulenceModels/turbulenceModels/RAS/kOmegaSST/kOmegaSST.C b/src/TurbulenceModels/turbulenceModels/RAS/kOmegaSST/kOmegaSST.C index 76fd9cf0c3..4353b5a748 100644 --- a/src/TurbulenceModels/turbulenceModels/RAS/kOmegaSST/kOmegaSST.C +++ b/src/TurbulenceModels/turbulenceModels/RAS/kOmegaSST/kOmegaSST.C @@ -24,8 +24,6 @@ License \*---------------------------------------------------------------------------*/ #include "kOmegaSST.H" -#include "bound.H" -#include "wallDist.H" // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // @@ -34,91 +32,22 @@ namespace Foam namespace RASModels { -// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * // - -template -tmp kOmegaSST::kOmegaSST::F1 -( - const volScalarField& CDkOmega -) const -{ - tmp CDkOmegaPlus = max - ( - CDkOmega, - dimensionedScalar("1.0e-10", dimless/sqr(dimTime), 1.0e-10) - ); - - tmp arg1 = min - ( - min - ( - max - ( - (scalar(1)/betaStar_)*sqrt(k_)/(omega_*y_), - scalar(500)*(this->mu()/this->rho_)/(sqr(y_)*omega_) - ), - (4*alphaOmega2_)*k_/(CDkOmegaPlus*sqr(y_)) - ), - scalar(10) - ); - - return tanh(pow4(arg1)); -} - -template -tmp kOmegaSST::kOmegaSST::F2() const -{ - tmp arg2 = min - ( - max - ( - (scalar(2)/betaStar_)*sqrt(k_)/(omega_*y_), - scalar(500)*(this->mu()/this->rho_)/(sqr(y_)*omega_) - ), - scalar(100) - ); - - return tanh(sqr(arg2)); -} - -template -tmp kOmegaSST::kOmegaSST::F3() const -{ - tmp arg3 = min - ( - 150*(this->mu()/this->rho_)/(omega_*sqr(y_)), - scalar(10) - ); - - return 1 - tanh(pow4(arg3)); -} - -template -tmp kOmegaSST::kOmegaSST::F23() const -{ - tmp f23(F2()); - - if (F3_) - { - f23() *= F3(); - } - - return f23; -} - +// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * // template void kOmegaSST::correctNut(const volScalarField& S2) { - this->nut_ = a1_*k_/max(a1_*omega_, b1_*F23()*sqrt(S2)); - this->nut_.correctBoundaryConditions(); + // Correct the turbulence viscosity + kOmegaSSTBase > >::correctNut + ( + S2 + ); + // Correct the turbulence thermal diffusivity BasicTurbulenceModel::correctNut(); } -// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * // - template void kOmegaSST::correctNut() { @@ -126,53 +55,6 @@ void kOmegaSST::correctNut() } -template -tmp kOmegaSST::kSource() const -{ - return tmp - ( - new fvScalarMatrix - ( - k_, - dimVolume*this->rho_.dimensions()*k_.dimensions()/dimTime - ) - ); -} - - -template -tmp kOmegaSST::omegaSource() const -{ - return tmp - ( - new fvScalarMatrix - ( - omega_, - dimVolume*this->rho_.dimensions()*omega_.dimensions()/dimTime - ) - ); -} - - -template -tmp kOmegaSST::Qsas -( - const volScalarField& S2, - const volScalarField& gamma, - const volScalarField& beta -) const -{ - return tmp - ( - new fvScalarMatrix - ( - omega_, - dimVolume*this->rho_.dimensions()*omega_.dimensions()/dimTime - ) - ); -} - - // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // template @@ -188,7 +70,7 @@ kOmegaSST::kOmegaSST const word& type ) : - eddyViscosity > + kOmegaSSTBase > > ( type, alpha, @@ -198,286 +80,17 @@ kOmegaSST::kOmegaSST phi, transport, propertiesName - ), - - alphaK1_ - ( - dimensioned::lookupOrAddToDict - ( - "alphaK1", - this->coeffDict_, - 0.85 - ) - ), - alphaK2_ - ( - dimensioned::lookupOrAddToDict - ( - "alphaK2", - this->coeffDict_, - 1.0 - ) - ), - alphaOmega1_ - ( - dimensioned::lookupOrAddToDict - ( - "alphaOmega1", - this->coeffDict_, - 0.5 - ) - ), - alphaOmega2_ - ( - dimensioned::lookupOrAddToDict - ( - "alphaOmega2", - this->coeffDict_, - 0.856 - ) - ), - gamma1_ - ( - dimensioned::lookupOrAddToDict - ( - "gamma1", - this->coeffDict_, - 5.0/9.0 - ) - ), - gamma2_ - ( - dimensioned::lookupOrAddToDict - ( - "gamma2", - this->coeffDict_, - 0.44 - ) - ), - beta1_ - ( - dimensioned::lookupOrAddToDict - ( - "beta1", - this->coeffDict_, - 0.075 - ) - ), - beta2_ - ( - dimensioned::lookupOrAddToDict - ( - "beta2", - this->coeffDict_, - 0.0828 - ) - ), - betaStar_ - ( - dimensioned::lookupOrAddToDict - ( - "betaStar", - this->coeffDict_, - 0.09 - ) - ), - a1_ - ( - dimensioned::lookupOrAddToDict - ( - "a1", - this->coeffDict_, - 0.31 - ) - ), - b1_ - ( - dimensioned::lookupOrAddToDict - ( - "b1", - this->coeffDict_, - 1.0 - ) - ), - c1_ - ( - dimensioned::lookupOrAddToDict - ( - "c1", - this->coeffDict_, - 10.0 - ) - ), - F3_ - ( - Switch::lookupOrAddToDict - ( - "F3", - this->coeffDict_, - false - ) - ), - - y_(wallDist::New(this->mesh_).y()), - - k_ - ( - IOobject - ( - IOobject::groupName("k", U.group()), - this->runTime_.timeName(), - this->mesh_, - IOobject::MUST_READ, - IOobject::AUTO_WRITE - ), - this->mesh_ - ), - omega_ - ( - IOobject - ( - IOobject::groupName("omega", U.group()), - this->runTime_.timeName(), - this->mesh_, - IOobject::MUST_READ, - IOobject::AUTO_WRITE - ), - this->mesh_ ) { - bound(k_, this->kMin_); - bound(omega_, this->omegaMin_); + correctNut(); if (type == typeName) { - correctNut(); this->printCoeffs(type); } } -// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // - -template -bool kOmegaSST::read() -{ - if (eddyViscosity >::read()) - { - alphaK1_.readIfPresent(this->coeffDict()); - alphaK2_.readIfPresent(this->coeffDict()); - alphaOmega1_.readIfPresent(this->coeffDict()); - alphaOmega2_.readIfPresent(this->coeffDict()); - gamma1_.readIfPresent(this->coeffDict()); - gamma2_.readIfPresent(this->coeffDict()); - beta1_.readIfPresent(this->coeffDict()); - beta2_.readIfPresent(this->coeffDict()); - betaStar_.readIfPresent(this->coeffDict()); - a1_.readIfPresent(this->coeffDict()); - b1_.readIfPresent(this->coeffDict()); - c1_.readIfPresent(this->coeffDict()); - F3_.readIfPresent("F3", this->coeffDict()); - - return true; - } - else - { - return false; - } -} - - -template -void kOmegaSST::correct() -{ - if (!this->turbulence_) - { - return; - } - - // Local references - const alphaField& alpha = this->alpha_; - const rhoField& rho = this->rho_; - const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_; - const volVectorField& U = this->U_; - volScalarField& nut = this->nut_; - - eddyViscosity >::correct(); - - volScalarField divU(fvc::div(fvc::absolute(this->phi(), U))); - - tmp tgradU = fvc::grad(U); - volScalarField S2(2*magSqr(symm(tgradU()))); - volScalarField GbyNu((tgradU() && dev(twoSymm(tgradU())))); - volScalarField G(this->GName(), nut*GbyNu); - tgradU.clear(); - - // Update omega and G at the wall - omega_.boundaryField().updateCoeffs(); - - volScalarField CDkOmega - ( - (2*alphaOmega2_)*(fvc::grad(k_) & fvc::grad(omega_))/omega_ - ); - - volScalarField F1(this->F1(CDkOmega)); - - { - volScalarField gamma(this->gamma(F1)); - volScalarField beta(this->beta(F1)); - - // Turbulent frequency equation - tmp omegaEqn - ( - fvm::ddt(alpha, rho, omega_) - + fvm::div(alphaRhoPhi, omega_) - - fvm::laplacian(alpha*rho*DomegaEff(F1), omega_) - == - alpha*rho*gamma - *min - ( - GbyNu, - (c1_/a1_)*betaStar_*omega_*max(a1_*omega_, b1_*F23()*sqrt(S2)) - ) - - fvm::SuSp((2.0/3.0)*alpha*rho*gamma*divU, omega_) - - fvm::Sp(alpha*rho*beta*omega_, omega_) - - fvm::SuSp - ( - alpha*rho*(F1 - scalar(1))*CDkOmega/omega_, - omega_ - ) - + Qsas(S2, gamma, beta) - + omegaSource() - ); - - omegaEqn().relax(); - - omegaEqn().boundaryManipulate(omega_.boundaryField()); - - solve(omegaEqn); - bound(omega_, this->omegaMin_); - } - - // Turbulent kinetic energy equation - tmp kEqn - ( - fvm::ddt(alpha, rho, k_) - + fvm::div(alphaRhoPhi, k_) - - fvm::laplacian(alpha*rho*DkEff(F1), k_) - == - min(alpha*rho*G, (c1_*betaStar_)*alpha*rho*k_*omega_) - - fvm::SuSp((2.0/3.0)*alpha*rho*divU, k_) - - fvm::Sp(alpha*rho*betaStar_*omega_, k_) - + kSource() - ); - - kEqn().relax(); - solve(kEqn); - bound(k_, this->kMin_); - - correctNut(S2); -} - - // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // } // End namespace RASModels diff --git a/src/TurbulenceModels/turbulenceModels/RAS/kOmegaSST/kOmegaSST.H b/src/TurbulenceModels/turbulenceModels/RAS/kOmegaSST/kOmegaSST.H index 1bbbedc91d..17327ad9dd 100644 --- a/src/TurbulenceModels/turbulenceModels/RAS/kOmegaSST/kOmegaSST.H +++ b/src/TurbulenceModels/turbulenceModels/RAS/kOmegaSST/kOmegaSST.H @@ -31,6 +31,8 @@ Description Implementation of the k-omega-SST turbulence model for incompressible and compressible flows. + Light wrapper around base class. + Turbulence model described in: \verbatim Menter, F. R. & Esch, T. (2001). @@ -98,11 +100,15 @@ Description SourceFiles kOmegaSST.C +SeeAlso + kOmegaSSTBase.H + \*---------------------------------------------------------------------------*/ #ifndef kOmegaSST_H #define kOmegaSST_H +#include "kOmegaSSTBase.H" #include "RASModel.H" #include "eddyViscosity.H" @@ -114,110 +120,21 @@ namespace RASModels { /*---------------------------------------------------------------------------*\ - Class kOmegaSST Declaration + Class kOmegaSST Declaration \*---------------------------------------------------------------------------*/ template class kOmegaSST : - public eddyViscosity > + public kOmegaSSTBase > > { - // Private Member Functions - - // Disallow default bitwise copy construct and assignment - kOmegaSST(const kOmegaSST&); - kOmegaSST& operator=(const kOmegaSST&); - protected: - // Protected data - - // Model coefficients - - dimensionedScalar alphaK1_; - dimensionedScalar alphaK2_; - - dimensionedScalar alphaOmega1_; - dimensionedScalar alphaOmega2_; - - dimensionedScalar gamma1_; - dimensionedScalar gamma2_; - - dimensionedScalar beta1_; - dimensionedScalar beta2_; - - dimensionedScalar betaStar_; - - dimensionedScalar a1_; - dimensionedScalar b1_; - dimensionedScalar c1_; - - Switch F3_; - - - // Fields - - //- Wall distance - // Note: different to wall distance in parent RASModel - // which is for near-wall cells only - const volScalarField& y_; - - volScalarField k_; - volScalarField omega_; - - - // Private Member Functions - - tmp F1(const volScalarField& CDkOmega) const; - tmp F2() const; - tmp F3() const; - tmp F23() const; - - tmp blend - ( - const volScalarField& F1, - const dimensionedScalar& psi1, - const dimensionedScalar& psi2 - ) const - { - return F1*(psi1 - psi2) + psi2; - } - - tmp alphaK(const volScalarField& F1) const - { - return blend(F1, alphaK1_, alphaK2_); - } - - tmp alphaOmega(const volScalarField& F1) const - { - return blend(F1, alphaOmega1_, alphaOmega2_); - } - - tmp beta(const volScalarField& F1) const - { - return blend(F1, beta1_, beta2_); - } - - tmp gamma(const volScalarField& F1) const - { - return blend(F1, gamma1_, gamma2_); - } - - void correctNut(const volScalarField& S2); - - - // Protected Member Functions + // Protecetd Member Functions + virtual void correctNut(const volScalarField& S2); virtual void correctNut(); - virtual tmp kSource() const; - virtual tmp omegaSource() const; - virtual tmp Qsas - ( - const volScalarField& S2, - const volScalarField& gamma, - const volScalarField& beta - ) const; public: @@ -250,68 +167,6 @@ public: //- Destructor virtual ~kOmegaSST() {} - - - // Member Functions - - //- Re-read model coefficients if they have changed - virtual bool read(); - - //- Return the effective diffusivity for k - tmp DkEff(const volScalarField& F1) const - { - return tmp - ( - new volScalarField("DkEff", alphaK(F1)*this->nut_ + this->nu()) - ); - } - - //- Return the effective diffusivity for omega - tmp DomegaEff(const volScalarField& F1) const - { - return tmp - ( - new volScalarField - ( - "DomegaEff", - alphaOmega(F1)*this->nut_ + this->nu() - ) - ); - } - - //- Return the turbulence kinetic energy - virtual tmp k() const - { - return k_; - } - - //- Return the turbulence kinetic energy dissipation rate - virtual tmp epsilon() const - { - return tmp - ( - new volScalarField - ( - IOobject - ( - "epsilon", - this->mesh_.time().timeName(), - this->mesh_ - ), - betaStar_*k_*omega_, - omega_.boundaryField().types() - ) - ); - } - - //- Return the turbulence kinetic energy dissipation rate - virtual tmp omega() const - { - return omega_; - } - - //- Solve the turbulence equations and correct the turbulence viscosity - virtual void correct(); }; diff --git a/src/TurbulenceModels/turbulenceModels/RAS/kOmegaSST/kOmegaSSTBase.C b/src/TurbulenceModels/turbulenceModels/RAS/kOmegaSST/kOmegaSSTBase.C new file mode 100644 index 0000000000..19943c0c1b --- /dev/null +++ b/src/TurbulenceModels/turbulenceModels/RAS/kOmegaSST/kOmegaSSTBase.C @@ -0,0 +1,480 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 . + +\*---------------------------------------------------------------------------*/ + +#include "kOmegaSSTBase.H" +#include "bound.H" +#include "wallDist.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +// * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * // + +template +tmp kOmegaSSTBase::F1 +( + const volScalarField& CDkOmega +) const +{ + tmp CDkOmegaPlus = max + ( + CDkOmega, + dimensionedScalar("1.0e-10", dimless/sqr(dimTime), 1.0e-10) + ); + + tmp arg1 = min + ( + min + ( + max + ( + (scalar(1)/betaStar_)*sqrt(k_)/(omega_*y_), + scalar(500)*(this->mu()/this->rho_)/(sqr(y_)*omega_) + ), + (4*alphaOmega2_)*k_/(CDkOmegaPlus*sqr(y_)) + ), + scalar(10) + ); + + return tanh(pow4(arg1)); +} + + +template +tmp kOmegaSSTBase::F2() const +{ + tmp arg2 = min + ( + max + ( + (scalar(2)/betaStar_)*sqrt(k_)/(omega_*y_), + scalar(500)*(this->mu()/this->rho_)/(sqr(y_)*omega_) + ), + scalar(100) + ); + + return tanh(sqr(arg2)); +} + + +template +tmp kOmegaSSTBase::F3() const +{ + tmp arg3 = min + ( + 150*(this->mu()/this->rho_)/(omega_*sqr(y_)), + scalar(10) + ); + + return 1 - tanh(pow4(arg3)); +} + + +template +tmp kOmegaSSTBase::F23() const +{ + tmp f23(F2()); + + if (F3_) + { + f23() *= F3(); + } + + return f23; +} + + +template +void kOmegaSSTBase::correctNut +( + const volScalarField& S2 +) +{ + // Correct the turbulence viscosity + this->nut_ = a1_*k_/max(a1_*omega_, b1_*F23()*sqrt(S2)); + this->nut_.correctBoundaryConditions(); +} + + +template +void kOmegaSSTBase::correctNut() +{ + correctNut(2*magSqr(symm(fvc::grad(this->U_)))); +} + + +template +tmp kOmegaSSTBase::kSource() const +{ + return tmp + ( + new fvScalarMatrix + ( + k_, + dimVolume*this->rho_.dimensions()*k_.dimensions()/dimTime + ) + ); +} + + +template +tmp kOmegaSSTBase::omegaSource() const +{ + return tmp + ( + new fvScalarMatrix + ( + omega_, + dimVolume*this->rho_.dimensions()*omega_.dimensions()/dimTime + ) + ); +} + + +template +tmp kOmegaSSTBase::Qsas +( + const volScalarField& S2, + const volScalarField& gamma, + const volScalarField& beta +) const +{ + return tmp + ( + new fvScalarMatrix + ( + omega_, + dimVolume*this->rho_.dimensions()*omega_.dimensions()/dimTime + ) + ); +} + + +// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * // + +template +kOmegaSSTBase::kOmegaSSTBase +( + const word& type, + const alphaField& alpha, + const rhoField& rho, + const volVectorField& U, + const surfaceScalarField& alphaRhoPhi, + const surfaceScalarField& phi, + const transportModel& transport, + const word& propertiesName +) +: + BasicEddyViscosityModel + ( + type, + alpha, + rho, + U, + alphaRhoPhi, + phi, + transport, + propertiesName + ), + + alphaK1_ + ( + dimensioned::lookupOrAddToDict + ( + "alphaK1", + this->coeffDict_, + 0.85 + ) + ), + alphaK2_ + ( + dimensioned::lookupOrAddToDict + ( + "alphaK2", + this->coeffDict_, + 1.0 + ) + ), + alphaOmega1_ + ( + dimensioned::lookupOrAddToDict + ( + "alphaOmega1", + this->coeffDict_, + 0.5 + ) + ), + alphaOmega2_ + ( + dimensioned::lookupOrAddToDict + ( + "alphaOmega2", + this->coeffDict_, + 0.856 + ) + ), + gamma1_ + ( + dimensioned::lookupOrAddToDict + ( + "gamma1", + this->coeffDict_, + 5.0/9.0 + ) + ), + gamma2_ + ( + dimensioned::lookupOrAddToDict + ( + "gamma2", + this->coeffDict_, + 0.44 + ) + ), + beta1_ + ( + dimensioned::lookupOrAddToDict + ( + "beta1", + this->coeffDict_, + 0.075 + ) + ), + beta2_ + ( + dimensioned::lookupOrAddToDict + ( + "beta2", + this->coeffDict_, + 0.0828 + ) + ), + betaStar_ + ( + dimensioned::lookupOrAddToDict + ( + "betaStar", + this->coeffDict_, + 0.09 + ) + ), + a1_ + ( + dimensioned::lookupOrAddToDict + ( + "a1", + this->coeffDict_, + 0.31 + ) + ), + b1_ + ( + dimensioned::lookupOrAddToDict + ( + "b1", + this->coeffDict_, + 1.0 + ) + ), + c1_ + ( + dimensioned::lookupOrAddToDict + ( + "c1", + this->coeffDict_, + 10.0 + ) + ), + F3_ + ( + Switch::lookupOrAddToDict + ( + "F3", + this->coeffDict_, + false + ) + ), + + y_(wallDist::New(this->mesh_).y()), + + k_ + ( + IOobject + ( + IOobject::groupName("k", U.group()), + this->runTime_.timeName(), + this->mesh_, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + this->mesh_ + ), + omega_ + ( + IOobject + ( + IOobject::groupName("omega", U.group()), + this->runTime_.timeName(), + this->mesh_, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + this->mesh_ + ) +{ + bound(k_, this->kMin_); + bound(omega_, this->omegaMin_); +} + + +// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * // + +template +bool kOmegaSSTBase::read() +{ + if (BasicEddyViscosityModel::read()) + { + alphaK1_.readIfPresent(this->coeffDict()); + alphaK2_.readIfPresent(this->coeffDict()); + alphaOmega1_.readIfPresent(this->coeffDict()); + alphaOmega2_.readIfPresent(this->coeffDict()); + gamma1_.readIfPresent(this->coeffDict()); + gamma2_.readIfPresent(this->coeffDict()); + beta1_.readIfPresent(this->coeffDict()); + beta2_.readIfPresent(this->coeffDict()); + betaStar_.readIfPresent(this->coeffDict()); + a1_.readIfPresent(this->coeffDict()); + b1_.readIfPresent(this->coeffDict()); + c1_.readIfPresent(this->coeffDict()); + F3_.readIfPresent("F3", this->coeffDict()); + + return true; + } + else + { + return false; + } +} + + +template +void kOmegaSSTBase::correct() +{ + if (!this->turbulence_) + { + return; + } + + // Local references + const alphaField& alpha = this->alpha_; + const rhoField& rho = this->rho_; + const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_; + const volVectorField& U = this->U_; + volScalarField& nut = this->nut_; + + BasicEddyViscosityModel::correct(); + + volScalarField divU(fvc::div(fvc::absolute(this->phi(), U))); + + tmp tgradU = fvc::grad(U); + volScalarField S2(2*magSqr(symm(tgradU()))); + volScalarField GbyNu((tgradU() && dev(twoSymm(tgradU())))); + volScalarField G(this->GName(), nut*GbyNu); + tgradU.clear(); + + // Update omega and G at the wall + omega_.boundaryField().updateCoeffs(); + + volScalarField CDkOmega + ( + (2*alphaOmega2_)*(fvc::grad(k_) & fvc::grad(omega_))/omega_ + ); + + volScalarField F1(this->F1(CDkOmega)); + + { + volScalarField gamma(this->gamma(F1)); + volScalarField beta(this->beta(F1)); + + // Turbulent frequency equation + tmp omegaEqn + ( + fvm::ddt(alpha, rho, omega_) + + fvm::div(alphaRhoPhi, omega_) + - fvm::laplacian(alpha*rho*DomegaEff(F1), omega_) + == + alpha*rho*gamma + *min + ( + GbyNu, + (c1_/a1_)*betaStar_*omega_*max(a1_*omega_, b1_*F23()*sqrt(S2)) + ) + - fvm::SuSp((2.0/3.0)*alpha*rho*gamma*divU, omega_) + - fvm::Sp(alpha*rho*beta*omega_, omega_) + - fvm::SuSp + ( + alpha*rho*(F1 - scalar(1))*CDkOmega/omega_, + omega_ + ) + + Qsas(S2, gamma, beta) + + omegaSource() + ); + + omegaEqn().relax(); + + omegaEqn().boundaryManipulate(omega_.boundaryField()); + + solve(omegaEqn); + bound(omega_, this->omegaMin_); + } + + // Turbulent kinetic energy equation + tmp kEqn + ( + fvm::ddt(alpha, rho, k_) + + fvm::div(alphaRhoPhi, k_) + - fvm::laplacian(alpha*rho*DkEff(F1), k_) + == + min(alpha*rho*G, (c1_*betaStar_)*alpha*rho*k_*omega_) + - fvm::SuSp((2.0/3.0)*alpha*rho*divU, k_) + - fvm::Sp(alpha*rho*betaStar_*omega_, k_) + + kSource() + ); + + kEqn().relax(); + solve(kEqn); + bound(k_, this->kMin_); + + correctNut(S2); +} + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// ************************************************************************* // diff --git a/src/TurbulenceModels/turbulenceModels/RAS/kOmegaSST/kOmegaSSTBase.H b/src/TurbulenceModels/turbulenceModels/RAS/kOmegaSST/kOmegaSSTBase.H new file mode 100644 index 0000000000..eaa4bd7652 --- /dev/null +++ b/src/TurbulenceModels/turbulenceModels/RAS/kOmegaSST/kOmegaSSTBase.H @@ -0,0 +1,320 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 . + +Class + Foam::kOmegaSSTBase + +Description + Base class implementation of the k-omega-SST turbulence model for + incompressible and compressible flows. + + Turbulence model described in: + \verbatim + Menter, F. R. & Esch, T. (2001). + Elements of Industrial Heat Transfer Prediction. + 16th Brazilian Congress of Mechanical Engineering (COBEM). + \endverbatim + + with updated coefficients from + \verbatim + Menter, F. R., Kuntz, M., and Langtry, R. (2003). + Ten Years of Industrial Experience with the SST Turbulence Model. + Turbulence, Heat and Mass Transfer 4, ed: K. Hanjalic, Y. Nagano, + & M. Tummers, Begell House, Inc., 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. (1998). + "Some Improvements in Menter’s k-omega-SST turbulence model" + 29th AIAA Fluid Dynamics Conference, AIAA-98-2554. + \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 + manner. The paper suggests that sigma is blended but this would not be + consistent with the blending of the k-epsilon and k-omega models. + + 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 if y+ becomes + sufficiently small blending u_tau in this manner clearly becomes nonsense. + + The default model coefficients are + \verbatim + kOmegaSSTBaseCoeffs + { + alphaK1 0.85; + alphaK2 1.0; + alphaOmega1 0.5; + alphaOmega2 0.856; + beta1 0.075; + beta2 0.0828; + betaStar 0.09; + gamma1 5/9; + gamma2 0.44; + a1 0.31; + b1 1.0; + c1 10.0; + F3 no; + } + \endverbatim + +SourceFiles + kOmegaSSTBase.C + +\*---------------------------------------------------------------------------*/ + +#ifndef kOmegaSSTBase_H +#define kOmegaSSTBase_H + +#include "RASModel.H" +#include "eddyViscosity.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +namespace Foam +{ + +/*---------------------------------------------------------------------------*\ + Class kOmegaSSTBase Declaration +\*---------------------------------------------------------------------------*/ + +template +class kOmegaSSTBase +: + public BasicEddyViscosityModel +{ + // Private Member Functions + + // Disallow default bitwise copy construct and assignment + kOmegaSSTBase(const kOmegaSSTBase&); + kOmegaSSTBase& operator=(const kOmegaSSTBase&); + + +protected: + + // Protected data + + // Model coefficients + + dimensionedScalar alphaK1_; + dimensionedScalar alphaK2_; + + dimensionedScalar alphaOmega1_; + dimensionedScalar alphaOmega2_; + + dimensionedScalar gamma1_; + dimensionedScalar gamma2_; + + dimensionedScalar beta1_; + dimensionedScalar beta2_; + + dimensionedScalar betaStar_; + + dimensionedScalar a1_; + dimensionedScalar b1_; + dimensionedScalar c1_; + + Switch F3_; + + + // Fields + + //- Wall distance + // Note: different to wall distance in parent RASModel + // which is for near-wall cells only + const volScalarField& y_; + + volScalarField k_; + volScalarField omega_; + + + // Protected Member Functions + + tmp F1(const volScalarField& CDkOmega) const; + tmp F2() const; + tmp F3() const; + tmp F23() const; + + tmp blend + ( + const volScalarField& F1, + const dimensionedScalar& psi1, + const dimensionedScalar& psi2 + ) const + { + return F1*(psi1 - psi2) + psi2; + } + + tmp alphaK(const volScalarField& F1) const + { + return blend(F1, alphaK1_, alphaK2_); + } + + tmp alphaOmega(const volScalarField& F1) const + { + return blend(F1, alphaOmega1_, alphaOmega2_); + } + + tmp beta(const volScalarField& F1) const + { + return blend(F1, beta1_, beta2_); + } + + tmp gamma(const volScalarField& F1) const + { + return blend(F1, gamma1_, gamma2_); + } + + void correctNut(const volScalarField& S2); + + virtual void correctNut(); + virtual tmp kSource() const; + virtual tmp omegaSource() const; + virtual tmp Qsas + ( + const volScalarField& S2, + const volScalarField& gamma, + const volScalarField& beta + ) const; + + +public: + + typedef typename BasicEddyViscosityModel::alphaField alphaField; + typedef typename BasicEddyViscosityModel::rhoField rhoField; + typedef typename BasicEddyViscosityModel::transportModel transportModel; + + + // Constructors + + //- Construct from components + kOmegaSSTBase + ( + const word& type, + const alphaField& alpha, + const rhoField& rho, + const volVectorField& U, + const surfaceScalarField& alphaRhoPhi, + const surfaceScalarField& phi, + const transportModel& transport, + const word& propertiesName = turbulenceModel::propertiesName + ); + + + //- Destructor + virtual ~kOmegaSSTBase() + {} + + + // Member Functions + + //- Re-read model coefficients if they have changed + virtual bool read(); + + //- Return the effective diffusivity for k + tmp DkEff(const volScalarField& F1) const + { + return tmp + ( + new volScalarField("DkEff", alphaK(F1)*this->nut_ + this->nu()) + ); + } + + //- Return the effective diffusivity for omega + tmp DomegaEff(const volScalarField& F1) const + { + return tmp + ( + new volScalarField + ( + "DomegaEff", + alphaOmega(F1)*this->nut_ + this->nu() + ) + ); + } + + //- Return the turbulence kinetic energy + virtual tmp k() const + { + return k_; + } + + //- Return the turbulence kinetic energy dissipation rate + virtual tmp epsilon() const + { + return tmp + ( + new volScalarField + ( + IOobject + ( + "epsilon", + this->mesh_.time().timeName(), + this->mesh_ + ), + betaStar_*k_*omega_, + omega_.boundaryField().types() + ) + ); + } + + //- Return the turbulence kinetic energy dissipation rate + virtual tmp omega() const + { + return omega_; + } + + //- Solve the turbulence equations and correct the turbulence viscosity + virtual void correct(); +}; + + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +} // End namespace Foam + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#ifdef NoRepository +# include "kOmegaSSTBase.C" +#endif + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +#endif + +// ************************************************************************* // diff --git a/src/TurbulenceModels/turbulenceModels/turbulenceModelDoc.H b/src/TurbulenceModels/turbulenceModels/doc/turbulenceModelDoc.H similarity index 100% rename from src/TurbulenceModels/turbulenceModels/turbulenceModelDoc.H rename to src/TurbulenceModels/turbulenceModels/doc/turbulenceModelDoc.H