From 78a0b623e76f0e997e87d4d087be9b3981fefcd9 Mon Sep 17 00:00:00 2001
From: Andrew Heather <a.heather@opencfd.co.uk>
Date: Thu, 21 Apr 2016 12:46:23 +0100
Subject: [PATCH] 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
---
.../applyBoundaryLayer/applyBoundaryLayer.C | 296 ++++++++----------
1 file changed, 123 insertions(+), 173 deletions(-)
diff --git a/applications/utilities/preProcessing/applyBoundaryLayer/applyBoundaryLayer.C b/applications/utilities/preProcessing/applyBoundaryLayer/applyBoundaryLayer.C
index 10d5d3b4a3..086c88f6eb 100644
--- a/applications/utilities/preProcessing/applyBoundaryLayer/applyBoundaryLayer.C
+++ b/applications/utilities/preProcessing/applyBoundaryLayer/applyBoundaryLayer.C
@@ -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;
--
GitLab