All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
phaseSystem.C
Go to the documentation of this file.
1 /*---------------------------------------------------------------------------*\
2  ========= |
3  \\ / F ield | OpenFOAM: The Open Source CFD Toolbox
4  \\ / O peration | Website: https://openfoam.org
5  \\ / A nd | Copyright (C) 2015-2021 OpenFOAM Foundation
6  \\/ M anipulation |
7 -------------------------------------------------------------------------------
8 License
9  This file is part of OpenFOAM.
10 
11  OpenFOAM is free software: you can redistribute it and/or modify it
12  under the terms of the GNU General Public License as published by
13  the Free Software Foundation, either version 3 of the License, or
14  (at your option) any later version.
15 
16  OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
17  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
19  for more details.
20 
21  You should have received a copy of the GNU General Public License
22  along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
23 
24 \*---------------------------------------------------------------------------*/
25 
26 #include "phaseSystem.H"
27 #include "surfaceTensionModel.H"
28 #include "aspectRatioModel.H"
29 #include "surfaceInterpolate.H"
30 #include "fvcDdt.H"
31 #include "localEulerDdtScheme.H"
32 #include "fvcDiv.H"
33 #include "fvcGrad.H"
34 #include "fvcSnGrad.H"
35 #include "CorrectPhi.H"
36 #include "fvcMeshPhi.H"
37 #include "alphaContactAngleFvPatchScalarField.H"
38 #include "unitConversion.H"
39 #include "dragModel.H"
42 #include "pimpleControl.H"
43 #include "pressureReference.H"
44 
45 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
46 
47 namespace Foam
48 {
49  defineTypeNameAndDebug(phaseSystem, 0);
50  defineRunTimeSelectionTable(phaseSystem, dictionary);
51 }
52 
53 const Foam::word Foam::phaseSystem::propertiesName("phaseProperties");
54 
55 
56 // * * * * * * * * * * * * Protected Member Functions * * * * * * * * * * * //
57 
59 (
60  const phaseModelList& phaseModels
61 ) const
62 {
63  tmp<surfaceScalarField> tmpPhi
64  (
66  (
67  "phi",
68  fvc::interpolate(phaseModels[0])*phaseModels[0].phi()
69  )
70  );
71 
72  for (label phasei=1; phasei<phaseModels.size(); phasei++)
73  {
74  tmpPhi.ref() +=
75  fvc::interpolate(phaseModels[phasei])*phaseModels[phasei].phi();
76  }
77 
78  return tmpPhi;
79 }
80 
81 
83 (
84  const dictTable& modelDicts
85 )
86 {
87  forAllConstIter(dictTable, modelDicts, iter)
88  {
89  const phasePairKey& key = iter.key();
90 
91  // pair already exists
92  if (phasePairs_.found(key))
93  {}
94 
95  // new ordered pair
96  else if (key.ordered())
97  {
98  phasePairs_.insert
99  (
100  key,
101  autoPtr<phasePair>
102  (
103  new orderedPhasePair
104  (
105  phaseModels_[key.first()],
106  phaseModels_[key.second()]
107  )
108  )
109  );
110  }
111 
112  // new unordered pair
113  else
114  {
115  phasePairs_.insert
116  (
117  key,
118  autoPtr<phasePair>
119  (
120  new phasePair
121  (
122  phaseModels_[key.first()],
123  phaseModels_[key.second()]
124  )
125  )
126  );
127  }
128  }
129 }
130 
131 
133 {
134  tmp<volScalarField> sumAlphaMoving
135  (
137  (
138  "sumAlphaMoving",
139  movingPhaseModels_[0],
140  calculatedFvPatchScalarField::typeName
141  )
142  );
143 
144  for
145  (
146  label movingPhasei=1;
147  movingPhasei<movingPhaseModels_.size();
148  movingPhasei++
149  )
150  {
151  sumAlphaMoving.ref() += movingPhaseModels_[movingPhasei];
152  }
153 
154  return sumAlphaMoving;
155 }
156 
157 
159 {
160  // Calculate the mean velocity difference with respect to Um0
161  // from the current velocity of the moving phases
162  volVectorField dUm(Um0);
163 
164  forAll(movingPhaseModels_, movingPhasei)
165  {
166  dUm -=
167  movingPhaseModels_[movingPhasei]
168  *movingPhaseModels_[movingPhasei].U();
169  }
170 
171  forAll(movingPhaseModels_, movingPhasei)
172  {
173  movingPhaseModels_[movingPhasei].URef() += dUm;
174  }
175 }
176 
177 
179 (
180  const PtrList<surfaceScalarField>& alphafs,
181  const surfaceScalarField& phim0
182 )
183 {
184  // Calculate the mean flux difference with respect to phim0
185  // from the current flux of the moving phases
186  surfaceScalarField dphim(phim0);
187 
188  forAll(movingPhaseModels_, movingPhasei)
189  {
190  dphim -=
191  alphafs[movingPhaseModels_[movingPhasei].index()]
192  *movingPhaseModels_[movingPhasei].phi();
193  }
194 
195  forAll(movingPhaseModels_, movingPhasei)
196  {
197  movingPhaseModels_[movingPhasei].phiRef() += dphim;
198  }
199 }
200 
201 
203 (
204  const volScalarField& alpha1,
205  const volScalarField& alpha2
206 ) const
207 {
208  /*
209  // Cell gradient of alpha
210  volVectorField gradAlpha =
211  alpha2*fvc::grad(alpha1) - alpha1*fvc::grad(alpha2);
212 
213  // Interpolated face-gradient of alpha
214  surfaceVectorField gradAlphaf = fvc::interpolate(gradAlpha);
215  */
216 
217  surfaceVectorField gradAlphaf
218  (
220  - fvc::interpolate(alpha1)*fvc::interpolate(fvc::grad(alpha2))
221  );
222 
223  // Face unit interface normal
224  return gradAlphaf/(mag(gradAlphaf) + deltaN_);
225 }
226 
227 
229 (
230  const volScalarField& alpha1,
231  const volScalarField& alpha2
232 ) const
233 {
234  // Face unit interface normal flux
235  return nHatfv(alpha1, alpha2) & mesh_.Sf();
236 }
237 
238 
240 (
241  const phaseModel& phase1,
242  const phaseModel& phase2,
243  surfaceVectorField::Boundary& nHatb
244 ) const
245 {
246  const volScalarField::Boundary& gbf
247  = phase1.boundaryField();
248 
249  const fvBoundaryMesh& boundary = mesh_.boundary();
250 
251  forAll(boundary, patchi)
252  {
253  if (isA<alphaContactAngleFvPatchScalarField>(gbf[patchi]))
254  {
255  const alphaContactAngleFvPatchScalarField& acap =
256  refCast<const alphaContactAngleFvPatchScalarField>(gbf[patchi]);
257 
258  vectorField& nHatPatch = nHatb[patchi];
259 
260  vectorField AfHatPatch
261  (
262  mesh_.Sf().boundaryField()[patchi]
263  /mesh_.magSf().boundaryField()[patchi]
264  );
265 
266  alphaContactAngleFvPatchScalarField::thetaPropsTable::
267  const_iterator tp =
268  acap.thetaProps()
269  .find(phasePairKey(phase1.name(), phase2.name()));
270 
271  if (tp == acap.thetaProps().end())
272  {
274  << "Cannot find interface "
275  << phasePairKey(phase1.name(), phase2.name())
276  << "\n in table of theta properties for patch "
277  << acap.patch().name()
278  << exit(FatalError);
279  }
280 
281  bool matched = (tp.key().first() == phase1.name());
282 
283  scalar theta0 = degToRad(tp().theta0(matched));
284  scalarField theta(boundary[patchi].size(), theta0);
285 
286  scalar uTheta = tp().uTheta();
287 
288  // Calculate the dynamic contact angle if required
289  if (uTheta > small)
290  {
291  scalar thetaA = degToRad(tp().thetaA(matched));
292  scalar thetaR = degToRad(tp().thetaR(matched));
293 
294  // Calculated the component of the velocity parallel to the wall
295  vectorField Uwall
296  (
297  phase1.U()().boundaryField()[patchi].patchInternalField()
298  - phase1.U()().boundaryField()[patchi]
299  );
300  Uwall -= (AfHatPatch & Uwall)*AfHatPatch;
301 
302  // Find the direction of the interface parallel to the wall
303  vectorField nWall
304  (
305  nHatPatch - (AfHatPatch & nHatPatch)*AfHatPatch
306  );
307 
308  // Normalise nWall
309  nWall /= (mag(nWall) + small);
310 
311  // Calculate Uwall resolved normal to the interface parallel to
312  // the interface
313  scalarField uwall(nWall & Uwall);
314 
315  theta += (thetaA - thetaR)*tanh(uwall/uTheta);
316  }
317 
318 
319  // Reset nHatPatch to correspond to the contact angle
320 
321  scalarField a12(nHatPatch & AfHatPatch);
322 
323  scalarField b1(cos(theta));
324 
325  scalarField b2(nHatPatch.size());
326 
327  forAll(b2, facei)
328  {
329  b2[facei] = cos(acos(a12[facei]) - theta[facei]);
330  }
331 
332  scalarField det(1 - a12*a12);
333 
334  scalarField a((b1 - a12*b2)/det);
335  scalarField b((b2 - a12*b1)/det);
336 
337  nHatPatch = a*AfHatPatch + b*nHatPatch;
338 
339  nHatPatch /= (mag(nHatPatch) + deltaN_.value());
340  }
341  }
342 }
343 
344 
346 (
347  const phaseModel& phase1,
348  const phaseModel& phase2
349 ) const
350 {
351  tmp<surfaceVectorField> tnHatfv = nHatfv(phase1, phase2);
352 
353  correctContactAngle(phase1, phase2, tnHatfv.ref().boundaryFieldRef());
354 
355  // Simple expression for curvature
356  return -fvc::div(tnHatfv & mesh_.Sf());
357 }
358 
359 
360 // * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
361 
363 (
364  const fvMesh& mesh
365 )
366 :
367  IOdictionary
368  (
369  IOobject
370  (
371  "phaseProperties",
372  mesh.time().constant(),
373  mesh,
374  IOobject::MUST_READ_IF_MODIFIED,
375  IOobject::NO_WRITE
376  )
377  ),
378 
379  mesh_(mesh),
380 
381  referencePhaseName_(lookupOrDefault("referencePhase", word::null)),
382 
383  phaseModels_
384  (
385  lookup("phases"),
386  phaseModel::iNew(*this, referencePhaseName_)
387  ),
388 
389  phi_(calcPhi(phaseModels_)),
390 
391  dpdt_
392  (
393  IOobject
394  (
395  "dpdt",
396  mesh.time().timeName(),
397  mesh
398  ),
399  mesh,
401  ),
402 
403  MRF_(mesh_),
404 
405  cAlphas_(lookupOrDefault("interfaceCompression", cAlphaTable())),
406 
407  deltaN_
408  (
409  "deltaN",
410  1e-8/pow(average(mesh_.V()), 1.0/3.0)
411  )
412 {
413  // Groupings
414  label movingPhasei = 0;
415  label stationaryPhasei = 0;
416  label anisothermalPhasei = 0;
417  label multiComponentPhasei = 0;
418  forAll(phaseModels_, phasei)
419  {
420  phaseModel& phase = phaseModels_[phasei];
421  movingPhasei += !phase.stationary();
422  stationaryPhasei += phase.stationary();
423  anisothermalPhasei += !phase.isothermal();
424  multiComponentPhasei += !phase.pure();
425  }
426  movingPhaseModels_.resize(movingPhasei);
427  stationaryPhaseModels_.resize(stationaryPhasei);
428  anisothermalPhaseModels_.resize(anisothermalPhasei);
429  multiComponentPhaseModels_.resize(multiComponentPhasei);
430 
431  movingPhasei = 0;
432  stationaryPhasei = 0;
433  anisothermalPhasei = 0;
434  multiComponentPhasei = 0;
435  forAll(phaseModels_, phasei)
436  {
437  phaseModel& phase = phaseModels_[phasei];
438  if (!phase.stationary())
439  {
440  movingPhaseModels_.set(movingPhasei++, &phase);
441  }
442  if (phase.stationary())
443  {
444  stationaryPhaseModels_.set(stationaryPhasei++, &phase);
445  }
446  if (!phase.isothermal())
447  {
448  anisothermalPhaseModels_.set(anisothermalPhasei++, &phase);
449  }
450  if (!phase.pure())
451  {
452  multiComponentPhaseModels_.set(multiComponentPhasei++, &phase);
453  }
454  }
455 
456  // Write phi
457  phi_.writeOpt() = IOobject::AUTO_WRITE;
458 
459  // Blending methods
460  forAllConstIter(dictionary, subDict("blending"), iter)
461  {
462  blendingMethods_.insert
463  (
464  iter().keyword(),
466  (
467  iter().keyword(),
468  iter().dict(),
469  phaseModels_.toc()
470  )
471  );
472  }
473 
474  // Sub-models
475  generatePairsAndSubModels("surfaceTension", surfaceTensionModels_);
476  generatePairsAndSubModels("aspectRatio", aspectRatioModels_);
477 
478  // Update motion fields
479  correctKinematics();
480 
481  // Set the optional reference phase fraction from the other phases
482  if (referencePhaseName_ != word::null)
483  {
484  phaseModel* referencePhasePtr = &phases()[referencePhaseName_];
485  volScalarField& referenceAlpha = *referencePhasePtr;
486 
487  referenceAlpha = 1;
488 
489  forAll(phaseModels_, phasei)
490  {
491  if (&phaseModels_[phasei] != referencePhasePtr)
492  {
493  referenceAlpha -= phaseModels_[phasei];
494  }
495  }
496  }
497 
498  forAll(phases(), phasei)
499  {
500  const volScalarField& alphai = phases()[phasei];
501  mesh_.setFluxRequired(alphai.name());
502  }
503 }
504 
505 
506 // * * * * * * * * * * * * * * * * Destructor * * * * * * * * * * * * * * * //
507 
509 {}
510 
511 
512 // * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
513 
515 {
516  tmp<volScalarField> rho(movingPhaseModels_[0]*movingPhaseModels_[0].rho());
517 
518  for
519  (
520  label movingPhasei=1;
521  movingPhasei<movingPhaseModels_.size();
522  movingPhasei++
523  )
524  {
525  rho.ref() +=
526  movingPhaseModels_[movingPhasei]
527  *movingPhaseModels_[movingPhasei].rho();
528  }
529 
530  if (stationaryPhaseModels_.empty())
531  {
532  return rho;
533  }
534  else
535  {
536  return rho/sumAlphaMoving();
537  }
538 }
539 
540 
542 {
543  tmp<volVectorField> U(movingPhaseModels_[0]*movingPhaseModels_[0].U());
544 
545  for
546  (
547  label movingPhasei=1;
548  movingPhasei<movingPhaseModels_.size();
549  movingPhasei++
550  )
551  {
552  U.ref() +=
553  movingPhaseModels_[movingPhasei]
554  *movingPhaseModels_[movingPhasei].U();
555  }
556 
557  if (stationaryPhaseModels_.empty())
558  {
559  return U;
560  }
561  else
562  {
563  return U/sumAlphaMoving();
564  }
565 }
566 
567 
569 Foam::phaseSystem::E(const phasePairKey& key) const
570 {
571  if (aspectRatioModels_.found(key))
572  {
573  return aspectRatioModels_[key]->E();
574  }
575  else
576  {
577  return volScalarField::New
578  (
579  aspectRatioModel::typeName + ":E",
580  mesh_,
582  );
583  }
584 }
585 
586 
588 Foam::phaseSystem::sigma(const phasePairKey& key) const
589 {
590  if (surfaceTensionModels_.found(key))
591  {
592  return surfaceTensionModels_[key]->sigma();
593  }
594  else
595  {
596  return volScalarField::New
597  (
598  surfaceTensionModel::typeName + ":sigma",
599  mesh_,
601  );
602  }
603 }
604 
605 
607 Foam::phaseSystem::sigma(const phasePairKey& key, label patchi) const
608 {
609  if (surfaceTensionModels_.found(key))
610  {
611  return surfaceTensionModels_[key]->sigma(patchi);
612  }
613  else
614  {
615  return tmp<scalarField>
616  (
617  new scalarField(mesh_.boundary()[patchi].size(), 0)
618  );
619  }
620 }
621 
622 
625 {
626  tmp<volScalarField> tnearInt
627  (
629  (
630  "nearInterface",
631  mesh_,
633  )
634  );
635 
636  forAll(phases(), phasei)
637  {
638  tnearInt.ref() = max
639  (
640  tnearInt(),
641  pos0(phases()[phasei] - 0.01)*pos0(0.99 - phases()[phasei])
642  );
643  }
644 
645  return tnearInt;
646 }
647 
648 
650 (
651  const phasePairKey& key
652 ) const
653 {
654  const phasePair pair
655  (
656  phaseModels_[key.first()],
657  phaseModels_[key.second()]
658  );
659 
660  return volScalarField::New
661  (
662  IOobject::groupName("dmdtf", pair.name()),
663  mesh(),
665  );
666 }
667 
668 
670 {
671  return PtrList<volScalarField>(phaseModels_.size());
672 }
673 
674 
676 {
677  return PtrList<volScalarField>(phaseModels_.size());
678 }
679 
680 
682 {
683  forAll(phaseModels_, phasei)
684  {
685  if (!phaseModels_[phasei].incompressible())
686  {
687  return false;
688  }
689  }
690 
691  return true;
692 }
693 
694 
695 bool Foam::phaseSystem::implicitPhasePressure(const phaseModel& phase) const
696 {
697  return false;
698 }
699 
700 
702 {
703  return false;
704 }
705 
706 
708 (
709  const phaseModel& phase1
710 ) const
711 {
712  tmp<surfaceScalarField> tSurfaceTension
713  (
715  (
716  "surfaceTension",
717  mesh_,
718  dimensionedScalar(dimensionSet(1, -2, -2, 0, 0), 0)
719  )
720  );
721 
722  forAll(phases(), phasej)
723  {
724  const phaseModel& phase2 = phases()[phasej];
725 
726  if (&phase2 != &phase1)
727  {
728  phasePairKey key12(phase1.name(), phase2.name());
729 
730  cAlphaTable::const_iterator cAlpha(cAlphas_.find(key12));
731 
732  if (cAlpha != cAlphas_.end())
733  {
734  tSurfaceTension.ref() +=
735  fvc::interpolate(sigma(key12)*K(phase1, phase2))
736  *(
737  fvc::interpolate(phase2)*fvc::snGrad(phase1)
738  - fvc::interpolate(phase1)*fvc::snGrad(phase2)
739  );
740  }
741  }
742  }
743 
744  return tSurfaceTension;
745 }
746 
747 
749 {
750  forAll(phaseModels_, phasei)
751  {
752  phaseModels_[phasei].correct();
753  }
754 }
755 
756 
758 {
759  const PtrList<volScalarField> dmdts = this->dmdts();
760 
761  forAll(movingPhaseModels_, movingPhasei)
762  {
763  phaseModel& phase = movingPhaseModels_[movingPhasei];
764  const volScalarField& alpha = phase;
765  volScalarField& rho = phase.thermoRef().rho();
766 
768  (
770  (
771  IOobject::groupName("source", phase.name()),
772  mesh_,
774  )
775  );
776 
777  if (fvModels().addsSupToField(rho.name()))
778  {
779  source += fvModels().source(alpha, rho)&rho;
780  }
781 
782  if (dmdts.set(phase.index()))
783  {
784  source += dmdts[phase.index()];
785  }
786 
787  phase.correctContinuityError(source);
788  }
789 }
790 
791 
793 {
794  bool updateDpdt = false;
795 
796  forAll(phaseModels_, phasei)
797  {
798  phaseModels_[phasei].correctKinematics();
799 
800  updateDpdt = updateDpdt || phaseModels_[phasei].thermo().dpdt();
801  }
802 
803  // Update the pressure time-derivative if required
804  if (updateDpdt)
805  {
806  dpdt_ = fvc::ddt(phaseModels_.begin()().thermo().p());
807  }
808 }
809 
810 
812 {
813  forAll(phaseModels_, phasei)
814  {
815  phaseModels_[phasei].correctThermo();
816  }
817 }
818 
819 
821 {
822  forAll(phaseModels_, phasei)
823  {
824  phaseModels_[phasei].correctReactions();
825  }
826 }
827 
828 
830 {
831  forAll(phaseModels_, phasei)
832  {
833  phaseModels_[phasei].correctSpecies();
834  }
835 }
836 
837 
839 {
840  forAll(phaseModels_, phasei)
841  {
842  phaseModels_[phasei].correctTurbulence();
843  }
844 }
845 
846 
848 {
849  forAll(phaseModels_, phasei)
850  {
851  phaseModels_[phasei].correctEnergyTransport();
852  }
853 }
854 
855 
857 {
858  if (mesh_.changing())
859  {
860  MRF_.update();
861 
862  // forAll(phaseModels_, phasei)
863  // {
864  // phaseModels_[phasei].meshUpdate();
865  // }
866  }
867 }
868 
869 
871 {
872  forAll(movingPhases(), movingPhasei)
873  {
874  phaseModel& phase = movingPhases()[movingPhasei];
875 
876  const volVectorField::Boundary& UBf = phase.U()().boundaryField();
877 
878  FieldField<fvsPatchField, scalar> phiRelBf
879  (
880  MRF_.relative(mesh_.Sf().boundaryField() & UBf)
881  );
882 
883  surfaceScalarField::Boundary& phiBf = phase.phiRef().boundaryFieldRef();
884 
885  forAll(mesh_.boundary(), patchi)
886  {
887  if
888  (
889  isA<fixedValueFvsPatchScalarField>(phiBf[patchi])
890  && !isA<movingWallVelocityFvPatchVectorField>(UBf[patchi])
891  )
892  {
893  phiBf[patchi] == phiRelBf[patchi];
894  }
895  }
896  }
897 }
898 
899 
901 (
902  const volScalarField& p_rgh,
903  const tmp<volScalarField>& divU,
905  nonOrthogonalSolutionControl& pimple
906 )
907 {
908  forAll(movingPhases(), movingPhasei)
909  {
910  phaseModel& phase = movingPhases()[movingPhasei];
911 
912  volVectorField::Boundary& Ubf = phase.URef().boundaryFieldRef();
913  surfaceVectorField::Boundary& UfBf = phase.UfRef().boundaryFieldRef();
914 
915  forAll(Ubf, patchi)
916  {
917  if (Ubf[patchi].fixesValue())
918  {
919  Ubf[patchi].initEvaluate();
920  }
921  }
922 
923  forAll(Ubf, patchi)
924  {
925  if (Ubf[patchi].fixesValue())
926  {
927  Ubf[patchi].evaluate();
928  UfBf[patchi] = Ubf[patchi];
929  }
930  }
931  }
932 
933  // Correct fixed-flux BCs to be consistent with the velocity BCs
934  correctBoundaryFlux();
935 
936  {
937  phi_ = Zero;
938  PtrList<surfaceScalarField> alphafs(phaseModels_.size());
939  forAll(movingPhases(), movingPhasei)
940  {
941  phaseModel& phase = movingPhases()[movingPhasei];
942  const label phasei = phase.index();
943  const volScalarField& alpha = phase;
944 
945  alphafs.set(phasei, fvc::interpolate(alpha).ptr());
946 
947  // Calculate absolute flux
948  // from the mapped surface velocity
949  phi_ += alphafs[phasei]*(mesh_.Sf() & phase.Uf());
950  }
951 
952  CorrectPhi
953  (
954  phi_,
955  movingPhases()[0].U(),
956  p_rgh,
957  // surfaceScalarField("rAUf", fvc::interpolate(rAU())),
959  divU(),
960  pressureReference,
961  pimple
962  );
963 
964  // Make the flux relative to the mesh motion
965  fvc::makeRelative(phi_, movingPhases()[0].U());
966 
967  setMixturePhi(alphafs, phi_);
968  }
969 }
970 
971 
973 {
974  if (regIOobject::read())
975  {
976  bool readOK = true;
977 
978  forAll(phaseModels_, phasei)
979  {
980  readOK &= phaseModels_[phasei].read();
981  }
982 
983  // models ...
984 
985  return readOK;
986  }
987  else
988  {
989  return false;
990  }
991 }
992 
993 
995 {
996  if (fv::localEulerDdt::enabled(vf.mesh()))
997  {
998  return fv::localEulerDdt::localRDeltaT(vf.mesh())*vf;
999  }
1000  else
1001  {
1002  return vf/vf.mesh().time().deltaT();
1003  }
1004 }
1005 
1006 
1008 {
1009  if (fv::localEulerDdt::enabled(sf.mesh()))
1010  {
1011  return fv::localEulerDdt::localRDeltaTf(sf.mesh())*sf;
1012  }
1013  else
1014  {
1015  return sf/sf.mesh().time().deltaT();
1016  }
1017 }
1018 
1019 
1020 // ************************************************************************* //
static const word propertiesName
Default name of the phase properties dictionary.
Definition: phaseSystem.H:341
pressureReference & pressureReference
volScalarField & p_rgh
dimensionedScalar tanh(const dimensionedScalar &ds)
dictionary dict
void generatePairs(const dictTable &modelDicts)
Generate pairs.
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:434
tmp< volScalarField > E(const phasePairKey &key) const
Return the aspect-ratio for a pair.
virtual void correctTurbulence()
Correct the turbulence.
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
fluidReactionThermo & thermo
Definition: createFields.H:28
pimpleNoLoopControl & pimple
static const surfaceScalarField & localRDeltaTf(const fvMesh &mesh)
Return the reciprocal of the local face time-step.
Definition: localEulerDdt.C:58
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:124
label phasei
Definition: pEqn.H:27
tmp< surfaceScalarField > nHatf(const volScalarField &alpha1, const volScalarField &alpha2) const
Normal to interface between two phases dotted with face areas.
virtual PtrList< volScalarField > d2mdtdps() const
Return the mass transfer pressure implicit coefficients.
virtual bool read()
Read object.
error FatalError
dimensioned< Type > max(const dimensioned< Type > &, const dimensioned< Type > &)
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:323
PtrList< surfaceScalarField > alphafs(phases.size())
tmp< GeometricField< Type, fvPatchField, volMesh > > div(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcDiv.C:47
Unit conversion functions.
const dimensionSet dimPressure
void setMixturePhi(const PtrList< surfaceScalarField > &alphafs, const surfaceScalarField &phim)
Re-normalise the flux of the phases.
volScalarField alpha(IOobject("alpha", runTime.timeName(), mesh, IOobject::READ_IF_PRESENT, IOobject::AUTO_WRITE), lambda *max(Ua &U, zeroSensitivity))
static tmp< GeometricField< scalar, fvsPatchField, surfaceMesh > > New(const word &name, const Internal &, const PtrList< fvsPatchField< scalar >> &)
Return a temporary field constructed from name,.
GeometricField< vector, fvsPatchField, surfaceMesh > surfaceVectorField
const dimensionedScalar b
Wien displacement law constant: default SI units: [m K].
Definition: createFields.H:27
tmp< volScalarField > sigma(const phasePairKey &key) const
Return the surface tension coefficient for a pair.
tmp< volScalarField > K(const phaseModel &alpha1, const phaseModel &alpha2) const
Curvature of interface between two phases.
Calculate the snGrad of the given volField.
void correctPhi(const volScalarField &p_rgh, const tmp< volScalarField > &divU, const pressureReference &pressureReference, nonOrthogonalSolutionControl &pimple)
const dimensionSet dimless
dimensionedScalar det(const dimensionedSphericalTensor &dt)
GeometricField< vector, fvPatchField, volMesh > volVectorField
Definition: volFieldsFwd.H:58
CGAL::Exact_predicates_exact_constructions_kernel K
tmp< volScalarField > nearInterface() const
Indicator of the proximity of the interface.
void correctBoundaryFlux()
Correct fixed-flux BCs to be consistent with the velocity BCs.
scalar degToRad(const scalar deg)
Conversion from degrees to radians.
void CorrectPhi(surfaceScalarField &phi, const volVectorField &U, const volScalarField &p, const RAUfType &rAUf, const DivUType &divU, const pressureReference &pressureReference, nonOrthogonalSolutionControl &pcorrControl)
tmp< GeometricField< Type, fvPatchField, volMesh > > ddt(const dimensioned< Type > dt, const fvMesh &mesh)
Definition: fvcDdt.C:45
tmp< surfaceScalarField > surfaceTension(const phaseModel &phase) const
Return the surface tension force.
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:57
const dimensionedScalar sigma
Stefan-Boltzmann constant: default SI units: [W/m^2/K^4].
Calculate the first temporal derivative.
virtual void correctContinuityError()
Correct the continuity errors.
const dimensionSet dimTime
virtual void correctKinematics()
Correct the kinematics.
stressControl lookup("compactNormalStress") >> compactNormalStress
dynamicFvMesh & mesh
static bool enabled(const fvMesh &mesh)
Return true if LTS is enabled.
Definition: localEulerDdt.C:37
dimensionedScalar cos(const dimensionedScalar &ds)
phaseSystem(const fvMesh &mesh)
Construct from fvMesh.
Calculate the gradient of the given field.
A class for handling words, derived from string.
Definition: word.H:59
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
static word groupName(Name name, const word &group)
tmp< volVectorField > U() const
Return the mixture velocity.
static autoPtr< blendingMethod > New(const word &modelName, const dictionary &dict, const wordList &phaseNames)
const dimensionSet dimDensity
virtual void correctThermo()
Correct the thermodynamics.
const Type & value() const
Return const reference to value.
static const word null
An empty word.
Definition: word.H:77
phi
Definition: correctPhi.H:3
word timeName
Definition: getTimeIndex.H:3
static const zero Zero
Definition: zero.H:97
static const volScalarField & localRDeltaT(const fvMesh &mesh)
Return the reciprocal of the local time-step.
Definition: localEulerDdt.C:46
dimensioned< Type > average(const DimensionedField< Type, GeoMesh > &df)
faceListList boundary(nPatches)
forAllConstIter(PtrDictionary< phaseModel >, mixture.phases(), phase)
Definition: pEqn.H:29
tmp< surfaceVectorField > nHatfv(const volScalarField &alpha1, const volScalarField &alpha2) const
Normal to interface between two phases.
Calculate the divergence of the given field.
dimensionedScalar pos0(const dimensionedScalar &ds)
virtual void correct()
Correct the fluid properties other than those listed below.
tmp< volScalarField > rho() const
Return the mixture density.
virtual PtrList< volScalarField > dmdts() const
Return the mass transfer rates for each phase.
defineRunTimeSelectionTable(reactionRateFlameArea, dictionary)
tmp< volScalarField > sumAlphaMoving() const
Return the sum of the phase fractions of the moving phases.
const Mesh & mesh() const
Return mesh.
defineTypeNameAndDebug(combustionModel, 0)
tmp< volScalarField > divU
Definition: alphaSuSp.H:9
volScalarField sf(fieldObject, mesh)
bool incompressible() const
Return incompressibility.
tmp< volScalarField > byDt(const volScalarField &vf)
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.
virtual void correctReactions()
Correct the reactions.
Foam::fvModels & fvModels
dimensionedScalar pow(const dimensionedScalar &ds, const dimensionedScalar &expt)
virtual void correctEnergyTransport()
Correct the energy transport e.g. alphat.
label patchi
tmp< surfaceScalarField > calcPhi(const phaseModelList &phaseModels) const
Calculate and return the mixture flux.
U
Definition: pEqn.H:72
fvModels source(alpha1, mixture.thermo1().rho())
static const dimensionSet dimSigma
Surface tension coefficient dimensions.
virtual ~phaseSystem()
Destructor.
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
A templated 1D list of pointers to objects of type <T>, where the size of the array is known and used...
Definition: List.H:70
virtual bool read()
Read base phaseProperties dictionary.
phasePairKey()
Construct null.
phaseSystem::phaseModelList & phases
Definition: createFields.H:12
dimensioned< scalar > mag(const dimensioned< Type > &)
Field< vector > vectorField
Specialisation of Field<T> for vector.
void setMixtureU(const volVectorField &Um)
Re-normalise the velocity of the phases.
const doubleScalar e
Elementary charge.
Definition: doubleScalar.H:105
void correctContactAngle(const phaseModel &alpha1, const phaseModel &alpha2, surfaceVectorField::Boundary &nHatb) const
Correction for the boundary condition on the unit normal nHat on.
A class for managing temporary objects.
Definition: PtrList.H:53
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
void makeRelative(surfaceScalarField &phi, const volVectorField &U)
Make the given flux relative.
Definition: fvcMeshPhi.C:75
virtual void meshUpdate()
Update the fluid properties for mesh changes.
virtual bool implicitPhasePressure() const
Returns true if the phase pressure is treated implicitly.
tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > snGrad(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvcSnGrad.C:45
virtual void correctSpecies()
Correct the species mass fractions.
Namespace for OpenFOAM.
virtual tmp< volScalarField > dmdtf(const phasePairKey &key) const
Return the mass transfer rate for an interface.