exaFOAM/foam-extend4.1-coherent-io
Archived
1
1
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity

FEATURE: Coupled kEpsilon and kOmegaSST turbulence models. Author: Hrvoje Jasak. Merge: Dominik Christ.

Browse source
This commit is contained in:
Dominik Christ 2015-05-12 16:17:36 +01:00
commit 74a0f0747d
5 changed files with 797 additions and 19 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

1
src/turbulenceModels/incompressible/RAS/Make/files
View file

@ -19,6 +19,7 @@ LienLeschzinerLowRe/LienLeschzinerLowRe.C
LamBremhorstKE/LamBremhorstKE.C
coupledKEpsilon/coupledKEpsilon.C
coupledKOmegaSST/coupledKOmegaSST.C
/* Wall functions */
wallFunctions = derivedFvPatchFields/wallFunctions

46
src/turbulenceModels/incompressible/RAS/coupledKEpsilon/coupledKEpsilon.C
View file

@ -127,11 +127,11 @@ coupledKEpsilon::coupledKEpsilon
),
autoCreateNut("nut", mesh_)
),
ke_
kEpsilon_
(
IOobject
(
"ke",
"kEpsilon",
runTime_.timeName(),
U_.db(),
IOobject::NO_READ,
@ -237,9 +237,9 @@ void coupledKEpsilon::correct()
return;
}
// Make coupled matrix
fvBlockMatrix<vector2> keEqn(kEpsilon_);
fvBlockMatrix<vector2> keEqn(ke_);
volScalarField G("RASModel::G", nut_*2*magSqr(symm(fvc::grad(U_))));
// Update epsilon and G at the wall
@ -253,15 +253,37 @@ void coupledKEpsilon::correct()
+ fvm::div(phi_, epsilon_)
+ fvm::SuSp(-fvc::div(phi_), epsilon_)
- fvm::laplacian(DepsilonEff(), epsilon_)
+ fvm::Sp(C2_*epsilon_/k_, epsilon_)
==
C1_*G*epsilon_/k_
- fvm::Sp(C2_*epsilon_/k_, epsilon_)
);
epsEqn.relax();
epsEqn.completeAssembly();
keEqn.insertEquation(1, epsEqn);
// Add coupling term:
// G_epsilon = C1*Cmu*(symm(grad(U))) k
// but with wall function corrections: must be calculated from G
// HJ, 27/Apr/2015
volScalarField coupling
(
"coupling",
-C1_*G*epsilon_/sqr(k_)
);
scalarField& couplingIn = coupling.internalField();
// Eliminate coupling in wall function cells
labelList eliminated = epsEqn.eliminatedEqns().toc();
forAll (eliminated, cellI)
{
couplingIn[eliminated[cellI]] *= 0;
}
// Insert coupling
keEqn.insertEquationCoupling(1, 0, coupling);
}
// Turbulent kinetic energy equation
@ -272,9 +294,7 @@ void coupledKEpsilon::correct()
+ fvm::div(phi_, k_)
+ fvm::SuSp(-fvc::div(phi_), k_)
- fvm::laplacian(DkEff(), k_)
==
G
- fvm::Sp(epsilon_/k_, k_)
+ fvm::SuSp((epsilon_ - G)/k_, k_)
);
kEqn.relax();
@ -282,16 +302,6 @@ void coupledKEpsilon::correct()
keEqn.insertEquation(0, kEqn);
}
// Add coupling term: C1*Cmu*(symm(grad(U))) k but with wall function
// corrections: must be calculated from G. HJ, 27/Apr/2015
// Add coupling term: epsilon source depends on k
// k, e sink terms cannot be changed because of boundedness
keEqn.insertEquationCoupling
(
1, 0, -C1_*G*epsilon_/sqr(k_)
);
// Update source coupling: coupling terms eliminated from source
keEqn.updateSourceCoupling();

5
src/turbulenceModels/incompressible/RAS/coupledKEpsilon/coupledKEpsilon.H
View file

