38 template<
class ParcelType>
41 Info<<
nl <<
"Constructing constant properties for" <<
endl;
42 constProps_.setSize(typeIdList_.size());
44 dictionary moleculeProperties
46 particleProperties_.subDict(
"moleculeProperties")
51 const word& id(typeIdList_[i]);
55 const dictionary& molDict(moleculeProperties.subDict(
id));
58 typename ParcelType::constantProperties::constantProperties(molDict);
63 template<
class ParcelType>
68 cellOccupancy_[cO].clear();
71 forAllIter(
typename DSMCCloud<ParcelType>, *
this, iter)
73 cellOccupancy_[iter().cell()].append(&iter());
78 template<
class ParcelType>
81 const IOdictionary& dsmcInitialiseDict
86 const scalar temperature
88 readScalar(dsmcInitialiseDict.lookup(
"temperature"))
91 const vector velocity(dsmcInitialiseDict.lookup(
"velocity"));
93 const dictionary& numberDensitiesDict
95 dsmcInitialiseDict.subDict(
"numberDensities")
98 List<word> molecules(numberDensitiesDict.toc());
100 Field<scalar> numberDensities(molecules.size());
106 numberDensitiesDict.lookup(molecules[i])
110 numberDensities /= nParticle_;
112 forAll(mesh_.cells(), celli)
122 const tetIndices& cellTetIs = cellTets[tetI];
124 scalar tetVolume = tet.mag();
128 const word& moleculeName(molecules[i]);
135 <<
"typeId " << moleculeName <<
"not defined." <<
nl 139 const typename ParcelType::constantProperties& cP =
142 scalar numberDensity = numberDensities[i];
145 scalar particlesRequired = numberDensity*tetVolume;
148 label nParticlesToInsert =
label(particlesRequired);
154 (particlesRequired - nParticlesToInsert)
158 nParticlesToInsert++;
161 for (
label pI = 0; pI < nParticlesToInsert; pI++)
163 point p = tet.randomPoint(rndGen_);
165 vector U = equipartitionLinearVelocity
171 scalar Ei = equipartitionInternalEnergy
174 cP.internalDegreesOfFreedom()
200 const typename ParcelType::constantProperties& cP = constProps
205 sigmaTcRMax_.primitiveFieldRef() = cP.sigmaT()*maxwellianMostProbableSpeed
211 sigmaTcRMax_.correctBoundaryConditions();
215 template<
class ParcelType>
218 if (!binaryCollision().active())
224 List<DynamicList<label>> subCells(8);
226 scalar deltaT =
mesh().time().deltaTValue();
228 label collisionCandidates = 0;
230 label collisions = 0;
232 forAll(cellOccupancy_, celli)
234 const DynamicList<ParcelType*>& cellParcels(cellOccupancy_[celli]);
236 label nC(cellParcels.size());
250 List<label> whichSubCell(cellParcels.size());
252 const point& cC = mesh_.cellCentres()[celli];
256 const ParcelType& p = *cellParcels[i];
257 vector relPos = p.position() - cC;
260 pos(relPos.x()) + 2*
pos(relPos.y()) + 4*
pos(relPos.z());
262 subCells[subCell].append(i);
263 whichSubCell[i] = subCell;
268 scalar sigmaTcRMax = sigmaTcRMax_[celli];
270 scalar selectedPairs =
271 collisionSelectionRemainder_[celli]
272 + 0.5*nC*(nC - 1)*nParticle_*sigmaTcRMax*deltaT
273 /mesh_.cellVolumes()[celli];
275 label nCandidates(selectedPairs);
276 collisionSelectionRemainder_[celli] = selectedPairs - nCandidates;
277 collisionCandidates += nCandidates;
279 for (
label c = 0;
c < nCandidates;
c++)
285 label candidateP = rndGen_.integer(0, nC - 1);
288 label candidateQ = -1;
290 List<label> subCellPs = subCells[whichSubCell[candidateP]];
291 label nSC = subCellPs.size();
301 candidateQ = subCellPs[rndGen_.integer(0, nSC - 1)];
302 }
while (candidateP == candidateQ);
312 candidateQ = rndGen_.integer(0, nC - 1);
313 }
while (candidateP == candidateQ);
333 ParcelType& parcelP = *cellParcels[candidateP];
334 ParcelType& parcelQ = *cellParcels[candidateQ];
336 scalar sigmaTcR = binaryCollision().sigmaTcR
346 if (sigmaTcR > sigmaTcRMax_[celli])
348 sigmaTcRMax_[celli] = sigmaTcR;
351 if ((sigmaTcR/sigmaTcRMax) > rndGen_.scalar01())
353 binaryCollision().collide
365 reduce(collisions, sumOp<label>());
367 reduce(collisionCandidates, sumOp<label>());
369 sigmaTcRMax_.correctBoundaryConditions();
371 if (collisionCandidates)
373 Info<<
" Collisions = " 375 <<
" Acceptance rate = " 376 << scalar(collisions)/scalar(collisionCandidates) <<
nl 386 template<
class ParcelType>
394 dimensionSet(1, -1, -2, 0, 0),
408 dimensionSet(1, -2, -1, 0, 0),
414 template<
class ParcelType>
419 scalarField& dsmcRhoN = dsmcRhoN_.primitiveFieldRef();
420 scalarField& linearKE = linearKE_.primitiveFieldRef();
421 scalarField& internalE = internalE_.primitiveFieldRef();
423 vectorField& momentum = momentum_.primitiveFieldRef();
427 const ParcelType& p = iter();
428 const label celli = p.cell();
431 rhoM[celli] += constProps(p.typeId()).mass();
433 linearKE[celli] += 0.5*constProps(p.typeId()).mass()*(p.U() & p.U());
434 internalE[celli] += p.Ei();
435 iDof[celli] += constProps(p.typeId()).internalDegreesOfFreedom();
436 momentum[celli] += constProps(p.typeId()).mass()*p.U();
439 rhoN *= nParticle_/
mesh().cellVolumes();
440 rhoN_.correctBoundaryConditions();
442 rhoM *= nParticle_/
mesh().cellVolumes();
443 rhoM_.correctBoundaryConditions();
445 dsmcRhoN_.correctBoundaryConditions();
447 linearKE *= nParticle_/
mesh().cellVolumes();
448 linearKE_.correctBoundaryConditions();
450 internalE *= nParticle_/
mesh().cellVolumes();
451 internalE_.correctBoundaryConditions();
453 iDof *= nParticle_/
mesh().cellVolumes();
454 iDof_.correctBoundaryConditions();
456 momentum *= nParticle_/
mesh().cellVolumes();
457 momentum_.correctBoundaryConditions();
463 template<
class ParcelType>
470 const label tetFacei,
475 ParcelType* pPtr =
new ParcelType
487 this->addParticle(pPtr);
493 template<
class ParcelType>
496 const word& cloudName,
503 cloudName_(cloudName),
509 cloudName +
"Properties",
516 typeIdList_(particleProperties_.lookup(
"typeIdList")),
517 nParticle_(
readScalar(particleProperties_.lookup(
"nEquivalentParticles"))),
518 cellOccupancy_(mesh_.nCells()),
523 this->
name() +
"SigmaTcRMax",
524 mesh_.time().timeName(),
531 collisionSelectionRemainder_
535 this->
name() +
":collisionSelectionRemainder",
536 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(),
619 mesh_.time().timeName(),
631 mesh_.time().timeName(),
643 mesh_.time().timeName(),
659 mesh_.time().timeName(),
674 mesh_.time().timeName(),
682 binaryCollisionModel_
690 wallInteractionModel_
708 buildCellOccupancy();
712 forAll(collisionSelectionRemainder_, i)
714 collisionSelectionRemainder_[i] = rndGen_.scalar01();
724 template<
class ParcelType>
727 const word& cloudName,
734 cloudName_(cloudName),
740 cloudName +
"Properties",
747 typeIdList_(particleProperties_.lookup(
"typeIdList")),
748 nParticle_(
readScalar(particleProperties_.lookup(
"nEquivalentParticles"))),
754 this->
name() +
"SigmaTcRMax",
755 mesh_.time().timeName(),
762 zeroGradientFvPatchScalarField::typeName
764 collisionSelectionRemainder_
768 this->
name() +
":collisionSelectionRemainder",
769 mesh_.time().timeName(),
780 mesh_.time().timeName(),
792 this->
name() +
"fD_",
793 mesh_.time().timeName(),
810 this->
name() +
"rhoN_",
811 mesh_.time().timeName(),
823 this->
name() +
"rhoM_",
824 mesh_.time().timeName(),
836 this->
name() +
"dsmcRhoN_",
837 mesh_.time().timeName(),
849 this->
name() +
"linearKE_",
850 mesh_.time().timeName(),
862 this->
name() +
"internalE_",
863 mesh_.time().timeName(),
875 this->
name() +
"iDof_",
876 mesh_.time().timeName(),
888 this->
name() +
"momentum_",
889 mesh_.time().timeName(),
911 mesh_.time().timeName(),
927 mesh_.time().timeName(),
941 binaryCollisionModel_(),
942 wallInteractionModel_(),
943 inflowBoundaryModel_()
947 initialise(dsmcInitialiseDict);
953 template<
class ParcelType>
960 template<
class ParcelType>
963 typename ParcelType::trackingData td(*
this);
970 this->dumpParticlePositions();
974 this->inflowBoundary().inflow();
980 buildCellOccupancy();
990 template<
class ParcelType>
993 label nDSMCParticles = this->size();
996 scalar nMol = nDSMCParticles*nParticle_;
998 vector linearMomentum = linearMomentumOfSystem();
1001 scalar linearKineticEnergy = linearKineticEnergyOfSystem();
1004 scalar internalEnergy = internalEnergyOfSystem();
1008 <<
" Number of dsmc particles = " 1014 Info<<
" Number of molecules = " 1016 <<
" Mass in system = " 1018 <<
" Average linear momentum = " 1019 << linearMomentum/nMol <<
nl 1020 <<
" |Average linear momentum| = " 1021 <<
mag(linearMomentum)/nMol <<
nl 1022 <<
" Average linear kinetic energy = " 1023 << linearKineticEnergy/nMol <<
nl 1024 <<
" Average internal energy = " 1025 << internalEnergy/nMol <<
nl 1026 <<
" Average total energy = " 1027 << (internalEnergy + linearKineticEnergy)/nMol
1033 template<
class ParcelType>
1044 rndGen_.GaussNormal(),
1045 rndGen_.GaussNormal(),
1046 rndGen_.GaussNormal()
1051 template<
class ParcelType>
1064 else if (iDof < 2.0 + SMALL && iDof > 2.0 - SMALL)
1071 scalar a = 0.5*iDof - 1;
1077 energyRatio = 10*rndGen_.scalar01();
1078 P =
pow((energyRatio/a), a)*
exp(a - energyRatio);
1079 }
while (P < rndGen_.scalar01());
1088 template<
class ParcelType>
1093 this->db().time().
path()/
"parcelPositions_" 1094 + this->
name() +
"_" 1095 + this->db().time().
timeName() +
".obj" 1100 const ParcelType& p = iter();
1102 pObj<<
"v " << p.position().x()
1103 <<
" " << p.position().y()
1104 <<
" " << p.position().z()
1112 template<
class ParcelType>
1115 typedef typename ParcelType::trackingData tdType;
1120 cellOccupancy_.setSize(mesh_.nCells());
1121 buildCellOccupancy();
1124 this->inflowBoundary().autoMap(mapper);
const Time & time() const
Return time.
void addNewParcel(const vector &position, const vector &U, const scalar Ei, const label celli, const label tetFacei, const label tetPtI, const label typeId)
Add new parcel.
virtual void autoMap(const mapPolyMesh &)
Remap the particles to the correct cells following mesh change.
#define forAll(list, i)
Loop across all elements in list.
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
dimensionedScalar log(const dimensionedScalar &ds)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
tetrahedron< point, const point & > tetPointRef
Templated inflow boundary model class.
#define forAllIter(Container, container, iter)
Iterate across all elements in the container object of type.
static int myProcNo(const label communicator=0)
Number of this process (starting from masterNo() = 0)
Templated DSMC particle collision class.
dimensioned< vector > dimensionedVector
Dimensioned vector obtained from generic dimensioned type.
dimensionedScalar sqrt(const dimensionedScalar &ds)
Ostream & endl(Ostream &os)
Add newline and flush stream.
const word cloudName(propsDict.lookup("cloudName"))
const Type & value() const
Return const reference to value.
void dumpParticlePositions() const
Dump particle positions to .obj file.
Vector< scalar > vector
A scalar version of the templated Vector.
Class containing mesh-to-mesh mapping information after a change in polyMesh topology.
void readFields(const typename GeoFieldType::Mesh &mesh, const IOobjectList &objects, const HashSet< word > &selectedFields, LIFOStack< regIOobject * > &storedObjects)
Read the selected GeometricFields of the specified type.
static List< tetIndices > cellTetIndices(const polyMesh &mesh, label cI)
Return the tet decomposition of the given cell, see.
tmp< DimensionedField< TypeR, GeoMesh > > New(const tmp< DimensionedField< TypeR, GeoMesh >> &tdf1, const word &name, const dimensionSet &dimensions)
IOdictionary is derived from dictionary and IOobject to give the dictionary automatic IO functionalit...
dimensionedScalar pos(const dimensionedScalar &ds)
Dimension set for the base types.
dimensionedScalar exp(const dimensionedScalar &ds)
A class for handling words, derived from string.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
virtual void flush()
Flush stream.
const word & constant() const
Return constant name.
bool readScalar(const char *buf, doubleScalar &s)
Read whole of buf as a scalar. Return true if succesful.
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
errorManip< error > abort(error &err)
scalar equipartitionInternalEnergy(scalar temperature, direction internalDegreesOfFreedom)
Generate a random internal energy, sampled from the.
void move(TrackData &td, const scalar trackTime)
Move the particles.
label findMax(const ListType &, const label start=0)
Find index of max element (and larger than given element).
fileName path(UMean.rootPath()/UMean.caseName()/"graphs"/UMean.instance())
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurence of given element and return index,.
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
word name(const complex &)
Return a string representation of a complex.
virtual ~DSMCCloud()
Destructor.
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
const dimensionedScalar k
Boltzmann constant.
const dimensionSet dimless(0, 0, 0, 0, 0, 0, 0)
vector point
Point is a vector.
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
Mesh data needed to do the Finite Volume discretisation.
const dimensionedScalar c
Speed of light in a vacuum.
void evolve()
Evolve the cloud (move, collide)
dimensioned< scalar > mag(const dimensioned< Type > &)
Field< vector > vectorField
Specialisation of Field<T> for vector.
Virtual abstract base class for templated DSMCCloud.
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
Templated wall interaction model class.
vector equipartitionLinearVelocity(scalar temperature, scalar mass)
Generate a random velocity sampled from the Maxwellian speed.
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Templated base class for dsmc cloud.
void info() const
Print cloud information.