30 template <
class CloudType>
41 particleFluxAccumulators_()
47 forAll(cloud.mesh().boundaryMesh(),
p)
49 const polyPatch& patch = cloud.mesh().boundaryMesh()[
p];
51 if (isType<polyPatch>(patch))
57 patches_.transfer(patches);
61 this->coeffDict().subDict(
"numberDensities")
67 particleFluxAccumulators_.
setSize(patches_.size());
71 const polyPatch& patch = cloud.mesh().boundaryMesh()[patches_[
p]];
80 moleculeTypeIds_.
setSize(molecules.size());
82 numberDensities_.setSize(molecules.size());
88 numberDensitiesDict.
lookup(molecules[i])
91 moleculeTypeIds_[i] =
findIndex(cloud.typeIdList(), molecules[i]);
93 if (moleculeTypeIds_[i] == -1)
97 "Foam::FreeStream<CloudType>::FreeStream"
102 ) <<
"typeId " << molecules[i] <<
"not defined in cloud." <<
nl
107 numberDensities_ /= cloud.nParticle();
114 template <
class CloudType>
121 template <
class CloudType>
124 CloudType& cloud(this->owner());
128 const scalar deltaT =
mesh.time().deltaTValue();
130 Random& rndGen(cloud.rndGen());
134 label particlesInserted = 0;
136 const volScalarField::GeometricBoundaryField& boundaryT
138 cloud.boundaryT().boundaryField()
141 const volVectorField::GeometricBoundaryField& boundaryU
143 cloud.boundaryU().boundaryField()
148 label patchI = patches_[
p];
160 label typeId = moleculeTypeIds_[i];
162 scalar mass = cloud.constProps(typeId).mass();
164 if (
min(boundaryT[patchI]) < SMALL)
167 <<
"Zero boundary temperature detected, "
168 <<
"check boundaryT condition." <<
nl
174 cloud.maxwellianMostProbableSpeed
196 exp(-
sqr(sCosTheta)) + sqrtPi*sCosTheta*(1 +
erf(sCosTheta))
208 label nVertices = faceVertices.
size();
210 label globalFaceIndex =
f + patch.
start();
212 label
cell =
mesh.faceOwner()[globalFaceIndex];
222 for (label v = 0; v < nVertices - 1; v++)
224 const point& vA =
mesh.points()[faceVertices[v]];
226 const point& vB =
mesh.points()[faceVertices[(v + 1)]];
229 0.5*
mag((vA - fC)^(vB - fC))/fA
230 + cTriAFracs[
max((v - 1), 0)];
235 cTriAFracs[nVertices - 1] = 1.0;
244 vector t1 = fC - (
mesh.points()[faceVertices[0]]);
252 scalar faceTemperature = boundaryT[patchI][
f];
254 const vector& faceVelocity = boundaryU[patchI][
f];
258 scalar& faceAccumulator = pFA[i][
f];
261 label nI =
max(label(faceAccumulator), 0);
265 if ((faceAccumulator - nI) > rndGen.scalar01())
270 faceAccumulator -= nI;
272 label typeId = moleculeTypeIds_[i];
274 scalar mass = cloud.constProps(typeId).mass();
276 for (label i = 0; i < nI; i++)
281 scalar triSelection = rndGen.scalar01();
290 if (cTriAFracs[tri] >= triSelection)
304 const point& B =
mesh.points()[faceVertices[sTri]];
306 mesh.points()[faceVertices[(sTri + 1) % nVertices]];
308 scalar s = rndGen.scalar01();
309 scalar t =
sqrt(rndGen.scalar01());
311 point p = (1 - t)*A + (1 - s)*t*B + s*t*C;
315 scalar mostProbableSpeed
317 cloud.maxwellianMostProbableSpeed
324 scalar sCosTheta = (faceVelocity & n)/mostProbableSpeed;
327 scalar uNormProbCoeffA =
328 sCosTheta +
sqrt(
sqr(sCosTheta) + 2.0);
330 scalar uNormProbCoeffB =
334 + sCosTheta*(sCosTheta -
sqrt(
sqr(sCosTheta) + 2.0))
338 scalar randomScaling = 3.0;
342 randomScaling =
mag(sCosTheta) + 1;
350 scalar uNormalThermal;
356 randomScaling*(2.0*rndGen.scalar01() - 1);
358 uNormal = uNormalThermal + sCosTheta;
366 P = 2.0*uNormal/uNormProbCoeffA
367 *
exp(uNormProbCoeffB -
sqr(uNormalThermal));
370 }
while (P < rndGen.scalar01());
373 sqrt(CloudType::kb*faceTemperature/mass)
375 rndGen.GaussNormal()*t1
376 + rndGen.GaussNormal()*t2
378 + (t1 & faceVelocity)*t1
379 + (t2 & faceVelocity)*t2
380 + mostProbableSpeed*uNormal*n;
382 scalar Ei = cloud.equipartitionInternalEnergy
385 cloud.constProps(typeId).internalDegreesOfFreedom()
405 Info<<
" Particles inserted = "
406 << particlesInserted <<
endl;