FreeStream.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-2017 OpenFOAM Foundation
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
26\*---------------------------------------------------------------------------*/
27
28#include "FreeStream.H"
29#include "constants.H"
30#include "triPointRef.H"
31#include "tetIndices.H"
32
33using namespace Foam::constant::mathematical;
34
35// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
36
37template<class CloudType>
39(
40 const dictionary& dict,
42)
43:
45 patches_(),
46 moleculeTypeIds_(),
47 numberDensities_(),
48 particleFluxAccumulators_()
49{
50 // Identify which patches to use
51
53
54 forAll(cloud.mesh().boundaryMesh(), p)
55 {
56 const polyPatch& patch = cloud.mesh().boundaryMesh()[p];
57
58 if (isType<polyPatch>(patch))
59 {
61 }
62 }
63
64 patches_.transfer(patches);
65
66 const dictionary& numberDensitiesDict
67 (
68 this->coeffDict().subDict("numberDensities")
69 );
70
71 List<word> molecules(numberDensitiesDict.toc());
72
73 // Initialise the particleFluxAccumulators_
74 particleFluxAccumulators_.setSize(patches_.size());
75
76 forAll(patches_, p)
77 {
78 const polyPatch& patch = cloud.mesh().boundaryMesh()[patches_[p]];
79
80 particleFluxAccumulators_[p] = List<Field<scalar>>
81 (
82 molecules.size(),
83 Field<scalar>(patch.size(), Zero)
84 );
85 }
86
87 moleculeTypeIds_.setSize(molecules.size());
88
89 numberDensities_.setSize(molecules.size());
90
91 forAll(molecules, i)
92 {
93 numberDensities_[i] = numberDensitiesDict.get<scalar>(molecules[i]);
94
95 moleculeTypeIds_[i] = cloud.typeIdList().find(molecules[i]);
96
97 if (moleculeTypeIds_[i] == -1)
98 {
100 << "typeId " << molecules[i] << "not defined in cloud." << nl
101 << abort(FatalError);
102 }
103 }
104
105 numberDensities_ /= cloud.nParticle();
106}
107
108
109
110// * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
111
112template<class CloudType>
114{}
115
116
117// * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
118
119template<class CloudType>
121{
122 CloudType& cloud(this->owner());
123 const polyMesh& mesh(cloud.mesh());
124
125 forAll(patches_, p)
126 {
127 label patchi = patches_[p];
128
129 const polyPatch& patch = mesh.boundaryMesh()[patchi];
130 List<Field<scalar>>& pFA = particleFluxAccumulators_[p];
131
132 forAll(pFA, facei)
133 {
134 pFA[facei].setSize(patch.size(), 0);
135 }
136 }
137}
138
139
140template<class CloudType>
142{
143 CloudType& cloud(this->owner());
144
145 const polyMesh& mesh(cloud.mesh());
146
147 const scalar deltaT = mesh.time().deltaTValue();
148
149 Random& rndGen = cloud.rndGen();
150
151 scalar sqrtPi = sqrt(pi);
152
153 label particlesInserted = 0;
154
155 const volScalarField::Boundary& boundaryT
156 (
157 cloud.boundaryT().boundaryField()
158 );
159
160 const volVectorField::Boundary& boundaryU
161 (
162 cloud.boundaryU().boundaryField()
163 );
164
165
166 forAll(patches_, p)
167 {
168 label patchi = patches_[p];
169
170 const polyPatch& patch = mesh.boundaryMesh()[patchi];
171
172 // Add mass to the accumulators. negative face area dotted with the
173 // velocity to point flux into the domain.
174
175 // Take a reference to the particleFluxAccumulator for this patch
176 List<Field<scalar>>& pFA = particleFluxAccumulators_[p];
177
178 forAll(pFA, i)
179 {
180 label typeId = moleculeTypeIds_[i];
181
182 scalar mass = cloud.constProps(typeId).mass();
183
184 if (min(boundaryT[patchi]) < SMALL)
185 {
187 << "Zero boundary temperature detected, check boundaryT "
188 << "condition." << nl
189 << nl << abort(FatalError);
190 }
191
192 scalarField mostProbableSpeed
193 (
194 cloud.maxwellianMostProbableSpeed
195 (
196 boundaryT[patchi],
197 mass
198 )
199 );
200
201 // Dotting boundary velocity with the face unit normal
202 // (which points out of the domain, so it must be
203 // negated), dividing by the most probable speed to form
204 // molecularSpeedRatio * cosTheta
205
206 scalarField sCosTheta
207 (
208 (boundaryU[patchi] & -patch.faceAreas()/mag(patch.faceAreas()))
209 / mostProbableSpeed
210 );
211
212 // From Bird eqn 4.22
213
214 pFA[i] +=
215 mag(patch.faceAreas())*numberDensities_[i]*deltaT
216 *mostProbableSpeed
217 *(
218 exp(-sqr(sCosTheta)) + sqrtPi*sCosTheta*(1 + erf(sCosTheta))
219 )
220 /(2.0*sqrtPi);
221 }
222
223 forAll(patch, pFI)
224 {
225 // Loop over all faces as the outer loop to avoid calculating
226 // geometrical properties multiple times.
227
228 const face& f = patch[pFI];
229
230 label globalFaceIndex = pFI + patch.start();
231
232 label celli = mesh.faceOwner()[globalFaceIndex];
233
234 const vector& fC = patch.faceCentres()[pFI];
235
236 scalar fA = mag(patch.faceAreas()[pFI]);
237
239 (
240 mesh,
241 globalFaceIndex,
242 celli
243 );
244
245 // Cumulative triangle area fractions
246 List<scalar> cTriAFracs(faceTets.size(), Zero);
247
248 scalar previousCummulativeSum = 0.0;
249
250 forAll(faceTets, triI)
251 {
252 const tetIndices& faceTetIs = faceTets[triI];
253
254 cTriAFracs[triI] =
255 faceTetIs.faceTri(mesh).mag()/fA
256 + previousCummulativeSum;
257
258 previousCummulativeSum = cTriAFracs[triI];
259 }
260
261 // Force the last area fraction value to 1.0 to avoid any
262 // rounding/non-flat face errors giving a value < 1.0
263 cTriAFracs.last() = 1.0;
264
265 // Normal unit vector *negative* so normal is pointing into the
266 // domain
267 const vector n = -normalised(patch.faceAreas()[pFI]);
268
269 // Wall tangential unit vector. Use the direction between the
270 // face centre and the first vertex in the list
271 const vector t1 = normalised(fC - (mesh.points()[f[0]]));
272
273 // Other tangential unit vector. Rescaling in case face is not
274 // flat and n and t1 aren't perfectly orthogonal
275 const vector t2 = normalised(n ^ t1);
276
277 scalar faceTemperature = boundaryT[patchi][pFI];
278
279 const vector& faceVelocity = boundaryU[patchi][pFI];
280
281 forAll(pFA, i)
282 {
283 scalar& faceAccumulator = pFA[i][pFI];
284
285 // Number of whole particles to insert
286 label nI = max(label(faceAccumulator), 0);
287
288 // Add another particle with a probability proportional to the
289 // remainder of taking the integer part of faceAccumulator
290 if ((faceAccumulator - nI) > rndGen.sample01<scalar>())
291 {
292 nI++;
293 }
294
295 faceAccumulator -= nI;
296
297 label typeId = moleculeTypeIds_[i];
298
299 scalar mass = cloud.constProps(typeId).mass();
300
301 for (label i = 0; i < nI; i++)
302 {
303 // Choose a triangle to insert on, based on their relative
304 // area
305
306 scalar triSelection = rndGen.sample01<scalar>();
307
308 // Selected triangle
309 label selectedTriI = -1;
310
311 forAll(cTriAFracs, triI)
312 {
313 selectedTriI = triI;
314
315 if (cTriAFracs[triI] >= triSelection)
316 {
317 break;
318 }
319 }
320
321 // Randomly distribute the points on the triangle.
322
323 const tetIndices& faceTetIs = faceTets[selectedTriI];
324
325 point p = faceTetIs.faceTri(mesh).randomPoint(rndGen);
326
327 // Velocity generation
328
329 scalar mostProbableSpeed
330 (
331 cloud.maxwellianMostProbableSpeed
332 (
333 faceTemperature,
334 mass
335 )
336 );
337
338 scalar sCosTheta = (faceVelocity & n)/mostProbableSpeed;
339
340 // Coefficients required for Bird eqn 12.5
341 scalar uNormProbCoeffA =
342 sCosTheta + sqrt(sqr(sCosTheta) + 2.0);
343
344 scalar uNormProbCoeffB =
345 0.5*
346 (
347 1.0
348 + sCosTheta*(sCosTheta - sqrt(sqr(sCosTheta) + 2.0))
349 );
350
351 // Equivalent to the QA value in Bird's DSMC3.FOR
352 scalar randomScaling = 3.0;
353
354 if (sCosTheta < -3)
355 {
356 randomScaling = mag(sCosTheta) + 1;
357 }
358
359 scalar P = -1;
360
361 // Normalised candidates for the normal direction velocity
362 // component
363 scalar uNormal;
364 scalar uNormalThermal;
365
366 // Select a velocity using Bird eqn 12.5
367 do
368 {
369 uNormalThermal =
370 randomScaling*(2.0*rndGen.sample01<scalar>() - 1);
371
372 uNormal = uNormalThermal + sCosTheta;
373
374 if (uNormal < 0.0)
375 {
376 P = -1;
377 }
378 else
379 {
380 P = 2.0*uNormal/uNormProbCoeffA
381 *exp(uNormProbCoeffB - sqr(uNormalThermal));
382 }
383
384 } while (P < rndGen.sample01<scalar>());
385
386 vector U =
387 sqrt(physicoChemical::k.value()*faceTemperature/mass)
388 *(
389 rndGen.GaussNormal<scalar>()*t1
390 + rndGen.GaussNormal<scalar>()*t2
391 )
392 + (t1 & faceVelocity)*t1
393 + (t2 & faceVelocity)*t2
394 + mostProbableSpeed*uNormal*n;
395
396 scalar Ei = cloud.equipartitionInternalEnergy
397 (
398 faceTemperature,
399 cloud.constProps(typeId).internalDegreesOfFreedom()
400 );
401
402 cloud.addNewParcel(p, celli, U, Ei, typeId);
403
404 particlesInserted++;
405 }
406 }
407 }
408 }
409
410 reduce(particlesInserted, sumOp<label>());
411
412 Info<< " Particles inserted = "
413 << particlesInserted << endl;
414
415}
416
417
418// ************************************************************************* //
label n
Templated base class for dsmc cloud.
Definition: DSMCCloud.H:75
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects.
Definition: DynamicList.H:72
Generic templated field type.
Definition: Field.H:82
Inserting new particles across the faces of a all patched of type "patch" for a free stream....
Definition: FreeStream.H:54
virtual void inflow()
Introduce particles.
Definition: FreeStream.C:141
virtual void autoMap(const mapPolyMesh &)
Remap the particles to the correct cells following mesh change.
Definition: FreeStream.C:120
virtual ~FreeStream()
Destructor.
Definition: FreeStream.C:113
iterator find(const Key &key)
Find and return an iterator set at the hashed entry.
Definition: HashTableI.H:114
Templated inflow boundary model class.
const dictionary & coeffDict() const
Return the coefficients dictionary.
A 1D array of objects of type <T>, where the size of the vector is known and used for subscript bound...
Definition: List.H:77
void transfer(List< T > &list)
Definition: List.C:447
void setSize(const label n)
Alias for resize()
Definition: List.H:218
void append(T *ptr)
Append an element to the end of the list.
Definition: PtrListI.H:113
Random number generator.
Definition: Random.H:60
Type GaussNormal()
Type sample01()
Return a sample whose components lie in the range [0,1].
scalar deltaTValue() const noexcept
Return time step value.
Definition: TimeStateI.H:43
void size(const label n)
Older name for setAddressableSize.
Definition: UList.H:114
T & last()
Return the last element of the list.
Definition: UListI.H:216
A cloud is a registry collection of lagrangian particles.
Definition: cloud.H:60
A list of keyword definitions, which are a keyword followed by a number of values (eg,...
Definition: dictionary.H:126
T get(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
wordList toc() const
Return the table of contents.
Definition: dictionary.C:602
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
const Time & time() const
Return the top-level database.
Definition: fvMesh.H:290
Class containing mesh-to-mesh mapping information after a change in polyMesh topology.
Definition: mapPolyMesh.H:162
static List< tetIndices > faceTetIndices(const polyMesh &mesh, label fI, label cI)
Return the tet decomposition of the given face, with.
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:81
virtual const labelList & faceOwner() const
Return face owner.
Definition: polyMesh.C:1121
const polyBoundaryMesh & boundaryMesh() const
Return boundary mesh.
Definition: polyMesh.H:456
virtual const pointField & points() const
Return raw points.
Definition: polyMesh.C:1083
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:75
Storage and named access for the indices of a tet which is part of the decomposition of a cell.
Definition: tetIndices.H:84
triPointRef faceTri(const polyMesh &mesh) const
Definition: tetIndicesI.H:149
scalar mag() const
Return scalar magnitude.
Definition: triangleI.H:128
Point randomPoint(Random &rndGen) const
Return a random point on the triangle from a uniform.
Definition: triangleI.H:254
U
Definition: pEqn.H:72
volScalarField & p
const polyBoundaryMesh & patches
dynamicFvMesh & mesh
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
Mathematical constants.
constexpr scalar pi(M_PI)
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
dimensionedScalar exp(const dimensionedScalar &ds)
dimensionedSymmTensor sqr(const dimensionedVector &dv)
quaternion normalised(const quaternion &q)
Return the normalised (unit) quaternion of the given quaternion.
Definition: quaternionI.H:680
dimensionedScalar erf(const dimensionedScalar &ds)
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)
void reduce(const List< UPstream::commsStruct > &comms, T &value, const BinaryOp &bop, const int tag, const label comm)
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
errorManip< error > abort(error &err)
Definition: errorManip.H:144
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
error FatalError
constexpr char nl
The newline '\n' character (0x0a)
Definition: Ostream.H:53
labelList f(nPoints)
dictionary dict
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333
Random rndGen
Definition: createFields.H:23