FreeFOAM The Cross-Platform CFD Toolkit
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) 2009-2010 OpenCFD Ltd.
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"
31 
32 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
33 
34 template<class ParcelType>
35 Foam::scalar Foam::DsmcCloud<ParcelType>::kb = 1.380650277e-23;
36 
37 
38 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
39 
40 template<class ParcelType>
42 {
43  Info<< nl << "Constructing constant properties for" << endl;
44  constProps_.setSize(typeIdList_.size());
45 
46  dictionary moleculeProperties
47  (
48  particleProperties_.subDict("moleculeProperties")
49  );
50 
51  forAll(typeIdList_, i)
52  {
53  const word& id(typeIdList_[i]);
54 
55  Info<< " " << id << endl;
56 
57  const dictionary& molDict(moleculeProperties.subDict(id));
58 
59  constProps_[i] =
60  typename ParcelType::constantProperties::constantProperties(molDict);
61  }
62 }
63 
64 
65 template<class ParcelType>
67 {
68  forAll(cellOccupancy_, cO)
69  {
70  cellOccupancy_[cO].clear();
71  }
72 
73  forAllIter(typename DsmcCloud<ParcelType>, *this, iter)
74  {
75  cellOccupancy_[iter().cell()].append(&iter());
76  }
77 }
78 
79 
80 template<class ParcelType>
82 (
83  const IOdictionary& dsmcInitialiseDict
84 )
85 {
86  Info<< nl << "Initialising particles" << endl;
87 
88  const scalar temperature
89  (
90  readScalar(dsmcInitialiseDict.lookup("temperature"))
91  );
92 
93  const vector velocity(dsmcInitialiseDict.lookup("velocity"));
94 
95  const dictionary& numberDensitiesDict
96  (
97  dsmcInitialiseDict.subDict("numberDensities")
98  );
99 
100  List<word> molecules(numberDensitiesDict.toc());
101 
102  Field<scalar> numberDensities(molecules.size());
103 
104  forAll(molecules, i)
105  {
106  numberDensities[i] = readScalar
107  (
108  numberDensitiesDict.lookup(molecules[i])
109  );
110  }
111 
112  numberDensities /= nParticle_;
113 
114  forAll(mesh_.cells(), cell)
115  {
116  const vector& cC = mesh_.cellCentres()[cell];
117  const labelList& cellFaces = mesh_.cells()[cell];
118  const scalar cV = mesh_.cellVolumes()[cell];
119 
120  label nTets = 0;
121 
122  // Each face is split into nEdges (or nVertices) - 2 tets.
123  forAll(cellFaces, face)
124  {
125  nTets += mesh_.faces()[cellFaces[face]].size() - 2;
126  }
127 
128  // Calculate the cumulative tet volumes circulating around the cell and
129  // record the vertex labels of each.
130  scalarList cTetVFracs(nTets, 0.0);
131 
132  List<labelList> tetPtIs(nTets, labelList(3,-1));
133 
134  // Keep track of which tet this is.
135  label tet = 0;
136 
137  forAll(cellFaces, face)
138  {
139  const labelList& facePoints = mesh_.faces()[cellFaces[face]];
140 
141  label pointI = 1;
142  while ((pointI + 1) < facePoints.size())
143  {
144 
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]];
148 
149  cTetVFracs[tet] =
150  mag(((pA - cC) ^ (pB - cC)) & (pC - cC))/(cV*6.0)
151  + cTetVFracs[max((tet - 1),0)];
152 
153  tetPtIs[tet][0] = facePoints[0];
154  tetPtIs[tet][1] = facePoints[pointI];
155  tetPtIs[tet][2] = facePoints[pointI + 1];
156 
157  pointI++;
158  tet++;
159  }
160  }
161 
162  // Force the last volume fraction value to 1.0 to avoid any
163  // rounding/non-flat face errors giving a value < 1.0
164  cTetVFracs[nTets - 1] = 1.0;
165 
166  forAll(molecules, i)
167  {
168  const word& moleculeName(molecules[i]);
169 
170  label typeId(findIndex(typeIdList_, moleculeName));
171 
172  if (typeId == -1)
173  {
174  FatalErrorIn("Foam::DsmcCloud<ParcelType>::initialise")
175  << "typeId " << moleculeName << "not defined." << nl
176  << abort(FatalError);
177  }
178 
179  const typename ParcelType::constantProperties& cP =
180  constProps(typeId);
181 
182  scalar numberDensity = numberDensities[i];
183 
184  // Calculate the number of particles required
185  scalar particlesRequired = numberDensity*mesh_.cellVolumes()[cell];
186 
187  // Only integer numbers of particles can be inserted
188  label nParticlesToInsert = label(particlesRequired);
189 
190  // Add another particle with a probability proportional to the
191  // remainder of taking the integer part of particlesRequired
192  if ((particlesRequired - nParticlesToInsert) > rndGen_.scalar01())
193  {
194  nParticlesToInsert++;
195  }
196 
197  for (label pI = 0; pI < nParticlesToInsert; pI++)
198  {
199  // Choose a random point in a generic tetrahedron
200 
201  scalar s = rndGen_.scalar01();
202  scalar t = rndGen_.scalar01();
203  scalar u = rndGen_.scalar01();
204 
205  if (s + t > 1.0)
206  {
207  s = 1.0 - s;
208  t = 1.0 - t;
209  }
210 
211  if (t + u > 1.0)
212  {
213  scalar tmp = u;
214  u = 1.0 - s - t;
215  t = 1.0 - tmp;
216  }
217  else if (s + t + u > 1.0)
218  {
219  scalar tmp = u;
220  u = s + t + u - 1.0;
221  s = 1.0 - t - tmp;
222  }
223 
224  // Choose a tetrahedron to insert in, based on their relative
225  // volumes
226  scalar tetSelection = rndGen_.scalar01();
227 
228  // Selected tetrahedron
229  label sTet = -1;
230 
231  forAll(cTetVFracs, tet)
232  {
233  sTet = tet;
234 
235  if (cTetVFracs[tet] >= tetSelection)
236  {
237  break;
238  }
239  }
240 
241  vector p =
242  (1 - s - t - u)*cC
243  + s*mesh_.points()[tetPtIs[sTet][0]]
244  + t*mesh_.points()[tetPtIs[sTet][1]]
245  + u*mesh_.points()[tetPtIs[sTet][2]];
246 
247  vector U = equipartitionLinearVelocity
248  (
249  temperature,
250  cP.mass()
251  );
252 
253  scalar Ei = equipartitionInternalEnergy
254  (
255  temperature,
256  cP.internalDegreesOfFreedom()
257  );
258 
259  U += velocity;
260 
261  addNewParcel
262  (
263  p,
264  U,
265  Ei,
266  cell,
267  typeId
268  );
269  }
270  }
271  }
272 
273  // Initialise the sigmaTcRMax_ field to the product of the cross section of
274  // the most abundant species and the most probable thermal speed (Bird,
275  // p222-223)
276 
277  label mostAbundantType(findMax(numberDensities));
278 
279  const typename ParcelType::constantProperties& cP = constProps
280  (
281  mostAbundantType
282  );
283 
284  sigmaTcRMax_.internalField() = cP.sigmaT()*maxwellianMostProbableSpeed
285  (
286  temperature,
287  cP.mass()
288  );
289 
290  sigmaTcRMax_.correctBoundaryConditions();
291 }
292 
293 
294 template<class ParcelType>
296 {
297  buildCellOccupancy();
298 
299  // Temporary storage for subCells
300  List<DynamicList<label> > subCells(8);
301 
302  scalar deltaT = mesh().time().deltaTValue();
303 
304  label collisionCandidates = 0;
305 
306  label collisions = 0;
307 
308  forAll(cellOccupancy_, celli)
309  {
310  const DynamicList<ParcelType*>& cellParcels(cellOccupancy_[celli]);
311 
312  label nC(cellParcels.size());
313 
314  if (nC > 1)
315  {
316 
317  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
318  // Assign particles to one of 8 Cartesian subCells
319 
320  // Clear temporary lists
321  forAll(subCells, i)
322  {
323  subCells[i].clear();
324  }
325 
326  // Inverse addressing specifying which subCell a parcel is in
327  List<label> whichSubCell(cellParcels.size());
328 
329  const point& cC = mesh_.cellCentres()[celli];
330 
331  forAll(cellParcels, i)
332  {
333  ParcelType* p = cellParcels[i];
334 
335  vector relPos = p->position() - cC;
336 
337  label subCell =
338  pos(relPos.x()) + 2*pos(relPos.y()) + 4*pos(relPos.z());
339 
340  subCells[subCell].append(i);
341 
342  whichSubCell[i] = subCell;
343  }
344 
345  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
346 
347  scalar sigmaTcRMax = sigmaTcRMax_[celli];
348 
349  scalar selectedPairs =
350  collisionSelectionRemainder_[celli]
351  + 0.5*nC*(nC - 1)*nParticle_*sigmaTcRMax*deltaT
352  /mesh_.cellVolumes()[celli];
353 
354  label nCandidates(selectedPairs);
355 
356  collisionSelectionRemainder_[celli] = selectedPairs - nCandidates;
357 
358  collisionCandidates += nCandidates;
359 
360  for (label c = 0; c < nCandidates; c++)
361  {
362  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
363  // subCell candidate selection procedure
364 
365  // Select the first collision candidate
366  label candidateP = rndGen_.integer(0, nC - 1);
367 
368  // Declare the second collision candidate
369  label candidateQ = -1;
370 
371  List<label> subCellPs = subCells[whichSubCell[candidateP]];
372 
373  label nSC = subCellPs.size();
374 
375  if (nSC > 1)
376  {
377  // If there are two or more particle in a subCell, choose
378  // another from the same cell. If the same candidate is
379  // chosen, choose again.
380 
381  do
382  {
383  candidateQ = subCellPs[rndGen_.integer(0, nSC - 1)];
384 
385  } while(candidateP == candidateQ);
386  }
387  else
388  {
389  // Select a possible second collision candidate from the
390  // whole cell. If the same candidate is chosen, choose
391  // again.
392 
393  do
394  {
395  candidateQ = rndGen_.integer(0, nC - 1);
396 
397  } while(candidateP == candidateQ);
398  }
399 
400  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
401  // uniform candidate selection procedure
402 
403  // // Select the first collision candidate
404  // label candidateP = rndGen_.integer(0, nC-1);
405 
406  // // Select a possible second collision candidate
407  // label candidateQ = rndGen_.integer(0, nC-1);
408 
409  // // If the same candidate is chosen, choose again
410  // while(candidateP == candidateQ)
411  // {
412  // candidateQ = rndGen_.integer(0, nC-1);
413  // }
414 
415  // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
416 
417  ParcelType* parcelP = cellParcels[candidateP];
418  ParcelType* parcelQ = cellParcels[candidateQ];
419 
420  label typeIdP = parcelP->typeId();
421  label typeIdQ = parcelQ->typeId();
422 
423  scalar sigmaTcR = binaryCollision().sigmaTcR
424  (
425  typeIdP,
426  typeIdQ,
427  parcelP->U(),
428  parcelQ->U()
429  );
430 
431  // Update the maximum value of sigmaTcR stored, but use the
432  // initial value in the acceptance-rejection criteria because
433  // the number of collision candidates selected was based on this
434 
435  if (sigmaTcR > sigmaTcRMax_[celli])
436  {
437  sigmaTcRMax_[celli] = sigmaTcR;
438  }
439 
440  if ((sigmaTcR/sigmaTcRMax) > rndGen_.scalar01())
441  {
442  binaryCollision().collide
443  (
444  typeIdP,
445  typeIdQ,
446  parcelP->U(),
447  parcelQ->U(),
448  parcelP->Ei(),
449  parcelQ->Ei()
450  );
451 
452  collisions++;
453  }
454  }
455  }
456  }
457 
458  reduce(collisions, sumOp<label>());
459 
460  reduce(collisionCandidates, sumOp<label>());
461 
462  sigmaTcRMax_.correctBoundaryConditions();
463 
464  if (collisionCandidates)
465  {
466  Info<< " Collisions = "
467  << collisions << nl
468  << " Acceptance rate = "
469  << scalar(collisions)/scalar(collisionCandidates) << nl
470  << endl;
471  }
472  else
473  {
474  Info<< " No collisions" << endl;
475  }
476 }
477 
478 
479 template<class ParcelType>
481 {
482  q_ = dimensionedScalar("zero", dimensionSet(1, 0, -3, 0, 0), 0.0);
483 
484  fD_ = dimensionedVector
485  (
486  "zero",
487  dimensionSet(1, -1, -2, 0, 0),
488  vector::zero
489  );
490 
491  rhoN_ = dimensionedScalar("zero", dimensionSet(0, -3, 0, 0, 0), VSMALL);
492 
493  rhoM_ = dimensionedScalar("zero", dimensionSet(1, -3, 0, 0, 0), VSMALL);
494 
495  dsmcRhoN_ = dimensionedScalar("zero", dimensionSet(0, -3, 0, 0, 0), 0.0);
496 
497  linearKE_ = dimensionedScalar("zero", dimensionSet(1, -1, -2, 0, 0), 0.0);
498 
499  internalE_ = dimensionedScalar("zero", dimensionSet(1, -1, -2, 0, 0), 0.0);
500 
501  iDof_ = dimensionedScalar("zero", dimensionSet(0, -3, 0, 0, 0), VSMALL);
502 
503  momentum_ = dimensionedVector
504  (
505  "zero",
506  dimensionSet(1, -2, -1, 0, 0),
507  vector::zero
508  );
509 }
510 
511 
512 template<class ParcelType>
514 {
515  scalarField& rhoN = rhoN_.internalField();
516 
517  scalarField& rhoM = rhoM_.internalField();
518 
519  scalarField& dsmcRhoN = dsmcRhoN_.internalField();
520 
521  scalarField& linearKE = linearKE_.internalField();
522 
523  scalarField& internalE = internalE_.internalField();
524 
525  scalarField& iDof = iDof_.internalField();
526 
527  vectorField& momentum = momentum_.internalField();
528 
529  forAllConstIter(typename DsmcCloud<ParcelType>, *this, iter)
530  {
531  const ParcelType& p = iter();
532  const label cellI = p.cell();
533 
534  rhoN[cellI]++;
535 
536  rhoM[cellI] += constProps(p.typeId()).mass();
537 
538  dsmcRhoN[cellI]++;
539 
540  linearKE[cellI] += 0.5*constProps(p.typeId()).mass()*(p.U() & p.U());
541 
542  internalE[cellI] += p.Ei();
543 
544  iDof[cellI] += constProps(p.typeId()).internalDegreesOfFreedom();
545 
546  momentum[cellI] += constProps(p.typeId()).mass()*p.U();
547  }
548 
549  rhoN *= nParticle_/mesh().cellVolumes();
550  rhoN_.correctBoundaryConditions();
551 
552  rhoM *= nParticle_/mesh().cellVolumes();
553  rhoM_.correctBoundaryConditions();
554 
555  dsmcRhoN_.correctBoundaryConditions();
556 
557  linearKE *= nParticle_/mesh().cellVolumes();
558  linearKE_.correctBoundaryConditions();
559 
560  internalE *= nParticle_/mesh().cellVolumes();
561  internalE_.correctBoundaryConditions();
562 
563  iDof *= nParticle_/mesh().cellVolumes();
564  iDof_.correctBoundaryConditions();
565 
566  momentum *= nParticle_/mesh().cellVolumes();
567  momentum_.correctBoundaryConditions();
568 }
569 
570 
571 // * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * * //
572 
573 template<class ParcelType>
575 (
576  const vector& position,
577  const vector& U,
578  const scalar Ei,
579  const label cellId,
580  const label typeId
581 )
582 {
583  ParcelType* pPtr = new ParcelType
584  (
585  *this,
586  position,
587  U,
588  Ei,
589  cellId,
590  typeId
591  );
592 
593  addParticle(pPtr);
594 }
595 
596 
597 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
598 
599 template<class ParcelType>
601 (
602  const word& cloudName,
603  const fvMesh& mesh,
604  bool readFields
605 )
606 :
608  DsmcBaseCloud(),
609  cloudName_(cloudName),
610  mesh_(mesh),
611  particleProperties_
612  (
613  IOobject
614  (
615  cloudName + "Properties",
616  mesh_.time().constant(),
617  mesh_,
618  IOobject::MUST_READ,
619  IOobject::NO_WRITE
620  )
621  ),
622  typeIdList_(particleProperties_.lookup("typeIdList")),
623  nParticle_(readScalar(particleProperties_.lookup("nEquivalentParticles"))),
624  cellOccupancy_(mesh_.nCells()),
625  sigmaTcRMax_
626  (
627  IOobject
628  (
629  this->name() + "SigmaTcRMax",
630  mesh_.time().timeName(),
631  mesh_,
632  IOobject::MUST_READ,
633  IOobject::AUTO_WRITE
634  ),
635  mesh_
636  ),
637  collisionSelectionRemainder_(mesh_.nCells(), 0),
638  q_
639  (
640  IOobject
641  (
642  "q",
643  mesh_.time().timeName(),
644  mesh_,
645  IOobject::MUST_READ,
646  IOobject::AUTO_WRITE
647  ),
648  mesh_
649  ),
650  fD_
651  (
652  IOobject
653  (
654  "fD",
655  mesh_.time().timeName(),
656  mesh_,
657  IOobject::MUST_READ,
658  IOobject::AUTO_WRITE
659  ),
660  mesh_
661  ),
662  rhoN_
663  (
664  IOobject
665  (
666  "rhoN",
667  mesh_.time().timeName(),
668  mesh_,
669  IOobject::MUST_READ,
670  IOobject::AUTO_WRITE
671  ),
672  mesh_
673  ),
674  rhoM_
675  (
676  IOobject
677  (
678  "rhoM",
679  mesh_.time().timeName(),
680  mesh_,
681  IOobject::MUST_READ,
682  IOobject::AUTO_WRITE
683  ),
684  mesh_
685  ),
686  dsmcRhoN_
687  (
688  IOobject
689  (
690  "dsmcRhoN",
691  mesh_.time().timeName(),
692  mesh_,
693  IOobject::MUST_READ,
694  IOobject::AUTO_WRITE
695  ),
696  mesh_
697  ),
698  linearKE_
699  (
700  IOobject
701  (
702  "linearKE",
703  mesh_.time().timeName(),
704  mesh_,
705  IOobject::MUST_READ,
706  IOobject::AUTO_WRITE
707  ),
708  mesh_
709  ),
710  internalE_
711  (
712  IOobject
713  (
714  "internalE",
715  mesh_.time().timeName(),
716  mesh_,
717  IOobject::MUST_READ,
718  IOobject::AUTO_WRITE
719  ),
720  mesh_
721  ),
722  iDof_
723  (
724  IOobject
725  (
726  "iDof",
727  mesh_.time().timeName(),
728  mesh_,
729  IOobject::MUST_READ,
730  IOobject::AUTO_WRITE
731  ),
732  mesh_
733  ),
734  momentum_
735  (
736  IOobject
737  (
738  "momentum",
739  mesh_.time().timeName(),
740  mesh_,
741  IOobject::MUST_READ,
742  IOobject::AUTO_WRITE
743  ),
744  mesh_
745  ),
746  constProps_(),
747  rndGen_(label(149382906) + 7183*Pstream::myProcNo()),
748  boundaryT_
749  (
751  (
752  IOobject
753  (
754  "boundaryT",
755  mesh_.time().timeName(),
756  mesh_,
757  IOobject::MUST_READ,
758  IOobject::AUTO_WRITE
759  ),
760  mesh_
761  )
762  ),
763  boundaryU_
764  (
766  (
767  IOobject
768  (
769  "boundaryU",
770  mesh_.time().timeName(),
771  mesh_,
772  IOobject::MUST_READ,
773  IOobject::AUTO_WRITE
774  ),
775  mesh_
776  )
777  ),
778  binaryCollisionModel_
779  (
781  (
782  particleProperties_,
783  *this
784  )
785  ),
786  wallInteractionModel_
787  (
789  (
790  particleProperties_,
791  *this
792  )
793  ),
794  inflowBoundaryModel_
795  (
797  (
798  particleProperties_,
799  *this
800  )
801  )
802 {
803  buildConstProps();
804 
805  buildCellOccupancy();
806 
807  // Initialise the collision selection remainder to a random value between 0
808  // and 1.
809  forAll(collisionSelectionRemainder_, i)
810  {
811  collisionSelectionRemainder_[i] = rndGen_.scalar01();
812  }
813 
814  if (readFields)
815  {
816  ParcelType::readFields(*this);
817  }
818 }
819 
820 
821 template<class ParcelType>
823 (
824  const word& cloudName,
825  const fvMesh& mesh,
826  const IOdictionary& dsmcInitialiseDict
827 )
828  :
830  DsmcBaseCloud(),
831  cloudName_(cloudName),
832  mesh_(mesh),
833  particleProperties_
834  (
835  IOobject
836  (
837  cloudName + "Properties",
838  mesh_.time().constant(),
839  mesh_,
840  IOobject::MUST_READ,
841  IOobject::NO_WRITE
842  )
843  ),
844  typeIdList_(particleProperties_.lookup("typeIdList")),
845  nParticle_(readScalar(particleProperties_.lookup("nEquivalentParticles"))),
846  cellOccupancy_(),
847  sigmaTcRMax_
848  (
849  IOobject
850  (
851  this->name() + "SigmaTcRMax",
852  mesh_.time().timeName(),
853  mesh_,
854  IOobject::NO_READ,
855  IOobject::AUTO_WRITE
856  ),
857  mesh_,
858  dimensionedScalar("zero", dimensionSet(0, 3, -1, 0, 0), 0.0),
859  zeroGradientFvPatchScalarField::typeName
860  ),
861  collisionSelectionRemainder_(),
862  q_
863  (
864  IOobject
865  (
866  this->name() + "q_",
867  mesh_.time().timeName(),
868  mesh_,
869  IOobject::NO_READ,
870  IOobject::NO_WRITE
871  ),
872  mesh_,
873  dimensionedScalar("zero", dimensionSet(1, 0, -3, 0, 0), 0.0)
874  ),
875  fD_
876  (
877  IOobject
878  (
879  this->name() + "fD_",
880  mesh_.time().timeName(),
881  mesh_,
882  IOobject::NO_READ,
883  IOobject::NO_WRITE
884  ),
885  mesh_,
887  (
888  "zero",
889  dimensionSet(1, -1, -2, 0, 0),
890  vector::zero
891  )
892  ),
893  rhoN_
894  (
895  IOobject
896  (
897  this->name() + "rhoN_",
898  mesh_.time().timeName(),
899  mesh_,
900  IOobject::NO_READ,
901  IOobject::NO_WRITE
902  ),
903  mesh_,
904  dimensionedScalar("zero", dimensionSet(0, -3, 0, 0, 0), VSMALL)
905  ),
906  rhoM_
907  (
908  IOobject
909  (
910  this->name() + "rhoM_",
911  mesh_.time().timeName(),
912  mesh_,
913  IOobject::NO_READ,
914  IOobject::NO_WRITE
915  ),
916  mesh_,
917  dimensionedScalar("zero", dimensionSet(1, -3, 0, 0, 0), VSMALL)
918  ),
919  dsmcRhoN_
920  (
921  IOobject
922  (
923  this->name() + "dsmcRhoN_",
924  mesh_.time().timeName(),
925  mesh_,
926  IOobject::NO_READ,
927  IOobject::NO_WRITE
928  ),
929  mesh_,
930  dimensionedScalar("zero", dimensionSet(0, -3, 0, 0, 0), 0.0)
931  ),
932  linearKE_
933  (
934  IOobject
935  (
936  this->name() + "linearKE_",
937  mesh_.time().timeName(),
938  mesh_,
939  IOobject::NO_READ,
940  IOobject::NO_WRITE
941  ),
942  mesh_,
943  dimensionedScalar("zero", dimensionSet(1, -1, -2, 0, 0), 0.0)
944  ),
945  internalE_
946  (
947  IOobject
948  (
949  this->name() + "internalE_",
950  mesh_.time().timeName(),
951  mesh_,
952  IOobject::NO_READ,
953  IOobject::NO_WRITE
954  ),
955  mesh_,
956  dimensionedScalar("zero", dimensionSet(1, -1, -2, 0, 0), 0.0)
957  ),
958  iDof_
959  (
960  IOobject
961  (
962  this->name() + "iDof_",
963  mesh_.time().timeName(),
964  mesh_,
965  IOobject::NO_READ,
966  IOobject::NO_WRITE
967  ),
968  mesh_,
969  dimensionedScalar("zero", dimensionSet(0, -3, 0, 0, 0), VSMALL)
970  ),
971  momentum_
972  (
973  IOobject
974  (
975  this->name() + "momentum_",
976  mesh_.time().timeName(),
977  mesh_,
978  IOobject::NO_READ,
979  IOobject::NO_WRITE
980  ),
981  mesh_,
983  (
984  "zero",
985  dimensionSet(1, -2, -1, 0, 0),
986  vector::zero
987  )
988  ),
989  constProps_(),
990  rndGen_(label(971501) + 1526*Pstream::myProcNo()),
991  boundaryT_
992  (
994  (
995  IOobject
996  (
997  "boundaryT",
998  mesh_.time().timeName(),
999  mesh_,
1000  IOobject::NO_READ,
1001  IOobject::NO_WRITE
1002  ),
1003  mesh_,
1004  dimensionedScalar("zero", dimensionSet(0, 0, 0, 1, 0), 0.0)
1005  )
1006  ),
1007  boundaryU_
1008  (
1010  (
1011  IOobject
1012  (
1013  "boundaryU",
1014  mesh_.time().timeName(),
1015  mesh_,
1016  IOobject::NO_READ,
1017  IOobject::NO_WRITE
1018  ),
1019  mesh_,
1021  (
1022  "zero",
1023  dimensionSet(0, 1, -1, 0, 0),
1024  vector::zero
1025  )
1026  )
1027  ),
1028  binaryCollisionModel_(),
1029  wallInteractionModel_(),
1030  inflowBoundaryModel_()
1031 {
1032  clear();
1033 
1034  buildConstProps();
1035 
1036  initialise(dsmcInitialiseDict);
1037 }
1038 
1039 
1040 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
1041 
1042 template<class ParcelType>
1044 {}
1045 
1046 
1047 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
1048 
1049 template<class ParcelType>
1051 {
1052  typename ParcelType::trackData td(*this);
1053 
1054  // Reset the data collection fields
1055  resetFields();
1056 
1057  if (debug)
1058  {
1059  this->dumpParticlePositions();
1060  }
1061 
1062  // Insert new particles from the inflow boundary
1063  this->inflowBoundary().inflow();
1064 
1065  // Move the particles ballistically with their current velocities
1067 
1068  // Calculate new velocities via stochastic collisions
1069  collisions();
1070 
1071  // Calculate the volume field data
1072  calculateFields();
1073 }
1074 
1075 
1076 template<class ParcelType>
1078 {
1079  label nDsmcParticles = this->size();
1080  reduce(nDsmcParticles, sumOp<label>());
1081 
1082  scalar nMol = nDsmcParticles*nParticle_;
1083 
1084  vector linearMomentum = linearMomentumOfSystem();
1085  reduce(linearMomentum, sumOp<vector>());
1086 
1087  scalar linearKineticEnergy = linearKineticEnergyOfSystem();
1088  reduce(linearKineticEnergy, sumOp<scalar>());
1089 
1090  scalar internalEnergy = internalEnergyOfSystem();
1091  reduce(internalEnergy, sumOp<scalar>());
1092 
1093  Info<< "Cloud name: " << this->name() << nl
1094  << " Number of dsmc particles = "
1095  << nDsmcParticles
1096  << endl;
1097 
1098  if (nDsmcParticles)
1099  {
1100  Info<< " Number of molecules = "
1101  << nMol << nl
1102  << " Mass in system = "
1103  << returnReduce(massInSystem(), sumOp<scalar>()) << nl
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
1114  << endl;
1115  }
1116 }
1117 
1118 
1119 template<class ParcelType>
1122  scalar temperature,
1123  scalar mass
1124 )
1125 {
1126  return
1127  sqrt(kb*temperature/mass)
1128  *vector
1129  (
1130  rndGen_.GaussNormal(),
1131  rndGen_.GaussNormal(),
1132  rndGen_.GaussNormal()
1133  );
1134 }
1135 
1136 
1137 template<class ParcelType>
1140  scalar temperature,
1141  scalar iDof
1142 )
1143 {
1144  scalar Ei = 0.0;
1145 
1146  if (iDof < SMALL)
1147  {
1148  return Ei;
1149  }
1150  else if (iDof < 2.0 + SMALL && iDof > 2.0 - SMALL)
1151  {
1152  // Special case for iDof = 2, i.e. diatomics;
1153  Ei = -log(rndGen_.scalar01())*kb*temperature;
1154  }
1155  else
1156  {
1157  scalar a = 0.5*iDof - 1;
1158 
1159  scalar energyRatio;
1160 
1161  scalar P = -1;
1162 
1163  do
1164  {
1165  energyRatio = 10*rndGen_.scalar01();
1166 
1167  P = pow((energyRatio/a), a)*exp(a - energyRatio);
1168 
1169  } while (P < rndGen_.scalar01());
1170 
1171  Ei = energyRatio*kb*temperature;
1172  }
1173 
1174  return Ei;
1175 }
1176 
1177 
1178 template<class ParcelType>
1180 {
1181  OFstream pObj
1182  (
1183  this->db().time().path()/"parcelPositions_"
1184  + this->name() + "_"
1185  + this->db().time().timeName() + ".obj"
1186  );
1187 
1188  forAllConstIter(typename DsmcCloud<ParcelType>, *this, iter)
1189  {
1190  const ParcelType& p = iter();
1191 
1192  pObj<< "v " << p.position().x()
1193  << " " << p.position().y()
1194  << " " << p.position().z()
1195  << nl;
1196  }
1197 
1198  pObj.flush();
1199 }
1200 
1201 
1202 // ************************ vim: set sw=4 sts=4 et: ************************ //