Negative Spalart Allmaras compressible model for OpenFOAM

Negative Spalart Allmaras model (labeled SA-neg)[1] was implemented for compressible Spalart Allmaras model in OpenFOAM.

New source codes are shown below. Spalart Shur Rotation/Curvature Correction (SA-RC), and Compressibility Correction are also included. Differences from the original codes are indicated in red for Compressibility Correction, blue for Spalart Shur Correction, and magenta for Negative SA model.

New Stilda modification[1] that prevents negative Stilda is also added in this version.

References
1. S.R. Allmaras, F.T. Johnson, P.R. Spalart, Modifications and Clarifications for the Implementation of the Spalart-Allmaras Turbulence Model, ICCFD7-1902, 7th International Conference on Computational Fluid Dynamics, 2012.

SpalartAllmaras.H

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

Class
Foam::compressible::RASModels::SpalartAllmaras

Group
grpCmpRASTurbulence

Description
Spalart-Allmaras one-eqn mixing-length model for compressible
external flows.

Reference:
\verbatim
"A One-Equation Turbulence Model for Aerodynamic Flows"
P.R. Spalart,
S.R. Allmaras,
La Recherche Aerospatiale, No. 1, 1994, pp. 5-21.
\endverbatim

Extended according to:
\verbatim
"An Unstructured Grid Generation and Adaptive Solution Technique
for High Reynolds Number Compressible Flows"
G.A. Ashford,
Ph.D. thesis, University of Michigan, 1996.
\endverbatim
by using the optional flag \c ashfordCorrection,

and
\verbatim
"Turbulence Modeling in Rotating and Curved Channels: Assessing the
Spalart-Shur Correction"
M.L. Shur,
M.K. Strelets,
A.K. Travin,
P.R. Spalart,
AIAA Journal Vol. 38, No. 5, 2000, pp. 784-792.
\endverbatim
using the optional flag \c spalartShurCorrection,

and
\verbatim
"Trends in Turbulence Treatments"
P.R. Spalart,
AIAA Paper No. 2000-2603, 2000.
\endverbatim
using the optional flag \c compressibilityCorrection,

and
\verbatim
"Modifications and Clarifications for the Implementation of the
Spalart-Allmaras Turbulence Model"
S.R. Allmaras,
F.T. Johnson,
P.R. Spalart,
ICCFD7-1902, 7th International Conference on
Computational Fluid Dynamics, 2012.
\endverbatim
using the optional flag \c StildaModification,
and \c negativeNuTilda

The default model coefficients correspond to the following:
\verbatim
SpalartAllmarasCoeffs
{
Cb1 0.1355;
Cb2 0.622;
Cw2 0.3;
Cw3 2.0;
Cv1 7.1;
Cv2 5.0;
sigmaNut 0.66666;
Prt 1.0; // only for compressible
kappa 0.41;
Cr1 1.0;
Cr2 12.0;
Cr3 1.0;
C5 3.5;
Cn1 16.0;
Cn2 0.7;
Cn3 0.9;
Ct3 1.2;
Ct4 0.5;
}
\endverbatim

SourceFiles
SpalartAllmaras.C

\*---------------------------------------------------------------------------*/

#ifndef compressibleSpalartAllmaras_H
#define compressibleSpalartAllmaras_H

#include "RASModel.H"
#include "wallDist.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

