
[Sponsors] 
April 5, 2011, 07:12 
TwoLayer kEpsilon Turbulence Model

#1 
New Member
Join Date: Apr 2011
Posts: 4
Rep Power: 6 
Hello everybody,
I'm modifying the standard kepsilon model included in OpenFOAM to include a twolayer 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 twolayer nearwall treatment in openfoam? Thanks a lot! James Code:
void Ketl::correct() { if (!turbulence_) { // Recalculate viscosity mut_ = rho_*Cmu_*sqr(k_)/(epsilon_ + epsilonSmall_); mut_.correctBoundaryConditions(); // Recalculate 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_); // Recalculate 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 kepsilon 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 * * * // Recalculate thermal diffusivity alphat_ = mut_/Prt_; alphat_.correctBoundaryConditions(); } [2] http://www.kxcad.net/STARCCM/online...ulence32.html 

April 6, 2011, 05:06 

#2 
Member

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 LowReynolds turbulence models. For the test case indicated you find a review of the near wall behaviour of several LowReynolds 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 you can find this article: 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
__________________
Cosimo Bianchini Ergon Research s.r.l. Via Panciatichi, 92 50127 Florence  ITALY Tel: +39 055 0763716 Mob: +39 320 9460153 email: cosimo.bianchini@ergonresearch.it URL: www.ergonresearch.it 

April 7, 2011, 05:11 

#3 
New Member
Join Date: Apr 2011
Posts: 4
Rep Power: 6 
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 

April 7, 2011, 10:56 

#4 
New Member
Join Date: Apr 2011
Posts: 4
Rep Power: 6 
I modified my twolayer 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 freestream region and the nearwall 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_ ) The viscosity is calculated as follows: Code:
// * * * N E W * * * // Recalculate 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 twolayer model of turbulence in calculation of a boundary layer with a pressure gradient", Journal of Engineering Physics and Thermodynamics, Volume 80 

February 11, 2013, 02:00 

#6 
New Member
OpenFoam
Join Date: Jul 2012
Posts: 24
Rep Power: 5 
Hi all,
Did anyone try and succeed in implementing 2 layer model.... i am trying to do the same into realizableKE model. 

May 4, 2013, 12:26 
I Did It..!

#7 
New Member
OpenFoam
Join Date: Jul 2012
Posts: 24
Rep Power: 5 
Finally i did it...!! thanks for help of all forum members..!!


February 21, 2014, 13:43 

#8 
Member
pooyan
Join Date: Nov 2011
Posts: 62
Rep Power: 6 
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, 

Thread Tools  
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
gammaReTheta turbulence model for predicting transitional flows  FelixL  OpenFOAM Programming & Development  79  October 11, 2015 15:30 
Turbulence Model and limitation to Reynolds number  qascapri  FLUENT  0  January 24, 2011 11:48 
Sato's model of bubble induced turbulence  ukbid  CFX  0  January 3, 2011 10:04 
how could i define a custom turbulence model  FredPacheo  FLUENT  0  July 24, 2008 11:06 
Wall turbulence, viscosity, boundary layer  Patrick Godon  Main CFD Forum  1  November 5, 2003 16:39 