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

Two-Layer k-Epsilon Turbulence Model

Register Blogs Members List Search Today's Posts Mark Forums Read

Like Tree9Likes
  • 1 Post By james.t
  • 2 Post By cosimobianchini
  • 2 Post By james.t
  • 4 Post By neeraj

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   April 5, 2011, 08:12
Default Two-Layer k-Epsilon Turbulence Model
  #1
New Member
 
Join Date: Apr 2011
Posts: 4
Rep Power: 14
james.t is on a distinguished road
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
[2] http://www.kxcad.net/STAR-CCM/online...ulence-32.html
lpz456 likes this.
james.t is offline   Reply With Quote

Old   April 6, 2011, 06:06
Default
  #2
Member
 
cosimo bianchini
Join Date: Mar 2009
Location: Florence, Tuscany, Italy
Posts: 88
Rep Power: 17
cosimobianchini is on a distinguished road
Send a message via Skype™ to cosimobianchini
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
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
calim_cfd and lpz456 like this.
__________________
Cosimo Bianchini

Ergon Research s.r.l.
Via Panciatichi, 92
50127 Florence - ITALY
Tel: +39 055 0763716
Mob: +39 320 9460153
e-mail: cosimo.bianchini@ergonresearch.it
URL: www.ergonresearch.it
cosimobianchini is offline   Reply With Quote

Old   April 7, 2011, 06:11
Default
  #3
New Member
 
Join Date: Apr 2011
Posts: 4
Rep Power: 14
james.t is on a distinguished road
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 is offline   Reply With Quote

Old   April 7, 2011, 11:56
Default
  #4
New Member
 
Join Date: Apr 2011
Posts: 4
Rep Power: 14
james.t is on a distinguished road
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
sam1364 and randolph like this.
james.t is offline   Reply With Quote

Old   January 25, 2012, 13:12
Default
  #5
Senior Member
 
lakeat's Avatar
 
Daniel WEI (老魏)
Join Date: Mar 2009
Location: Beijing, China
Posts: 689
Blog Entries: 9
Rep Power: 21
lakeat is on a distinguished road
Send a message via Skype™ to lakeat
Hi James,

Did you make it?
__________________
~
Daniel WEI
-------------
Boeing Research & Technology - China
Beijing, China
Email
lakeat is offline   Reply With Quote

Old   February 11, 2013, 02:00
Default
  #6
New Member
 
neeraj's Avatar
 
OpenFoam
Join Date: Jul 2012
Posts: 24
Rep Power: 13
neeraj is on a distinguished road
Hi all,
Did anyone try and succeed in implementing 2 layer model.... i am trying to do the same into realizableKE model.
neeraj is offline   Reply With Quote

Old   May 4, 2013, 13:26
Default I Did It..!
  #7
New Member
 
neeraj's Avatar
 
OpenFoam
Join Date: Jul 2012
Posts: 24
Rep Power: 13
neeraj is on a distinguished road
Finally i did it...!! thanks for help of all forum members..!!
neeraj is offline   Reply With Quote

Old   February 21, 2014, 13:43
Default
  #8
Member
 
pooyan
Join Date: Nov 2011
Posts: 62
Rep Power: 14
sam1364 is on a distinguished road
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,
sam1364 is offline   Reply With Quote

Old   March 25, 2016, 02:34
Default
  #9
New Member
 
M Kashif Tehseen
Join Date: Feb 2016
Location: Islamabad
Posts: 4
Rep Power: 10
Tehseen is on a distinguished road
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.

Thanks in advance
Tehseen is offline   Reply With Quote

Reply

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

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
gamma-ReTheta turbulence model for predicting transitional flows FelixL OpenFOAM Programming & Development 123 August 30, 2022 12:50
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 12:06
Wall turbulence, viscosity, boundary layer Patrick Godon Main CFD Forum 1 November 5, 2003 16:39


All times are GMT -4. The time now is 06:57.