namespace Foam
{
namespace compressible
{
namespace RASModels
{

/*---------------------------------------------------------------------------*\
Class SpalartAllmaras Declaration
\*---------------------------------------------------------------------------*/

class SpalartAllmaras
:
public RASModel
{

protected:

// Protected data

// Model coefficients

dimensionedScalar sigmaNut_;
dimensionedScalar kappa_;
dimensionedScalar Prt_;

dimensionedScalar Cb1_;
dimensionedScalar Cb2_;
dimensionedScalar Cw1_;
dimensionedScalar Cw2_;
dimensionedScalar Cw3_;
dimensionedScalar Cv1_;
dimensionedScalar Cv2_;
dimensionedScalar Cr1_;
dimensionedScalar Cr2_;
dimensionedScalar Cr3_;
dimensionedScalar C5_;
dimensionedScalar Cn1_;
dimensionedScalar Cn2_;
dimensionedScalar Cn3_;
dimensionedScalar Ct3_;
dimensionedScalar Ct4_;

//- Optional flag to activate the Ashford correction
Switch ashfordCorrection_;

//- Optional flag to activate the compressibility correction
Switch compressibilityCorrection_;

//- Optional flag to activate the Rotation/Curvature correction
Switch spalartShurCorrection_;

//- Optional flag to activate the new STilda correction
Switch StildaModification_;

//- Optional flag to activate the negative nuTilda model
Switch negativeNuTilda_;

// Fields

volScalarField nuTilda_;
volScalarField mut_;
volScalarField alphat_;

volScalarField Sxx_;
volScalarField Sxy_;
volScalarField Sxz_;
volScalarField Syx_;
volScalarField Syy_;
volScalarField Syz_;
volScalarField Szx_;
volScalarField Szy_;
volScalarField Szz_;

//- Wall distance
wallDist d_;

// Protected Member Functions

tmp<volScalarField> chi() const;
tmp<volScalarField> fv1(const volScalarField& chi) const;
tmp<volScalarField> fv2
(
const volScalarField& chi,
const volScalarField& fv1
) const;
tmp<volScalarField> fv3
(
const volScalarField& chi,
const volScalarField& fv1
) const;
tmp<volScalarField> fw(const volScalarField& Stilda) const;

tmp<volScalarField> fr1
(
const volSymmTensorField& S,
const volTensorField& W
) const;

tmp<volScalarField> ft2(const volScalarField& chi) const;

public:

//- Runtime type information
TypeName("SpalartAllmaras");

// Constructors

//- Construct from components
SpalartAllmaras
(
const volScalarField& rho,
const volVectorField& U,
const surfaceScalarField& phi,
const fluidThermo& thermophysicalModel,
const word& turbulenceModelName = turbulenceModel::typeName,
const word& modelName = typeName
);

//- Destructor
virtual ~SpalartAllmaras()
{}

// Member Functions

//- Return the effective diffusivity for nuTilda
virtual tmp<volScalarField> DnuTildaEff(const volScalarField& chi) const;
/*
{
return tmp<volScalarField>
(
new volScalarField
(
"DnuTildaEff",
(rho_*nuTilda_ + mu())/sigmaNut_
)
);
}
*/
//- Return the turbulence viscosity
virtual tmp<volScalarField> mut() const
{
return mut_;
}

//- Return the turbulence thermal diffusivity
virtual tmp<volScalarField> alphat() const
{
return alphat_;
}

//- Return the turbulence kinetic energy
virtual tmp<volScalarField> k() const;

//- Return the turbulence kinetic energy dissipation rate
virtual tmp<volScalarField> epsilon() const;

//- Return the Reynolds stress tensor
virtual tmp<volSymmTensorField> R() const;

//- Return the effective stress tensor including the laminar stress
virtual tmp<volSymmTensorField> devRhoReff() const;

//- Return the source term for the momentum equation
virtual tmp<fvVectorMatrix> divDevRhoReff(volVectorField& U) const;

//- Solve the turbulence equations and correct the turbulence viscosity
virtual void correct();

//- Read RASProperties dictionary
virtual bool read();
};

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

} // End namespace RASModels
} // End namespace compressible
} // End namespace Foam

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

#endif

// ************************************************************************* //

SpalartAllmaras.C
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
\\ / O peration |
\\ / A nd | Copyright (C) 2011-2013 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 "SpalartAllmaras.H"
#include "addToRunTimeSelectionTable.H"

#include "backwardsCompatibilityWallFunctions.H"
#include "fvcDdt.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

