CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM

Some confusions about SST-DES (code and formulas)

Register Blogs Community New Posts Updated Threads Search

Like Tree7Likes
  • 7 Post By fertinaz

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   April 21, 2018, 08:45
Question Some confusions about SST-DES (code and formulas)
  #1
Member
 
王莹
Join Date: May 2017
Posts: 51
Rep Power: 9
Alisa_W is on a distinguished road
Hello,
I'm trying to modify the SST-DES in OpenFOAM. I have found the main formulas of SST-DES and its code in OpenFOAM. However, when I compared the two parts so as to match the code with formulas one to one, it really confused me.....I can't even locate the key variables in the code.

I wonder if I found the wrong formulas of SST-DES, so I enclose the formulas and code. Please, any suggestion is welcome.

p.s. So sorry that I can't post the .doc file......Size of the file is limited.
Attached Images
File Type: png formulas.png (117.7 KB, 117 views)
Attached Files
File Type: txt code SST-DES main body.txt (5.6 KB, 22 views)
Alisa_W is offline   Reply With Quote

Old   April 22, 2018, 02:27
Default
  #2
Member
 
Fatih Ertinaz
Join Date: Feb 2011
Location: Istanbul
Posts: 64
Rep Power: 15
fertinaz is on a distinguished road
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

Last edited by fertinaz; April 22, 2018 at 02:52. Reason: Added references
fertinaz is offline   Reply With Quote

Old   April 24, 2018, 10:30
Default
  #3
Member
 
王莹
Join Date: May 2017
Posts: 51
Rep Power: 9
Alisa_W is on a distinguished road
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 View Post
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
Alisa_W is offline   Reply With Quote

Reply

Tags
code, des, modify, sst, turbulence model


Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Some confusions about SST-DES (code and formulas) Alisa_W OpenFOAM Programming & Development 1 April 12, 2018 23:29


All times are GMT -4. The time now is 12:04.