eddy.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) 2015 OpenFOAM Foundation
9 Copyright (C) 2016-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
27\*---------------------------------------------------------------------------*/
28
29#include "eddy.H"
31
32using namespace Foam::constant;
33
34// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
35
36Foam::label Foam::eddy::Gamma2Values[] = {1, 2, 3, 4, 5, 6, 7, 8};
37Foam::UList<Foam::label> Foam::eddy::Gamma2(&Gamma2Values[0], 8);
38int Foam::eddy::debug = 0;
39
40
41// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
42
43bool Foam::eddy::setScales
44(
45 const scalar sigmaX,
46 const label gamma2,
47 const vector& e,
48 const vector& lambda,
49 vector& sigma,
51) const
52{
53 // Static array of gamma^2 vs c2 coefficient (PCR:Table 1)
54 static const scalar gamma2VsC2[8] =
55 {2, 1.875, 1.737, 1.75, 0.91, 0.825, 0.806, 1.5};
56
57 const scalar gamma = Foam::sqrt(scalar(gamma2));
58
59 // c2 coefficient retrieved from array
60 const scalar c2 = gamma2VsC2[gamma2 - 1];
61
62 // Length scale in the largest eigenvalue direction
63 const label d1 = dir1_;
64 const label d2 = (d1 + 1) % 3;
65 const label d3 = (d1 + 2) % 3;
66
67 sigma[d1] = sigmaX;
68
69 // Note: sigma_average = 1/3*(sigma_x + sigma_y + sigma_z)
70 // Substituting for sigma_y = sigma_x/gamma and sigma_z = sigma_y
71 // Other length scales equal, as function of major axis length and gamma
72 sigma[d2] = sigma[d1]/gamma;
73 sigma[d3] = sigma[d2];
74
75 // (PCR:Eq. 13)
76 const vector sigma2(cmptMultiply(sigma, sigma));
77 const scalar slos2 = cmptSum(cmptDivide(lambda, sigma2));
78
79 bool ok = true;
80
81 for (label beta = 0; beta < vector::nComponents; ++beta)
82 {
83 const scalar x = slos2 - 2*lambda[beta]/sigma2[beta];
84
85 if (x < 0)
86 {
87 alpha[beta] = 0;
88 ok = false;
89 }
90 else
91 {
92 // (SST:Eq. 23)
93 alpha[beta] = e[beta]*sqrt(x/(2*c2));
94 }
95 }
96
97 if (debug > 1)
98 {
99 Pout<< "c2:" << c2
100 << ", gamma2:" << gamma2
101 << ", gamma:" << gamma
102 << ", lambda:" << lambda
103 << ", sigma2: " << sigma2
104 << ", slos2: " << slos2
105 << ", sigmaX:" << sigmaX
106 << ", sigma:" << sigma
107 << ", alpha:" << alpha
108 << endl;
109 }
110
111 return ok;
112}
113
114
115// * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * * * //
116
118:
119 patchFaceI_(-1),
120 position0_(Zero),
121 x_(0),
122 sigma_(Zero),
123 alpha_(Zero),
124 Rpg_(tensor::I),
125 c1_(-1),
126 dir1_(0)
127{}
128
129
131(
132 const label patchFaceI,
133 const point& position0,
134 const scalar x,
135 const scalar sigmaX,
136 const symmTensor& R,
138)
139:
140 patchFaceI_(patchFaceI),
141 position0_(position0),
142 x_(x),
143 sigma_(Zero),
144 alpha_(Zero),
145 Rpg_(tensor::I),
146 c1_(-1),
147 dir1_(0)
148{
149 // Principal stresses - eigenvalues returned in ascending order
150 const vector lambda(eigenValues(R));
151
152 // Eddy rotation from principal-to-global axes
153 // - given by the 3 eigenvectors of the Reynold stress tensor as rows in
154 // the result tensor (transposed transformation tensor)
155 // - returned in ascending eigenvalue order
156 Rpg_ = eigenVectors(R, lambda).T();
157
158 if (debug)
159 {
160 // Global->Principal transform = Rgp = Rpg.T()
161 // Rgp & R & Rgp.T() should have eigenvalues on its diagonal and
162 // zeros for all other components
163 Pout<< "Rpg.T() & R & Rpg: " << (Rpg_.T() & R & Rpg_) << endl;
164 }
165
166 // Set the eddy orientation to position of max eigenvalue
167 // (direction of eddy major axis, sigma_x in reference)
168 dir1_ = 2;
169
170 // Random vector of 1's and -1's
171 const vector e(epsilon(rndGen));
172
173 // Set intensities and length scales
174 bool found = false;
175 forAll(Gamma2, i)
176 {
177 // Random length scale ratio, gamma = sigmax/sigmay = sigmax/sigmaz
178 // - using gamma^2 to ease lookup of c2 coefficient
179 const label gamma2 = Gamma2[i];
180
181 if (setScales(sigmaX, gamma2, e, lambda, sigma_, alpha_))
182 {
183 found = true;
184 break;
185 }
186 }
187
188 // Normalisation coefficient (PCR:Eq. 11)
189 // Note: sqrt(10*V)/sqrt(nEddy) applied outside when computing uPrime
190 c1_ = cmptAv(sigma_)/cmptProduct(sigma_)*cmptMin(sigma_);
191
192 if (found)
193 {
194 // Shuffle the gamma^2 values
195 rndGen.shuffle(Gamma2);
196 }
197 else
198 {
199 if (debug)
200 {
201 // If not found typically means that the stress has a repeated
202 // eigenvalue/not covered by the selection of Gamma values, e.g.
203 // as seen by range of applicability on Lumley diagram
205 << "Unable to set eddy intensity for eddy: " << *this
206 << endl;
207 }
208
209 // Remove the influence of this eddy/indicate that its initialisation
210 // failed
211 patchFaceI_ = -1;
212 }
213}
214
215
217:
218 patchFaceI_(e.patchFaceI_),
219 position0_(e.position0_),
220 x_(e.x_),
221 sigma_(e.sigma_),
222 alpha_(e.alpha_),
223 Rpg_(e.Rpg_),
224 c1_(e.c1_),
225 dir1_(e.dir1_)
226{}
227
228
229// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
230
232{
233 // Relative position inside eddy (global system) (PCR:p. 524)
234 const vector r(cmptDivide(xp - position(n), sigma_));
235
236 if (mag(r) >= scalar(1))
237 {
238 return vector::zero;
239 }
240
241 // Relative position inside eddy (eddy principal system)
242 const vector rp(Rpg_.T() & r);
243
244 // Shape function (eddy principal system)
245 const vector q(cmptMultiply(sigma_, vector::one - cmptMultiply(rp, rp)));
246
247 // Fluctuating velocity (eddy principal system) (PCR:Eq. 8)
248 const vector uPrimep(cmptMultiply(q, rp^alpha_));
249
250 // Convert into global system (PCR:Eq. 10)
251 return c1_*(Rpg_ & uPrimep);
252}
253
254
256(
257 const vector& n,
258 Ostream& os
259) const
260{
261 const point p(position(n));
262 os << "v " << p.x() << " " << p.y() << " " << p.z() << nl;
263}
264
265
267(
268 const label pointOffset,
269 const vector& n,
270 Ostream& os
271) const
272{
273 if (patchFaceI_ < 0)
274 {
275 // Invalid eddy
276 return 0;
277 }
278
279 static const label nFaceAxis = 20;
280 static const label nFaceTheta = 22;
281 static const label nEddyPoints = (nFaceAxis - 1)*nFaceTheta + 2;
283
284 static scalar dTheta = mathematical::twoPi/nFaceTheta;
285 static scalar dPhi = mathematical::pi/scalar(nFaceAxis);
286
287 label pointI = pointOffset;
288
289 const vector& s = sigma_;
290
291 // Unit vector
292 vector axisDir(Zero);
293 axisDir[dir1_] = 1;
294
295 const label dir2 = (dir1_ + 1) % 3;
296 const label dir3 = (dir1_ + 2) % 3;
297
298 // Calculate the point positions
299 x[0] = axisDir*s[dir1_];
300 x[nEddyPoints - 1] = - axisDir*s[dir1_];
301
302 label eddyPtI = 1;
303 for (label axisI = 1; axisI < nFaceAxis; ++axisI)
304 {
305 const scalar z = s[dir1_]*cos(axisI*dPhi);
306 const scalar r = sqrt(sqr(s[dir2])*(1 - sqr(z)/sqr(s[dir1_])));
307
308 for (label thetaI = 0; thetaI < nFaceTheta; ++thetaI)
309 {
310 const scalar theta = thetaI*dTheta;
311 point& p = x[eddyPtI++];
312 p[dir1_] = z;
313 p[dir2] = r*sin(theta);
314 p[dir3] = r*cos(theta);
315 }
316 }
317
318 // Write points
319 forAll(x, i)
320 {
321 const point p = position(n) + (Rpg_ & x[i]);
322 os << "v " << p.x() << " " << p.y() << " " << p.z() << nl;
323 }
324
325 // Write the end cap tri faces
326 for (label faceI = 0; faceI < nFaceTheta; ++faceI)
327 {
328 const label p1 = pointI + 1;
329 const label p2 = p1 + faceI + 1;
330 label p3 = p2 + 1;
331 if (faceI == nFaceTheta - 1) p3 -= nFaceTheta;
332 os << "f " << p1 << " " << p2 << " " << p3 << nl;
333
334 const label q1 = pointI + nEddyPoints;
335 const label q2 = q1 - faceI - 1;
336 label q3 = q2 - 1;
337 if (faceI == nFaceTheta - 1) q3 += nFaceTheta;
338 os << "f " << q1 << " " << q2 << " " << q3 << nl;
339 }
340
341 // Write quad faces
342 for (label axisI = 1; axisI < nFaceAxis - 1; ++axisI)
343 {
344 for (label thetaI = 0; thetaI < nFaceTheta; ++thetaI)
345 {
346 const label p1 = pointI + 1 + (axisI - 1)*nFaceTheta + thetaI + 1;
347 const label p2 = p1 + nFaceTheta;
348 label p3 = p2 + 1;
349 label p4 = p1 + 1;
350
351 if (thetaI == nFaceTheta - 1)
352 {
353 p3 -= nFaceTheta;
354 p4 -= nFaceTheta;
355 }
356 os << "f " << p1 << " " << p2 << " " << p3 << " " << p4 << nl;
357 }
358 }
359
360 return nEddyPoints;
361}
362
363
364// ************************************************************************* //
#define R(A, B, C, D, E, F, K, M)
bool found
label n
A 1D vector of objects of type <T> with a fixed length <N>.
Definition: FixedList.H:81
An Ostream is an abstract base class for all output systems (streams, files, token lists,...
Definition: Ostream.H:62
Random number generator.
Definition: Random.H:60
void shuffle(UList< Type > &values)
Shuffle the values in the list.
A 1D vector of objects of type <T>, where the size of the vector is known and can be used for subscri...
Definition: UList.H:94
dimensioned< Type > T() const
Return transpose.
Class to describe eddies for the turbulentDFSEMInletFvPatchVectorField boundary condition.
Definition: eddy.H:70
eddy()
Construct null.
Definition: eddy.C:117
void writeCentreOBJ(const vector &n, Ostream &os) const
Write the eddy centre in OBJ format.
Definition: eddy.C:256
const vector & alpha() const noexcept
Return the time-averaged intensity.
Definition: eddyI.H:66
vector uPrime(const point &xp, const vector &n) const
Return the fluctuating velocity contribution at local point xp.
Definition: eddy.C:231
label writeSurfaceOBJ(const label pointOffset, const vector &n, Ostream &os) const
Write the eddy surface in OBJ format.
Definition: eddy.C:267
static int debug
Flag to activate debug statements.
Definition: eddy.H:147
scalar x() const noexcept
Return the distance from the reference position.
Definition: eddyI.H:54
const vector & sigma() const noexcept
Return the length scales in 3-D space.
Definition: eddyI.H:60
A class representing the concept of 1 (one) that can be used to avoid manipulating objects known to b...
Definition: one.H:62
static constexpr direction nComponents
Number of components in bool is 1.
Definition: bool.H:98
A Vector of values with scalar precision, where scalar is float/double depending on the compilation f...
volScalarField & p
const scalar gamma
Definition: EEqn.H:9
regionProperties rp(runTime)
Foam::UList< Foam::label > Foam::eddy::Gamma2 &[0] Gamma2Values
Definition: eddy.C:37
scalar epsilon
OBJstream os(runTime.globalPath()/outputName)
gmvFile<< "tracers "<< particles.size()<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().x()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().y()<< " ";}gmvFile<< nl;for(const passiveParticle &p :particles){ gmvFile<< p.position().z()<< " ";}gmvFile<< nl;forAll(lagrangianScalarNames, i){ word name=lagrangianScalarNames[i];IOField< scalar > s(IOobject(name, runTime.timeName(), cloud::prefix, mesh, IOobject::MUST_READ, IOobject::NO_WRITE))
#define WarningInFunction
Report a warning using Foam::Warning.
constexpr scalar pi(M_PI)
constexpr scalar twoPi(2 *M_PI)
const dimensionedScalar c2
Second radiation constant: default SI units: [m.K].
Different types of constants.
dimensionedTensor eigenVectors(const dimensionedSymmTensor &dt)
dimensionedSymmTensor sqr(const dimensionedVector &dv)
dimensionedVector eigenValues(const dimensionedSymmTensor &dt)
dimensionedScalar sin(const dimensionedScalar &ds)
void cmptMin(FieldField< Field, typename FieldField< Field, Type >::cmptType > &cf, const FieldField< Field, Type > &f)
dimensioned< Type > cmptDivide(const dimensioned< Type > &, const dimensioned< Type > &)
static const Identity< scalar > I
Definition: Identity.H:94
Cmpt cmptProduct(const VectorSpace< Form, Cmpt, Ncmpts > &vs)
Definition: VectorSpaceI.H:598
Cmpt cmptSum(const SphericalTensor< Cmpt > &st)
Return the sum of components of a SphericalTensor.
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)
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
dimensioned< Type > cmptMultiply(const dimensioned< Type > &, const dimensioned< Type > &)
prefixOSstream Pout
OSstream wrapped stdout (std::cout) with parallel prefix.
tmp< DimensionedField< typename DimensionedField< Type, GeoMesh >::cmptType, GeoMesh > > cmptAv(const DimensionedField< Type, GeoMesh > &df)
dimensionedScalar cos(const dimensionedScalar &ds)
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
volScalarField & alpha
volScalarField & e
Definition: createFields.H:11
dimensionedScalar beta("beta", dimless/dimTemperature, laminarTransport)
dimensionedScalar lambda("lambda", dimTime/sqr(dimLength), laminarTransport)
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333
A non-counting (dummy) refCount.
Definition: refCount.H:59
Random rndGen
Definition: createFields.H:23