572 lines
14 KiB
C
572 lines
14 KiB
C
/*---------------------------------------------------------------------------*\
|
||
========= |
|
||
\\ / F ield | OpenFOAM: The Open Source CFD Toolbox
|
||
\\ / O peration |
|
||
\\ / A nd | www.openfoam.com
|
||
\\/ M anipulation |
|
||
-------------------------------------------------------------------------------
|
||
Copyright (C) 2022 OpenCFD Ltd.
|
||
-------------------------------------------------------------------------------
|
||
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/>.
|
||
|
||
Application
|
||
setTurbulenceFields
|
||
|
||
Group
|
||
grpPreProcessingUtilities
|
||
|
||
Description
|
||
Initialises turbulence fields according to
|
||
various empirical governing equations.
|
||
|
||
References:
|
||
\verbatim
|
||
Initialisation method (tag:M):
|
||
Manceau, R. (n.d.).
|
||
A two-step automatic initialization
|
||
procedure for RANS computations.
|
||
(Unpublished).
|
||
|
||
Previous-version of the initialisation model (tag:LM):
|
||
Lardeau, S., & Manceau, R. (2014).
|
||
Computations of complex flow configurations using
|
||
a modified elliptic-blending Reynolds-stress model.
|
||
10th International ERCOFTAC Symposium on Engineering
|
||
Turbulence Modelling and Measurements. Marbella, Spain.
|
||
https://hal.archives-ouvertes.fr/hal-01051799
|
||
\endverbatim
|
||
|
||
Usage
|
||
Minimal example by using \c system/setTurbulenceFieldsDict:
|
||
\verbatim
|
||
// Mandatory entries
|
||
uRef <scalar>;
|
||
|
||
// Optional entries
|
||
initialiseU <bool>;
|
||
initialiseEpsilon <bool>;
|
||
initialiseK <bool>;
|
||
initialiseOmega <bool>;
|
||
initialiseR <bool>;
|
||
writeF <bool>;
|
||
|
||
kappa <scalar>;
|
||
Cmu <scalar>;
|
||
dPlusRef <scalar>;
|
||
|
||
f <word>;
|
||
U <word>;
|
||
epsilon <word>;
|
||
k <word>;
|
||
omega <word>;
|
||
R <word>;
|
||
\endverbatim
|
||
|
||
where the entries mean:
|
||
\table
|
||
Property | Description | Type | Reqd | Deflt
|
||
uRef | Reference speed | scalar | yes | -
|
||
initialiseU | Flag to initialise U | bool | no | false
|
||
initialiseEpsilon | Flag to initialise epsilon | bool | no | false
|
||
initialiseK | Flag to initialise k | bool | no | false
|
||
initialiseOmega | Flag to initialise omega | bool | no | false
|
||
initialiseR | Flag to initialise R | bool | no | false
|
||
writeF | Flag to write elliptic-blending field, f | bool | no | false
|
||
kappa | von Karman constant | scalar | no | 0.41
|
||
Cmu | Empirical constant | scalar | no | 0.09
|
||
dPlusRef | Reference dPlus | scalar | no | 15
|
||
f | Name of operand f field | word | no | f
|
||
U | Name of operand U field | word | no | U
|
||
epsilon | Name of operand epsilon field | word | no | epsilon
|
||
k | Name of operand k field | word | no | k
|
||
omega | Name of operand omega field | word | no | omega
|
||
R | Name of operand R field | word | no | R
|
||
\endtable
|
||
|
||
Note
|
||
- Time that the utility applies to is determined by the
|
||
\c controlDict.startFrom and \c controlDict.startTime entries.
|
||
- The utility modifies near-wall fields, hence
|
||
can be more effective for low-Re mesh cases.
|
||
|
||
\*---------------------------------------------------------------------------*/
|
||
|
||
#include "fvCFD.H"
|
||
#include "singlePhaseTransportModel.H"
|
||
#include "turbulentTransportModel.H"
|
||
#include "turbulentFluidThermoModel.H"
|
||
#include "processorFvPatchField.H"
|
||
#include "wallFvPatch.H"
|
||
#include "fixedValueFvPatchFields.H"
|
||
|
||
using namespace Foam;
|
||
|
||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||
|
||
void InfoField(const word& fldName)
|
||
{
|
||
Info<< "Writing field: " << fldName << nl << endl;
|
||
}
|
||
|
||
|
||
template<class Type>
|
||
void correctProcessorPatches
|
||
(
|
||
GeometricField<Type, fvPatchField, volMesh>& vf
|
||
)
|
||
{
|
||
if (!Pstream::parRun())
|
||
{
|
||
return;
|
||
}
|
||
|
||
// Not possible to use correctBoundaryConditions on fields as they may
|
||
// use local info as opposed to the constraint values employed here,
|
||
// but still need to update processor patches
|
||
auto& bf = vf.boundaryFieldRef();
|
||
|
||
forAll(bf, patchi)
|
||
{
|
||
if (isA<processorFvPatchField<Type>>(bf[patchi]))
|
||
{
|
||
bf[patchi].initEvaluate();
|
||
}
|
||
}
|
||
|
||
forAll(bf, patchi)
|
||
{
|
||
if (isA<processorFvPatchField<Type>>(bf[patchi]))
|
||
{
|
||
bf[patchi].evaluate();
|
||
}
|
||
}
|
||
}
|
||
|
||
|
||
IOobject createIOobject
|
||
(
|
||
const fvMesh& mesh,
|
||
const word& name,
|
||
IOobject::readOption rOpt = IOobject::READ_IF_PRESENT
|
||
)
|
||
{
|
||
return
|
||
IOobject
|
||
(
|
||
name,
|
||
mesh.time().timeName(),
|
||
mesh,
|
||
rOpt,
|
||
IOobject::NO_WRITE,
|
||
false // do not register
|
||
);
|
||
}
|
||
|
||
|
||
tmp<volScalarField> nu
|
||
(
|
||
const fvMesh& mesh,
|
||
const volVectorField& U
|
||
)
|
||
{
|
||
const Time& runTime = mesh.time();
|
||
|
||
if
|
||
(
|
||
IOobject
|
||
(
|
||
basicThermo::dictName,
|
||
runTime.constant(),
|
||
mesh
|
||
).typeHeaderOk<IOdictionary>(true)
|
||
)
|
||
{
|
||
// Compressible
|
||
autoPtr<fluidThermo> pThermo(fluidThermo::New(mesh));
|
||
fluidThermo& thermo = pThermo();
|
||
volScalarField rho(thermo.rho());
|
||
|
||
// Create phi
|
||
#include "compressibleCreatePhi.H"
|
||
|
||
autoPtr<compressible::turbulenceModel> turbulence
|
||
(
|
||
compressible::turbulenceModel::New
|
||
(
|
||
rho,
|
||
U,
|
||
phi,
|
||
thermo
|
||
)
|
||
);
|
||
|
||
turbulence->validate();
|
||
|
||
return tmp<volScalarField>::New(turbulence->nu());
|
||
}
|
||
else
|
||
{
|
||
// Incompressible
|
||
|
||
// Create phi
|
||
#include "createPhi.H"
|
||
|
||
singlePhaseTransportModel laminarTransport(U, phi);
|
||
|
||
autoPtr<incompressible::turbulenceModel> turbulence
|
||
(
|
||
incompressible::turbulenceModel::New(U, phi, laminarTransport)
|
||
);
|
||
|
||
turbulence->validate();
|
||
|
||
return tmp<volScalarField>::New(turbulence->nu());
|
||
}
|
||
}
|
||
|
||
|
||
// Calculate elliptic blending function
|
||
// between near-wall and weakly-inhomogeneous regions
|
||
void calcF
|
||
(
|
||
const volScalarField& L,
|
||
volScalarField& f
|
||
)
|
||
{
|
||
tmp<volScalarField> tinvLsqr = scalar(1)/sqr(L);
|
||
const volScalarField& invLsqr = tinvLsqr.cref();
|
||
|
||
// (M:Eq. 6)
|
||
tmp<fvScalarMatrix> fEqn
|
||
(
|
||
fvm::Sp(invLsqr, f)
|
||
- fvm::laplacian(f)
|
||
==
|
||
invLsqr
|
||
);
|
||
|
||
tinvLsqr.clear();
|
||
|
||
fEqn.ref().relax();
|
||
solve(fEqn);
|
||
|
||
// (M:p. 2)
|
||
const dimensioned<scalarMinMax> fMinMax
|
||
(
|
||
dimless,
|
||
scalarMinMax(Zero, scalar(1) - Foam::exp(-scalar(400)/scalar(50)))
|
||
);
|
||
|
||
f.clip(fMinMax);
|
||
}
|
||
|
||
|
||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||
// Main program:
|
||
int main(int argc, char *argv[])
|
||
{
|
||
argList::addNote
|
||
(
|
||
"Sets initial turbulence fields based on"
|
||
" various empirical equations"
|
||
);
|
||
argList::noFunctionObjects();
|
||
argList::addOption("dict", "file", "Alternative setTurbulenceFieldsDict");
|
||
|
||
#include "addRegionOption.H"
|
||
|
||
#include "setRootCase.H"
|
||
#include "createTime.H"
|
||
|
||
const word dictName("setTurbulenceFieldsDict");
|
||
#include "setSystemRunTimeDictionaryIO.H"
|
||
Info<< "Reading " << dictIO.name() << nl << endl;
|
||
IOdictionary dict(dictIO);
|
||
|
||
#include "createNamedMesh.H"
|
||
|
||
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
|
||
|
||
IOstream::defaultPrecision(15);
|
||
|
||
|
||
// Dictionary input (M:p. 2)
|
||
|
||
const scalar uRef = dict.getCheck<scalar>("uRef", scalarMinMax::ge(SMALL));
|
||
const dimensionedScalar uTau(dimVelocity, 0.05*uRef);
|
||
const scalar kappa =
|
||
dict.getCheckOrDefault<scalar>("kappa", 0.41, scalarMinMax::ge(SMALL));
|
||
const scalar Cmu =
|
||
dict.getCheckOrDefault<scalar>("Cmu", 0.09, scalarMinMax::ge(SMALL));
|
||
const scalar dPlusRef =
|
||
dict.getCheckOrDefault<scalar>("dPlusRef", 15, scalarMinMax::ge(SMALL));
|
||
const word fName = dict.getOrDefault<word>("f", "f");
|
||
const word UName = dict.getOrDefault<word>("U", "U");
|
||
const word epsilonName = dict.getOrDefault<word>("epsilon", "epsilon");
|
||
const word kName = dict.getOrDefault<word>("k", "k");
|
||
const word omegaName = dict.getOrDefault<word>("omega", "omega");
|
||
const word RName = dict.getOrDefault<word>("R", "R");
|
||
const bool initU = dict.getOrDefault<bool>("initialiseU", false);
|
||
const bool initEpsilon =
|
||
dict.getOrDefault<bool>("initialiseEpsilon", false);
|
||
const bool initK = dict.getOrDefault<bool>("initialiseK", false);
|
||
const bool initOmega = dict.getOrDefault<bool>("initialiseOmega", false);
|
||
const bool initR = dict.getOrDefault<bool>("initialiseR", false);
|
||
const bool writeF = dict.getOrDefault<bool>("writeF", false);
|
||
|
||
|
||
// Start initialising the operand fields
|
||
|
||
// Read operand fields
|
||
|
||
auto tU = tmp<volVectorField>::New
|
||
(
|
||
createIOobject(mesh, UName, IOobject::MUST_READ),
|
||
mesh
|
||
);
|
||
|
||
// Infer the initial BCs from the velocity
|
||
const wordList bcTypes
|
||
(
|
||
tU.cref().boundaryField().size(),
|
||
fixedValueFvPatchScalarField::typeName
|
||
);
|
||
|
||
tmp<volScalarField> tepsilon;
|
||
tmp<volScalarField> tk;
|
||
tmp<volScalarField> tomega;
|
||
tmp<volSymmTensorField> tR;
|
||
|
||
if (initEpsilon)
|
||
{
|
||
tepsilon = tmp<volScalarField>::New
|
||
(
|
||
createIOobject(mesh, epsilonName),
|
||
mesh,
|
||
dimensionedScalar(sqr(dimLength)/pow3(dimTime), SMALL),
|
||
bcTypes
|
||
);
|
||
}
|
||
|
||
if (initK)
|
||
{
|
||
tk = tmp<volScalarField>::New
|
||
(
|
||
createIOobject(mesh, kName),
|
||
mesh,
|
||
dimensionedScalar(sqr(dimLength/dimTime), SMALL),
|
||
bcTypes
|
||
);
|
||
}
|
||
|
||
if (initOmega)
|
||
{
|
||
tomega = tmp<volScalarField>::New
|
||
(
|
||
createIOobject(mesh, omegaName),
|
||
mesh,
|
||
dimensionedScalar(dimless/dimTime, SMALL),
|
||
bcTypes
|
||
);
|
||
}
|
||
|
||
if (initR)
|
||
{
|
||
tR = tmp<volSymmTensorField>::New
|
||
(
|
||
createIOobject(mesh, RName),
|
||
mesh,
|
||
dimensionedSymmTensor(sqr(dimLength/dimTime), Zero),
|
||
bcTypes
|
||
);
|
||
}
|
||
|
||
|
||
// Create elliptic blending factor field
|
||
|
||
volScalarField f
|
||
(
|
||
IOobject
|
||
(
|
||
fName,
|
||
runTime.timeName(),
|
||
mesh,
|
||
IOobject::NO_READ,
|
||
IOobject::NO_WRITE,
|
||
false // do not register
|
||
),
|
||
mesh,
|
||
dimensionedScalar(dimless, scalar(1)),
|
||
fixedValueFvPatchField<scalar>::typeName
|
||
);
|
||
|
||
for (fvPatchScalarField& pfld : f.boundaryFieldRef())
|
||
{
|
||
if (isA<wallFvPatch>(pfld.patch()))
|
||
{
|
||
pfld == Zero;
|
||
}
|
||
}
|
||
|
||
|
||
// Create auxillary fields for the initialisation
|
||
|
||
tmp<volScalarField> tnu = nu(mesh, tU.cref());
|
||
const volScalarField& nu = tnu.cref();
|
||
|
||
// (M:p. 2)
|
||
tmp<volScalarField> tL = scalar(50)*nu/uTau;
|
||
const volScalarField& L = tL.cref();
|
||
|
||
calcF(L, f);
|
||
|
||
// (M:Eq. 8)
|
||
tmp<volScalarField> td = -tL*Foam::log(scalar(1) - f);
|
||
const volScalarField& d = td.cref();
|
||
|
||
// (M:p. 2)
|
||
const volScalarField dPlus(d*uTau/nu);
|
||
|
||
// (M:Eq. 11)
|
||
const volScalarField epsilon
|
||
(
|
||
pow4(uTau)/(kappa*nu*max(dPlus, dPlusRef))
|
||
);
|
||
|
||
// (M:Eq. 13)
|
||
const volScalarField fk(Foam::exp(-dPlus/scalar(25)));
|
||
|
||
// (M:Eq. 12)
|
||
const volScalarField k
|
||
(
|
||
(epsilon*sqr(td)*sqr(fk))/(2*nu)
|
||
+ sqr(uTau)*sqr(scalar(1) - fk)/Foam::sqrt(Cmu)
|
||
);
|
||
|
||
|
||
// Initialise operand fields
|
||
|
||
if (initU)
|
||
{
|
||
volVectorField& U = tU.ref();
|
||
|
||
// Reichardt’s law (M:Eq. 10)
|
||
const scalar C = 7.8;
|
||
const scalar B1 = 11;
|
||
const scalar B2 = 3;
|
||
const volScalarField fRei
|
||
(
|
||
Foam::log(scalar(1) + kappa*dPlus)/kappa
|
||
+ C*
|
||
(
|
||
scalar(1)
|
||
- Foam::exp(-dPlus/B1)
|
||
- dPlus/B1*Foam::exp(-dPlus/B2)
|
||
)
|
||
);
|
||
|
||
// (M:Eq. 9)
|
||
const dimensionedScalar maxU(dimVelocity, SMALL);
|
||
U *= min(scalar(1), fRei*uTau/max(mag(U), maxU));
|
||
correctProcessorPatches<vector>(U);
|
||
}
|
||
|
||
if (tepsilon.valid())
|
||
{
|
||
tepsilon.ref() = epsilon;
|
||
correctProcessorPatches<scalar>(tepsilon.ref());
|
||
}
|
||
|
||
if (tk.valid())
|
||
{
|
||
tk.ref() = k;
|
||
correctProcessorPatches<scalar>(tk.ref());
|
||
}
|
||
|
||
if (tomega.valid())
|
||
{
|
||
const dimensionedScalar k0(sqr(dimLength/dimTime), SMALL);
|
||
tomega.ref() = Cmu*epsilon/(k + k0);
|
||
correctProcessorPatches<scalar>(tomega.ref());
|
||
}
|
||
|
||
if (tR.valid())
|
||
{
|
||
volSymmTensorField& R = tR.ref();
|
||
|
||
// (M:Eq. 3)
|
||
const volSphericalTensorField Rdiag(k*twoThirdsI);
|
||
forAll(R, celli)
|
||
{
|
||
R[celli] = Rdiag[celli];
|
||
}
|
||
correctProcessorPatches<symmTensor>(R);
|
||
}
|
||
|
||
|
||
// Write operand fields
|
||
|
||
Info<< endl;
|
||
|
||
if (initU)
|
||
{
|
||
InfoField(tU->name());
|
||
tU->write();
|
||
}
|
||
|
||
if (tepsilon.valid())
|
||
{
|
||
InfoField(tepsilon->name());
|
||
tepsilon->write();
|
||
}
|
||
|
||
if (tk.valid())
|
||
{
|
||
InfoField(tk->name());
|
||
tk->write();
|
||
}
|
||
|
||
if (tomega.valid())
|
||
{
|
||
InfoField(tomega->name());
|
||
tomega->write();
|
||
}
|
||
|
||
if (tR.valid())
|
||
{
|
||
InfoField(tR->name());
|
||
tR->write();
|
||
}
|
||
|
||
if (writeF)
|
||
{
|
||
InfoField(f.name());
|
||
f.write();
|
||
}
|
||
|
||
|
||
Info<< nl;
|
||
runTime.printExecutionTime(Info);
|
||
|
||
Info<< "End\n" << endl;
|
||
|
||
return 0;
|
||
}
|
||
|
||
|
||
// ************************************************************************* //
|