Commit 6d380ddd authored by Andrew Heather's avatar Andrew Heather
Browse files

ENH: Updated kOmegaSST family of models

parent 3dbd3914
......@@ -142,14 +142,31 @@ template<class BasicEddyViscosityModel>
tmp<volScalarField::Internal>
kOmegaSSTBase<BasicEddyViscosityModel>::epsilonByk
(
const volScalarField::Internal& F1,
const volScalarField::Internal& F2
const volScalarField& F1,
const volTensorField& gradU
) const
{
return betaStar_*omega_();
}
template<class BasicEddyViscosityModel>
tmp<volScalarField::Internal> kOmegaSSTBase<BasicEddyViscosityModel>::GbyNu
(
const volScalarField::Internal& GbyNu0,
const volScalarField::Internal& F2,
const volScalarField::Internal& S2
) const
{
return min
(
GbyNu0,
(c1_/a1_)*betaStar_*omega_()
*max(a1_*omega_(), b1_*F2*sqrt(S2))
);
}
template<class BasicEddyViscosityModel>
tmp<fvScalarMatrix> kOmegaSSTBase<BasicEddyViscosityModel>::kSource() const
{
......@@ -181,9 +198,9 @@ tmp<fvScalarMatrix> kOmegaSSTBase<BasicEddyViscosityModel>::omegaSource() const
template<class BasicEddyViscosityModel>
tmp<fvScalarMatrix> kOmegaSSTBase<BasicEddyViscosityModel>::Qsas
(
const volScalarField& S2,
const volScalarField& gamma,
const volScalarField& beta
const volScalarField::Internal& S2,
const volScalarField::Internal& gamma,
const volScalarField::Internal& beta
) const
{
return tmp<fvScalarMatrix>
......@@ -422,13 +439,12 @@ void kOmegaSSTBase<BasicEddyViscosityModel>::correct()
BasicEddyViscosityModel::correct();
volScalarField divU(fvc::div(fvc::absolute(this->phi(), U)));
volScalarField::Internal divU(fvc::div(fvc::absolute(this->phi(), U)));
tmp<volTensorField> tgradU = fvc::grad(U);
volScalarField S2(2*magSqr(symm(tgradU())));
volScalarField GbyNu((tgradU() && dev(twoSymm(tgradU()))));
volScalarField G(this->GName(), nut*GbyNu);
tgradU.clear();
volScalarField::Internal GbyNu0((tgradU() && dev(twoSymm(tgradU()))));
volScalarField::Internal G(this->GName(), nut*GbyNu0);
// Update omega and G at the wall
omega_.boundaryFieldRef().updateCoeffs();
......@@ -439,10 +455,11 @@ void kOmegaSSTBase<BasicEddyViscosityModel>::correct()
);
volScalarField F1(this->F1(CDkOmega));
volScalarField F23(this->F23());
{
volScalarField gamma(this->gamma(F1));
volScalarField beta(this->beta(F1));
volScalarField::Internal gamma(this->gamma(F1));
volScalarField::Internal beta(this->beta(F1));
// Turbulent frequency equation
tmp<fvScalarMatrix> omegaEqn
......@@ -451,20 +468,15 @@ void kOmegaSSTBase<BasicEddyViscosityModel>::correct()
+ 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_)
alpha*rho*gamma*GbyNu(GbyNu0, F23(), 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_,
alpha()*rho()*(F1() - scalar(1))*CDkOmega()/omega_(),
omega_
)
+ Qsas(S2, gamma, beta)
+ Qsas(S2(), gamma, beta)
+ omegaSource()
+ fvOptions(alpha, rho, omega_)
);
......@@ -484,13 +496,15 @@ void kOmegaSSTBase<BasicEddyViscosityModel>::correct()
+ 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_)
alpha()*rho()*Pk(G)
- fvm::SuSp((2.0/3.0)*alpha()*rho()*divU, k_)
- fvm::Sp(alpha()*rho()*epsilonByk(F1, tgradU()), k_)
+ kSource()
+ fvOptions(alpha, rho, k_)
);
tgradU.clear();
kEqn.ref().relax();
fvOptions.constrain(kEqn.ref());
solve(kEqn);
......
......@@ -100,9 +100,6 @@ SourceFiles
#ifndef kOmegaSSTBase_H
#define kOmegaSSTBase_H
#include "RASModel.H"
#include "eddyViscosity.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
......@@ -121,7 +118,7 @@ class kOmegaSSTBase
// Disallow default bitwise copy construct and assignment
kOmegaSSTBase(const kOmegaSSTBase&);
kOmegaSSTBase& operator=(const kOmegaSSTBase&);
void operator=(const kOmegaSSTBase&);
protected:
......@@ -164,10 +161,10 @@ protected:
// Protected Member Functions
tmp<volScalarField> F1(const volScalarField& CDkOmega) const;
tmp<volScalarField> F2() const;
tmp<volScalarField> F3() const;
tmp<volScalarField> F23() const;
virtual tmp<volScalarField> F1(const volScalarField& CDkOmega) const;
virtual tmp<volScalarField> F2() const;
virtual tmp<volScalarField> F3() const;
virtual tmp<volScalarField> F23() const;
tmp<volScalarField> blend
(
......@@ -179,6 +176,16 @@ protected:
return F1*(psi1 - psi2) + psi2;
}
tmp<volScalarField::Internal> blend
(
const volScalarField::Internal& F1,
const dimensionedScalar& psi1,
const dimensionedScalar& psi2
) const
{
return F1*(psi1 - psi2) + psi2;
}
tmp<volScalarField> alphaK(const volScalarField& F1) const
{
return blend(F1, alphaK1_, alphaK2_);
......@@ -189,12 +196,18 @@ protected:
return blend(F1, alphaOmega1_, alphaOmega2_);
}
tmp<volScalarField> beta(const volScalarField& F1) const
tmp<volScalarField::Internal> beta
(
const volScalarField::Internal& F1
) const
{
return blend(F1, beta1_, beta2_);
}
tmp<volScalarField> gamma(const volScalarField& F1) const
tmp<volScalarField::Internal> gamma
(
const volScalarField::Internal& F1
) const
{
return blend(F1, gamma1_, gamma2_);
}
......@@ -212,8 +225,16 @@ protected:
//- Return epsilon/k which for standard RAS is betaStar*omega
virtual tmp<volScalarField::Internal> epsilonByk
(
const volScalarField::Internal& F1,
const volScalarField::Internal& F2
const volScalarField& F1,
const volTensorField& gradU
) const;
//- Return G/nu
virtual tmp<volScalarField::Internal> GbyNu
(
const volScalarField::Internal& GbyNu0,
const volScalarField::Internal& F2,
const volScalarField::Internal& S2
) const;
virtual tmp<fvScalarMatrix> kSource() const;
......@@ -222,9 +243,9 @@ protected:
virtual tmp<fvScalarMatrix> Qsas
(
const volScalarField& S2,
const volScalarField& gamma,
const volScalarField& beta
const volScalarField::Internal& S2,
const volScalarField::Internal& gamma,
const volScalarField::Internal& beta
) const;
......@@ -326,7 +347,7 @@ public:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "kOmegaSSTBase.C"
#include "kOmegaSSTBase.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
......
......@@ -72,7 +72,7 @@ class kOmegaSSTDDES
// Disallow default bitwise copy construct and assignment
kOmegaSSTDDES(const kOmegaSSTDDES&);
kOmegaSSTDDES& operator=(const kOmegaSSTDDES&);
void operator=(const kOmegaSSTDDES&);
protected:
......
......@@ -66,6 +66,30 @@ tmp<volScalarField> kOmegaSSTDES<BasicTurbulenceModel>::dTilda
}
template<class BasicTurbulenceModel>
tmp<volScalarField::Internal> kOmegaSSTDES<BasicTurbulenceModel>::epsilonByk
(
const volScalarField& F1,
const volTensorField& gradU
) const
{
volScalarField CDES(this->CDES(F1));
return sqrt(this->k_())/dTilda(mag(gradU), CDES)()();
}
template<class BasicTurbulenceModel>
tmp<volScalarField::Internal> kOmegaSSTDES<BasicTurbulenceModel>::GbyNu
(
const volScalarField::Internal& GbyNu0,
const volScalarField::Internal& F2,
const volScalarField::Internal& S2
) const
{
return GbyNu0; // Unlimited
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class BasicTurbulenceModel>
......@@ -150,96 +174,6 @@ bool kOmegaSSTDES<BasicTurbulenceModel>::read()
}
template<class BasicTurbulenceModel>
void kOmegaSSTDES<BasicTurbulenceModel>::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& k = this->k_;
volScalarField& omega = this->omega_;
volScalarField& nut = this->nut_;
DESModel<BasicTurbulenceModel>::correct();
volScalarField divU(fvc::div(fvc::absolute(this->phi(), U)));
tmp<volTensorField> tgradU = fvc::grad(U);
volScalarField magGradU(mag(tgradU()));
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.boundaryFieldRef().updateCoeffs();
volScalarField CDkOmega
(
(2*this->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<fvScalarMatrix> omegaEqn
(
fvm::ddt(alpha, rho, omega)
+ fvm::div(alphaRhoPhi, omega)
- fvm::laplacian(alpha*rho*this->DomegaEff(F1), omega)
==
alpha*rho*gamma*GbyNu // Using unlimited GybNu
- 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)
+ this->omegaSource()
);
omegaEqn.ref().relax();
omegaEqn.ref().boundaryManipulate(omega.boundaryFieldRef());
solve(omegaEqn);
bound(omega, this->omegaMin_);
}
{
volScalarField CDES(this->CDES(F1));
volScalarField dTilda(this->dTilda(magGradU, CDES));
// Turbulent kinetic energy equation
tmp<fvScalarMatrix> kEqn
(
fvm::ddt(alpha, rho, k)
+ fvm::div(alphaRhoPhi, k)
- fvm::laplacian(alpha*rho*this->DkEff(F1), k)
==
min(alpha*rho*G, (this->c1_*this->betaStar_)*alpha*rho*k*omega)
- fvm::SuSp((2.0/3.0)*alpha*rho*divU, k)
- fvm::Sp(alpha*rho*sqrt(k)/dTilda, k) // modified for DES
+ this->kSource()
);
kEqn.ref().relax();
solve(kEqn);
bound(k, this->kMin_);
}
this->correctNut(S2);
}
template<class BasicTurbulenceModel>
tmp<volScalarField> kOmegaSSTDES<BasicTurbulenceModel>::LESRegion() const
{
......
......@@ -74,7 +74,7 @@ class kOmegaSSTDES
// Disallow default bitwise copy construct and assignment
kOmegaSSTDES(const kOmegaSSTDES&);
kOmegaSSTDES& operator=(const kOmegaSSTDES&);
void operator=(const kOmegaSSTDES&);
protected:
......@@ -106,6 +106,21 @@ protected:
const volScalarField& CDES
) const;
//- Return epsilon/k
virtual tmp<volScalarField::Internal> epsilonByk
(
const volScalarField& F1,
const volTensorField& gradU
) const;
//- Return G/nu
virtual tmp<volScalarField::Internal> GbyNu
(
const volScalarField::Internal& GbyNu0,
const volScalarField::Internal& F2,
const volScalarField::Internal& S2
) const;
public:
......@@ -144,9 +159,6 @@ public:
//- Re-read model coefficients if they have changed
virtual bool read();
//- Solve the turbulence equations and correct the turbulence viscosity
virtual void correct();
//- Return the LES field indicator
virtual tmp<volScalarField> LESRegion() const;
};
......@@ -159,7 +171,7 @@ public:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "kOmegaSSTDES.C"
#include "kOmegaSSTDES.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
......
......@@ -85,7 +85,7 @@ class kOmegaSSTIDDES
// Disallow default bitwise copy construct and assignment
kOmegaSSTIDDES(const kOmegaSSTIDDES&);
kOmegaSSTIDDES& operator=(const kOmegaSSTIDDES&);
void operator=(const kOmegaSSTIDDES&);
protected:
......@@ -160,7 +160,7 @@ public:
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#ifdef NoRepository
# include "kOmegaSSTIDDES.C"
#include "kOmegaSSTIDDES.C"
#endif
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
......
......@@ -60,13 +60,13 @@ tmp<volScalarField::Internal> kOmegaSSTLM<BasicTurbulenceModel>::Pk
template<class BasicTurbulenceModel>
tmp<volScalarField::Internal> kOmegaSSTLM<BasicTurbulenceModel>::epsilonByk
(
const volScalarField::Internal& F1,
const volScalarField::Internal& F2
const volScalarField& F1,
const volTensorField& gradU
) const
{
return
min(max(gammaIntEff_, scalar(0.1)), scalar(1))
*kOmegaSST<BasicTurbulenceModel>::epsilonByk(F1, F2);
*kOmegaSST<BasicTurbulenceModel>::epsilonByk(F1, gradU);
}
......
......@@ -164,8 +164,8 @@ protected:
//- Modified form of the k-omega SST epsilon/k
virtual tmp<volScalarField::Internal> epsilonByk
(
const volScalarField::Internal& F1,
const volScalarField::Internal& F2
const volScalarField& F1,
const volTensorField& gradU
) const;
//- Freestream blending-function
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment