From a371f1792b1a23167cd91bc3e48179b5e1d867a9 Mon Sep 17 00:00:00 2001 From: Kutalmis Bercin Date: Thu, 7 Apr 2022 11:29:12 +0100 Subject: [PATCH] ENH: setTurbulenceFields: new automatic initialisation method for turbulence fields --- .../setTurbulenceFields/Make/files | 3 + .../setTurbulenceFields/Make/options | 22 + .../setTurbulenceFields/setTurbulenceFields.C | 531 ++++++++++++++++++ 3 files changed, 556 insertions(+) create mode 100644 applications/utilities/preProcessing/setTurbulenceFields/Make/files create mode 100644 applications/utilities/preProcessing/setTurbulenceFields/Make/options create mode 100644 applications/utilities/preProcessing/setTurbulenceFields/setTurbulenceFields.C diff --git a/applications/utilities/preProcessing/setTurbulenceFields/Make/files b/applications/utilities/preProcessing/setTurbulenceFields/Make/files new file mode 100644 index 0000000000..d4a1527323 --- /dev/null +++ b/applications/utilities/preProcessing/setTurbulenceFields/Make/files @@ -0,0 +1,3 @@ +setTurbulenceFields.C + +EXE = $(FOAM_APPBIN)/setTurbulenceFields diff --git a/applications/utilities/preProcessing/setTurbulenceFields/Make/options b/applications/utilities/preProcessing/setTurbulenceFields/Make/options new file mode 100644 index 0000000000..a887778b7b --- /dev/null +++ b/applications/utilities/preProcessing/setTurbulenceFields/Make/options @@ -0,0 +1,22 @@ +EXE_INC = \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude \ + -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \ + -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \ + -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \ + -I$(LIB_SRC)/transportModels \ + -I$(LIB_SRC)/transportModels/compressible/lnInclude \ + -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel + +EXE_LIBS = \ + -lfiniteVolume \ + -lfvOptions \ + -lmeshTools \ + -lturbulenceModels \ + -lincompressibleTurbulenceModels \ + -lcompressibleTurbulenceModels \ + -lfluidThermophysicalModels \ + -lincompressibleTransportModels \ + -lcompressibleTransportModels \ + -lgenericPatchFields diff --git a/applications/utilities/preProcessing/setTurbulenceFields/setTurbulenceFields.C b/applications/utilities/preProcessing/setTurbulenceFields/setTurbulenceFields.C new file mode 100644 index 0000000000..6e69f83457 --- /dev/null +++ b/applications/utilities/preProcessing/setTurbulenceFields/setTurbulenceFields.C @@ -0,0 +1,531 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 . + +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 ; + + // Optional entries + initialiseU ; + initialiseEpsilon ; + initialiseK ; + initialiseOmega ; + initialiseR ; + writeF ; + + kappa ; + Cmu ; + dPlusRef ; + + f ; + U ; + epsilon ; + k ; + omega ; + R ; + \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 "wallFvPatch.H" +#include "fixedValueFvPatchFields.H" + +using namespace Foam; + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +void InfoField(const word& fldName) +{ + Info<< "Writing field: " << fldName << nl << endl; +} + + +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 nu +( + const fvMesh& mesh, + const volVectorField& U +) +{ + const Time& runTime = mesh.time(); + + if + ( + IOobject + ( + basicThermo::dictName, + runTime.constant(), + mesh + ).typeHeaderOk(true) + ) + { + // Compressible + autoPtr pThermo(fluidThermo::New(mesh)); + fluidThermo& thermo = pThermo(); + volScalarField rho(thermo.rho()); + + // Create phi + #include "compressibleCreatePhi.H" + + autoPtr turbulence + ( + compressible::turbulenceModel::New + ( + rho, + U, + phi, + thermo + ) + ); + + turbulence->validate(); + + return tmp::New(turbulence->nu()); + } + else + { + // Incompressible + + // Create phi + #include "createPhi.H" + + singlePhaseTransportModel laminarTransport(U, phi); + + autoPtr turbulence + ( + incompressible::turbulenceModel::New(U, phi, laminarTransport) + ); + + turbulence->validate(); + + return tmp::New(turbulence->nu()); + } +} + + +// Calculate elliptic blending function +// between near-wall and weakly-inhomogeneous regions +void calcF +( + const volScalarField& L, + volScalarField& f +) +{ + tmp tinvLsqr = scalar(1)/sqr(L); + const volScalarField& invLsqr = tinvLsqr.cref(); + + // (M:Eq. 6) + tmp fEqn + ( + fvm::Sp(invLsqr, f) + - fvm::laplacian(f) + == + invLsqr + ); + + tinvLsqr.clear(); + + fEqn.ref().relax(); + solve(fEqn); + + // (M:p. 2) + const dimensioned 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("uRef", scalarMinMax::ge(SMALL)); + const dimensionedScalar uTau(dimVelocity, 0.05*uRef); + const scalar kappa = + dict.getCheckOrDefault("kappa", 0.41, scalarMinMax::ge(SMALL)); + const scalar Cmu = + dict.getCheckOrDefault("Cmu", 0.09, scalarMinMax::ge(SMALL)); + const scalar dPlusRef = + dict.getCheckOrDefault("dPlusRef", 15, scalarMinMax::ge(SMALL)); + const word fName = dict.getOrDefault("f", "f"); + const word UName = dict.getOrDefault("U", "U"); + const word epsilonName = dict.getOrDefault("epsilon", "epsilon"); + const word kName = dict.getOrDefault("k", "k"); + const word omegaName = dict.getOrDefault("omega", "omega"); + const word RName = dict.getOrDefault("R", "R"); + const bool initU = dict.getOrDefault("initialiseU", false); + const bool initEpsilon = + dict.getOrDefault("initialiseEpsilon", false); + const bool initK = dict.getOrDefault("initialiseK", false); + const bool initOmega = dict.getOrDefault("initialiseOmega", false); + const bool initR = dict.getOrDefault("initialiseR", false); + const bool writeF = dict.getOrDefault("writeF", false); + + + // Start initialising the operand fields + + // Read operand fields + + auto tU = tmp::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 tepsilon; + tmp tk; + tmp tomega; + tmp tR; + + if (initEpsilon) + { + tepsilon = tmp::New + ( + createIOobject(mesh, epsilonName), + mesh, + dimensionedScalar(sqr(dimLength)/pow3(dimTime), SMALL), + bcTypes + ); + } + + if (initK) + { + tk = tmp::New + ( + createIOobject(mesh, kName), + mesh, + dimensionedScalar(sqr(dimLength/dimTime), SMALL), + bcTypes + ); + } + + if (initOmega) + { + tomega = tmp::New + ( + createIOobject(mesh, omegaName), + mesh, + dimensionedScalar(dimless/dimTime, SMALL), + bcTypes + ); + } + + if (initR) + { + tR = tmp::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::typeName + ); + + for (fvPatchScalarField& pfld : f.boundaryFieldRef()) + { + if (isA(pfld.patch())) + { + pfld == Zero; + } + } + + + // Create auxillary fields for the initialisation + + tmp tnu = nu(mesh, tU.cref()); + const volScalarField& nu = tnu.cref(); + + // (M:p. 2) + tmp tL = scalar(50)*nu/uTau; + const volScalarField& L = tL.cref(); + + calcF(L, f); + + // (M:Eq. 8) + tmp 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)); + } + + if (tepsilon.valid()) + { + tepsilon.ref() = epsilon; + } + + if (tk.valid()) + { + tk.ref() = k; + } + + if (tomega.valid()) + { + const dimensionedScalar k0(sqr(dimLength/dimTime), SMALL); + tomega.ref() = Cmu*epsilon/(k + k0); + } + + if (tR.valid()) + { + volSymmTensorField& R = tR.ref(); + + // (M:Eq. 3) + const volSphericalTensorField Rdiag(k*twoThirdsI); + forAll(R, celli) + { + R[celli] = Rdiag[celli]; + } + } + + + // 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; +} + + +// ************************************************************************* //