FreeStream.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2011-2022 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "FreeStream.H"
27 #include "constants.H"
28 #include "triPointRef.H"
29 #include "tetIndices.H"
30 
31 using namespace Foam::constant::mathematical;
32 
33 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
34 
35 template<class CloudType>
37 (
38  const dictionary& dict,
39  CloudType& cloud
40 )
41 :
42  InflowBoundaryModel<CloudType>(dict, cloud, typeName),
43  patches_(),
44  moleculeTypeIds_(),
45  numberDensities_(),
46  particleFluxAccumulators_()
47 {
48  // Identify which patches to use
49 
51 
52  forAll(cloud.mesh().boundaryMesh(), p)
53  {
54  const polyPatch& patch = cloud.mesh().boundaryMesh()[p];
55 
56  if (isType<polyPatch>(patch))
57  {
58  patches.append(p);
59  }
60  }
61 
62  patches_.transfer(patches);
63 
64  const dictionary& numberDensitiesDict
65  (
66  this->coeffDict().subDict("numberDensities")
67  );
68 
69  List<word> molecules(numberDensitiesDict.toc());
70 
71  // Initialise the particleFluxAccumulators_
72  particleFluxAccumulators_.setSize(patches_.size());
73 
74  forAll(patches_, p)
75  {
76  const polyPatch& patch = cloud.mesh().boundaryMesh()[patches_[p]];
77 
78  particleFluxAccumulators_[p] = List<Field<scalar>>
79  (
80  molecules.size(),
81  Field<scalar>(patch.size(), 0.0)
82  );
83  }
84 
85  moleculeTypeIds_.setSize(molecules.size());
86 
87  numberDensities_.setSize(molecules.size());
88 
89  forAll(molecules, i)
90  {
91  numberDensities_[i] =
92  numberDensitiesDict.lookup<scalar>(molecules[i]);
93 
94  moleculeTypeIds_[i] = findIndex(cloud.typeIdList(), molecules[i]);
95 
96  if (moleculeTypeIds_[i] == -1)
97  {
99  << "typeId " << molecules[i] << "not defined in cloud." << nl
100  << abort(FatalError);
101  }
102  }
103 
104  numberDensities_ /= cloud.nParticle();
105 }
106 
107 
108 
109 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
110 
111 template<class CloudType>
113 {}
114 
115 
116 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
117 
118 template<class CloudType>
120 {
121  CloudType& cloud(this->owner());
122  const polyMesh& mesh(cloud.mesh());
123 
124  forAll(patches_, p)
125  {
126  label patchi = patches_[p];
127 
128  const polyPatch& patch = mesh.boundaryMesh()[patchi];
129  List<Field<scalar>>& pFA = particleFluxAccumulators_[p];
130 
131  forAll(pFA, facei)
132  {
133  pFA[facei].setSize(patch.size(), 0);
134  }
135  }
136 }
137 
138 
139 template<class CloudType>
141 {
142  CloudType& cloud(this->owner());
143 
144  const polyMesh& mesh(cloud.mesh());
145 
146  const scalar deltaT = mesh.time().deltaTValue();
147 
148  Random& rndGen(cloud.rndGen());
149 
150  scalar sqrtPi = sqrt(pi);
151 
152  label particlesInserted = 0;
153 
154  const volScalarField::Boundary& boundaryT
155  (
156  cloud.boundaryT().boundaryField()
157  );
158 
159  const volVectorField::Boundary& boundaryU
160  (
161  cloud.boundaryU().boundaryField()
162  );
163 
164 
165  forAll(patches_, p)
166  {
167  label patchi = patches_[p];
168 
169  const polyPatch& patch = mesh.boundaryMesh()[patchi];
170 
171  // Add mass to the accumulators. negative face area dotted with the
172  // velocity to point flux into the domain.
173 
174  // Take a reference to the particleFluxAccumulator for this patch
175  List<Field<scalar>>& pFA = particleFluxAccumulators_[p];
176 
177  forAll(pFA, i)
178  {
179  label typeId = moleculeTypeIds_[i];
180 
181  scalar mass = cloud.constProps(typeId).mass();
182 
183  if (min(boundaryT[patchi]) < small)
184  {
186  << "Zero boundary temperature detected, check boundaryT "
187  << "condition." << nl
188  << nl << abort(FatalError);
189  }
190 
191  scalarField mostProbableSpeed
192  (
194  (
195  boundaryT[patchi],
196  mass
197  )
198  );
199 
200  // Dotting boundary velocity with the face unit normal
201  // (which points out of the domain, so it must be
202  // negated), dividing by the most probable speed to form
203  // molecularSpeedRatio * cosTheta
204 
205  scalarField sCosTheta
206  (
207  (boundaryU[patchi] & -patch.faceAreas()/patch.magFaceAreas())
208  / mostProbableSpeed
209  );
210 
211  // From Bird eqn 4.22
212 
213  pFA[i] +=
214  patch.magFaceAreas()*numberDensities_[i]*deltaT
215  *mostProbableSpeed
216  *(
217  exp(-sqr(sCosTheta)) + sqrtPi*sCosTheta*(1 + erf(sCosTheta))
218  )
219  /(2.0*sqrtPi);
220  }
221 
222  forAll(patch, pFI)
223  {
224  // Loop over all faces as the outer loop to avoid calculating
225  // geometrical properties multiple times.
226 
227  const face& f = patch[pFI];
228 
229  label globalFaceIndex = pFI + patch.start();
230 
231  label celli = mesh.faceOwner()[globalFaceIndex];
232 
233  const vector& fC = patch.faceCentres()[pFI];
234 
235  scalar fA = patch.magFaceAreas()[pFI];
236 
238  (
239  mesh,
240  globalFaceIndex,
241  celli
242  );
243 
244  // Cumulative triangle area fractions
245  List<scalar> cTriAFracs(faceTets.size(), 0.0);
246 
247  scalar previousCumulativeSum = 0.0;
248 
249  forAll(faceTets, triI)
250  {
251  const tetIndices& faceTetIs = faceTets[triI];
252 
253  cTriAFracs[triI] =
254  faceTetIs.faceTri(mesh).mag()/fA
255  + previousCumulativeSum;
256 
257  previousCumulativeSum = cTriAFracs[triI];
258  }
259 
260  // Force the last area fraction value to 1.0 to avoid any
261  // rounding/non-flat face errors giving a value < 1.0
262  cTriAFracs.last() = 1.0;
263 
264  // Normal unit vector *negative* so normal is pointing into the
265  // domain
266  vector n = patch.faceAreas()[pFI];
267  n /= -mag(n);
268 
269  // Wall tangential unit vector. Use the direction between the
270  // face centre and the first vertex in the list
271  vector t1 = fC - (mesh.points()[f[0]]);
272  t1 /= mag(t1);
273 
274  // Other tangential unit vector. Rescaling in case face is not
275  // flat and n and t1 aren't perfectly orthogonal
276  vector t2 = n^t1;
277  t2 /= mag(t2);
278 
279  scalar faceTemperature = boundaryT[patchi][pFI];
280 
281  const vector& faceVelocity = boundaryU[patchi][pFI];
282 
283  forAll(pFA, i)
284  {
285  scalar& faceAccumulator = pFA[i][pFI];
286 
287  // Number of whole particles to insert
288  label nI = max(label(faceAccumulator), 0);
289 
290  // Add another particle with a probability proportional to the
291  // remainder of taking the integer part of faceAccumulator
292  if ((faceAccumulator - nI) > rndGen.scalar01())
293  {
294  nI++;
295  }
296 
297  faceAccumulator -= nI;
298 
299  label typeId = moleculeTypeIds_[i];
300 
301  scalar mass = cloud.constProps(typeId).mass();
302 
303  for (label i = 0; i < nI; i++)
304  {
305  // Choose a triangle to insert on, based on their relative
306  // area
307 
308  scalar triSelection = rndGen.scalar01();
309 
310  // Selected triangle
311  label selectedTriI = -1;
312 
313  forAll(cTriAFracs, triI)
314  {
315  selectedTriI = triI;
316 
317  if (cTriAFracs[triI] >= triSelection)
318  {
319  break;
320  }
321  }
322 
323  // Randomly distribute the points on the triangle.
324 
325  const tetIndices& faceTetIs = faceTets[selectedTriI];
326 
327  point p = faceTetIs.faceTri(mesh).randomPoint(rndGen);
328 
329  // Velocity generation
330 
331  scalar mostProbableSpeed
332  (
334  (
335  faceTemperature,
336  mass
337  )
338  );
339 
340  scalar sCosTheta = (faceVelocity & n)/mostProbableSpeed;
341 
342  // Coefficients required for Bird eqn 12.5
343  scalar uNormProbCoeffA =
344  sCosTheta + sqrt(sqr(sCosTheta) + 2.0);
345 
346  scalar uNormProbCoeffB =
347  0.5*
348  (
349  1.0
350  + sCosTheta*(sCosTheta - sqrt(sqr(sCosTheta) + 2.0))
351  );
352 
353  // Equivalent to the QA value in Bird's DSMC3.FOR
354  scalar randomScaling = 3.0;
355 
356  if (sCosTheta < -3)
357  {
358  randomScaling = mag(sCosTheta) + 1;
359  }
360 
361  scalar P = -1;
362 
363  // Normalised candidates for the normal direction velocity
364  // component
365  scalar uNormal;
366  scalar uNormalThermal;
367 
368  // Select a velocity using Bird eqn 12.5
369  do
370  {
371  uNormalThermal =
372  randomScaling*(2.0*rndGen.scalar01() - 1);
373 
374  uNormal = uNormalThermal + sCosTheta;
375 
376  if (uNormal < 0.0)
377  {
378  P = -1;
379  }
380  else
381  {
382  P = 2.0*uNormal/uNormProbCoeffA
383  *exp(uNormProbCoeffB - sqr(uNormalThermal));
384  }
385 
386  } while (P < rndGen.scalar01());
387 
388  vector U =
389  sqrt(physicoChemical::k.value()*faceTemperature/mass)
390  *(
391  rndGen.scalarNormal()*t1
392  + rndGen.scalarNormal()*t2
393  )
394  + (t1 & faceVelocity)*t1
395  + (t2 & faceVelocity)*t2
396  + mostProbableSpeed*uNormal*n;
397 
398  scalar Ei = cloud.equipartitionInternalEnergy
399  (
400  faceTemperature,
401  cloud.constProps(typeId).internalDegreesOfFreedom()
402  );
403 
404  cloud.addNewParcel(p, celli, U, Ei, typeId);
405 
406  particlesInserted++;
407  }
408  }
409  }
410  }
411 
412  reduce(particlesInserted, sumOp<label>());
413 
414  Info<< " Particles inserted = "
415  << particlesInserted << endl;
416 
417 }
418 
419 
420 // ************************************************************************* //
const fvPatchList & patches
const volScalarField & boundaryT() const
Return macroscopic temperature.
Definition: DSMCCloudI.H:193
dictionary dict
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
layerAndWeight max(const layerAndWeight &a, const layerAndWeight &b)
FvWallInfoData< WallInfo, label > label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
virtual void autoMap(const polyTopoChangeMap &)
Remap the particles to the correct cells following mesh change.
Definition: FreeStream.C:119
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
U
Definition: pEqn.H:72
scalar mag() const
Return scalar magnitude.
Definition: triangleI.H:103
error FatalError
A list of keyword definitions, which are a keyword followed by any number of values (e...
Definition: dictionary.H:156
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:306
const Boundary & boundaryField() const
Return const-reference to the boundary field.
triPointRef faceTri(const polyMesh &mesh) const
Return the geometry corresponding to the tri on the face for.
Definition: tetIndicesI.H:123
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Templated inflow boundary model class.
Definition: DSMCCloud.H:61
void size(const label)
Override size to be inconsistent with allocated storage.
Definition: ListI.H:164
dimensionedScalar sqrt(const dimensionedScalar &ds)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
fvMesh & mesh
wordList toc() const
Return the table of contents.
Definition: dictionary.C:1115
Random rndGen(label(0))
virtual void inflow()
Introduce particles.
Definition: FreeStream.C:140
FreeStream(const dictionary &dict, CloudType &cloud)
Construct from dictionary.
Definition: FreeStream.C:37
dimensionedScalar exp(const dimensionedScalar &ds)
A 1D vector of objects of type <T> that resizes itself as necessary to accept the new objects...
Definition: DynamicList.H:56
const volVectorField & boundaryU() const
Return macroscopic velocity.
Definition: DSMCCloudI.H:201
mathematical constants.
Random & rndGen()
Return references to the random object.
Definition: DSMCCloudI.H:121
Class containing mesh-to-mesh mapping information after a change in polyMesh topology.
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
Definition: DynamicListI.H:296
layerAndWeight min(const layerAndWeight &a, const layerAndWeight &b)
Storage and named access for the indices of a tet which is part of the decomposition of a cell...
Definition: tetIndices.H:81
errorManip< error > abort(error &err)
Definition: errorManip.H:131
scalar maxwellianMostProbableSpeed(scalar temperature, scalar mass) const
Most probable speed.
Definition: DSMCCloudI.H:388
Random number generator.
Definition: Random.H:57
scalar equipartitionInternalEnergy(scalar temperature, direction internalDegreesOfFreedom)
Generate a random internal energy, sampled from the.
Definition: DSMCCloud.C:1027
const vectorField::subField faceAreas() const
Return face areas.
Definition: polyPatch.C:309
static const char nl
Definition: Ostream.H:260
dimensionedScalar erf(const dimensionedScalar &ds)
labelList f(nPoints)
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurrence of given element and return index,.
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
void addNewParcel(const vector &position, const label celli, const vector &U, const scalar Ei, const label typeId)
Add new parcel.
Definition: DSMCCloud.C:454
const fvMesh & mesh() const
Return references to the mesh.
Definition: DSMCCloudI.H:41
void setSize(const label)
Reset size of List.
Definition: List.C:281
label patchi
const scalarField::subField magFaceAreas() const
Return face area magnitudes.
Definition: polyPatch.C:315
const dimensionedScalar k
Boltzmann constant.
label start() const
Return start label of this patch in the polyMesh face list.
Definition: polyPatch.H:306
messageStream Info
Point randomPoint(Random &rndGen) const
Return a random point on the triangle from a uniform.
Definition: triangleI.H:241
dimensioned< scalar > mag(const dimensioned< Type > &)
label n
Mesh consisting of general polyhedral cells.
Definition: polyMesh.H:76
volScalarField & p
A patch is a list of labels that address the faces in the global face list.
Definition: polyPatch.H:66
T & last()
Return the last element of the list.
Definition: UListI.H:128
virtual ~FreeStream()
Destructor.
Definition: FreeStream.C:112
Templated base class for dsmc cloud.
Definition: DSMCCloud.H:75
const List< typename ParcelType::constantProperties > & constProps() const
Return all of the constant properties.
Definition: DSMCCloudI.H:95
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:864
static List< tetIndices > faceTetIndices(const polyMesh &mesh, label fI, label cI)
Return the tet decomposition of the given face, with.
const vectorField::subField faceCentres() const
Return face centres.
Definition: polyPatch.C:303