Merge branch 'feature-appyBoundaryLayer-writeTurbulenceFields' into 'develop'

ENH: applyBoundaryLayer - optionally write turbulence fields

See merge request Development/openfoam!313
This commit is contained in:
Sergio Ferraris
2019-12-16 22:18:38 +00:00

View File

@ -6,7 +6,7 @@
\\/ M anipulation | \\/ M anipulation |
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
Copyright (C) 2011-2016 OpenFOAM Foundation Copyright (C) 2011-2016 OpenFOAM Foundation
Copyright (C) 2015-2016 OpenCFD Ltd. Copyright (C) 2015-2019 OpenCFD Ltd.
------------------------------------------------------------------------------- -------------------------------------------------------------------------------
License License
This file is part of OpenFOAM. This file is part of OpenFOAM.
@ -71,12 +71,10 @@ void correctProcessorPatches
return; return;
} }
typedef GeometricField<Type, fvPatchField, volMesh> volFieldType;
// Not possible to use correctBoundaryConditions on fields as they may // Not possible to use correctBoundaryConditions on fields as they may
// use local info as opposed to the constraint values employed here, // use local info as opposed to the constraint values employed here,
// but still need to update processor patches // but still need to update processor patches
typename volFieldType::Boundary& bf = vf.boundaryFieldRef(); auto& bf = vf.boundaryFieldRef();
forAll(bf, patchi) forAll(bf, patchi)
{ {
@ -296,8 +294,12 @@ int main(int argc, char *argv[])
); );
argList::addBoolOption argList::addBoolOption
( (
"write-nut", "writeTurbulenceFields", // (until 1906 was write-nut)
"Write the turbulence viscosity field" "Write the turbulence fields"
);
argList::addOptionCompat
(
"writeTurbulenceFields", {"write-nut", 1906}
); );
#include "setRootCase.H" #include "setRootCase.H"
@ -319,7 +321,7 @@ int main(int argc, char *argv[])
<< exit(FatalError); << exit(FatalError);
} }
const bool writeNut = args.found("write-nut"); const bool writeTurbulenceFields = args.found("writeTurbulenceFields");
#include "createTime.H" #include "createTime.H"
#include "createNamedMesh.H" #include "createNamedMesh.H"
@ -343,36 +345,37 @@ int main(int argc, char *argv[])
mask.correctBoundaryConditions(); mask.correctBoundaryConditions();
correctProcessorPatches<vector>(U); correctProcessorPatches<vector>(U);
// Retrieve nut from turbulence model if (writeTurbulenceFields)
volScalarField nut(calcNut(mesh, U));
// Blend nut using boundary layer profile
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<scalar>(nut);
if (writeNut)
{ {
// Retrieve nut from turbulence model
volScalarField nut(calcNut(mesh, U));
// Blend nut using boundary layer profile
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<scalar>(nut);
Info<< "Writing nut\n" << endl; Info<< "Writing nut\n" << endl;
nut.write(); 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);
} }
// 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);
if (writeNut) setField(mesh, "nuTilda", nut);
// Write the updated U field // Write the updated U field
Info<< "Writing U\n" << endl; Info<< "Writing U\n" << endl;
U.write(); U.write();