Newer
Older
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
{
if (pimple.nCorrPIMPLE() == 1)
{
p =
(
rho
- (1.0 - gamma)*rhol0
- ((gamma*psiv + (1.0 - gamma)*psil) - psi)*pSat
)/psi;
}
surfaceScalarField rhof("rhof", fvc::interpolate(rho));
volScalarField rAU(1.0/UEqn.A());
surfaceScalarField rAUf("Dp", rhof*fvc::interpolate(rAU));
volVectorField HbyA("HbyA", U);
HbyA = rAU*UEqn.H();
phiv = (fvc::interpolate(HbyA) & mesh.Sf())
+ fvc::ddtPhiCorr(rAU, rho, U, phivAbs);
fvc::makeRelative(phiv, U);
surfaceScalarField phiGradp(rAUf*mesh.magSf()*fvc::snGrad(p));
phiv -= phiGradp/rhof;
volScalarField rho0(rho - psi*p);
while (pimple.correctNonOrthogonal())
{
fvScalarMatrix pEqn
(
fvc::ddt(rho)
+ psi*correction(fvm::ddt(p))
+ fvc::div(phiv, rho)
+ fvc::div(phiGradp)
- fvm::laplacian(rAUf, p)
);
pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
if (pimple.finalNonOrthogonalIter())
{
phiv += (phiGradp + pEqn.flux())/rhof;
phivAbs = fvc::absolute(phiv, U);
}
}
Info<< "Predicted p max-min : " << max(p).value()
<< " " << min(p).value() << endl;
rho == max(rho0 + psi*p, rhoMin);
#include "gammaPsi.H"
p =
(
rho
- (1.0 - gamma)*rhol0
- ((gamma*psiv + (1.0 - gamma)*psil) - psi)*pSat
)/psi;
p.correctBoundaryConditions();
Info<< "Phase-change corrected p max-min : " << max(p).value()
<< " " << min(p).value() << endl;
// Correct velocity
U = HbyA - rAU*fvc::grad(p);
// Remove the swirl component of velocity for "wedge" cases
if (pimple.dict().found("removeSwirl"))
{
label swirlCmpt(readLabel(pimple.dict().lookup("removeSwirl")));
Info<< "Removing swirl component-" << swirlCmpt << " of U" << endl;
U.field().replace(swirlCmpt, 0.0);
}
U.correctBoundaryConditions();
Info<< "max(U) " << max(mag(U)).value() << endl;
}