namespace Foam
{
namespace compressible
{
namespace RASModels
{

// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //

defineTypeNameAndDebug(SpalartAllmaras, 0);
addToRunTimeSelectionTable(RASModel, SpalartAllmaras, dictionary);

// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //

tmp<volScalarField> SpalartAllmaras::chi() const
{
return rho_*nuTilda_/mu();
}

tmp<volScalarField> SpalartAllmaras::fv1(const volScalarField& chi) const
{
const volScalarField chi3(pow3(chi));
return chi3/(chi3 + pow3(Cv1_));
}

tmp<volScalarField> SpalartAllmaras::fv2
(
const volScalarField& chi,
const volScalarField& fv1
) const
{
if (ashfordCorrection_)
{
return 1.0/pow3(scalar(1) + chi/Cv2_);
}
else
{
return 1.0 - chi/(1.0 + chi*fv1);
}
}

tmp<volScalarField> SpalartAllmaras::fv3
(
const volScalarField& chi,
const volScalarField& fv1
) const
{
if (ashfordCorrection_)
{
const volScalarField chiByCv2((1/Cv2_)*chi);

return
(scalar(1) + chi*fv1)
*(1/Cv2_)
*(3*(scalar(1) + chiByCv2) + sqr(chiByCv2))
/pow3(scalar(1) + chiByCv2);
}
else
{
return tmp<volScalarField>
(
new volScalarField
(
IOobject
(
"fv3",
mesh_.time().timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh_,
dimensionedScalar("fv3", dimless, 1),
zeroGradientFvPatchScalarField::typeName
)
);
}
}

tmp<volScalarField> SpalartAllmaras::fw(const volScalarField& Stilda) const
{
volScalarField r
(
min
(
nuTilda_
/(
max
(
Stilda,
dimensionedScalar("SMALL", Stilda.dimensions(), SMALL)
)
*sqr(kappa_*d_)
),
scalar(10.0)
)
);
r.boundaryField() == 0.0;

const volScalarField g(r + Cw2_*(pow6(r) - r));

return g*pow((1.0 + pow6(Cw3_))/(pow6(g) + pow6(Cw3_)), 1.0/6.0);
}

/*
* calculation of rotation/curvature correction
* returns 1.0 when spalartShurCorrection = false
*/
tmp<volScalarField> SpalartAllmaras::fr1(const volSymmTensorField& S, const volTensorField& W) const
{
if (spalartShurCorrection_)
{
volScalarField sqrS(2.0*magSqr(S));
volScalarField sqrW(2.0*magSqr(W));
volScalarField sqrD(0.5*(sqrS + sqrW));

volScalarField rStar(sqrt(sqrS/max(sqrW, dimensionedScalar("SMALL", sqrW.dimensions(), SMALL))));

volScalarField Wxx(W.component(tensor::XX));
volScalarField Wxy(W.component(tensor::XY));
volScalarField Wxz(W.component(tensor::XZ));
volScalarField Wyx(W.component(tensor::YX));
volScalarField Wyy(W.component(tensor::YY));
volScalarField Wyz(W.component(tensor::YZ));
volScalarField Wzx(W.component(tensor::ZX));
volScalarField Wzy(W.component(tensor::ZY));
volScalarField Wzz(W.component(tensor::ZZ));

volScalarField rTilda
(
2.0/(rho_*sqr(max(sqrD, dimensionedScalar("SMALL", sqrD.dimensions(), SMALL))))
*(
(Wxx*Sxx_ + Wxy*Sxy_ + Wxz*Sxz_)
*(fvc::ddt(rho_, Sxx_) + fvc::div(phi_, Sxx_)) //i,j=1,1
+ (Wxx*Syx_ + Wxy*Syy_ + Wxz*Syz_)
*(fvc::ddt(rho_, Sxy_) + fvc::div(phi_, Sxy_)) //i,j=1,2
+ (Wxx*Szx_ + Wxy*Szy_ + Wxz*Szz_)
*(fvc::ddt(rho_, Sxz_) + fvc::div(phi_, Sxz_)) //i,j=1,3
+ (Wyx*Sxx_ + Wyy*Sxy_ + Wyz*Sxz_)
*(fvc::ddt(rho_, Syx_) + fvc::div(phi_, Syx_)) //i,j=2,1
+ (Wyx*Syx_ + Wyy*Syy_ + Wyz*Syz_)
*(fvc::ddt(rho_, Syy_) + fvc::div(phi_, Syy_)) //i,j=2,2
+ (Wyx*Szx_ + Wyy*Szy_ + Wyz*Szz_)
*(fvc::ddt(rho_, Syz_) + fvc::div(phi_, Syz_)) //i,j=2,3
+ (Wzx*Sxx_ + Wzy*Sxy_ + Wzz*Sxz_)
*(fvc::ddt(rho_, Szx_) + fvc::div(phi_, Szx_)) //i,j=3,1
+ (Wzx*Syx_ + Wzy*Syy_ + Wzz*Syz_)
*(fvc::ddt(rho_, Szy_) + fvc::div(phi_, Szy_)) //i,j=3,2
+ (Wzx*Szx_ + Wzy*Szy_ + Wzz*Szz_)
*(fvc::ddt(rho_, Szz_) + fvc::div(phi_, Szz_)) //i,j=3,3
)
);

return
(1 + Cr1_)*2.0*rStar/(1.0 + rStar)*(1.0 - Cr3_*atan(Cr2_*rTilda))
- Cr1_;
}
else
{
return tmp<volScalarField>
(
new volScalarField
(
IOobject
(
"fr1",
mesh_.time().timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh_,
dimensionedScalar("fr1", dimless, 1.0),
zeroGradientFvPatchScalarField::typeName
)
);
}
}

/*
* calculation of ft2
* returns 0.0 when negativeNuTilda = false
*/
tmp<volScalarField> SpalartAllmaras::ft2(const volScalarField& chi) const
{
if(negativeNuTilda_)
{
return Ct3_*exp(-Ct4_*sqr(chi));
}
else
{
return tmp<volScalarField>
(
new volScalarField
(
IOobject
(
"ft2",
mesh_.time().timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh_,
dimensionedScalar("ft2", dimless, 0),
zeroGradientFvPatchScalarField::typeName
)
);
}
}

// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //

SpalartAllmaras::SpalartAllmaras
(
const volScalarField& rho,
const volVectorField& U,
const surfaceScalarField& phi,
const fluidThermo& thermophysicalModel,
const word& turbulenceModelName,
const word& modelName
)
:
RASModel(modelName, rho, U, phi, thermophysicalModel, turbulenceModelName),

sigmaNut_
(
dimensioned<scalar>::lookupOrAddToDict
(
"sigmaNut",
coeffDict_,
0.66666
)
),
kappa_
(
dimensioned<scalar>::lookupOrAddToDict
(
"kappa",
coeffDict_,
0.41
)
),
Prt_
(
dimensioned<scalar>::lookupOrAddToDict
(
"Prt",
coeffDict_,
1.0
)
),

Cb1_
(
dimensioned<scalar>::lookupOrAddToDict
(
"Cb1",
coeffDict_,
0.1355
)
),
Cb2_
(
dimensioned<scalar>::lookupOrAddToDict
(
"Cb2",
coeffDict_,
0.622
)
),
Cw1_(Cb1_/sqr(kappa_) + (1.0 + Cb2_)/sigmaNut_),
Cw2_
(
dimensioned<scalar>::lookupOrAddToDict
(
"Cw2",
coeffDict_,
0.3
)
),
Cw3_
(
dimensioned<scalar>::lookupOrAddToDict
(
"Cw3",
coeffDict_,
2.0
)
),
Cv1_
(
dimensioned<scalar>::lookupOrAddToDict
(
"Cv1",
coeffDict_,
7.1
)
),
Cv2_
(
dimensioned<scalar>::lookupOrAddToDict
(
"Cv2",
coeffDict_,
5.0
)
),
Cr1_
(
dimensioned<scalar>::lookupOrAddToDict
(
"Cr1",
coeffDict_,
1.0
)
),
Cr2_
(
dimensioned<scalar>::lookupOrAddToDict
(
"Cr2",
coeffDict_,
12.0
)
),
Cr3_
(
dimensioned<scalar>::lookupOrAddToDict
(
"Cr3",
coeffDict_,
1.0
)
),
C5_
(
dimensioned<scalar>::lookupOrAddToDict
(
"C5",
coeffDict_,
3.5
)
),
Cn1_
(
dimensioned<scalar>::lookupOrAddToDict
(
"Cn1",
coeffDict_,
16.0
)
),
Cn2_
(
dimensioned<scalar>::lookupOrAddToDict
(
"Cn2",
coeffDict_,
0.7
)
),
Cn3_
(
dimensioned<scalar>::lookupOrAddToDict
(
"Cn3",
coeffDict_,
0.9
)
),
Ct3_
(
dimensioned<scalar>::lookupOrAddToDict
(
"Ct3",
coeffDict_,
1.2
)
),
Ct4_
(
dimensioned<scalar>::lookupOrAddToDict
(
"Ct4",
coeffDict_,
0.5
)
),

ashfordCorrection_(coeffDict_.lookupOrDefault("ashfordCorrection", true)),

compressibilityCorrection_
(
coeffDict_.lookupOrDefault("compressibilityCorrection", false)
),

spalartShurCorrection_
(
coeffDict_.lookupOrDefault("spalartShurCorrection", false)
),

StildaModification_
(
coeffDict_.lookupOrDefault("StildaModification", false)
),

negativeNuTilda_
(
coeffDict_.lookupOrDefault("negativeNuTilda", false)
),

nuTilda_
(
IOobject
(
"nuTilda",
runTime_.timeName(),
mesh_,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh_
),

mut_
(
IOobject
(
"mut",
runTime_.timeName(),
mesh_,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh_
),

alphat_
(
IOobject
(
"alphat",
runTime_.timeName(),
mesh_,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
autoCreateAlphat("alphat", mesh_)
),

Sxx_
(
IOobject
(
"Sxx",
runTime_.timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh_,
dimensionedScalar("Sxx", dimensionSet(0, 0, -1, 0, 0, 0, 0), 0),
zeroGradientFvPatchScalarField::typeName
),

Sxy_
(
IOobject
(
"Sxy",
runTime_.timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh_,
dimensionedScalar("Sxy", dimensionSet(0, 0, -1, 0, 0, 0, 0), 0),
zeroGradientFvPatchScalarField::typeName
),

Sxz_
(
IOobject
(
"Sxz",
runTime_.timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh_,
dimensionedScalar("Sxz", dimensionSet(0, 0, -1, 0, 0, 0, 0), 0),
zeroGradientFvPatchScalarField::typeName
),

Syx_
(
IOobject
(
"Syx",
runTime_.timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh_,
dimensionedScalar("Syx", dimensionSet(0, 0, -1, 0, 0, 0, 0), 0),
zeroGradientFvPatchScalarField::typeName
),

Syy_
(
IOobject
(
"Syy",
runTime_.timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh_,
dimensionedScalar("Syy", dimensionSet(0, 0, -1, 0, 0, 0, 0), 0),
zeroGradientFvPatchScalarField::typeName
),

Syz_
(
IOobject
(
"Syz",
runTime_.timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh_,
dimensionedScalar("Syz", dimensionSet(0, 0, -1, 0, 0, 0, 0), 0),
zeroGradientFvPatchScalarField::typeName
),

Szx_
(
IOobject
(
"Szx",
runTime_.timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh_,
dimensionedScalar("Szx", dimensionSet(0, 0, -1, 0, 0, 0, 0), 0),
zeroGradientFvPatchScalarField::typeName
),

Szy_
(
IOobject
(
"Szy",
runTime_.timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh_,
dimensionedScalar("Szy", dimensionSet(0, 0, -1, 0, 0, 0, 0), 0),
zeroGradientFvPatchScalarField::typeName
),

Szz_
(
IOobject
(
"Szz",
runTime_.timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh_,
dimensionedScalar("Szz", dimensionSet(0, 0, -1, 0, 0, 0, 0), 0),
zeroGradientFvPatchScalarField::typeName
),

d_(mesh_)
{
alphat_ = mut_/Prt_;
alphat_.correctBoundaryConditions();

volSymmTensorField S(symm(fvc::grad(U_)));
Sxx_ = S.component(tensor::XX);
Sxy_ = S.component(tensor::XY);
Sxz_ = S.component(tensor::XZ);
Syx_ = S.component(tensor::YX);
Syy_ = S.component(tensor::YY);
Syz_ = S.component(tensor::YZ);
Szx_ = S.component(tensor::ZX);
Szy_ = S.component(tensor::ZY);
Szz_ = S.component(tensor::ZZ);

printCoeffs();

if (ashfordCorrection_)
{
Info<< " Employing Ashford correction" << endl;
}
if (compressibilityCorrection_)
{
Info<< " Employing compressibility correction" << endl;
}
if (spalartShurCorrection_)
{
Info<< " Employing Rotation/Curvature correction" << endl;
}
if (StildaModification_)
{
Info<< " Enabling new Stilda modification" << endl;
}
if (negativeNuTilda_)
{
Info<< " Enabling negative nuTilda" << endl;
}
}

// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //

/*
* calculation of DnuTildaEff
* modifies DnuTildaEff value when chi < 0
*/
tmp<volScalarField> SpalartAllmaras::DnuTildaEff(const volScalarField& chi) const
{
volScalarField pow3chi(pow(chi, 3));
volScalarField fn
(
pos(chi) + neg(chi)*(Cn1_ + pow3chi)/(Cn1_ - pow3chi)
);

return tmp<volScalarField>
(
new volScalarField
(
"DnuTildaEff",
(
(rho_*nuTilda_*fn + mu())/sigmaNut_
)
)
);
}

tmp<volScalarField> SpalartAllmaras::k() const
{
WarningIn("tmp<volScalarField> SpalartAllmaras::k() const")
<< "Turbulence kinetic energy not defined for Spalart-Allmaras model. "
<< "Returning zero field" << endl;

return tmp<volScalarField>
(
new volScalarField
(
IOobject
(
"k",
runTime_.timeName(),
mesh_
),
mesh_,
dimensionedScalar("0", dimensionSet(0, 2, -2, 0, 0), 0)
)
);
}

tmp<volScalarField> SpalartAllmaras::epsilon() const
{
WarningIn("tmp<volScalarField> SpalartAllmaras::epsilon() const")
<< "Turbulence kinetic energy dissipation rate not defined for "
<< "Spalart-Allmaras model. Returning zero field"
<< endl;

return tmp<volScalarField>
(
new volScalarField
(
IOobject
(
"epsilon",
runTime_.timeName(),
mesh_
),
mesh_,
dimensionedScalar("0", dimensionSet(0, 2, -3, 0, 0), 0)
)
);
}

tmp<volSymmTensorField> SpalartAllmaras::R() const
{
return tmp<volSymmTensorField>
(
new volSymmTensorField
(
IOobject
(
"R",
runTime_.timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE
),
((2.0/3.0)*I)*k() - (mut_/rho_)*dev(twoSymm(fvc::grad(U_)))
)
);
}

tmp<volSymmTensorField> SpalartAllmaras::devRhoReff() const
{
return tmp<volSymmTensorField>
(
new volSymmTensorField
(
IOobject
(
"devRhoReff",
runTime_.timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE
),
-muEff()*dev(twoSymm(fvc::grad(U_)))
)
);
}

tmp<fvVectorMatrix> SpalartAllmaras::divDevRhoReff(volVectorField& U) const
{
const volScalarField muEff_(muEff());

return
(
- fvm::laplacian(muEff_, U)
- fvc::div(muEff_*dev2(T(fvc::grad(U))))
);
}

bool SpalartAllmaras::read()
{
if (RASModel::read())
{
sigmaNut_.readIfPresent(coeffDict());
kappa_.readIfPresent(coeffDict());
Prt_.readIfPresent(coeffDict());

Cb1_.readIfPresent(coeffDict());
Cb2_.readIfPresent(coeffDict());
Cw1_ = Cb1_/sqr(kappa_) + (1.0 + Cb2_)/sigmaNut_;
Cw2_.readIfPresent(coeffDict());
Cw3_.readIfPresent(coeffDict());
Cv1_.readIfPresent(coeffDict());
Cv2_.readIfPresent(coeffDict());
Cr1_.readIfPresent(coeffDict());
Cr2_.readIfPresent(coeffDict());
Cr3_.readIfPresent(coeffDict());
C5_.readIfPresent(coeffDict());
Cn1_.readIfPresent(coeffDict());
Cn2_.readIfPresent(coeffDict());
Cn3_.readIfPresent(coeffDict());
Ct3_.readIfPresent(coeffDict());
Ct4_.readIfPresent(coeffDict());

ashfordCorrection_.readIfPresent("ashfordCorrection", coeffDict());
compressibilityCorrection_.readIfPresent
(
"compressibilityCorrection", coeffDict()
);
spalartShurCorrection_.readIfPresent
(
"spalartShurCorrection", coeffDict()
);
StildaModification_.readIfPresent
(
"StildaModification", coeffDict()
);
negativeNuTilda_.readIfPresent
(
"negativeNuTilda", coeffDict()
);

return true;
}
else
{
return false;
}
}

void SpalartAllmaras::correct()
{
if (!turbulence_)
{
// Re-calculate viscosity
mut_ = rho_*nuTilda_*fv1(chi());
mut_.correctBoundaryConditions();

// Re-calculate thermal diffusivity
alphat_ = mut_/Prt_;
alphat_.correctBoundaryConditions();

return;
}

RASModel::correct();

if (mesh_.changing())
{
d_.correct();
}

const volScalarField chi(this->chi());
const volScalarField fv1(this->fv1(chi));

volTensorField gradU(fvc::grad(U_));
volSymmTensorField S(symm(gradU));
volTensorField W(skew(gradU));
volScalarField Omega(sqrt(2.0)*mag(W));

// new Stilda modification used when StildaModification = true
volScalarField Stilda(
IOobject
(
"Stilda",
mesh_.time().timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh_,
dimensionedScalar("Stilda", Omega.dimensions(), 0.0)
);
if (StildaModification_)
{
volScalarField Sbar
(
fv2(chi, fv1)*nuTilda_/sqr(kappa_*d_)
);
Stilda = Omega
+ pos(Cn2_*Omega + Sbar)*Sbar
+ neg(Cn2_*Omega + Sbar)
*(Omega*(sqr(Cn2_)*Omega + Cn3_*Sbar))
/max(((Cn3_ - 2.0*Cn2_)*Omega - Sbar),
dimensionedScalar("SMALL", Omega.dimensions(), SMALL));
}
else
{
Stilda = fv3(chi, fv1)*::sqrt(2.0)*mag(W)
+ fv2(chi, fv1)*nuTilda_/sqr(kappa_*d_);
}
/*
const volScalarField Stilda
(
fv3(chi, fv1)*::sqrt(2.0)*mag(skew(fvc::grad(U_)))
+ fv2(chi, fv1)*nuTilda_/sqr(kappa_*d_)
);
*/

Sxx_ = S.component(tensor::XX);
Sxy_ = S.component(tensor::XY);
Sxz_ = S.component(tensor::XZ);
Syx_ = S.component(tensor::YX);
Syy_ = S.component(tensor::YY);
Syz_ = S.component(tensor::YZ);
Szx_ = S.component(tensor::ZX);
Szy_ = S.component(tensor::ZY);
Szz_ = S.component(tensor::ZZ);

// calculation of compressibility correction
volScalarField compCorr
(
IOobject
(
"compCorr",
mesh_.time().timeName(),
mesh_,
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh_,
dimensionedScalar("compCorr", dimensionSet(1,-1,-2,0,0,0,0), 0.0)
);
if(compressibilityCorrection_)
{
volScalarField Cp(thermo().Cp());
volScalarField Cv(thermo().Cv());
volScalarField a(sqrt((Cp/Cv)*(Cp - Cv)*thermo().T()));
compCorr = -C5_*rho_*sqr(nuTilda_/a)*magSqr(gradU);
}

tmp<fvScalarMatrix> nuTildaEqn
(
fvm::ddt(rho_, nuTilda_)
+ fvm::div(phi_, nuTilda_)
// now DnuTildaEff is dependent on chi
- fvm::laplacian(DnuTildaEff(chi), nuTilda_)
- Cb2_/sigmaNut_*rho_*magSqr(fvc::grad(nuTilda_))
==
pos(nuTilda_)
*(
// RHS for nuTilda >= 0, rotation/curvature correction factor is multiplied
Cb1_*(1.0 - ft2(chi))*rho_*Stilda*nuTilda_*fr1(S, W)
- fvm::Sp((Cw1_*fw(Stilda) - Cb1_/sqr(kappa_)*ft2(chi))*nuTilda_*rho_/sqr(d_), nuTilda_)
)
+ neg(nuTilda_)
*(
// RHS for nuTilda < 0
Cb1_*(1.0 - Ct3_)*rho_*Omega*nuTilda_
+ fvm::Sp(Cw1_*nuTilda_*rho_/sqr(d_), nuTilda_)
)
// add compressibility correction term
+ compCorr
);

nuTildaEqn().relax();
solve(nuTildaEqn);
if(!negativeNuTilda_)
{
// bound nuTilda when negativeNuTilda = false
bound(nuTilda_, dimensionedScalar("0", nuTilda_.dimensions(), 0.0));
}
nuTilda_.correctBoundaryConditions();

// Re-calculate viscosity
mut_.internalField() = mag(fv1)*nuTilda_.internalField()*rho_.internalField();
if(negativeNuTilda_)
{
// bound mut when negativeNuTilda = true
bound(mut_, dimensionedScalar("0", mut_.dimensions(), 0.0));
}
mut_.correctBoundaryConditions();

// Re-calculate thermal diffusivity
alphat_ = mut_/Prt_;
alphat_.correctBoundaryConditions();
}

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

} // End namespace RASModels
} // End namespace compressible
} // End namespace Foam

// ************************************************************************* //

hobbyFoam

Search This Blog

Negative Spalart Allmaras compressible model for OpenFOAM

Comments

Post a Comment