1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
|
{
word alphaScheme("div(phi,alpha)");
word alpharScheme("div(phirb,alpha)");
surfaceScalarField phir("phir", phic*interface.nHatf());
Pair<tmp<volScalarField>> vDotAlphal =
mixture->vDotAlphal();
const volScalarField& vDotcAlphal = vDotAlphal[0]();
const volScalarField& vDotvAlphal = vDotAlphal[1]();
const volScalarField vDotvmcAlphal(vDotvAlphal - vDotcAlphal);
tmp<surfaceScalarField> talphaPhi;
if (MULESCorr)
{
fvScalarMatrix alpha1Eqn
(
fv::EulerDdtScheme<scalar>(mesh).fvmDdt(alpha1)
+ fv::gaussConvectionScheme<scalar>
(
mesh,
phi,
upwind<scalar>(mesh, phi)
).fvmDiv(phi, alpha1)
- fvm::Sp(divU, alpha1)
==
fvm::Sp(vDotvmcAlphal, alpha1)
+ vDotcAlphal
);
alpha1Eqn.solve();
Info<< "Phase-1 volume fraction = "
<< alpha1.weightedAverage(mesh.Vsc()).value()
<< " Min(" << alpha1.name() << ") = " << min(alpha1).value()
<< " Max(" << alpha1.name() << ") = " << max(alpha1).value()
<< endl;
talphaPhi = alpha1Eqn.flux();
}
volScalarField alpha10("alpha10", alpha1);
for (int aCorr=0; aCorr<nAlphaCorr; aCorr++)
{
tmp<surfaceScalarField> talphaPhiCorr
(
fvc::flux
(
phi,
alpha1,
alphaScheme
)
+ fvc::flux
(
-fvc::flux(-phir, alpha2, alpharScheme),
alpha1,
alpharScheme
)
);
if (MULESCorr)
{
talphaPhiCorr.ref() -= talphaPhi();
volScalarField alpha100("alpha100", alpha10);
alpha10 = alpha1;
MULES::correct
(
geometricOneField(),
alpha1,
talphaPhi(),
talphaPhiCorr.ref(),
vDotvmcAlphal,
(
divU*(alpha10 - alpha100)
- vDotvmcAlphal*alpha10
)(),
1,
0
);
// Under-relax the correction for all but the 1st corrector
if (aCorr == 0)
{
talphaPhi.ref() += talphaPhiCorr();
}
else
{
alpha1 = 0.5*alpha1 + 0.5*alpha10;
talphaPhi.ref() += 0.5*talphaPhiCorr();
}
}
else
{
MULES::explicitSolve
(
geometricOneField(),
alpha1,
phi,
talphaPhiCorr.ref(),
vDotvmcAlphal,
(divU*alpha1 + vDotcAlphal)(),
1,
0
);
talphaPhi = talphaPhiCorr;
}
alpha2 = 1.0 - alpha1;
}
rhoPhi = talphaPhi()*(rho1 - rho2) + phi*rho2;
Info<< "Liquid phase volume fraction = "
<< alpha1.weightedAverage(mesh.V()).value()
<< " Min(" << alpha1.name() << ") = " << min(alpha1).value()
<< " Max(" << alpha1.name() << ") = " << max(alpha1).value()
<< endl;
}
|
