src/lagrangian/dsmc/clouds/Templates/DsmcCloud/DsmcCloudTemplate.C

   1 /*---------------------------------------------------------------------------*\
   2   =========                 |
   3   \\      /  F ield         | foam-extend: Open Source CFD
   4    \\    /   O peration     | Version:     3.2
   5     \\  /    A nd           | Web:         http://www.foam-extend.org
   6      \\/     M anipulation  | For copyright notice see file Copyright
   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 "DsmcCloudTemplate.H"
  27 #include "BinaryCollisionModel.H"
  28 #include "WallInteractionModel.H"
  29 #include "InflowBoundaryModel.H"
  30
  31 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
  32
  33 template<class ParcelType>
  34 Foam::scalar Foam::DsmcCloud<ParcelType>::kb = 1.380650277e-23;
  35
  36
  37 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
  38
  39 template<class ParcelType>
  40 void Foam::DsmcCloud<ParcelType>::buildConstProps()
  41 {
  42     Info<< nl << "Constructing constant properties for" << endl;
  43     constProps_.setSize(typeIdList_.size());
  44
  45     dictionary moleculeProperties
  46     (
  47         particleProperties_.subDict("moleculeProperties")
  48     );
  49
  50     forAll(typeIdList_, i)
  51     {
  52         const word& id(typeIdList_[i]);
  53
  54         Info<< "    " << id << endl;
  55
  56         const dictionary& molDict(moleculeProperties.subDict(id));
  57
  58         constProps_[i] =
  59         typename ParcelType::constantProperties::constantProperties(molDict);
  60     }
  61 }
  62
  63
  64 template<class ParcelType>
  65 void Foam::DsmcCloud<ParcelType>::buildCellOccupancy()
  66 {
  67     forAll(cellOccupancy_, cO)
  68     {
  69         cellOccupancy_[cO].clear();
  70     }
  71
  72     forAllIter(typename DsmcCloud<ParcelType>, *this, iter)
  73     {
  74         cellOccupancy_[iter().cell()].append(&iter());
  75     }
  76 }
  77
  78
  79 template<class ParcelType>
  80 void Foam::DsmcCloud<ParcelType>::initialise
  81 (
  82     const IOdictionary& dsmcInitialiseDict
  83 )
  84 {
  85     Info<< nl << "Initialising particles" << endl;
  86
  87     const scalar temperature
  88     (
  89         readScalar(dsmcInitialiseDict.lookup("temperature"))
  90     );
  91
  92     const vector velocity(dsmcInitialiseDict.lookup("velocity"));
  93
  94     const dictionary& numberDensitiesDict
  95     (
  96         dsmcInitialiseDict.subDict("numberDensities")
  97     );
  98
  99     List<word> molecules(numberDensitiesDict.toc());
 100
 101     Field<scalar> numberDensities(molecules.size());
 102
 103     forAll(molecules, i)
 104     {
 105         numberDensities[i] = readScalar
 106         (
 107             numberDensitiesDict.lookup(molecules[i])
 108         );
 109     }
 110
 111     numberDensities /= nParticle_;
 112
 113     forAll(mesh_.cells(), cell)
 114     {
 115         const vector& cC = mesh_.cellCentres()[cell];
 116         const labelList& cellFaces = mesh_.cells()[cell];
 117         const scalar cV = mesh_.cellVolumes()[cell];
 118
 119         label nTets = 0;
 120
 121         // Each face is split into nEdges (or nVertices) - 2 tets.
 122         forAll(cellFaces, face)
 123         {
 124             nTets += mesh_.faces()[cellFaces[face]].size() - 2;
 125         }
 126
 127         // Calculate the cumulative tet volumes circulating around the cell and
 128         // record the vertex labels of each.
 129         scalarList cTetVFracs(nTets, 0.0);
 130
 131         List<labelList> tetPtIs(nTets, labelList(3,-1));
 132
 133         // Keep track of which tet this is.
 134         label tet = 0;
 135
 136         forAll(cellFaces, face)
 137         {
 138             const labelList& facePoints = mesh_.faces()[cellFaces[face]];
 139
 140             label pointI = 1;
 141             while ((pointI + 1) < facePoints.size())
 142             {
 143
 144                 const vector& pA = mesh_.points()[facePoints[0]];
 145                 const vector& pB = mesh_.points()[facePoints[pointI]];
 146                 const vector& pC = mesh_.points()[facePoints[pointI + 1]];
 147
 148                 cTetVFracs[tet] =
 149                     mag(((pA - cC) ^ (pB - cC)) & (pC - cC))/(cV*6.0)
 150                   + cTetVFracs[max((tet - 1),0)];
 151
 152                 tetPtIs[tet][0] = facePoints[0];
 153                 tetPtIs[tet][1] = facePoints[pointI];
 154                 tetPtIs[tet][2] = facePoints[pointI + 1];
 155
 156                 pointI++;
 157                 tet++;
 158             }
 159         }
 160
 161         // Force the last volume fraction value to 1.0 to avoid any
 162         // rounding/non-flat face errors giving a value < 1.0
 163         cTetVFracs[nTets - 1] = 1.0;
 164
 165         forAll(molecules, i)
 166         {
 167             const word& moleculeName(molecules[i]);
 168
 169             label typeId(findIndex(typeIdList_, moleculeName));
 170
 171             if (typeId == -1)
 172             {
 173                 FatalErrorIn("Foam::DsmcCloud<ParcelType>::initialise")
 174                 << "typeId " << moleculeName << "not defined." << nl
 175                     << abort(FatalError);
 176             }
 177
 178             const typename ParcelType::constantProperties& cP =
 179                 constProps(typeId);
 180
 181             scalar numberDensity = numberDensities[i];
 182
 183             // Calculate the number of particles required
 184             scalar particlesRequired = numberDensity*mesh_.cellVolumes()[cell];
 185
 186             // Only integer numbers of particles can be inserted
 187             label nParticlesToInsert = label(particlesRequired);
 188
 189             // Add another particle with a probability proportional to the
 190             // remainder of taking the integer part of particlesRequired
 191             if ((particlesRequired - nParticlesToInsert) > rndGen_.scalar01())
 192             {
 193                 nParticlesToInsert++;
 194             }
 195
 196             for (label pI = 0; pI < nParticlesToInsert; pI++)
 197             {
 198                 // Choose a random point in a generic tetrahedron
 199
 200                 scalar s = rndGen_.scalar01();
 201                 scalar t = rndGen_.scalar01();
 202                 scalar u = rndGen_.scalar01();
 203
 204                 if (s + t > 1.0)
 205                 {
 206                     s = 1.0 - s;
 207                     t = 1.0 - t;
 208                 }
 209
 210                 if (t + u > 1.0)
 211                 {
 212                     scalar tmp = u;
 213                     u = 1.0 - s - t;
 214                     t = 1.0 - tmp;
 215                 }
 216                 else if (s + t + u > 1.0)
 217                 {
 218                     scalar tmp = u;
 219                     u = s + t + u - 1.0;
 220                     s = 1.0 - t - tmp;
 221                 }
 222
 223                 // Choose a tetrahedron to insert in, based on their relative
 224                 // volumes
 225                 scalar tetSelection = rndGen_.scalar01();
 226
 227                 // Selected tetrahedron
 228                 label sTet = -1;
 229
 230                 forAll(cTetVFracs, tet)
 231                 {
 232                     sTet = tet;
 233
 234                     if (cTetVFracs[tet] >= tetSelection)
 235                     {
 236                         break;
 237                     }
 238                 }
 239
 240                 vector p =
 241                     (1 - s - t - u)*cC
 242                   + s*mesh_.points()[tetPtIs[sTet][0]]
 243                   + t*mesh_.points()[tetPtIs[sTet][1]]
 244                   + u*mesh_.points()[tetPtIs[sTet][2]];
 245
 246                 vector U = equipartitionLinearVelocity
 247                 (
 248                     temperature,
 249                     cP.mass()
 250                 );
 251
 252                 scalar Ei = equipartitionInternalEnergy
 253                 (
 254                     temperature,
 255                     cP.internalDegreesOfFreedom()
 256                 );
 257
 258                 U += velocity;
 259
 260                 addNewParcel
 261                 (
 262                     p,
 263                     U,
 264                     Ei,
 265                     cell,
 266                     typeId
 267                 );
 268             }
 269         }
 270     }
 271
 272     // Initialise the sigmaTcRMax_ field to the product of the cross section of
 273     // the most abundant species and the most probable thermal speed (Bird,
 274     // p222-223)
 275
 276     label mostAbundantType(findMax(numberDensities));
 277
 278     const typename ParcelType::constantProperties& cP = constProps
 279     (
 280         mostAbundantType
 281     );
 282
 283     sigmaTcRMax_.internalField() = cP.sigmaT()*maxwellianMostProbableSpeed
 284     (
 285         temperature,
 286         cP.mass()
 287     );
 288
 289     sigmaTcRMax_.correctBoundaryConditions();
 290 }
 291
 292
 293 template<class ParcelType>
 294 void Foam::DsmcCloud<ParcelType>::collisions()
 295 {
 296     buildCellOccupancy();
 297
 298     // Temporary storage for subCells
 299     List<dynamicLabelList > subCells(8);
 300
 301     scalar deltaT = mesh().time().deltaTValue();
 302
 303     label collisionCandidates = 0;
 304
 305     label collisions = 0;
 306
 307     forAll(cellOccupancy_, celli)
 308     {
 309         const DynamicList<ParcelType*>& cellParcels(cellOccupancy_[celli]);
 310
 311         label nC(cellParcels.size());
 312
 313         if (nC > 1)
 314         {
 315
 316             // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 317             // Assign particles to one of 8 Cartesian subCells
 318
 319             // Clear temporary lists
 320             forAll(subCells, i)
 321             {
 322                 subCells[i].clear();
 323             }
 324
 325             // Inverse addressing specifying which subCell a parcel is in
 326             List<label> whichSubCell(cellParcels.size());
 327
 328             const point& cC = mesh_.cellCentres()[celli];
 329
 330             forAll(cellParcels, i)
 331             {
 332                 ParcelType* p = cellParcels[i];
 333
 334                 vector relPos = p->position() - cC;
 335
 336                 label subCell =
 337                     pos(relPos.x()) + 2*pos(relPos.y()) + 4*pos(relPos.z());
 338
 339                 subCells[subCell].append(i);
 340
 341                 whichSubCell[i] = subCell;
 342             }
 343
 344             // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 345
 346             scalar sigmaTcRMax = sigmaTcRMax_[celli];
 347
 348             scalar selectedPairs =
 349                 collisionSelectionRemainder_[celli]
 350               + 0.5*nC*(nC - 1)*nParticle_*sigmaTcRMax*deltaT
 351                /mesh_.cellVolumes()[celli];
 352
 353             label nCandidates(selectedPairs);
 354
 355             collisionSelectionRemainder_[celli] = selectedPairs - nCandidates;
 356
 357             collisionCandidates += nCandidates;
 358
 359             for (label c = 0; c < nCandidates; c++)
 360             {
 361                 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 362                 // subCell candidate selection procedure
 363
 364                 // Select the first collision candidate
 365                 label candidateP = rndGen_.integer(0, nC - 1);
 366
 367                 // Declare the second collision candidate
 368                 label candidateQ = -1;
 369
 370                 List<label> subCellPs = subCells[whichSubCell[candidateP]];
 371
 372                 label nSC = subCellPs.size();
 373
 374                 if (nSC > 1)
 375                 {
 376                     // If there are two or more particle in a subCell, choose
 377                     // another from the same cell.  If the same candidate is
 378                     // chosen, choose again.
 379
 380                     do
 381                     {
 382                         candidateQ = subCellPs[rndGen_.integer(0, nSC - 1)];
 383
 384                     } while(candidateP == candidateQ);
 385                 }
 386                 else
 387                 {
 388                     // Select a possible second collision candidate from the
 389                     // whole cell.  If the same candidate is chosen, choose
 390                     // again.
 391
 392                     do
 393                     {
 394                         candidateQ = rndGen_.integer(0, nC - 1);
 395
 396                     } while(candidateP == candidateQ);
 397                 }
 398
 399                 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 400                 // uniform candidate selection procedure
 401
 402                 // // Select the first collision candidate
 403                 // label candidateP = rndGen_.integer(0, nC-1);
 404
 405                 // // Select a possible second collision candidate
 406                 // label candidateQ = rndGen_.integer(0, nC-1);
 407
 408                 // // If the same candidate is chosen, choose again
 409                 // while(candidateP == candidateQ)
 410                 // {
 411                 //     candidateQ = rndGen_.integer(0, nC-1);
 412                 // }
 413
 414                 // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 415
 416                 ParcelType* parcelP = cellParcels[candidateP];
 417                 ParcelType* parcelQ = cellParcels[candidateQ];
 418
 419                 label typeIdP = parcelP->typeId();
 420                 label typeIdQ = parcelQ->typeId();
 421
 422                 scalar sigmaTcR = binaryCollision().sigmaTcR
 423                 (
 424                     typeIdP,
 425                     typeIdQ,
 426                     parcelP->U(),
 427                     parcelQ->U()
 428                 );
 429
 430                 // Update the maximum value of sigmaTcR stored, but use the
 431                 // initial value in the acceptance-rejection criteria because
 432                 // the number of collision candidates selected was based on this
 433
 434                 if (sigmaTcR > sigmaTcRMax_[celli])
 435                 {
 436                     sigmaTcRMax_[celli] = sigmaTcR;
 437                 }
 438
 439                 if ((sigmaTcR/sigmaTcRMax) > rndGen_.scalar01())
 440                 {
 441                     binaryCollision().collide
 442                     (
 443                         typeIdP,
 444                         typeIdQ,
 445                         parcelP->U(),
 446                         parcelQ->U(),
 447                         parcelP->Ei(),
 448                         parcelQ->Ei()
 449                     );
 450
 451                     collisions++;
 452                 }
 453             }
 454         }
 455     }
 456
 457     reduce(collisions, sumOp<label>());
 458
 459     reduce(collisionCandidates, sumOp<label>());
 460
 461     if (collisionCandidates)
 462     {
 463         Info<< "    Collisions                      = "
 464             << collisions << nl
 465             << "    Acceptance rate                 = "
 466             << scalar(collisions)/scalar(collisionCandidates) << nl
 467             << endl;
 468     }
 469     else
 470     {
 471         Info<< "    No collisions" << endl;
 472     }
 473 }
 474
 475
 476 template<class ParcelType>
 477 void Foam::DsmcCloud<ParcelType>::resetFields()
 478 {
 479     q_ = dimensionedScalar("zero",  dimensionSet(1, 0, -3, 0, 0), 0.0);
 480
 481     fD_ = dimensionedVector
 482     (
 483         "zero",
 484         dimensionSet(1, -1, -2, 0, 0),
 485         vector::zero
 486     );
 487
 488     rhoN_ = dimensionedScalar("zero",  dimensionSet(0, -3, 0, 0, 0), VSMALL);
 489
 490     rhoM_ =  dimensionedScalar("zero",  dimensionSet(1, -3, 0, 0, 0), VSMALL);
 491
 492     dsmcRhoN_ = dimensionedScalar("zero",  dimensionSet(0, -3, 0, 0, 0), 0.0);
 493
 494     linearKE_ = dimensionedScalar("zero",  dimensionSet(1, -1, -2, 0, 0), 0.0);
 495
 496     internalE_ = dimensionedScalar("zero",  dimensionSet(1, -1, -2, 0, 0), 0.0);
 497
 498     iDof_ = dimensionedScalar("zero",  dimensionSet(0, -3, 0, 0, 0), VSMALL);
 499
 500     momentum_ = dimensionedVector
 501     (
 502         "zero",
 503         dimensionSet(1, -2, -1, 0, 0),
 504         vector::zero
 505     );
 506 }
 507
 508
 509 template<class ParcelType>
 510 void Foam::DsmcCloud<ParcelType>::calculateFields()
 511 {
 512     scalarField& rhoN = rhoN_.internalField();
 513
 514     scalarField& rhoM = rhoM_.internalField();
 515
 516     scalarField& dsmcRhoN = dsmcRhoN_.internalField();
 517
 518     scalarField& linearKE = linearKE_.internalField();
 519
 520     scalarField& internalE = internalE_.internalField();
 521
 522     scalarField& iDof = iDof_.internalField();
 523
 524     vectorField& momentum = momentum_.internalField();
 525
 526     forAllConstIter(typename DsmcCloud<ParcelType>, *this, iter)
 527     {
 528         const ParcelType& p = iter();
 529         const label cellI = p.cell();
 530
 531         rhoN[cellI]++;
 532
 533         rhoM[cellI] += constProps(p.typeId()).mass();
 534
 535         dsmcRhoN[cellI]++;
 536
 537         linearKE[cellI] += 0.5*constProps(p.typeId()).mass()*(p.U() & p.U());
 538
 539         internalE[cellI] += p.Ei();
 540
 541         iDof[cellI] += constProps(p.typeId()).internalDegreesOfFreedom();
 542
 543         momentum[cellI] += constProps(p.typeId()).mass()*p.U();
 544     }
 545
 546     rhoN *= nParticle_/mesh().cellVolumes();
 547     rhoN_.correctBoundaryConditions();
 548
 549     rhoM *= nParticle_/mesh().cellVolumes();
 550     rhoM_.correctBoundaryConditions();
 551
 552     linearKE *= nParticle_/mesh().cellVolumes();
 553     linearKE_.correctBoundaryConditions();
 554
 555     internalE *= nParticle_/mesh().cellVolumes();
 556     internalE_.correctBoundaryConditions();
 557
 558     iDof *= nParticle_/mesh().cellVolumes();
 559     iDof_.correctBoundaryConditions();
 560
 561     momentum *= nParticle_/mesh().cellVolumes();
 562     momentum_.correctBoundaryConditions();
 563 }
 564
 565
 566 // * * * * * * * * * * * Protected Member Functions  * * * * * * * * * * * * //
 567
 568 template<class ParcelType>
 569 void Foam::DsmcCloud<ParcelType>::addNewParcel
 570 (
 571     const vector& position,
 572     const vector& U,
 573     const scalar Ei,
 574     const label cellId,
 575     const label typeId
 576 )
 577 {
 578     ParcelType* pPtr = new ParcelType
 579     (
 580         *this,
 581         position,
 582         U,
 583         Ei,
 584         cellId,
 585         typeId
 586     );
 587
 588     this->addParticle(pPtr);
 589 }
 590
 591
 592 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 593
 594 template<class ParcelType>
 595 Foam::DsmcCloud<ParcelType>::DsmcCloud
 596 (
 597     const word& cloudName,
 598     const fvMesh& mesh,
 599     bool readFields
 600 )
 601 :
 602     Cloud<ParcelType>(mesh, cloudName, false),
 603     DsmcBaseCloud(),
 604     cloudName_(cloudName),
 605     mesh_(mesh),
 606     particleProperties_
 607     (
 608         IOobject
 609         (
 610             cloudName + "Properties",
 611             mesh_.time().constant(),
 612             mesh_,
 613             IOobject::MUST_READ,
 614             IOobject::NO_WRITE
 615         )
 616     ),
 617     typeIdList_(particleProperties_.lookup("typeIdList")),
 618     nParticle_(readScalar(particleProperties_.lookup("nEquivalentParticles"))),
 619     cellOccupancy_(mesh_.nCells()),
 620     sigmaTcRMax_
 621     (
 622         IOobject
 623         (
 624             this->name() + "SigmaTcRMax",
 625             mesh_.time().timeName(),
 626             mesh_,
 627             IOobject::MUST_READ,
 628             IOobject::AUTO_WRITE
 629         ),
 630         mesh_
 631     ),
 632     collisionSelectionRemainder_(mesh_.nCells(), 0),
 633     q_
 634     (
 635         IOobject
 636         (
 637             "q",
 638             mesh_.time().timeName(),
 639             mesh_,
 640             IOobject::MUST_READ,
 641             IOobject::AUTO_WRITE
 642         ),
 643         mesh_
 644     ),
 645     fD_
 646     (
 647         IOobject
 648         (
 649             "fD",
 650             mesh_.time().timeName(),
 651             mesh_,
 652             IOobject::MUST_READ,
 653             IOobject::AUTO_WRITE
 654         ),
 655         mesh_
 656     ),
 657     rhoN_
 658     (
 659         IOobject
 660         (
 661             "rhoN",
 662             mesh_.time().timeName(),
 663             mesh_,
 664             IOobject::MUST_READ,
 665             IOobject::AUTO_WRITE
 666         ),
 667         mesh_
 668     ),
 669     rhoM_
 670     (
 671         IOobject
 672         (
 673             "rhoM",
 674             mesh_.time().timeName(),
 675             mesh_,
 676             IOobject::MUST_READ,
 677             IOobject::AUTO_WRITE
 678         ),
 679         mesh_
 680     ),
 681     dsmcRhoN_
 682     (
 683         IOobject
 684         (
 685             "dsmcRhoN",
 686             mesh_.time().timeName(),
 687             mesh_,
 688             IOobject::MUST_READ,
 689             IOobject::AUTO_WRITE
 690         ),
 691         mesh_
 692     ),
 693     linearKE_
 694     (
 695         IOobject
 696         (
 697             "linearKE",
 698             mesh_.time().timeName(),
 699             mesh_,
 700             IOobject::MUST_READ,
 701             IOobject::AUTO_WRITE
 702         ),
 703         mesh_
 704     ),
 705     internalE_
 706     (
 707         IOobject
 708         (
 709             "internalE",
 710             mesh_.time().timeName(),
 711             mesh_,
 712             IOobject::MUST_READ,
 713             IOobject::AUTO_WRITE
 714         ),
 715         mesh_
 716     ),
 717     iDof_
 718     (
 719         IOobject
 720         (
 721             "iDof",
 722             mesh_.time().timeName(),
 723             mesh_,
 724             IOobject::MUST_READ,
 725             IOobject::AUTO_WRITE
 726         ),
 727         mesh_
 728     ),
 729     momentum_
 730     (
 731         IOobject
 732         (
 733             "momentum",
 734             mesh_.time().timeName(),
 735             mesh_,
 736             IOobject::MUST_READ,
 737             IOobject::AUTO_WRITE
 738         ),
 739         mesh_
 740     ),
 741     constProps_(),
 742     rndGen_(label(149382906) + 7183*Pstream::myProcNo()),
 743     boundaryT_
 744     (
 745         volScalarField
 746         (
 747             IOobject
 748             (
 749                 "boundaryT",
 750                 mesh_.time().timeName(),
 751                 mesh_,
 752                 IOobject::MUST_READ,
 753                 IOobject::AUTO_WRITE
 754             ),
 755             mesh_
 756         )
 757     ),
 758     boundaryU_
 759     (
 760         volVectorField
 761         (
 762             IOobject
 763             (
 764                 "boundaryU",
 765                 mesh_.time().timeName(),
 766                 mesh_,
 767                 IOobject::MUST_READ,
 768                 IOobject::AUTO_WRITE
 769             ),
 770             mesh_
 771         )
 772     ),
 773     binaryCollisionModel_
 774     (
 775         BinaryCollisionModel<DsmcCloud<ParcelType> >::New
 776         (
 777             particleProperties_,
 778             *this
 779         )
 780     ),
 781     wallInteractionModel_
 782     (
 783         WallInteractionModel<DsmcCloud<ParcelType> >::New
 784         (
 785             particleProperties_,
 786             *this
 787         )
 788     ),
 789     inflowBoundaryModel_
 790     (
 791         InflowBoundaryModel<DsmcCloud<ParcelType> >::New
 792         (
 793             particleProperties_,
 794             *this
 795         )
 796     )
 797 {
 798     buildConstProps();
 799
 800     buildCellOccupancy();
 801
 802     // Initialise the collision selection remainder to a random value between 0
 803     // and 1.
 804     forAll(collisionSelectionRemainder_, i)
 805     {
 806         collisionSelectionRemainder_[i] = rndGen_.scalar01();
 807     }
 808
 809     if (readFields)
 810     {
 811         ParcelType::readFields(*this);
 812     }
 813 }
 814
 815
 816 template<class ParcelType>
 817 Foam::DsmcCloud<ParcelType>::DsmcCloud
 818 (
 819     const word& cloudName,
 820     const fvMesh& mesh,
 821     const IOdictionary& dsmcInitialiseDict
 822 )
 823     :
 824     Cloud<ParcelType>(mesh, cloudName, false),
 825     DsmcBaseCloud(),
 826     cloudName_(cloudName),
 827     mesh_(mesh),
 828     particleProperties_
 829     (
 830         IOobject
 831         (
 832             cloudName + "Properties",
 833             mesh_.time().constant(),
 834             mesh_,
 835             IOobject::MUST_READ,
 836             IOobject::NO_WRITE
 837         )
 838     ),
 839     typeIdList_(particleProperties_.lookup("typeIdList")),
 840     nParticle_(readScalar(particleProperties_.lookup("nEquivalentParticles"))),
 841     cellOccupancy_(),
 842     sigmaTcRMax_
 843     (
 844         IOobject
 845         (
 846             this->name() + "SigmaTcRMax",
 847             mesh_.time().timeName(),
 848             mesh_,
 849             IOobject::NO_READ,
 850             IOobject::AUTO_WRITE
 851         ),
 852         mesh_,
 853         dimensionedScalar("zero",  dimensionSet(0, 3, -1, 0, 0), 0.0)
 854     ),
 855     collisionSelectionRemainder_(),
 856     q_
 857     (
 858         IOobject
 859         (
 860             this->name() + "q_",
 861             mesh_.time().timeName(),
 862             mesh_,
 863             IOobject::NO_READ,
 864             IOobject::NO_WRITE
 865         ),
 866         mesh_,
 867         dimensionedScalar("zero",  dimensionSet(1, 0, -3, 0, 0), 0.0)
 868     ),
 869     fD_
 870     (
 871         IOobject
 872         (
 873             this->name() + "fD_",
 874             mesh_.time().timeName(),
 875             mesh_,
 876             IOobject::NO_READ,
 877             IOobject::NO_WRITE
 878         ),
 879         mesh_,
 880         dimensionedVector
 881         (
 882             "zero",
 883             dimensionSet(1, -1, -2, 0, 0),
 884             vector::zero
 885         )
 886     ),
 887     rhoN_
 888     (
 889         IOobject
 890         (
 891             this->name() + "rhoN_",
 892             mesh_.time().timeName(),
 893             mesh_,
 894             IOobject::NO_READ,
 895             IOobject::NO_WRITE
 896         ),
 897         mesh_,
 898         dimensionedScalar("zero",  dimensionSet(0, -3, 0, 0, 0), VSMALL)
 899     ),
 900     rhoM_
 901     (
 902         IOobject
 903         (
 904             this->name() + "rhoM_",
 905             mesh_.time().timeName(),
 906             mesh_,
 907             IOobject::NO_READ,
 908             IOobject::NO_WRITE
 909         ),
 910         mesh_,
 911         dimensionedScalar("zero",  dimensionSet(1, -3, 0, 0, 0), VSMALL)
 912     ),
 913     dsmcRhoN_
 914     (
 915         IOobject
 916         (
 917             this->name() + "dsmcRhoN_",
 918             mesh_.time().timeName(),
 919             mesh_,
 920             IOobject::NO_READ,
 921             IOobject::NO_WRITE
 922         ),
 923         mesh_,
 924         dimensionedScalar("zero",  dimensionSet(0, -3, 0, 0, 0), 0.0)
 925     ),
 926     linearKE_
 927     (
 928         IOobject
 929         (
 930             this->name() + "linearKE_",
 931             mesh_.time().timeName(),
 932             mesh_,
 933             IOobject::NO_READ,
 934             IOobject::NO_WRITE
 935         ),
 936         mesh_,
 937         dimensionedScalar("zero",  dimensionSet(1, -1, -2, 0, 0), 0.0)
 938     ),
 939     internalE_
 940     (
 941         IOobject
 942         (
 943             this->name() + "internalE_",
 944             mesh_.time().timeName(),
 945             mesh_,
 946             IOobject::NO_READ,
 947             IOobject::NO_WRITE
 948         ),
 949         mesh_,
 950         dimensionedScalar("zero",  dimensionSet(1, -1, -2, 0, 0), 0.0)
 951     ),
 952     iDof_
 953     (
 954         IOobject
 955         (
 956             this->name() + "iDof_",
 957             mesh_.time().timeName(),
 958             mesh_,
 959             IOobject::NO_READ,
 960             IOobject::NO_WRITE
 961         ),
 962         mesh_,
 963         dimensionedScalar("zero",  dimensionSet(0, -3, 0, 0, 0), VSMALL)
 964     ),
 965     momentum_
 966     (
 967         IOobject
 968         (
 969             this->name() + "momentum_",
 970             mesh_.time().timeName(),
 971             mesh_,
 972             IOobject::NO_READ,
 973             IOobject::NO_WRITE
 974         ),
 975         mesh_,
 976         dimensionedVector
 977         (
 978             "zero",
 979             dimensionSet(1, -2, -1, 0, 0),
 980             vector::zero
 981         )
 982     ),
 983     constProps_(),
 984     rndGen_(label(971501) + 1526*Pstream::myProcNo()),
 985     boundaryT_
 986     (
 987         volScalarField
 988         (
 989             IOobject
 990             (
 991                 "boundaryT",
 992                 mesh_.time().timeName(),
 993                 mesh_,
 994                 IOobject::NO_READ,
 995                 IOobject::NO_WRITE
 996             ),
 997             mesh_,
 998             dimensionedScalar("zero",  dimensionSet(0, 0, 0, 1, 0), 0.0)
 999         )
1000     ),
1001     boundaryU_
1002     (
1003         volVectorField
1004         (
1005             IOobject
1006             (
1007                 "boundaryU",
1008                 mesh_.time().timeName(),
1009                 mesh_,
1010                 IOobject::NO_READ,
1011                 IOobject::NO_WRITE
1012             ),
1013             mesh_,
1014             dimensionedVector
1015             (
1016                 "zero",
1017                 dimensionSet(0, 1, -1, 0, 0),
1018                 vector::zero
1019             )
1020         )
1021     ),
1022     binaryCollisionModel_(),
1023     wallInteractionModel_(),
1024     inflowBoundaryModel_()
1025 {
1026     clear();
1027
1028     buildConstProps();
1029
1030     initialise(dsmcInitialiseDict);
1031 }
1032
1033
1034 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
1035
1036 template<class ParcelType>
1037 Foam::DsmcCloud<ParcelType>::~DsmcCloud()
1038 {}
1039
1040
1041 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
1042
1043 template<class ParcelType>
1044 void Foam::DsmcCloud<ParcelType>::evolve()
1045 {
1046     typename ParcelType::trackData td(*this);
1047
1048     // Reset the data collection fields
1049     resetFields();
1050
1051     if (debug)
1052     {
1053         this->dumpParticlePositions();
1054     }
1055
1056     // Insert new particles from the inflow boundary
1057     this->inflowBoundary().inflow();
1058
1059     // Move the particles ballistically with their current velocities
1060     Cloud<ParcelType>::move(td);
1061
1062     // Calculate new velocities via stochastic collisions
1063     collisions();
1064
1065     // Calculate the volume field data
1066     calculateFields();
1067 }
1068
1069
1070 template<class ParcelType>
1071 void Foam::DsmcCloud<ParcelType>::info() const
1072 {
1073     label nDsmcParticles = this->size();
1074     reduce(nDsmcParticles, sumOp<label>());
1075
1076     scalar nMol = nDsmcParticles*nParticle_;
1077
1078     vector linearMomentum = linearMomentumOfSystem();
1079     reduce(linearMomentum, sumOp<vector>());
1080
1081     scalar linearKineticEnergy = linearKineticEnergyOfSystem();
1082     reduce(linearKineticEnergy, sumOp<scalar>());
1083
1084     scalar internalEnergy = internalEnergyOfSystem();
1085     reduce(internalEnergy, sumOp<scalar>());
1086
1087     Info<< "Cloud name: " << this->name() << nl
1088         << "    Number of dsmc particles        = "
1089         << nDsmcParticles
1090         << endl;
1091
1092     if (nDsmcParticles)
1093     {
1094         Info<< "    Number of molecules             = "
1095             << nMol << nl
1096             << "    Mass in system                  = "
1097             << returnReduce(massInSystem(), sumOp<scalar>()) << nl
1098             << "    Average linear momentum         = "
1099             << linearMomentum/nMol << nl
1100             << "   |Average linear momentum|        = "
1101             << mag(linearMomentum)/nMol << nl
1102             << "    Average linear kinetic energy   = "
1103             << linearKineticEnergy/nMol << nl
1104             << "    Average internal energy         = "
1105             << internalEnergy/nMol << nl
1106             << "    Average total energy            = "
1107             << (internalEnergy + linearKineticEnergy)/nMol
1108             << endl;
1109     }
1110 }
1111
1112
1113 template<class ParcelType>
1114 Foam::vector Foam::DsmcCloud<ParcelType>::equipartitionLinearVelocity
1115 (
1116     scalar temperature,
1117     scalar mass
1118 )
1119 {
1120     return
1121         sqrt(kb*temperature/mass)
1122        *vector
1123         (
1124             rndGen_.GaussNormal(),
1125             rndGen_.GaussNormal(),
1126             rndGen_.GaussNormal()
1127         );
1128 }
1129
1130
1131 template<class ParcelType>
1132 Foam::scalar Foam::DsmcCloud<ParcelType>::equipartitionInternalEnergy
1133 (
1134     scalar temperature,
1135     scalar iDof
1136 )
1137 {
1138     scalar Ei = 0.0;
1139
1140     if (iDof < SMALL)
1141     {
1142         return Ei;
1143     }
1144     else if (iDof < 2.0 + SMALL && iDof > 2.0 - SMALL)
1145     {
1146         // Special case for iDof = 2, i.e. diatomics;
1147         Ei = -log(rndGen_.scalar01())*kb*temperature;
1148     }
1149     else
1150     {
1151         scalar a = 0.5*iDof - 1;
1152
1153         scalar energyRatio;
1154
1155         scalar P = -1;
1156
1157         do
1158         {
1159             energyRatio = 10*rndGen_.scalar01();
1160
1161             P = pow((energyRatio/a), a)*exp(a - energyRatio);
1162
1163         } while (P < rndGen_.scalar01());
1164
1165         Ei = energyRatio*kb*temperature;
1166     }
1167
1168     return Ei;
1169 }
1170
1171
1172 template<class ParcelType>
1173 void Foam::DsmcCloud<ParcelType>::dumpParticlePositions() const
1174 {
1175     OFstream pObj
1176     (
1177         this->db().time().path()/"parcelPositions_"
1178       + this->name() + "_"
1179       + this->db().time().timeName() + ".obj"
1180     );
1181
1182     forAllConstIter(typename DsmcCloud<ParcelType>, *this, iter)
1183     {
1184         const ParcelType& p = iter();
1185
1186         pObj<< "v " << p.position().x()
1187             << " "  << p.position().y()
1188             << " "  << p.position().z()
1189             << nl;
1190     }
1191
1192     pObj.flush();
1193 }
1194
1195
1196 // ************************************************************************* //