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()
179 addNewParcel(p, celli, U, Ei, typeId);
191 const typename ParcelType::constantProperties& cP = constProps
196 sigmaTcRMax_.primitiveFieldRef() = cP.sigmaT()*maxwellianMostProbableSpeed
202 sigmaTcRMax_.correctBoundaryConditions();
206 template<
class ParcelType>
209 if (!binaryCollision().active())
215 List<DynamicList<label>> subCells(8);
217 scalar deltaT =
mesh().time().deltaTValue();
219 label collisionCandidates = 0;
221 label collisions = 0;
223 forAll(cellOccupancy_, celli)
225 const DynamicList<ParcelType*>& cellParcels(cellOccupancy_[celli]);
227 label nC(cellParcels.size());
241 List<label> whichSubCell(cellParcels.size());
243 const point& cC = mesh_.cellCentres()[celli];
247 const ParcelType& p = *cellParcels[i];
248 vector relPos = p.position() - cC;
251 pos0(relPos.x()) + 2*
pos0(relPos.y()) + 4*
pos0(relPos.z());
253 subCells[subCell].append(i);
254 whichSubCell[i] = subCell;
259 scalar sigmaTcRMax = sigmaTcRMax_[celli];
261 scalar selectedPairs =
262 collisionSelectionRemainder_[celli]
263 + 0.5*nC*(nC - 1)*nParticle_*sigmaTcRMax*deltaT
264 /mesh_.cellVolumes()[celli];
266 label nCandidates(selectedPairs);
267 collisionSelectionRemainder_[celli] = selectedPairs - nCandidates;
268 collisionCandidates += nCandidates;
270 for (
label c = 0;
c < nCandidates;
c++)
276 label candidateP = rndGen_.integer(0, nC - 1);
279 label candidateQ = -1;
281 List<label> subCellPs = subCells[whichSubCell[candidateP]];
282 label nSC = subCellPs.size();
292 candidateQ = subCellPs[rndGen_.integer(0, nSC - 1)];
293 }
while (candidateP == candidateQ);
303 candidateQ = rndGen_.integer(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_.scalar01())
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>
385 dimensionSet(1, -1, -2, 0, 0),
399 dimensionSet(1, -2, -1, 0, 0),
405 template<
class ParcelType>
410 scalarField& dsmcRhoN = dsmcRhoN_.primitiveFieldRef();
411 scalarField& linearKE = linearKE_.primitiveFieldRef();
412 scalarField& internalE = internalE_.primitiveFieldRef();
414 vectorField& momentum = momentum_.primitiveFieldRef();
418 const ParcelType& p = iter();
419 const label celli = p.cell();
422 rhoM[celli] += constProps(p.typeId()).mass();
424 linearKE[celli] += 0.5*constProps(p.typeId()).mass()*(p.U() & p.U());
425 internalE[celli] += p.Ei();
426 iDof[celli] += constProps(p.typeId()).internalDegreesOfFreedom();
427 momentum[celli] += constProps(p.typeId()).mass()*p.U();
430 rhoN *= nParticle_/
mesh().cellVolumes();
431 rhoN_.correctBoundaryConditions();
433 rhoM *= nParticle_/
mesh().cellVolumes();
434 rhoM_.correctBoundaryConditions();
436 dsmcRhoN_.correctBoundaryConditions();
438 linearKE *= nParticle_/
mesh().cellVolumes();
439 linearKE_.correctBoundaryConditions();
441 internalE *= nParticle_/
mesh().cellVolumes();
442 internalE_.correctBoundaryConditions();
444 iDof *= nParticle_/
mesh().cellVolumes();
445 iDof_.correctBoundaryConditions();
447 momentum *= nParticle_/
mesh().cellVolumes();
448 momentum_.correctBoundaryConditions();
454 template<
class ParcelType>
464 this->addParticle(
new ParcelType(mesh_, position, celli, U, Ei, typeId));
470 template<
class ParcelType>
473 const word& cloudName,
480 cloudName_(cloudName),
486 cloudName +
"Properties",
493 typeIdList_(particleProperties_.lookup(
"typeIdList")),
494 nParticle_(
readScalar(particleProperties_.lookup(
"nEquivalentParticles"))),
495 cellOccupancy_(mesh_.nCells()),
500 this->
name() +
"SigmaTcRMax",
501 mesh_.time().timeName(),
508 collisionSelectionRemainder_
512 this->
name() +
":collisionSelectionRemainder",
513 mesh_.time().timeName(),
524 mesh_.time().timeName(),
536 mesh_.time().timeName(),
548 mesh_.time().timeName(),
560 mesh_.time().timeName(),
572 mesh_.time().timeName(),
584 mesh_.time().timeName(),
596 mesh_.time().timeName(),
608 mesh_.time().timeName(),
620 mesh_.time().timeName(),
636 mesh_.time().timeName(),
651 mesh_.time().timeName(),
659 binaryCollisionModel_
667 wallInteractionModel_
685 buildCellOccupancy();
689 forAll(collisionSelectionRemainder_, i)
691 collisionSelectionRemainder_[i] = rndGen_.scalar01();
701 template<
class ParcelType>
704 const word& cloudName,
711 cloudName_(cloudName),
717 cloudName +
"Properties",
724 typeIdList_(particleProperties_.lookup(
"typeIdList")),
725 nParticle_(
readScalar(particleProperties_.lookup(
"nEquivalentParticles"))),
731 this->
name() +
"SigmaTcRMax",
732 mesh_.time().timeName(),
739 zeroGradientFvPatchScalarField::typeName
741 collisionSelectionRemainder_
745 this->
name() +
":collisionSelectionRemainder",
746 mesh_.time().timeName(),
757 mesh_.time().timeName(),
769 this->
name() +
"fD_",
770 mesh_.time().timeName(),
787 this->
name() +
"rhoN_",
788 mesh_.time().timeName(),
800 this->
name() +
"rhoM_",
801 mesh_.time().timeName(),
813 this->
name() +
"dsmcRhoN_",
814 mesh_.time().timeName(),
826 this->
name() +
"linearKE_",
827 mesh_.time().timeName(),
839 this->
name() +
"internalE_",
840 mesh_.time().timeName(),
852 this->
name() +
"iDof_",
853 mesh_.time().timeName(),
865 this->
name() +
"momentum_",
866 mesh_.time().timeName(),
888 mesh_.time().timeName(),
904 mesh_.time().timeName(),
918 binaryCollisionModel_(),
919 wallInteractionModel_(),
920 inflowBoundaryModel_()
924 initialise(dsmcInitialiseDict);
930 template<
class ParcelType>
937 template<
class ParcelType>
940 typename ParcelType::trackingData td(*
this);
947 this->dumpParticlePositions();
951 this->inflowBoundary().inflow();
957 buildCellOccupancy();
967 template<
class ParcelType>
970 label nDSMCParticles = this->size();
973 scalar nMol = nDSMCParticles*nParticle_;
975 vector linearMomentum = linearMomentumOfSystem();
978 scalar linearKineticEnergy = linearKineticEnergyOfSystem();
981 scalar internalEnergy = internalEnergyOfSystem();
985 <<
" Number of dsmc particles = " 991 Info<<
" Number of molecules = " 993 <<
" Mass in system = " 995 <<
" Average linear momentum = " 996 << linearMomentum/nMol <<
nl 997 <<
" |Average linear momentum| = " 998 <<
mag(linearMomentum)/nMol <<
nl 999 <<
" Average linear kinetic energy = " 1000 << linearKineticEnergy/nMol <<
nl 1001 <<
" Average internal energy = " 1002 << internalEnergy/nMol <<
nl 1003 <<
" Average total energy = " 1004 << (internalEnergy + linearKineticEnergy)/nMol
1010 template<
class ParcelType>
1021 rndGen_.GaussNormal(),
1022 rndGen_.GaussNormal(),
1023 rndGen_.GaussNormal()
1028 template<
class ParcelType>
1041 else if (iDof < 2.0 + SMALL && iDof > 2.0 - SMALL)
1048 scalar a = 0.5*iDof - 1;
1054 energyRatio = 10*rndGen_.scalar01();
1055 P =
pow((energyRatio/a), a)*
exp(a - energyRatio);
1056 }
while (P < rndGen_.scalar01());
1065 template<
class ParcelType>
1070 this->db().time().
path()/
"parcelPositions_" 1071 + this->
name() +
"_" 1072 + this->db().time().
timeName() +
".obj" 1077 const ParcelType& p = iter();
1079 pObj<<
"v " << p.position().x()
1080 <<
" " << p.position().y()
1081 <<
" " << p.position().z()
1089 template<
class ParcelType>
1095 cellOccupancy_.setSize(mesh_.nCells());
1096 buildCellOccupancy();
1099 this->inflowBoundary().autoMap(mapper);
tetrahedron< point, const point & > tetPointRef
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)
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.
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
void info() const
Print cloud information.
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.
Vector< scalar > vector
A scalar version of the templated Vector.
void autoMap(const mapPolyMesh &)
Remap the cells of particles corresponding to the.
Class containing mesh-to-mesh mapping information after a change in polyMesh topology.
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...
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.
const Type & value() const
Return const reference to value.
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).
dimensionedScalar pos0(const dimensionedScalar &ds)
const Time & time() const
Return time.
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.
void addNewParcel(const vector &position, const label celli, const vector &U, const scalar Ei, const label typeId)
Add new parcel.
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)
const word cloudName(propsDict.lookup("cloudName"))
Templated wall interaction model class.
void dumpParticlePositions() const
Dump particle positions to .obj file.
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.