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-2018 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] = readScalar
92  (
93  numberDensitiesDict.lookup(molecules[i])
94  );
95 
96  moleculeTypeIds_[i] = findIndex(cloud.typeIdList(), molecules[i]);
97 
98  if (moleculeTypeIds_[i] == -1)
99  {
101  << "typeId " << molecules[i] << "not defined in cloud." << nl
102  << abort(FatalError);
103  }
104  }
105 
106  numberDensities_ /= cloud.nParticle();
107 }
108 
109 
110 
111 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
112 
113 template<class CloudType>
115 {}
116 
117 
118 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
119 
120 template<class CloudType>
122 {
123  CloudType& cloud(this->owner());
124  const polyMesh& mesh(cloud.mesh());
125 
126  forAll(patches_, p)
127  {
128  label patchi = patches_[p];
129 
130  const polyPatch& patch = mesh.boundaryMesh()[patchi];
131  List<Field<scalar>>& pFA = particleFluxAccumulators_[p];
132 
133  forAll(pFA, facei)
134  {
135  pFA[facei].setSize(patch.size(), 0);
136  }
137  }
138 }
139 
140 
141 template<class CloudType>
143 {
144  CloudType& cloud(this->owner());
145 
146  const polyMesh& mesh(cloud.mesh());
147 
148  const scalar deltaT = mesh.time().deltaTValue();
149 
150  Random& rndGen(cloud.rndGen());
151 
152  scalar sqrtPi = sqrt(pi);
153 
154  label particlesInserted = 0;
155 
156  const volScalarField::Boundary& boundaryT
157  (
158  cloud.boundaryT().boundaryField()
159  );
160 
161  const volVectorField::Boundary& boundaryU
162  (
163  cloud.boundaryU().boundaryField()
164  );
165 
166 
167  forAll(patches_, p)
168  {
169  label patchi = patches_[p];
170 
171  const polyPatch& patch = mesh.boundaryMesh()[patchi];
172 
173  // Add mass to the accumulators. negative face area dotted with the
174  // velocity to point flux into the domain.
175 
176  // Take a reference to the particleFluxAccumulator for this patch
177  List<Field<scalar>>& pFA = particleFluxAccumulators_[p];
178 
179  forAll(pFA, i)
180  {
181  label typeId = moleculeTypeIds_[i];
182 
183  scalar mass = cloud.constProps(typeId).mass();
184 
185  if (min(boundaryT[patchi]) < small)
186  {
188  << "Zero boundary temperature detected, check boundaryT "
189  << "condition." << nl
190  << nl << abort(FatalError);
191  }
192 
193  scalarField mostProbableSpeed
194  (
196  (
197  boundaryT[patchi],
198  mass
199  )
200  );
201 
202  // Dotting boundary velocity with the face unit normal
203  // (which points out of the domain, so it must be
204  // negated), dividing by the most probable speed to form
205  // molecularSpeedRatio * cosTheta
206 
207  scalarField sCosTheta
208  (
209  (boundaryU[patchi] & -patch.faceAreas()/mag(patch.faceAreas()))
210  / mostProbableSpeed
211  );
212 
213  // From Bird eqn 4.22
214 
215  pFA[i] +=
216  mag(patch.faceAreas())*numberDensities_[i]*deltaT
217  *mostProbableSpeed
218  *(
219  exp(-sqr(sCosTheta)) + sqrtPi*sCosTheta*(1 + erf(sCosTheta))
220  )
221  /(2.0*sqrtPi);
222  }
223 
224  forAll(patch, pFI)
225  {
226  // Loop over all faces as the outer loop to avoid calculating
227  // geometrical properties multiple times.
228 
229  const face& f = patch[pFI];
230 
231  label globalFaceIndex = pFI + patch.start();
232 
233  label celli = mesh.faceOwner()[globalFaceIndex];
234 
235  const vector& fC = patch.faceCentres()[pFI];
236 
237  scalar fA = mag(patch.faceAreas()[pFI]);
238 
240  (
241  mesh,
242  globalFaceIndex,
243  celli
244  );
245 
246  // Cumulative triangle area fractions
247  List<scalar> cTriAFracs(faceTets.size(), 0.0);
248 
249  scalar previousCummulativeSum = 0.0;
250 
251  forAll(faceTets, triI)
252  {
253  const tetIndices& faceTetIs = faceTets[triI];
254 
255  cTriAFracs[triI] =
256  faceTetIs.faceTri(mesh).mag()/fA
257  + previousCummulativeSum;
258 
259  previousCummulativeSum = cTriAFracs[triI];
260  }
261 
262  // Force the last area fraction value to 1.0 to avoid any
263  // rounding/non-flat face errors giving a value < 1.0
264  cTriAFracs.last() = 1.0;
265 
266  // Normal unit vector *negative* so normal is pointing into the
267  // domain
268  vector n = patch.faceAreas()[pFI];
269  n /= -mag(n);
270 
271  // Wall tangential unit vector. Use the direction between the
272  // face centre and the first vertex in the list
273  vector t1 = fC - (mesh.points()[f[0]]);
274  t1 /= mag(t1);
275 
276  // Other tangential unit vector. Rescaling in case face is not
277  // flat and n and t1 aren't perfectly orthogonal
278  vector t2 = n^t1;
279  t2 /= mag(t2);
280 
281  scalar faceTemperature = boundaryT[patchi][pFI];
282 
283  const vector& faceVelocity = boundaryU[patchi][pFI];
284 
285  forAll(pFA, i)
286  {
287  scalar& faceAccumulator = pFA[i][pFI];
288 
289  // Number of whole particles to insert
290  label nI = max(label(faceAccumulator), 0);
291 
292  // Add another particle with a probability proportional to the
293  // remainder of taking the integer part of faceAccumulator
294  if ((faceAccumulator - nI) > rndGen.scalar01())
295  {
296  nI++;
297  }
298 
299  faceAccumulator -= nI;
300 
301  label typeId = moleculeTypeIds_[i];
302 
303  scalar mass = cloud.constProps(typeId).mass();
304 
305  for (label i = 0; i < nI; i++)
306  {
307  // Choose a triangle to insert on, based on their relative
308  // area
309 
310  scalar triSelection = rndGen.scalar01();
311 
312  // Selected triangle
313  label selectedTriI = -1;
314 
315  forAll(cTriAFracs, triI)
316  {
317  selectedTriI = triI;
318 
319  if (cTriAFracs[triI] >= triSelection)
320  {
321  break;
322  }
323  }
324 
325  // Randomly distribute the points on the triangle.
326 
327  const tetIndices& faceTetIs = faceTets[selectedTriI];
328 
329  point p = faceTetIs.faceTri(mesh).randomPoint(rndGen);
330 
331  // Velocity generation
332 
333  scalar mostProbableSpeed
334  (
336  (
337  faceTemperature,
338  mass
339  )
340  );
341 
342  scalar sCosTheta = (faceVelocity & n)/mostProbableSpeed;
343 
344  // Coefficients required for Bird eqn 12.5
345  scalar uNormProbCoeffA =
346  sCosTheta + sqrt(sqr(sCosTheta) + 2.0);
347 
348  scalar uNormProbCoeffB =
349  0.5*
350  (
351  1.0
352  + sCosTheta*(sCosTheta - sqrt(sqr(sCosTheta) + 2.0))
353  );
354 
355  // Equivalent to the QA value in Bird's DSMC3.FOR
356  scalar randomScaling = 3.0;
357 
358  if (sCosTheta < -3)
359  {
360  randomScaling = mag(sCosTheta) + 1;
361  }
362 
363  scalar P = -1;
364 
365  // Normalised candidates for the normal direction velocity
366  // component
367  scalar uNormal;
368  scalar uNormalThermal;
369 
370  // Select a velocity using Bird eqn 12.5
371  do
372  {
373  uNormalThermal =
374  randomScaling*(2.0*rndGen.scalar01() - 1);
375 
376  uNormal = uNormalThermal + sCosTheta;
377 
378  if (uNormal < 0.0)
379  {
380  P = -1;
381  }
382  else
383  {
384  P = 2.0*uNormal/uNormProbCoeffA
385  *exp(uNormProbCoeffB - sqr(uNormalThermal));
386  }
387 
388  } while (P < rndGen.scalar01());
389 
390  vector U =
391  sqrt(physicoChemical::k.value()*faceTemperature/mass)
392  *(
393  rndGen.scalarNormal()*t1
394  + rndGen.scalarNormal()*t2
395  )
396  + (t1 & faceVelocity)*t1
397  + (t2 & faceVelocity)*t2
398  + mostProbableSpeed*uNormal*n;
399 
400  scalar Ei = cloud.equipartitionInternalEnergy
401  (
402  faceTemperature,
403  cloud.constProps(typeId).internalDegreesOfFreedom()
404  );
405 
406  cloud.addNewParcel(p, celli, U, Ei, typeId);
407 
408  particlesInserted++;
409  }
410  }
411  }
412  }
413 
414  reduce(particlesInserted, sumOp<label>());
415 
416  Info<< " Particles inserted = "
417  << particlesInserted << endl;
418 
419 }
420 
421 
422 // ************************************************************************* //
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
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
A face is a list of labels corresponding to mesh vertices.
Definition: face.H:75
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:158
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
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:113
dimensionedSymmTensor sqr(const dimensionedVector &dv)
Templated inflow boundary model class.
Definition: DSMCCloud.H:62
dimensionedScalar sqrt(const dimensionedScalar &ds)
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:256
virtual void autoMap(const mapPolyMesh &)
Remap the particles to the correct cells following mesh change.
Definition: FreeStream.C:121
patches[0]
Class containing mesh-to-mesh mapping information after a change in polyMesh topology.
Definition: mapPolyMesh.H:158
wordList toc() const
Return the table of contents.
Definition: dictionary.C:783
Random rndGen(label(0))
virtual void inflow()
Introduce particles.
Definition: FreeStream.C:142
FreeStream(const dictionary &dict, CloudType &cloud)
Construct from dictionary.
Definition: FreeStream.C:37
dynamicFvMesh & mesh
dimensionedScalar exp(const dimensionedScalar &ds)
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
DynamicList< T, SizeInc, SizeMult, SizeDiv > & append(const T &)
Append an element at the end of the list.
Definition: DynamicListI.H:296
Storage and named access for the indices of a tet which is part of the decomposition of a cell...
Definition: tetIndices.H:81
bool readScalar(const char *buf, doubleScalar &s)
Read whole of buf as a scalar. Return true if successful.
Definition: doubleScalar.H:68
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:1025
const vectorField::subField faceAreas() const
Return face normals.
Definition: polyPatch.C:290
static const char nl
Definition: Ostream.H:265
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,.
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
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:456
const fvMesh & mesh() const
Return references to the mesh.
Definition: DSMCCloudI.H:41
const dimensionedScalar k
Boltzmann constant.
void setSize(const label)
Reset size of List.
Definition: List.C:281
label patchi
U
Definition: pEqn.H:72
label start() const
Return start label of this patch in the polyMesh face list.
Definition: polyPatch.H:303
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:74
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:114
Templated base class for dsmc cloud.
Definition: DSMCCloud.H:69
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:583
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:284