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 | www.openfoam.com
6 \\/ M anipulation |
7-------------------------------------------------------------------------------
8 Copyright (C) 2013-2017 OpenFOAM Foundation
9 Copyright (C) 2020 OpenCFD Ltd.
10-------------------------------------------------------------------------------
11License
12 This file is part of OpenFOAM.
13
14 OpenFOAM is free software: you can redistribute it and/or modify it
15 under the terms of the GNU General Public License as published by
16 the Free Software Foundation, either version 3 of the License, or
17 (at your option) any later version.
18
19 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
20 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
21 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
22 for more details.
23
24 You should have received a copy of the GNU General Public License
25 along with OpenFOAM. If not, see <http://www.gnu.org/licenses/>.
26
27\*---------------------------------------------------------------------------*/
28
30#include "alphaContactAngleFvPatchScalarField.H"
31#include "Time.H"
32#include "subCycle.H"
33#include "MULES.H"
34#include "fvcDiv.H"
35#include "fvcGrad.H"
36#include "fvcSnGrad.H"
37#include "fvcFlux.H"
38#include "fvcMeshPhi.H"
39#include "surfaceInterpolate.H"
40#include "unitConversion.H"
41
42// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
43
44namespace Foam
45{
46 defineTypeNameAndDebug(multiphaseMixtureThermo, 0);
47}
48
49
50// * * * * * * * * * * * * * Private Member Functions * * * * * * * * * * * //
51
52void Foam::multiphaseMixtureThermo::calcAlphas()
53{
54 scalar level = 0.0;
55 alphas_ == 0.0;
56
57 for (const phaseModel& phase : phases_)
58 {
59 alphas_ += level * phase;
60 level += 1.0;
61 }
62}
63
64
65// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
66
68(
69 const volVectorField& U,
71)
72:
73 psiThermo(U.mesh(), word::null),
74 phases_(lookup("phases"), phaseModel::iNew(p_, T_)),
75
76 mesh_(U.mesh()),
77 U_(U),
78 phi_(phi),
79
80 rhoPhi_
81 (
82 IOobject
83 (
84 "rhoPhi",
85 mesh_.time().timeName(),
86 mesh_,
87 IOobject::NO_READ,
88 IOobject::NO_WRITE
89 ),
90 mesh_,
92 ),
93
94 alphas_
95 (
96 IOobject
97 (
98 "alphas",
99 mesh_.time().timeName(),
100 mesh_,
101 IOobject::NO_READ,
102 IOobject::AUTO_WRITE
103 ),
104 mesh_,
106 ),
107
108 sigmas_(lookup("sigmas")),
109 dimSigma_(1, 0, -2, 0, 0),
110 deltaN_
111 (
112 "deltaN",
113 1e-8/cbrt(average(mesh_.V()))
114 )
115{
116 rhoPhi_.setOriented();
117 calcAlphas();
118 alphas_.write();
119 correct();
120}
121
122
123// * * * * * * * * * * * * * * Member Functions * * * * * * * * * * * * * * //
124
126{
127 for (phaseModel& phase : phases_)
128 {
129 phase.correct();
130 }
131
132 auto phasei = phases_.cbegin();
133
134 psi_ = phasei()*phasei().thermo().psi();
135 mu_ = phasei()*phasei().thermo().mu();
136 alpha_ = phasei()*phasei().thermo().alpha();
137
138 for (++phasei; phasei != phases_.cend(); ++phasei)
139 {
140 psi_ += phasei()*phasei().thermo().psi();
141 mu_ += phasei()*phasei().thermo().mu();
142 alpha_ += phasei()*phasei().thermo().alpha();
143 }
144}
145
146
148{
149 for (phaseModel& phase : phases_)
150 {
151 phase.thermo().rho() += phase.thermo().psi()*dp;
152 }
153}
154
155
157{
158 auto phasei = phases_.cbegin();
159
160 word name = phasei().thermo().thermoName();
161
162 for (++phasei; phasei != phases_.cend(); ++phasei)
163 {
164 name += ',' + phasei().thermo().thermoName();
165 }
166
167 return name;
168}
169
170
172{
173 for (const phaseModel& phase : phases_)
174 {
175 if (!phase.thermo().incompressible())
176 {
177 return false;
178 }
179 }
180
181 return true;
182}
183
184
186{
187 for (const phaseModel& phase : phases_)
188 {
189 if (!phase.thermo().isochoric())
190 {
191 return false;
192 }
193 }
194
195 return true;
196}
197
198
200(
201 const volScalarField& p,
202 const volScalarField& T
203) const
204{
205 auto phasei = phases_.cbegin();
206
207 tmp<volScalarField> the(phasei()*phasei().thermo().he(p, T));
208
209 for (++phasei; phasei != phases_.cend(); ++phasei)
210 {
211 the.ref() += phasei()*phasei().thermo().he(p, T);
212 }
213
214 return the;
215}
216
217
219(
220 const scalarField& p,
221 const scalarField& T,
222 const labelList& cells
223) const
224{
225 auto phasei = phases_.cbegin();
226
227 tmp<scalarField> the
228 (
230 );
231
232 for (++phasei; phasei != phases_.cend(); ++phasei)
233 {
234 the.ref() +=
235 scalarField(phasei(), cells)*phasei().thermo().he(p, T, cells);
236 }
237
238 return the;
239}
240
241
243(
244 const scalarField& p,
245 const scalarField& T,
246 const label patchi
247) const
248{
249 auto phasei = phases_.cbegin();
250
251 tmp<scalarField> the
252 (
253 phasei().boundaryField()[patchi]*phasei().thermo().he(p, T, patchi)
254 );
255
256 for (++phasei; phasei != phases_.cend(); ++phasei)
257 {
258 the.ref() +=
259 phasei().boundaryField()[patchi]*phasei().thermo().he(p, T, patchi);
260 }
261
262 return the;
263}
264
265
267{
268 auto phasei = phases_.cbegin();
269
270 tmp<volScalarField> thc(phasei()*phasei().thermo().hc());
271
272 for (++phasei; phasei != phases_.cend(); ++phasei)
273 {
274 thc.ref() += phasei()*phasei().thermo().hc();
275 }
276
277 return thc;
278}
279
280
282(
283 const scalarField& h,
284 const scalarField& p,
285 const scalarField& T0,
286 const labelList& cells
287) const
288{
290 return T0;
291}
292
293
295(
296 const scalarField& h,
297 const scalarField& p,
298 const scalarField& T0,
299 const label patchi
300) const
301{
303 return T0;
304}
305
306
308{
309 auto phasei = phases_.cbegin();
310
311 tmp<volScalarField> trho(phasei()*phasei().thermo().rho());
312
313 for (++phasei; phasei != phases_.cend(); ++phasei)
314 {
315 trho.ref() += phasei()*phasei().thermo().rho();
316 }
317
318 return trho;
319}
320
321
323(
324 const label patchi
325) const
326{
327 auto phasei = phases_.cbegin();
328
329 tmp<scalarField> trho
330 (
331 phasei().boundaryField()[patchi]*phasei().thermo().rho(patchi)
332 );
333
334 for (++phasei; phasei != phases_.cend(); ++phasei)
335 {
336 trho.ref() +=
337 phasei().boundaryField()[patchi]*phasei().thermo().rho(patchi);
338 }
339
340 return trho;
341}
342
343
345{
346 auto phasei = phases_.cbegin();
347
348 tmp<volScalarField> tCp(phasei()*phasei().thermo().Cp());
349
350 for (++phasei; phasei != phases_.cend(); ++phasei)
351 {
352 tCp.ref() += phasei()*phasei().thermo().Cp();
353 }
354
355 return tCp;
356}
357
358
360(
361 const scalarField& p,
362 const scalarField& T,
363 const label patchi
364) const
365{
366 auto phasei = phases_.cbegin();
367
368 tmp<scalarField> tCp
369 (
370 phasei().boundaryField()[patchi]*phasei().thermo().Cp(p, T, patchi)
371 );
372
373 for (++phasei; phasei != phases_.cend(); ++phasei)
374 {
375 tCp.ref() +=
376 phasei().boundaryField()[patchi]*phasei().thermo().Cp(p, T, patchi);
377 }
378
379 return tCp;
380}
381
382
384{
385 auto phasei = phases_.cbegin();
386
387 tmp<volScalarField> tCv(phasei()*phasei().thermo().Cv());
388
389 for (++phasei; phasei != phases_.cend(); ++phasei)
390 {
391 tCv.ref() += phasei()*phasei().thermo().Cv();
392 }
393
394 return tCv;
395}
396
397
399(
400 const scalarField& p,
401 const scalarField& T,
402 const label patchi
403) const
404{
405 auto phasei = phases_.cbegin();
406
407 tmp<scalarField> tCv
408 (
409 phasei().boundaryField()[patchi]*phasei().thermo().Cv(p, T, patchi)
410 );
411
412 for (++phasei; phasei != phases_.cend(); ++phasei)
413 {
414 tCv.ref() +=
415 phasei().boundaryField()[patchi]*phasei().thermo().Cv(p, T, patchi);
416 }
417
418 return tCv;
419}
420
421
423{
424 auto phasei = phases_.cbegin();
425
426 tmp<volScalarField> tgamma(phasei()*phasei().thermo().gamma());
427
428 for (++phasei; phasei != phases_.cend(); ++phasei)
429 {
430 tgamma.ref() += phasei()*phasei().thermo().gamma();
431 }
432
433 return tgamma;
434}
435
436
438(
439 const scalarField& p,
440 const scalarField& T,
441 const label patchi
442) const
443{
444 auto phasei = phases_.cbegin();
445
446 tmp<scalarField> tgamma
447 (
448 phasei().boundaryField()[patchi]*phasei().thermo().gamma(p, T, patchi)
449 );
450
451 for (++phasei; phasei != phases_.cend(); ++phasei)
452 {
453 tgamma.ref() +=
454 phasei().boundaryField()[patchi]
455 *phasei().thermo().gamma(p, T, patchi);
456 }
457
458 return tgamma;
459}
460
461
463{
464 auto phasei = phases_.cbegin();
465
466 tmp<volScalarField> tCpv(phasei()*phasei().thermo().Cpv());
467
468 for (++phasei; phasei != phases_.cend(); ++phasei)
469 {
470 tCpv.ref() += phasei()*phasei().thermo().Cpv();
471 }
472
473 return tCpv;
474}
475
476
478(
479 const scalarField& p,
480 const scalarField& T,
481 const label patchi
482) const
483{
484 auto phasei = phases_.cbegin();
485
486 tmp<scalarField> tCpv
487 (
488 phasei().boundaryField()[patchi]*phasei().thermo().Cpv(p, T, patchi)
489 );
490
491 for (++phasei; phasei != phases_.cend(); ++phasei)
492 {
493 tCpv.ref() +=
494 phasei().boundaryField()[patchi]
495 *phasei().thermo().Cpv(p, T, patchi);
496 }
497
498 return tCpv;
499}
500
501
503{
504 auto phasei = phases_.cbegin();
505
506 tmp<volScalarField> tCpByCpv(phasei()*phasei().thermo().CpByCpv());
507
508 for (++phasei; phasei != phases_.cend(); ++phasei)
509 {
510 tCpByCpv.ref() += phasei()*phasei().thermo().CpByCpv();
511 }
512
513 return tCpByCpv;
514}
515
516
518(
519 const scalarField& p,
520 const scalarField& T,
521 const label patchi
522) const
523{
524 auto phasei = phases_.cbegin();
525
526 tmp<scalarField> tCpByCpv
527 (
528 phasei().boundaryField()[patchi]*phasei().thermo().CpByCpv(p, T, patchi)
529 );
530
531 for (++phasei; phasei != phases_.cend(); ++phasei)
532 {
533 tCpByCpv.ref() +=
534 phasei().boundaryField()[patchi]
535 *phasei().thermo().CpByCpv(p, T, patchi);
536 }
537
538 return tCpByCpv;
539}
540
541
543{
544 auto phasei = phases_.cbegin();
545
546 tmp<volScalarField> tW(phasei()*phasei().thermo().W());
547
548 for (++phasei; phasei != phases_.cend(); ++phasei)
549 {
550 tW.ref() += phasei()*phasei().thermo().W();
551 }
552
553 return tW;
554}
555
556
558{
559 return mu()/rho();
560}
561
562
564(
565 const label patchi
566) const
567{
568 return mu(patchi)/rho(patchi);
569}
570
571
573{
574 auto phasei = phases_.cbegin();
575
576 tmp<volScalarField> tkappa(phasei()*phasei().thermo().kappa());
577
578 for (++phasei; phasei != phases_.cend(); ++phasei)
579 {
580 tkappa.ref() += phasei()*phasei().thermo().kappa();
581 }
582
583 return tkappa;
584}
585
586
588(
589 const label patchi
590) const
591{
592 auto phasei = phases_.cbegin();
593
594 tmp<scalarField> tkappa
595 (
596 phasei().boundaryField()[patchi]*phasei().thermo().kappa(patchi)
597 );
598
599 for (++phasei; phasei != phases_.cend(); ++phasei)
600 {
601 tkappa.ref() +=
602 phasei().boundaryField()[patchi]*phasei().thermo().kappa(patchi);
603 }
604
605 return tkappa;
606}
607
608
610{
611 PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
612
613 tmp<volScalarField> talphaEff(phasei()*phasei().thermo().alphahe());
614
615 for (++phasei; phasei != phases_.end(); ++phasei)
616 {
617 talphaEff.ref() += phasei()*phasei().thermo().alphahe();
618 }
619
620 return talphaEff;
621}
622
623
625(
626 const label patchi
627) const
628{
629 PtrDictionary<phaseModel>::const_iterator phasei = phases_.begin();
630
631 tmp<scalarField> talphaEff
632 (
633 phasei().boundaryField()[patchi]
634 *phasei().thermo().alphahe(patchi)
635 );
636
637 for (++phasei; phasei != phases_.end(); ++phasei)
638 {
639 talphaEff.ref() +=
640 phasei().boundaryField()[patchi]
641 *phasei().thermo().alphahe(patchi);
642 }
643
644 return talphaEff;
645}
646
647
649(
650 const volScalarField& alphat
651) const
652{
653 auto phasei = phases_.cbegin();
654
655 tmp<volScalarField> tkappaEff(phasei()*phasei().thermo().kappaEff(alphat));
656
657 for (++phasei; phasei != phases_.cend(); ++phasei)
658 {
659 tkappaEff.ref() += phasei()*phasei().thermo().kappaEff(alphat);
660 }
661
662 return tkappaEff;
663}
664
665
667(
668 const scalarField& alphat,
669 const label patchi
670) const
671{
672 auto phasei = phases_.cbegin();
673
674 tmp<scalarField> tkappaEff
675 (
676 phasei().boundaryField()[patchi]
677 *phasei().thermo().kappaEff(alphat, patchi)
678 );
679
680 for (++phasei; phasei != phases_.cend(); ++phasei)
681 {
682 tkappaEff.ref() +=
683 phasei().boundaryField()[patchi]
684 *phasei().thermo().kappaEff(alphat, patchi);
685 }
686
687 return tkappaEff;
688}
689
690
692(
693 const volScalarField& alphat
694) const
695{
696 auto phasei = phases_.cbegin();
697
698 tmp<volScalarField> talphaEff(phasei()*phasei().thermo().alphaEff(alphat));
699
700 for (++phasei; phasei != phases_.cend(); ++phasei)
701 {
702 talphaEff.ref() += phasei()*phasei().thermo().alphaEff(alphat);
703 }
704
705 return talphaEff;
706}
707
708
710(
711 const scalarField& alphat,
712 const label patchi
713) const
714{
715 auto phasei = phases_.cbegin();
716
717 tmp<scalarField> talphaEff
718 (
719 phasei().boundaryField()[patchi]
720 *phasei().thermo().alphaEff(alphat, patchi)
721 );
722
723 for (++phasei; phasei != phases_.cend(); ++phasei)
724 {
725 talphaEff.ref() +=
726 phasei().boundaryField()[patchi]
727 *phasei().thermo().alphaEff(alphat, patchi);
728 }
729
730 return talphaEff;
731}
732
733
735{
736 auto phasei = phases_.cbegin();
737
738 tmp<volScalarField> trCv(phasei()/phasei().thermo().Cv());
739
740 for (++phasei; phasei != phases_.cend(); ++phasei)
741 {
742 trCv.ref() += phasei()/phasei().thermo().Cv();
743 }
744
745 return trCv;
746}
747
748
751{
752 tmp<surfaceScalarField> tstf
753 (
755 (
756 IOobject
757 (
758 "surfaceTensionForce",
759 mesh_.time().timeName(),
760 mesh_
761 ),
762 mesh_,
763 dimensionedScalar(dimensionSet(1, -2, -2, 0, 0), Zero)
764 )
765 );
766
767 surfaceScalarField& stf = tstf.ref();
768 stf.setOriented();
769
770 forAllConstIters(phases_, phase1)
771 {
772 const phaseModel& alpha1 = *phase1;
773
774 auto phase2 = phase1;
775
776 for (++phase2; phase2 != phases_.cend(); ++phase2)
777 {
778 const phaseModel& alpha2 = *phase2;
779
780 auto sigma = sigmas_.cfind(interfacePair(alpha1, alpha2));
781
782 if (!sigma.found())
783 {
785 << "Cannot find interface " << interfacePair(alpha1, alpha2)
786 << " in list of sigma values"
787 << exit(FatalError);
788 }
789
790 stf += dimensionedScalar("sigma", dimSigma_, *sigma)
792 (
795 );
796 }
797 }
798
799 return tstf;
800}
801
802
804{
805 const Time& runTime = mesh_.time();
806
807 const dictionary& alphaControls = mesh_.solverDict("alpha");
808 label nAlphaSubCycles(alphaControls.get<label>("nAlphaSubCycles"));
809 scalar cAlpha(alphaControls.get<scalar>("cAlpha"));
810
811 volScalarField& alpha = phases_.first();
812
813 if (nAlphaSubCycles > 1)
814 {
815 surfaceScalarField rhoPhiSum(0.0*rhoPhi_);
816 dimensionedScalar totalDeltaT = runTime.deltaT();
817
818 for
819 (
820 subCycle<volScalarField> alphaSubCycle(alpha, nAlphaSubCycles);
821 !(++alphaSubCycle).end();
822 )
823 {
824 solveAlphas(cAlpha);
825 rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi_;
826 }
827
828 rhoPhi_ = rhoPhiSum;
829 }
830 else
831 {
832 solveAlphas(cAlpha);
833 }
834}
835
836
837Foam::tmp<Foam::surfaceVectorField> Foam::multiphaseMixtureThermo::nHatfv
838(
839 const volScalarField& alpha1,
841) const
842{
843 /*
844 // Cell gradient of alpha
845 volVectorField gradAlpha =
846 alpha2*fvc::grad(alpha1) - alpha1*fvc::grad(alpha2);
847
848 // Interpolated face-gradient of alpha
849 surfaceVectorField gradAlphaf = fvc::interpolate(gradAlpha);
850 */
851
852 surfaceVectorField gradAlphaf
853 (
856 );
857
858 // Face unit interface normal
859 return gradAlphaf/(mag(gradAlphaf) + deltaN_);
860}
861
862
864(
865 const volScalarField& alpha1,
867) const
868{
869 // Face unit interface normal flux
870 return nHatfv(alpha1, alpha2) & mesh_.Sf();
871}
872
873
874// Correction for the boundary condition on the unit normal nHat on
875// walls to produce the correct contact angle.
876
877// The dynamic contact angle is calculated from the component of the
878// velocity on the direction of the interface, parallel to the wall.
879
880void Foam::multiphaseMixtureThermo::correctContactAngle
881(
882 const phaseModel& alpha1,
883 const phaseModel& alpha2,
885) const
886{
887 const volScalarField::Boundary& gbf
889
890 const fvBoundaryMesh& boundary = mesh_.boundary();
891
892 forAll(boundary, patchi)
893 {
894 if (isA<alphaContactAngleFvPatchScalarField>(gbf[patchi]))
895 {
896 const alphaContactAngleFvPatchScalarField& acap =
897 refCast<const alphaContactAngleFvPatchScalarField>(gbf[patchi]);
898
899 vectorField& nHatPatch = nHatb[patchi];
900
901 vectorField AfHatPatch
902 (
903 mesh_.Sf().boundaryField()[patchi]
904 /mesh_.magSf().boundaryField()[patchi]
905 );
906
907 const auto tp =
908 acap.thetaProps().cfind(interfacePair(alpha1, alpha2));
909
910 if (!tp.found())
911 {
913 << "Cannot find interface " << interfacePair(alpha1, alpha2)
914 << "\n in table of theta properties for patch "
915 << acap.patch().name()
916 << exit(FatalError);
917 }
918
919 const bool matched = (tp.key().first() == alpha1.name());
920
921 const scalar theta0 = degToRad(tp().theta0(matched));
922 scalarField theta(boundary[patchi].size(), theta0);
923
924 const scalar uTheta = tp().uTheta();
925
926 // Calculate the dynamic contact angle if required
927 if (uTheta > SMALL)
928 {
929 const scalar thetaA = degToRad(tp().thetaA(matched));
930 const scalar thetaR = degToRad(tp().thetaR(matched));
931
932 // Calculated the component of the velocity parallel to the wall
933 vectorField Uwall
934 (
935 U_.boundaryField()[patchi].patchInternalField()
936 - U_.boundaryField()[patchi]
937 );
938 Uwall -= (AfHatPatch & Uwall)*AfHatPatch;
939
940 // Find the direction of the interface parallel to the wall
941 vectorField nWall
942 (
943 nHatPatch - (AfHatPatch & nHatPatch)*AfHatPatch
944 );
945
946 // Normalise nWall
947 nWall /= (mag(nWall) + SMALL);
948
949 // Calculate Uwall resolved normal to the interface parallel to
950 // the interface
951 scalarField uwall(nWall & Uwall);
952
953 theta += (thetaA - thetaR)*tanh(uwall/uTheta);
954 }
955
956
957 // Reset nHatPatch to correspond to the contact angle
958
959 scalarField a12(nHatPatch & AfHatPatch);
960
961 scalarField b1(cos(theta));
962
963 scalarField b2(nHatPatch.size());
964
965 forAll(b2, facei)
966 {
967 b2[facei] = cos(acos(a12[facei]) - theta[facei]);
968 }
969
970 scalarField det(1.0 - a12*a12);
971
972 scalarField a((b1 - a12*b2)/det);
973 scalarField b((b2 - a12*b1)/det);
974
975 nHatPatch = a*AfHatPatch + b*nHatPatch;
976
977 nHatPatch /= (mag(nHatPatch) + deltaN_.value());
978 }
979 }
980}
981
982
984(
985 const phaseModel& alpha1,
986 const phaseModel& alpha2
987) const
988{
989 tmp<surfaceVectorField> tnHatfv = nHatfv(alpha1, alpha2);
990
991 correctContactAngle(alpha1, alpha2, tnHatfv.ref().boundaryFieldRef());
992
993 // Simple expression for curvature
994 return -fvc::div(tnHatfv & mesh_.Sf());
995}
996
997
1000{
1001 tmp<volScalarField> tnearInt
1002 (
1003 new volScalarField
1004 (
1005 IOobject
1006 (
1007 "nearInterface",
1008 mesh_.time().timeName(),
1009 mesh_
1010 ),
1011 mesh_,
1013 )
1014 );
1015
1016 for (const phaseModel& phase : phases_)
1017 {
1018 tnearInt.ref() =
1019 max(tnearInt(), pos0(phase - 0.01)*pos0(0.99 - phase));
1020 }
1021
1022 return tnearInt;
1023}
1024
1025
1026void Foam::multiphaseMixtureThermo::solveAlphas
1027(
1028 const scalar cAlpha
1029)
1030{
1031 static label nSolves(-1);
1032 ++nSolves;
1033
1034 const word alphaScheme("div(phi,alpha)");
1035 const word alpharScheme("div(phirb,alpha)");
1036
1037 surfaceScalarField phic(mag(phi_/mesh_.magSf()));
1038 phic = min(cAlpha*phic, max(phic));
1039
1040 PtrList<surfaceScalarField> alphaPhiCorrs(phases_.size());
1041
1042 int phasei = 0;
1043 for (phaseModel& alpha : phases_)
1044 {
1045 alphaPhiCorrs.set
1046 (
1047 phasei,
1049 (
1050 phi_.name() + alpha.name(),
1051 fvc::flux
1052 (
1053 phi_,
1054 alpha,
1055 alphaScheme
1056 )
1057 )
1058 );
1059
1060 surfaceScalarField& alphaPhiCorr = alphaPhiCorrs[phasei];
1061
1062 for (phaseModel& alpha2 : phases_)
1063 {
1064 if (&alpha2 == &alpha) continue;
1065
1067
1068 alphaPhiCorr += fvc::flux
1069 (
1071 alpha,
1073 );
1074 }
1075
1077 (
1078 1.0/mesh_.time().deltaT().value(),
1079 geometricOneField(),
1080 alpha,
1081 phi_,
1082 alphaPhiCorr,
1083 zeroField(),
1084 zeroField(),
1085 oneField(),
1086 zeroField(),
1087 true
1088 );
1089
1090 ++phasei;
1091 }
1092
1093 MULES::limitSum(alphaPhiCorrs);
1094
1095 rhoPhi_ = dimensionedScalar(dimensionSet(1, 0, -1, 0, 0), Zero);
1096
1097 volScalarField sumAlpha
1098 (
1099 IOobject
1100 (
1101 "sumAlpha",
1102 mesh_.time().timeName(),
1103 mesh_
1104 ),
1105 mesh_,
1107 );
1108
1109
1111
1112
1113 phasei = 0;
1114
1115 for (phaseModel& alpha : phases_)
1116 {
1117 surfaceScalarField& alphaPhi = alphaPhiCorrs[phasei];
1118 alphaPhi += upwind<scalar>(mesh_, phi_).flux(alpha);
1119
1121 (
1122 IOobject
1123 (
1124 "Sp",
1125 mesh_.time().timeName(),
1126 mesh_
1127 ),
1128 mesh_,
1130 );
1131
1133 (
1134 IOobject
1135 (
1136 "Su",
1137 mesh_.time().timeName(),
1138 mesh_
1139 ),
1140 // Divergence term is handled explicitly to be
1141 // consistent with the explicit transport solution
1142 divU*min(alpha, scalar(1))
1143 );
1144
1145 {
1146 const scalarField& dgdt = alpha.dgdt();
1147
1148 forAll(dgdt, celli)
1149 {
1150 if (dgdt[celli] < 0.0 && alpha[celli] > 0.0)
1151 {
1152 Sp[celli] += dgdt[celli]*alpha[celli];
1153 Su[celli] -= dgdt[celli]*alpha[celli];
1154 }
1155 else if (dgdt[celli] > 0.0 && alpha[celli] < 1.0)
1156 {
1157 Sp[celli] -= dgdt[celli]*(1.0 - alpha[celli]);
1158 }
1159 }
1160 }
1161
1162 for (const phaseModel& alpha2 : phases_)
1163 {
1164 if (&alpha2 == &alpha) continue;
1165
1166 const scalarField& dgdt2 = alpha2.dgdt();
1167
1168 forAll(dgdt2, celli)
1169 {
1170 if (dgdt2[celli] > 0.0 && alpha2[celli] < 1.0)
1171 {
1172 Sp[celli] -= dgdt2[celli]*(1.0 - alpha2[celli]);
1173 Su[celli] += dgdt2[celli]*alpha[celli];
1174 }
1175 else if (dgdt2[celli] < 0.0 && alpha2[celli] > 0.0)
1176 {
1177 Sp[celli] += dgdt2[celli]*alpha2[celli];
1178 }
1179 }
1180 }
1181
1183 (
1184 geometricOneField(),
1185 alpha,
1186 alphaPhi,
1187 Sp,
1188 Su
1189 );
1190
1191 rhoPhi_ += fvc::interpolate(alpha.thermo().rho())*alphaPhi;
1192
1193 Info<< alpha.name() << " volume fraction, min, max = "
1194 << alpha.weightedAverage(mesh_.V()).value()
1195 << ' ' << min(alpha).value()
1196 << ' ' << max(alpha).value()
1197 << endl;
1198
1199 sumAlpha += alpha;
1200
1201 ++phasei;
1202 }
1203
1204 Info<< "Phase-sum volume fraction, min, max = "
1205 << sumAlpha.weightedAverage(mesh_.V()).value()
1206 << ' ' << min(sumAlpha).value()
1207 << ' ' << max(sumAlpha).value()
1208 << endl;
1209
1210 calcAlphas();
1211}
1212
1213
1214// ************************************************************************* //
CGAL::Exact_predicates_exact_constructions_kernel K
surfaceScalarField phic(mixture.cAlpha() *mag(alphaPhic/mesh.magSf()))
MULES: Multidimensional universal limiter for explicit solution.
volScalarField & he
Definition: YEEqn.H:52
surfaceScalarField & phi
const volScalarField & alpha1
phaseModel & phase1
const volScalarField & alpha2
phaseModel & phase2
const dimensionSet & dimensions() const
Return dimensions.
dimensioned< Type > weightedAverage(const DimensionedField< scalar, GeoMesh > &weightField) const
Calculate and return weighted average.
void setOriented(const bool oriented=true) noexcept
Set the oriented flag.
DimensionedField< scalar, volMesh > Internal
The internal field type from which this GeometricField is derived.
Internal & ref(const bool updateAccessTime=true)
Return a reference to the dimensioned internal field.
GeometricBoundaryField< vector, fvsPatchField, surfaceMesh > Boundary
Type of boundary fields.
const Boundary & boundaryField() const
Return const-reference to the boundary field.
const word & name() const noexcept
Return the object name.
Definition: IOobjectI.H:65
dimensionedScalar deltaT() const
Return time step.
Definition: TimeStateI.H:55
T & first()
Return the first element of the list.
Definition: UListI.H:202
const_iterator cend() const noexcept
Return const_iterator to end traversing the constant UList.
Definition: UListI.H:364
volScalarField alpha_
Laminar thermal diffusivity [kg/m/s].
Definition: basicThermo.H:142
T get(const word &keyword, enum keyType::option matchOpt=keyType::REGEX) const
const surfaceScalarField & nHatf() const
tmp< volScalarField > alphaEff() const
Effective thermal turbulent diffusivity of mixture [kg/m/s].
tmp< volScalarField > nearInterface() const
Indicator of the proximity of the interface.
virtual word thermoName() const
Return the name of the thermo physics.
tmp< surfaceScalarField > surfaceTensionForce() const
virtual tmp< volScalarField > Cv() const
Heat capacity at constant volume [J/kg/K].
virtual tmp< volScalarField > rho() const
Density [kg/m^3].
virtual tmp< volScalarField > W() const
Molecular weight [kg/kmol].
virtual tmp< volScalarField > gamma() const
Gamma = Cp/Cv [].
virtual bool incompressible() const
Return true if the equation of state is incompressible.
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 > kappa() const
Thermal diffusivity for temperature of mixture [J/m/s/K].
tmp< volScalarField > rCv() const
Return the phase-averaged reciprocal Cv.
virtual tmp< volScalarField > Cp() const
Heat capacity at constant pressure [J/kg/K].
virtual void correct()
Update properties.
virtual volScalarField & he()
Enthalpy/Internal energy [J/kg].
virtual tmp< volScalarField > hc() const
Chemical enthalpy [J/kg].
virtual tmp< volScalarField > nu() const
Kinematic viscosity of mixture [m^2/s].
void correctRho(const volScalarField &dp)
Update densities for given pressure change.
virtual tmp< volScalarField > alphahe() const
Thermal diffusivity for energy of mixture [kg/m/s].
virtual tmp< volScalarField > CpByCpv() const
Heat capacity ratio [].
void solve()
Solve for the mixture phase-fractions.
virtual bool isochoric() const
Return true if the equation of state is isochoric.
const Time & time() const noexcept
Return time registry.
volScalarField mu_
Dynamic viscosity [kg/m/s].
Definition: psiThermo.H:68
volScalarField psi_
Compressibility [s^2/m^2].
Definition: psiThermo.H:65
Basic thermodynamics type based on the use of fitting functions for cp, h, s obtained from the templa...
A class for managing temporary objects.
Definition: tmp.H:65
A class for handling words, derived from Foam::string.
Definition: word.H:68
#define defineTypeNameAndDebug(Type, DebugSwitch)
Define the typeName and debug information.
Definition: className.H:121
U
Definition: pEqn.H:72
volScalarField & p
thermo correct()
const volScalarField & T
const volScalarField & mu
const tmp< volScalarField > & tCv
Definition: EEqn.H:5
const volScalarField & Cv
Definition: EEqn.H:8
const scalar gamma
Definition: EEqn.H:9
const tmp< volScalarField > & tCp
Definition: EEqn.H:4
const volScalarField & Cp
Definition: EEqn.H:7
faceListList boundary
dynamicFvMesh & mesh
engineTime & runTime
surfaceScalarField phir(fvc::flux(UdmModel.Udm()))
word alpharScheme("div(phirb,alpha)")
#define NotImplemented
Issue a FatalErrorIn for a function not currently implemented.
Definition: error.H:517
#define FatalErrorInFunction
Report an error message using Foam::FatalError.
Definition: error.H:453
Calculate the divergence of the given field.
Calculate the face-flux of the given field.
Calculate the gradient of the given field.
Calculate the mesh motion flux and convert fluxes from absolute to relative and back.
Calculate the snGrad of the given volField.
const cellShapeList & cells
volScalarField alphaEff("alphaEff", turbulence->nu()/Pr+alphat)
zeroField divU
Definition: alphaSuSp.H:3
zeroField Su
Definition: alphaSuSp.H:1
zeroField Sp
Definition: alphaSuSp.H:2
word timeName
Definition: getTimeIndex.H:3
label phasei
Definition: pEqn.H:27
kappaEff
Definition: TEqn.H:10
void explicitSolve(const RdeltaTType &rDeltaT, const RhoType &rho, volScalarField &psi, const surfaceScalarField &phiPsi, const SpType &Sp, const SuType &Su)
void limit(const RdeltaTType &rDeltaT, const RhoType &rho, const volScalarField &psi, const surfaceScalarField &phi, surfaceScalarField &phiPsi, const SpType &Sp, const SuType &Su, const PsiMaxType &psiMax, const PsiMinType &psiMin, const bool returnCorr)
void limitSum(UPtrList< scalarField > &phiPsiCorrs)
Definition: MULES.C:34
const dimensionedScalar kappa
Coulomb constant: default SI units: [N.m2/C2].
const dimensionedScalar sigma
Stefan-Boltzmann constant: default SI units: [W/m2/K4].
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.
tmp< GeometricField< typename outerProduct< vector, Type >::type, fvPatchField, volMesh > > grad(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcGrad.C:54
tmp< GeometricField< Type, fvsPatchField, surfaceMesh > > snGrad(const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
Definition: fvcSnGrad.C:47
tmp< GeometricField< Type, fvPatchField, volMesh > > div(const GeometricField< Type, fvsPatchField, surfaceMesh > &ssf)
Definition: fvcDiv.C:49
tmp< surfaceScalarField > flux(const volVectorField &vvf)
Return the face-flux field obtained from the given volVectorField.
Definition: fvcFlux.C:34
tmp< surfaceScalarField > absolute(const tmp< surfaceScalarField > &tphi, const volVectorField &U)
Return the given relative flux in absolute form.
Definition: fvcMeshPhi.C:190
vtkSmartPointer< vtkFloatArray > zeroField(const word &name, const label size)
Create named field initialized to zero.
Namespace for OpenFOAM.
GeometricField< vector, fvsPatchField, surfaceMesh > surfaceVectorField
dimensionedScalar det(const dimensionedSphericalTensor &dt)
dimensioned< scalar > dimensionedScalar
Dimensioned scalar obtained from generic dimensioned type.
label max(const labelHashSet &set, label maxValue=labelMin)
Find the max value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:47
GeometricField< vector, fvPatchField, volMesh > volVectorField
Definition: volFieldsFwd.H:83
const dimensionSet dimless
Dimensionless.
dimensionedScalar pos0(const dimensionedScalar &ds)
List< label > labelList
A List of labels.
Definition: List.H:66
const dimensionSet dimTime(0, 0, 1, 0, 0, 0, 0)
Definition: dimensionSets.H:53
GeometricField< scalar, fvPatchField, volMesh > volScalarField
Definition: volFieldsFwd.H:82
dimensionedScalar tanh(const dimensionedScalar &ds)
messageStream Info
Information stream (stdout output on master, null elsewhere)
dimensioned< Type > average(const DimensionedField< Type, GeoMesh > &df)
GeometricField< scalar, fvsPatchField, surfaceMesh > surfaceScalarField
Field< scalar > scalarField
Specialisation of Field<T> for scalar.
Ostream & endl(Ostream &os)
Add newline and flush stream.
Definition: Ostream.H:372
constexpr scalar degToRad() noexcept
Multiplication factor for degrees to radians conversion.
dimensioned< typename typeOfMag< Type >::type > mag(const dimensioned< Type > &dt)
label min(const labelHashSet &set, label minValue=labelMax)
Find the min value in labelHashSet, optionally limited by second argument.
Definition: hashSets.C:33
Field< vector > vectorField
Specialisation of Field<T> for vector.
static constexpr const zero Zero
Global zero (0)
Definition: zero.H:131
error FatalError
dimensionedScalar cbrt(const dimensionedScalar &ds)
word name(const expressions::valueTypeCode typeCode)
A word representation of a valueTypeCode. Empty for INVALID.
Definition: exprTraits.C:59
errorManipArg< error, int > exit(error &err, const int errNo=1)
Definition: errorManip.H:130
dimensionedScalar cos(const dimensionedScalar &ds)
dimensionedScalar acos(const dimensionedScalar &ds)
const dimensionSet dimMass(1, 0, 0, 0, 0, 0, 0)
Definition: dimensionSets.H:51
constexpr auto end(C &c) -> decltype(c.end())
Return iterator to the end of the container c.
Definition: stdFoam.H:158
volScalarField & alpha
tmp< volScalarField > trho
volScalarField & h
volScalarField & b
Definition: createFields.H:27
scalar T0
Definition: createFields.H:22
volScalarField & e
Definition: createFields.H:11
const dictionary & alphaControls
Definition: alphaControls.H:1
label nAlphaSubCycles(alphaControls.get< label >("nAlphaSubCycles"))
#define forAll(list, i)
Loop across all elements in list.
Definition: stdFoam.H:333
#define forAllConstIters(container, iter)
Iterate across all elements of the container object with const access.
Definition: stdFoam.H:278
surfaceScalarField alphaPhi(phi.name()+alpha1.name(), fvc::flux(phi, alpha1, alphaScheme))
Unit conversion functions.