multiphaseSystem.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-2016 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 "multiphaseSystem.H"
27 #include "alphaContactAngleFvPatchScalarField.H"
29 #include "Time.H"
30 #include "subCycle.H"
31 #include "MULES.H"
32 #include "surfaceInterpolate.H"
33 #include "fvcGrad.H"
34 #include "fvcSnGrad.H"
35 #include "fvcDiv.H"
36 #include "fvcFlux.H"
37 #include "fvcAverage.H"
38 
39 // * * * * * * * * * * * * * * * Static Member Data * * * * * * * * * * * * //
40 
41 const Foam::scalar Foam::multiphaseSystem::convertToRad =
43 
44 
45 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
46 
47 void Foam::multiphaseSystem::calcAlphas()
48 {
49  scalar level = 0.0;
50  alphas_ == 0.0;
51 
52  forAllIter(PtrDictionary<phaseModel>, phases_, iter)
53  {
54  alphas_ += level*iter();
55  level += 1.0;
56  }
57 }
58 
59 
60 void Foam::multiphaseSystem::solveAlphas()
61 {
62  PtrList<surfaceScalarField> alphaPhiCorrs(phases_.size());
63  int phasei = 0;
64 
65  forAllIter(PtrDictionary<phaseModel>, phases_, iter)
66  {
67  phaseModel& phase1 = iter();
68  volScalarField& alpha1 = phase1;
69 
70  phase1.alphaPhi() =
71  dimensionedScalar("0", dimensionSet(0, 3, -1, 0, 0), 0);
72 
73  alphaPhiCorrs.set
74  (
75  phasei,
77  (
78  "phi" + alpha1.name() + "Corr",
79  fvc::flux
80  (
81  phi_,
82  phase1,
83  "div(phi," + alpha1.name() + ')'
84  )
85  )
86  );
87 
88  surfaceScalarField& alphaPhiCorr = alphaPhiCorrs[phasei];
89 
90  forAllIter(PtrDictionary<phaseModel>, phases_, iter2)
91  {
92  phaseModel& phase2 = iter2();
93  volScalarField& alpha2 = phase2;
94 
95  if (&phase2 == &phase1) continue;
96 
97  surfaceScalarField phir(phase1.phi() - phase2.phi());
98 
100  (
101  cAlphas_.find(interfacePair(phase1, phase2))
102  );
103 
104  if (cAlpha != cAlphas_.end())
105  {
107  (
108  (mag(phi_) + mag(phir))/mesh_.magSf()
109  );
110 
111  phir += min(cAlpha()*phic, max(phic))*nHatf(phase1, phase2);
112  }
113 
114  word phirScheme
115  (
116  "div(phir," + alpha2.name() + ',' + alpha1.name() + ')'
117  );
118 
119  alphaPhiCorr += fvc::flux
120  (
121  -fvc::flux(-phir, phase2, phirScheme),
122  phase1,
123  phirScheme
124  );
125  }
126 
127  surfaceScalarField::Boundary& alphaPhiCorrBf =
128  alphaPhiCorr.boundaryFieldRef();
129 
130  // Ensure that the flux at inflow BCs is preserved
131  forAll(alphaPhiCorrBf, patchi)
132  {
133  fvsPatchScalarField& alphaPhiCorrp = alphaPhiCorrBf[patchi];
134 
135  if (!alphaPhiCorrp.coupled())
136  {
137  const scalarField& phi1p = phase1.phi().boundaryField()[patchi];
138  const scalarField& alpha1p = alpha1.boundaryField()[patchi];
139 
140  forAll(alphaPhiCorrp, facei)
141  {
142  if (phi1p[facei] < 0)
143  {
144  alphaPhiCorrp[facei] = alpha1p[facei]*phi1p[facei];
145  }
146  }
147  }
148  }
149 
151  (
152  1.0/mesh_.time().deltaT().value(),
153  geometricOneField(),
154  phase1,
155  phi_,
156  alphaPhiCorr,
157  zeroField(),
158  zeroField(),
159  1,
160  0,
161  true
162  );
163 
164  phasei++;
165  }
166 
167  MULES::limitSum(alphaPhiCorrs);
168 
169  volScalarField sumAlpha
170  (
171  IOobject
172  (
173  "sumAlpha",
174  mesh_.time().timeName(),
175  mesh_
176  ),
177  mesh_,
178  dimensionedScalar("sumAlpha", dimless, 0)
179  );
180 
181  phasei = 0;
182 
183  forAllIter(PtrDictionary<phaseModel>, phases_, iter)
184  {
185  phaseModel& phase1 = iter();
186 
187  surfaceScalarField& alphaPhi = alphaPhiCorrs[phasei];
188  alphaPhi += upwind<scalar>(mesh_, phi_).flux(phase1);
189 
191  (
192  geometricOneField(),
193  phase1,
194  alphaPhi,
195  zeroField(),
196  zeroField()
197  );
198 
199  phase1.alphaPhi() += alphaPhi;
200 
201  Info<< phase1.name() << " volume fraction, min, max = "
202  << phase1.weightedAverage(mesh_.V()).value()
203  << ' ' << min(phase1).value()
204  << ' ' << max(phase1).value()
205  << endl;
206 
207  sumAlpha += phase1;
208 
209  phasei++;
210  }
211 
212  Info<< "Phase-sum volume fraction, min, max = "
213  << sumAlpha.weightedAverage(mesh_.V()).value()
214  << ' ' << min(sumAlpha).value()
215  << ' ' << max(sumAlpha).value()
216  << endl;
217 
218  calcAlphas();
219 }
220 
221 
222 Foam::tmp<Foam::surfaceVectorField> Foam::multiphaseSystem::nHatfv
223 (
224  const volScalarField& alpha1,
225  const volScalarField& alpha2
226 ) const
227 {
228  /*
229  // Cell gradient of alpha
230  volVectorField gradAlpha =
231  alpha2*fvc::grad(alpha1) - alpha1*fvc::grad(alpha2);
232 
233  // Interpolated face-gradient of alpha
234  surfaceVectorField gradAlphaf = fvc::interpolate(gradAlpha);
235  */
236 
237  surfaceVectorField gradAlphaf
238  (
240  - fvc::interpolate(alpha1)*fvc::interpolate(fvc::grad(alpha2))
241  );
242 
243  // Face unit interface normal
244  return gradAlphaf/(mag(gradAlphaf) + deltaN_);
245 }
246 
247 
248 Foam::tmp<Foam::surfaceScalarField> Foam::multiphaseSystem::nHatf
249 (
250  const volScalarField& alpha1,
251  const volScalarField& alpha2
252 ) const
253 {
254  // Face unit interface normal flux
255  return nHatfv(alpha1, alpha2) & mesh_.Sf();
256 }
257 
258 
259 // Correction for the boundary condition on the unit normal nHat on
260 // walls to produce the correct contact angle.
261 
262 // The dynamic contact angle is calculated from the component of the
263 // velocity on the direction of the interface, parallel to the wall.
264 
265 void Foam::multiphaseSystem::correctContactAngle
266 (
267  const phaseModel& phase1,
268  const phaseModel& phase2,
269  surfaceVectorField::Boundary& nHatb
270 ) const
271 {
272  const volScalarField::Boundary& gbf
273  = phase1.boundaryField();
274 
275  const fvBoundaryMesh& boundary = mesh_.boundary();
276 
277  forAll(boundary, patchi)
278  {
279  if (isA<alphaContactAngleFvPatchScalarField>(gbf[patchi]))
280  {
281  const alphaContactAngleFvPatchScalarField& acap =
282  refCast<const alphaContactAngleFvPatchScalarField>(gbf[patchi]);
283 
284  vectorField& nHatPatch = nHatb[patchi];
285 
286  vectorField AfHatPatch
287  (
288  mesh_.Sf().boundaryField()[patchi]
289  /mesh_.magSf().boundaryField()[patchi]
290  );
291 
293  const_iterator tp =
294  acap.thetaProps().find(interfacePair(phase1, phase2));
295 
296  if (tp == acap.thetaProps().end())
297  {
299  << "Cannot find interface " << interfacePair(phase1, phase2)
300  << "\n in table of theta properties for patch "
301  << acap.patch().name()
302  << exit(FatalError);
303  }
304 
305  bool matched = (tp.key().first() == phase1.name());
306 
307  scalar theta0 = convertToRad*tp().theta0(matched);
308  scalarField theta(boundary[patchi].size(), theta0);
309 
310  scalar uTheta = tp().uTheta();
311 
312  // Calculate the dynamic contact angle if required
313  if (uTheta > SMALL)
314  {
315  scalar thetaA = convertToRad*tp().thetaA(matched);
316  scalar thetaR = convertToRad*tp().thetaR(matched);
317 
318  // Calculated the component of the velocity parallel to the wall
319  vectorField Uwall
320  (
321  phase1.U().boundaryField()[patchi].patchInternalField()
322  - phase1.U().boundaryField()[patchi]
323  );
324  Uwall -= (AfHatPatch & Uwall)*AfHatPatch;
325 
326  // Find the direction of the interface parallel to the wall
327  vectorField nWall
328  (
329  nHatPatch - (AfHatPatch & nHatPatch)*AfHatPatch
330  );
331 
332  // Normalise nWall
333  nWall /= (mag(nWall) + SMALL);
334 
335  // Calculate Uwall resolved normal to the interface parallel to
336  // the interface
337  scalarField uwall(nWall & Uwall);
338 
339  theta += (thetaA - thetaR)*tanh(uwall/uTheta);
340  }
341 
342 
343  // Reset nHatPatch to correspond to the contact angle
344 
345  scalarField a12(nHatPatch & AfHatPatch);
346 
347  scalarField b1(cos(theta));
348 
349  scalarField b2(nHatPatch.size());
350 
351  forAll(b2, facei)
352  {
353  b2[facei] = cos(acos(a12[facei]) - theta[facei]);
354  }
355 
356  scalarField det(1.0 - a12*a12);
357 
358  scalarField a((b1 - a12*b2)/det);
359  scalarField b((b2 - a12*b1)/det);
360 
361  nHatPatch = a*AfHatPatch + b*nHatPatch;
362 
363  nHatPatch /= (mag(nHatPatch) + deltaN_.value());
364  }
365  }
366 }
367 
368 
369 Foam::tmp<Foam::volScalarField> Foam::multiphaseSystem::K
370 (
371  const phaseModel& phase1,
372  const phaseModel& phase2
373 ) const
374 {
375  tmp<surfaceVectorField> tnHatfv = nHatfv(phase1, phase2);
376 
377  correctContactAngle(phase1, phase2, tnHatfv.ref().boundaryFieldRef());
378 
379  // Simple expression for curvature
380  return -fvc::div(tnHatfv & mesh_.Sf());
381 }
382 
383 
384 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
385 
387 (
388  const volVectorField& U,
389  const surfaceScalarField& phi
390 )
391 :
393  (
394  IOobject
395  (
396  "transportProperties",
397  U.time().constant(),
398  U.db(),
401  )
402  ),
403 
404  phases_(lookup("phases"), phaseModel::iNew(U.mesh())),
405 
406  mesh_(U.mesh()),
407  phi_(phi),
408 
409  alphas_
410  (
411  IOobject
412  (
413  "alphas",
414  mesh_.time().timeName(),
415  mesh_,
416  IOobject::NO_READ,
418  ),
419  mesh_,
420  dimensionedScalar("alphas", dimless, 0.0)
421  ),
422 
423  sigmas_(lookup("sigmas")),
424  dimSigma_(1, 0, -2, 0, 0),
425  cAlphas_(lookup("interfaceCompression")),
426  Cvms_(lookup("virtualMass")),
427  deltaN_
428  (
429  "deltaN",
430  1e-8/pow(average(mesh_.V()), 1.0/3.0)
431  )
432 {
433  calcAlphas();
434  alphas_.write();
435 
436  interfaceDictTable dragModelsDict(lookup("drag"));
437 
438  forAllConstIter(interfaceDictTable, dragModelsDict, iter)
439  {
440  dragModels_.insert
441  (
442  iter.key(),
444  (
445  iter(),
446  *phases_.lookup(iter.key().first()),
447  *phases_.lookup(iter.key().second())
448  ).ptr()
449  );
450  }
451 
452  forAllConstIter(PtrDictionary<phaseModel>, phases_, iter1)
453  {
454  const phaseModel& phase1 = iter1();
455 
456  forAllConstIter(PtrDictionary<phaseModel>, phases_, iter2)
457  {
458  const phaseModel& phase2 = iter2();
459 
460  if (&phase2 != &phase1)
461  {
462  scalarCoeffSymmTable::const_iterator sigma
463  (
464  sigmas_.find(interfacePair(phase1, phase2))
465  );
466 
467  if (sigma != sigmas_.end())
468  {
469  scalarCoeffSymmTable::const_iterator cAlpha
470  (
471  cAlphas_.find(interfacePair(phase1, phase2))
472  );
473 
474  if (cAlpha == cAlphas_.end())
475  {
477  << "Compression coefficient not specified for "
478  "phase pair ("
479  << phase1.name() << ' ' << phase2.name()
480  << ") for which a surface tension "
481  "coefficient is specified"
482  << endl;
483  }
484  }
485  }
486  }
487  }
488 }
489 
490 
491 // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
492 
494 {
495  PtrDictionary<phaseModel>::const_iterator iter = phases_.begin();
496 
497  tmp<volScalarField> trho = iter()*iter().rho();
498  volScalarField& rho = trho.ref();
499 
500  for (++iter; iter != phases_.end(); ++iter)
501  {
502  rho += iter()*iter().rho();
503  }
504 
505  return trho;
506 }
507 
508 
510 Foam::multiphaseSystem::rho(const label patchi) const
511 {
512  PtrDictionary<phaseModel>::const_iterator iter = phases_.begin();
513 
514  tmp<scalarField> trho = iter().boundaryField()[patchi]*iter().rho().value();
515  scalarField& rho = trho.ref();
516 
517  for (++iter; iter != phases_.end(); ++iter)
518  {
519  rho += iter().boundaryField()[patchi]*iter().rho().value();
520  }
521 
522  return trho;
523 }
524 
525 
527 {
528  PtrDictionary<phaseModel>::const_iterator iter = phases_.begin();
529 
530  tmp<volScalarField> tmu = iter()*(iter().rho()*iter().nu());
531  volScalarField& mu = tmu.ref();
532 
533  for (++iter; iter != phases_.end(); ++iter)
534  {
535  mu += iter()*(iter().rho()*iter().nu());
536  }
537 
538  return tmu/rho();
539 }
540 
541 
543 Foam::multiphaseSystem::nu(const label patchi) const
544 {
545  PtrDictionary<phaseModel>::const_iterator iter = phases_.begin();
546 
547  tmp<scalarField> tmu =
548  iter().boundaryField()[patchi]
549  *(iter().rho().value()*iter().nu().value());
550  scalarField& mu = tmu.ref();
551 
552  for (++iter; iter != phases_.end(); ++iter)
553  {
554  mu +=
555  iter().boundaryField()[patchi]
556  *(iter().rho().value()*iter().nu().value());
557  }
558 
559  return tmu/rho(patchi);
560 }
561 
562 
564 (
565  const phaseModel& phase
566 ) const
567 {
568  tmp<volScalarField> tCvm
569  (
570  new volScalarField
571  (
572  IOobject
573  (
574  "Cvm",
575  mesh_.time().timeName(),
576  mesh_
577  ),
578  mesh_,
580  (
581  "Cvm",
582  dimensionSet(1, -3, 0, 0, 0),
583  0
584  )
585  )
586  );
587 
588  forAllConstIter(PtrDictionary<phaseModel>, phases_, iter)
589  {
590  const phaseModel& phase2 = iter();
591 
592  if (&phase2 != &phase)
593  {
594  scalarCoeffTable::const_iterator Cvm
595  (
596  Cvms_.find(interfacePair(phase, phase2))
597  );
598 
599  if (Cvm != Cvms_.end())
600  {
601  tCvm.ref() += Cvm()*phase2.rho()*phase2;
602  }
603  else
604  {
605  Cvm = Cvms_.find(interfacePair(phase2, phase));
606 
607  if (Cvm != Cvms_.end())
608  {
609  tCvm.ref() += Cvm()*phase.rho()*phase2;
610  }
611  }
612  }
613  }
614 
615  return tCvm;
616 }
617 
618 
620 (
621  const phaseModel& phase
622 ) const
623 {
624  tmp<volVectorField> tSvm
625  (
626  new volVectorField
627  (
628  IOobject
629  (
630  "Svm",
631  mesh_.time().timeName(),
632  mesh_
633  ),
634  mesh_,
636  (
637  "Svm",
638  dimensionSet(1, -2, -2, 0, 0),
639  Zero
640  )
641  )
642  );
643 
644  forAllConstIter(PtrDictionary<phaseModel>, phases_, iter)
645  {
646  const phaseModel& phase2 = iter();
647 
648  if (&phase2 != &phase)
649  {
650  scalarCoeffTable::const_iterator Cvm
651  (
652  Cvms_.find(interfacePair(phase, phase2))
653  );
654 
655  if (Cvm != Cvms_.end())
656  {
657  tSvm.ref() += Cvm()*phase2.rho()*phase2*phase2.DDtU();
658  }
659  else
660  {
661  Cvm = Cvms_.find(interfacePair(phase2, phase));
662 
663  if (Cvm != Cvms_.end())
664  {
665  tSvm.ref() += Cvm()*phase.rho()*phase2*phase2.DDtU();
666  }
667  }
668  }
669  }
670 
671  volVectorField::Boundary& SvmBf =
672  tSvm.ref().boundaryFieldRef();
673 
674  // Remove virtual mass at fixed-flux boundaries
675  forAll(phase.phi().boundaryField(), patchi)
676  {
677  if
678  (
679  isA<fixedValueFvsPatchScalarField>
680  (
681  phase.phi().boundaryField()[patchi]
682  )
683  )
684  {
685  SvmBf[patchi] = Zero;
686  }
687  }
688 
689  return tSvm;
690 }
691 
692 
695 {
696  autoPtr<dragCoeffFields> dragCoeffsPtr(new dragCoeffFields);
697 
698  forAllConstIter(dragModelTable, dragModels_, iter)
699  {
700  const dragModel& dm = *iter();
701 
702  volScalarField* Kptr =
703  (
704  max
705  (
706  //fvc::average(dm.phase1()*dm.phase2()),
707  //fvc::average(dm.phase1())*fvc::average(dm.phase2()),
708  dm.phase1()*dm.phase2(),
709  dm.residualPhaseFraction()
710  )
711  *dm.K
712  (
713  max
714  (
715  mag(dm.phase1().U() - dm.phase2().U()),
716  dm.residualSlip()
717  )
718  )
719  ).ptr();
720 
721  volScalarField::Boundary& Kbf = Kptr->boundaryFieldRef();
722 
723  // Remove drag at fixed-flux boundaries
724  forAll(dm.phase1().phi().boundaryField(), patchi)
725  {
726  if
727  (
728  isA<fixedValueFvsPatchScalarField>
729  (
730  dm.phase1().phi().boundaryField()[patchi]
731  )
732  )
733  {
734  Kbf[patchi] = 0.0;
735  }
736  }
737 
738  dragCoeffsPtr().insert(iter.key(), Kptr);
739  }
740 
741  return dragCoeffsPtr;
742 }
743 
744 
746 (
747  const phaseModel& phase,
748  const dragCoeffFields& dragCoeffs
749 ) const
750 {
751  tmp<volScalarField> tdragCoeff
752  (
753  new volScalarField
754  (
755  IOobject
756  (
757  "dragCoeff",
758  mesh_.time().timeName(),
759  mesh_
760  ),
761  mesh_,
763  (
764  "dragCoeff",
765  dimensionSet(1, -3, -1, 0, 0),
766  0
767  )
768  )
769  );
770 
771  dragModelTable::const_iterator dmIter = dragModels_.begin();
772  dragCoeffFields::const_iterator dcIter = dragCoeffs.begin();
773  for
774  (
775  ;
776  dmIter != dragModels_.end() && dcIter != dragCoeffs.end();
777  ++dmIter, ++dcIter
778  )
779  {
780  if
781  (
782  &phase == &dmIter()->phase1()
783  || &phase == &dmIter()->phase2()
784  )
785  {
786  tdragCoeff.ref() += *dcIter();
787  }
788  }
789 
790  return tdragCoeff;
791 }
792 
793 
795 (
796  const phaseModel& phase1
797 ) const
798 {
799  tmp<surfaceScalarField> tSurfaceTension
800  (
802  (
803  IOobject
804  (
805  "surfaceTension",
806  mesh_.time().timeName(),
807  mesh_
808  ),
809  mesh_,
811  (
812  "surfaceTension",
813  dimensionSet(1, -2, -2, 0, 0),
814  0
815  )
816  )
817  );
818 
819  forAllConstIter(PtrDictionary<phaseModel>, phases_, iter)
820  {
821  const phaseModel& phase2 = iter();
822 
823  if (&phase2 != &phase1)
824  {
825  scalarCoeffSymmTable::const_iterator sigma
826  (
827  sigmas_.find(interfacePair(phase1, phase2))
828  );
829 
830  if (sigma != sigmas_.end())
831  {
832  tSurfaceTension.ref() +=
833  dimensionedScalar("sigma", dimSigma_, sigma())
834  *fvc::interpolate(K(phase1, phase2))*
835  (
836  fvc::interpolate(phase2)*fvc::snGrad(phase1)
837  - fvc::interpolate(phase1)*fvc::snGrad(phase2)
838  );
839  }
840  }
841  }
842 
843  return tSurfaceTension;
844 }
845 
846 
849 {
850  tmp<volScalarField> tnearInt
851  (
852  new volScalarField
853  (
854  IOobject
855  (
856  "nearInterface",
857  mesh_.time().timeName(),
858  mesh_
859  ),
860  mesh_,
861  dimensionedScalar("nearInterface", dimless, 0.0)
862  )
863  );
864 
865  forAllConstIter(PtrDictionary<phaseModel>, phases_, iter)
866  {
867  tnearInt.ref() = max(tnearInt(), pos(iter() - 0.01)*pos(0.99 - iter()));
868  }
869 
870  return tnearInt;
871 }
872 
873 
875 {
876  forAllIter(PtrDictionary<phaseModel>, phases_, iter)
877  {
878  iter().correct();
879  }
880 
881  const Time& runTime = mesh_.time();
882 
883  const dictionary& alphaControls = mesh_.solverDict("alpha");
884  label nAlphaSubCycles(readLabel(alphaControls.lookup("nAlphaSubCycles")));
885 
886  if (nAlphaSubCycles > 1)
887  {
888  dimensionedScalar totalDeltaT = runTime.deltaT();
889 
890  PtrList<volScalarField> alpha0s(phases_.size());
891  PtrList<surfaceScalarField> alphaPhiSums(phases_.size());
892 
893  int phasei = 0;
894  forAllIter(PtrDictionary<phaseModel>, phases_, iter)
895  {
896  phaseModel& phase = iter();
898 
899  alpha0s.set
900  (
901  phasei,
902  new volScalarField(alpha.oldTime())
903  );
904 
905  alphaPhiSums.set
906  (
907  phasei,
909  (
910  IOobject
911  (
912  "phiSum" + alpha.name(),
913  runTime.timeName(),
914  mesh_
915  ),
916  mesh_,
917  dimensionedScalar("0", dimensionSet(0, 3, -1, 0, 0), 0)
918  )
919  );
920 
921  phasei++;
922  }
923 
924  for
925  (
926  subCycleTime alphaSubCycle
927  (
928  const_cast<Time&>(runTime),
930  );
931  !(++alphaSubCycle).end();
932  )
933  {
934  solveAlphas();
935 
936  int phasei = 0;
937  forAllIter(PtrDictionary<phaseModel>, phases_, iter)
938  {
939  alphaPhiSums[phasei] += iter().alphaPhi()/nAlphaSubCycles;
940  phasei++;
941  }
942  }
943 
944  phasei = 0;
945  forAllIter(PtrDictionary<phaseModel>, phases_, iter)
946  {
947  phaseModel& phase = iter();
948  volScalarField& alpha = phase;
949 
950  phase.alphaPhi() = alphaPhiSums[phasei];
951 
952  // Correct the time index of the field
953  // to correspond to the global time
954  alpha.timeIndex() = runTime.timeIndex();
955 
956  // Reset the old-time field value
957  alpha.oldTime() = alpha0s[phasei];
958  alpha.oldTime().timeIndex() = runTime.timeIndex();
959 
960  phasei++;
961  }
962  }
963  else
964  {
965  solveAlphas();
966  }
967 }
968 
969 
971 {
972  if (regIOobject::read())
973  {
974  bool readOK = true;
975 
976  PtrList<entry> phaseData(lookup("phases"));
977  label phasei = 0;
978 
979  forAllIter(PtrDictionary<phaseModel>, phases_, iter)
980  {
981  readOK &= iter().read(phaseData[phasei++].dict());
982  }
983 
984  lookup("sigmas") >> sigmas_;
985  lookup("interfaceCompression") >> cAlphas_;
986  lookup("virtualMass") >> Cvms_;
987 
988  return readOK;
989  }
990  else
991  {
992  return false;
993  }
994 }
995 
996 
997 // ************************************************************************* //
tmp< volScalarField > nu() const
Return the mixture laminar viscosity.
fvsPatchField< scalar > fvsPatchScalarField
dimensionedScalar tanh(const dimensionedScalar &ds)
dictionary dict
tmp< GeometricField< typename outerProduct< vector, Type >::type, fvPatchField, volMesh >> grad(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcGrad.C:52
dimensionedScalar acos(const dimensionedScalar &ds)
#define forAll(list, i)
Loop across all elements in list.
Definition: UList.H:428
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
tmp< surfaceScalarField > surfaceTension(const phaseModel &phase) const
const objectRegistry & db() const
Return the local objectRegistry.
Definition: IOobject.C:221
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
label phasei
Definition: pEqn.H:27
static iteratorEnd end()
iteratorEnd set to beyond the end of any HashTable
Definition: HashTable.H:106
const double e
Elementary charge.
Definition: doubleFloat.H:78
virtual bool read()
Read object.
error FatalError
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
tmp< volVectorField > U() const
Return the mixture velocity.
tmp< GeometricField< Type, fvPatchField, volMesh > > div(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcDiv.C:47
#define forAllIter(Container, container, iter)
Iterate across all elements in the container object of type.
Definition: UList.H:453
IOobject(const word &name, const fileName &instance, const objectRegistry &registry, readOption r=NO_READ, writeOption w=NO_WRITE, bool registerObject=true)
Construct from name, instance, registry, io options.
Definition: IOobject.C:114
GeometricField< vector, fvsPatchField, surfaceMesh > surfaceVectorField
dimensioned< vector > dimensionedVector
Dimensioned vector obtained from generic dimensioned type.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
tmp< volScalarField > trho
const dimensionedScalar sigma
Stefan-Boltzmann constant: default SI units: [W/m2/K4].
autoPtr< dragCoeffFields > dragCoeffs() const
Return the drag coefficients for all of the interfaces.
autoPtr< multiphaseSystem::dragCoeffFields > dragCoeffs(fluid.dragCoeffs())
const Type & value() const
Return const reference to value.
tmp< volScalarField > Cvm(const phaseModel &phase) const
Return the virtual-mass coefficient for the given phase.
Calculate the snGrad of the given volField.
dimensionedScalar det(const dimensionedSphericalTensor &dt)
label size() const
Return number of elements in list.
Definition: DLListBaseI.H:77
GeometricField< vector, fvPatchField, volMesh > volVectorField
Definition: volFieldsFwd.H:55
CGAL::Exact_predicates_exact_constructions_kernel K
Area-weighted average a surfaceField creating a volField.
surfaceScalarField phir(phic *interface.nHatf())
iterator find(const Key &)
Find and return an iterator set at the hashedEntry.
Definition: HashTable.C:138
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:52
const surfaceScalarField & phi() const
Constant access the mixture flux.
Definition: phaseSystemI.H:55
tmp< volScalarField > dragCoeff(const phaseModel &phase, const dragCoeffFields &dragCoeffs) const
Return the sum of the drag coefficients for the given phase.
dimensionedScalar pos(const dimensionedScalar &ds)
surfaceScalarField alphaPhi(IOobject("alphaPhi", runTime.timeName(), mesh, IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE), phi *fvc::interpolate(alpha1))
stressControl lookup("compactNormalStress") >> compactNormalStress
static autoPtr< dragModel > New(const dictionary &interfaceDict, const phaseModel &phase1, const phaseModel &phase2)
dimensionedScalar cos(const dimensionedScalar &ds)
const dictionary & alphaControls
Definition: alphaControls.H:1
void limitSum(UPtrList< scalarField > &phiPsiCorrs)
Definition: MULES.C:51
const fvMesh & mesh() const
Constant access the mesh.
Definition: phaseSystemI.H:28
const dimensionedScalar b
Wien displacement law constant: default SI units: [m.K].
Definition: createFields.H:27
Calculate the gradient of the given field.
void explicitSolve(const RdeltaTType &rDeltaT, const RhoType &rho, volScalarField &psi, const surfaceScalarField &phiPsi, const SpType &Sp, const SuType &Su)
const Time & time() const
Return time.
Definition: IOobject.C:227
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Calculate the face-flux of the given field.
friend class const_iterator
Declare friendship with the const_iterator.
Definition: HashTable.H:192
label timeIndex() const
Return the time index of the field.
word timeName
Definition: getTimeIndex.H:3
static const zero Zero
Definition: zero.H:91
dimensioned< Type > average(const DimensionedField< Type, GeoMesh > &df)
faceListList boundary(nPatches)
void solve()
Solve for the mixture phase-fractions.
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
Definition: pEqn.H:29
label readLabel(Istream &is)
Definition: label.H:64
phaseModel & phase1
tmp< volScalarField > rho() const
Return the mixture density.
Calculate the divergence of the given field.
tmp< volVectorField > Svm(const phaseModel &phase) const
Return the virtual-mass source for the given phase.
bool read()
Read base transportProperties dictionary.
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
static tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > interpolate(const GeometricField< Type, fvPatchField, volMesh > &tvf, const surfaceScalarField &faceFlux, Istream &schemeData)
Interpolate field onto faces using scheme given by Istream.
const dimensionedScalar mu
Atomic mass unit.
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
IOdictionary(const IOobject &)
Construct given an IOobject.
Definition: IOdictionary.C:45
label patchi
const dimensionSet dimless(0, 0, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:47
#define WarningInFunction
Report a warning using Foam::Warning.
Boundary & boundaryFieldRef()
Return a reference to the boundary field.
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
const phaseModel & phase() const
Return the phase.
messageStream Info
dimensioned< scalar > mag(const dimensioned< Type > &)
Field< vector > vectorField
Specialisation of Field<T> for vector.
MULES: Multidimensional universal limiter for explicit solution.
const surfaceScalarField & alphaPhi() const
Definition: phaseModel.H:200
An auto-pointer similar to the STL auto_ptr but with automatic casting to a reference to the type and...
Definition: PtrList.H:53
tmp< volScalarField > nearInterface() const
Indicator of the proximity of the interface.
void limit(const RdeltaTType &rDeltaT, const RhoType &rho, const volScalarField &psi, const surfaceScalarField &phi, surfaceScalarField &phiPsi, const SpType &Sp, const SuType &Su, const scalar psiMax, const scalar psiMin, const bool returnCorr)
phaseModel & phase2
A class for managing temporary objects.
Definition: PtrList.H:54
tmp< surfaceScalarField > flux(const volVectorField &vvf)
Return the face-flux field obtained from the given volVectorField.
Definition: fvcFlux.C:32
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
multiphaseSystem(const volVectorField &U, const surfaceScalarField &phi)
Construct from components.
const dimensionedScalar alpha
Fine-structure constant: default SI units: [].
tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > snGrad(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvcSnGrad.C:45
label nAlphaSubCycles(readLabel(alphaControls.lookup("nAlphaSubCycles")))
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:451