applyBoundaryLayer.C
Go to the documentation of this file.
1/*---------------------------------------------------------------------------*\
2 ========= |
3 \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4 \\ / O peration |
5 \\ / A nd | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2011-2016 OpenFOAM Foundation
9 Copyright (C) 2015-2021 OpenCFD Ltd.
10-------------------------------------------------------------------------------
11License
12 This file is part of OpenFOAM.
13
14 OpenFOAM is free software: you can redistribute it and/or modify it
15 under the terms of the GNU General Public License as published by
16 the Free Software Foundation, either version 3 of the License, or
17 (at your option) any later version.
18
19 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22 for more details.
23
24 You should have received a copy of the GNU General Public License
25 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26
27Application
28 applyBoundaryLayer
29
30Group
31 grpPreProcessingUtilities
32
33Description
34 Apply a simplified boundary-layer model to the velocity and
35 turbulence fields based on the 1/7th power-law.
36
37 The uniform boundary-layer thickness is either provided via the -ybl option
38 or calculated as the average of the distance to the wall scaled with
39 the thickness coefficient supplied via the option -Cbl. If both options
40 are provided -ybl is used.
41
42 Compressible modes is automatically selected based on the existence of the
43 "thermophysicalProperties" dictionary required to construct the
44 thermodynamics package.
45
46\*---------------------------------------------------------------------------*/
47
48#include "fvCFD.H"
52#include "wallDist.H"
55
56// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
57
58// Turbulence constants - file-scope
59static const scalar Cmu(0.09);
60static const scalar kappa(0.41);
61
62
63template<class Type>
64void correctProcessorPatches
65(
66 GeometricField<Type, fvPatchField, volMesh>& vf
67)
68{
69 if (!Pstream::parRun())
70 {
71 return;
72 }
73
74 // Not possible to use correctBoundaryConditions on fields as they may
75 // use local info as opposed to the constraint values employed here,
76 // but still need to update processor patches
77 auto& bf = vf.boundaryFieldRef();
78
79 forAll(bf, patchi)
80 {
81 if (isA<processorFvPatchField<Type>>(bf[patchi]))
82 {
83 bf[patchi].initEvaluate();
84 }
85 }
86
87 forAll(bf, patchi)
88 {
89 if (isA<processorFvPatchField<Type>>(bf[patchi]))
90 {
91 bf[patchi].evaluate();
92 }
93 }
94}
95
96
97void blendField
98(
99 const word& fieldName,
100 const fvMesh& mesh,
101 const scalarField& mask,
102 const scalarField& boundaryLayerField
103)
104{
105 IOobject fieldHeader
106 (
107 fieldName,
108 mesh.time().timeName(),
109 mesh,
110 IOobject::MUST_READ,
111 IOobject::NO_WRITE,
112 false
113 );
114
115 if (fieldHeader.typeHeaderOk<volScalarField>(true))
116 {
117 volScalarField fld(fieldHeader, mesh);
118 scalarField& pf = fld.primitiveFieldRef();
119 pf = (1 - mask)*pf + mask*boundaryLayerField;
120 fld.max(SMALL);
121
122 // Do not correct BC
123 // - operation may use inconsistent fields wrt these local
124 // manipulations
125 //fld.correctBoundaryConditions();
126 correctProcessorPatches<scalar>(fld);
127
128 Info<< "Writing " << fieldName << nl << endl;
129 fld.write();
130 }
131}
132
133
134void calcOmegaField
135(
136 const fvMesh& mesh,
137 const scalarField& mask,
138 const scalarField& kBL,
139 const scalarField& epsilonBL
140)
141{
142 // Turbulence omega
143 IOobject omegaHeader
144 (
145 "omega",
146 mesh.time().timeName(),
147 mesh,
148 IOobject::MUST_READ,
149 IOobject::NO_WRITE,
150 false
151 );
152
153 if (omegaHeader.typeHeaderOk<volScalarField>(true))
154 {
155 volScalarField omega(omegaHeader, mesh);
156 scalarField& pf = omega.primitiveFieldRef();
157
158 pf = (1 - mask)*pf + mask*epsilonBL/(Cmu*kBL + SMALL);
159 omega.max(SMALL);
160
161 // Do not correct BC
162 // - operation may use inconsistent fields wrt these local
163 // manipulations
164 // omega.correctBoundaryConditions();
165 correctProcessorPatches<scalar>(omega);
166
167 Info<< "Writing omega\n" << endl;
168 omega.write();
169 }
170}
171
172
173void setField
174(
175 const fvMesh& mesh,
176 const word& fieldName,
177 const volScalarField& value
178)
179{
180 IOobject fldHeader
181 (
182 fieldName,
183 mesh.time().timeName(),
184 mesh,
185 IOobject::MUST_READ,
186 IOobject::NO_WRITE,
187 false
188 );
189
190 if (fldHeader.typeHeaderOk<volScalarField>(true))
191 {
192 volScalarField fld(fldHeader, mesh);
193 fld = value;
194
195 // Do not correct BC
196 // - operation may use inconsistent fields wrt these local
197 // manipulations
198 // fld.correctBoundaryConditions();
199 correctProcessorPatches<scalar>(fld);
200
201 Info<< "Writing " << fieldName << nl << endl;
202 fld.write();
203 }
204}
205
206
207tmp<volScalarField> calcNut
208(
209 const fvMesh& mesh,
210 const volVectorField& U
211)
212{
213 const Time& runTime = mesh.time();
214
215 if
216 (
217 IOobject
218 (
219 basicThermo::dictName,
220 runTime.constant(),
221 mesh
222 ).typeHeaderOk<IOdictionary>(true)
223 )
224 {
225 // Compressible
226 autoPtr<fluidThermo> pThermo(fluidThermo::New(mesh));
227 fluidThermo& thermo = pThermo();
229
230 // Update/re-write phi
231 #include "compressibleCreatePhi.H"
232 phi.write();
233
234 autoPtr<compressible::turbulenceModel> turbulence
235 (
236 compressible::turbulenceModel::New
237 (
238 rho,
239 U,
240 phi,
241 thermo
242 )
243 );
244
245 // Correct nut
246 turbulence->validate();
247
248 return tmp<volScalarField>::New(turbulence->nut());
249 }
250 else
251 {
252 // Incompressible
253
254 // Update/re-write phi
255 #include "createPhi.H"
256 phi.write();
257
258 singlePhaseTransportModel laminarTransport(U, phi);
259
260 autoPtr<incompressible::turbulenceModel> turbulence
261 (
262 incompressible::turbulenceModel::New(U, phi, laminarTransport)
263 );
264
265 // Correct nut
266 turbulence->validate();
267
268 return tmp<volScalarField>::New(turbulence->nut());
269 }
270}
271
272
273int main(int argc, char *argv[])
274{
275 argList::addNote
276 (
277 "Apply a simplified boundary-layer model to the velocity and"
278 " turbulence fields based on the 1/7th power-law."
279 );
280
281 #include "addRegionOption.H"
282
283 argList::addOption
284 (
285 "ybl",
286 "scalar",
287 "Specify the boundary-layer thickness"
288 );
289 argList::addOption
290 (
291 "Cbl",
292 "scalar",
293 "Boundary-layer thickness as Cbl * mean distance to wall"
294 );
295 argList::addBoolOption
296 (
297 "writeTurbulenceFields", // (until 1906 was write-nut)
298 "Write the turbulence fields"
299 );
300 argList::addOptionCompat
301 (
302 "writeTurbulenceFields", {"write-nut", 1906}
303 );
304
305 #include "setRootCase.H"
306
307 if (!args.found("ybl") && !args.found("Cbl"))
308 {
310 << "Neither option 'ybl' or 'Cbl' have been provided to calculate "
311 << "the boundary-layer thickness.\n"
312 << "Please choose either 'ybl' OR 'Cbl'."
313 << exit(FatalError);
314 }
315 else if (args.found("ybl") && args.found("Cbl"))
316 {
318 << "Both 'ybl' and 'Cbl' have been provided to calculate "
319 << "the boundary-layer thickness.\n"
320 << "Please choose either 'ybl' OR 'Cbl'."
321 << exit(FatalError);
322 }
323
324 const bool writeTurbulenceFields = args.found("writeTurbulenceFields");
325
326 #include "createTime.H"
327 #include "createNamedMesh.H"
328 #include "createFields.H"
329
330 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
331
332 // Modify velocity by applying a 1/7th power law boundary-layer
333 // u/U0 = (y/ybl)^(1/7)
334 // assumes U0 is the same as the current cell velocity
335 Info<< "Setting boundary layer velocity" << nl << endl;
336 const scalar yblv = ybl.value();
337 forAll(U, celli)
338 {
339 if ((y[celli] > 0) && (y[celli] <= yblv))
340 {
341 mask[celli] = 1;
342 U[celli] *= ::pow(y[celli]/yblv, (1.0/7.0));
343 }
344 }
345 mask.correctBoundaryConditions();
346 correctProcessorPatches<vector>(U);
347
348 if (writeTurbulenceFields)
349 {
350 // Retrieve nut from turbulence model
351 volScalarField nut(calcNut(mesh, U));
352
353 // Blend nut using boundary layer profile
354 volScalarField S("S", mag(dev(symm(fvc::grad(U)))));
355 nut = (1 - mask)*nut + mask*sqr(kappa*min(y, ybl))*::sqrt(2)*S;
356
357 // Do not correct BC - wall functions will 'undo' manipulation above
358 // by using nut from turbulence model
359 correctProcessorPatches<scalar>(nut);
360
361 Info<< "Writing nut\n" << endl;
362 nut.write();
363
364 // Boundary layer turbulence kinetic energy
365 scalar ck0 = pow025(Cmu)*kappa;
366 scalarField kBL(sqr(nut/(ck0*min(y, ybl))));
367
368 // Boundary layer turbulence dissipation
369 scalar ce0 = ::pow(Cmu, 0.75)/kappa;
370 scalarField epsilonBL(ce0*kBL*sqrt(kBL)/min(y, ybl));
371
372 // Process fields if they are present
373 blendField("k", mesh, mask, kBL);
374 blendField("epsilon", mesh, mask, epsilonBL);
375 calcOmegaField(mesh, mask, kBL, epsilonBL);
376 setField(mesh, "nuTilda", nut);
377 }
378
379 // Write the updated U field
380 Info<< "Writing U\n" << endl;
381 U.write();
382
383 Info<< nl;
384 runTime.printExecutionTime(Info);
385
386 Info<< "End\n" << endl;
387
388 return 0;
389}
390
391
392// ************************************************************************* //
scalar y
Info<< nl<< "Wrote faMesh in vtk format: "<< writer.output().name()<< nl;}{ vtk::lineWriter writer(aMesh.points(), aMesh.edges(), fileName(aMesh.mesh().time().globalPath()/"finiteArea-edges"));writer.writeGeometry();writer.beginCellData(4);writer.writeProcIDs();{ Field< scalar > fld(faMeshTools::flattenEdgeField(aMesh.magLe(), true))
surfaceScalarField & phi
virtual bool write(const token &tok)=0
Write token to stream or otherwise handle it.
bool found(const word &optName) const
Return true if the named option is found.
Definition: argListI.H:178
Basic thermodynamics type based on the use of fitting functions for cp, h, s obtained from the templa...
U
Definition: pEqn.H:72
Creates and initialises the face-flux field phi.
dynamicFvMesh & mesh
engineTime & runTime
Required Variables.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
scalar nut
compressible::turbulenceModel & turbulence
const dimensionedScalar kappa
Coulomb constant: default SI units: [N.m2/C2].
dimensionedSymmTensor dev(const dimensionedSymmTensor &dt)
dimensionedSymmTensor symm(const dimensionedSymmTensor &dt)
dimensionedSymmTensor sqr(const dimensionedVector &dv)
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:82
messageStream Info
Information stream (stdout output on master, null elsewhere)
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
dimensionedScalar sqrt(const dimensionedScalar &ds)
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:130
dimensionedScalar pow025(const dimensionedScalar &ds)
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
surfacesMesh setField(triSurfaceToAgglom)
Foam::argList args(argc, argv)
Info<< "Reading thermophysical properties\n"<< endl;autoPtr< psiReactionThermo > pThermo(psiReactionThermo::New(mesh))
singlePhaseTransportModel laminarTransport(U, phi)
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333