Skip to content
Snippets Groups Projects
Commit b846422e authored by Andrew Heather's avatar Andrew Heather
Browse files

ENH: Updates to applyBoundaryLayer utility 

Old:
- Previous versions created k and epsilon fields by default, and
  then processed omega and nuTilda fields if present.
- Depending on the choice of turbulence model, not all of these fields
  would be used, and could lead to errors when running some utilities
  due to erroneous values.
- If the omega field did not exist, it would be derived from the epsilon
  field, and also inherit the epsilon boundary conditions (wall
  functions)

New:
- This version will only update fields that already exist on file, i.e.
  will not generate any new fields, and will preserve the boundary
  conditions
parent a2c69ae4
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
applications/utilities/preProcessing/applyBoundaryLayer/applyBoundaryLayer.C
+ 123
− 173
View file @ b846422e
......@@ -49,7 +49,7 @@ Description
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
// turbulence constants - file-scope
// Turbulence constants - file-scope
static const scalar Cmu(0.09);
static const scalar kappa(0.41);
......@@ -110,75 +110,49 @@ void correctProcessorPatches(volScalarField& vf)
}
template<class TurbulenceModel>
Foam::tmp<Foam::volScalarField> calcK
void blendField
(
TurbulenceModel& turbulence,
const volScalarField& mask,
const volScalarField& nut,
const volScalarField& y,
const dimensionedScalar& ybl,
const scalar Cmu,
const scalar kappa
)
{
// Turbulence k
tmp<volScalarField> tk = turbulence->k();
volScalarField& k = tk();
scalar ck0 = pow025(Cmu)*kappa;
k = (1 - mask)*k + mask*sqr(nut/(ck0*min(y, ybl)));
k.rename("k");
// Do not correct BC
// - operation may use inconsistent fields wrt these local manipulations
//k.correctBoundaryConditions();
correctProcessorPatches(k);
Info<< "Writing k\n" << endl;
k.write();
return tk;
}
template<class TurbulenceModel>
Foam::tmp<Foam::volScalarField> calcEpsilon
(
TurbulenceModel& turbulence,
const volScalarField& mask,
const volScalarField& k,
const volScalarField& y,
const dimensionedScalar& ybl,
const scalar Cmu,
const scalar kappa
const word& fieldName,
const fvMesh& mesh,
const scalarField& mask,
const scalarField& boundaryLayerField
)
{
// Turbulence epsilon
tmp<volScalarField> tepsilon = turbulence->epsilon();
volScalarField& epsilon = tepsilon();
scalar ce0 = ::pow(Cmu, 0.75)/kappa;
epsilon = (1 - mask)*epsilon + mask*ce0*k*sqrt(k)/min(y, ybl);
epsilon.max(SMALL);
epsilon.rename("epsilon");
IOobject fieldHeader
(
fieldName,
mesh.time().timeName(),
mesh,
IOobject::MUST_READ,
IOobject::NO_WRITE,
false
);
// Do not correct BC
// - operation may use inconsistent fields wrt these local manipulations
// epsilon.correctBoundaryConditions();
correctProcessorPatches(epsilon);
if (fieldHeader.headerOk())
{
volScalarField fld(fieldHeader, mesh);
scalarField& internalField = fld.internalField();
internalField = (1 - mask)*internalField + mask*boundaryLayerField;
fld.max(SMALL);
Info<< "Writing epsilon\n" << endl;
epsilon.write();
// Do not correct BC
// - operation may use inconsistent fields wrt these local
// manipulations
//fld.correctBoundaryConditions();
correctProcessorPatches(fld);
return tepsilon;
Info<< "Writing " << fieldName << nl << endl;
fld.write();
}
}
void calcOmega
void calcOmegaField
(
const fvMesh& mesh,
const volScalarField& mask,
const volScalarField& k,
const volScalarField& epsilon
const fvMesh& mesh,
const scalarField& mask,
const scalarField& kBL,
const scalarField& epsilonBL
)
{
// Turbulence omega
......@@ -195,9 +169,10 @@ void calcOmega
if (omegaHeader.headerOk())
{
volScalarField omega(omegaHeader, mesh);
dimensionedScalar k0("SMALL", k.dimensions(), SMALL);
scalarField& internalField = omega.internalField();
omega = (1 - mask)*omega + mask*epsilon/(Cmu*k + k0);
internalField =
(1 - mask)*internalField + mask*epsilonBL/(Cmu*kBL + SMALL);
omega.max(SMALL);
// Do not correct BC
......@@ -246,114 +221,79 @@ void setField
}
void calcCompressible
tmp<volScalarField> calcNut
(
const fvMesh& mesh,
const volScalarField& mask,
const volVectorField& U,
const volScalarField& y,
const dimensionedScalar& ybl
const volVectorField& U
)
{
const Time& runTime = mesh.time();
autoPtr<fluidThermo> pThermo(fluidThermo::New(mesh));
fluidThermo& thermo = pThermo();
volScalarField rho(thermo.rho());
// Update/re-write phi
#include "compressibleCreatePhi.H"
phi.write();
autoPtr<compressible::turbulenceModel> turbulence
if
(
compressible::turbulenceModel::New
IOobject
(
rho,
U,
phi,
thermo
)
);
// Hack to correct nut
// Note: in previous versions of the code, nut was initialised on
// construction of the turbulence model. This is no longer the
// case for the Templated Turbulence models. The call to correct
// below will evolve the turbulence model equations and update nut,
// whereas only nut update is required. Need to revisit.
turbulence->correct();
tmp<volScalarField> tnut = turbulence->nut();
volScalarField& nut = tnut();
volScalarField S(mag(dev(symm(fvc::grad(U)))));
nut = (1 - mask)*nut + mask*sqr(kappa*min(y, ybl))*::sqrt(2)*S;
// Do not correct BC - wall functions will 'undo' manipulation above
// by using nut from turbulence model
correctProcessorPatches(nut);
nut.write();
tmp<volScalarField> k =
calcK(turbulence, mask, nut, y, ybl, Cmu, kappa);
tmp<volScalarField> epsilon =
calcEpsilon(turbulence, mask, k, y, ybl, Cmu, kappa);
calcOmega(mesh, mask, k, epsilon);
setField(mesh, "nuTilda", nut);
}
void calcIncompressible
(
const fvMesh& mesh,
const volScalarField& mask,
const volVectorField& U,
const volScalarField& y,
const dimensionedScalar& ybl
)
{
const Time& runTime = mesh.time();
// Update/re-write phi
#include "createPhi.H"
phi.write();
singlePhaseTransportModel laminarTransport(U, phi);
autoPtr<incompressible::turbulenceModel> turbulence
(
incompressible::turbulenceModel::New(U, phi, laminarTransport)
);
// Hack to correct nut
// Note: in previous versions of the code, nut was initialised on
// construction of the turbulence model. This is no longer the
// case for the Templated Turbulence models. The call to correct
// below will evolve the turbulence model equations and update nut,
// whereas only nut update is required. Need to revisit.
turbulence->correct();
basicThermo::dictName,
runTime.constant(),
mesh
).headerOk()
)
{
// Compressible
autoPtr<fluidThermo> pThermo(fluidThermo::New(mesh));
fluidThermo& thermo = pThermo();
volScalarField rho(thermo.rho());
tmp<volScalarField> tnut = turbulence->nut();
// Update/re-write phi
#include "compressibleCreatePhi.H"
phi.write();
// Calculate nut - reference nut is calculated by the turbulence model
// on its construction
volScalarField& nut = tnut();
autoPtr<compressible::turbulenceModel> turbulence
(
compressible::turbulenceModel::New
(
rho,
U,
phi,
thermo
)
);
// Hack to correct nut
// Note: in previous versions of the code, nut was initialised on
// construction of the turbulence model. This is no longer the
// case for the Templated Turbulence models. The call to correct
// below will evolve the turbulence model equations and update nut,
// whereas only nut update is required. Need to revisit.
turbulence->correct();
return tmp<volScalarField>(new volScalarField(turbulence->nut()));
}
else
{
// Incompressible
volScalarField S("S", mag(dev(symm(fvc::grad(U)))));
nut = (1 - mask)*nut + mask*sqr(kappa*min(y, ybl))*::sqrt(2)*S;
// Update/re-write phi
#include "createPhi.H"
phi.write();
// Do not correct BC - wall functions will 'undo' manipulation above
// by using nut from turbulence model
correctProcessorPatches(nut);
nut.write();
singlePhaseTransportModel laminarTransport(U, phi);
tmp<volScalarField> k =
calcK(turbulence, mask, nut, y, ybl, Cmu, kappa);
tmp<volScalarField> epsilon =
calcEpsilon(turbulence, mask, k, y, ybl, Cmu, kappa);
calcOmega(mesh, mask, k, epsilon);
setField(mesh, "nuTilda", nut);
autoPtr<incompressible::turbulenceModel> turbulence
(
incompressible::turbulenceModel::New(U, phi, laminarTransport)
);
// Hack to correct nut
// Note: in previous versions of the code, nut was initialised on
// construction of the turbulence model. This is no longer the
// case for the Templated Turbulence models. The call to correct
// below will evolve the turbulence model equations and update nut,
// whereas only nut update is required. Need to revisit.
turbulence->correct();
return tmp<volScalarField>(new volScalarField(turbulence->nut()));
}
}
......@@ -427,22 +367,32 @@ int main(int argc, char *argv[])
mask.correctBoundaryConditions();
Udash.correctBoundaryConditions();
if
(
IOobject
(
basicThermo::dictName,
runTime.constant(),
mesh
).headerOk()
)
{
calcCompressible(mesh, mask, Udash, y, ybl);
}
else
{
calcIncompressible(mesh, mask, Udash, y, ybl);
}
// Retrieve nut from turbulence model
volScalarField nut(calcNut(mesh, Udash));
// Blend nut using boundary layer profile
volScalarField S("S", mag(dev(symm(fvc::grad(Udash)))));
nut = (1 - mask)*nut + mask*sqr(kappa*min(y, ybl))*::sqrt(2)*S;
// Do not correct BC - wall functions will 'undo' manipulation above
// by using nut from turbulence model
correctProcessorPatches(nut);
nut.write();
// Boundary layer turbulence kinetic energy
scalar ck0 = pow025(Cmu)*kappa;
scalarField kBL(sqr(nut/(ck0*min(y, ybl))));
// Boundary layer turbulence dissipation
scalar ce0 = ::pow(Cmu, 0.75)/kappa;
scalarField epsilonBL(ce0*kBL*sqrt(kBL)/min(y, ybl));
// Process fields if they are present
blendField("k", mesh, mask, kBL);
blendField("epsilon", mesh, mask, epsilonBL);
calcOmegaField(mesh, mask, kBL, epsilonBL);
setField(mesh, "nuTilda", nut);
// Copy internal field Udash into U before writing
Info<< "Writing U\n" << endl;
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment