CFD Online Discussion Forums - Some confusions about SST-DES (code and formulas)

Alisa

Short answer: Check

Code:

$FOAM_SRC/turbulenceModels/Base/kOmegaSST/kOmegaSSTBase.C

Follow below for long answer.

Find the constructor in the code you attached:

Code:

template<class BasicTurbulenceModel>

kOmegaSSTDES<BasicTurbulenceModel>::kOmegaSSTDES

(

    someParametersHere

)

    kOmegaSST

    <

        LESeddyViscosity<BasicTurbulenceModel>,

        BasicTurbulenceModel

    >

This basically means that the class kOmegaSSTDES is a sub-class of the class kOmegaSST. Parameters in <> are parameters of the general class kOmegaSST.

You can find the constructor implementation of the template class in $FOAM_SRC/TurbulenceModels/ turbulenceModels/*

So now it is crucial to figure out how kOmegaSST is implemented. That information is in

Code:

$FOAM_SRC/turbulenceModels/RAS/kOmegaSST/kOmegaSST.C

Code:

template<class BasicTurbulenceModel>

kOmegaSST<BasicTurbulenceModel>::kOmegaSST

(

    some parameters

)

:

    Foam::kOmegaSST

    <

        eddyViscosity<RASModel<BasicTurbulenceModel>>,

        BasicTurbulenceModel

    >

This file in fact doesn't provide much information except that it is an eddyViscosity model. However when you check the .H file, you will see that it includes kOmegaSSTBase.H which involves whole mathematical representation of the k-omega SST model.

This is the mathematical model for k-omega SST:

${{\partial k} \over {\partial t}} + U_j {{\partial k} \over {\partial x_j }} = P_k - \beta ^* k\omega + {\partial \over {\partial x_j }}\left[ {\left( {\nu + \sigma_k \nu _T } \right){{\partial k} \over {\partial x_j }}} \right]$
${{\partial \omega } \over {\partial t}} + U_j {{\partial \omega } \over {\partial x_j }} = \alpha S^2 - \beta \omega ^2 + {\partial \over {\partial x_j }}\left[ {\left( {\nu + \sigma_{\omega} \nu _T } \right){{\partial \omega } \over {\partial x_j }}} \right] + 2( 1 - F_1 ) \sigma_{\omega 2} {1 \over \omega} {{\partial k } \over {\partial x_i}} {{\partial \omega } \over {\partial x_i}}$

$F_2=\mbox{tanh} \left[ \left[ \mbox{max} \left( { 2 \sqrt{k} \over \beta^* \omega y } , { 500 \nu \over y^2 \omega } \right) \right]^2 \right]$
$P_k=\mbox{min} \left(\tau _{ij} {{\partial U_i } \over {\partial x_j }} , 10\beta^* k \omega \right)$
$F_1=\mbox{tanh} \left\{ \left\{ \mbox{min} \left[ \mbox{max} \left( {\sqrt{k} \over \beta ^* \omega y}, {500 \nu \over y^2 \omega} \right) , {4 \sigma_{\omega 2} k \over CD_{k\omega} y^2} \right] \right\} ^4 \right\}$
$CD_{k\omega}=\mbox{max} \left( 2\rho\sigma_{\omega 2} {1 \over \omega} {{\partial k} \over {\partial x_i}} {{\partial \omega} \over {\partial x_i}}, 10 ^{-10} \right )$

So let's start with turbulence kinetic energy:

Code:

    // Turbulent kinetic energy equation

    tmp<fvScalarMatrix> kEqn

    (

        fvm::ddt(alpha, rho, k_)

      + fvm::div(alphaRhoPhi, k_)

      - fvm::laplacian(alpha*rho*DkEff(F1), k_)

     ==

        alpha()*rho()*Pk(G)

      - fvm::SuSp((2.0/3.0)*alpha()*rho()*divU, k_)

      - fvm::Sp(alpha()*rho()*epsilonByk(F1, F23), k_)

      + kSource()

      + fvOptions(alpha, rho, k_)

    );

I assume that first two terms (time derivation and convection) are clear. To understand the third term, one needs to clarify DkEff(F1).

F1 is the blending factor which is in charge of the transition between boundary layer (=1) and far field (=0). This is implemented as a protected member function and code below clearly corresponds to the F1 formula written above:

Code:

    tmp<volScalarField> CDkOmegaPlus = max

    (

        CDkOmega,

        dimensionedScalar("1.0e-10", dimless/sqr(dimTime), 1.0e-10)

    );



    tmp<volScalarField> 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));

where

Code:

    volScalarField CDkOmega

    (

        (2*alphaOmega2_)*(fvc::grad(k_) & fvc::grad(omega_))/omega_

    );

So what does DkEff(F1), also called effective diffusivity do? See kOmegaSSTBase.H:

Code:

        //- Return the effective diffusivity for k

        tmp<volScalarField> DkEff(const volScalarField& F1) const

        {

            return tmp<volScalarField>

            (

                new volScalarField("DkEff", alphaK(F1)*this->nut_ + this->nu())

            );

        }

I think this one looks very much like the term $\left( {\nu + \sigma_k \nu _T } \right)$ . Except we need to figure out alphaK function:

Code:

        tmp<volScalarField> alphaK(const volScalarField& F1) const

        {

            return blend(F1, alphaK1_, alphaK2_);

        }

So it calls blending function. Nothing fancy. Keep in mind that alpha = 1/sigma by the way. This is what blend function does:

Code:

        tmp<volScalarField::Internal> blend

        (

            const volScalarField::Internal& F1,

            const dimensionedScalar& psi1,

            const dimensionedScalar& psi2

        ) const

        {

            return F1*(psi1 - psi2) + psi2;

        }

where Psi is: $\phi = \phi_1 F_1 + \phi_2 (1 - F_1)$ . Looks like terms are re-arranged but code still represents the correct formula.

So now we should take a look at the definition of this->nut_. It is called kinematic eddy viscosity and is defined as:

$\nu _T = {a_1 k \over \mbox{max}(a_1 \omega, S F_2) }$

and the corresponding code is:

Code:

this->nut_ = a1_*k_/max(a1_*omega_, b1_*F2*sqrt(S2));

F2 is given above. See its code:

Code:

    tmp<volScalarField> 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));

and S2 is:

Code:

    

    volScalarField S2(2*magSqr(symm(tgradU())));

where $S = \left | \frac{1}{2} \left( \partial_j u_i + \partial_i u_j \right) \right |$ and S = sqrt(S2).

I think this explains the third term. In general what you see in the code is hidden behind the DkEff function. When you dig deeper it becomes clearer.

The first term on the RHS P_k is also explicitly implemented. You can find in the same file:

Code:

$FOAM_SRC/turbulenceModels/Base/kOmegaSST/kOmegaSSTBase.C

where return value for P_k is

Code:

return min(G, (c1_*betaStar_)*this->k_()*this->omega_());

where G is

Code:

    volScalarField::Internal GbyNu(dev(twoSymm(tgradU()())) && tgradU()());

    volScalarField::Internal G(this->GName(), nut()*GbyNu);

This is the mathematical definition of the G term: $G = \nu_t \frac{\partial u_i}{\partial x_j} \left( \frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i} \right)$

Now remember the definition of the tgradU that is written above. I don't have enough energy and time to look at the details of the GbyNu but it corresponds to the shear tensor.

So the final term $- \beta^* k \omega$ :

Code:

      - fvm::SuSp((2.0/3.0)*alpha()*rho()*divU, k_)

      - fvm::Sp(alpha()*rho()*epsilonByk(F1, F23), k_)

These two functions, SuSp and Sp are two methods of implicit source term implementation in OpenFOAM. You can refer to the Users Guide to understand their mechanism.

See epsilonByk function:

Code:

return betaStar_*omega_();

I hope this clarifies implementation of k-omega SST model a bit more for you. To understand DES expansion you need to figure out how length scale is implemented which is I guess is more clear.

References:
+ https://pingpong.chalmers.se/public/...o?item=3256524
+ https://pingpong.chalmers.se/public/...OmegaSST-1.pdf

Dear Fertinaz,
Thank you so much for your detailed explanation and I am also sorry for my late reply.

After reading your reply, I have really known more about the frame of OF, which contains so much inheritance and deriving that I sometimes can not find the root of certain variable. Of course, more important reason is that my knowledge on CFD is very poor.

Again, thank you so much. I greatly appreciate your generosity and patience.

Quote:

Originally Posted by fertinaz (Post 689794)

Alisa

Short answer: Check

Code:

$FOAM_SRC/turbulenceModels/Base/kOmegaSST/kOmegaSSTBase.C

Follow below for long answer.

Find the constructor in the code you attached:

Code:

template<class BasicTurbulenceModel>

kOmegaSSTDES<BasicTurbulenceModel>::kOmegaSSTDES

(

    someParametersHere

)

    kOmegaSST

    <

        LESeddyViscosity<BasicTurbulenceModel>,

        BasicTurbulenceModel

    >

Code:

$FOAM_SRC/turbulenceModels/RAS/kOmegaSST/kOmegaSST.C

Code:

template<class BasicTurbulenceModel>

kOmegaSST<BasicTurbulenceModel>::kOmegaSST

(

    some parameters

)

:

    Foam::kOmegaSST

    <

        eddyViscosity<RASModel<BasicTurbulenceModel>>,

        BasicTurbulenceModel

    >

${{\partial k} \over {\partial t}} + U_j {{\partial k} \over {\partial x_j }} = P_k - \beta ^* k\omega + {\partial \over {\partial x_j }}\left[ {\left( {\nu + \sigma_k \nu _T } \right){{\partial k} \over {\partial x_j }}} \right]$
${{\partial \omega } \over {\partial t}} + U_j {{\partial \omega } \over {\partial x_j }} = \alpha S^2 - \beta \omega ^2 + {\partial \over {\partial x_j }}\left[ {\left( {\nu + \sigma_{\omega} \nu _T } \right){{\partial \omega } \over {\partial x_j }}} \right] + 2( 1 - F_1 ) \sigma_{\omega 2} {1 \over \omega} {{\partial k } \over {\partial x_i}} {{\partial \omega } \over {\partial x_i}}$

$F_2=\mbox{tanh} \left[ \left[ \mbox{max} \left( { 2 \sqrt{k} \over \beta^* \omega y } , { 500 \nu \over y^2 \omega } \right) \right]^2 \right]$
$P_k=\mbox{min} \left(\tau _{ij} {{\partial U_i } \over {\partial x_j }} , 10\beta^* k \omega \right)$
$F_1=\mbox{tanh} \left\{ \left\{ \mbox{min} \left[ \mbox{max} \left( {\sqrt{k} \over \beta ^* \omega y}, {500 \nu \over y^2 \omega} \right) , {4 \sigma_{\omega 2} k \over CD_{k\omega} y^2} \right] \right\} ^4 \right\}$
$CD_{k\omega}=\mbox{max} \left( 2\rho\sigma_{\omega 2} {1 \over \omega} {{\partial k} \over {\partial x_i}} {{\partial \omega} \over {\partial x_i}}, 10 ^{-10} \right )$

So let's start with turbulence kinetic energy:

Code:

    // Turbulent kinetic energy equation

    tmp<fvScalarMatrix> kEqn

    (

        fvm::ddt(alpha, rho, k_)

      + fvm::div(alphaRhoPhi, k_)

      - fvm::laplacian(alpha*rho*DkEff(F1), k_)

     ==

        alpha()*rho()*Pk(G)

      - fvm::SuSp((2.0/3.0)*alpha()*rho()*divU, k_)

      - fvm::Sp(alpha()*rho()*epsilonByk(F1, F23), k_)

      + kSource()

      + fvOptions(alpha, rho, k_)

    );

Code:

    tmp<volScalarField> CDkOmegaPlus = max

    (

        CDkOmega,

        dimensionedScalar("1.0e-10", dimless/sqr(dimTime), 1.0e-10)

    );



    tmp<volScalarField> 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));

where

Code:

    volScalarField CDkOmega

    (

        (2*alphaOmega2_)*(fvc::grad(k_) & fvc::grad(omega_))/omega_

    );

So what does DkEff(F1), also called effective diffusivity do? See kOmegaSSTBase.H:

Code:

        //- Return the effective diffusivity for k

        tmp<volScalarField> DkEff(const volScalarField& F1) const

        {

            return tmp<volScalarField>

            (

                new volScalarField("DkEff", alphaK(F1)*this->nut_ + this->nu())

            );

        }

I think this one looks very much like the term $\left( {\nu + \sigma_k \nu _T } \right)$ . Except we need to figure out alphaK function:

Code:

        tmp<volScalarField> alphaK(const volScalarField& F1) const

        {

            return blend(F1, alphaK1_, alphaK2_);

        }

So it calls blending function. Nothing fancy. Keep in mind that alpha = 1/sigma by the way. This is what blend function does:

Code:

        tmp<volScalarField::Internal> blend

        (

            const volScalarField::Internal& F1,

            const dimensionedScalar& psi1,

            const dimensionedScalar& psi2

        ) const

        {

            return F1*(psi1 - psi2) + psi2;

        }

$\nu _T = {a_1 k \over \mbox{max}(a_1 \omega, S F_2) }$

and the corresponding code is:

Code:

this->nut_ = a1_*k_/max(a1_*omega_, b1_*F2*sqrt(S2));

F2 is given above. See its code:

Code:

    tmp<volScalarField> 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));

and S2 is:

Code:

    

    volScalarField S2(2*magSqr(symm(tgradU())));

Code:

$FOAM_SRC/turbulenceModels/Base/kOmegaSST/kOmegaSSTBase.C

where return value for P_k is

Code:

return min(G, (c1_*betaStar_)*this->k_()*this->omega_());

where G is

Code:

    volScalarField::Internal GbyNu(dev(twoSymm(tgradU()())) && tgradU()());

    volScalarField::Internal G(this->GName(), nut()*GbyNu);

Code:

      - fvm::SuSp((2.0/3.0)*alpha()*rho()*divU, k_)

      - fvm::Sp(alpha()*rho()*epsilonByk(F1, F23), k_)

These two functions, SuSp and Sp are two methods of implicit source term implementation in OpenFOAM. You can refer to the Users Guide to understand their mechanism.

See epsilonByk function:

Code:

return betaStar_*omega_();