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

Underprediction Issues (Turbulence Model Implementation)

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

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   June 28, 2021, 16:17
Default Underprediction Issues (Turbulence Model Implementation)
  #1
New Member
 
Join Date: Jun 2020
Posts: 2
Rep Power: 0
muycfd is on a distinguished road
Dear Foamers,


I am a fairly beginner user, trying to implement a RANS turbulence model (k-e-t) that I found in the literature, which has multiple time scales. As far as I know, this model and its coefficients are tuned for plane jet case (for shear-free flow). It is an incompressible flow case. Up to this point, I have tried many things but couldn’t match it neither with the experimental results nor with published simulation results of it. My simulations show slight underprediction in turbulence quantities (i.e. k/Uc^2, uv/Uc^2). However, velocity distribution is almost the same. Any suggestion or corrections would be invaluable. Model equations is attached as an image. And the code can be found below:



Code:
                        
 /*---------------------------------------------------------------------------*\
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
     \\  /    A nd           | Copyright (C) 2011-2020 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
     This file is part of OpenFOAM.
 
 
     OpenFOAM is free software: you can redistribute it and/or modify it
     under the terms of the GNU General Public License as published by
     the Free Software Foundation, either version 3 of the License, or
     (at your option) any later version.
 
 
     OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
     ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
     FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
     for more details.
 
 
     You should have received a copy of the GNU General Public License
     along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
 
 
 \*---------------------------------------------------------------------------*/
 
 
 #include "kET.H"
 #include "fvOptions.H"
 #include "bound.H"
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 
 namespace Foam
 {
 namespace RASModels
 {
 
 
 // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
 
 
 template<class BasicMomentumTransportModel>
 void kET<BasicMomentumTransportModel>::correctNut()
 {
     this->nut_ = Cmu_*sqr(k_)/epsilon_;
     this->nut_.correctBoundaryConditions();
     fv::options::New(this->mesh_).correct(this->nut_);
 }
 
 
 
 
 template<class BasicMomentumTransportModel>
 tmp<fvScalarMatrix> kET<BasicMomentumTransportModel>::kSource() const
 {
     return tmp<fvScalarMatrix>
     (
         new fvScalarMatrix
         (
             k_,
             dimVolume*this->rho_.dimensions()*k_.dimensions()
             /dimTime
         )
     );
 }
 
 
 
 
 template<class BasicMomentumTransportModel>
 tmp<fvScalarMatrix> kET<BasicMomentumTransportModel>::epsilonSource() const
 {
     return tmp<fvScalarMatrix>
     (
         new fvScalarMatrix
         (
             epsilon_,
             dimVolume*this->rho_.dimensions()*epsilon_.dimensions()
             /dimTime
         )
     );
 }
 
 
 template<class BasicMomentumTransportModel>
 tmp<fvScalarMatrix> kET<BasicMomentumTransportModel>::taoSource() const
 {
     return tmp<fvScalarMatrix>
     (
         new fvScalarMatrix
         (
             tao_,
             dimVolume*this->rho_.dimensions()*tao_.dimensions()
             /dimTime
         )
     );
 }
 
 
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 
 template<class BasicMomentumTransportModel>
 kET<BasicMomentumTransportModel>::kET
 (
     const alphaField& alpha,
     const rhoField& rho,
     const volVectorField& U,
     const surfaceScalarField& alphaRhoPhi,
     const surfaceScalarField& phi,
     const transportModel& transport,
     const word& type
 )
 :
     eddyViscosity<RASModel<BasicMomentumTransportModel>>
     (
         type,
         alpha,
         rho,
         U,
         alphaRhoPhi,
         phi,
         transport
     ),
 
 
     Cmu_
     (
         dimensioned<scalar>::lookupOrAddToDict
         (
             "Cmu",
             this->coeffDict_,
             0.09
         )
     ),
     Ceps1_
     (
         dimensioned<scalar>::lookupOrAddToDict
         (
             "Ceps1",
             this->coeffDict_,
             0.75
         )
     ),
     Ceps2_
     (
         dimensioned<scalar>::lookupOrAddToDict
         (
             "Ceps2",
             this->coeffDict_,
             1.0
         )
     ),
     Ceps3_
     (
         dimensioned<scalar>::lookupOrAddToDict
         (
             "Ceps3",
             this->coeffDict_,
             0.67
         )
     ),
     Ctao0_
     (
         dimensioned<scalar>::lookupOrAddToDict
         (
             "Ctao0",
             this->coeffDict_,
             1.054
         )
     ),
     Ctao1_
     (
         dimensioned<scalar>::lookupOrAddToDict
         (
             "Ctao1",
             this->coeffDict_,
             1.1
         )
     ),
     Ctao2_
     (
         dimensioned<scalar>::lookupOrAddToDict
         (
             "Ctao2",
             this->coeffDict_,
             0.59
         )
     ),
     Ctao3_
     (
         dimensioned<scalar>::lookupOrAddToDict
         (
             "Ctao3",
             this->coeffDict_,
             0.83
         )
     ),
     sigmak_
     (
         dimensioned<scalar>::lookupOrAddToDict
         (
             "sigmak",
             this->coeffDict_,
             1.0
         )
     ),
     sigmaEps_
     (
         dimensioned<scalar>::lookupOrAddToDict
         (
             "sigmaEps",
             this->coeffDict_,
             1.2
         )
     ),
     sigmaTao_
     (
         dimensioned<scalar>::lookupOrAddToDict
         (
             "sigmaTao",
             this->coeffDict_,
             1.1
         )
     ),
 
 
     k_
     (
         IOobject
         (
             IOobject::groupName("k", alphaRhoPhi.group()),
             this->runTime_.timeName(),
             this->mesh_,
             IOobject::MUST_READ,
             IOobject::AUTO_WRITE
         ),
         this->mesh_
     ),
     epsilon_
     (
         IOobject
         (
             IOobject::groupName("epsilon", alphaRhoPhi.group()),
             this->runTime_.timeName(),
             this->mesh_,
             IOobject::MUST_READ,
             IOobject::AUTO_WRITE
         ),
         this->mesh_
     ),
     tao_
     (
         IOobject
         (
             IOobject::groupName("tao", alphaRhoPhi.group()),
             this->runTime_.timeName(),
             this->mesh_,
             IOobject::MUST_READ,
             IOobject::AUTO_WRITE
         ),
         this->mesh_
     )
 {
     bound(k_, this->kMin_);
     bound(epsilon_, this->epsilonMin_);
     bound(tao_, this->taoMin_);
 
 
     if (type == typeName)
     {
         this->printCoeffs(type);
     }
 }
 
 
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 
 template<class BasicMomentumTransportModel>
 bool kET<BasicMomentumTransportModel>::read()
 {
     if (eddyViscosity<RASModel<BasicMomentumTransportModel>>::read())
     {
         Cmu_.readIfPresent(this->coeffDict());
         Ceps1_.readIfPresent(this->coeffDict());
         Ceps2_.readIfPresent(this->coeffDict());
         Ceps3_.readIfPresent(this->coeffDict());
         Ctao0_.readIfPresent(this->coeffDict());
         Ctao1_.readIfPresent(this->coeffDict());
         Ctao2_.readIfPresent(this->coeffDict());
         Ctao3_.readIfPresent(this->coeffDict());
         sigmak_.readIfPresent(this->coeffDict());
         sigmaEps_.readIfPresent(this->coeffDict());
         sigmaTao_.readIfPresent(this->coeffDict());
 
 
         return true;
     }
     else
     {
         return false;
     }
 }
 
 
 
 
 template<class BasicMomentumTransportModel>
 void kET<BasicMomentumTransportModel>::correct()
 {
     if (!this->turbulence_)
     {
         return;
     }
 
 
     // Local references
     const alphaField& alpha = this->alpha_;
     const rhoField& rho = this->rho_;
     const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_;
     const volVectorField& U = this->U_;
     volScalarField& nut = this->nut_;
     fv::options& fvOptions(fv::options::New(this->mesh_));
 
 
     eddyViscosity<RASModel<BasicMomentumTransportModel>>::correct();
 
 
     volScalarField::Internal divU
     (
         fvc::div(fvc::absolute(this->phi(), U))()
     );
 
 
     tmp<volTensorField> tgradU = fvc::grad(U);
      
     volScalarField Str(2*magSqr(dev(symm(tgradU()))));
     volScalarField magStr(sqrt(Str));
      
     volScalarField::Internal G
     (
         this->GName(),
         nut()*(dev(twoSymm(tgradU().v())) && tgradU().v())
     );
     tgradU.clear();
 
 
     // Update epsilon and G at the wall
     epsilon_.boundaryFieldRef().updateCoeffs();
 
 
     // Dissipation equation
     tmp<fvScalarMatrix> epsEqn
     (
         fvm::ddt(alpha, rho, epsilon_)
       + fvm::div(alphaRhoPhi, epsilon_)
       - fvm::laplacian(alpha*rho*DepsilonEff(), epsilon_)
      ==
         Ceps1_*alpha()*rho()*G/tao_()
       - fvm::SuSp(((2.0/3.0)*Ceps1_)*alpha()*rho()*k_()*divU/tao_()/epsilon_(), epsilon_)
       - fvm::Sp(Ceps2_*alpha()*rho()/tao_(), epsilon_)
       - fvm::SuSp(Ceps3_*alpha()*rho()*magStr(), epsilon_)
       + epsilonSource()
       + fvOptions(alpha, rho, epsilon_)
     );
 
 
     epsEqn.ref().relax();
     fvOptions.constrain(epsEqn.ref());
     epsEqn.ref().boundaryManipulate(epsilon_.boundaryFieldRef());
     solve(epsEqn);
     fvOptions.correct(epsilon_);
     bound(epsilon_, this->epsilonMin_);
 
 
     // Turbulent kinetic energy equation
     tmp<fvScalarMatrix> kEqn
     (
         fvm::ddt(alpha, rho, k_)
       + fvm::div(alphaRhoPhi, k_)
       - fvm::laplacian(alpha*rho*DkEff(), k_)
      ==
         alpha()*rho()*G
       - fvm::SuSp((2.0/3.0)*alpha()*rho()*divU, k_)
       - fvm::Sp(alpha()*rho()*epsilon_()/k_(), k_)
       + kSource()
       + fvOptions(alpha, rho, k_)
     );
 
 
     kEqn.ref().relax();
     fvOptions.constrain(kEqn.ref());
     solve(kEqn);
     fvOptions.correct(k_);
     bound(k_, this->kMin_);
      
     // Tao equation
     tmp<fvScalarMatrix> taoEqn
     (
         fvm::ddt(alpha, rho, tao_)
       + fvm::div(alphaRhoPhi, tao_)
       - fvm::laplacian(alpha*rho*DtaoEff(), tao_)
      ==
         Ctao0_*alpha()*rho()
       - fvm::Sp(Ctao1_*alpha()*rho()*epsilon_()/k_(),tao_)
       - fvm::Sp(Ctao2_*alpha()*rho()*G/k_(),tao_)
       + fvm::SuSp(((2.0/3.0)*Ctao2_)*alpha()*rho()*divU, tao_)
       - fvm::SuSp(Ctao3_*alpha()*rho()*magStr(), tao_)
       + taoSource()
       + fvOptions(alpha, rho, tao_)
     );
 
 
     taoEqn.ref().relax();
     fvOptions.constrain(taoEqn.ref());
     solve(taoEqn);
     fvOptions.correct(tao_);
     bound(tao_, this->taoMin_);
 
 
     correctNut();
 }
 
 
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 
 } // End namespace RASModels
 } // End namespace Foam
 
 
 // ************************************************************************* //
