[muycfd]

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:

                        
 /*---------------------------------------------------------------------------*\
   =========                 |
   \\      /  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.

[klausb]

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?...

[muycfd]

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.

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:

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 }