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