@ -39,6 +39,9 @@ Description
}
@endverbatim
Author
Hrvoje Jasak, Wikki Ltd. All rights reserved.
SourceFiles
coupledKEpsilon.C
@ -83,7 +86,7 @@ class coupledKEpsilon
volScalarField nut_;
// Block vector field for k-epsilon
volVector2Field ke_;
volVector2Field kEpsilon_;
public:

469
src/turbulenceModels/incompressible/RAS/coupledKOmegaSST/coupledKOmegaSST.C Normal file
View file

@ -0,0 +1,469 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration |
\\ / A nd | For copyright notice see file Copyright
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of foam-extend.
foam-extend 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.
foam-extend 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 foam-extend. If not, see <http://www.gnu.org/licenses/>.
\*---------------------------------------------------------------------------*/
#include "coupledKOmegaSST.H"
#include "fvBlockMatrix.H"
#include "addToRunTimeSelectionTable.H"
#include "backwardsCompatibilityWallFunctions.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace incompressible
{
namespace RASModels
{
// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
defineTypeNameAndDebug(coupledKOmegaSST, 0);
addToRunTimeSelectionTable(RASModel, coupledKOmegaSST, dictionary);
// * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
tmp<volScalarField> coupledKOmegaSST::F1(const volScalarField& CDkOmega) const
{
volScalarField CDkOmegaPlus = max
(
CDkOmega,
dimensionedScalar("1.0e-10", dimless/sqr(dimTime), 1.0e-10)
);
volScalarField arg1 = min
(
min
(
max
(
(scalar(1)/betaStar_)*sqrt(k_)/(omega_*y_),
scalar(500)*nu()/(sqr(y_)*omega_)
),
(4*alphaOmega2_)*k_/(CDkOmegaPlus*sqr(y_))
),
scalar(10)
);
return tanh(pow4(arg1));
}
tmp<volScalarField> coupledKOmegaSST::F2() const
{
volScalarField arg2 = min
(
max
(
(scalar(2)/betaStar_)*sqrt(k_)/(omega_*y_),
scalar(500)*nu()/(sqr(y_)*omega_)
),
scalar(100)
);
return tanh(sqr(arg2));
}
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
coupledKOmegaSST::coupledKOmegaSST
(
const volVectorField& U,
const surfaceScalarField& phi,
transportModel& lamTransportModel
)
:
RASModel(typeName, U, phi, lamTransportModel),
alphaK1_
(
dimensioned<scalar>::lookupOrAddToDict
(
"alphaK1",
coeffDict_,
0.85034
)
),
alphaK2_
(
dimensioned<scalar>::lookupOrAddToDict
(
"alphaK2",
coeffDict_,
1.0
)
),
alphaOmega1_
(
dimensioned<scalar>::lookupOrAddToDict
(
"alphaOmega1",
coeffDict_,
0.5
)
),
alphaOmega2_
(
dimensioned<scalar>::lookupOrAddToDict
(
"alphaOmega2",
coeffDict_,
0.85616
)
),
gamma1_
(
dimensioned<scalar>::lookupOrAddToDict
(
"gamma1",
coeffDict_,
0.5532
)
),
gamma2_
(
dimensioned<scalar>::lookupOrAddToDict
(
"gamma2",
coeffDict_,
0.4403
)
),
beta1_
(
dimensioned<scalar>::lookupOrAddToDict
(
"beta1",
coeffDict_,
0.075
)
),
beta2_
(
dimensioned<scalar>::lookupOrAddToDict
(
"beta2",
coeffDict_,
0.0828
)
),
betaStar_
(
dimensioned<scalar>::lookupOrAddToDict
(
"betaStar",
coeffDict_,
0.09
)
),
a1_
(
dimensioned<scalar>::lookupOrAddToDict
(
"a1",
coeffDict_,
0.31
)
),
c1_
(
dimensioned<scalar>::lookupOrAddToDict
(
"c1",
coeffDict_,
10.0
)
),
y_(mesh_),
k_
(
IOobject
(
"k",
runTime_.timeName(),
U_.db(),
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
autoCreateK("k", mesh_, U_.db())
),
omega_
(
IOobject
(
"omega",
runTime_.timeName(),
U_.db(),
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
autoCreateOmega("omega", mesh_, U_.db())
),
nut_
(
IOobject
(
"nut",
runTime_.timeName(),
U_.db(),
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
autoCreateNut("nut", mesh_, U_.db())
),
kOmega_
(
IOobject
(
"kOmega",
runTime_.timeName(),
U_.db(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
mesh_,
dimensionedVector2("zero", dimless, vector2::zero)
)
{
bound(omega_, omega0_);
nut_ = a1_*k_/max(a1_*omega_, F2()*sqrt(2.0)*mag(symm(fvc::grad(U_))));
nut_.correctBoundaryConditions();
printCoeffs();
}
// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
tmp<volSymmTensorField> coupledKOmegaSST::R() const
{
return tmp<volSymmTensorField>
(
new volSymmTensorField
(
IOobject
(
"R",
runTime_.timeName(),
U_.db(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
((2.0/3.0)*I)*k_ - nut_*twoSymm(fvc::grad(U_)),
k_.boundaryField().types()
)
);
}
tmp<volSymmTensorField> coupledKOmegaSST::devReff() const
{
return tmp<volSymmTensorField>
(
new volSymmTensorField
(
IOobject
(
"devRhoReff",
runTime_.timeName(),
U_.db(),
IOobject::NO_READ,
IOobject::NO_WRITE
),
-nuEff()*dev(twoSymm(fvc::grad(U_)))
)
);
}
tmp<fvVectorMatrix> coupledKOmegaSST::divDevReff(volVectorField& U) const
{
return
(
- fvm::laplacian(nuEff(), U)
- fvc::div(nuEff()*dev(fvc::grad(U)().T()))
);
}
bool coupledKOmegaSST::read()
{
if (RASModel::read())
{
alphaK1_.readIfPresent(coeffDict());
alphaK2_.readIfPresent(coeffDict());
alphaOmega1_.readIfPresent(coeffDict());
alphaOmega2_.readIfPresent(coeffDict());
gamma1_.readIfPresent(coeffDict());
gamma2_.readIfPresent(coeffDict());
beta1_.readIfPresent(coeffDict());
beta2_.readIfPresent(coeffDict());
betaStar_.readIfPresent(coeffDict());
a1_.readIfPresent(coeffDict());
c1_.readIfPresent(coeffDict());
return true;
}
else
{
return false;
}
}
void coupledKOmegaSST::correct()
{
// Bound in case of topological change
// HJ, 22/Aug/2007
if (mesh_.changing())
{
bound(k_, k0_);
bound(omega_, omega0_);
}
RASModel::correct();
if (!turbulence_)
{
return;
}
if (mesh_.changing())
{
y_.correct();
}
volScalarField S2 = magSqr(symm(fvc::grad(U_)));
volScalarField G("RASModel::G", nut_*2*S2);
// Make coupled matrix
fvBlockMatrix<vector2> kOmegaEqn(kOmega_);
// Update omega and G at the wall
omega_.boundaryField().updateCoeffs();
volScalarField CDkOmega =
(2*alphaOmega2_)*(fvc::grad(k_) & fvc::grad(omega_))/omega_;
volScalarField F1 = this->F1(CDkOmega);
Info<< "Coupled k-omega" << endl;
// Turbulent frequency equation
{
fvScalarMatrix omegaEqn
(
fvm::ddt(omega_)
+ fvm::div(phi_, omega_)
+ fvm::SuSp(-fvc::div(phi_), omega_)
- fvm::laplacian(DomegaEff(F1), omega_)
+ fvm::SuSp
(
beta(F1)*omega_
+ (F1 - scalar(1))*CDkOmega/omega_,
omega_
)
==
gamma(F1)*2*S2
);
omegaEqn.relax();
omegaEqn.completeAssembly();
kOmegaEqn.insertEquation(1, omegaEqn);
// Add coupling term
volScalarField coupling
(
"coupling",
-gamma(F1)*2*S2/k_
);
scalarField& couplingIn = coupling.internalField();
// Eliminate coupling in wall function cells
labelList eliminated = omegaEqn.eliminatedEqns().toc();
forAll (eliminated, cellI)
{
couplingIn[eliminated[cellI]] *= 0;
}
// Insert coupling
kOmegaEqn.insertEquationCoupling(1, 0, coupling);
}
// Turbulent kinetic energy equation
{
fvScalarMatrix kEqn
(
fvm::ddt(k_)
+ fvm::div(phi_, k_)
+ fvm::SuSp(-fvc::div(phi_), k_)
- fvm::laplacian(DkEff(F1), k_)
+ fvm::SuSp
(
betaStar_*omega_
- min(G/k_, c1_*betaStar_*omega_),
k_
)
);
kEqn.relax();
kOmegaEqn.insertEquation(0, kEqn);
}
// Update source coupling: coupling terms eliminated from source
kOmegaEqn.updateSourceCoupling();
kOmegaEqn.solve();
// Retrieve solution
kOmegaEqn.retrieveSolution(0, k_.internalField());
kOmegaEqn.retrieveSolution(1, omega_.internalField());
k_.correctBoundaryConditions();
omega_.correctBoundaryConditions();
bound(k_, k0_);
bound(omega_, omega0_);
// Re-calculate viscosity
nut_ = a1_*k_/max(a1_*omega_, F2()*sqrt(2*S2));
nut_ = min(nut_, nuRatio()*nu());
nut_.correctBoundaryConditions();
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace RASModels
} // End namespace incompressible
} // End namespace Foam
// ************************************************************************* //

295
src/turbulenceModels/incompressible/RAS/coupledKOmegaSST/coupledKOmegaSST.H Normal file
View file

@ -0,0 +1,295 @@
/*---------------------------------------------------------------------------*\
========= |
\\ / F ield | foam-extend: Open Source CFD
\\ / O peration |
\\ / A nd | For copyright notice see file Copyright
\\/ M anipulation |
-------------------------------------------------------------------------------
License
This file is part of foam-extend.
foam-extend 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.
foam-extend 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 foam-extend. If not, see <http://www.gnu.org/licenses/>.
Class
Foam::incompressible::RASModels::coupledKOmegaSST
Description
Implementation of the k-omega-SST turbulence model for incompressible
flows using a block-coupled solution.
Turbulence model described in:
@verbatim
Menter, F., Esch, T.
"Elements of Industrial Heat Transfer Prediction"
16th Brazilian Congress of Mechanical Engineering (COBEM),
Nov. 2001
@endverbatim
Note that this implementation is written in terms of alpha diffusion
coefficients rather than the more traditional sigma (alpha = 1/sigma) so
that the blending can be applied to all coefficuients in a consistent
manner. The paper suggests that sigma is blended but this would not be
consistent with the blending of the k-epsilon and k-omega models.
Also note that the error in the last term of equation (2) relating to
sigma has been corrected.
Wall-functions are applied in this implementation by using equations (14)
to specify the near-wall omega as appropriate.
The blending functions (15) and (16) are not currently used because of the
uncertainty in their origin, range of applicability and that is y+ becomes
sufficiently small blending u_tau in this manner clearly becomes nonsense.
The default model coefficients correspond to the following:
@verbatim
coupledKOmegaSSTCoeffs
{
alphaK1 0.85034;
alphaK2 1.0;
alphaOmega1 0.5;
alphaOmega2 0.85616;
beta1 0.075;
beta2 0.0828;
betaStar 0.09;
gamma1 0.5532;
gamma2 0.4403;
a1 0.31;
c1 10.0;
}
@endverbatim
Author
Hrvoje Jasak, Wikki Ltd. All rights reserved.
SourceFiles
coupledKOmegaSST.C
\*---------------------------------------------------------------------------*/
#ifndef coupledKOmegaSST_H
#define coupledKOmegaSST_H
#include "RASModel.H"
#include "wallDist.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
namespace Foam
{
namespace incompressible
{
namespace RASModels
{
/*---------------------------------------------------------------------------*\
Class kOmega Declaration
\*---------------------------------------------------------------------------*/
class coupledKOmegaSST
:
public RASModel
{
// Private data
// Model coefficients
dimensionedScalar alphaK1_;
dimensionedScalar alphaK2_;
dimensionedScalar alphaOmega1_;
dimensionedScalar alphaOmega2_;
dimensionedScalar gamma1_;
dimensionedScalar gamma2_;
dimensionedScalar beta1_;
dimensionedScalar beta2_;
dimensionedScalar betaStar_;
dimensionedScalar a1_;
dimensionedScalar c1_;
//- Wall distance field
// Note: different to wall distance in parent RASModel
wallDist y_;
// Fields
volScalarField k_;
volScalarField omega_;
volScalarField nut_;
// Block vector field for k-omega
volVector2Field kOmega_;
// Private member functions
tmp<volScalarField> F1(const volScalarField& CDkOmega) const;
tmp<volScalarField> F2() const;
tmp<volScalarField> blend
(
const volScalarField& F1,
const dimensionedScalar& psi1,
const dimensionedScalar& psi2
) const
{
return F1*(psi1 - psi2) + psi2;
}
tmp<volScalarField> alphaK
(
const volScalarField& F1
) const
{
return blend(F1, alphaK1_, alphaK2_);
}
tmp<volScalarField> alphaOmega
(
const volScalarField& F1
) const
{
return blend(F1, alphaOmega1_, alphaOmega2_);
}
tmp<volScalarField> beta
(
const volScalarField& F1
) const
{
return blend(F1, beta1_, beta2_);
}
tmp<volScalarField> gamma
(
const volScalarField& F1
) const
{
return blend(F1, gamma1_, gamma2_);
}
public:
//- Runtime type information
TypeName("coupledKOmegaSST");
// Constructors
//- Construct from components
coupledKOmegaSST
(
const volVectorField& U,
const surfaceScalarField& phi,
transportModel& transport
);
//- Destructor
virtual ~coupledKOmegaSST()
{}
// Member Functions
//- Return the turbulence viscosity
virtual tmp<volScalarField> nut() const
{
return nut_;
}
//- Return the effective diffusivity for k
tmp<volScalarField> DkEff(const volScalarField& F1) const
{
return tmp<volScalarField>
(
new volScalarField("DkEff", alphaK(F1)*nut_ + nu())
);
}
//- Return the effective diffusivity for omega
tmp<volScalarField> DomegaEff(const volScalarField& F1) const
{
return tmp<volScalarField>
(
new volScalarField("DomegaEff", alphaOmega(F1)*nut_ + nu())
);
}
//- Return the turbulence kinetic energy
virtual tmp<volScalarField> k() const
{
return k_;
}
//- Return the turbulence specific dissipation rate
virtual tmp<volScalarField> omega() const
{
return omega_;
}
//- Return the turbulence kinetic energy dissipation rate
virtual tmp<volScalarField> epsilon() const
{
return tmp<volScalarField>
(
new volScalarField
(
IOobject
(
"epsilon",
mesh_.time().timeName(),
mesh_
),
betaStar_*k_*omega_,
omega_.boundaryField().types()
)
);
}
//- Return the Reynolds stress tensor
virtual tmp<volSymmTensorField> R() const;
//- Return the effective stress tensor including the laminar stress
virtual tmp<volSymmTensorField> devReff() const;
//- Return the source term for the momentum equation
virtual tmp<fvVectorMatrix> divDevReff(volVectorField& U) const;
//- Solve the turbulence equations and correct the turbulence viscosity
virtual void correct();
//- Read RASProperties dictionary
virtual bool read();
};
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
} // End namespace RASModels
} // namespace incompressible
} // End namespace Foam
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
#endif
// ************************************************************************* //