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
|
{
surfaceScalarField alphaPhi
(
IOobject
(
"alphaPhi",
runTime.timeName(),
mesh
),
mesh,
dimensionedScalar("0", phi.dimensions(), 0)
);
surfaceScalarField phir(fvc::flux(UdmModel.Udm()));
if (nAlphaSubCycles > 1)
{
dimensionedScalar totalDeltaT = runTime.deltaT();
surfaceScalarField alphaPhiSum
(
IOobject
(
"alphaPhiSum",
runTime.timeName(),
mesh
),
mesh,
dimensionedScalar("0", phi.dimensions(), 0)
);
for
(
subCycle<volScalarField> alphaSubCycle(alpha1, nAlphaSubCycles);
!(++alphaSubCycle).end();
)
{
#include "alphaEqn.H"
alphaPhiSum += (runTime.deltaT()/totalDeltaT)*alphaPhi;
}
alphaPhi = alphaPhiSum;
}
else
{
#include "alphaEqn.H"
}
// Apply the diffusion term separately to allow implicit solution
// and boundedness of the explicit advection
{
fvScalarMatrix alpha1Eqn
(
fvm::ddt(alpha1) - fvc::ddt(alpha1)
- fvm::laplacian(turbulence->nut(), alpha1)
);
alpha1Eqn.solve(mesh.solver("alpha1Diffusion"));
alphaPhi += alpha1Eqn.flux();
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;
}
rhoPhi = alphaPhi*(rho1 - rho2) + phi*rho2;
rho = mixture.rho();
}
|
