pEqn.H 1.52 KB
Newer Older
1 2
volScalarField rAU(1.0/UEqn.A());
volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
3 4 5 6 7 8
surfaceScalarField phiHbyA("phiHbyA", fvc::flux(HbyA));

if (pimple.ddtCorr())
{
    phiHbyA += MRF.zeroFilter(fvc::interpolate(rAU)*fvc::ddtCorr(U, phi, Uf));
}
9 10 11 12
else
{
    phiHbyA += MRF.zeroFilter(fvc::interpolate(rAU));
}
henry's avatar
henry committed
13

Henry's avatar
Henry committed
14
MRF.makeRelative(phiHbyA);
15

16 17 18 19 20 21
if (p.needReference())
{
    fvc::makeRelative(phiHbyA, U);
    adjustPhi(phiHbyA, U, p);
    fvc::makeAbsolute(phiHbyA, U);
}
22

23 24 25 26
tmp<volScalarField> rAtU(rAU);

if (pimple.consistent())
{
27
    rAtU = 1.0/max(1.0/rAU - UEqn.H1(), 0.1/rAU);
28 29 30 31 32 33 34
    phiHbyA +=
        fvc::interpolate(rAtU() - rAU)*fvc::snGrad(p)*mesh.magSf();
    HbyA -= (rAU - rAtU())*fvc::grad(p);
}

if (pimple.nCorrPISO() <= 1)
{
35
    tUEqn.clear();
36 37
}

38 39
// Update the pressure BCs to ensure flux consistency
constrainPressure(p, U, phiHbyA, rAtU(), MRF);
40

henry's avatar
henry committed
41
// Non-orthogonal pressure corrector loop
42
while (pimple.correctNonOrthogonal())
henry's avatar
henry committed
43 44 45
{
    fvScalarMatrix pEqn
    (
46
        fvm::laplacian(rAtU(), p) == fvc::div(phiHbyA)
henry's avatar
henry committed
47 48 49 50
    );

    pEqn.setReference(pRefCell, pRefValue);

51
    pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
henry's avatar
henry committed
52

53
    if (pimple.finalNonOrthogonalIter())
henry's avatar
henry committed
54
    {
Henry's avatar
Henry committed
55
        phi = phiHbyA - pEqn.flux();
henry's avatar
henry committed
56 57 58 59 60
    }
}

#include "continuityErrs.H"

61 62
// Explicitly relax pressure for momentum corrector
p.relax();
henry's avatar
henry committed
63

64
U = HbyA - rAtU*fvc::grad(p);
henry's avatar
henry committed
65
U.correctBoundaryConditions();
66
fvOptions.correct(U);
67 68 69 70 71 72

// Correct Uf if the mesh is moving
fvc::correctUf(Uf, U, phi);

// Make the fluxes relative to the mesh motion
fvc::makeRelative(phi, U);