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
|
{
word alphaScheme("div(phi,alpha)");
word alpharScheme("div(phirb,alpha)");
if (MULESCorr)
{
fvScalarMatrix alpha1Eqn
(
fv::EulerDdtScheme<scalar>(mesh).fvmDdt(alpha1)
+ fv::gaussConvectionScheme<scalar>
(
mesh,
phi,
upwind<scalar>(mesh, phi)
).fvmDiv(phi, alpha1)
);
solve(alpha1Eqn);
Info<< "Phase-1 volume fraction = "
<< alpha1.weightedAverage(mesh.Vsc()).value()
<< " Min(" << alpha1.name() << ") = " << min(alpha1).value()
<< " Max(" << alpha1.name() << ") = " << max(alpha1).value()
<< endl;
tmp<surfaceScalarField> talphaPhiUD(alpha1Eqn.flux());
alphaPhi = talphaPhiUD();
if (alphaApplyPrevCorr && talphaPhiCorr0.valid())
{
Info<< "Applying the previous iteration correction flux" << endl;
MULES::correct
(
alpha1,
alphaPhi,
talphaPhiCorr0.ref(),
mixture.alphaMax(),
0
);
alphaPhi += talphaPhiCorr0();
}
// Cache the upwind-flux
talphaPhiCorr0 = talphaPhiUD;
}
for (int aCorr=0; aCorr<nAlphaCorr; aCorr++)
{
tmp<surfaceScalarField> talphaPhiUn
(
fvc::flux
(
phi,
alpha1,
alphaScheme
)
+ fvc::flux
(
phir,
alpha1,
alpharScheme
)
);
if (MULESCorr)
{
tmp<surfaceScalarField> talphaPhiCorr(talphaPhiUn() - alphaPhi);
volScalarField alpha10("alpha10", alpha1);
MULES::correct
(
alpha1,
talphaPhiUn(),
talphaPhiCorr.ref(),
mixture.alphaMax(),
0
);
// Under-relax the correction for all but the 1st corrector
if (aCorr == 0)
{
alphaPhi += talphaPhiCorr();
}
else
{
alpha1 = 0.5*alpha1 + 0.5*alpha10;
alphaPhi += 0.5*talphaPhiCorr();
}
}
else
{
alphaPhi = talphaPhiUn;
MULES::explicitSolve
(
alpha1,
phi,
alphaPhi,
mixture.alphaMax(),
0
);
}
}
if (alphaApplyPrevCorr && MULESCorr)
{
talphaPhiCorr0 = alphaPhi - talphaPhiCorr0;
}
alpha2 = 1.0 - alpha1;
Info<< "Phase-1 volume fraction = "
<< alpha1.weightedAverage(mesh.Vsc()).value()
<< " Min(" << alpha1.name() << ") = " << min(alpha1).value()
<< " Max(" << alpha1.name() << ") = " << max(alpha1).value()
<< endl;
}
|
