DSMCCloud.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration |
5  \\ / A nd | Copyright (C) 2011-2015 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "DSMCCloud.H"
27 #include "BinaryCollisionModel.H"
28 #include "WallInteractionModel.H"
29 #include "InflowBoundaryModel.H"
30 #include "constants.H"
33 
34 using namespace Foam::constant;
35 
36 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
37 
38 template<class ParcelType>
40 {
41  Info<< nl << "Constructing constant properties for" << endl;
42  constProps_.setSize(typeIdList_.size());
43 
44  dictionary moleculeProperties
45  (
46  particleProperties_.subDict("moleculeProperties")
47  );
48 
49  forAll(typeIdList_, i)
50  {
51  const word& id(typeIdList_[i]);
52 
53  Info<< " " << id << endl;
54 
55  const dictionary& molDict(moleculeProperties.subDict(id));
56 
57  constProps_[i] =
58  typename ParcelType::constantProperties::constantProperties(molDict);
59  }
60 }
61 
62 
63 template<class ParcelType>
65 {
66  forAll(cellOccupancy_, cO)
67  {
68  cellOccupancy_[cO].clear();
69  }
70 
71  forAllIter(typename DSMCCloud<ParcelType>, *this, iter)
72  {
73  cellOccupancy_[iter().cell()].append(&iter());
74  }
75 }
76 
77 
78 template<class ParcelType>
80 (
81  const IOdictionary& dsmcInitialiseDict
82 )
83 {
84  Info<< nl << "Initialising particles" << endl;
85 
86  const scalar temperature
87  (
88  readScalar(dsmcInitialiseDict.lookup("temperature"))
89  );
90 
91  const vector velocity(dsmcInitialiseDict.lookup("velocity"));
92 
93  const dictionary& numberDensitiesDict
94  (
95  dsmcInitialiseDict.subDict("numberDensities")
96  );
97 
98  List<word> molecules(numberDensitiesDict.toc());
99 
100  Field<scalar> numberDensities(molecules.size());
101 
102  forAll(molecules, i)
103  {
104  numberDensities[i] = readScalar
105  (
106  numberDensitiesDict.lookup(molecules[i])
107  );
108  }
109 
110  numberDensities /= nParticle_;
111 
112  forAll(mesh_.cells(), cellI)
113  {
114  List<tetIndices> cellTets = polyMeshTetDecomposition::cellTetIndices
115  (
116  mesh_,
117  cellI
118  );
119 
120  forAll(cellTets, tetI)
121  {
122  const tetIndices& cellTetIs = cellTets[tetI];
123  tetPointRef tet = cellTetIs.tet(mesh_);
124  scalar tetVolume = tet.mag();
125 
126  forAll(molecules, i)
127  {
128  const word& moleculeName(molecules[i]);
129 
130  label typeId(findIndex(typeIdList_, moleculeName));
131 
132  if (typeId == -1)
133  {
134  FatalErrorIn("Foam::DSMCCloud<ParcelType>::initialise")
135  << "typeId " << moleculeName << "not defined." << nl
136  << abort(FatalError);
137  }
138 
139  const typename ParcelType::constantProperties& cP =
140  constProps(typeId);
141 
142  scalar numberDensity = numberDensities[i];
143 
144  // Calculate the number of particles required
145  scalar particlesRequired = numberDensity*tetVolume;
146 
147  // Only integer numbers of particles can be inserted
148  label nParticlesToInsert = label(particlesRequired);
149 
150  // Add another particle with a probability proportional to the
151  // remainder of taking the integer part of particlesRequired
152  if
153  (
154  (particlesRequired - nParticlesToInsert)
155  > rndGen_.scalar01()
156  )
157  {
158  nParticlesToInsert++;
159  }
160 
161  for (label pI = 0; pI < nParticlesToInsert; pI++)
162  {
163  point p = tet.randomPoint(rndGen_);
164 
165  vector U = equipartitionLinearVelocity
166  (
167  temperature,
168  cP.mass()
169  );
170 
171  scalar Ei = equipartitionInternalEnergy
172  (
173  temperature,
174  cP.internalDegreesOfFreedom()
175  );
176 
177  U += velocity;
178 
179  addNewParcel
180  (
181  p,
182  U,
183  Ei,
184  cellI,
185  cellTetIs.face(),
186  cellTetIs.tetPt(),
187  typeId
188  );
189  }
190  }
191  }
192  }
193 
194  // Initialise the sigmaTcRMax_ field to the product of the cross section of
195  // the most abundant species and the most probable thermal speed (Bird,
196  // p222-223)
197 
198  label mostAbundantType(findMax(numberDensities));
199 
200  const typename ParcelType::constantProperties& cP = constProps
201  (
202  mostAbundantType
203  );
204 
205  sigmaTcRMax_.internalField() = cP.sigmaT()*maxwellianMostProbableSpeed
206  (
207  temperature,
208  cP.mass()
209  );
210 
211  sigmaTcRMax_.correctBoundaryConditions();
212 }
213 
214 
215 template<class ParcelType>
217 {
218  if (!binaryCollision().active())
219  {
220  return;
221  }
222 
223  // Temporary storage for subCells
224  List<DynamicList<label> > subCells(8);
225 
226  scalar deltaT = mesh().time().deltaTValue();
227 
228  label collisionCandidates = 0;
229 
230  label collisions = 0;
231 
232  forAll(cellOccupancy_, cellI)
233  {
234  const DynamicList<ParcelType*>& cellParcels(cellOccupancy_[cellI]);
235 
236  label nC(cellParcels.size());
237 
238  if (nC > 1)
239  {
240  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
241  // Assign particles to one of 8 Cartesian subCells
242 
243  // Clear temporary lists
244  forAll(subCells, i)
245  {
246  subCells[i].clear();
247  }
248 
249  // Inverse addressing specifying which subCell a parcel is in
250  List<label> whichSubCell(cellParcels.size());
251 
252  const point& cC = mesh_.cellCentres()[cellI];
253 
254  forAll(cellParcels, i)
255  {
256  const ParcelType& p = *cellParcels[i];
257  vector relPos = p.position() - cC;
258 
259  label subCell =
260  pos(relPos.x()) + 2*pos(relPos.y()) + 4*pos(relPos.z());
261 
262  subCells[subCell].append(i);
263  whichSubCell[i] = subCell;
264  }
265 
266  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
267 
268  scalar sigmaTcRMax = sigmaTcRMax_[cellI];
269 
270  scalar selectedPairs =
271  collisionSelectionRemainder_[cellI]
272  + 0.5*nC*(nC - 1)*nParticle_*sigmaTcRMax*deltaT
273  /mesh_.cellVolumes()[cellI];
274 
275  label nCandidates(selectedPairs);
276  collisionSelectionRemainder_[cellI] = selectedPairs - nCandidates;
277  collisionCandidates += nCandidates;
278 
279  for (label c = 0; c < nCandidates; c++)
280  {
281  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
282  // subCell candidate selection procedure
283 
284  // Select the first collision candidate
285  label candidateP = rndGen_.integer(0, nC - 1);
286 
287  // Declare the second collision candidate
288  label candidateQ = -1;
289 
290  List<label> subCellPs = subCells[whichSubCell[candidateP]];
291  label nSC = subCellPs.size();
292 
293  if (nSC > 1)
294  {
295  // If there are two or more particle in a subCell, choose
296  // another from the same cell. If the same candidate is
297  // chosen, choose again.
298 
299  do
300  {
301  candidateQ = subCellPs[rndGen_.integer(0, nSC - 1)];
302  } while (candidateP == candidateQ);
303  }
304  else
305  {
306  // Select a possible second collision candidate from the
307  // whole cell. If the same candidate is chosen, choose
308  // again.
309 
310  do
311  {
312  candidateQ = rndGen_.integer(0, nC - 1);
313  } while (candidateP == candidateQ);
314  }
315 
316  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
317  // uniform candidate selection procedure
318 
319  // // Select the first collision candidate
320  // label candidateP = rndGen_.integer(0, nC-1);
321 
322  // // Select a possible second collision candidate
323  // label candidateQ = rndGen_.integer(0, nC-1);
324 
325  // // If the same candidate is chosen, choose again
326  // while (candidateP == candidateQ)
327  // {
328  // candidateQ = rndGen_.integer(0, nC-1);
329  // }
330 
331  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
332 
333  ParcelType& parcelP = *cellParcels[candidateP];
334  ParcelType& parcelQ = *cellParcels[candidateQ];
335 
336  scalar sigmaTcR = binaryCollision().sigmaTcR
337  (
338  parcelP,
339  parcelQ
340  );
341 
342  // Update the maximum value of sigmaTcR stored, but use the
343  // initial value in the acceptance-rejection criteria because
344  // the number of collision candidates selected was based on this
345 
346  if (sigmaTcR > sigmaTcRMax_[cellI])
347  {
348  sigmaTcRMax_[cellI] = sigmaTcR;
349  }
350 
351  if ((sigmaTcR/sigmaTcRMax) > rndGen_.scalar01())
352  {
353  binaryCollision().collide
354  (
355  parcelP,
356  parcelQ
357  );
358 
359  collisions++;
360  }
361  }
362  }
363  }
364 
365  reduce(collisions, sumOp<label>());
366 
367  reduce(collisionCandidates, sumOp<label>());
368 
369  sigmaTcRMax_.correctBoundaryConditions();
370 
371  if (collisionCandidates)
372  {
373  Info<< " Collisions = "
374  << collisions << nl
375  << " Acceptance rate = "
376  << scalar(collisions)/scalar(collisionCandidates) << nl
377  << endl;
378  }
379  else
380  {
381  Info<< " No collisions" << endl;
382  }
383 }
384 
385 
386 template<class ParcelType>
388 {
389  q_ = dimensionedScalar("zero", dimensionSet(1, 0, -3, 0, 0), 0.0);
390 
391  fD_ = dimensionedVector
392  (
393  "zero",
394  dimensionSet(1, -1, -2, 0, 0),
396  );
397 
398  rhoN_ = dimensionedScalar("zero", dimensionSet(0, -3, 0, 0, 0), VSMALL);
399  rhoM_ = dimensionedScalar("zero", dimensionSet(1, -3, 0, 0, 0), VSMALL);
400  dsmcRhoN_ = dimensionedScalar("zero", dimensionSet(0, -3, 0, 0, 0), 0.0);
401  linearKE_ = dimensionedScalar("zero", dimensionSet(1, -1, -2, 0, 0), 0.0);
402  internalE_ = dimensionedScalar("zero", dimensionSet(1, -1, -2, 0, 0), 0.0);
403  iDof_ = dimensionedScalar("zero", dimensionSet(0, -3, 0, 0, 0), VSMALL);
404 
405  momentum_ = dimensionedVector
406  (
407  "zero",
408  dimensionSet(1, -2, -1, 0, 0),
410  );
411 }
412 
413 
414 template<class ParcelType>
416 {
417  scalarField& rhoN = rhoN_.internalField();
418  scalarField& rhoM = rhoM_.internalField();
419  scalarField& dsmcRhoN = dsmcRhoN_.internalField();
420  scalarField& linearKE = linearKE_.internalField();
421  scalarField& internalE = internalE_.internalField();
422  scalarField& iDof = iDof_.internalField();
423  vectorField& momentum = momentum_.internalField();
424 
425  forAllConstIter(typename DSMCCloud<ParcelType>, *this, iter)
426  {
427  const ParcelType& p = iter();
428  const label cellI = p.cell();
429 
430  rhoN[cellI]++;
431  rhoM[cellI] += constProps(p.typeId()).mass();
432  dsmcRhoN[cellI]++;
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();
437  }
438 
439  rhoN *= nParticle_/mesh().cellVolumes();
440  rhoN_.correctBoundaryConditions();
441 
442  rhoM *= nParticle_/mesh().cellVolumes();
443  rhoM_.correctBoundaryConditions();
444 
445  dsmcRhoN_.correctBoundaryConditions();
446 
447  linearKE *= nParticle_/mesh().cellVolumes();
448  linearKE_.correctBoundaryConditions();
449 
450  internalE *= nParticle_/mesh().cellVolumes();
451  internalE_.correctBoundaryConditions();
452 
453  iDof *= nParticle_/mesh().cellVolumes();
454  iDof_.correctBoundaryConditions();
455 
456  momentum *= nParticle_/mesh().cellVolumes();
457  momentum_.correctBoundaryConditions();
458 }
459 
460 
461 // * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * * //
462 
463 template<class ParcelType>
465 (
466  const vector& position,
467  const vector& U,
468  const scalar Ei,
469  const label cellI,
470  const label tetFaceI,
471  const label tetPtI,
472  const label typeId
473 )
474 {
475  ParcelType* pPtr = new ParcelType
476  (
477  mesh_,
478  position,
479  U,
480  Ei,
481  cellI,
482  tetFaceI,
483  tetPtI,
484  typeId
485  );
486 
487  this->addParticle(pPtr);
488 }
489 
490 
491 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
492 
493 template<class ParcelType>
495 (
496  const word& cloudName,
497  const fvMesh& mesh,
498  bool readFields
499 )
500 :
502  DSMCBaseCloud(),
503  cloudName_(cloudName),
504  mesh_(mesh),
505  particleProperties_
506  (
507  IOobject
508  (
509  cloudName + "Properties",
510  mesh_.time().constant(),
511  mesh_,
514  )
515  ),
516  typeIdList_(particleProperties_.lookup("typeIdList")),
517  nParticle_(readScalar(particleProperties_.lookup("nEquivalentParticles"))),
518  cellOccupancy_(mesh_.nCells()),
519  sigmaTcRMax_
520  (
521  IOobject
522  (
523  this->name() + "SigmaTcRMax",
524  mesh_.time().timeName(),
525  mesh_,
528  ),
529  mesh_
530  ),
531  collisionSelectionRemainder_
532  (
533  IOobject
534  (
535  this->name() + ":collisionSelectionRemainder",
536  mesh_.time().timeName(),
537  mesh_
538  ),
539  mesh_,
540  dimensionedScalar("collisionSelectionRemainder", dimless, 0)
541  ),
542  q_
543  (
544  IOobject
545  (
546  "q",
547  mesh_.time().timeName(),
548  mesh_,
551  ),
552  mesh_
553  ),
554  fD_
555  (
556  IOobject
557  (
558  "fD",
559  mesh_.time().timeName(),
560  mesh_,
563  ),
564  mesh_
565  ),
566  rhoN_
567  (
568  IOobject
569  (
570  "rhoN",
571  mesh_.time().timeName(),
572  mesh_,
575  ),
576  mesh_
577  ),
578  rhoM_
579  (
580  IOobject
581  (
582  "rhoM",
583  mesh_.time().timeName(),
584  mesh_,
587  ),
588  mesh_
589  ),
590  dsmcRhoN_
591  (
592  IOobject
593  (
594  "dsmcRhoN",
595  mesh_.time().timeName(),
596  mesh_,
599  ),
600  mesh_
601  ),
602  linearKE_
603  (
604  IOobject
605  (
606  "linearKE",
607  mesh_.time().timeName(),
608  mesh_,
611  ),
612  mesh_
613  ),
614  internalE_
615  (
616  IOobject
617  (
618  "internalE",
619  mesh_.time().timeName(),
620  mesh_,
623  ),
624  mesh_
625  ),
626  iDof_
627  (
628  IOobject
629  (
630  "iDof",
631  mesh_.time().timeName(),
632  mesh_,
635  ),
636  mesh_
637  ),
638  momentum_
639  (
640  IOobject
641  (
642  "momentum",
643  mesh_.time().timeName(),
644  mesh_,
647  ),
648  mesh_
649  ),
650  constProps_(),
651  rndGen_(label(149382906) + 7183*Pstream::myProcNo()),
652  boundaryT_
653  (
655  (
656  IOobject
657  (
658  "boundaryT",
659  mesh_.time().timeName(),
660  mesh_,
663  ),
664  mesh_
665  )
666  ),
667  boundaryU_
668  (
670  (
671  IOobject
672  (
673  "boundaryU",
674  mesh_.time().timeName(),
675  mesh_,
678  ),
679  mesh_
680  )
681  ),
682  binaryCollisionModel_
683  (
685  (
686  particleProperties_,
687  *this
688  )
689  ),
690  wallInteractionModel_
691  (
693  (
694  particleProperties_,
695  *this
696  )
697  ),
698  inflowBoundaryModel_
699  (
701  (
702  particleProperties_,
703  *this
704  )
705  )
706 {
707  buildConstProps();
708  buildCellOccupancy();
709 
710  // Initialise the collision selection remainder to a random value between 0
711  // and 1.
712  forAll(collisionSelectionRemainder_, i)
713  {
714  collisionSelectionRemainder_[i] = rndGen_.scalar01();
715  }
716 
717  if (readFields)
718  {
719  ParcelType::readFields(*this);
720  }
721 }
722 
723 
724 template<class ParcelType>
726 (
727  const word& cloudName,
728  const fvMesh& mesh,
729  const IOdictionary& dsmcInitialiseDict
730 )
731  :
733  DSMCBaseCloud(),
734  cloudName_(cloudName),
735  mesh_(mesh),
736  particleProperties_
737  (
738  IOobject
739  (
740  cloudName + "Properties",
741  mesh_.time().constant(),
742  mesh_,
745  )
746  ),
747  typeIdList_(particleProperties_.lookup("typeIdList")),
748  nParticle_(readScalar(particleProperties_.lookup("nEquivalentParticles"))),
749  cellOccupancy_(),
750  sigmaTcRMax_
751  (
752  IOobject
753  (
754  this->name() + "SigmaTcRMax",
755  mesh_.time().timeName(),
756  mesh_,
759  ),
760  mesh_,
761  dimensionedScalar("zero", dimensionSet(0, 3, -1, 0, 0), 0.0),
762  zeroGradientFvPatchScalarField::typeName
763  ),
764  collisionSelectionRemainder_
765  (
766  IOobject
767  (
768  this->name() + ":collisionSelectionRemainder",
769  mesh_.time().timeName(),
770  mesh_
771  ),
772  mesh_,
773  dimensionedScalar("collisionSelectionRemainder", dimless, 0)
774  ),
775  q_
776  (
777  IOobject
778  (
779  this->name() + "q_",
780  mesh_.time().timeName(),
781  mesh_,
784  ),
785  mesh_,
786  dimensionedScalar("zero", dimensionSet(1, 0, -3, 0, 0), 0.0)
787  ),
788  fD_
789  (
790  IOobject
791  (
792  this->name() + "fD_",
793  mesh_.time().timeName(),
794  mesh_,
797  ),
798  mesh_,
800  (
801  "zero",
802  dimensionSet(1, -1, -2, 0, 0),
804  )
805  ),
806  rhoN_
807  (
808  IOobject
809  (
810  this->name() + "rhoN_",
811  mesh_.time().timeName(),
812  mesh_,
815  ),
816  mesh_,
817  dimensionedScalar("zero", dimensionSet(0, -3, 0, 0, 0), VSMALL)
818  ),
819  rhoM_
820  (
821  IOobject
822  (
823  this->name() + "rhoM_",
824  mesh_.time().timeName(),
825  mesh_,
828  ),
829  mesh_,
830  dimensionedScalar("zero", dimensionSet(1, -3, 0, 0, 0), VSMALL)
831  ),
832  dsmcRhoN_
833  (
834  IOobject
835  (
836  this->name() + "dsmcRhoN_",
837  mesh_.time().timeName(),
838  mesh_,
841  ),
842  mesh_,
843  dimensionedScalar("zero", dimensionSet(0, -3, 0, 0, 0), 0.0)
844  ),
845  linearKE_
846  (
847  IOobject
848  (
849  this->name() + "linearKE_",
850  mesh_.time().timeName(),
851  mesh_,
854  ),
855  mesh_,
856  dimensionedScalar("zero", dimensionSet(1, -1, -2, 0, 0), 0.0)
857  ),
858  internalE_
859  (
860  IOobject
861  (
862  this->name() + "internalE_",
863  mesh_.time().timeName(),
864  mesh_,
867  ),
868  mesh_,
869  dimensionedScalar("zero", dimensionSet(1, -1, -2, 0, 0), 0.0)
870  ),
871  iDof_
872  (
873  IOobject
874  (
875  this->name() + "iDof_",
876  mesh_.time().timeName(),
877  mesh_,
880  ),
881  mesh_,
882  dimensionedScalar("zero", dimensionSet(0, -3, 0, 0, 0), VSMALL)
883  ),
884  momentum_
885  (
886  IOobject
887  (
888  this->name() + "momentum_",
889  mesh_.time().timeName(),
890  mesh_,
893  ),
894  mesh_,
896  (
897  "zero",
898  dimensionSet(1, -2, -1, 0, 0),
900  )
901  ),
902  constProps_(),
903  rndGen_(label(971501) + 1526*Pstream::myProcNo()),
904  boundaryT_
905  (
907  (
908  IOobject
909  (
910  "boundaryT",
911  mesh_.time().timeName(),
912  mesh_,
915  ),
916  mesh_,
917  dimensionedScalar("zero", dimensionSet(0, 0, 0, 1, 0), 0.0)
918  )
919  ),
920  boundaryU_
921  (
923  (
924  IOobject
925  (
926  "boundaryU",
927  mesh_.time().timeName(),
928  mesh_,
931  ),
932  mesh_,
934  (
935  "zero",
936  dimensionSet(0, 1, -1, 0, 0),
938  )
939  )
940  ),
941  binaryCollisionModel_(),
942  wallInteractionModel_(),
943  inflowBoundaryModel_()
944 {
945  clear();
946  buildConstProps();
947  initialise(dsmcInitialiseDict);
948 }
949 
950 
951 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
952 
953 template<class ParcelType>
955 {}
956 
957 
958 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
959 
960 template<class ParcelType>
962 {
963  typename ParcelType::trackingData td(*this);
964 
965  // Reset the data collection fields
966  resetFields();
967 
968  if (debug)
969  {
970  this->dumpParticlePositions();
971  }
972 
973  // Insert new particles from the inflow boundary
974  this->inflowBoundary().inflow();
975 
976  // Move the particles ballistically with their current velocities
977  Cloud<ParcelType>::move(td, mesh_.time().deltaTValue());
978 
979  // Update cell occupancy
980  buildCellOccupancy();
981 
982  // Calculate new velocities via stochastic collisions
983  collisions();
984 
985  // Calculate the volume field data
986  calculateFields();
987 }
988 
989 
990 template<class ParcelType>
992 {
993  label nDSMCParticles = this->size();
994  reduce(nDSMCParticles, sumOp<label>());
995 
996  scalar nMol = nDSMCParticles*nParticle_;
997 
998  vector linearMomentum = linearMomentumOfSystem();
999  reduce(linearMomentum, sumOp<vector>());
1000 
1001  scalar linearKineticEnergy = linearKineticEnergyOfSystem();
1002  reduce(linearKineticEnergy, sumOp<scalar>());
1003 
1004  scalar internalEnergy = internalEnergyOfSystem();
1005  reduce(internalEnergy, sumOp<scalar>());
1006 
1007  Info<< "Cloud name: " << this->name() << nl
1008  << " Number of dsmc particles = "
1009  << nDSMCParticles
1010  << endl;
1011 
1012  if (nDSMCParticles)
1013  {
1014  Info<< " Number of molecules = "
1015  << nMol << nl
1016  << " Mass in system = "
1017  << returnReduce(massInSystem(), sumOp<scalar>()) << nl
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
1028  << endl;
1029  }
1030 }
1031 
1032 
1033 template<class ParcelType>
1036  scalar temperature,
1037  scalar mass
1038 )
1039 {
1040  return
1041  sqrt(physicoChemical::k.value()*temperature/mass)
1042  *vector
1043  (
1044  rndGen_.GaussNormal(),
1045  rndGen_.GaussNormal(),
1046  rndGen_.GaussNormal()
1047  );
1048 }
1049 
1050 
1051 template<class ParcelType>
1054  scalar temperature,
1055  direction iDof
1056 )
1057 {
1058  scalar Ei = 0.0;
1059 
1060  if (iDof < SMALL)
1061  {
1062  return Ei;
1063  }
1064  else if (iDof < 2.0 + SMALL && iDof > 2.0 - SMALL)
1065  {
1066  // Special case for iDof = 2, i.e. diatomics;
1067  Ei = -log(rndGen_.scalar01())*physicoChemical::k.value()*temperature;
1068  }
1069  else
1070  {
1071  scalar a = 0.5*iDof - 1;
1072  scalar energyRatio;
1073  scalar P = -1;
1074 
1075  do
1076  {
1077  energyRatio = 10*rndGen_.scalar01();
1078  P = pow((energyRatio/a), a)*exp(a - energyRatio);
1079  } while (P < rndGen_.scalar01());
1080 
1081  Ei = energyRatio*physicoChemical::k.value()*temperature;
1082  }
1083 
1084  return Ei;
1085 }
1086 
1087 
1088 template<class ParcelType>
1090 {
1091  OFstream pObj
1092  (
1093  this->db().time().path()/"parcelPositions_"
1094  + this->name() + "_"
1095  + this->db().time().timeName() + ".obj"
1096  );
1097 
1098  forAllConstIter(typename DSMCCloud<ParcelType>, *this, iter)
1099  {
1100  const ParcelType& p = iter();
1101 
1102  pObj<< "v " << p.position().x()
1103  << " " << p.position().y()
1104  << " " << p.position().z()
1105  << nl;
1106  }
1107 
1108  pObj.flush();
1109 }
1110 
1111 
1112 template<class ParcelType>
1114 {
1115  typedef typename ParcelType::trackingData tdType;
1116  tdType td(*this);
1117  Cloud<ParcelType>::template autoMap<tdType>(td, mapper);
1118 
1119  // Update the cell occupancy field
1120  cellOccupancy_.setSize(mesh_.nCells());
1121  buildCellOccupancy();
1122 
1123  // Update the inflow BCs
1124  this->inflowBoundary().autoMap(mapper);
1125 }
1126 
1127 
1128 // ************************************************************************* //
void reduce(const List< UPstream::commsStruct > &comms, T &Value, const BinaryOp &bop, const int tag, const label comm)
dimensionedScalar sqrt(const dimensionedScalar &ds)
Output to file stream.
Definition: OFstream.H:81
unsigned char direction
Definition: direction.H:43
virtual ~DSMCCloud()
Destructor.
Definition: DSMCCloud.C:954
Mesh data needed to do the Finite Volume discretisation.
Definition: fvMesh.H:78
const word cloudName(propsDict.lookup("cloudName"))
T returnReduce(const T &Value, const BinaryOp &bop, const int tag=Pstream::msgType(), const label comm=UPstream::worldComm)
vector point
Point is a vector.
Definition: point.H:41
label findIndex(const ListType &, typename ListType::const_reference, const label start=0)
Find first occurence of given element and return index,.
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
dimensioned< scalar > mag(const dimensioned< Type > &)
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
Collection of constants.
Class containing mesh-to-mesh mapping information after a change in polyMesh topology.
Definition: mapPolyMesh.H:158
Virtual abstract base class for templated DSMCCloud.
Definition: DSMCBaseCloud.H:48
label findMax(const ListType &, const label start=0)
Find index of max element (and larger than given element).
#define forAllIter(Container, container, iter)
Definition: UList.H:440
A class for handling words, derived from string.
Definition: word.H:59
intWM_LABEL_SIZE_t label
A label is an int32_t or int64_t as specified by the pre-processor macro WM_LABEL_SIZE.
Definition: label.H:59
dimensionedScalar exp(const dimensionedScalar &ds)
messageStream Info
void move(TrackData &td, const scalar trackTime)
Move the particles.
Definition: Cloud.C:187
dynamicFvMesh & mesh
const word & constant() const
Return constant name.
Definition: TimePaths.H:124
autoPtr< BasicCompressibleTurbulenceModel > New(const volScalarField &rho, const volVectorField &U, const surfaceScalarField &phi, const typename BasicCompressibleTurbulenceModel::transportModel &transport, const word &propertiesName)
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.
Definition: DSMCCloud.C:465
Templated wall interaction model class.
Definition: DSMCCloud.H:58
tetrahedron< point, const point & > tetPointRef
Definition: tetrahedron.H:78
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
scalar equipartitionInternalEnergy(scalar temperature, direction internalDegreesOfFreedom)
Generate a random internal energy, sampled from the.
Definition: DSMCCloud.C:1053
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
Definition: pEqn.H:39
dimensionedScalar log(const dimensionedScalar &ds)
static const char nl
Definition: Ostream.H:260
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
IOobject defines the attributes of an object for which implicit objectRegistry management is supporte...
Definition: IOobject.H:91
virtual void flush()
Flush stream.
Definition: OSstream.C:246
void readFields(const Mesh &mesh, const IOobjectList &objects, PtrList< GeoField > &fields)
void dumpParticlePositions() const
Dump particle positions to .obj file.
Definition: DSMCCloud.C:1089
Field< vector > vectorField
Specialisation of Field<T> for vector.
Dimension set for the base types.
Definition: dimensionSet.H:116
#define forAll(list, i)
Definition: UList.H:421
bool readScalar(const char *buf, doubleScalar &s)
Read whole of buf as a scalar. Return true if succesful.
Definition: doubleScalar.H:63
IOdictionary is derived from dictionary and IOobject to give the dictionary automatic IO functionalit...
Definition: IOdictionary.H:53
void evolve()
Evolve the cloud (move, collide)
Definition: DSMCCloud.C:961
vector equipartitionLinearVelocity(scalar temperature, scalar mass)
Generate a random velocity sampled from the Maxwellian speed.
Definition: DSMCCloud.C:1035
dimensioned< vector > dimensionedVector
Dimensioned vector obtained from generic dimensioned type.
virtual void autoMap(const mapPolyMesh &)
Remap the particles to the correct cells following mesh change.
Definition: DSMCCloud.C:1113
errorManip< error > abort(error &err)
Definition: errorManip.H:131
static int myProcNo(const label communicator=0)
Number of this process (starting from masterNo() = 0)
Definition: UPstream.H:404
fileName path(UMean.rootPath()/UMean.caseName()/"graphs"/UMean.instance())
#define FatalErrorIn(functionName)
Report an error message using Foam::FatalError.
Definition: error.H:314
void info() const
Print cloud information.
Definition: DSMCCloud.C:991
error FatalError
static const Vector zero
Definition: Vector.H:80
const dimensionedScalar c
Speed of light in a vacuum.
Vector< scalar > vector
A scalar version of the templated Vector.
Definition: vector.H:49
Templated DSMC particle collision class.
Definition: DSMCCloud.H:55
const dimensionSet dimless(0, 0, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:47
const Time & time() const
Return time.
Templated base class for dsmc cloud.
Definition: DSMCCloud.H:68
dimensionedScalar pos(const dimensionedScalar &ds)
Templated inflow boundary model class.
Definition: DSMCCloud.H:61
UEqn clear()
const dimensionedScalar k
Boltzmann constant.
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
static List< tetIndices > cellTetIndices(const polyMesh &mesh, label cI)
Return the tet decomposition of the given cell, see.
const Type & value() const
Return const reference to value.
word timeName
Definition: getTimeIndex.H:3