34 template<
class ParcelType>
40 template<
class ParcelType>
43 Info<<
nl <<
"Constructing constant properties for" <<
endl;
44 constProps_.setSize(typeIdList_.size());
46 dictionary moleculeProperties
48 particleProperties_.subDict(
"moleculeProperties")
53 const word& id(typeIdList_[i]);
57 const dictionary& molDict(moleculeProperties.subDict(
id));
60 typename ParcelType::constantProperties::constantProperties(molDict);
65 template<
class ParcelType>
70 cellOccupancy_[cO].clear();
73 forAllIter(
typename DsmcCloud<ParcelType>, *
this, iter)
75 cellOccupancy_[iter().cell()].append(&iter());
80 template<
class ParcelType>
83 const IOdictionary& dsmcInitialiseDict
88 const scalar temperature
90 readScalar(dsmcInitialiseDict.lookup(
"temperature"))
93 const vector velocity(dsmcInitialiseDict.lookup(
"velocity"));
95 const dictionary& numberDensitiesDict
97 dsmcInitialiseDict.subDict(
"numberDensities")
100 List<word> molecules(numberDensitiesDict.toc());
102 Field<scalar> numberDensities(molecules.size());
108 numberDensitiesDict.lookup(molecules[i])
112 numberDensities /= nParticle_;
114 forAll(mesh_.cells(), cell)
116 const vector& cC = mesh_.cellCentres()[cell];
117 const labelList& cellFaces = mesh_.cells()[cell];
118 const scalar cV = mesh_.cellVolumes()[cell];
125 nTets += mesh_.faces()[cellFaces[face]].size() - 2;
132 List<labelList> tetPtIs(nTets,
labelList(3,-1));
139 const labelList& facePoints = mesh_.faces()[cellFaces[face]];
142 while ((pointI + 1) < facePoints.size())
145 const vector& pA = mesh_.points()[facePoints[0]];
146 const vector& pB = mesh_.points()[facePoints[pointI]];
147 const vector& pC = mesh_.points()[facePoints[pointI + 1]];
150 mag(((pA - cC) ^ (pB - cC)) & (pC - cC))/(cV*6.0)
151 + cTetVFracs[
max((tet - 1),0)];
153 tetPtIs[tet][0] = facePoints[0];
154 tetPtIs[tet][1] = facePoints[pointI];
155 tetPtIs[tet][2] = facePoints[pointI + 1];
164 cTetVFracs[nTets - 1] = 1.0;
168 const word& moleculeName(molecules[i]);
170 label typeId(
findIndex(typeIdList_, moleculeName));
175 <<
"typeId " << moleculeName <<
"not defined." <<
nl
179 const typename ParcelType::constantProperties& cP =
182 scalar numberDensity = numberDensities[i];
185 scalar particlesRequired = numberDensity*mesh_.cellVolumes()[cell];
188 label nParticlesToInsert = label(particlesRequired);
192 if ((particlesRequired - nParticlesToInsert) > rndGen_.scalar01())
194 nParticlesToInsert++;
197 for (label pI = 0; pI < nParticlesToInsert; pI++)
201 scalar s = rndGen_.scalar01();
202 scalar t = rndGen_.scalar01();
203 scalar u = rndGen_.scalar01();
217 else if (s + t + u > 1.0)
226 scalar tetSelection = rndGen_.scalar01();
235 if (cTetVFracs[tet] >= tetSelection)
243 + s*mesh_.points()[tetPtIs[sTet][0]]
244 + t*mesh_.points()[tetPtIs[sTet][1]]
245 + u*mesh_.points()[tetPtIs[sTet][2]];
247 vector U = equipartitionLinearVelocity
253 scalar Ei = equipartitionInternalEnergy
256 cP.internalDegreesOfFreedom()
277 label mostAbundantType(
findMax(numberDensities));
279 const typename ParcelType::constantProperties& cP = constProps
284 sigmaTcRMax_.internalField() = cP.sigmaT()*maxwellianMostProbableSpeed
290 sigmaTcRMax_.correctBoundaryConditions();
294 template<
class ParcelType>
297 buildCellOccupancy();
300 List<DynamicList<label> > subCells(8);
302 scalar deltaT =
mesh().time().deltaTValue();
304 label collisionCandidates = 0;
306 label collisions = 0;
308 forAll(cellOccupancy_, celli)
310 const DynamicList<ParcelType*>& cellParcels(cellOccupancy_[celli]);
312 label nC(cellParcels.size());
327 List<label> whichSubCell(cellParcels.size());
329 const point& cC = mesh_.cellCentres()[celli];
333 ParcelType* p = cellParcels[i];
335 vector relPos = p->position() - cC;
338 pos(relPos.x()) + 2*
pos(relPos.y()) + 4*
pos(relPos.z());
340 subCells[subCell].append(i);
342 whichSubCell[i] = subCell;
347 scalar sigmaTcRMax = sigmaTcRMax_[celli];
349 scalar selectedPairs =
350 collisionSelectionRemainder_[celli]
351 + 0.5*nC*(nC - 1)*nParticle_*sigmaTcRMax*deltaT
352 /mesh_.cellVolumes()[celli];
354 label nCandidates(selectedPairs);
356 collisionSelectionRemainder_[celli] = selectedPairs - nCandidates;
358 collisionCandidates += nCandidates;
360 for (label c = 0; c < nCandidates; c++)
366 label candidateP = rndGen_.integer(0, nC - 1);
369 label candidateQ = -1;
371 List<label> subCellPs = subCells[whichSubCell[candidateP]];
373 label nSC = subCellPs.size();
383 candidateQ = subCellPs[rndGen_.integer(0, nSC - 1)];
385 }
while(candidateP == candidateQ);
395 candidateQ = rndGen_.integer(0, nC - 1);
397 }
while(candidateP == candidateQ);
417 ParcelType* parcelP = cellParcels[candidateP];
418 ParcelType* parcelQ = cellParcels[candidateQ];
420 label typeIdP = parcelP->typeId();
421 label typeIdQ = parcelQ->typeId();
423 scalar sigmaTcR = binaryCollision().sigmaTcR
435 if (sigmaTcR > sigmaTcRMax_[celli])
437 sigmaTcRMax_[celli] = sigmaTcR;
440 if ((sigmaTcR/sigmaTcRMax) > rndGen_.scalar01())
442 binaryCollision().collide
458 reduce(collisions, sumOp<label>());
460 reduce(collisionCandidates, sumOp<label>());
462 sigmaTcRMax_.correctBoundaryConditions();
464 if (collisionCandidates)
466 Info<<
" Collisions = "
468 <<
" Acceptance rate = "
469 << scalar(collisions)/scalar(collisionCandidates) <<
nl
479 template<
class ParcelType>
487 dimensionSet(1, -1, -2, 0, 0),
506 dimensionSet(1, -2, -1, 0, 0),
512 template<
class ParcelType>
523 scalarField& internalE = internalE_.internalField();
531 const ParcelType& p = iter();
532 const label cellI = p.cell();
536 rhoM[cellI] += constProps(p.typeId()).mass();
540 linearKE[cellI] += 0.5*constProps(p.typeId()).mass()*(p.U() & p.U());
542 internalE[cellI] += p.Ei();
544 iDof[cellI] += constProps(p.typeId()).internalDegreesOfFreedom();
546 momentum[cellI] += constProps(p.typeId()).mass()*p.U();
549 rhoN *= nParticle_/
mesh().cellVolumes();
550 rhoN_.correctBoundaryConditions();
552 rhoM *= nParticle_/
mesh().cellVolumes();
553 rhoM_.correctBoundaryConditions();
555 dsmcRhoN_.correctBoundaryConditions();
557 linearKE *= nParticle_/
mesh().cellVolumes();
558 linearKE_.correctBoundaryConditions();
560 internalE *= nParticle_/
mesh().cellVolumes();
561 internalE_.correctBoundaryConditions();
563 iDof *= nParticle_/
mesh().cellVolumes();
564 iDof_.correctBoundaryConditions();
566 momentum *= nParticle_/
mesh().cellVolumes();
567 momentum_.correctBoundaryConditions();
573 template<
class ParcelType>
583 ParcelType* pPtr =
new ParcelType
599 template<
class ParcelType>
609 cloudName_(cloudName),
615 cloudName +
"Properties",
622 typeIdList_(particleProperties_.lookup(
"typeIdList")),
623 nParticle_(
readScalar(particleProperties_.lookup(
"nEquivalentParticles"))),
624 cellOccupancy_(mesh_.nCells()),
629 this->
name() +
"SigmaTcRMax",
630 mesh_.time().timeName(),
637 collisionSelectionRemainder_(mesh_.nCells(), 0),
643 mesh_.time().timeName(),
655 mesh_.time().timeName(),
667 mesh_.time().timeName(),
679 mesh_.time().timeName(),
691 mesh_.time().timeName(),
703 mesh_.time().timeName(),
715 mesh_.time().timeName(),
727 mesh_.time().timeName(),
739 mesh_.time().timeName(),
747 rndGen_(label(149382906) + 7183*Pstream::myProcNo()),
755 mesh_.time().timeName(),
770 mesh_.time().timeName(),
778 binaryCollisionModel_
786 wallInteractionModel_
805 buildCellOccupancy();
809 forAll(collisionSelectionRemainder_, i)
811 collisionSelectionRemainder_[i] = rndGen_.scalar01();
821 template<
class ParcelType>
824 const word& cloudName,
831 cloudName_(cloudName),
837 cloudName +
"Properties",
844 typeIdList_(particleProperties_.lookup(
"typeIdList")),
845 nParticle_(
readScalar(particleProperties_.lookup(
"nEquivalentParticles"))),
851 this->
name() +
"SigmaTcRMax",
852 mesh_.time().timeName(),
859 zeroGradientFvPatchScalarField::typeName
861 collisionSelectionRemainder_(),
867 mesh_.time().timeName(),
879 this->
name() +
"fD_",
880 mesh_.time().timeName(),
897 this->
name() +
"rhoN_",
898 mesh_.time().timeName(),
910 this->
name() +
"rhoM_",
911 mesh_.time().timeName(),
923 this->
name() +
"dsmcRhoN_",
924 mesh_.time().timeName(),
936 this->
name() +
"linearKE_",
937 mesh_.time().timeName(),
949 this->
name() +
"internalE_",
950 mesh_.time().timeName(),
962 this->
name() +
"iDof_",
963 mesh_.time().timeName(),
975 this->
name() +
"momentum_",
976 mesh_.time().timeName(),
990 rndGen_(label(971501) + 1526*Pstream::myProcNo()),
998 mesh_.time().timeName(),
1014 mesh_.time().timeName(),
1028 binaryCollisionModel_(),
1029 wallInteractionModel_(),
1030 inflowBoundaryModel_()
1036 initialise(dsmcInitialiseDict);
1042 template<
class ParcelType>
1049 template<
class ParcelType>
1052 typename ParcelType::trackData td(*
this);
1059 this->dumpParticlePositions();
1063 this->inflowBoundary().inflow();
1076 template<
class ParcelType>
1079 label nDsmcParticles = this->size();
1082 scalar nMol = nDsmcParticles*nParticle_;
1084 vector linearMomentum = linearMomentumOfSystem();
1087 scalar linearKineticEnergy = linearKineticEnergyOfSystem();
1090 scalar internalEnergy = internalEnergyOfSystem();
1094 <<
" Number of dsmc particles = "
1100 Info<<
" Number of molecules = "
1102 <<
" Mass in system = "
1104 <<
" Average linear momentum = "
1105 << linearMomentum/nMol <<
nl
1106 <<
" |Average linear momentum| = "
1107 <<
mag(linearMomentum)/nMol <<
nl
1108 <<
" Average linear kinetic energy = "
1109 << linearKineticEnergy/nMol <<
nl
1110 <<
" Average internal energy = "
1111 << internalEnergy/nMol <<
nl
1112 <<
" Average total energy = "
1113 << (internalEnergy + linearKineticEnergy)/nMol
1119 template<
class ParcelType>
1127 sqrt(kb*temperature/mass)
1130 rndGen_.GaussNormal(),
1131 rndGen_.GaussNormal(),
1132 rndGen_.GaussNormal()
1137 template<
class ParcelType>
1150 else if (iDof < 2.0 + SMALL && iDof > 2.0 - SMALL)
1153 Ei = -
log(rndGen_.scalar01())*kb*temperature;
1157 scalar a = 0.5*iDof - 1;
1165 energyRatio = 10*rndGen_.scalar01();
1167 P =
pow((energyRatio/a), a)*
exp(a - energyRatio);
1169 }
while (P < rndGen_.scalar01());
1171 Ei = energyRatio*kb*temperature;
1178 template<
class ParcelType>
1183 this->db().time().path()/
"parcelPositions_"
1184 + this->
name() +
"_"
1185 + this->db().time().
timeName() +
".obj"
1190 const ParcelType& p = iter();
1192 pObj<<
"v " << p.position().x()
1193 <<
" " << p.position().y()
1194 <<
" " << p.position().z()