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 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  PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
309 
310  tmp<volScalarField> tCp(phasei()*phasei().thermo().Cp());
311 
312  for (++phasei; phasei != phases_.end(); ++phasei)
313  {
314  tCp() += phasei()*phasei().thermo().Cp();
315  }
316 
317  return tCp;
318 }
319 
320 
322 (
323  const scalarField& p,
324  const scalarField& T,
325  const label patchi
326 ) const
327 {
328  PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
329 
330  tmp<scalarField> tCp
331  (
332  phasei().boundaryField()[patchi]*phasei().thermo().Cp(p, T, patchi)
333  );
334 
335  for (++phasei; phasei != phases_.end(); ++phasei)
336  {
337  tCp() +=
338  phasei().boundaryField()[patchi]*phasei().thermo().Cp(p, T, patchi);
339  }
340 
341  return tCp;
342 }
343 
344 
346 {
347  PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
348 
349  tmp<volScalarField> tCv(phasei()*phasei().thermo().Cv());
350 
351  for (++phasei; phasei != phases_.end(); ++phasei)
352  {
353  tCv() += phasei()*phasei().thermo().Cv();
354  }
355 
356  return tCv;
357 }
358 
359 
361 (
362  const scalarField& p,
363  const scalarField& T,
364  const label patchi
365 ) const
366 {
367  PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
368 
369  tmp<scalarField> tCv
370  (
371  phasei().boundaryField()[patchi]*phasei().thermo().Cv(p, T, patchi)
372  );
373 
374  for (++phasei; phasei != phases_.end(); ++phasei)
375  {
376  tCv() +=
377  phasei().boundaryField()[patchi]*phasei().thermo().Cv(p, T, patchi);
378  }
379 
380  return tCv;
381 }
382 
383 
385 {
386  PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
387 
388  tmp<volScalarField> tgamma(phasei()*phasei().thermo().gamma());
389 
390  for (++phasei; phasei != phases_.end(); ++phasei)
391  {
392  tgamma() += phasei()*phasei().thermo().gamma();
393  }
394 
395  return tgamma;
396 }
397 
398 
400 (
401  const scalarField& p,
402  const scalarField& T,
403  const label patchi
404 ) const
405 {
406  PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
407 
408  tmp<scalarField> tgamma
409  (
410  phasei().boundaryField()[patchi]*phasei().thermo().gamma(p, T, patchi)
411  );
412 
413  for (++phasei; phasei != phases_.end(); ++phasei)
414  {
415  tgamma() +=
416  phasei().boundaryField()[patchi]
417  *phasei().thermo().gamma(p, T, patchi);
418  }
419 
420  return tgamma;
421 }
422 
423 
425 {
426  PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
427 
428  tmp<volScalarField> tCpv(phasei()*phasei().thermo().Cpv());
429 
430  for (++phasei; phasei != phases_.end(); ++phasei)
431  {
432  tCpv() += phasei()*phasei().thermo().Cpv();
433  }
434 
435  return tCpv;
436 }
437 
438 
440 (
441  const scalarField& p,
442  const scalarField& T,
443  const label patchi
444 ) const
445 {
446  PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
447 
448  tmp<scalarField> tCpv
449  (
450  phasei().boundaryField()[patchi]*phasei().thermo().Cpv(p, T, patchi)
451  );
452 
453  for (++phasei; phasei != phases_.end(); ++phasei)
454  {
455  tCpv() +=
456  phasei().boundaryField()[patchi]
457  *phasei().thermo().Cpv(p, T, patchi);
458  }
459 
460  return tCpv;
461 }
462 
463 
465 {
466  PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
467 
468  tmp<volScalarField> tCpByCpv(phasei()*phasei().thermo().CpByCpv());
469 
470  for (++phasei; phasei != phases_.end(); ++phasei)
471  {
472  tCpByCpv() += phasei()*phasei().thermo().CpByCpv();
473  }
474 
475  return tCpByCpv;
476 }
477 
478 
480 (
481  const scalarField& p,
482  const scalarField& T,
483  const label patchi
484 ) const
485 {
486  PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
487 
488  tmp<scalarField> tCpByCpv
489  (
490  phasei().boundaryField()[patchi]*phasei().thermo().CpByCpv(p, T, patchi)
491  );
492 
493  for (++phasei; phasei != phases_.end(); ++phasei)
494  {
495  tCpByCpv() +=
496  phasei().boundaryField()[patchi]
497  *phasei().thermo().CpByCpv(p, T, patchi);
498  }
499 
500  return tCpByCpv;
501 }
502 
503 
505 {
506  PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
507 
508  tmp<volScalarField> tkappa(phasei()*phasei().thermo().kappa());
509 
510  for (++phasei; phasei != phases_.end(); ++phasei)
511  {
512  tkappa() += phasei()*phasei().thermo().kappa();
513  }
514 
515  return tkappa;
516 }
517 
518 
520 (
521  const label patchi
522 ) const
523 {
524  PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
525 
526  tmp<scalarField> tkappa
527  (
528  phasei().boundaryField()[patchi]*phasei().thermo().kappa(patchi)
529  );
530 
531  for (++phasei; phasei != phases_.end(); ++phasei)
532  {
533  tkappa() +=
534  phasei().boundaryField()[patchi]*phasei().thermo().kappa(patchi);
535  }
536 
537  return tkappa;
538 }
539 
540 
542 (
543  const volScalarField& alphat
544 ) const
545 {
546  PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
547 
548  tmp<volScalarField> tkappaEff(phasei()*phasei().thermo().kappaEff(alphat));
549 
550  for (++phasei; phasei != phases_.end(); ++phasei)
551  {
552  tkappaEff() += phasei()*phasei().thermo().kappaEff(alphat);
553  }
554 
555  return tkappaEff;
556 }
557 
558 
560 (
561  const scalarField& alphat,
562  const label patchi
563 ) const
564 {
565  PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
566 
567  tmp<scalarField> tkappaEff
568  (
569  phasei().boundaryField()[patchi]
570  *phasei().thermo().kappaEff(alphat, patchi)
571  );
572 
573  for (++phasei; phasei != phases_.end(); ++phasei)
574  {
575  tkappaEff() +=
576  phasei().boundaryField()[patchi]
577  *phasei().thermo().kappaEff(alphat, patchi);
578  }
579 
580  return tkappaEff;
581 }
582 
583 
585 (
586  const volScalarField& alphat
587 ) const
588 {
589  PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
590 
591  tmp<volScalarField> talphaEff(phasei()*phasei().thermo().alphaEff(alphat));
592 
593  for (++phasei; phasei != phases_.end(); ++phasei)
594  {
595  talphaEff() += phasei()*phasei().thermo().alphaEff(alphat);
596  }
597 
598  return talphaEff;
599 }
600 
601 
603 (
604  const scalarField& alphat,
605  const label patchi
606 ) const
607 {
608  PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
609 
610  tmp<scalarField> talphaEff
611  (
612  phasei().boundaryField()[patchi]
613  *phasei().thermo().alphaEff(alphat, patchi)
614  );
615 
616  for (++phasei; phasei != phases_.end(); ++phasei)
617  {
618  talphaEff() +=
619  phasei().boundaryField()[patchi]
620  *phasei().thermo().alphaEff(alphat, patchi);
621  }
622 
623  return talphaEff;
624 }
625 
626 
628 {
629  PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
630 
631  tmp<volScalarField> trCv(phasei()/phasei().thermo().Cv());
632 
633  for (++phasei; phasei != phases_.end(); ++phasei)
634  {
635  trCv() += phasei()/phasei().thermo().Cv();
636  }
637 
638  return trCv;
639 }
640 
641 
644 {
645  tmp<surfaceScalarField> tstf
646  (
648  (
649  IOobject
650  (
651  "surfaceTensionForce",
652  mesh_.time().timeName(),
653  mesh_
654  ),
655  mesh_,
657  (
658  "surfaceTensionForce",
659  dimensionSet(1, -2, -2, 0, 0),
660  0.0
661  )
662  )
663  );
664 
665  surfaceScalarField& stf = tstf();
666 
667  forAllConstIter(PtrDictionary<phaseModel>, phases_, phase1)
668  {
669  const phaseModel& alpha1 = phase1();
670 
671  PtrDictionary<phaseModel>::const_iterator phase2 = phase1;
672  ++phase2;
673 
674  for (; phase2 != phases_.end(); ++phase2)
675  {
676  const phaseModel& alpha2 = phase2();
677 
678  sigmaTable::const_iterator sigma =
679  sigmas_.find(interfacePair(alpha1, alpha2));
680 
681  if (sigma == sigmas_.end())
682  {
684  (
685  "multiphaseMixtureThermo::surfaceTensionForce() const"
686  ) << "Cannot find interface " << interfacePair(alpha1, alpha2)
687  << " in list of sigma values"
688  << exit(FatalError);
689  }
690 
691  stf += dimensionedScalar("sigma", dimSigma_, sigma())
692  *fvc::interpolate(K(alpha1, alpha2))*
693  (
694  fvc::interpolate(alpha2)*fvc::snGrad(alpha1)
695  - fvc::interpolate(alpha1)*fvc::snGrad(alpha2)
696  );
697  }
698  }
699 
700  return tstf;
701 }
702 
703 
705 {
706  const Time& runTime = mesh_.time();
707 
708  const dictionary& pimpleDict = mesh_.solutionDict().subDict("PIMPLE");
709 
710  label nAlphaSubCycles(readLabel(pimpleDict.lookup("nAlphaSubCycles")));
711 
712  scalar cAlpha(readScalar(pimpleDict.lookup("cAlpha")));
713 
714 
715  volScalarField& alpha = phases_.first();
716 
717  if (nAlphaSubCycles > 1)
718  {
719  surfaceScalarField rhoPhiSum(0.0*rhoPhi_);
720  dimensionedScalar totalDeltaT = runTime.deltaT();
721 
722  for
723  (
724  subCycle<volScalarField> alphaSubCycle(alpha, nAlphaSubCycles);
725  !(++alphaSubCycle).end();
726  )
727  {
728  solveAlphas(cAlpha);
729  rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi_;
730  }
731 
732  rhoPhi_ = rhoPhiSum;
733  }
734  else
735  {
736  solveAlphas(cAlpha);
737  }
738 }
739 
740 
741 Foam::tmp<Foam::surfaceVectorField> Foam::multiphaseMixtureThermo::nHatfv
742 (
743  const volScalarField& alpha1,
744  const volScalarField& alpha2
745 ) const
746 {
747  /*
748  // Cell gradient of alpha
749  volVectorField gradAlpha =
750  alpha2*fvc::grad(alpha1) - alpha1*fvc::grad(alpha2);
751 
752  // Interpolated face-gradient of alpha
753  surfaceVectorField gradAlphaf = fvc::interpolate(gradAlpha);
754  */
755 
756  surfaceVectorField gradAlphaf
757  (
759  - fvc::interpolate(alpha1)*fvc::interpolate(fvc::grad(alpha2))
760  );
761 
762  // Face unit interface normal
763  return gradAlphaf/(mag(gradAlphaf) + deltaN_);
764 }
765 
766 
767 Foam::tmp<Foam::surfaceScalarField> Foam::multiphaseMixtureThermo::nHatf
768 (
769  const volScalarField& alpha1,
770  const volScalarField& alpha2
771 ) const
772 {
773  // Face unit interface normal flux
774  return nHatfv(alpha1, alpha2) & mesh_.Sf();
775 }
776 
777 
778 // Correction for the boundary condition on the unit normal nHat on
779 // walls to produce the correct contact angle.
780 
781 // The dynamic contact angle is calculated from the component of the
782 // velocity on the direction of the interface, parallel to the wall.
783 
784 void Foam::multiphaseMixtureThermo::correctContactAngle
785 (
786  const phaseModel& alpha1,
787  const phaseModel& alpha2,
788  surfaceVectorField::GeometricBoundaryField& nHatb
789 ) const
790 {
791  const volScalarField::GeometricBoundaryField& gbf
792  = alpha1.boundaryField();
793 
794  const fvBoundaryMesh& boundary = mesh_.boundary();
795 
796  forAll(boundary, patchi)
797  {
798  if (isA<alphaContactAngleFvPatchScalarField>(gbf[patchi]))
799  {
800  const alphaContactAngleFvPatchScalarField& acap =
801  refCast<const alphaContactAngleFvPatchScalarField>(gbf[patchi]);
802 
803  vectorField& nHatPatch = nHatb[patchi];
804 
805  vectorField AfHatPatch
806  (
807  mesh_.Sf().boundaryField()[patchi]
808  /mesh_.magSf().boundaryField()[patchi]
809  );
810 
811  alphaContactAngleFvPatchScalarField::thetaPropsTable::
812  const_iterator tp =
813  acap.thetaProps().find(interfacePair(alpha1, alpha2));
814 
815  if (tp == acap.thetaProps().end())
816  {
818  (
819  "multiphaseMixtureThermo::correctContactAngle"
820  "(const phaseModel& alpha1, const phaseModel& alpha2, "
821  "fvPatchVectorFieldField& nHatb) const"
822  ) << "Cannot find interface " << interfacePair(alpha1, alpha2)
823  << "\n in table of theta properties for patch "
824  << acap.patch().name()
825  << exit(FatalError);
826  }
827 
828  bool matched = (tp.key().first() == alpha1.name());
829 
830  scalar theta0 = convertToRad*tp().theta0(matched);
831  scalarField theta(boundary[patchi].size(), theta0);
832 
833  scalar uTheta = tp().uTheta();
834 
835  // Calculate the dynamic contact angle if required
836  if (uTheta > SMALL)
837  {
838  scalar thetaA = convertToRad*tp().thetaA(matched);
839  scalar thetaR = convertToRad*tp().thetaR(matched);
840 
841  // Calculated the component of the velocity parallel to the wall
842  vectorField Uwall
843  (
844  U_.boundaryField()[patchi].patchInternalField()
845  - U_.boundaryField()[patchi]
846  );
847  Uwall -= (AfHatPatch & Uwall)*AfHatPatch;
848 
849  // Find the direction of the interface parallel to the wall
850  vectorField nWall
851  (
852  nHatPatch - (AfHatPatch & nHatPatch)*AfHatPatch
853  );
854 
855  // Normalise nWall
856  nWall /= (mag(nWall) + SMALL);
857 
858  // Calculate Uwall resolved normal to the interface parallel to
859  // the interface
860  scalarField uwall(nWall & Uwall);
861 
862  theta += (thetaA - thetaR)*tanh(uwall/uTheta);
863  }
864 
865 
866  // Reset nHatPatch to correspond to the contact angle
867 
868  scalarField a12(nHatPatch & AfHatPatch);
869 
870  scalarField b1(cos(theta));
871 
872  scalarField b2(nHatPatch.size());
873 
874  forAll(b2, facei)
875  {
876  b2[facei] = cos(acos(a12[facei]) - theta[facei]);
877  }
878 
879  scalarField det(1.0 - a12*a12);
880 
881  scalarField a((b1 - a12*b2)/det);
882  scalarField b((b2 - a12*b1)/det);
883 
884  nHatPatch = a*AfHatPatch + b*nHatPatch;
885 
886  nHatPatch /= (mag(nHatPatch) + deltaN_.value());
887  }
888  }
889 }
890 
891 
892 Foam::tmp<Foam::volScalarField> Foam::multiphaseMixtureThermo::K
893 (
894  const phaseModel& alpha1,
895  const phaseModel& alpha2
896 ) const
897 {
898  tmp<surfaceVectorField> tnHatfv = nHatfv(alpha1, alpha2);
899 
900  correctContactAngle(alpha1, alpha2, tnHatfv().boundaryField());
901 
902  // Simple expression for curvature
903  return -fvc::div(tnHatfv & mesh_.Sf());
904 }
905 
906 
909 {
910  tmp<volScalarField> tnearInt
911  (
912  new volScalarField
913  (
914  IOobject
915  (
916  "nearInterface",
917  mesh_.time().timeName(),
918  mesh_
919  ),
920  mesh_,
921  dimensionedScalar("nearInterface", dimless, 0.0)
922  )
923  );
924 
925  forAllConstIter(PtrDictionary<phaseModel>, phases_, phase)
926  {
927  tnearInt() = max(tnearInt(), pos(phase() - 0.01)*pos(0.99 - phase()));
928  }
929 
930  return tnearInt;
931 }
932 
933 
934 void Foam::multiphaseMixtureThermo::solveAlphas
935 (
936  const scalar cAlpha
937 )
938 {
939  static label nSolves=-1;
940  nSolves++;
941 
942  word alphaScheme("div(phi,alpha)");
943  word alpharScheme("div(phirb,alpha)");
944 
945  surfaceScalarField phic(mag(phi_/mesh_.magSf()));
946  phic = min(cAlpha*phic, max(phic));
947 
948  PtrList<surfaceScalarField> phiAlphaCorrs(phases_.size());
949  int phasei = 0;
950 
951  forAllIter(PtrDictionary<phaseModel>, phases_, phase)
952  {
953  phaseModel& alpha = phase();
954 
955  phiAlphaCorrs.set
956  (
957  phasei,
959  (
960  fvc::flux
961  (
962  phi_,
963  alpha,
964  alphaScheme
965  )
966  )
967  );
968 
969  surfaceScalarField& phiAlphaCorr = phiAlphaCorrs[phasei];
970 
971  forAllIter(PtrDictionary<phaseModel>, phases_, phase2)
972  {
973  phaseModel& alpha2 = phase2();
974 
975  if (&alpha2 == &alpha) continue;
976 
977  surfaceScalarField phir(phic*nHatf(alpha, alpha2));
978 
979  phiAlphaCorr += fvc::flux
980  (
981  -fvc::flux(-phir, alpha2, alpharScheme),
982  alpha,
984  );
985  }
986 
988  (
989  1.0/mesh_.time().deltaT().value(),
990  geometricOneField(),
991  alpha,
992  phi_,
993  phiAlphaCorr,
994  zeroField(),
995  zeroField(),
996  1,
997  0,
998  3,
999  true
1000  );
1001 
1002  phasei++;
1003  }
1004 
1005  MULES::limitSum(phiAlphaCorrs);
1006 
1007  rhoPhi_ = dimensionedScalar("0", dimensionSet(1, 0, -1, 0, 0), 0);
1008 
1009  volScalarField sumAlpha
1010  (
1011  IOobject
1012  (
1013  "sumAlpha",
1014  mesh_.time().timeName(),
1015  mesh_
1016  ),
1017  mesh_,
1018  dimensionedScalar("sumAlpha", dimless, 0)
1019  );
1020 
1021 
1023 
1024 
1025  phasei = 0;
1026 
1027  forAllIter(PtrDictionary<phaseModel>, phases_, phase)
1028  {
1029  phaseModel& alpha = phase();
1030 
1031  surfaceScalarField& phiAlpha = phiAlphaCorrs[phasei];
1032  phiAlpha += upwind<scalar>(mesh_, phi_).flux(alpha);
1033 
1035  (
1036  IOobject
1037  (
1038  "Sp",
1039  mesh_.time().timeName(),
1040  mesh_
1041  ),
1042  mesh_,
1043  dimensionedScalar("Sp", alpha.dgdt().dimensions(), 0.0)
1044  );
1045 
1047  (
1048  IOobject
1049  (
1050  "Su",
1051  mesh_.time().timeName(),
1052  mesh_
1053  ),
1054  // Divergence term is handled explicitly to be
1055  // consistent with the explicit transport solution
1056  divU*min(alpha, scalar(1))
1057  );
1058 
1059  {
1060  const scalarField& dgdt = alpha.dgdt();
1061 
1062  forAll(dgdt, celli)
1063  {
1064  if (dgdt[celli] < 0.0 && alpha[celli] > 0.0)
1065  {
1066  Sp[celli] += dgdt[celli]*alpha[celli];
1067  Su[celli] -= dgdt[celli]*alpha[celli];
1068  }
1069  else if (dgdt[celli] > 0.0 && alpha[celli] < 1.0)
1070  {
1071  Sp[celli] -= dgdt[celli]*(1.0 - alpha[celli]);
1072  }
1073  }
1074  }
1075 
1076  forAllConstIter(PtrDictionary<phaseModel>, phases_, phase2)
1077  {
1078  const phaseModel& alpha2 = phase2();
1079 
1080  if (&alpha2 == &alpha) continue;
1081 
1082  const scalarField& dgdt2 = alpha2.dgdt();
1083 
1084  forAll(dgdt2, celli)
1085  {
1086  if (dgdt2[celli] > 0.0 && alpha2[celli] < 1.0)
1087  {
1088  Sp[celli] -= dgdt2[celli]*(1.0 - alpha2[celli]);
1089  Su[celli] += dgdt2[celli]*alpha[celli];
1090  }
1091  else if (dgdt2[celli] < 0.0 && alpha2[celli] > 0.0)
1092  {
1093  Sp[celli] += dgdt2[celli]*alpha2[celli];
1094  }
1095  }
1096  }
1097 
1099  (
1100  geometricOneField(),
1101  alpha,
1102  phiAlpha,
1103  Sp,
1104  Su
1105  );
1106 
1107  rhoPhi_ += fvc::interpolate(alpha.thermo().rho())*phiAlpha;
1108 
1109  Info<< alpha.name() << " volume fraction, min, max = "
1110  << alpha.weightedAverage(mesh_.V()).value()
1111  << ' ' << min(alpha).value()
1112  << ' ' << max(alpha).value()
1113  << endl;
1114 
1115  sumAlpha += alpha;
1116 
1117  phasei++;
1118  }
1119 
1120  Info<< "Phase-sum volume fraction, min, max = "
1121  << sumAlpha.weightedAverage(mesh_.V()).value()
1122  << ' ' << min(sumAlpha).value()
1123  << ' ' << max(sumAlpha).value()
1124  << endl;
1125 
1126  calcAlphas();
1127 }
1128 
1129 
1130 // ************************************************************************* //
messageStream Info
surfaceScalarField phir(phic *interface.nHatf())
virtual bool incompressible() const
Return true if the equation of state is incompressible.
#define forAllIter(Container, container, iter)
Definition: UList.H:440
CGAL::Exact_predicates_exact_constructions_kernel K
surfaceScalarField phiAlpha(fvc::flux(phi, alpha1, alphaScheme))
dimensionedScalar det(const dimensionedSphericalTensor &dt)
dimensionedScalar pos(const dimensionedScalar &ds)
dimensioned< Type > min(const dimensioned< Type > &, const dimensioned< Type > &)
#define forAllConstIter(Container, container, iter)
Definition: UList.H:461
dynamicFvMesh & mesh
A class for managing temporary objects.
Definition: PtrList.H:90
virtual tmp< volScalarField > hc() const
Chemical enthalpy [J/kg].
label patchi
const dimensionedScalar sigma
Stefan-Boltzmann constant: default SI units: [W/m2/K4].
GeometricField< vector, fvPatchField, volMesh > volVectorField
Definition: volFieldsFwd.H:55
Calculate the face-flux of the given field.
const double e
Definition: doubleFloat.H:78
const dimensionSet dimless(0, 0, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:47
virtual bool isochoric() const
Return true if the equation of state is isochoric.
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
#define notImplemented(fn)
Definition: error.H:351
void explicitSolve(const RdeltaTType &rDeltaT, const RhoType &rho, volScalarField &psi, const surfaceScalarField &phiPsi, const SpType &Sp, const SuType &Su)
virtual tmp< volScalarField > kappa() const
Thermal diffusivity for temperature of mixture [J/m/s/K].
List< label > labelList
A List of labels.
Definition: labelList.H:56
dimensioned< scalar > mag(const dimensioned< Type > &)
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
Definition: dimensionSets.H:51
void limitSum(UPtrList< scalarField > &phiPsiCorrs)
Definition: MULES.C:72
virtual tmp< volScalarField > Cpv() const
Heat capacity at constant pressure/volume [J/kg/K].
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.
virtual tmp< volScalarField > Cp() const
Heat capacity at constant pressure [J/kg/K].
Calculate the snGrad of the given volField.
Calculate the mesh motion flux and convert fluxes from absolute to relative and back.
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
psiReactionThermo & thermo
Definition: createFields.H:32
tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > interpolate(const GeometricField< Type, fvPatchField, volMesh > &vf, const surfaceScalarField &faceFlux, Istream &schemeData)
Interpolate field onto faces using scheme given by Istream.
virtual tmp< volScalarField > alphaEff(const volScalarField &alphat) const
Effective thermal diffusivity of mixture [J/m/s/K].
tmp< volScalarField > nearInterface() const
Indicator of the proximity of the interface.
#define forAll(list, i)
Definition: UList.H:421
tmp< surfaceScalarField > absolute(const tmp< surfaceScalarField > &tphi, const volVectorField &U)
Return the given relative flux in absolute form.
Definition: fvcMeshPhi.C:187
virtual tmp< volScalarField > gamma() const
gamma = Cp/Cv []
tmp< GeometricField< Type, fvPatchField, volMesh > > div(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcDiv.C:47
const dimensionedScalar alpha
Fine-structure constant: default SI units: [].
multiphaseMixtureThermo(const volVectorField &U, const surfaceScalarField &phi)
Construct from components.
GeometricField< vector, fvsPatchField, surfaceMesh > surfaceVectorField
MULES: Multidimensional universal limiter for explicit solution.
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:125
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:249
label nAlphaSubCycles(readLabel(alphaControls.lookup("nAlphaSubCycles")))
const dimensionedScalar b
Wien displacement law constant: default SI units: [m.K].
Definition: createFields.H:28
virtual volScalarField & he()
Enthalpy/Internal energy [J/kg].
Calculate the gradient of the given field.
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
virtual tmp< volScalarField > CpByCpv() const
Heat capacity ratio [].
virtual void correct()
Update properties.
iterator begin()
Return an iterator to begin traversing the UList.
Definition: UListI.H:216
virtual tmp< volScalarField > rho() const
Density [kg/m^3].
error FatalError
Definition: error.H:303
word timeName
Definition: getTimeIndex.H:3
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
virtual tmp< volScalarField > kappaEff(const volScalarField &alphat) const
Effective thermal diffusivity of mixture [J/m/s/K].
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
tmp< volScalarField > rCv() const
Return the phase-averaged reciprocal Cv.
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
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 label nLimiterIter, const bool returnCorr)
tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > flux(const surfaceScalarField &phi, const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvcFlux.C:45
tmp< GeometricField< Type, fvPatchField, volMesh > > Sp(const volScalarField &sp, const GeometricField< Type, fvPatchField, volMesh > &vf)
Definition: fvcSup.C:67
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:52
void correctRho(const volScalarField &dp)
Update densities for given pressure change.
word alpharScheme("div(phirb,alpha)")
const volScalarField & alphaEff
Definition: setAlphaEff.H:93
const dimensionedScalar h
Planck constant.
Definition: createFields.H:6
dimensionedScalar tanh(const dimensionedScalar &ds)
volScalarField & he
Definition: YEEqn.H:56
dimensionedScalar cos(const dimensionedScalar &ds)
label phasei
Definition: pEqn.H:37
tmp< GeometricField< typename outerProduct< vector, Type >::type, fvPatchField, volMesh >> grad(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcGrad.C:52
Calculate the divergence of the given field.
#define FatalErrorIn(fn)
Definition: error.H:318
Field< vector > vectorField
Specialisation of Field<T> for vector.
defineTypeNameAndDebug(combustionModel, 0)
const dimensionedScalar kappa
Coulomb constant: default SI units: [N.m2/C2].
volScalarField divU(fvc::div(fvc::absolute(phi, U)))
dimensionedScalar acos(const dimensionedScalar &ds)
const dimensionSet dimMass(1, 0, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:49
phiP1 boundaryField()
tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > snGrad(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvcSnGrad.C:45
virtual void correct()
Update properties.
Surface Interpolation.
tmp< volScalarField > trho
DimensionedField< Type, GeoMesh > DimensionedInternalField
dimensioned< Type > average(const DimensionedField< Type, GeoMesh > &df)
bool readScalar(const char *buf, doubleScalar &s)
Read whole of buf as a scalar. Return true if succesful.
Definition: doubleScalar.H:62
void solve()
Solve for the mixture phase-fractions.
Info<< "Creating turbulence model\n"<< endl;tmp< volScalarField > talphaEff
Definition: setAlphaEff.H:2
tmp< surfaceScalarField > surfaceTensionForce() const
faceListList boundary(nPatches)