CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   OpenFOAM Programming & Development (http://www.cfd-online.com/Forums/openfoam-programming-development/)
-   -   Two-Layer k-Epsilon Turbulence Model (http://www.cfd-online.com/Forums/openfoam-programming-development/86888-two-layer-k-epsilon-turbulence-model.html)

 james.t April 5, 2011 07:12

Two-Layer k-Epsilon Turbulence Model

Hello everybody,

I'm modifying the standard k-epsilon model included in OpenFOAM to include a two-layer boundary treatment according to [1, 2]. My goal is to achieve a grid independent solution. My reference case is the simulation of a turbulent boundary layer and comparison to the results of Whtie and Wieghardt.

I directly changed the code in the kEpsilon.C to solve the transport equation for k and epsilon for the whole field (no modification so far). For lower Reynolds numbers Re_y < 200, one equation for epsilon and one for mut are calculated to overwrite the results of the previous epsilon transport equation (see source code).

However, directly at the wall Re_y goes to zero, so mut would be zero and epsilon infinity. How should I treat these values directly at the wall? Has someone experience with the two-layer near-wall treatment in openfoam?

Thanks a lot!

James

Code:

``` void Ketl::correct() {         if (!turbulence_)         {                 // Re-calculate viscosity                 mut_ = rho_*Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_);                 mut_.correctBoundaryConditions();                 // Re-calculate thermal diffusivity                 alphat_ = mut_/Prt_;                 alphat_.correctBoundaryConditions();                 return;         }         RASModel::correct();         volScalarField divU = fvc::div(phi_/fvc::interpolate(rho_));         if (mesh_.moving())         {                 divU += fvc::div(mesh_.phi());         }         tmp<volTensorField> tgradU = fvc::grad(U_);         volScalarField G("RASModel::G", mut_*(tgradU() && dev(twoSymm(tgradU()))));         tgradU.clear();         // Update espsilon and G at the wall         epsilon_.boundaryField().updateCoeffs();         // Dissipation equation         tmp<fvScalarMatrix> epsEqn         (                 fvm::ddt(rho_, epsilon_)             + fvm::div(phi_, epsilon_)             - fvm::laplacian(DepsilonEff(), epsilon_)         ==                 C1_*G*epsilon_/k_             - fvm::SuSp(((2.0/3.0)*C1_ + C3_)*rho_*divU, epsilon_)             - fvm::Sp(C2_*rho_*epsilon_/k_, epsilon_)         );         epsEqn().relax();         epsEqn().boundaryManipulate(epsilon_.boundaryField());         solve(epsEqn);         bound(epsilon_, epsilon0_);         // * * * N E W * * *         // wall distance         volScalarField y_ = wallDist(mesh_).y();         // reynolds number based on wall distance         Rey_ = rho_ * y_ * sqrt(k_) / mu();         // constants         scalar Cmu75 = pow(Cmu_.value(), 0.75);         scalar kappa_ = 0.42;         scalar Aeps_ = 2 * kappa_ / Cmu75;         // loop over all cells         forAll(Rey_, cellI)         {                 if(Rey_[cellI] < 200)                 {                         // length scale                         scalar Leps = y_[cellI] * kappa_ / Cmu75 * (1 - exp( -Rey_[cellI] / Aeps_ ));                         // dissipation                         epsilon_[cellI] = pow(k_[cellI], 1.5) / Leps;                 }         }         // * * * N E W * * *         // Turbulent kinetic energy equation         tmp<fvScalarMatrix> kEqn         (                 fvm::ddt(rho_, k_)             + fvm::div(phi_, k_)             - fvm::laplacian(DkEff(), k_)         ==                 G             - fvm::SuSp((2.0/3.0)*rho_*divU, k_)             - fvm::Sp(rho_*epsilon_/k_, k_)         );         kEqn().relax();         solve(kEqn);         bound(k_, k0_);         // Re-calculate viscosity         mut_ = rho_*Cmu_*sqr(k_)/epsilon_;         mut_.correctBoundaryConditions();                 // * * * N E W * * *         // constants         scalar Amu_ = 70;         scalar A_ = 10 / tanh(0.98);         // loop over all cells         forAll(Rey_, cellI)         {                 if(Rey_[cellI] < 200)                 {                         // length scale                         scalar Lmu = y_[cellI] * kappa_ / Cmu75 * (1 - exp( -Rey_[cellI] / Amu_ ) );                         // viscosity according to the standard k-epsilon model                         scalar mutKE = rho_[cellI] * Cmu_.value() * sqr(k_[cellI]) / (epsilon_[cellI] + epsilonSmall_.value());                         // viscosity according to wolfstein                         scalar mutTL = rho_[cellI] * Cmu_.value() * Lmu * sqrt(k_[cellI]);                         // blending function                         scalar lambda_ = 0.5* ( 1 + tanh( (Rey_[cellI] - 200) / A_ ) );                         // blended viscosity                         mut_[cellI] = lambda_ * mutKE + (1 - lambda_) * mutTL;                 }         }         // * * * N E W * * *         // Re-calculate thermal diffusivity         alphat_ = mut_/Prt_;         alphat_.correctBoundaryConditions(); }```
[1] http://my.fit.edu/itresources/manual...ug/node514.htm

 cosimobianchini April 6, 2011 05:06

Hi James,
simply updating epsilon field after solving the transport equation would not be enough to obtain solution with smooth epsilon and mut field. You should consider exploit the setValues member function of fvMatrix to ensure that you satisfy at the same time the extended wall function in the near wall zone as well as the transport equation in the free jet.
Concerning the limiting behavior at the wall for the turbulent quantities your are right that y and so mut tend towards 0, however k is going to zero as well and the value of epsilon should be finite. This behavior is common to all Low-Reynolds turbulence models. For the test case indicated you find a review of the near wall behaviour of several Low-Reynolds turbulence model in:

V. C. Patel, W. Rodi, and G. Sheuerer. Turbulence models
for near wall and low reynolds number flows: a review.
AIAA Journal, 26:1308–1319, 1993

Here:http://www.opensourcecfd.com/conference2008/2007/index.php
Heat Transfer Applications in Turbomachinery - L. Mangani, C. Bianchini
dealing partially with the test case you are referring to.
More results could be find in: http://powerlab.fsb.hr/ped/kturbo/Op...aniPhD2008.pdfhttp://powerlab.fsb.hr/ped/kturbo/Op...aniPhD2008.pdf
Hope you find this interesting,
Cosimo

 james.t April 7, 2011 05:11

Hi cosimo,

thank you very much for your help!

I just found out, that I forgot to implement the blending function for epsilon! One of the reaons, my results weren't that good. Also thanks for the setValues command, I already included it and it seems to work perfectly :-)

Thanks again for the papers! I hadn't good reference so far, so that's invaluable for me!

Regards

James

 james.t April 7, 2011 10:56

I modified my two-layer model according to the paper of Volkov [1]. Now a blending function lambda is directly included in the dissipation transport equation to distinguish between the free-stream region and the near-wall region. (Limit: Re_y = 200)

The transport equation for epsilon is as follows:
Code:

```        // Update espsilon and G at the wall         epsilon_.boundaryField().updateCoeffs();```
Code:

```         // Dissipation equation         tmp<fvScalarMatrix> epsEqn         (                 fvm::ddt(rho_, epsilon_)             + lambda_ * fvm::div(phi_, epsilon_)             - lambda_ * fvm::laplacian(DepsilonEff(), epsilon_)         ==               lambda_            * C1_*G*epsilon_/k_             - lambda_            * fvm::SuSp(((2.0/3.0)*C1_ + C3_)*rho_*divU, epsilon_)             - lambda_            * fvm::Sp(C2_*rho_*epsilon_/k_, epsilon_)             + (1 - lambda_)  * alpha * rho_ * ( pow(k_, 1.5)/Leps_ - epsilon_ )         );         epsEqn().relax();         epsEqn().boundaryManipulate(epsilon_.boundaryField());         solve(epsEqn);         bound(epsilon_, epsilon0_); ```

where lambda_ is the blending function defined as
Code:

```        // wall distance         volScalarField y_ = wallDist(mesh_).y();         // reynolds number based on wall distance         Rey_ = rho_ * y_ * sqrt(k_) / mu();         // blending function         scalar A_ = 10 / tanh(0.98);         volScalarField lambda_ = 0.5* ( 1 + tanh( (Rey_ - 200) / A_ ) );```

lambda is zero at the wall, so that just
Code:

``` fvm::ddt(rho_, epsilon_) == + alpha * rho_ * ( pow(k_, 1.5)/Leps_ - epsilon_ ) ```
is solved.

The viscosity is calculated as follows:
Code:

```        // * * * N E W * * *         // Re-calculate viscosity         scalar Amu_ = 70;         volScalarField Lmu = y_ * kappa_ / Cmu75 * (1 - exp( -Rey_ / Amu_ ) );         volScalarField mut_ke = rho_*Cmu_*sqr(k_)/epsilon_;         volScalarField mut_tl = rho_*Cmu_*Lmu*sqrt(k_);         mut_ = lambda_ * mut_ke + (1 - lambda_) * mut_tl;         // * * * N E W * * *         mut_.correctBoundaryConditions();```

However, the simulations always show some regions with very high turbulent viscosity outside the boundary layer (Re_y > 200) with values of mut = 10 mio and more.

Can someone give me an advise, where the mistake might be?

Thanks a lot!

Regards

James

[1]
K.N. Volkov,
"Application of a two-layer model of turbulence in calculation of a boundary layer with a pressure gradient", Journal of Engineering Physics and Thermodynamics, Volume 80

 lakeat January 25, 2012 13:12

Hi James,

Did you make it?

 neeraj February 11, 2013 02:00

Hi all,
Did anyone try and succeed in implementing 2 layer model.... i am trying to do the same into realizableKE model.

 neeraj May 4, 2013 12:26

I Did It..!

Finally i did it...!! thanks for help of all forum members..!!

 sam1364 February 21, 2014 13:43

Hey Neeraj,

can you please post your code here. I wanna do the same thing and I would like to check on that with you.

Thanks,

 Tehseen March 25, 2016 02:34

Dear Neeraj, will you please like to share the code as it will be much valuable for me in my thesis writing. You can contact me at mktehseen@gmail.com too.