multiphaseMixtureThermo.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) 2013-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 
27 #include "alphaContactAngleFvPatchScalarField.H"
28 #include "Time.H"
29 #include "subCycle.H"
30 #include "MULES.H"
31 #include "fvcDiv.H"
32 #include "fvcGrad.H"
33 #include "fvcSnGrad.H"
34 #include "fvcFlux.H"
35 #include "fvcMeshPhi.H"
36 #include "surfaceInterpolate.H"
37 
38 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
39 
40 namespace Foam
41 {
42  defineTypeNameAndDebug(multiphaseMixtureThermo, 0);
43 }
44 
45 
46 const Foam::scalar Foam::multiphaseMixtureThermo::convertToRad =
48 
49 
50 // * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
51 
52 void Foam::multiphaseMixtureThermo::calcAlphas()
53 {
54  scalar level = 0.0;
55  alphas_ == 0.0;
56 
57  forAllIter(PtrDictionary<phaseModel>, phases_, phase)
58  {
59  alphas_ += level*phase();
60  level += 1.0;
61  }
62 }
63 
64 
65 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
66 
68 (
69  const volVectorField& U,
70  const surfaceScalarField& phi
71 )
72 :
73  psiThermo(U.mesh(), word::null),
74  phases_(lookup("phases"), phaseModel::iNew(p_, T_)),
75 
76  mesh_(U.mesh()),
77  U_(U),
78  phi_(phi),
79 
80  rhoPhi_
81  (
82  IOobject
83  (
84  "rhoPhi",
85  mesh_.time().timeName(),
86  mesh_,
87  IOobject::NO_READ,
88  IOobject::NO_WRITE
89  ),
90  mesh_,
91  dimensionedScalar("rhoPhi", dimMass/dimTime, 0.0)
92  ),
93 
94  alphas_
95  (
96  IOobject
97  (
98  "alphas",
99  mesh_.time().timeName(),
100  mesh_,
101  IOobject::NO_READ,
102  IOobject::AUTO_WRITE
103  ),
104  mesh_,
105  dimensionedScalar("alphas", dimless, 0.0)
106  ),
107 
108  sigmas_(lookup("sigmas")),
109  dimSigma_(1, 0, -2, 0, 0),
110  deltaN_
111  (
112  "deltaN",
113  1e-8/pow(average(mesh_.V()), 1.0/3.0)
114  )
115 {
116  calcAlphas();
117  alphas_.write();
118  correct();
119 }
120 
121 
122 // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
123 
125 {
126  forAllIter(PtrDictionary<phaseModel>, phases_, phasei)
127  {
128  phasei().correct();
129  }
130 
131  PtrDictionary<phaseModel>::iterator phasei = phases_.begin();
132 
133  psi_ = phasei()*phasei().thermo().psi();
134  mu_ = phasei()*phasei().thermo().mu();
135  alpha_ = phasei()*phasei().thermo().alpha();
136 
137  for (++phasei; phasei != phases_.end(); ++phasei)
138  {
139  psi_ += phasei()*phasei().thermo().psi();
140  mu_ += phasei()*phasei().thermo().mu();
141  alpha_ += phasei()*phasei().thermo().alpha();
142  }
143 }
144 
145 
147 {
148  forAllIter(PtrDictionary<phaseModel>, phases_, phasei)
149  {
150  phasei().thermo().rho() += phasei().thermo().psi()*dp;
151  }
152 }
153 
154 
156 {
157  bool ico = true;
158 
159  forAllConstIter(PtrDictionary<phaseModel>, phases_, phase)
160  {
161  ico &= phase().thermo().incompressible();
162  }
163 
164  return ico;
165 }
166 
167 
169 {
170  bool iso = true;
171 
172  forAllConstIter(PtrDictionary<phaseModel>, phases_, phase)
173  {
174  iso &= phase().thermo().incompressible();
175  }
176 
177  return iso;
178 }
179 
180 
182 (
183  const volScalarField& p,
184  const volScalarField& T
185 ) const
186 {
187  PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
188 
189  tmp<volScalarField> the(phasei()*phasei().thermo().he(p, T));
190 
191  for (++phasei; phasei != phases_.end(); ++phasei)
192  {
193  the.ref() += phasei()*phasei().thermo().he(p, T);
194  }
195 
196  return the;
197 }
198 
199 
201 (
202  const scalarField& p,
203  const scalarField& T,
204  const labelList& cells
205 ) const
206 {
207  PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
208 
209  tmp<scalarField> the
210  (
211  scalarField(phasei(), cells)*phasei().thermo().he(p, T, cells)
212  );
213 
214  for (++phasei; phasei != phases_.end(); ++phasei)
215  {
216  the.ref() +=
217  scalarField(phasei(), cells)*phasei().thermo().he(p, T, cells);
218  }
219 
220  return the;
221 }
222 
223 
225 (
226  const scalarField& p,
227  const scalarField& T,
228  const label patchi
229 ) const
230 {
231  PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
232 
233  tmp<scalarField> the
234  (
235  phasei().boundaryField()[patchi]*phasei().thermo().he(p, T, patchi)
236  );
237 
238  for (++phasei; phasei != phases_.end(); ++phasei)
239  {
240  the.ref() +=
241  phasei().boundaryField()[patchi]*phasei().thermo().he(p, T, patchi);
242  }
243 
244  return the;
245 }
246 
247 
249 {
250  PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
251 
252  tmp<volScalarField> thc(phasei()*phasei().thermo().hc());
253 
254  for (++phasei; phasei != phases_.end(); ++phasei)
255  {
256  thc.ref() += phasei()*phasei().thermo().hc();
257  }
258 
259  return thc;
260 }
261 
262 
264 (
265  const scalarField& h,
266  const scalarField& p,
267  const scalarField& T0,
268  const labelList& cells
269 ) const
270 {
272  return T0;
273 }
274 
275 
277 (
278  const scalarField& h,
279  const scalarField& p,
280  const scalarField& T0,
281  const label patchi
282 ) const
283 {
285  return T0;
286 }
287 
288 
290 {
291  PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
292 
293  tmp<volScalarField> trho(phasei()*phasei().thermo().rho());
294 
295  for (++phasei; phasei != phases_.end(); ++phasei)
296  {
297  trho.ref() += phasei()*phasei().thermo().rho();
298  }
299 
300  return trho;
301 }
302 
303 
305 (
306  const label patchi
307 ) const
308 {
309  PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
310 
311  tmp<scalarField> trho
312  (
313  phasei().boundaryField()[patchi]*phasei().thermo().rho(patchi)
314  );
315 
316  for (++phasei; phasei != phases_.end(); ++phasei)
317  {
318  trho.ref() +=
319  phasei().boundaryField()[patchi]*phasei().thermo().rho(patchi);
320  }
321 
322  return trho;
323 }
324 
325 
327 {
328  PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
329 
330  tmp<volScalarField> tCp(phasei()*phasei().thermo().Cp());
331 
332  for (++phasei; phasei != phases_.end(); ++phasei)
333  {
334  tCp.ref() += phasei()*phasei().thermo().Cp();
335  }
336 
337  return tCp;
338 }
339 
340 
342 (
343  const scalarField& p,
344  const scalarField& T,
345  const label patchi
346 ) const
347 {
348  PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
349 
350  tmp<scalarField> tCp
351  (
352  phasei().boundaryField()[patchi]*phasei().thermo().Cp(p, T, patchi)
353  );
354 
355  for (++phasei; phasei != phases_.end(); ++phasei)
356  {
357  tCp.ref() +=
358  phasei().boundaryField()[patchi]*phasei().thermo().Cp(p, T, patchi);
359  }
360 
361  return tCp;
362 }
363 
364 
366 {
367  PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
368 
369  tmp<volScalarField> tCv(phasei()*phasei().thermo().Cv());
370 
371  for (++phasei; phasei != phases_.end(); ++phasei)
372  {
373  tCv.ref() += phasei()*phasei().thermo().Cv();
374  }
375 
376  return tCv;
377 }
378 
379 
381 (
382  const scalarField& p,
383  const scalarField& T,
384  const label patchi
385 ) const
386 {
387  PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
388 
389  tmp<scalarField> tCv
390  (
391  phasei().boundaryField()[patchi]*phasei().thermo().Cv(p, T, patchi)
392  );
393 
394  for (++phasei; phasei != phases_.end(); ++phasei)
395  {
396  tCv.ref() +=
397  phasei().boundaryField()[patchi]*phasei().thermo().Cv(p, T, patchi);
398  }
399 
400  return tCv;
401 }
402 
403 
405 {
406  PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
407 
408  tmp<volScalarField> tgamma(phasei()*phasei().thermo().gamma());
409 
410  for (++phasei; phasei != phases_.end(); ++phasei)
411  {
412  tgamma.ref() += phasei()*phasei().thermo().gamma();
413  }
414 
415  return tgamma;
416 }
417 
418 
420 (
421  const scalarField& p,
422  const scalarField& T,
423  const label patchi
424 ) const
425 {
426  PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
427 
428  tmp<scalarField> tgamma
429  (
430  phasei().boundaryField()[patchi]*phasei().thermo().gamma(p, T, patchi)
431  );
432 
433  for (++phasei; phasei != phases_.end(); ++phasei)
434  {
435  tgamma.ref() +=
436  phasei().boundaryField()[patchi]
437  *phasei().thermo().gamma(p, T, patchi);
438  }
439 
440  return tgamma;
441 }
442 
443 
445 {
446  PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
447 
448  tmp<volScalarField> tCpv(phasei()*phasei().thermo().Cpv());
449 
450  for (++phasei; phasei != phases_.end(); ++phasei)
451  {
452  tCpv.ref() += phasei()*phasei().thermo().Cpv();
453  }
454 
455  return tCpv;
456 }
457 
458 
460 (
461  const scalarField& p,
462  const scalarField& T,
463  const label patchi
464 ) const
465 {
466  PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
467 
468  tmp<scalarField> tCpv
469  (
470  phasei().boundaryField()[patchi]*phasei().thermo().Cpv(p, T, patchi)
471  );
472 
473  for (++phasei; phasei != phases_.end(); ++phasei)
474  {
475  tCpv.ref() +=
476  phasei().boundaryField()[patchi]
477  *phasei().thermo().Cpv(p, T, patchi);
478  }
479 
480  return tCpv;
481 }
482 
483 
485 {
486  PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
487 
488  tmp<volScalarField> tCpByCpv(phasei()*phasei().thermo().CpByCpv());
489 
490  for (++phasei; phasei != phases_.end(); ++phasei)
491  {
492  tCpByCpv.ref() += phasei()*phasei().thermo().CpByCpv();
493  }
494 
495  return tCpByCpv;
496 }
497 
498 
500 (
501  const scalarField& p,
502  const scalarField& T,
503  const label patchi
504 ) const
505 {
506  PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
507 
508  tmp<scalarField> tCpByCpv
509  (
510  phasei().boundaryField()[patchi]*phasei().thermo().CpByCpv(p, T, patchi)
511  );
512 
513  for (++phasei; phasei != phases_.end(); ++phasei)
514  {
515  tCpByCpv.ref() +=
516  phasei().boundaryField()[patchi]
517  *phasei().thermo().CpByCpv(p, T, patchi);
518  }
519 
520  return tCpByCpv;
521 }
522 
523 
525 {
526  return mu()/rho();
527 }
528 
529 
531 (
532  const label patchi
533 ) const
534 {
535  return mu(patchi)/rho(patchi);
536 }
537 
538 
540 {
541  PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
542 
543  tmp<volScalarField> tkappa(phasei()*phasei().thermo().kappa());
544 
545  for (++phasei; phasei != phases_.end(); ++phasei)
546  {
547  tkappa.ref() += phasei()*phasei().thermo().kappa();
548  }
549 
550  return tkappa;
551 }
552 
553 
555 (
556  const label patchi
557 ) const
558 {
559  PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
560 
561  tmp<scalarField> tkappa
562  (
563  phasei().boundaryField()[patchi]*phasei().thermo().kappa(patchi)
564  );
565 
566  for (++phasei; phasei != phases_.end(); ++phasei)
567  {
568  tkappa.ref() +=
569  phasei().boundaryField()[patchi]*phasei().thermo().kappa(patchi);
570  }
571 
572  return tkappa;
573 }
574 
575 
577 (
578  const volScalarField& alphat
579 ) const
580 {
581  PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
582 
583  tmp<volScalarField> tkappaEff(phasei()*phasei().thermo().kappaEff(alphat));
584 
585  for (++phasei; phasei != phases_.end(); ++phasei)
586  {
587  tkappaEff.ref() += phasei()*phasei().thermo().kappaEff(alphat);
588  }
589 
590  return tkappaEff;
591 }
592 
593 
595 (
596  const scalarField& alphat,
597  const label patchi
598 ) const
599 {
600  PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
601 
602  tmp<scalarField> tkappaEff
603  (
604  phasei().boundaryField()[patchi]
605  *phasei().thermo().kappaEff(alphat, patchi)
606  );
607 
608  for (++phasei; phasei != phases_.end(); ++phasei)
609  {
610  tkappaEff.ref() +=
611  phasei().boundaryField()[patchi]
612  *phasei().thermo().kappaEff(alphat, patchi);
613  }
614 
615  return tkappaEff;
616 }
617 
618 
620 (
621  const volScalarField& alphat
622 ) const
623 {
624  PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
625 
626  tmp<volScalarField> talphaEff(phasei()*phasei().thermo().alphaEff(alphat));
627 
628  for (++phasei; phasei != phases_.end(); ++phasei)
629  {
630  talphaEff.ref() += phasei()*phasei().thermo().alphaEff(alphat);
631  }
632 
633  return talphaEff;
634 }
635 
636 
638 (
639  const scalarField& alphat,
640  const label patchi
641 ) const
642 {
643  PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
644 
645  tmp<scalarField> talphaEff
646  (
647  phasei().boundaryField()[patchi]
648  *phasei().thermo().alphaEff(alphat, patchi)
649  );
650 
651  for (++phasei; phasei != phases_.end(); ++phasei)
652  {
653  talphaEff.ref() +=
654  phasei().boundaryField()[patchi]
655  *phasei().thermo().alphaEff(alphat, patchi);
656  }
657 
658  return talphaEff;
659 }
660 
661 
663 {
664  PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
665 
666  tmp<volScalarField> trCv(phasei()/phasei().thermo().Cv());
667 
668  for (++phasei; phasei != phases_.end(); ++phasei)
669  {
670  trCv.ref() += phasei()/phasei().thermo().Cv();
671  }
672 
673  return trCv;
674 }
675 
676 
679 {
680  tmp<surfaceScalarField> tstf
681  (
683  (
684  IOobject
685  (
686  "surfaceTensionForce",
687  mesh_.time().timeName(),
688  mesh_
689  ),
690  mesh_,
692  (
693  "surfaceTensionForce",
694  dimensionSet(1, -2, -2, 0, 0),
695  0.0
696  )
697  )
698  );
699 
700  surfaceScalarField& stf = tstf.ref();
701 
702  forAllConstIter(PtrDictionary<phaseModel>, phases_, phase1)
703  {
704  const phaseModel& alpha1 = phase1();
705 
707  ++phase2;
708 
709  for (; phase2 != phases_.end(); ++phase2)
710  {
711  const phaseModel& alpha2 = phase2();
712 
713  sigmaTable::const_iterator sigma =
714  sigmas_.find(interfacePair(alpha1, alpha2));
715 
716  if (sigma == sigmas_.end())
717  {
719  << "Cannot find interface " << interfacePair(alpha1, alpha2)
720  << " in list of sigma values"
721  << exit(FatalError);
722  }
723 
724  stf += dimensionedScalar("sigma", dimSigma_, sigma())
725  *fvc::interpolate(K(alpha1, alpha2))*
726  (
727  fvc::interpolate(alpha2)*fvc::snGrad(alpha1)
728  - fvc::interpolate(alpha1)*fvc::snGrad(alpha2)
729  );
730  }
731  }
732 
733  return tstf;
734 }
735 
736 
738 {
739  const Time& runTime = mesh_.time();
740 
741  const dictionary& alphaControls = mesh_.solverDict("alpha");
742  label nAlphaSubCycles(readLabel(alphaControls.lookup("nAlphaSubCycles")));
743  scalar cAlpha(readScalar(alphaControls.lookup("cAlpha")));
744 
745  volScalarField& alpha = phases_.first();
746 
747  if (nAlphaSubCycles > 1)
748  {
749  surfaceScalarField rhoPhiSum(0.0*rhoPhi_);
750  dimensionedScalar totalDeltaT = runTime.deltaT();
751 
752  for
753  (
754  subCycle<volScalarField> alphaSubCycle(alpha, nAlphaSubCycles);
755  !(++alphaSubCycle).end();
756  )
757  {
758  solveAlphas(cAlpha);
759  rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi_;
760  }
761 
762  rhoPhi_ = rhoPhiSum;
763  }
764  else
765  {
766  solveAlphas(cAlpha);
767  }
768 }
769 
770 
771 Foam::tmp<Foam::surfaceVectorField> Foam::multiphaseMixtureThermo::nHatfv
772 (
773  const volScalarField& alpha1,
774  const volScalarField& alpha2
775 ) const
776 {
777  /*
778  // Cell gradient of alpha
779  volVectorField gradAlpha =
780  alpha2*fvc::grad(alpha1) - alpha1*fvc::grad(alpha2);
781 
782  // Interpolated face-gradient of alpha
783  surfaceVectorField gradAlphaf = fvc::interpolate(gradAlpha);
784  */
785 
786  surfaceVectorField gradAlphaf
787  (
789  - fvc::interpolate(alpha1)*fvc::interpolate(fvc::grad(alpha2))
790  );
791 
792  // Face unit interface normal
793  return gradAlphaf/(mag(gradAlphaf) + deltaN_);
794 }
795 
796 
797 Foam::tmp<Foam::surfaceScalarField> Foam::multiphaseMixtureThermo::nHatf
798 (
799  const volScalarField& alpha1,
800  const volScalarField& alpha2
801 ) const
802 {
803  // Face unit interface normal flux
804  return nHatfv(alpha1, alpha2) & mesh_.Sf();
805 }
806 
807 
808 // Correction for the boundary condition on the unit normal nHat on
809 // walls to produce the correct contact angle.
810 
811 // The dynamic contact angle is calculated from the component of the
812 // velocity on the direction of the interface, parallel to the wall.
813 
814 void Foam::multiphaseMixtureThermo::correctContactAngle
815 (
816  const phaseModel& alpha1,
817  const phaseModel& alpha2,
818  surfaceVectorField::Boundary& nHatb
819 ) const
820 {
821  const volScalarField::Boundary& gbf
822  = alpha1.boundaryField();
823 
824  const fvBoundaryMesh& boundary = mesh_.boundary();
825 
826  forAll(boundary, patchi)
827  {
828  if (isA<alphaContactAngleFvPatchScalarField>(gbf[patchi]))
829  {
830  const alphaContactAngleFvPatchScalarField& acap =
831  refCast<const alphaContactAngleFvPatchScalarField>(gbf[patchi]);
832 
833  vectorField& nHatPatch = nHatb[patchi];
834 
835  vectorField AfHatPatch
836  (
837  mesh_.Sf().boundaryField()[patchi]
838  /mesh_.magSf().boundaryField()[patchi]
839  );
840 
842  const_iterator tp =
843  acap.thetaProps().find(interfacePair(alpha1, alpha2));
844 
845  if (tp == acap.thetaProps().end())
846  {
848  << "Cannot find interface " << interfacePair(alpha1, alpha2)
849  << "\n in table of theta properties for patch "
850  << acap.patch().name()
851  << exit(FatalError);
852  }
853 
854  bool matched = (tp.key().first() == alpha1.name());
855 
856  scalar theta0 = convertToRad*tp().theta0(matched);
857  scalarField theta(boundary[patchi].size(), theta0);
858 
859  scalar uTheta = tp().uTheta();
860 
861  // Calculate the dynamic contact angle if required
862  if (uTheta > SMALL)
863  {
864  scalar thetaA = convertToRad*tp().thetaA(matched);
865  scalar thetaR = convertToRad*tp().thetaR(matched);
866 
867  // Calculated the component of the velocity parallel to the wall
868  vectorField Uwall
869  (
870  U_.boundaryField()[patchi].patchInternalField()
871  - U_.boundaryField()[patchi]
872  );
873  Uwall -= (AfHatPatch & Uwall)*AfHatPatch;
874 
875  // Find the direction of the interface parallel to the wall
876  vectorField nWall
877  (
878  nHatPatch - (AfHatPatch & nHatPatch)*AfHatPatch
879  );
880 
881  // Normalise nWall
882  nWall /= (mag(nWall) + SMALL);
883 
884  // Calculate Uwall resolved normal to the interface parallel to
885  // the interface
886  scalarField uwall(nWall & Uwall);
887 
888  theta += (thetaA - thetaR)*tanh(uwall/uTheta);
889  }
890 
891 
892  // Reset nHatPatch to correspond to the contact angle
893 
894  scalarField a12(nHatPatch & AfHatPatch);
895 
896  scalarField b1(cos(theta));
897 
898  scalarField b2(nHatPatch.size());
899 
900  forAll(b2, facei)
901  {
902  b2[facei] = cos(acos(a12[facei]) - theta[facei]);
903  }
904 
905  scalarField det(1.0 - a12*a12);
906 
907  scalarField a((b1 - a12*b2)/det);
908  scalarField b((b2 - a12*b1)/det);
909 
910  nHatPatch = a*AfHatPatch + b*nHatPatch;
911 
912  nHatPatch /= (mag(nHatPatch) + deltaN_.value());
913  }
914  }
915 }
916 
917 
918 Foam::tmp<Foam::volScalarField> Foam::multiphaseMixtureThermo::K
919 (
920  const phaseModel& alpha1,
921  const phaseModel& alpha2
922 ) const
923 {
924  tmp<surfaceVectorField> tnHatfv = nHatfv(alpha1, alpha2);
925 
926  correctContactAngle(alpha1, alpha2, tnHatfv.ref().boundaryFieldRef());
927 
928  // Simple expression for curvature
929  return -fvc::div(tnHatfv & mesh_.Sf());
930 }
931 
932 
935 {
936  tmp<volScalarField> tnearInt
937  (
938  new volScalarField
939  (
940  IOobject
941  (
942  "nearInterface",
943  mesh_.time().timeName(),
944  mesh_
945  ),
946  mesh_,
947  dimensionedScalar("nearInterface", dimless, 0.0)
948  )
949  );
950 
951  forAllConstIter(PtrDictionary<phaseModel>, phases_, phase)
952  {
953  tnearInt.ref() =
954  max(tnearInt(), pos(phase() - 0.01)*pos(0.99 - phase()));
955  }
956 
957  return tnearInt;
958 }
959 
960 
961 void Foam::multiphaseMixtureThermo::solveAlphas
962 (
963  const scalar cAlpha
964 )
965 {
966  static label nSolves=-1;
967  nSolves++;
968 
969  word alphaScheme("div(phi,alpha)");
970  word alpharScheme("div(phirb,alpha)");
971 
972  surfaceScalarField phic(mag(phi_/mesh_.magSf()));
973  phic = min(cAlpha*phic, max(phic));
974 
975  PtrList<surfaceScalarField> alphaPhiCorrs(phases_.size());
976  int phasei = 0;
977 
978  forAllIter(PtrDictionary<phaseModel>, phases_, phase)
979  {
980  phaseModel& alpha = phase();
981 
982  alphaPhiCorrs.set
983  (
984  phasei,
986  (
987  phi_.name() + alpha.name(),
988  fvc::flux
989  (
990  phi_,
991  alpha,
992  alphaScheme
993  )
994  )
995  );
996 
997  surfaceScalarField& alphaPhiCorr = alphaPhiCorrs[phasei];
998 
999  forAllIter(PtrDictionary<phaseModel>, phases_, phase2)
1000  {
1001  phaseModel& alpha2 = phase2();
1002 
1003  if (&alpha2 == &alpha) continue;
1004 
1005  surfaceScalarField phir(phic*nHatf(alpha, alpha2));
1006 
1007  alphaPhiCorr += fvc::flux
1008  (
1009  -fvc::flux(-phir, alpha2, alpharScheme),
1010  alpha,
1011  alpharScheme
1012  );
1013  }
1014 
1015  MULES::limit
1016  (
1017  1.0/mesh_.time().deltaT().value(),
1018  geometricOneField(),
1019  alpha,
1020  phi_,
1021  alphaPhiCorr,
1022  zeroField(),
1023  zeroField(),
1024  1,
1025  0,
1026  true
1027  );
1028 
1029  phasei++;
1030  }
1031 
1032  MULES::limitSum(alphaPhiCorrs);
1033 
1034  rhoPhi_ = dimensionedScalar("0", dimensionSet(1, 0, -1, 0, 0), 0);
1035 
1036  volScalarField sumAlpha
1037  (
1038  IOobject
1039  (
1040  "sumAlpha",
1041  mesh_.time().timeName(),
1042  mesh_
1043  ),
1044  mesh_,
1045  dimensionedScalar("sumAlpha", dimless, 0)
1046  );
1047 
1048 
1050 
1051 
1052  phasei = 0;
1053 
1054  forAllIter(PtrDictionary<phaseModel>, phases_, phase)
1055  {
1056  phaseModel& alpha = phase();
1057 
1058  surfaceScalarField& alphaPhi = alphaPhiCorrs[phasei];
1059  alphaPhi += upwind<scalar>(mesh_, phi_).flux(alpha);
1060 
1062  (
1063  IOobject
1064  (
1065  "Sp",
1066  mesh_.time().timeName(),
1067  mesh_
1068  ),
1069  mesh_,
1070  dimensionedScalar("Sp", alpha.dgdt().dimensions(), 0.0)
1071  );
1072 
1074  (
1075  IOobject
1076  (
1077  "Su",
1078  mesh_.time().timeName(),
1079  mesh_
1080  ),
1081  // Divergence term is handled explicitly to be
1082  // consistent with the explicit transport solution
1083  divU*min(alpha, scalar(1))
1084  );
1085 
1086  {
1087  const scalarField& dgdt = alpha.dgdt();
1088 
1089  forAll(dgdt, celli)
1090  {
1091  if (dgdt[celli] < 0.0 && alpha[celli] > 0.0)
1092  {
1093  Sp[celli] += dgdt[celli]*alpha[celli];
1094  Su[celli] -= dgdt[celli]*alpha[celli];
1095  }
1096  else if (dgdt[celli] > 0.0 && alpha[celli] < 1.0)
1097  {
1098  Sp[celli] -= dgdt[celli]*(1.0 - alpha[celli]);
1099  }
1100  }
1101  }
1102 
1103  forAllConstIter(PtrDictionary<phaseModel>, phases_, phase2)
1104  {
1105  const phaseModel& alpha2 = phase2();
1106 
1107  if (&alpha2 == &alpha) continue;
1108 
1109  const scalarField& dgdt2 = alpha2.dgdt();
1110 
1111  forAll(dgdt2, celli)
1112  {
1113  if (dgdt2[celli] > 0.0 && alpha2[celli] < 1.0)
1114  {
1115  Sp[celli] -= dgdt2[celli]*(1.0 - alpha2[celli]);
1116  Su[celli] += dgdt2[celli]*alpha[celli];
1117  }
1118  else if (dgdt2[celli] < 0.0 && alpha2[celli] > 0.0)
1119  {
1120  Sp[celli] += dgdt2[celli]*alpha2[celli];
1121  }
1122  }
1123  }
1124 
1126  (
1127  geometricOneField(),
1128  alpha,
1129  alphaPhi,
1130  Sp,
1131  Su
1132  );
1133 
1134  rhoPhi_ += fvc::interpolate(alpha.thermo().rho())*alphaPhi;
1135 
1136  Info<< alpha.name() << " volume fraction, min, max = "
1137  << alpha.weightedAverage(mesh_.V()).value()
1138  << ' ' << min(alpha).value()
1139  << ' ' << max(alpha).value()
1140  << endl;
1141 
1142  sumAlpha += alpha;
1143 
1144  phasei++;
1145  }
1146 
1147  Info<< "Phase-sum volume fraction, min, max = "
1148  << sumAlpha.weightedAverage(mesh_.V()).value()
1149  << ' ' << min(sumAlpha).value()
1150  << ' ' << max(sumAlpha).value()
1151  << endl;
1152 
1153  calcAlphas();
1154 }
1155 
1156 
1157 // ************************************************************************* //
tmp< surfaceScalarField > surfaceTensionForce() const
dimensionedScalar tanh(const dimensionedScalar &ds)
virtual tmp< volScalarField > rho() const
Density [kg/m^3] - uses current value of pressure.
Definition: psiThermo.C:93
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
const double e
Elementary charge.
Definition: doubleFloat.H:78
virtual tmp< volScalarField > Cpv() const
Heat capacity at constant pressure/volume [J/kg/K].
virtual tmp< volScalarField > Cp() const
Heat capacity at constant pressure [J/kg/K].
error FatalError
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:319
Info<< "Creating turbulence model\n"<< endl;tmp< volScalarField > talphaEff
Definition: setAlphaEff.H:2
tmp< GeometricField< Type, fvPatchField, volMesh > > div(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcDiv.C:47
virtual tmp< volScalarField > Cpv() const
Heat capacity at constant pressure/volume [J/kg/K].
virtual tmp< volScalarField > Cv() const
Heat capacity at constant volume [J/kg/K].
#define forAllIter(Container, container, iter)
Iterate across all elements in the container object of type.
Definition: UList.H:453
virtual tmp< volScalarField > kappa() const
Thermal diffusivity for temperature of mixture [J/m/s/K].
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
virtual tmp< volScalarField > Cp() const
Heat capacity at constant pressure [J/kg/K].
word alpharScheme("div(phirb,alpha)")
GeometricField< vector, fvsPatchField, surfaceMesh > surfaceVectorField
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:253
virtual tmp< volScalarField > rho() const
Density [kg/m^3].
tmp< volScalarField > trho
const dimensionedScalar sigma
Stefan-Boltzmann constant: default SI units: [W/m2/K4].
virtual tmp< volScalarField > CpByCpv() const
Heat capacity ratio [].
tmp< GeometricField< Type, fvPatchField, volMesh > > Sp(const volScalarField &sp, const GeometricField< Type, fvPatchField, volMesh > &vf)
Definition: fvcSup.C:67
const Type & value() const
Return const reference to value.
Calculate the snGrad of the given volField.
dimensionedScalar det(const dimensionedSphericalTensor &dt)
virtual tmp< volScalarField > kappaEff(const volScalarField &alphat) const
Effective thermal diffusivity of mixture [J/m/s/K].
label size() const
Return number of elements in list.
Definition: DLListBaseI.H:77
GeometricField< vector, fvPatchField, volMesh > volVectorField
Definition: volFieldsFwd.H:55
T & first()
Return the first element of the list.
Definition: UListI.H:114
DimensionedField< scalar, volMesh > Internal
Type of the internal field from which this GeometricField is derived.
CGAL::Exact_predicates_exact_constructions_kernel K
virtual tmp< volScalarField > kappaEff(const volScalarField &alphat) const
Effective thermal diffusivity of mixture [J/m/s/K].
virtual tmp< volScalarField > CpByCpv() const
Heat capacity ratio [].
surfaceScalarField phir(phic *interface.nHatf())
virtual tmp< volScalarField > mu() const
Dynamic viscosity of mixture [kg/m/s].
Definition: psiThermo.C:111
dictionary()
Construct top-level dictionary null.
Definition: dictionary.C:113
tmp< volScalarField > rCv() const
Return the phase-averaged reciprocal Cv.
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:52
void correctRho(const volScalarField &dp)
Update densities for given pressure change.
psiReactionThermo & thermo
Definition: createFields.H:31
dimensionedScalar pos(const dimensionedScalar &ds)
virtual tmp< volScalarField > hc() const
Chemical enthalpy [J/kg].
void solve()
Solve for the mixture phase-fractions.
surfaceScalarField alphaPhi(IOobject("alphaPhi", runTime.timeName(), mesh, IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE), phi *fvc::interpolate(alpha1))
stressControl lookup("compactNormalStress") >> compactNormalStress
dynamicFvMesh & mesh
virtual tmp< volScalarField > kappa() const
Thermal diffusivity for temperature of mixture [J/m/s/K].
dimensionedScalar cos(const dimensionedScalar &ds)
void limitSum(UPtrList< scalarField > &phiPsiCorrs)
Definition: MULES.C:51
virtual tmp< volScalarField > gamma() const
Gamma = Cp/Cv [].
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)
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Calculate the face-flux of the given field.
virtual tmp< volScalarField > Cv() const
Heat capacity at constant volume [J/kg/K].
word timeName
Definition: getTimeIndex.H:3
List< label > labelList
A List of labels.
Definition: labelList.H:56
virtual tmp< volScalarField > gamma() const
Gamma = Cp/Cv [].
dimensioned< Type > average(const DimensionedField< Type, GeoMesh > &df)
faceListList boundary(nPatches)
bool readScalar(const char *buf, doubleScalar &s)
Read whole of buf as a scalar. Return true if succesful.
Definition: doubleScalar.H:63
virtual tmp< volScalarField > hc() const
Chemical enthalpy [J/kg].
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
Definition: pEqn.H:29
label readLabel(Istream &is)
Definition: label.H:64
phaseModel & phase1
Calculate the divergence of the given field.
virtual bool incompressible() const
Return true if the equation of state is incompressible.
virtual tmp< scalarField > THE(const scalarField &h, const scalarField &p, const scalarField &T0, const labelList &cells) const
Temperature from enthalpy/internal energy for cell-set.
defineTypeNameAndDebug(combustionModel, 0)
virtual tmp< volScalarField > nu() const
Kinematic viscosity of mixture [m^2/s].
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
virtual volScalarField & he()
Enthalpy/Internal energy [J/kg].
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.
Calculate the mesh motion flux and convert fluxes from absolute to relative and back.
const dimensionedScalar h
Planck constant.
Definition: createFields.H:6
Internal & ref()
Return a reference to the dimensioned internal field.
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
multiphaseMixtureThermo(const volVectorField &U, const surfaceScalarField &phi)
Construct from components.
volScalarField alpha_
Laminar thermal diffusuvity [kg/m/s].
Definition: basicThermo.H:74
virtual tmp< volScalarField > alphaEff(const volScalarField &alphat) const
Effective thermal diffusivity of mixture [J/m/s/K].
label patchi
tmp< surfaceScalarField > absolute(const tmp< surfaceScalarField > &tphi, const volVectorField &U)
Return the given relative flux in absolute form.
Definition: fvcMeshPhi.C:188
const dimensionSet dimless(0, 0, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:47
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
virtual void correct()
Update properties.
virtual const volScalarField & alpha() const
Thermal diffusivity for enthalpy of mixture [kg/m/s].
Definition: basicThermo.C:490
tmp< volScalarField > nearInterface() const
Indicator of the proximity of the interface.
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
Definition: dimensionSets.H:51
messageStream Info
dimensioned< scalar > mag(const dimensioned< Type > &)
volScalarField divU(fvc::div(fvc::absolute(phi, U)))
Field< vector > vectorField
Specialisation of Field<T> for vector.
MULES: Multidimensional universal limiter for explicit solution.
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)
volScalarField mu_
Dynamic viscosity [kg/m/s].
Definition: psiThermo.H:62
const dimensionSet dimMass(1, 0, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:49
virtual tmp< volScalarField > alphaEff(const volScalarField &alphat) const
Effective thermal diffusivity of mixture [J/m/s/K].
volScalarField psi_
Compressibility [s^2/m^2].
Definition: psiThermo.H:59
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
virtual bool isochoric() const
Return true if the equation of state is isochoric.
tmp< GeometricField< Type, fvPatchField, volMesh > > Su(const GeometricField< Type, fvPatchField, volMesh > &su, const GeometricField< Type, fvPatchField, volMesh > &vf)
Definition: fvcSup.C:44
#define NotImplemented
Issue a FatalErrorIn for a function not currently implemented.
Definition: error.H:366
virtual volScalarField & he()
Enthalpy/Internal energy [J/kg].
tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > snGrad(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvcSnGrad.C:45
Namespace for OpenFOAM.
label nAlphaSubCycles(readLabel(alphaControls.lookup("nAlphaSubCycles")))
virtual void correct()
Update properties.