alphaEqn.H 3.19 KB
Newer Older
1
2
3
4
5
6
{
    word alphaScheme("div(phi,alpha)");
    word alpharScheme("div(phirb,alpha)");

    surfaceScalarField phir("phir", phic*interface.nHatf());

7
    Pair<tmp<volScalarField>> vDotAlphal =
8
        mixture->vDotAlphal();
9
10
11
12
    const volScalarField& vDotcAlphal = vDotAlphal[0]();
    const volScalarField& vDotvAlphal = vDotAlphal[1]();
    const volScalarField vDotvmcAlphal(vDotvAlphal - vDotcAlphal);

13
    tmp<surfaceScalarField> talphaPhi;
14
15

    if (MULESCorr)
16
    {
17
18
        fvScalarMatrix alpha1Eqn
        (
19
20
21
22
23
24
25
26
            fv::EulerDdtScheme<scalar>(mesh).fvmDdt(alpha1)
          + fv::gaussConvectionScheme<scalar>
            (
                mesh,
                phi,
                upwind<scalar>(mesh, phi)
            ).fvmDiv(phi, alpha1)
          - fvm::Sp(divU, alpha1)
27
28
29
30
31
32
33
34
35
         ==
            fvm::Sp(vDotvmcAlphal, alpha1)
          + vDotcAlphal
        );

        alpha1Eqn.solve();

        Info<< "Phase-1 volume fraction = "
            << alpha1.weightedAverage(mesh.Vsc()).value()
36
37
            << "  Min(" << alpha1.name() << ") = " << min(alpha1).value()
            << "  Max(" << alpha1.name() << ") = " << max(alpha1).value()
38
39
            << endl;

40
        talphaPhi = alpha1Eqn.flux();
41
42
43
44
45
46
    }

    volScalarField alpha10("alpha10", alpha1);

    for (int aCorr=0; aCorr<nAlphaCorr; aCorr++)
    {
47
        tmp<surfaceScalarField> talphaPhiCorr
48
        (
49
50
51
52
53
54
55
56
            fvc::flux
            (
                phi,
                alpha1,
                alphaScheme
            )
          + fvc::flux
            (
57
                -fvc::flux(-phir, alpha2, alpharScheme),
58
59
                alpha1,
                alpharScheme
60
61
            )
        );
62

63
64
        if (MULESCorr)
        {
65
            talphaPhiCorr() -= talphaPhi();
66

67
68
69
70
            volScalarField alpha100("alpha100", alpha10);
            alpha10 = alpha1;

            MULES::correct
71
            (
72
73
                geometricOneField(),
                alpha1,
74
75
                talphaPhi(),
                talphaPhiCorr(),
76
77
78
79
80
81
82
83
                vDotvmcAlphal,
                (
                    divU*(alpha10 - alpha100)
                  - vDotvmcAlphal*alpha10
                )(),
                1,
                0
            );
84

85
86
87
            // Under-relax the correction for all but the 1st corrector
            if (aCorr == 0)
            {
88
                talphaPhi() += talphaPhiCorr();
89
90
91
92
            }
            else
            {
                alpha1 = 0.5*alpha1 + 0.5*alpha10;
93
                talphaPhi() += 0.5*talphaPhiCorr();
94
            }
95
96
97
98
        }
        else
        {
            MULES::explicitSolve
99
            (
100
101
102
                geometricOneField(),
                alpha1,
                phi,
103
                talphaPhiCorr(),
104
105
106
107
108
                vDotvmcAlphal,
                (divU*alpha1 + vDotcAlphal)(),
                1,
                0
            );
109

110
            talphaPhi = talphaPhiCorr;
111
        }
112

113
        alpha2 = 1.0 - alpha1;
114
115
    }

116
    rhoPhi = talphaPhi()*(rho1 - rho2) + phi*rho2;
117

118
119
    Info<< "Liquid phase volume fraction = "
        << alpha1.weightedAverage(mesh.V()).value()
120
121
        << "  Min(" << alpha1.name() << ") = " << min(alpha1).value()
        << "  Max(" << alpha1.name() << ") = " << max(alpha1).value()
122
123
        << endl;
}