FreeFOAM The Cross-Platform CFD Toolkit
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) 2009-2010 OpenCFD Ltd.
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 
28 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
29 
30 template <class CloudType>
32 (
33  const dictionary& dict,
34  CloudType& cloud
35 )
36 :
37  InflowBoundaryModel<CloudType>(dict, cloud, typeName),
38  patches_(),
39  moleculeTypeIds_(),
40  numberDensities_(),
41  particleFluxAccumulators_()
42 {
43  // Identify which patches to use
44 
46 
47  forAll(cloud.mesh().boundaryMesh(), p)
48  {
49  const polyPatch& patch = cloud.mesh().boundaryMesh()[p];
50 
51  if (isType<polyPatch>(patch))
52  {
53  patches.append(p);
54  }
55  }
56 
57  patches_.transfer(patches);
58 
59  const dictionary& numberDensitiesDict
60  (
61  this->coeffDict().subDict("numberDensities")
62  );
63 
64  List<word> molecules(numberDensitiesDict.toc());
65 
66  // Initialise the particleFluxAccumulators_
67  particleFluxAccumulators_.setSize(patches_.size());
68 
69  forAll(patches_, p)
70  {
71  const polyPatch& patch = cloud.mesh().boundaryMesh()[patches_[p]];
72 
73  particleFluxAccumulators_[p] = List<Field<scalar> >
74  (
75  molecules.size(),
76  Field<scalar>(patch.size(), 0.0)
77  );
78  }
79 
80  moleculeTypeIds_.setSize(molecules.size());
81 
82  numberDensities_.setSize(molecules.size());
83 
84  forAll(molecules, i)
85  {
86  numberDensities_[i] = readScalar
87  (
88  numberDensitiesDict.lookup(molecules[i])
89  );
90 
91  moleculeTypeIds_[i] = findIndex(cloud.typeIdList(), molecules[i]);
92 
93  if (moleculeTypeIds_[i] == -1)
94  {
96  (
97  "Foam::FreeStream<CloudType>::FreeStream"
98  "("
99  "const dictionary&, "
100  "CloudType&"
101  ")"
102  ) << "typeId " << molecules[i] << "not defined in cloud." << nl
103  << abort(FatalError);
104  }
105  }
106 
107  numberDensities_ /= cloud.nParticle();
108 }
109 
110 
111 
112 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
113 
114 template <class CloudType>
116 {}
117 
118 
119 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
120 
121 template <class CloudType>
123 {
124  CloudType& cloud(this->owner());
125 
126  const polyMesh& mesh(cloud.mesh());
127 
128  const scalar deltaT = mesh.time().deltaTValue();
129 
130  Random& rndGen(cloud.rndGen());
131 
132  scalar sqrtPi = sqrt(mathematicalConstant::pi);
133 
134  label particlesInserted = 0;
135 
136  const volScalarField::GeometricBoundaryField& boundaryT
137  (
138  cloud.boundaryT().boundaryField()
139  );
140 
141  const volVectorField::GeometricBoundaryField& boundaryU
142  (
143  cloud.boundaryU().boundaryField()
144  );
145 
146  forAll(patches_, p)
147  {
148  label patchI = patches_[p];
149 
150  const polyPatch& patch = mesh.boundaryMesh()[patchI];
151 
152  // Add mass to the accumulators. negative face area dotted with the
153  // velocity to point flux into the domain.
154 
155  // Take a reference to the particleFluxAccumulator for this patch
156  List<Field<scalar> >& pFA = particleFluxAccumulators_[p];
157 
158  forAll(pFA, i)
159  {
160  label typeId = moleculeTypeIds_[i];
161 
162  scalar mass = cloud.constProps(typeId).mass();
163 
164  if (min(boundaryT[patchI]) < SMALL)
165  {
166  FatalErrorIn ("Foam::FreeStream<CloudType>::inflow()")
167  << "Zero boundary temperature detected, "
168  << "check boundaryT condition." << nl
169  << nl << abort(FatalError);
170  }
171 
172  scalarField mostProbableSpeed
173  (
174  cloud.maxwellianMostProbableSpeed
175  (
176  boundaryT[patchI],
177  mass
178  )
179  );
180 
181  // Dotting boundary velocity with the face unit normal (which points
182  // out of the domain, so it must be negated), dividing by the most
183  // probable speed to form molecularSpeedRatio * cosTheta
184 
185  scalarField sCosTheta =
186  boundaryU[patchI]
187  & -patch.faceAreas()/mag(patch.faceAreas())
188  /mostProbableSpeed;
189 
190  // From Bird eqn 4.22
191 
192  pFA[i] +=
193  mag(patch.faceAreas())*numberDensities_[i]*deltaT
194  *mostProbableSpeed
195  *(
196  exp(-sqr(sCosTheta)) + sqrtPi*sCosTheta*(1 + erf(sCosTheta))
197  )
198  /(2.0*sqrtPi);
199  }
200 
201  forAll(patch, f)
202  {
203  // Loop over all faces as the outer loop to avoid calculating
204  // geometrical properties multiple times.
205 
206  labelList faceVertices = patch[f];
207 
208  label nVertices = faceVertices.size();
209 
210  label globalFaceIndex = f + patch.start();
211 
212  label cell = mesh.faceOwner()[globalFaceIndex];
213 
214  const vector& fC = patch.faceCentres()[f];
215 
216  scalar fA = mag(patch.faceAreas()[f]);
217 
218  // Cumulative triangle area fractions
219  List<scalar> cTriAFracs(nVertices);
220  cTriAFracs[0] = 0.0;
221 
222  for (label v = 0; v < nVertices - 1; v++)
223  {
224  const point& vA = mesh.points()[faceVertices[v]];
225 
226  const point& vB = mesh.points()[faceVertices[(v + 1)]];
227 
228  cTriAFracs[v] =
229  0.5*mag((vA - fC)^(vB - fC))/fA
230  + cTriAFracs[max((v - 1), 0)];
231  }
232 
233  // Force the last area fraction value to 1.0 to avoid any
234  // rounding/non-flat face errors giving a value < 1.0
235  cTriAFracs[nVertices - 1] = 1.0;
236 
237  // Normal unit vector *negative* so normal is pointing into the
238  // domain
239  vector n = patch.faceAreas()[f];
240  n /= -mag(n);
241 
242  // Wall tangential unit vector. Use the direction between the
243  // face centre and the first vertex in the list
244  vector t1 = fC - (mesh.points()[faceVertices[0]]);
245  t1 /= mag(t1);
246 
247  // Other tangential unit vector. Rescaling in case face is not
248  // flat and n and t1 aren't perfectly orthogonal
249  vector t2 = n^t1;
250  t2 /= mag(t2);
251 
252  scalar faceTemperature = boundaryT[patchI][f];
253 
254  const vector& faceVelocity = boundaryU[patchI][f];
255 
256  forAll(pFA, i)
257  {
258  scalar& faceAccumulator = pFA[i][f];
259 
260  // Number of whole particles to insert
261  label nI = max(label(faceAccumulator), 0);
262 
263  // Add another particle with a probability proportional to the
264  // remainder of taking the integer part of faceAccumulator
265  if ((faceAccumulator - nI) > rndGen.scalar01())
266  {
267  nI++;
268  }
269 
270  faceAccumulator -= nI;
271 
272  label typeId = moleculeTypeIds_[i];
273 
274  scalar mass = cloud.constProps(typeId).mass();
275 
276  for (label i = 0; i < nI; i++)
277  {
278  // Choose a triangle to insert on, based on their relative
279  // area
280 
281  scalar triSelection = rndGen.scalar01();
282 
283  // Selected triangle
284  label sTri = -1;
285 
286  forAll(cTriAFracs, tri)
287  {
288  sTri = tri;
289 
290  if (cTriAFracs[tri] >= triSelection)
291  {
292  break;
293  }
294  }
295 
296  // Randomly distribute the points on the triangle, using the
297  // method from:
298  // Generating Random Points in Triangles
299  // by Greg Turk
300  // from "Graphics Gems", Academic Press, 1990
301  // http://tog.acm.org/GraphicsGems/gems/TriPoints.c
302 
303  const point& A = fC;
304  const point& B = mesh.points()[faceVertices[sTri]];
305  const point& C =
306  mesh.points()[faceVertices[(sTri + 1) % nVertices]];
307 
308  scalar s = rndGen.scalar01();
309  scalar t = sqrt(rndGen.scalar01());
310 
311  point p = (1 - t)*A + (1 - s)*t*B + s*t*C;
312 
313  // Velocity generation
314 
315  scalar mostProbableSpeed
316  (
317  cloud.maxwellianMostProbableSpeed
318  (
319  faceTemperature,
320  mass
321  )
322  );
323 
324  scalar sCosTheta = (faceVelocity & n)/mostProbableSpeed;
325 
326  // Coefficients required for Bird eqn 12.5
327  scalar uNormProbCoeffA =
328  sCosTheta + sqrt(sqr(sCosTheta) + 2.0);
329 
330  scalar uNormProbCoeffB =
331  0.5*
332  (
333  1.0
334  + sCosTheta*(sCosTheta - sqrt(sqr(sCosTheta) + 2.0))
335  );
336 
337  // Equivalent to the QA value in Bird's DSMC3.FOR
338  scalar randomScaling = 3.0;
339 
340  if (sCosTheta < -3)
341  {
342  randomScaling = mag(sCosTheta) + 1;
343  }
344 
345  scalar P = -1;
346 
347  // Normalised candidates for the normal direction velocity
348  // component
349  scalar uNormal;
350  scalar uNormalThermal;
351 
352  // Select a velocity using Bird eqn 12.5
353  do
354  {
355  uNormalThermal =
356  randomScaling*(2.0*rndGen.scalar01() - 1);
357 
358  uNormal = uNormalThermal + sCosTheta;
359 
360  if (uNormal < 0.0)
361  {
362  P = -1;
363  }
364  else
365  {
366  P = 2.0*uNormal/uNormProbCoeffA
367  *exp(uNormProbCoeffB - sqr(uNormalThermal));
368  }
369 
370  } while (P < rndGen.scalar01());
371 
372  vector U =
373  sqrt(CloudType::kb*faceTemperature/mass)
374  *(
375  rndGen.GaussNormal()*t1
376  + rndGen.GaussNormal()*t2
377  )
378  + (t1 & faceVelocity)*t1
379  + (t2 & faceVelocity)*t2
380  + mostProbableSpeed*uNormal*n;
381 
382  scalar Ei = cloud.equipartitionInternalEnergy
383  (
384  faceTemperature,
385  cloud.constProps(typeId).internalDegreesOfFreedom()
386  );
387 
388  cloud.addNewParcel
389  (
390  p,
391  U,
392  Ei,
393  cell,
394  typeId
395  );
396 
397  particlesInserted++;
398  }
399  }
400  }
401  }
402 
403  reduce(particlesInserted, sumOp<label>());
404 
405  Info<< " Particles inserted = "
406  << particlesInserted << endl;
407 
408 }
409 
410 
411 // ************************ vim: set sw=4 sts=4 et: ************************ //