Thanks,


P.S.: I just used tao instead of tau to prevent any unforeseen match of variables.
Attached Images
File Type: jpg IMG_2407.jpg (152.7 KB, 10 views)
muycfd is offline   Reply With Quote

Old   June 28, 2021, 17:47
Default
  #2
Senior Member
 
Klaus
Join Date: Mar 2009
Posts: 207
Rep Power: 19
klausb will become famous soon enough
please provide the full documentation of the k-e-tau turbulence model and reference cases

What does you mesh look like? Did you double check your turbulence properties?...
klausb is offline   Reply With Quote

Old   July 9, 2021, 07:33
Default
  #3
New Member
 
Join Date: Jun 2020
Posts: 2
Rep Power: 0
muycfd is on a distinguished road
Thank you for your suggestions Klaus. And sorry for the delay, multiple important stuff came up. I tried k-epsilon and other well-known RANS models before starting k-e-t and results of the k-epsilon fit perfectly for the same case when I compare it with experimental studies. I used the same domain for all the simulations. So, I think the grid may not be the reason here.


I have forgotten to add that z is equal to (epsilon*tau)/k.


In addition, I found that for the incompressible cases, the terms at the end of epsilon and tau equations were omitted (steps to get magStr too). The terms are given below to be clear.


Code:
fvm::SuSp(Ceps3_*alpha()*rho()*magStr(), epsilon_)
 
 
 fvm::SuSp(Ctao3_*alpha()*rho()*magStr(), tao_)
