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