All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
populationBalanceModel.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2017-2021 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 "populationBalanceModel.H"
27 #include "coalescenceModel.H"
28 #include "breakupModel.H"
29 #include "binaryBreakupModel.H"
30 #include "driftModel.H"
31 #include "nucleationModel.H"
32 #include "phaseSystem.H"
33 #include "surfaceTensionModel.H"
34 #include "fvm.H"
35 #include "fvcDdt.H"
37 #include "shapeModel.H"
38 
39 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
40 
41 namespace Foam
42 {
43 namespace diameterModels
44 {
45  defineTypeNameAndDebug(populationBalanceModel, 0);
46 }
47 }
48 
49 
50 // * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * * //
51 
52 void Foam::diameterModels::populationBalanceModel::registerVelocityGroups()
53 {
54  forAll(fluid_.phases(), phasei)
55  {
56  if (isA<velocityGroup>(fluid_.phases()[phasei].dPtr()()))
57  {
58  const velocityGroup& velGroup =
59  refCast<const velocityGroup>(fluid_.phases()[phasei].dPtr()());
60 
61  if (velGroup.popBalName() == this->name())
62  {
63  velocityGroups_.resize(velocityGroups_.size() + 1);
64 
65  velocityGroups_.set
66  (
67  velocityGroups_.size() - 1,
68  &const_cast<velocityGroup&>(velGroup)
69  );
70 
71  forAll(velGroup.sizeGroups(), i)
72  {
73  this->registerSizeGroups
74  (
75  const_cast<sizeGroup&>(velGroup.sizeGroups()[i])
76  );
77  }
78  }
79  }
80  }
81 }
82 
83 
84 void Foam::diameterModels::populationBalanceModel::registerSizeGroups
85 (
86  sizeGroup& group
87 )
88 {
89  if
90  (
91  sizeGroups().size() != 0
92  &&
93  group.x().value() <= sizeGroups().last().x().value()
94  )
95  {
97  << "Size groups must be entered according to their representative"
98  << " size"
99  << exit(FatalError);
100  }
101 
102  sizeGroups_.resize(sizeGroups().size() + 1);
103  sizeGroups_.set(sizeGroups().size() - 1, &group);
104 
105  if (sizeGroups().size() == 1)
106  {
107  v_.append
108  (
110  (
111  "v",
112  sizeGroups().last().x()
113  )
114  );
115 
116  v_.append
117  (
119  (
120  "v",
121  sizeGroups().last().x()
122  )
123  );
124  }
125  else
126  {
127  v_.last() =
128  0.5
129  *(
130  sizeGroups()[sizeGroups().size()-2].x()
131  + sizeGroups().last().x()
132  );
133 
134  v_.append
135  (
137  (
138  "v",
139  sizeGroups().last().x()
140  )
141  );
142  }
143 
144  delta_.append(new PtrList<dimensionedScalar>());
145 
146  Su_.append
147  (
148  new volScalarField
149  (
150  IOobject
151  (
152  "Su",
153  fluid_.time().timeName(),
154  mesh_
155  ),
156  mesh_,
157  dimensionedScalar(inv(dimTime), 0)
158  )
159  );
160 
161  Sp_.append
162  (
163  new volScalarField
164  (
165  IOobject
166  (
167  "Sp",
168  fluid_.time().timeName(),
169  mesh_
170  ),
171  mesh_,
172  dimensionedScalar(inv(dimTime), 0)
173  )
174  );
175 }
176 
177 
178 void Foam::diameterModels::populationBalanceModel::createPhasePairs()
179 {
180  forAll(velocityGroups_, i)
181  {
182  const phaseModel& phasei = velocityGroups_[i].phase();
183 
184  forAll(velocityGroups_, j)
185  {
186  const phaseModel& phasej = velocityGroups_[j].phase();
187 
188  if (&phasei != &phasej)
189  {
190  const phasePairKey key
191  (
192  phasei.name(),
193  phasej.name(),
194  false
195  );
196 
197  if (!phasePairs_.found(key))
198  {
199  phasePairs_.insert
200  (
201  key,
202  autoPtr<phasePair>
203  (
204  new phasePair
205  (
206  phasei,
207  phasej
208  )
209  )
210  );
211  }
212  }
213  }
214  }
215 }
216 
217 
218 void Foam::diameterModels::populationBalanceModel::precompute()
219 {
220  forAll(coalescenceModels_, model)
221  {
222  coalescenceModels_[model].precompute();
223  }
224 
225  forAll(breakupModels_, model)
226  {
227  breakupModels_[model].precompute();
228 
229  breakupModels_[model].dsdPtr()->precompute();
230  }
231 
232  forAll(binaryBreakupModels_, model)
233  {
234  binaryBreakupModels_[model].precompute();
235  }
236 
237  forAll(drift_, model)
238  {
239  drift_[model].precompute();
240  }
241 
242  forAll(nucleation_, model)
243  {
244  nucleation_[model].precompute();
245  }
246 }
247 
248 
249 void Foam::diameterModels::populationBalanceModel::birthByCoalescence
250 (
251  const label j,
252  const label k
253 )
254 {
255  const sizeGroup& fj = sizeGroups()[j];
256  const sizeGroup& fk = sizeGroups()[k];
257 
258  dimensionedScalar Eta;
259  dimensionedScalar v = fj.x() + fk.x();
260 
261  for (label i = j; i < sizeGroups().size(); i++)
262  {
263  Eta = eta(i, v);
264 
265  if (Eta.value() == 0) continue;
266 
267  const sizeGroup& fi = sizeGroups()[i];
268 
269  if (j == k)
270  {
271  Sui_ =
272  0.5*fi.x()/(fj.x()*fk.x())*Eta
273  *coalescenceRate_()*fj*fj.phase()*fk*fk.phase();
274  }
275  else
276  {
277  Sui_ =
278  fi.x()/(fj.x()*fk.x())*Eta
279  *coalescenceRate_()*fj*fj.phase()*fk*fk.phase();
280  }
281 
282  Su_[i] += Sui_;
283 
284  const phasePairKey pairij
285  (
286  fi.phase().name(),
287  fj.phase().name()
288  );
289 
290  if (pDmdt_.found(pairij))
291  {
292  const scalar dmdtSign
293  (
294  Pair<word>::compare(pDmdt_.find(pairij).key(), pairij)
295  );
296 
297  *pDmdt_[pairij] += dmdtSign*fj.x()/v*Sui_*fj.phase().rho();
298  }
299 
300  const phasePairKey pairik
301  (
302  fi.phase().name(),
303  fk.phase().name()
304  );
305 
306  if (pDmdt_.found(pairik))
307  {
308  const scalar dmdtSign
309  (
310  Pair<word>::compare(pDmdt_.find(pairik).key(), pairik)
311  );
312 
313  *pDmdt_[pairik] += dmdtSign*fk.x()/v*Sui_*fk.phase().rho();
314  }
315 
316  sizeGroups_[i].shapeModelPtr()->addCoalescence(Sui_, fj, fk);
317  }
318 }
319 
320 
321 void Foam::diameterModels::populationBalanceModel::deathByCoalescence
322 (
323  const label i,
324  const label j
325 )
326 {
327  const sizeGroup& fi = sizeGroups()[i];
328  const sizeGroup& fj = sizeGroups()[j];
329 
330  Sp_[i] += coalescenceRate_()*fi.phase()*fj*fj.phase()/fj.x();
331 
332  if (i != j)
333  {
334  Sp_[j] += coalescenceRate_()*fj.phase()*fi*fi.phase()/fi.x();
335  }
336 }
337 
338 
339 void Foam::diameterModels::populationBalanceModel::birthByBreakup
340 (
341  const label k,
342  const label model
343 )
344 {
345  const sizeGroup& fk = sizeGroups()[k];
346 
347  for (label i = 0; i <= k; i++)
348  {
349  const sizeGroup& fi = sizeGroups()[i];
350 
351  Sui_ =
352  fi.x()*breakupModels_[model].dsdPtr()().nik(i, k)/fk.x()
353  *breakupRate_()*fk*fk.phase();
354 
355  Su_[i] += Sui_;
356 
357  const phasePairKey pair
358  (
359  fi.phase().name(),
360  fk.phase().name()
361  );
362 
363  if (pDmdt_.found(pair))
364  {
365  const scalar dmdtSign
366  (
367  Pair<word>::compare(pDmdt_.find(pair).key(), pair)
368  );
369 
370  *pDmdt_[pair] += dmdtSign*Sui_*fk.phase().rho();
371  }
372 
373  sizeGroups_[i].shapeModelPtr()->addBreakup(Sui_, fk);
374  }
375 }
376 
377 
378 void Foam::diameterModels::populationBalanceModel::deathByBreakup(const label i)
379 {
380  Sp_[i] += breakupRate_()*sizeGroups()[i].phase();
381 }
382 
383 
384 void Foam::diameterModels::populationBalanceModel::calcDeltas()
385 {
386  forAll(sizeGroups(), i)
387  {
388  if (delta_[i].empty())
389  {
390  for (label j = 0; j <= sizeGroups().size() - 1; j++)
391  {
392  const sizeGroup& fj = sizeGroups()[j];
393 
394  delta_[i].append
395  (
397  (
398  "delta",
399  dimVolume,
400  v_[i+1].value() - v_[i].value()
401  )
402  );
403 
404  if
405  (
406  v_[i].value() < 0.5*fj.x().value()
407  &&
408  0.5*fj.x().value() < v_[i+1].value()
409  )
410  {
411  delta_[i][j] = mag(0.5*fj.x() - v_[i]);
412  }
413  else if (0.5*fj.x().value() <= v_[i].value())
414  {
415  delta_[i][j].value() = 0;
416  }
417  }
418  }
419  }
420 }
421 
422 
423 void Foam::diameterModels::populationBalanceModel::birthByBinaryBreakup
424 (
425  const label i,
426  const label j
427 )
428 {
429  const sizeGroup& fi = sizeGroups()[i];
430  const sizeGroup& fj = sizeGroups()[j];
431 
432  const volScalarField Su(binaryBreakupRate_()*fj*fj.phase());
433 
434  Sui_ = fi.x()*delta_[i][j]/fj.x()*Su;
435 
436  Su_[i] += Sui_;
437 
438  sizeGroups_[i].shapeModelPtr()->addBreakup(Sui_, fj);
439 
440  const phasePairKey pairij
441  (
442  fi.phase().name(),
443  fj.phase().name()
444  );
445 
446  if (pDmdt_.found(pairij))
447  {
448  const scalar dmdtSign
449  (
450  Pair<word>::compare(pDmdt_.find(pairij).key(), pairij)
451  );
452 
453  *pDmdt_[pairij] += dmdtSign*Sui_*fj.phase().rho();
454  }
455 
456  dimensionedScalar Eta;
457  dimensionedScalar v = fj.x() - fi.x();
458 
459  for (label k = 0; k <= j; k++)
460  {
461  Eta = eta(k, v);
462 
463  if (Eta.value() == 0) continue;
464 
465  const sizeGroup& fk = sizeGroups()[k];
466 
467  volScalarField& Suk = Sui_;
468 
469  Suk = fk.x()*delta_[i][j]*Eta/fj.x()*Su;
470 
471  Su_[k] += Suk;
472 
473  const phasePairKey pairkj
474  (
475  fk.phase().name(),
476  fj.phase().name()
477  );
478 
479  if (pDmdt_.found(pairkj))
480  {
481  const scalar dmdtSign
482  (
484  (
485  pDmdt_.find(pairkj).key(),
486  pairkj
487  )
488  );
489 
490  *pDmdt_[pairkj] += dmdtSign*Suk*fj.phase().rho();
491  }
492 
493  sizeGroups_[k].shapeModelPtr()->addBreakup(Suk, fj);
494  }
495 }
496 
497 
498 void Foam::diameterModels::populationBalanceModel::deathByBinaryBreakup
499 (
500  const label j,
501  const label i
502 )
503 {
504  Sp_[i] += sizeGroups()[i].phase()*binaryBreakupRate_()*delta_[j][i];
505 }
506 
507 
508 void Foam::diameterModels::populationBalanceModel::drift
509 (
510  const label i,
511  driftModel& model
512 )
513 {
514  model.addToDriftRate(driftRate_(), i);
515 
516  const sizeGroup& fp = sizeGroups()[i];
517 
518  if (i < sizeGroups().size() - 1)
519  {
520  const sizeGroup& fe = sizeGroups()[i+1];
521  volScalarField& Sue = Sui_;
522 
523  Sp_[i] += 1/(fe.x() - fp.x())*pos(driftRate_())*driftRate_()*fp.phase();
524 
525  Sue =
526  fe.x()/(fp.x()*(fe.x() - fp.x()))
527  *pos(driftRate_())*driftRate_()*fp*fp.phase();
528 
529  Su_[i+1] += Sue;
530 
531  const phasePairKey pairij
532  (
533  fp.phase().name(),
534  fe.phase().name()
535  );
536 
537  if (pDmdt_.found(pairij))
538  {
539  const scalar dmdtSign
540  (
541  Pair<word>::compare(pDmdt_.find(pairij).key(), pairij)
542  );
543 
544  *pDmdt_[pairij] -= dmdtSign*Sue*fp.phase().rho();
545  }
546 
547  sizeGroups_[i+1].shapeModelPtr()->addDrift(Sue, fp, model);
548  }
549 
550  if (i == sizeGroups().size() - 1)
551  {
552  Sp_[i] -= pos(driftRate_())*driftRate_()*fp.phase()/fp.x();
553  }
554 
555  if (i > 0)
556  {
557  const sizeGroup& fw = sizeGroups()[i-1];
558  volScalarField& Suw = Sui_;
559 
560  Sp_[i] += 1/(fw.x() - fp.x())*neg(driftRate_())*driftRate_()*fp.phase();
561 
562  Suw =
563  fw.x()/(fp.x()*(fw.x() - fp.x()))
564  *neg(driftRate_())*driftRate_()*fp*fp.phase();
565 
566  Su_[i-1] += Suw;
567 
568  const phasePairKey pairih
569  (
570  fp.phase().name(),
571  fw.phase().name()
572  );
573 
574  if (pDmdt_.found(pairih))
575  {
576  const scalar dmdtSign
577  (
578  Pair<word>::compare(pDmdt_.find(pairih).key(), pairih)
579  );
580 
581  *pDmdt_[pairih] -= dmdtSign*Suw*fp.phase().rho();
582  }
583 
584  sizeGroups_[i-1].shapeModelPtr()->addDrift(Suw, fp, model);
585  }
586 
587  if (i == 0)
588  {
589  Sp_[i] -= neg(driftRate_())*driftRate_()*fp.phase()/fp.x();
590  }
591 }
592 
593 
594 void Foam::diameterModels::populationBalanceModel::nucleation
595 (
596  const label i,
597  nucleationModel& model
598 )
599 {
600  const sizeGroup& fi = sizeGroups()[i];
601 
602  model.addToNucleationRate(nucleationRate_(), i);
603 
604  Sui_ = fi.x()*nucleationRate_();
605 
606  Su_[i] += Sui_;
607 
608  sizeGroups_[i].shapeModelPtr()->addNucleation(Sui_, fi, model);
609 }
610 
611 
612 void Foam::diameterModels::populationBalanceModel::sources()
613 {
614  forAll(sizeGroups(), i)
615  {
616  sizeGroups_[i].shapeModelPtr()->reset();
617  Su_[i] = Zero;
618  Sp_[i] = Zero;
619  }
620 
622  (
623  phasePairTable,
624  phasePairs(),
625  phasePairIter
626  )
627  {
628  *pDmdt_(phasePairIter()) = Zero;
629  }
630 
631  forAll(coalescencePairs_, coalescencePairi)
632  {
633  label i = coalescencePairs_[coalescencePairi].first();
634  label j = coalescencePairs_[coalescencePairi].second();
635 
636  coalescenceRate_() = Zero;
637 
638  forAll(coalescenceModels_, model)
639  {
640  coalescenceModels_[model].addToCoalescenceRate
641  (
642  coalescenceRate_(),
643  i,
644  j
645  );
646  }
647 
648  birthByCoalescence(i, j);
649 
650  deathByCoalescence(i, j);
651  }
652 
653  forAll(binaryBreakupPairs_, binaryBreakupPairi)
654  {
655  label i = binaryBreakupPairs_[binaryBreakupPairi].first();
656  label j = binaryBreakupPairs_[binaryBreakupPairi].second();
657 
658  binaryBreakupRate_() = Zero;
659 
660  forAll(binaryBreakupModels_, model)
661  {
662  binaryBreakupModels_[model].addToBinaryBreakupRate
663  (
664  binaryBreakupRate_(),
665  j,
666  i
667  );
668  }
669 
670  birthByBinaryBreakup(j, i);
671 
672  deathByBinaryBreakup(j, i);
673  }
674 
675  forAll(sizeGroups(), i)
676  {
677  forAll(breakupModels_, model)
678  {
679  breakupModels_[model].setBreakupRate(breakupRate_(), i);
680 
681  birthByBreakup(i, model);
682 
683  deathByBreakup(i);
684  }
685 
686  forAll(drift_, model)
687  {
688  driftRate_() = Zero;
689 
690  drift(i, drift_[model]);
691  }
692 
693  forAll(nucleation_, model)
694  {
695  nucleationRate_() = Zero;
696 
697  nucleation(i, nucleation_[model]);
698  }
699  }
700 }
701 
702 
703 void Foam::diameterModels::populationBalanceModel::calcAlphas()
704 {
705  alphas_() = Zero;
706 
707  forAll(velocityGroups_, v)
708  {
709  const phaseModel& phase = velocityGroups_[v].phase();
710 
711  alphas_() += max(phase, phase.residualAlpha());
712  }
713 }
714 
715 
717 Foam::diameterModels::populationBalanceModel::calcDsm()
718 {
719  tmp<volScalarField> tInvDsm
720  (
722  (
723  "invDsm",
724  mesh_,
726  )
727  );
728 
729  volScalarField& invDsm = tInvDsm.ref();
730 
731  forAll(velocityGroups_, v)
732  {
733  const phaseModel& phase = velocityGroups_[v].phase();
734 
735  invDsm += max(phase, phase.residualAlpha())/(phase.d()*alphas_());
736  }
737 
738  return 1/tInvDsm;
739 }
740 
741 
742 void Foam::diameterModels::populationBalanceModel::calcVelocity()
743 {
744  U_() = Zero;
745 
746  forAll(velocityGroups_, v)
747  {
748  const phaseModel& phase = velocityGroups_[v].phase();
749 
750  U_() += max(phase, phase.residualAlpha())*phase.U()/alphas_();
751  }
752 }
753 
754 bool Foam::diameterModels::populationBalanceModel::updateSources()
755 {
756  const bool result = sourceUpdateCounter_ % sourceUpdateInterval() == 0;
757 
758  ++ sourceUpdateCounter_;
759 
760  return result;
761 }
762 
763 
764 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
765 
767 (
768  const phaseSystem& fluid,
769  const word& name,
770  HashPtrTable<volScalarField, phasePairKey, phasePairKey::hash>& pDmdt
771 )
772 :
773  regIOobject
774  (
775  IOobject
776  (
777  name,
778  fluid.time().constant(),
779  fluid.mesh()
780  )
781  ),
782  fluid_(fluid),
783  pDmdt_(pDmdt),
784  mesh_(fluid_.mesh()),
785  name_(name),
786  dict_
787  (
788  fluid_.subDict("populationBalanceCoeffs").subDict(name_)
789  ),
790  pimple_(mesh_.lookupObject<pimpleControl>("solutionControl")),
791  continuousPhase_
792  (
793  mesh_.lookupObject<phaseModel>
794  (
795  IOobject::groupName("alpha", dict_.lookup("continuousPhase"))
796  )
797  ),
798  velocityGroups_(),
799  sizeGroups_(),
800  v_(),
801  delta_(),
802  Su_(),
803  Sp_(),
804  Sui_
805  (
806  IOobject
807  (
808  "Sui",
809  mesh_.time().timeName(),
810  mesh_
811  ),
812  mesh_,
813  dimensionedScalar(inv(dimTime), Zero)
814  ),
815  coalescenceModels_
816  (
817  dict_.lookup("coalescenceModels"),
818  coalescenceModel::iNew(*this)
819  ),
820  coalescenceRate_(),
821  coalescencePairs_(),
822  breakupModels_
823  (
824  dict_.lookup("breakupModels"),
825  breakupModel::iNew(*this)
826  ),
827  breakupRate_(),
828  binaryBreakupModels_
829  (
830  dict_.lookup("binaryBreakupModels"),
831  binaryBreakupModel::iNew(*this)
832  ),
833  binaryBreakupRate_(),
834  binaryBreakupPairs_(),
835  drift_
836  (
837  dict_.lookup("driftModels"),
838  driftModel::iNew(*this)
839  ),
840  driftRate_(),
841  nucleation_
842  (
843  dict_.lookup("nucleationModels"),
844  nucleationModel::iNew(*this)
845  ),
846  nucleationRate_(),
847  alphas_(),
848  dsm_(),
849  U_(),
850  sourceUpdateCounter_(0)
851 {
852  this->registerVelocityGroups();
853 
854  this->createPhasePairs();
855 
856  if (sizeGroups().size() < 3)
857  {
859  << "The populationBalance " << name_
860  << " requires a minimum number of three sizeGroups to be specified."
861  << exit(FatalError);
862  }
863 
864 
865  if (coalescenceModels_.size() != 0)
866  {
867  coalescenceRate_.set
868  (
869  new volScalarField
870  (
871  IOobject
872  (
873  "coalescenceRate",
874  mesh_.time().timeName(),
875  mesh_
876  ),
877  mesh_,
879  )
880  );
881 
882  forAll(sizeGroups(), i)
883  {
884  for (label j = 0; j <= i; j++)
885  {
886  coalescencePairs_.append(labelPair(i, j));
887  }
888  }
889  }
890 
891  if (breakupModels_.size() != 0)
892  {
893  breakupRate_.set
894  (
895  new volScalarField
896  (
897  IOobject
898  (
899  "breakupRate",
900  fluid_.time().timeName(),
901  mesh_
902  ),
903  mesh_,
904  dimensionedScalar(inv(dimTime), Zero)
905  )
906  );
907  }
908 
909  if (binaryBreakupModels_.size() != 0)
910  {
911  binaryBreakupRate_.set
912  (
913  new volScalarField
914  (
915  IOobject
916  (
917  "binaryBreakupRate",
918  fluid_.time().timeName(),
919  mesh_
920  ),
921  mesh_,
923  (
924  "binaryBreakupRate",
925  inv(dimVolume*dimTime),
926  Zero
927  )
928  )
929  );
930 
931  calcDeltas();
932 
933  forAll(sizeGroups(), i)
934  {
935  label j = 0;
936 
937  while (delta_[j][i].value() != 0)
938  {
939  binaryBreakupPairs_.append(labelPair(i, j));
940  j++;
941  }
942  }
943  }
944 
945  if (drift_.size() != 0)
946  {
947  driftRate_.set
948  (
949  new volScalarField
950  (
951  IOobject
952  (
953  "driftRate",
954  fluid_.time().timeName(),
955  mesh_
956  ),
957  mesh_,
959  )
960  );
961  }
962 
963  if (nucleation_.size() != 0)
964  {
965  nucleationRate_.set
966  (
967  new volScalarField
968  (
969  IOobject
970  (
971  "nucleationRate",
972  fluid_.time().timeName(),
973  mesh_
974  ),
975  mesh_,
977  (
978  "nucleationRate",
979  inv(dimTime*dimVolume),
980  Zero
981  )
982  )
983  );
984  }
985 
986  if (velocityGroups_.size() > 1)
987  {
988  alphas_.set
989  (
990  new volScalarField
991  (
992  IOobject
993  (
994  IOobject::groupName("alpha", this->name()),
995  fluid_.time().timeName(),
996  mesh_,
999  ),
1000  mesh_,
1002  )
1003  );
1004 
1005  dsm_.set
1006  (
1007  new volScalarField
1008  (
1009  IOobject
1010  (
1011  IOobject::groupName("d", this->name()),
1012  fluid_.time().timeName(),
1013  mesh_,
1016  ),
1017  mesh_,
1019  )
1020  );
1021 
1022  U_.set
1023  (
1024  new volVectorField
1025  (
1026  IOobject
1027  (
1028  IOobject::groupName("U", this->name()),
1029  fluid_.time().timeName(),
1030  mesh_,
1033  ),
1034  mesh_,
1036  )
1037  );
1038  }
1039 }
1040 
1041 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
1042 
1044 {}
1045 
1046 // * * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * //
1047 
1050 {
1051  notImplemented("populationBalance::clone() const");
1052  return autoPtr<populationBalanceModel>(nullptr);
1053 }
1054 
1055 
1057 {
1058  return os.good();
1059 }
1060 
1061 
1064 (
1065  const label i,
1066  const dimensionedScalar& v
1067 ) const
1068 {
1069  const dimensionedScalar& x0 = sizeGroups()[0].x();
1070  const dimensionedScalar& xi = sizeGroups()[i].x();
1071  const dimensionedScalar& xm = sizeGroups().last().x();
1072  dimensionedScalar lowerBoundary(x0);
1073  dimensionedScalar upperBoundary(xm);
1074 
1075  if (i != 0) lowerBoundary = sizeGroups()[i-1].x();
1076 
1077  if (i != sizeGroups().size() - 1) upperBoundary = sizeGroups()[i+1].x();
1078 
1079  if ((i == 0 && v < x0) || (i == sizeGroups().size() - 1 && v > xm))
1080  {
1081  return v/xi;
1082  }
1083  else if (v < lowerBoundary || v > upperBoundary)
1084  {
1085  return 0;
1086  }
1087  else if (v.value() == xi.value())
1088  {
1089  return 1;
1090  }
1091  else if (v > xi)
1092  {
1093  return (upperBoundary - v)/(upperBoundary - xi);
1094  }
1095  else
1096  {
1097  return (v - lowerBoundary)/(xi - lowerBoundary);
1098  }
1099 }
1100 
1101 
1104 (
1105  const phaseModel& dispersedPhase
1106 ) const
1107 {
1108  return
1109  fluid_.lookupSubModel<surfaceTensionModel>
1110  (
1111  phasePair(dispersedPhase, continuousPhase_)
1112  ).sigma();
1113 }
1114 
1115 
1118 {
1119  return
1120  mesh_.lookupObject<phaseCompressible::momentumTransportModel>
1121  (
1123  (
1124  momentumTransportModel::typeName,
1125  continuousPhase_.name()
1126  )
1127  );
1128 }
1129 
1130 
1132 {
1133  const dictionary& solutionControls = mesh_.solverDict(name_);
1134  bool solveOnFinalIterOnly =
1135  solutionControls.lookupOrDefault<bool>("solveOnFinalIterOnly", false);
1136 
1137  if (!solveOnFinalIterOnly || pimple_.finalPimpleIter())
1138  {
1139  const label nCorr = this->nCorr();
1140  const scalar tolerance =
1141  solutionControls.lookup<scalar>("tolerance");
1142 
1143  if (nCorr > 0)
1144  {
1145  precompute();
1146  }
1147 
1148  int iCorr = 0;
1149  scalar maxInitialResidual = 1;
1150 
1151  while (++iCorr <= nCorr && maxInitialResidual > tolerance)
1152  {
1153  Info<< "populationBalance "
1154  << this->name()
1155  << ": Iteration "
1156  << iCorr
1157  << endl;
1158 
1159  if (updateSources())
1160  {
1161  sources();
1162  }
1163 
1164  maxInitialResidual = 0;
1165 
1166  forAll(sizeGroups(), i)
1167  {
1168  sizeGroup& fi = sizeGroups_[i];
1169  const phaseModel& phase = fi.phase();
1170  const volScalarField& alpha = phase;
1171  const dimensionedScalar& residualAlpha = phase.residualAlpha();
1172  const volScalarField& rho = phase.thermo().rho();
1173 
1174  fvScalarMatrix sizeGroupEqn
1175  (
1176  fvm::ddt(alpha, fi)
1177  + fvm::div(phase.alphaPhi(), fi)
1178  ==
1179  Su_[i]
1180  - fvm::Sp(Sp_[i], fi)
1181  + fluid_.fvModels().source(alpha, rho, fi)/rho
1182  + fvc::ddt(residualAlpha, fi)
1183  - fvm::ddt(residualAlpha, fi)
1184  );
1185 
1186  sizeGroupEqn.relax();
1187  fluid_.fvConstraints().constrain(sizeGroupEqn);
1188 
1189  maxInitialResidual = max
1190  (
1191  sizeGroupEqn.solve().initialResidual(),
1192  maxInitialResidual
1193  );
1194 
1195  fluid_.fvConstraints().constrain(fi);
1196  }
1197  }
1198 
1199  volScalarField fAlpha0
1200  (
1201  sizeGroups().first()*sizeGroups().first().phase()
1202  );
1203 
1204  volScalarField fAlphaN
1205  (
1206  sizeGroups().last()*sizeGroups().last().phase()
1207  );
1208 
1209  Info<< this->name() << " sizeGroup phase fraction first, last = "
1210  << fAlpha0.weightedAverage(this->mesh().V()).value()
1211  << ' ' << fAlphaN.weightedAverage(this->mesh().V()).value()
1212  << endl;
1213  }
1214 }
1215 
1216 
1218 {
1219  if (velocityGroups_.size() > 1)
1220  {
1221  calcAlphas();
1222  dsm_() = calcDsm();
1223  calcVelocity();
1224  }
1225 }
1226 
1227 
1228 // ************************************************************************* //
Templated abstract base class for multiphase compressible turbulence models.
fvMatrix< scalar > fvScalarMatrix
Definition: fvMatricesFwd.H:42
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:434
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
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
label phasei
Definition: pEqn.H:27
error FatalError
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:323
populationBalanceModel(const phaseSystem &fluid, const word &name, HashPtrTable< volScalarField, phasePairKey, phasePairKey::hash > &pDmdt)
dimensionedSphericalTensor inv(const dimensionedSphericalTensor &dt)
static int compare(const Pair< word > &a, const Pair< word > &b)
Compare Pairs.
Definition: Pair.H:153
static tmp< GeometricField< scalar, fvPatchField, volMesh > > New(const word &name, const Internal &, const PtrList< fvPatchField< scalar >> &)
Return a temporary field constructed from name,.
dimensioned< vector > dimensionedVector
Dimensioned vector obtained from generic dimensioned type.
const phaseCompressible::momentumTransportModel & continuousTurbulence() const
Return reference to turbulence model of the continuous phase.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
const dimensionSet dimless
label k
Boltzmann constant.
Generic dimensioned Type class.
bool writeData(Ostream &) const
Dummy write for regIOobject.
tmp< fvMatrix< Type > > Sp(const volScalarField::Internal &, const GeometricField< Type, fvPatchField, volMesh > &)
GeometricField< vector, fvPatchField, volMesh > volVectorField
Definition: volFieldsFwd.H:58
dimensionedScalar neg(const dimensionedScalar &ds)
const dimensionedScalar eta(const label i, const dimensionedScalar &v) const
Return allocation coefficient.
const dimensionSet dimLength
tmp< GeometricField< Type, fvPatchField, volMesh > > ddt(const dimensioned< Type > dt, const fvMesh &mesh)
Definition: fvcDdt.C:45
autoPtr< populationBalanceModel > clone() const
Return clone.
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:57
const dimensionedScalar sigma
Stefan-Boltzmann constant: default SI units: [W/m^2/K^4].
Calculate the first temporal derivative.
dimensionedScalar pos(const dimensionedScalar &ds)
stressControl lookup("compactNormalStress") >> compactNormalStress
dynamicFvMesh & mesh
static word groupName(Name name, const word &group)
virtual void precompute()
Precompute diameter independent expressions.
tmp< fvMatrix< Type > > ddt(const GeometricField< Type, fvPatchField, volMesh > &vf)
Definition: fvmDdt.C:46
const Type & value() const
Return const reference to value.
void correct()
Correct derived quantities.
Pair< label > labelPair
Label pair.
Definition: labelPair.H:48
word timeName
Definition: getTimeIndex.H:3
static const zero Zero
Definition: zero.H:97
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
Definition: pEqn.H:29
tmp< fvMatrix< Type > > div(const surfaceScalarField &flux, const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvmDiv.C:46
const dimensionSet dimVelocity
defineTypeNameAndDebug(combustionModel, 0)
word name(const complex &)
Return a string representation of a complex.
Definition: complex.C:47
Internal & ref()
Return a reference to the dimensioned internal field.
#define notImplemented(functionName)
Issue a FatalErrorIn for a function not currently implemented.
Definition: error.H:360
PhaseCompressibleMomentumTransportModel< dynamicTransportModel > momentumTransportModel
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
int nCorr
Definition: createControls.H:3
const dimensionSet dimVolume
messageStream Info
dimensioned< scalar > mag(const dimensioned< Type > &)
An auto-pointer similar to the STL auto_ptr but with automatic casting to a reference to the type and...
Definition: PtrList.H:52
const tmp< volScalarField > sigmaWithContinuousPhase(const phaseModel &dispersedPhase) const
Return the surface tension coefficient between a given dispersed.
const tmp< volScalarField::Internal > & Su
Definition: alphaSuSp.H:6
A class for managing temporary objects.
Definition: PtrList.H:53
void solve()
Solve the population balance equation.
Namespace for OpenFOAM.