41 template<
class ParcelType>
44 Info<<
nl <<
"Constructing constant properties for" <<
endl;
45 constProps_.setSize(typeIdList_.size());
47 dictionary moleculeProperties
49 particleProperties_.subDict(
"moleculeProperties")
54 const word&
id = typeIdList_[i];
58 const dictionary& molDict(moleculeProperties.subDict(
id));
60 constProps_[i] =
typename ParcelType::constantProperties(molDict);
65 template<
class ParcelType>
68 for (
auto& list : cellOccupancy_)
73 for (ParcelType&
p : *
this)
75 cellOccupancy_[
p.cell()].append(&
p);
80 template<
class ParcelType>
83 const IOdictionary& dsmcInitialiseDict
88 const scalar temperature
90 dsmcInitialiseDict.get<scalar>(
"temperature")
93 const vector velocity(dsmcInitialiseDict.get<
vector>(
"velocity"));
95 const dictionary& numberDensitiesDict
97 dsmcInitialiseDict.subDict(
"numberDensities")
100 List<word> molecules(numberDensitiesDict.toc());
102 Field<scalar> numberDensities(molecules.size());
106 numberDensities[i] = numberDensitiesDict.get<scalar>(molecules[i]);
109 numberDensities /= nParticle_;
111 forAll(mesh_.cells(), celli)
121 const tetIndices& cellTetIs = cellTets[tetI];
123 scalar tetVolume = tet.mag();
127 const word& moleculeName(molecules[i]);
129 label typeId = typeIdList_.find(moleculeName);
134 <<
"typeId " << moleculeName <<
"not defined." <<
nl
138 const typename ParcelType::constantProperties& cP =
141 scalar numberDensity = numberDensities[i];
144 scalar particlesRequired = numberDensity*tetVolume;
147 label nParticlesToInsert = label(particlesRequired);
153 (particlesRequired - nParticlesToInsert)
154 > rndGen_.sample01<scalar>()
157 nParticlesToInsert++;
160 for (label pI = 0; pI < nParticlesToInsert; pI++)
162 point p = tet.randomPoint(rndGen_);
164 vector U = equipartitionLinearVelocity
170 scalar Ei = equipartitionInternalEnergy
173 cP.internalDegreesOfFreedom()
178 addNewParcel(
p, celli,
U, Ei, typeId);
188 label mostAbundantType(
findMax(numberDensities));
190 const typename ParcelType::constantProperties& cP = constProps
195 sigmaTcRMax_.primitiveFieldRef() = cP.sigmaT()*maxwellianMostProbableSpeed
201 sigmaTcRMax_.correctBoundaryConditions();
205 template<
class ParcelType>
208 if (!binaryCollision().active())
214 List<DynamicList<label>> subCells(8);
216 scalar deltaT =
mesh().time().deltaTValue();
218 label collisionCandidates = 0;
220 label collisions = 0;
222 forAll(cellOccupancy_, celli)
224 const DynamicList<ParcelType*>& cellParcels(cellOccupancy_[celli]);
226 label nC(cellParcels.size());
240 List<label> whichSubCell(cellParcels.size());
242 const point& cC = mesh_.cellCentres()[celli];
246 const ParcelType&
p = *cellParcels[i];
247 vector relPos =
p.position() - cC;
250 pos0(relPos.x()) + 2*
pos0(relPos.y()) + 4*
pos0(relPos.z());
252 subCells[subCell].append(i);
253 whichSubCell[i] = subCell;
258 scalar sigmaTcRMax = sigmaTcRMax_[celli];
260 scalar selectedPairs =
261 collisionSelectionRemainder_[celli]
262 + 0.5*nC*(nC - 1)*nParticle_*sigmaTcRMax*deltaT
263 /mesh_.cellVolumes()[celli];
265 label nCandidates(selectedPairs);
266 collisionSelectionRemainder_[celli] = selectedPairs - nCandidates;
267 collisionCandidates += nCandidates;
269 for (label
c = 0;
c < nCandidates;
c++)
275 label candidateP = rndGen_.position<label>(0, nC - 1);
278 label candidateQ = -1;
280 List<label> subCellPs = subCells[whichSubCell[candidateP]];
281 label nSC = subCellPs.size();
291 label i = rndGen_.position<label>(0, nSC - 1);
292 candidateQ = subCellPs[i];
293 }
while (candidateP == candidateQ);
303 candidateQ = rndGen_.position<label>(0, nC - 1);
304 }
while (candidateP == candidateQ);
324 ParcelType& parcelP = *cellParcels[candidateP];
325 ParcelType& parcelQ = *cellParcels[candidateQ];
327 scalar sigmaTcR = binaryCollision().sigmaTcR
337 if (sigmaTcR > sigmaTcRMax_[celli])
339 sigmaTcRMax_[celli] = sigmaTcR;
342 if ((sigmaTcR/sigmaTcRMax) > rndGen_.sample01<scalar>())
344 binaryCollision().collide
356 reduce(collisions, sumOp<label>());
358 reduce(collisionCandidates, sumOp<label>());
360 sigmaTcRMax_.correctBoundaryConditions();
362 if (collisionCandidates)
364 Info<<
" Collisions = "
366 <<
" Acceptance rate = "
367 << scalar(collisions)/scalar(collisionCandidates) <<
nl
377 template<
class ParcelType>
393 template<
class ParcelType>
398 scalarField& dsmcRhoN = dsmcRhoN_.primitiveFieldRef();
399 scalarField& linearKE = linearKE_.primitiveFieldRef();
400 scalarField& internalE = internalE_.primitiveFieldRef();
402 vectorField& momentum = momentum_.primitiveFieldRef();
404 for (
const ParcelType&
p : *
this)
406 const label celli =
p.cell();
409 rhoM[celli] += constProps(
p.typeId()).mass();
411 linearKE[celli] += 0.5*constProps(
p.typeId()).mass()*(
p.U() &
p.U());
412 internalE[celli] +=
p.Ei();
413 iDof[celli] += constProps(
p.typeId()).internalDegreesOfFreedom();
414 momentum[celli] += constProps(
p.typeId()).mass()*
p.U();
417 rhoN *= nParticle_/
mesh().cellVolumes();
418 rhoN_.correctBoundaryConditions();
420 rhoM *= nParticle_/
mesh().cellVolumes();
421 rhoM_.correctBoundaryConditions();
423 dsmcRhoN_.correctBoundaryConditions();
425 linearKE *= nParticle_/
mesh().cellVolumes();
426 linearKE_.correctBoundaryConditions();
428 internalE *= nParticle_/
mesh().cellVolumes();
429 internalE_.correctBoundaryConditions();
431 iDof *= nParticle_/
mesh().cellVolumes();
432 iDof_.correctBoundaryConditions();
434 momentum *= nParticle_/
mesh().cellVolumes();
435 momentum_.correctBoundaryConditions();
441 template<
class ParcelType>
451 this->addParticle(
new ParcelType(mesh_, position, celli,
U, Ei, typeId));
457 template<
class ParcelType>
474 mesh_.time().constant(),
480 typeIdList_(particleProperties_.lookup(
"typeIdList")),
481 nParticle_(particleProperties_.get<scalar>(
"nEquivalentParticles")),
482 cellOccupancy_(mesh_.nCells()),
487 this->
name() +
"SigmaTcRMax",
488 mesh_.time().timeName(),
495 collisionSelectionRemainder_
499 this->
name() +
":collisionSelectionRemainder",
500 mesh_.time().timeName(),
511 mesh_.time().timeName(),
523 mesh_.time().timeName(),
535 mesh_.time().timeName(),
547 mesh_.time().timeName(),
559 mesh_.time().timeName(),
571 mesh_.time().timeName(),
583 mesh_.time().timeName(),
595 mesh_.time().timeName(),
607 mesh_.time().timeName(),
623 mesh_.time().timeName(),
638 mesh_.time().timeName(),
646 binaryCollisionModel_
654 wallInteractionModel_
672 buildCellOccupancy();
676 forAll(collisionSelectionRemainder_, i)
678 collisionSelectionRemainder_[i] = rndGen_.sample01<scalar>();
688 template<
class ParcelType>
705 mesh_.time().constant(),
711 typeIdList_(particleProperties_.lookup(
"typeIdList")),
712 nParticle_(particleProperties_.get<scalar>(
"nEquivalentParticles")),
718 this->
name() +
"SigmaTcRMax",
719 mesh_.time().timeName(),
726 zeroGradientFvPatchScalarField::typeName
728 collisionSelectionRemainder_
732 this->
name() +
":collisionSelectionRemainder",
733 mesh_.time().timeName(),
744 mesh_.time().timeName(),
756 this->
name() +
"fD_",
757 mesh_.time().timeName(),
769 this->
name() +
"rhoN_",
770 mesh_.time().timeName(),
782 this->
name() +
"rhoM_",
783 mesh_.time().timeName(),
795 this->
name() +
"dsmcRhoN_",
796 mesh_.time().timeName(),
808 this->
name() +
"linearKE_",
809 mesh_.time().timeName(),
821 this->
name() +
"internalE_",
822 mesh_.time().timeName(),
834 this->
name() +
"iDof_",
835 mesh_.time().timeName(),
847 this->
name() +
"momentum_",
848 mesh_.time().timeName(),
865 mesh_.time().timeName(),
881 mesh_.time().timeName(),
890 binaryCollisionModel_(),
891 wallInteractionModel_(),
892 inflowBoundaryModel_()
896 initialise(dsmcInitialiseDict);
902 template<
class ParcelType>
909 template<
class ParcelType>
912 typename ParcelType::trackingData td(*
this);
919 this->dumpParticlePositions();
923 this->inflowBoundary().inflow();
929 buildCellOccupancy();
939 template<
class ParcelType>
942 label nDSMCParticles = this->size();
945 scalar nMol = nDSMCParticles*nParticle_;
947 vector linearMomentum = linearMomentumOfSystem();
950 scalar linearKineticEnergy = linearKineticEnergyOfSystem();
953 scalar internalEnergy = internalEnergyOfSystem();
957 <<
" Number of dsmc particles = "
963 Info<<
" Number of molecules = "
965 <<
" Mass in system = "
967 <<
" Average linear momentum = "
968 << linearMomentum/nMol <<
nl
969 <<
" |Average linear momentum| = "
970 <<
mag(linearMomentum)/nMol <<
nl
971 <<
" Average linear kinetic energy = "
972 << linearKineticEnergy/nMol <<
nl
973 <<
" Average internal energy = "
974 << internalEnergy/nMol <<
nl
975 <<
" Average total energy = "
976 << (internalEnergy + linearKineticEnergy)/nMol
982 template<
class ParcelType>
991 *rndGen_.GaussNormal<
vector>();
995 template<
class ParcelType>
1011 -
log(rndGen_.sample01<scalar>())
1017 const scalar a = 0.5*iDof - 1;
1018 scalar energyRatio = 0;
1023 energyRatio = 10*rndGen_.sample01<scalar>();
1024 P =
pow((energyRatio/a), a)*
exp(a - energyRatio);
1025 }
while (P < rndGen_.sample01<scalar>());
1031 template<
class ParcelType>
1036 this->db().time().
path()/
"parcelPositions_"
1037 + this->
name() +
"_"
1038 + this->db().time().
timeName() +
".obj"
1041 for (
const ParcelType&
p : *
this)
1043 pObj<<
"v " <<
p.position().x()
1044 <<
' ' <<
p.position().y()
1045 <<
' ' <<
p.position().z()
1053 template<
class ParcelType>
1059 cellOccupancy_.setSize(mesh_.nCells());
1060 buildCellOccupancy();
1063 this->inflowBoundary().autoMap(mapper);