createBoxTurb.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) 2018 OpenCFD Ltd.
9-------------------------------------------------------------------------------
10License
11 This file is part of OpenFOAM.
12
13 OpenFOAM is free software: you can redistribute it and/or modify it
14 under the terms of the GNU General Public License as published by
15 the Free Software Foundation, either version 3 of the License, or
16 (at your option) any later version.
17
18 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
19 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
20 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
21 for more details.
22
23 You should have received a copy of the GNU General Public License
24 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
25
26Application
27 createBoxTurb
28
29Description
30 Create a box of isotropic turbulence based on a user-specified
31 energy spectrum.
32
33 Based on the reference
34 \verbatim
35 Saad, T., Cline, D., Stoll, R., Sutherland, J.C.
36 "Scalable Tools for Generating Synthetic Isotropic Turbulence with
37 Arbitrary Spectra"
38 AIAA Journal, Vol. 55, No. 1 (2017), pp. 327-331.
39 \endverbatim
40
41 The \c -createBlockMesh option creates a block mesh and exits, which
42 can then be decomposed and the utility run in parallel.
43
44\*---------------------------------------------------------------------------*/
45
46#include "fvCFD.H"
47#include "block.H"
49
50using namespace Foam::constant;
51
52// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
53
54Foam::vector randomUnitVector(Random& rndGen)
55{
56 // Sample point on a sphere
57 scalar t = rndGen.globalPosition<scalar>(-1, 1);
58 scalar phim = rndGen.globalSample01<scalar>()*mathematical::twoPi;
59 scalar thetam = Foam::acos(t);
60
61 return vector
62 (
63 Foam::sin(thetam*Foam::cos(phim)),
64 Foam::sin(thetam*Foam::sin(phim)),
65 Foam::cos(thetam)
66 );
67}
68
69
70int main(int argc, char *argv[])
71{
72 argList::addNote
73 (
74 "Create a box of isotropic turbulence based on a user-specified"
75 " energy spectrum."
76 );
77
78 argList::addBoolOption
79 (
80 "createBlockMesh",
81 "create the block mesh and exit"
82 );
83
84 #include "setRootCase.H"
85
86 #include "createTime.H"
87 #include "createFields.H"
88
89 if (args.found("createBlockMesh"))
90 {
91 // Create a box block mesh with cyclic patches
92 #include "createBlockMesh.H"
93 Info<< "\nEnd\n" << endl;
94 return 0;
95 }
96
97 #include "createMesh.H"
98
99 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
100
101 // Minimum wave number
102 scalar kappa0 = mathematical::twoPi/cmptMin(L);
103
104 // Maximum wave number
105 scalar kappaMax = mathematical::pi/cmptMin(delta);
106
107 Info<< "Wave number min/max = " << kappa0 << ", " << kappaMax << endl;
108
109 Info<< "Generating velocity field" << endl;
110
112 (
114 (
115 "U",
116 runTime.timeName(),
117 mesh,
118 IOobject::NO_READ,
119 IOobject::NO_WRITE
120 ),
121 mesh,
122 dimensionedVector(dimVelocity, Zero)
123 );
124
125 vectorField& Uc = U.primitiveFieldRef();
126 const scalar deltaKappa = (kappaMax - kappa0)/scalar(nModes - 1);
127 const vectorField& C(mesh.C());
128 for (label modei = 1; modei <= nModes; ++modei)
129 {
130 // Equidistant wave mode
131 scalar kappaM = kappa0 + deltaKappa*(modei-1);
132
133 Info<< "Processing mode:" << modei << " kappaM:" << kappaM << endl;
134
135 // Energy
136 scalar E = Ek->value(kappaM);
137
138 // Wave amplitude
139 scalar qm = Foam::sqrt(E*deltaKappa);
140
141 // Wave number unit vector
142 const vector kappaHatm(randomUnitVector(rndGen));
143
144 vector kappaTildem(0.5*kappaM*cmptMultiply(kappaHatm, delta));
145 for (direction i = 0; i < 3; ++i)
146 {
147 kappaTildem[i] = 2/delta[i]*Foam::sin(kappaTildem[i]);
148 }
149
150 // Intermediate unit vector zeta
151 const vector zetaHatm(randomUnitVector(rndGen));
152
153 // Unit vector sigma
154 vector sigmaHatm(zetaHatm^kappaTildem);
155 sigmaHatm /= mag(kappaTildem);
156
157 // Phase angle
158 scalar psim = 0.5*rndGen.position(-mathematical::pi, mathematical::pi);
159
160 // Add the velocity contribution per mode
161 Uc += 2*qm*cos(kappaM*(kappaHatm & C) + psim)*sigmaHatm;
162 }
163
164 U.write();
165
166 {
167 Info<< "Generating kinetic energy field" << endl;
168 volScalarField k("k", 0.5*magSqr(U));
169 k.write();
170 Info<< "min/max/average k = "
171 << gMin(k) << ", " << gMax(k) << ", " << gAverage(k)
172 << endl;
173 }
174
175 {
176 Info<< "Generating div(U) field" << endl;
177 volScalarField divU(fvc::div(U));
178 divU.write();
179 Info<< "min/max/average div(U) = "
180 << gMin(divU) << ", " << gMax(divU) << ", " << gAverage(divU)
181 << endl;
182 }
183
184 Info<< nl;
185 runTime.printExecutionTime(Info);
186
187 Info<< "End\n" << endl;
188 return 0;
189}
190
191
192// ************************************************************************* //
static const Foam::dimensionedScalar C("", Foam::dimTemperature, 234.5)
label k
scalar delta
Graphite solid properties.
Definition: C.H:53
Defines the attributes of an object for which implicit objectRegistry management is supported,...
Definition: IOobject.H:170
Random number generator.
Definition: Random.H:60
bool found(const word &optName) const
Return true if the named option is found.
Definition: argListI.H:178
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
U
Definition: pEqn.H:72
dynamicFvMesh & mesh
engineTime & runTime
Required Variables.
zeroField divU
Definition: alphaSuSp.H:3
Different types of constants.
dimensionedScalar sin(const dimensionedScalar &ds)
void cmptMin(FieldField< Field, typename FieldField< Field, Type >::cmptType > &cf, const FieldField< Field, Type > &f)
messageStream Info
Information stream (stdout output on master, null elsewhere)
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)
uint8_t direction
Definition: direction.H:56
Type gAverage(const FieldField< Field, Type > &f)
dimensioned< Type > cmptMultiply(const dimensioned< Type > &, const dimensioned< Type > &)
Type gMin(const FieldField< Field, Type > &f)
Type gMax(const FieldField< Field, Type > &f)
dimensioned< typename typeOfMag< Type >::type > magSqr(const dimensioned< Type > &dt)
dimensionedScalar cos(const dimensionedScalar &ds)
dimensionedScalar acos(const dimensionedScalar &ds)
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
Foam::argList args(argc, argv)
Random rndGen
Definition: createFields.H:23
const label nModes
Definition: createFields.H:22
const vector L(dict.get< vector >("L"))
autoPtr< Function1< scalar > > Ek(Function1< scalar >::New("Ek", dict, &runTime))