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