I tried both cases with and without these terms, closer results were for without conditions.


For the boundary conditions, I used the same k and epsilon values (as it is calculated for the k epsilon model) at first and calculated tau using the given relations but then I found that they used different equations such as below:


Code:
k=.005*Ui^2
 epsilon=Cmu*k^2/(7*nu)
 tau=k/(epsilon+SMALL)/1.92
None of them were successful.


A professor told me that it could be related to source term linearization. However, I don’t know how to change the coefficients of Sc or Sp. I have scanned the src folder with grep -R but couldn’t reach any related file. From what I discover on that subject is that it affects only the stability and convergence.


For the reference cases, I used the same experimental studies and one numerical study (s-thesis.pdf).



I also added related information about the k-e-t model.



I hope these help for you to make recommendations.


Thanks!




p { margin-bottom: 0.1in; line-height: 115%; background: transparent }
muycfd is offline   Reply With Quote

Reply

Tags
implementation, rans, turbulence model, underprediction

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
Antti Hellsten model: New Advanced k–ω Turbulence Model purnp2 OpenFOAM Programming & Development 0 May 9, 2019 20:16
Only two turbulence options available in CFX Pre Jack001 CFX 5 March 30, 2016 03:47
implementation of Spalart-Allmaras Turbulence Model zhengjg Main CFD Forum 0 July 24, 2013 04:43
SA Turbulence model implementation ganesh Main CFD Forum 0 March 6, 2006 13:23
k-w Turbulence model implementation suneesh Main CFD Forum 4 November 23, 2005 18:35


All times are GMT -4. The time now is 20:26.