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) 2011 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>
  36 Foam::FreeStream<CloudType>::FreeStream
  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
  50     DynamicList<label> patches;
  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         {
 100             FatalErrorIn
 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>
 120 Foam::FreeStream<CloudType>::~FreeStream()
 121 {}
 122
 123
 124 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 125
 126 template <class CloudType>
 127 void Foam::FreeStream<CloudType>::inflow()
 128 {
 129     CloudType& cloud(this->owner());
 130
 131     const polyMesh& mesh(cloud.mesh());
 132
 133     const scalar deltaT = mesh.time().deltaTValue();
 134
 135     Random& rndGen(cloud.rndGen());
 136
 137     scalar sqrtPi = sqrt(pi);
 138
 139     label particlesInserted = 0;
 140
 141     const volScalarField::GeometricBoundaryField& boundaryT
 142     (
 143         cloud.boundaryT().boundaryField()
 144     );
 145
 146     const volVectorField::GeometricBoundaryField& boundaryU
 147     (
 148         cloud.boundaryU().boundaryField()
 149     );
 150
 151     forAll(patches_, p)
 152     {
 153         label patchI = patches_[p];
 154
 155         const polyPatch& patch = mesh.boundaryMesh()[patchI];
 156
 157         // Add mass to the accumulators.  negative face area dotted with the
 158         // velocity to point flux into the domain.
 159
 160         // Take a reference to the particleFluxAccumulator for this patch
 161         List<Field<scalar> >& pFA = particleFluxAccumulators_[p];
 162
 163         forAll(pFA, i)
 164         {
 165             label typeId = moleculeTypeIds_[i];
 166
 167             scalar mass = cloud.constProps(typeId).mass();
 168
 169             if (min(boundaryT[patchI]) < SMALL)
 170             {
 171                 FatalErrorIn ("Foam::FreeStream<CloudType>::inflow()")
 172                     << "Zero boundary temperature detected, check boundaryT "
 173                     << "condition." << nl
 174                     << nl << abort(FatalError);
 175             }
 176
 177             scalarField mostProbableSpeed
 178             (
 179                 cloud.maxwellianMostProbableSpeed
 180                 (
 181                     boundaryT[patchI],
 182                     mass
 183                 )
 184             );
 185
 186             // Dotting boundary velocity with the face unit normal
 187             // (which points out of the domain, so it must be
 188             // negated), dividing by the most probable speed to form
 189             // molecularSpeedRatio * cosTheta
 190
 191             scalarField sCosTheta
 192             (
 193                 (boundaryU[patchI] & -patch.faceAreas()/mag(patch.faceAreas()))
 194               / mostProbableSpeed
 195             );
 196
 197             // From Bird eqn 4.22
 198
 199             pFA[i] +=
 200                 mag(patch.faceAreas())*numberDensities_[i]*deltaT
 201                *mostProbableSpeed
 202                *(
 203                    exp(-sqr(sCosTheta)) + sqrtPi*sCosTheta*(1 + erf(sCosTheta))
 204                 )
 205                /(2.0*sqrtPi);
 206         }
 207
 208         forAll(patch, pFI)
 209         {
 210             // Loop over all faces as the outer loop to avoid calculating
 211             // geometrical properties multiple times.
 212
 213             const face& f = patch[pFI];
 214
 215             label globalFaceIndex = pFI + patch.start();
 216
 217             label cellI = mesh.faceOwner()[globalFaceIndex];
 218
 219             const vector& fC = patch.faceCentres()[pFI];
 220
 221             scalar fA = mag(patch.faceAreas()[pFI]);
 222
 223             List<tetIndices> faceTets = polyMeshTetDecomposition::faceTetIndices
 224             (
 225                 mesh,
 226                 globalFaceIndex,
 227                 cellI
 228             );
 229
 230             // Cumulative triangle area fractions
 231             List<scalar> cTriAFracs(faceTets.size(), 0.0);
 232
 233             scalar previousCummulativeSum = 0.0;
 234
 235             forAll(faceTets, triI)
 236             {
 237                 const tetIndices& faceTetIs = faceTets[triI];
 238
 239                 cTriAFracs[triI] =
 240                     faceTetIs.faceTri(mesh).mag()/fA
 241                   + previousCummulativeSum;
 242
 243                 previousCummulativeSum = cTriAFracs[triI];
 244             }
 245
 246             // Force the last area fraction value to 1.0 to avoid any
 247             // rounding/non-flat face errors giving a value < 1.0
 248             cTriAFracs.last() = 1.0;
 249
 250             // Normal unit vector *negative* so normal is pointing into the
 251             // domain
 252             vector n = patch.faceAreas()[pFI];
 253             n /= -mag(n);
 254
 255             // Wall tangential unit vector. Use the direction between the
 256             // face centre and the first vertex in the list
 257             vector t1 = fC - (mesh.points()[f[0]]);
 258             t1 /= mag(t1);
 259
 260             // Other tangential unit vector.  Rescaling in case face is not
 261             // flat and n and t1 aren't perfectly orthogonal
 262             vector t2 = n^t1;
 263             t2 /= mag(t2);
 264
 265             scalar faceTemperature = boundaryT[patchI][pFI];
 266
 267             const vector& faceVelocity = boundaryU[patchI][pFI];
 268
 269             forAll(pFA, i)
 270             {
 271                 scalar& faceAccumulator = pFA[i][pFI];
 272
 273                 // Number of whole particles to insert
 274                 label nI = max(label(faceAccumulator), 0);
 275
 276                 // Add another particle with a probability proportional to the
 277                 // remainder of taking the integer part of faceAccumulator
 278                 if ((faceAccumulator - nI) > rndGen.scalar01())
 279                 {
 280                     nI++;
 281                 }
 282
 283                 faceAccumulator -= nI;
 284
 285                 label typeId = moleculeTypeIds_[i];
 286
 287                 scalar mass = cloud.constProps(typeId).mass();
 288
 289                 for (label i = 0; i < nI; i++)
 290                 {
 291                     // Choose a triangle to insert on, based on their relative
 292                     // area
 293
 294                     scalar triSelection = rndGen.scalar01();
 295
 296                     // Selected triangle
 297                     label selectedTriI = -1;
 298
 299                     forAll(cTriAFracs, triI)
 300                     {
 301                         selectedTriI = triI;
 302
 303                         if (cTriAFracs[triI] >= triSelection)
 304                         {
 305                             break;
 306                         }
 307                     }
 308
 309                     // Randomly distribute the points on the triangle.
 310
 311                     const tetIndices& faceTetIs = faceTets[selectedTriI];
 312
 313                     point p = faceTetIs.faceTri(mesh).randomPoint(rndGen);
 314
 315                     // Velocity generation
 316
 317                     scalar mostProbableSpeed
 318                     (
 319                         cloud.maxwellianMostProbableSpeed
 320                         (
 321                             faceTemperature,
 322                             mass
 323                         )
 324                     );
 325
 326                     scalar sCosTheta = (faceVelocity & n)/mostProbableSpeed;
 327
 328                     // Coefficients required for Bird eqn 12.5
 329                     scalar uNormProbCoeffA =
 330                         sCosTheta + sqrt(sqr(sCosTheta) + 2.0);
 331
 332                     scalar uNormProbCoeffB =
 333                         0.5*
 334                         (
 335                             1.0
 336                           + sCosTheta*(sCosTheta - sqrt(sqr(sCosTheta) + 2.0))
 337                         );
 338
 339                     // Equivalent to the QA value in Bird's DSMC3.FOR
 340                     scalar randomScaling = 3.0;
 341
 342                     if (sCosTheta < -3)
 343                     {
 344                         randomScaling = mag(sCosTheta) + 1;
 345                     }
 346
 347                     scalar P = -1;
 348
 349                     // Normalised candidates for the normal direction velocity
 350                     // component
 351                     scalar uNormal;
 352                     scalar uNormalThermal;
 353
 354                     // Select a velocity using Bird eqn 12.5
 355                     do
 356                     {
 357                         uNormalThermal =
 358                             randomScaling*(2.0*rndGen.scalar01() - 1);
 359
 360                         uNormal = uNormalThermal + sCosTheta;
 361
 362                         if (uNormal < 0.0)
 363                         {
 364                             P = -1;
 365                         }
 366                         else
 367                         {
 368                             P = 2.0*uNormal/uNormProbCoeffA
 369                                *exp(uNormProbCoeffB - sqr(uNormalThermal));
 370                         }
 371
 372                     } while (P < rndGen.scalar01());
 373
 374                     vector U =
 375                         sqrt(physicoChemical::k.value()*faceTemperature/mass)
 376                        *(
 377                             rndGen.GaussNormal()*t1
 378                           + rndGen.GaussNormal()*t2
 379                         )
 380                       + (t1 & faceVelocity)*t1
 381                       + (t2 & faceVelocity)*t2
 382                       + mostProbableSpeed*uNormal*n;
 383
 384                     scalar Ei = cloud.equipartitionInternalEnergy
 385                     (
 386                         faceTemperature,
 387                         cloud.constProps(typeId).internalDegreesOfFreedom()
 388                     );
 389
 390                     cloud.addNewParcel
 391                     (
 392                         p,
 393                         U,
 394                         Ei,
 395                         cellI,
 396                         globalFaceIndex,
 397                         faceTetIs.tetPt(),
 398                         typeId
 399                     );
 400
 401                     particlesInserted++;
 402                 }
 403             }
 404         }
 405     }
 406
 407     reduce(particlesInserted, sumOp<label>());
 408
 409     Info<< "    Particles inserted              = "
 410         << particlesInserted << endl;
 411
 412 }
 413
 414
 415 // ************************************************************************* //