1

Previous: Description of turbulence model Up: Turbulence Modeling in OpenFOAM Next: Example of creating new

This is an automatically generated documentation by LaTeX2HTML utility. In case of any issue, please, contact us at info@cfdsupport.com.

View of the source code of turbulence model

  • Take a look at the two equation turbulence model $ k-\epsilon$ for incompressible fluid flow

     

  • Source code can be found in:

    $FOAM_SRC/TurbulenceModels/turbulenceModels/RAS/kEpsilon

     

  • In file kEpsilon.H there is defined class of $ k-\epsilon$ turbulence model and headers of class functions

     

  • In file kEpsilon.C the model itself is defined

     

  • Let us have a look at the file:

    # less $FOAM_SRC/TurbulenceModels/turbulenceModels/RAS/kEpsilon/kEpsilon.C

     

    /*---------------------------------------------------------------------------*\
      =========                 |
      \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
       \\    /   O peration     |
        \\  /    A nd           | Copyright (C) 2011-2017 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/>.
    
    \*---------------------------------------------------------------------------*/
    
  • OpenFOAM header

     

    #include "kEpsilon.H"
    #include "fvOptions.H"
    #include "bound.H"
    
    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
    
  • Include of header files

     

    namespace Foam
    {
    namespace RASModels
    {
    

     

  • Entering namespaces Foam, RASModels

     

    // * * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * //
    
    template<class BasicTurbulenceModel>
    void kEpsilon<BasicTurbulenceModel>::correctNut()
    {
        this->nut_ = Cmu_*sqr(k_)/epsilon_;
        this->nut_.correctBoundaryConditions();
        fv::options::New(this->mesh_).correct(this->nut_);
    
        BasicTurbulenceModel::correctNut();
    }
    
    
    template<class BasicTurbulenceModel>
    tmp<fvScalarMatrix> kEpsilon<BasicTurbulenceModel>::kSource() const
    {
        return tmp<fvScalarMatrix>
        (
            new fvScalarMatrix
            (
                k_,
                dimVolume*this->rho_.dimensions()*k_.dimensions()
                /dimTime
            )
        );
    }
    
    
    template<class BasicTurbulenceModel>
    tmp<fvScalarMatrix> kEpsilon<BasicTurbulenceModel>::epsilonSource() const
    {
        return tmp<fvScalarMatrix>
        (
            new fvScalarMatrix
            (
                epsilon_,
                dimVolume*this->rho_.dimensions()*epsilon_.dimensions()
                /dimTime
            )
        );
    }
    
    
    // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
    
    template<class BasicTurbulenceModel>
    kEpsilon<BasicTurbulenceModel>::kEpsilon
    (
        const alphaField& alpha,
        const rhoField& rho,
        const volVectorField& U,
        const surfaceScalarField& alphaRhoPhi,
        const surfaceScalarField& phi,
        const transportModel& transport,
        const word& propertiesName,
        const word& type
    )
    :
        eddyViscosity<RASModel<BasicTurbulenceModel>>
        (
            type,
            alpha,
            rho,
            U,
            alphaRhoPhi,
            phi,
            transport,
            propertiesName
        ),
    
        Cmu_
        (
            dimensioned<scalar>::lookupOrAddToDict
            (
                "Cmu",
                this->coeffDict_,
                0.09
            )
        ),
        C1_
        (
            dimensioned<scalar>::lookupOrAddToDict
            (
                "C1",
                this->coeffDict_,
                1.44
            )
        ),
        C2_
        (
            dimensioned<scalar>::lookupOrAddToDict
            (
                "C2",
                this->coeffDict_,
                1.92
            )
        ),
        C3_
        (
            dimensioned<scalar>::lookupOrAddToDict
            (
                "C3",
                this->coeffDict_,
                0
            )
        ),
        sigmak_
        (
            dimensioned<scalar>::lookupOrAddToDict
            (
                "sigmak",
                this->coeffDict_,
                1.0
            )
        ),
        sigmaEps_
        (
            dimensioned<scalar>::lookupOrAddToDict
            (
                "sigmaEps",
                this->coeffDict_,
                1.3
            )
        ),
    
        k_
        (
            IOobject
            (
                IOobject::groupName("k", U.group()),
                this->runTime_.timeName(),
                this->mesh_,
                IOobject::MUST_READ,
                IOobject::AUTO_WRITE
            ),
            this->mesh_
        ),
        epsilon_
        (
            IOobject
            (
                IOobject::groupName("epsilon", U.group()),
                this->runTime_.timeName(),
                this->mesh_,
                IOobject::MUST_READ,
                IOobject::AUTO_WRITE
            ),
            this->mesh_
        )
    {
        bound(k_, this->kMin_);
        bound(epsilon_, this->epsilonMin_);
    
        if (type == typeName)
        {
            this->printCoeffs(type);
        }
    }
    

     

  • Class $ k-\epsilon$ constructor
  • Constant definition
  • Definition of eddy viscosity
  • Reading of turbulent quantities

     

    // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
    
    template<class BasicTurbulenceModel>
    bool kEpsilon<BasicTurbulenceModel>::read()
    {
        if (eddyViscosity<RASModel<BasicTurbulenceModel>>::read())
        {
            Cmu_.readIfPresent(this->coeffDict());
            C1_.readIfPresent(this->coeffDict());
            C2_.readIfPresent(this->coeffDict());
            C3_.readIfPresent(this->coeffDict());
            sigmak_.readIfPresent(this->coeffDict());
            sigmaEps_.readIfPresent(this->coeffDict());
    
            return true;
        }
        else
        {
            return false;
        }
    }
    

     

  • Function read() reads model constants from turbulenceProperties

     

  • If coefficients not present in turbulenceProperties default coefficients are used

     

    template<class BasicTurbulenceModel>
    void kEpsilon<BasicTurbulenceModel>::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<BasicTurbulenceModel>>::correct();
    
        volScalarField::Internal divU
        (
            fvc::div(fvc::absolute(this->phi(), U))().v()
        );
    
        tmp<volTensorField> tgradU = fvc::grad(U);
        volScalarField::Internal G
        (
            this->GName(),
            nut.v()*(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_)
         ==
            C1_*alpha()*rho()*G*epsilon_()/k_()
          - fvm::SuSp(((2.0/3.0)*C1_ - C3_)*alpha()*rho()*divU, epsilon_)
          - fvm::Sp(C2_*alpha()*rho()*epsilon_()/k_(), 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_);
    
        correctNut();
    }
    
  • Function correct() solves differential equations of model:
  • $ G$ is production of turbulence
  • Equation for rate of dissipation of turbulent kinetic energy $ \epsilon$

     

    $\displaystyle \frac{\partial \epsilon}{\partial t} + u_{j} \sum_{j=1}^{3} \frac{\partial \epsilon}{\partial x_{j}} - \sum_{j=1}^{3} \frac{\partial}{\partial x_{j}} \left[ \left(\alpha_{\epsilon} \nu_{t} + \nu \right) \frac{\partial \epsilon}{\partial x_{j}} \right] = + C_{1\epsilon} \frac{\epsilon}{k} G - C_{2} \frac{\epsilon^{2}}{k}$(24.1)

     

     

  • Equation for turbulent kinetic energy $ k$:

     

    $\displaystyle \frac{\partial k}{\partial t} + u_{j} \sum_{j=1}^{3} \frac{\partial k}{\partial x_{j}} - \sum_{j=1}^{3} \frac{\partial}{\partial x_{j}} \left[ \left(\alpha_{k} \nu_{t} + \nu \right) \frac{\partial k}{\partial x_{j}} \right] = G - \epsilon$(24.2)

     

     

    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
    
    } // End namespace RASModels
    } // End namespace Foam
    
    // ************************************************************************* //
    
  • Closing namespaces
  •