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 -Foam::tmp calcK +void blendField ( - TurbulenceModel& turbulence, - const volScalarField& mask, - const volScalarField& nut, - 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 k - tmp 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"); + 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 - //k.correctBoundaryConditions(); - correctProcessorPatches(k); + if (fieldHeader.headerOk()) + { + volScalarField fld(fieldHeader, mesh); + scalarField& internalField = fld.internalField(); + internalField = (1 - mask)*internalField + mask*boundaryLayerField; + fld.max(SMALL); - Info<< "Writing k\n" << endl; - k.write(); + // Do not correct BC + // - operation may use inconsistent fields wrt these local + // manipulations + //fld.correctBoundaryConditions(); + correctProcessorPatches(fld); - return tk; + Info<< "Writing " << fieldName << nl << endl; + fld.write(); + } } -template -Foam::tmp calcEpsilon +void calcOmegaField ( - TurbulenceModel& turbulence, - const volScalarField& mask, - const volScalarField& k, - const volScalarField& y, - const dimensionedScalar& ybl, - const scalar Cmu, - const scalar kappa -) -{ - // Turbulence epsilon - tmp 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"); - - // Do not correct BC - // - operation may use inconsistent fields wrt these local manipulations - // epsilon.correctBoundaryConditions(); - correctProcessorPatches(epsilon); - - Info<< "Writing epsilon\n" << endl; - epsilon.write(); - - return tepsilon; -} - - -void calcOmega -( - 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 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 pThermo(fluidThermo::New(mesh)); - fluidThermo& thermo = pThermo(); - volScalarField rho(thermo.rho()); - - // Update/re-write phi - #include "compressibleCreatePhi.H" - phi.write(); - - autoPtr turbulence + if ( - compressible::turbulenceModel::New + IOobject ( - rho, - U, - phi, - thermo - ) - ); + basicThermo::dictName, + runTime.constant(), + mesh + ).headerOk() + ) + { + // Compressible + autoPtr pThermo(fluidThermo::New(mesh)); + fluidThermo& thermo = pThermo(); + volScalarField rho(thermo.rho()); - // 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(); + // Update/re-write phi + #include "compressibleCreatePhi.H" + phi.write(); - tmp tnut = turbulence->nut(); + autoPtr turbulence + ( + compressible::turbulenceModel::New + ( + rho, + U, + phi, + thermo + ) + ); - volScalarField& nut = tnut(); - volScalarField S(mag(dev(symm(fvc::grad(U))))); - nut = (1 - mask)*nut + mask*sqr(kappa*min(y, ybl))*::sqrt(2)*S; + // 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(); - // Do not correct BC - wall functions will 'undo' manipulation above - // by using nut from turbulence model - correctProcessorPatches(nut); - nut.write(); + return tmp(new volScalarField(turbulence->nut())); + } + else + { + // Incompressible - tmp k = - calcK(turbulence, mask, nut, y, ybl, Cmu, kappa); - tmp epsilon = - calcEpsilon(turbulence, mask, k, y, ybl, Cmu, kappa); - calcOmega(mesh, mask, k, epsilon); - setField(mesh, "nuTilda", nut); -} + // Update/re-write phi + #include "createPhi.H" + phi.write(); + singlePhaseTransportModel laminarTransport(U, phi); -void calcIncompressible -( - const fvMesh& mesh, - const volScalarField& mask, - const volVectorField& U, - const volScalarField& y, - const dimensionedScalar& ybl -) -{ - const Time& runTime = mesh.time(); + autoPtr turbulence + ( + incompressible::turbulenceModel::New(U, phi, laminarTransport) + ); - // Update/re-write phi - #include "createPhi.H" - phi.write(); + // 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(); - singlePhaseTransportModel laminarTransport(U, phi); - - autoPtr 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(); - - tmp tnut = turbulence->nut(); - - // Calculate nut - reference nut is calculated by the turbulence model - // on its construction - volScalarField& nut = tnut(); - - volScalarField S("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 k = - calcK(turbulence, mask, nut, y, ybl, Cmu, kappa); - tmp epsilon = - calcEpsilon(turbulence, mask, k, y, ybl, Cmu, kappa); - calcOmega(mesh, mask, k, epsilon); - setField(mesh, "nuTilda", nut); + return tmp(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;