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

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   4    \\    /   O peration     |
   5     \\  /    A nd           | Copyright (C) 2009-2009 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 the
  13     Free Software Foundation; either version 2 of the License, or (at your
  14     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, write to the Free Software Foundation,
  23     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
  24
  25 \*---------------------------------------------------------------------------*/
  26
  27 #include "FreeStream.H"
  28
  29 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
  30
  31 template <class CloudType>
  32 Foam::FreeStream<CloudType>::FreeStream
  33 (
  34     const dictionary& dict,
  35     CloudType& cloud
  36 )
  37 :
  38     InflowBoundaryModel<CloudType>(dict, cloud, typeName),
  39     patches_(),
  40     moleculeTypeIds_(),
  41     numberDensities_(),
  42     particleFluxAccumulators_()
  43 {
  44     // Identify which patches to use
  45
  46     DynamicList<label> patches;
  47
  48     forAll(cloud.mesh().boundaryMesh(), p)
  49     {
  50         const polyPatch& patch = cloud.mesh().boundaryMesh()[p];
  51
  52         if (patch.type() == polyPatch::typeName)
  53         {
  54             patches.append(p);
  55         }
  56     }
  57
  58     patches_.transfer(patches);
  59
  60     const dictionary& numberDensitiesDict
  61     (
  62         this->coeffDict().subDict("numberDensities")
  63     );
  64
  65     List<word> molecules(numberDensitiesDict.toc());
  66
  67     // Initialise the particleFluxAccumulators_
  68     particleFluxAccumulators_.setSize(patches_.size());
  69
  70     forAll(patches_, p)
  71     {
  72         const polyPatch& patch = cloud.mesh().boundaryMesh()[patches_[p]];
  73
  74         particleFluxAccumulators_[p] = List<Field<scalar> >
  75         (
  76             molecules.size(),
  77             Field<scalar>(patch.size(), 0.0)
  78         );
  79     }
  80
  81     moleculeTypeIds_.setSize(molecules.size());
  82
  83     numberDensities_.setSize(molecules.size());
  84
  85     forAll(molecules, i)
  86     {
  87         numberDensities_[i] = readScalar
  88         (
  89             numberDensitiesDict.lookup(molecules[i])
  90         );
  91
  92         moleculeTypeIds_[i] = findIndex(cloud.typeIdList(), molecules[i]);
  93
  94         if (moleculeTypeIds_[i] == -1)
  95         {
  96             FatalErrorIn
  97             (
  98                 "Foam::FreeStream<CloudType>::FreeStream"
  99                 "("
 100                     "const dictionary&, "
 101                     "CloudType&"
 102                 ")"
 103             )   << "typeId " << molecules[i] << "not defined in cloud." << nl
 104                 << abort(FatalError);
 105         }
 106     }
 107
 108     numberDensities_ /= cloud.nParticle();
 109 }
 110
 111
 112
 113 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 114
 115 template <class CloudType>
 116 Foam::FreeStream<CloudType>::~FreeStream()
 117 {}
 118
 119
 120 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 121
 122 template <class CloudType>
 123 void Foam::FreeStream<CloudType>::inflow()
 124 {
 125     CloudType& cloud(this->owner());
 126
 127     const polyMesh& mesh(cloud.mesh());
 128
 129     const scalar deltaT = cloud.cachedDeltaT();
 130
 131     Random& rndGen(cloud.rndGen());
 132
 133     scalar sqrtPi = sqrt(mathematicalConstant::pi);
 134
 135     label particlesInserted = 0;
 136
 137     const volScalarField::GeometricBoundaryField& boundaryT
 138     (
 139         cloud.T().boundaryField()
 140     );
 141
 142     const volVectorField::GeometricBoundaryField& boundaryU
 143     (
 144         cloud.U().boundaryField()
 145     );
 146
 147     forAll(patches_, p)
 148     {
 149         label patchI = patches_[p];
 150
 151         const polyPatch& patch = mesh.boundaryMesh()[patchI];
 152
 153         // Add mass to the accumulators.  negative face area dotted with the
 154         // velocity to point flux into the domain.
 155
 156         // Take a reference to the particleFluxAccumulator for this patch
 157         List<Field<scalar> >& pFA = particleFluxAccumulators_[p];
 158
 159         forAll(pFA, i)
 160         {
 161             label typeId = moleculeTypeIds_[i];
 162
 163             scalar mass = cloud.constProps(typeId).mass();
 164
 165             if (min(boundaryT[patchI]) < SMALL)
 166             {
 167                 FatalErrorIn ("Foam::FreeStream<CloudType>::inflow()")
 168                     << "Zero boundary temperature detected, 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             for (label v = 0; v < nVertices - 1; v++)
 222             {
 223                 const point& vA = mesh.points()[faceVertices[v]];
 224
 225                 const point& vB = mesh.points()[faceVertices[(v + 1)]];
 226
 227                 cTriAFracs[v] =
 228                     0.5*mag((vA - fC)^(vB - fC))/fA
 229                   + cTriAFracs[max((v - 1), 0)];
 230             }
 231
 232             // Force the last area fraction value to 1.0 to avoid any
 233             // rounding/non-flat face errors giving a value < 1.0
 234             cTriAFracs[nVertices - 1] = 1.0;
 235
 236             // Normal unit vector *negative* so normal is pointing into the
 237             // domain
 238             vector n = patch.faceAreas()[f];
 239             n /= -mag(n);
 240
 241             // Wall tangential unit vector. Use the direction between the
 242             // face centre and the first vertex in the list
 243             vector t1 = fC - (mesh.points()[faceVertices[0]]);
 244             t1 /= mag(t1);
 245
 246             // Other tangential unit vector.  Rescaling in case face is not
 247             // flat and n and t1 aren't perfectly orthogonal
 248             vector t2 = n^t1;
 249             t2 /= mag(t2);
 250
 251             scalar faceTemperature = boundaryT[patchI][f];
 252
 253             const vector& faceVelocity = boundaryU[patchI][f];
 254
 255             forAll(pFA, i)
 256             {
 257                 scalar& faceAccumulator = pFA[i][f];
 258
 259                 // Number of whole particles to insert
 260                 label nI = max(label(faceAccumulator), 0);
 261
 262                 // Add another particle with a probability proportional to the
 263                 // remainder of taking the integer part of faceAccumulator
 264                 if ((faceAccumulator - nI) > rndGen.scalar01())
 265                 {
 266                     nI++;
 267                 }
 268
 269                 faceAccumulator -= nI;
 270
 271                 label typeId = moleculeTypeIds_[i];
 272
 273                 scalar mass = cloud.constProps(typeId).mass();
 274
 275                 for (label i = 0; i < nI; i++)
 276                 {
 277                     // Choose a triangle to insert on, based on their relative
 278                     // area
 279
 280                     scalar triSelection = rndGen.scalar01();
 281
 282                     // Selected triangle
 283                     label sTri = -1;
 284
 285                     forAll(cTriAFracs, tri)
 286                     {
 287                         sTri = tri;
 288
 289                         if (cTriAFracs[tri] >= triSelection)
 290                         {
 291                             break;
 292                         }
 293                     }
 294
 295                     // Randomly distribute the points on the triangle, using the
 296                     // method from:
 297                     // Generating Random Points in Triangles
 298                     // by Greg Turk
 299                     // from "Graphics Gems", Academic Press, 1990
 300                     // http://tog.acm.org/GraphicsGems/gems/TriPoints.c
 301
 302                     const point& A = fC;
 303                     const point& B = mesh.points()[faceVertices[sTri]];
 304                     const point& C =
 305                     mesh.points()[faceVertices[(sTri + 1) % nVertices]];
 306
 307                     scalar s = rndGen.scalar01();
 308                     scalar t = sqrt(rndGen.scalar01());
 309
 310                     point p = (1 - t)*A + (1 - s)*t*B + s*t*C;
 311
 312                     // Velocity generation
 313
 314                     scalar mostProbableSpeed
 315                     (
 316                         cloud.maxwellianMostProbableSpeed
 317                         (
 318                             faceTemperature,
 319                             mass
 320                         )
 321                     );
 322
 323                     scalar sCosTheta = (faceVelocity & n)/mostProbableSpeed;
 324
 325                     // Coefficients required for Bird eqn 12.5
 326                     scalar uNormProbCoeffA =
 327                         sCosTheta + sqrt(sqr(sCosTheta) + 2.0);
 328
 329                     scalar uNormProbCoeffB =
 330                         0.5*
 331                         (
 332                             1.0
 333                           + sCosTheta*(sCosTheta - sqrt(sqr(sCosTheta) + 2.0))
 334                         );
 335
 336                     // Equivalent to the QA value in Bird's DSMC3.FOR
 337                     scalar randomScaling = 3.0;
 338
 339                     if (sCosTheta < -3)
 340                     {
 341                         randomScaling = mag(sCosTheta) + 1;
 342                     }
 343
 344                     scalar P = -1;
 345
 346                     // Normalised candidates for the normal direction velocity
 347                     // component
 348                     scalar uNormal;
 349                     scalar uNormalThermal;
 350
 351                     // Select a velocity using Bird eqn 12.5
 352                     do
 353                     {
 354                         uNormalThermal =
 355                             randomScaling*(2.0*rndGen.scalar01() - 1);
 356
 357                         uNormal = uNormalThermal + sCosTheta;
 358
 359                         if (uNormal < 0.0)
 360                         {
 361                             P = -1;
 362                         }
 363                         else
 364                         {
 365                             P = 2.0*uNormal/uNormProbCoeffA
 366                                *exp(uNormProbCoeffB - sqr(uNormalThermal));
 367                         }
 368
 369                     } while (P < rndGen.scalar01());
 370
 371                     vector U =
 372                         sqrt(CloudType::kb*faceTemperature/mass)
 373                        *(
 374                             rndGen.GaussNormal()*t1
 375                           + rndGen.GaussNormal()*t2
 376                         )
 377                       + (t1 & faceVelocity)*t1
 378                       + (t2 & faceVelocity)*t2
 379                       + mostProbableSpeed*uNormal*n;
 380
 381                     scalar Ei = cloud.equipartitionInternalEnergy
 382                     (
 383                         faceTemperature,
 384                         cloud.constProps(typeId).internalDegreesOfFreedom()
 385                     );
 386
 387                     cloud.addNewParcel
 388                     (
 389                         p,
 390                         U,
 391                         Ei,
 392                         cell,
 393                         typeId
 394                     );
 395
 396                     particlesInserted++;
 397                 }
 398             }
 399         }
 400     }
 401
 402     reduce(particlesInserted, sumOp<label>());
 403
 404     Info<< "    Particles inserted              = "
 405         << particlesInserted << endl;
 406
 407 }
 408
 409
 410 // ************************************************************************* //