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-2015 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "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  alphas_.correctBoundaryConditions();
59 }
60 
61 
62 void Foam::multiphaseSystem::solveAlphas()
63 {
64  PtrList<surfaceScalarField> alphaPhiCorrs(phases_.size());
65  int phasei = 0;
66 
67  forAllIter(PtrDictionary<phaseModel>, phases_, iter)
68  {
69  phaseModel& phase1 = iter();
70  volScalarField& alpha1 = phase1;
71 
72  phase1.alphaPhi() =
73  dimensionedScalar("0", dimensionSet(0, 3, -1, 0, 0), 0);
74 
75  alphaPhiCorrs.set
76  (
77  phasei,
79  (
80  "phi" + alpha1.name() + "Corr",
81  fvc::flux
82  (
83  phi_,
84  phase1,
85  "div(phi," + alpha1.name() + ')'
86  )
87  )
88  );
89 
90  surfaceScalarField& alphaPhiCorr = alphaPhiCorrs[phasei];
91 
92  forAllIter(PtrDictionary<phaseModel>, phases_, iter2)
93  {
94  phaseModel& phase2 = iter2();
95  volScalarField& alpha2 = phase2;
96 
97  if (&phase2 == &phase1) continue;
98 
99  surfaceScalarField phir(phase1.phi() - phase2.phi());
100 
102  (
103  cAlphas_.find(interfacePair(phase1, phase2))
104  );
105 
106  if (cAlpha != cAlphas_.end())
107  {
109  (
110  (mag(phi_) + mag(phir))/mesh_.magSf()
111  );
112 
113  phir += min(cAlpha()*phic, max(phic))*nHatf(phase1, phase2);
114  }
115 
116  word phirScheme
117  (
118  "div(phir," + alpha2.name() + ',' + alpha1.name() + ')'
119  );
120 
121  alphaPhiCorr += fvc::flux
122  (
123  -fvc::flux(-phir, phase2, phirScheme),
124  phase1,
125  phirScheme
126  );
127  }
128 
129  // Ensure that the flux at inflow BCs is preserved
130  forAll(alphaPhiCorr.boundaryField(), patchi)
131  {
132  fvsPatchScalarField& alphaPhiCorrp =
133  alphaPhiCorr.boundaryField()[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::GeometricBoundaryField& nHatb
270 ) const
271 {
272  const volScalarField::GeometricBoundaryField& 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 
292  alphaContactAngleFvPatchScalarField::thetaPropsTable::
293  const_iterator tp =
294  acap.thetaProps().find(interfacePair(phase1, phase2));
295 
296  if (tp == acap.thetaProps().end())
297  {
299  (
300  "multiphaseSystem::correctContactAngle"
301  "(const phaseModel& phase1, const phaseModel& phase2, "
302  "fvPatchVectorFieldField& nHatb) const"
303  ) << "Cannot find interface " << interfacePair(phase1, phase2)
304  << "\n in table of theta properties for patch "
305  << acap.patch().name()
306  << exit(FatalError);
307  }
308 
309  bool matched = (tp.key().first() == phase1.name());
310 
311  scalar theta0 = convertToRad*tp().theta0(matched);
312  scalarField theta(boundary[patchi].size(), theta0);
313 
314  scalar uTheta = tp().uTheta();
315 
316  // Calculate the dynamic contact angle if required
317  if (uTheta > SMALL)
318  {
319  scalar thetaA = convertToRad*tp().thetaA(matched);
320  scalar thetaR = convertToRad*tp().thetaR(matched);
321 
322  // Calculated the component of the velocity parallel to the wall
323  vectorField Uwall
324  (
325  phase1.U().boundaryField()[patchi].patchInternalField()
326  - phase1.U().boundaryField()[patchi]
327  );
328  Uwall -= (AfHatPatch & Uwall)*AfHatPatch;
329 
330  // Find the direction of the interface parallel to the wall
331  vectorField nWall
332  (
333  nHatPatch - (AfHatPatch & nHatPatch)*AfHatPatch
334  );
335 
336  // Normalise nWall
337  nWall /= (mag(nWall) + SMALL);
338 
339  // Calculate Uwall resolved normal to the interface parallel to
340  // the interface
341  scalarField uwall(nWall & Uwall);
342 
343  theta += (thetaA - thetaR)*tanh(uwall/uTheta);
344  }
345 
346 
347  // Reset nHatPatch to correspond to the contact angle
348 
349  scalarField a12(nHatPatch & AfHatPatch);
350 
351  scalarField b1(cos(theta));
352 
353  scalarField b2(nHatPatch.size());
354 
355  forAll(b2, facei)
356  {
357  b2[facei] = cos(acos(a12[facei]) - theta[facei]);
358  }
359 
360  scalarField det(1.0 - a12*a12);
361 
362  scalarField a((b1 - a12*b2)/det);
363  scalarField b((b2 - a12*b1)/det);
364 
365  nHatPatch = a*AfHatPatch + b*nHatPatch;
366 
367  nHatPatch /= (mag(nHatPatch) + deltaN_.value());
368  }
369  }
370 }
371 
372 
373 Foam::tmp<Foam::volScalarField> Foam::multiphaseSystem::K
374 (
375  const phaseModel& phase1,
376  const phaseModel& phase2
377 ) const
378 {
379  tmp<surfaceVectorField> tnHatfv = nHatfv(phase1, phase2);
380 
381  correctContactAngle(phase1, phase2, tnHatfv().boundaryField());
382 
383  // Simple expression for curvature
384  return -fvc::div(tnHatfv & mesh_.Sf());
385 }
386 
387 
388 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
389 
391 (
392  const volVectorField& U,
393  const surfaceScalarField& phi
394 )
395 :
397  (
398  IOobject
399  (
400  "transportProperties",
401  U.time().constant(),
402  U.db(),
405  )
406  ),
407 
408  phases_(lookup("phases"), phaseModel::iNew(U.mesh())),
409 
410  mesh_(U.mesh()),
411  phi_(phi),
412 
413  alphas_
414  (
415  IOobject
416  (
417  "alphas",
418  mesh_.time().timeName(),
419  mesh_,
420  IOobject::NO_READ,
422  ),
423  mesh_,
424  dimensionedScalar("alphas", dimless, 0.0),
425  zeroGradientFvPatchScalarField::typeName
426  ),
427 
428  sigmas_(lookup("sigmas")),
429  dimSigma_(1, 0, -2, 0, 0),
430  cAlphas_(lookup("interfaceCompression")),
431  Cvms_(lookup("virtualMass")),
432  deltaN_
433  (
434  "deltaN",
435  1e-8/pow(average(mesh_.V()), 1.0/3.0)
436  )
437 {
438  calcAlphas();
439  alphas_.write();
440 
441  interfaceDictTable dragModelsDict(lookup("drag"));
442 
443  forAllConstIter(interfaceDictTable, dragModelsDict, iter)
444  {
445  dragModels_.insert
446  (
447  iter.key(),
449  (
450  iter(),
451  *phases_.lookup(iter.key().first()),
452  *phases_.lookup(iter.key().second())
453  ).ptr()
454  );
455  }
456 
457  forAllConstIter(PtrDictionary<phaseModel>, phases_, iter1)
458  {
459  const phaseModel& phase1 = iter1();
460 
461  forAllConstIter(PtrDictionary<phaseModel>, phases_, iter2)
462  {
463  const phaseModel& phase2 = iter2();
464 
465  if (&phase2 != &phase1)
466  {
467  scalarCoeffSymmTable::const_iterator sigma
468  (
469  sigmas_.find(interfacePair(phase1, phase2))
470  );
471 
472  if (sigma != sigmas_.end())
473  {
474  scalarCoeffSymmTable::const_iterator cAlpha
475  (
476  cAlphas_.find(interfacePair(phase1, phase2))
477  );
478 
479  if (cAlpha == cAlphas_.end())
480  {
481  WarningIn
482  (
483  "multiphaseSystem::multiphaseSystem"
484  "(const volVectorField& U,"
485  "const surfaceScalarField& phi)"
486  ) << "Compression coefficient not specified for "
487  "phase pair ("
488  << phase1.name() << ' ' << phase2.name()
489  << ") for which a surface tension "
490  "coefficient is specified"
491  << endl;
492  }
493  }
494  }
495  }
496  }
497 }
498 
499 
500 // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
501 
503 {
504  PtrDictionary<phaseModel>::const_iterator iter = phases_.begin();
505 
506  tmp<volScalarField> trho = iter()*iter().rho();
507 
508  for (++iter; iter != phases_.end(); ++iter)
509  {
510  trho() += iter()*iter().rho();
511  }
512 
513  return trho;
514 }
515 
516 
518 Foam::multiphaseSystem::rho(const label patchi) const
519 {
520  PtrDictionary<phaseModel>::const_iterator iter = phases_.begin();
521 
522  tmp<scalarField> trho = iter().boundaryField()[patchi]*iter().rho().value();
523 
524  for (++iter; iter != phases_.end(); ++iter)
525  {
526  trho() += iter().boundaryField()[patchi]*iter().rho().value();
527  }
528 
529  return trho;
530 }
531 
532 
534 {
535  PtrDictionary<phaseModel>::const_iterator iter = phases_.begin();
536 
537  tmp<volScalarField> tmu = iter()*(iter().rho()*iter().nu());
538 
539  for (++iter; iter != phases_.end(); ++iter)
540  {
541  tmu() += iter()*(iter().rho()*iter().nu());
542  }
543 
544  return tmu/rho();
545 }
546 
547 
549 Foam::multiphaseSystem::nu(const label patchi) const
550 {
551  PtrDictionary<phaseModel>::const_iterator iter = phases_.begin();
552 
553  tmp<scalarField> tmu =
554  iter().boundaryField()[patchi]
555  *(iter().rho().value()*iter().nu().value());
556 
557  for (++iter; iter != phases_.end(); ++iter)
558  {
559  tmu() +=
560  iter().boundaryField()[patchi]
561  *(iter().rho().value()*iter().nu().value());
562  }
563 
564  return tmu/rho(patchi);
565 }
566 
567 
569 (
570  const phaseModel& phase
571 ) const
572 {
573  tmp<volScalarField> tCvm
574  (
575  new volScalarField
576  (
577  IOobject
578  (
579  "Cvm",
580  mesh_.time().timeName(),
581  mesh_
582  ),
583  mesh_,
585  (
586  "Cvm",
587  dimensionSet(1, -3, 0, 0, 0),
588  0
589  )
590  )
591  );
592 
593  forAllConstIter(PtrDictionary<phaseModel>, phases_, iter)
594  {
595  const phaseModel& phase2 = iter();
596 
597  if (&phase2 != &phase)
598  {
599  scalarCoeffTable::const_iterator Cvm
600  (
601  Cvms_.find(interfacePair(phase, phase2))
602  );
603 
604  if (Cvm != Cvms_.end())
605  {
606  tCvm() += Cvm()*phase2.rho()*phase2;
607  }
608  else
609  {
610  Cvm = Cvms_.find(interfacePair(phase2, phase));
611 
612  if (Cvm != Cvms_.end())
613  {
614  tCvm() += Cvm()*phase.rho()*phase2;
615  }
616  }
617  }
618  }
619 
620  return tCvm;
621 }
622 
623 
625 (
626  const phaseModel& phase
627 ) const
628 {
629  tmp<volVectorField> tSvm
630  (
631  new volVectorField
632  (
633  IOobject
634  (
635  "Svm",
636  mesh_.time().timeName(),
637  mesh_
638  ),
639  mesh_,
641  (
642  "Svm",
643  dimensionSet(1, -2, -2, 0, 0),
645  )
646  )
647  );
648 
649  forAllConstIter(PtrDictionary<phaseModel>, phases_, iter)
650  {
651  const phaseModel& phase2 = iter();
652 
653  if (&phase2 != &phase)
654  {
655  scalarCoeffTable::const_iterator Cvm
656  (
657  Cvms_.find(interfacePair(phase, phase2))
658  );
659 
660  if (Cvm != Cvms_.end())
661  {
662  tSvm() += Cvm()*phase2.rho()*phase2*phase2.DDtU();
663  }
664  else
665  {
666  Cvm = Cvms_.find(interfacePair(phase2, phase));
667 
668  if (Cvm != Cvms_.end())
669  {
670  tSvm() += Cvm()*phase.rho()*phase2*phase2.DDtU();
671  }
672  }
673  }
674  }
675 
676  // Remove virtual mass at fixed-flux boundaries
677  forAll(phase.phi().boundaryField(), patchi)
678  {
679  if
680  (
681  isA<fixedValueFvsPatchScalarField>
682  (
683  phase.phi().boundaryField()[patchi]
684  )
685  )
686  {
687  tSvm().boundaryField()[patchi] = vector::zero;
688  }
689  }
690 
691  return tSvm;
692 }
693 
694 
697 {
698  autoPtr<dragCoeffFields> dragCoeffsPtr(new dragCoeffFields);
699 
700  forAllConstIter(dragModelTable, dragModels_, iter)
701  {
702  const dragModel& dm = *iter();
703 
704  volScalarField* Kptr =
705  (
706  max
707  (
708  //fvc::average(dm.phase1()*dm.phase2()),
709  //fvc::average(dm.phase1())*fvc::average(dm.phase2()),
710  dm.phase1()*dm.phase2(),
711  dm.residualPhaseFraction()
712  )
713  *dm.K
714  (
715  max
716  (
717  mag(dm.phase1().U() - dm.phase2().U()),
718  dm.residualSlip()
719  )
720  )
721  ).ptr();
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  Kptr->boundaryField()[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() += *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() +=
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() = 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< GeometricField< Type, fvsPatchField, surfaceMesh > > interpolate(const GeometricField< Type, fvPatchField, volMesh > &vf, const surfaceScalarField &faceFlux, Istream &schemeData)
phaseModel & phase1
Definition: createFields.H:12
dimensioned< Type > average(const DimensionedField< Type, GeoMesh > &df)
void solve()
Solve for the mixture phase-fractions.
const phaseModel & phase() const
Return the phase.
void explicitSolve(const RdeltaTType &rDeltaT, const RhoType &rho, volScalarField &psi, const surfaceScalarField &phiPsi, const SpType &Sp, const SuType &Su)
phaseModel & phase2
Definition: createFields.H:13
faceListList boundary(nPatches)
Calculate the snGrad of the given volField.
tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > flux(const surfaceScalarField &phi, const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvcFlux.C:45
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
const dimensionedScalar b
Wien displacement law constant: default SI units: [m.K].
Definition: createFields.H:28
dimensioned< scalar > mag(const dimensioned< Type > &)
label phasei
Definition: pEqn.H:37
CGAL::Exact_predicates_exact_constructions_kernel K
const dimensionedScalar sigma
Stefan-Boltzmann constant: default SI units: [W/m2/K4].
surfaceScalarField phir(phic *interface.nHatf())
const surfaceScalarField & phi() const
Constant access the mixture flux.
Definition: phaseSystemI.H:55
static autoPtr< dragModel > New(const dictionary &interfaceDict, const phaseModel &phase1, const phaseModel &phase2)
autoPtr< multiphaseSystem::dragCoeffFields > dragCoeffs(fluid.dragCoeffs())
#define forAllIter(Container, container, iter)
Definition: UList.H:440
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
const surfaceScalarField & alphaPhi() const
Definition: phaseModel.H:200
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
Surface Interpolation.
friend class const_iterator
Declare friendship with the const_iterator.
Definition: HashTable.H:192
tmp< volVectorField > Svm(const phaseModel &phase) const
Return the virtual-mass source for the given phase.
messageStream Info
multiphaseSystem(const volVectorField &U, const surfaceScalarField &phi)
Construct from components.
tmp< volVectorField > U() const
Return the mixture velocity.
label size() const
Return number of elements in list.
Definition: DLListBaseI.H:77
const dictionary & alphaControls
Definition: alphaControls.H:1
MULES: Multidimensional universal limiter for explicit solution.
fvsPatchField< scalar > fvsPatchScalarField
label readLabel(Istream &is)
Definition: label.H:64
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Calculate the gradient of the given field.
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
Definition: pEqn.H:39
dictionary dict
const double e
Elementary charge.
Definition: doubleFloat.H:78
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:251
stressControl lookup("compactNormalStress") >> compactNormalStress
#define WarningIn(functionName)
Report a warning using Foam::Warning.
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:52
dimensionedScalar det(const dimensionedSphericalTensor &dt)
tmp< volScalarField > nu() const
Return the mixture laminar viscosity.
void limitSum(UPtrList< scalarField > &phiPsiCorrs)
Definition: MULES.C:51
const dimensionedScalar alpha
Fine-structure constant: default SI units: [].
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
Field< vector > vectorField
Specialisation of Field<T> for vector.
#define forAll(list, i)
Definition: UList.H:421
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
bool read()
Read base transportProperties dictionary.
tmp< volScalarField > nearInterface() const
Indicator of the proximity of the interface.
dimensionedScalar cos(const dimensionedScalar &ds)
label patchi
virtual bool read()
Read object.
dimensioned< vector > dimensionedVector
Dimensioned vector obtained from generic dimensioned type.
tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > snGrad(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvcSnGrad.C:45
surfaceScalarField alphaPhi(IOobject( "alphaPhi", runTime.timeName(), mesh, IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE ), phi *fvc::interpolate(alpha1))
dimensionedScalar acos(const dimensionedScalar &ds)
dimensionedScalar tanh(const dimensionedScalar &ds)
tmp< volScalarField > dragCoeff(const phaseModel &phase, const dragCoeffFields &dragCoeffs) const
Return the sum of the drag coefficients for the given phase.
#define FatalErrorIn(functionName)
Report an error message using Foam::FatalError.
Definition: error.H:314
tmp< surfaceScalarField > surfaceTension(const phaseModel &phase) const
tmp< volScalarField > rho() const
Return the mixture density.
GeometricField< vector, fvPatchField, volMesh > volVectorField
Definition: volFieldsFwd.H:55
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)
autoPtr< dragCoeffFields > dragCoeffs() const
Return the drag coefficients for all of the interfaces.
tmp< volScalarField > Cvm(const phaseModel &phase) const
Return the virtual-mass coefficient for the given phase.
label nAlphaSubCycles(readLabel(alphaControls.lookup("nAlphaSubCycles")))
Calculate the divergence of the given field.
ITstream & lookup(const word &, bool recursive=false, bool patternMatch=true) const
Find and return an entry data stream.
Definition: dictionary.C:452
Area-weighted average a surfaceField creating a volField.
error FatalError
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:131
void correctBoundaryConditions()
Correct boundary field.
static const Vector zero
Definition: Vector.H:80
const Time & time() const
Return time.
Definition: IOobject.C:251
label timeIndex() const
Return the time index of the field.
static iteratorEnd end()
iteratorEnd set to beyond the end of any HashTable
Definition: HashTable.H:106
const dimensionSet dimless(0, 0, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:47
tmp< GeometricField< Type, fvPatchField, volMesh > > div(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcDiv.C:47
tmp< GeometricField< typename outerProduct< vector, Type >::type, fvPatchField, volMesh >> grad(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcGrad.C:52
dimensionedScalar pos(const dimensionedScalar &ds)
iterator find(const Key &)
Find and return an iterator set at the hashedEntry.
Definition: HashTable.C:139
GeometricField< vector, fvsPatchField, surfaceMesh > surfaceVectorField
IOdictionary(const IOobject &)
Construct given an IOobject.
Definition: IOdictionary.C:45
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
const objectRegistry & db() const
Return the local objectRegistry.
Definition: IOobject.C:245
const fvMesh & mesh() const
Constant access the mesh.
Definition: phaseSystemI.H:28
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic 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:117
A class for managing temporary objects.
Definition: PtrList.H:118
const Type & value() const
Return const reference to value.
word timeName
Definition: getTimeIndex.H:3
tmp< volScalarField > trho
Calculate the face-flux of the given field.