src/lagrangian/dsmc/submodels/InflowBoundaryModel/FreeStream/FreeStream.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     |
   5     \\  /    A nd           | For copyright notice see file Copyright
   6      \\/     M anipulation  |
   7 -------------------------------------------------------------------------------
   8 License
   9     This file is part of foam-extend.
  10
  11     foam-extend is free software: you can redistribute it and/or modify it
  12     under the terms of the GNU General Public License as published by the
  13     Free Software Foundation, either version 3 of the License, or (at your
  14     option) any later version.
  15
  16     foam-extend is distributed in the hope that it will be useful, but
  17     WITHOUT ANY WARRANTY; without even the implied warranty of
  18     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19     General Public License for more details.
  20
  21     You should have received a copy of the GNU General Public License
  22     along with foam-extend.  If not, see <http://www.gnu.org/licenses/>.
  23
  24 \*---------------------------------------------------------------------------*/
  25
  26 #include "FreeStream.H"
  27
  28 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  29
  30 template <class CloudType>
  31 Foam::FreeStream<CloudType>::FreeStream
  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
  45     DynamicList<label> patches;
  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         {
  95             FatalErrorIn
  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>
 115 Foam::FreeStream<CloudType>::~FreeStream()
 116 {}
 117
 118
 119 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 120
 121 template <class CloudType>
 122 void Foam::FreeStream<CloudType>::inflow()
 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
 221             scalar previousCummulativeSum = 0.0;
 222
 223             for (label v = 0; v < nVertices - 1; v++)
 224             {
 225                 const point& vA = mesh.points()[faceVertices[v]];
 226
 227                 const point& vB = mesh.points()[faceVertices[(v + 1)]];
 228
 229                 cTriAFracs[v] =
 230                     0.5*mag((vA - fC)^(vB - fC))/fA
 231                   + previousCummulativeSum;
 232
 233                 previousCummulativeSum = cTriAFracs[v];
 234             }
 235
 236             // Force the last area fraction value to 1.0 to avoid any
 237             // rounding/non-flat face errors giving a value < 1.0
 238             cTriAFracs[nVertices - 1] = 1.0;
 239
 240             // Normal unit vector *negative* so normal is pointing into the
 241             // domain
 242             vector n = patch.faceAreas()[f];
 243             n /= -mag(n);
 244
 245             // Wall tangential unit vector. Use the direction between the
 246             // face centre and the first vertex in the list
 247             vector t1 = fC - (mesh.points()[faceVertices[0]]);
 248             t1 /= mag(t1);
 249
 250             // Other tangential unit vector.  Rescaling in case face is not
 251             // flat and n and t1 aren't perfectly orthogonal
 252             vector t2 = n^t1;
 253             t2 /= mag(t2);
 254
 255             scalar faceTemperature = boundaryT[patchI][f];
 256
 257             const vector& faceVelocity = boundaryU[patchI][f];
 258
 259             forAll(pFA, i)
 260             {
 261                 scalar& faceAccumulator = pFA[i][f];
 262
 263                 // Number of whole particles to insert
 264                 label nI = max(label(faceAccumulator), 0);
 265
 266                 // Add another particle with a probability proportional to the
 267                 // remainder of taking the integer part of faceAccumulator
 268                 if ((faceAccumulator - nI) > rndGen.scalar01())
 269                 {
 270                     nI++;
 271                 }
 272
 273                 faceAccumulator -= nI;
 274
 275                 label typeId = moleculeTypeIds_[i];
 276
 277                 scalar mass = cloud.constProps(typeId).mass();
 278
 279                 for (label i = 0; i < nI; i++)
 280                 {
 281                     // Choose a triangle to insert on, based on their relative
 282                     // area
 283
 284                     scalar triSelection = rndGen.scalar01();
 285
 286                     // Selected triangle
 287                     label sTri = -1;
 288
 289                     forAll(cTriAFracs, tri)
 290                     {
 291                         sTri = tri;
 292
 293                         if (cTriAFracs[tri] >= triSelection)
 294                         {
 295                             break;
 296                         }
 297                     }
 298
 299                     // Randomly distribute the points on the triangle, using the
 300                     // method from:
 301                     // Generating Random Points in Triangles
 302                     // by Greg Turk
 303                     // from "Graphics Gems", Academic Press, 1990
 304                     // http://tog.acm.org/GraphicsGems/gems/TriPoints.c
 305
 306                     const point& A = fC;
 307                     const point& B = mesh.points()[faceVertices[sTri]];
 308                     const point& C =
 309                     mesh.points()[faceVertices[(sTri + 1) % nVertices]];
 310
 311                     scalar s = rndGen.scalar01();
 312                     scalar t = sqrt(rndGen.scalar01());
 313
 314                     point p = (1 - t)*A + (1 - s)*t*B + s*t*C;
 315
 316                     // Velocity generation
 317
 318                     scalar mostProbableSpeed
 319                     (
 320                         cloud.maxwellianMostProbableSpeed
 321                         (
 322                             faceTemperature,
 323                             mass
 324                         )
 325                     );
 326
 327                     scalar sCosTheta = (faceVelocity & n)/mostProbableSpeed;
 328
 329                     // Coefficients required for Bird eqn 12.5
 330                     scalar uNormProbCoeffA =
 331                         sCosTheta + sqrt(sqr(sCosTheta) + 2.0);
 332
 333                     scalar uNormProbCoeffB =
 334                         0.5*
 335                         (
 336                             1.0
 337                           + sCosTheta*(sCosTheta - sqrt(sqr(sCosTheta) + 2.0))
 338                         );
 339
 340                     // Equivalent to the QA value in Bird's DSMC3.FOR
 341                     scalar randomScaling = 3.0;
 342
 343                     if (sCosTheta < -3)
 344                     {
 345                         randomScaling = mag(sCosTheta) + 1;
 346                     }
 347
 348                     scalar P = -1;
 349
 350                     // Normalised candidates for the normal direction velocity
 351                     // component
 352                     scalar uNormal;
 353                     scalar uNormalThermal;
 354
 355                     // Select a velocity using Bird eqn 12.5
 356                     do
 357                     {
 358                         uNormalThermal =
 359                             randomScaling*(2.0*rndGen.scalar01() - 1);
 360
 361                         uNormal = uNormalThermal + sCosTheta;
 362
 363                         if (uNormal < 0.0)
 364                         {
 365                             P = -1;
 366                         }
 367                         else
 368                         {
 369                             P = 2.0*uNormal/uNormProbCoeffA
 370                                *exp(uNormProbCoeffB - sqr(uNormalThermal));
 371                         }
 372
 373                     } while (P < rndGen.scalar01());
 374
 375                     vector U =
 376                         sqrt(CloudType::kb*faceTemperature/mass)
 377                        *(
 378                             rndGen.GaussNormal()*t1
 379                           + rndGen.GaussNormal()*t2
 380                         )
 381                       + (t1 & faceVelocity)*t1
 382                       + (t2 & faceVelocity)*t2
 383                       + mostProbableSpeed*uNormal*n;
 384
 385                     scalar Ei = cloud.equipartitionInternalEnergy
 386                     (
 387                         faceTemperature,
 388                         cloud.constProps(typeId).internalDegreesOfFreedom()
 389                     );
 390
 391                     cloud.addNewParcel
 392                     (
 393                         p,
 394                         U,
 395                         Ei,
 396                         cell,
 397                         typeId
 398                     );
 399
 400                     particlesInserted++;
 401                 }
 402             }
 403         }
 404     }
 405
 406     reduce(particlesInserted, sumOp<label>());
 407
 408     Info<< "    Particles inserted              = "
 409         << particlesInserted << endl;
 410
 411 }
 412
 413
 414 // ************************************************************************* //