multiphase: Changed phiU to phiHbyA and cache HbyA rather than store the velocity predictor in U

2025-11-28 03:28:01 +00:00 · 2012-04-19 09:39:56 +01:00
parent a6fac2d8bd
commit 59f63768d6
14 changed files with 207 additions and 136 deletions
--- a/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/compressibleInterFoam/pEqn.H
@ -2,8 +2,8 @@
    rho1 = eos1->rho(p, T);
    rho2 = eos2->rho(p, T);

-    volScalarField rAU = 1.0/UEqn.A();
-    surfaceScalarField rAUf = fvc::interpolate(rAU);
+    volScalarField rAU("rAU", 1.0/UEqn.A());
+    surfaceScalarField rAUf("Dp", fvc::interpolate(rAU));

    tmp<fvScalarMatrix> p_rghEqnComp1;
    tmp<fvScalarMatrix> p_rghEqnComp2;
@ -27,20 +27,25 @@
          - fvc::Sp(fvc::div(phid2), p_rgh);
    }

-    U = rAU*UEqn.H();
+    volVectorField HbyA("HbyA", U);
+    HbyA = rAU*UEqn.H();

-    surfaceScalarField phiU
+    surfaceScalarField phiHbyA
    (
-        "phiU",
-        (fvc::interpolate(U) & mesh.Sf())
+        "phiHbyA",
+        (fvc::interpolate(HbyA) & mesh.Sf())
      + fvc::ddtPhiCorr(rAU, rho, U, phi)
    );

-    phi = phiU +
+    surfaceScalarField phig
+    (
        (
            fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
          - ghf*fvc::snGrad(rho)
-        )*rAUf*mesh.magSf();
+        )*rAUf*mesh.magSf()
+    );
+
+    phiHbyA += phig;

    // Thermodynamic density needs to be updated by psi*d(p) after the
    // pressure solution - done in 2 parts. Part 1:
@ -50,7 +55,7 @@
    {
        fvScalarMatrix p_rghEqnIncomp
        (
-            fvc::div(phi)
+            fvc::div(phiHbyA)
          - fvm::laplacian(rAUf, p_rgh)
        );

@ -74,13 +79,15 @@
                pos(alpha2)*(p_rghEqnComp2 & p_rgh)/rho2
              - pos(alpha1)*(p_rghEqnComp1 & p_rgh)/rho1
            );
-            phi += p_rghEqnIncomp.flux();
+
+            phi = phiHbyA + p_rghEqnIncomp.flux();
+
+            U = HbyA
+              + rAU*fvc::reconstruct((phig + p_rghEqnIncomp.flux())/rAUf);
+            U.correctBoundaryConditions();
        }
    }

-    U += rAU*fvc::reconstruct((phi - phiU)/rAUf);
-    U.correctBoundaryConditions();
-
    p = max(p_rgh + (alpha1*rho1 + alpha2*rho2)*gh, pMin);

    rho1 = eos1->rho(p, T);
--- a/applications/solvers/multiphase/interFoam/MRFInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/interFoam/MRFInterFoam/pEqn.H
@ -1,6 +1,6 @@
 {
-    volScalarField rAU("Dp", 1.0/UEqn.A());
-    surfaceScalarField rAUf("Dpf", fvc::interpolate(rAU));
+    volScalarField rAU("rAU", 1.0/UEqn.A());
+    surfaceScalarField rAUf("Dp", fvc::interpolate(rAU));

    volVectorField HbyA("HbyA", U);
    HbyA = rAU*UEqn.H();
@ -11,7 +11,6 @@
        (fvc::interpolate(HbyA) & mesh.Sf())
      + fvc::ddtPhiCorr(rAU, rho, U, phi)
    );
-
    mrfZones.relativeFlux(phiHbyA);

    adjustPhi(phiHbyA, U, p_rgh);
--- a/applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H
+++ b/applications/solvers/multiphase/interFoam/interDyMFoam/pEqn.H
@ -1,6 +1,6 @@
 {
-    volScalarField rAU("Dp", 1.0/UEqn.A());
-    surfaceScalarField rAUf("Dpf", fvc::interpolate(rAU));
+    volScalarField rAU("rAU", 1.0/UEqn.A());
+    surfaceScalarField rAUf("Dp", fvc::interpolate(rAU));

    volVectorField HbyA("HbyA", U);
    HbyA = rAU*UEqn.H();
--- a/applications/solvers/multiphase/interFoam/pEqn.H
+++ b/applications/solvers/multiphase/interFoam/pEqn.H
@ -1,6 +1,6 @@
 {
-    volScalarField rAU("Dp", 1.0/UEqn.A());
-    surfaceScalarField rAUf("Dpf", fvc::interpolate(rAU));
+    volScalarField rAU("rAU", 1.0/UEqn.A());
+    surfaceScalarField rAUf("Dp", fvc::interpolate(rAU));

    volVectorField HbyA("HbyA", U);
    HbyA = rAU*UEqn.H();
--- a/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H
+++ b/applications/solvers/multiphase/interPhaseChangeFoam/pEqn.H
@ -1,24 +1,28 @@
 {
-    volScalarField rAU(1.0/UEqn.A());
-    surfaceScalarField rAUf(fvc::interpolate(rAU));
+    volScalarField rAU("rAU", 1.0/UEqn.A());
+    surfaceScalarField rAUf("Dp", fvc::interpolate(rAU));

-    U = rAU*UEqn.H();
+    volVectorField HbyA("HbyA", U);
+    HbyA = rAU*UEqn.H();

-    surfaceScalarField phiU
+    surfaceScalarField phiHbyA
    (
-        "phiU",
-        (fvc::interpolate(U) & mesh.Sf())
+        "phiHbyA",
+        (fvc::interpolate(HbyA) & mesh.Sf())
      + fvc::ddtPhiCorr(rAU, rho, U, phi)
    );

-    adjustPhi(phiU, U, p_rgh);
+    adjustPhi(phiHbyA, U, p_rgh);

-    phi =
-        phiU
-      + (
+    surfaceScalarField phig
+    (
+        (
            fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
          - ghf*fvc::snGrad(rho)
-        )*rAUf*mesh.magSf();
+        )*rAUf*mesh.magSf()
+    );
+
+    phiHbyA += phig;

    Pair<tmp<volScalarField> > vDotP = twoPhaseProperties->vDotP();
    const volScalarField& vDotcP = vDotP[0]();
@ -28,7 +32,7 @@
    {
        fvScalarMatrix p_rghEqn
        (
-            fvc::div(phi) - fvm::laplacian(rAUf, p_rgh)
+            fvc::div(phiHbyA) - fvm::laplacian(rAUf, p_rgh)
          - (vDotvP - vDotcP)*(pSat - rho*gh) + fvm::Sp(vDotvP - vDotcP, p_rgh)
        );

@ -38,13 +42,13 @@

        if (pimple.finalNonOrthogonalIter())
        {
-            phi += p_rghEqn.flux();
+            phi = phiHbyA + p_rghEqn.flux();
+
+            U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rAUf);
+            U.correctBoundaryConditions();
        }
    }

-    U += rAU*fvc::reconstruct((phi - phiU)/rAUf);
-    U.correctBoundaryConditions();
-
    #include "continuityErrs.H"

    p == p_rgh + rho*gh;
--- a/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/multiphaseInterFoam/MRFMultiphaseInterFoam/pEqn.H
@ -1,31 +1,32 @@
 {
-    volScalarField rAU(1.0/UEqn.A());
-    surfaceScalarField rAUf(fvc::interpolate(rAU));
+    volScalarField rAU("rAU", 1.0/UEqn.A());
+    surfaceScalarField rAUf("Dp", fvc::interpolate(rAU));

-    U = rAU*UEqn.H();
+    volVectorField HbyA("HbyA", U);
+    HbyA = rAU*UEqn.H();

-    surfaceScalarField phiU
+    surfaceScalarField phiHbyA
    (
-        "phiU",
-        (fvc::interpolate(U) & mesh.Sf())
-      //+ fvc::ddtPhiCorr(rAU, rho, U, phi)
+        "phiHbyA",
+        (fvc::interpolate(HbyA) & mesh.Sf())
+      + fvc::ddtPhiCorr(rAU, rho, U, phi)
    );
-    mrfZones.relativeFlux(phiU);
+    mrfZones.relativeFlux(phiHbyA);

-    adjustPhi(phiU, U, p_rgh);
+    adjustPhi(phiHbyA, U, p_rgh);

-    phi =
-        phiU
-      + (
-            mixture.surfaceTensionForce()
-          - ghf*fvc::snGrad(rho)
-        )*rAUf*mesh.magSf();
+    surfaceScalarField phig
+    (
+        - ghf*fvc::snGrad(rho)*rAUf*mesh.magSf()
+    );
+
+    phiHbyA += phig;

    while (pimple.correctNonOrthogonal())
    {
        fvScalarMatrix p_rghEqn
        (
-            fvm::laplacian(rAUf, p_rgh) == fvc::div(phi)
+            fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
        );

        p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@ -34,13 +35,13 @@

        if (pimple.finalNonOrthogonalIter())
        {
-            phi -= p_rghEqn.flux();
+            phi = phiHbyA - p_rghEqn.flux();
+
+            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
+            U.correctBoundaryConditions();
        }
    }

-    U += rAU*fvc::reconstruct((phi - phiU)/rAUf);
-    U.correctBoundaryConditions();
-
    #include "continuityErrs.H"

    p == p_rgh + rho*gh;
--- a/applications/solvers/multiphase/multiphaseInterFoam/pEqn.H
+++ b/applications/solvers/multiphase/multiphaseInterFoam/pEqn.H
@ -1,29 +1,34 @@
 {
-    volScalarField rAU(1.0/UEqn.A());
-    surfaceScalarField rAUf(fvc::interpolate(rAU));
+    volScalarField rAU("rAU", 1.0/UEqn.A());
+    surfaceScalarField rAUf("Dp", fvc::interpolate(rAU));

-    U = rAU*UEqn.H();
+    volVectorField HbyA("HbyA", U);
+    HbyA = rAU*UEqn.H();

-    surfaceScalarField phiU
+    surfaceScalarField phiHbyA
    (
-        "phiU",
-        (fvc::interpolate(U) & mesh.Sf())
+        "phiHbyA",
+        (fvc::interpolate(HbyA) & mesh.Sf())
      + fvc::ddtPhiCorr(rAU, rho, U, phi)
    );

-    adjustPhi(phiU, U, p_rgh);
+    adjustPhi(phiHbyA, U, p_rgh);

-    phi = phiU +
+    surfaceScalarField phig
    (
-        mixture.surfaceTensionForce()
-      - ghf*fvc::snGrad(rho)
-    )*rAUf*mesh.magSf();
+        (
+            mixture.surfaceTensionForce()
+          - ghf*fvc::snGrad(rho)
+        )*rAUf*mesh.magSf()
+    );
+
+    phiHbyA += phig;

    while (pimple.correctNonOrthogonal())
    {
        fvScalarMatrix p_rghEqn
        (
-            fvm::laplacian(rAUf, p_rgh) == fvc::div(phi)
+            fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
        );

        p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@ -32,13 +37,13 @@

        if (pimple.finalNonOrthogonalIter())
        {
-            phi -= p_rghEqn.flux();
+            phi = phiHbyA - p_rghEqn.flux();
+
+            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
+            U.correctBoundaryConditions();
        }
    }

-    U += rAU*fvc::reconstruct((phi - phiU)/rAUf);
-    U.correctBoundaryConditions();
-
    #include "continuityErrs.H"

    p == p_rgh + rho*gh;
--- a/applications/solvers/multiphase/settlingFoam/pEqn.H
+++ b/applications/solvers/multiphase/settlingFoam/pEqn.H
@ -1,54 +1,60 @@
-volScalarField rAU(1.0/UEqn.A());
-
-surfaceScalarField rAUf
-(
-    "(rho*(1|A(U)))",
-    fvc::interpolate(rho)*fvc::interpolate(rAU)
-);
-
-U = rAU*UEqn.H();
-phi =
-    fvc::interpolate(rho)
-   *(
-       (fvc::interpolate(U) & mesh.Sf())
-     + fvc::ddtPhiCorr(rAU, rho, U, phi)
-    );
-
-surfaceScalarField phiU("phiU", phi);
-phi -= ghf*fvc::snGrad(rho)*rAUf*mesh.magSf();
-
-while (pimple.correctNonOrthogonal())
 {
-    fvScalarMatrix p_rghEqn
+    volScalarField rAU("rAU", 1.0/UEqn.A());
+
+    surfaceScalarField rAUf("Dp", fvc::interpolate(rho*rAU));
+
+    volVectorField HbyA("HbyA", U);
+    HbyA = rAU*UEqn.H();
+
+    surfaceScalarField phiHbyA
    (
-        fvm::laplacian(rAUf, p_rgh) == fvc::ddt(rho) + fvc::div(phi)
+        fvc::interpolate(rho)
+       *(
+           (fvc::interpolate(HbyA) & mesh.Sf())
+         + fvc::ddtPhiCorr(rAU, rho, U, phi)
+        )
    );

-    p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
+    surfaceScalarField phig
+    (
+        - ghf*fvc::snGrad(rho)*rAUf*mesh.magSf()
+    );

-    p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
+    phiHbyA += phig;

-    if (pimple.finalNonOrthogonalIter())
+    while (pimple.correctNonOrthogonal())
    {
-        phi -= p_rghEqn.flux();
+        fvScalarMatrix p_rghEqn
+        (
+            fvm::laplacian(rAUf, p_rgh) == fvc::ddt(rho) + fvc::div(phiHbyA)
+        );
+
+        p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
+
+        p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));
+
+        if (pimple.finalNonOrthogonalIter())
+        {
+            phi = phiHbyA - p_rghEqn.flux();
+
+            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
+            U.correctBoundaryConditions();
+        }
    }
+
+    p == p_rgh + rho*gh;
+
+    if (p_rgh.needReference())
+    {
+        p += dimensionedScalar
+        (
+            "p",
+            p.dimensions(),
+            pRefValue - getRefCellValue(p, pRefCell)
+        );
+        p_rgh = p - rho*gh;
+    }
+
+    #include "rhoEqn.H"
+    #include "compressibleContinuityErrs.H"
 }
-
-p == p_rgh + rho*gh;
-
-if (p_rgh.needReference())
-{
-    p += dimensionedScalar
-    (
-        "p",
-        p.dimensions(),
-        pRefValue - getRefCellValue(p, pRefCell)
-    );
-    p_rgh = p - rho*gh;
-}
-
-#include "rhoEqn.H"
-#include "compressibleContinuityErrs.H"
-
-U += rAU*fvc::reconstruct((phi - phiU)/rAUf);
-U.correctBoundaryConditions();
--- a/applications/solvers/multiphase/twoLiquidMixingFoam/pEqn.H
+++ b/applications/solvers/multiphase/twoLiquidMixingFoam/pEqn.H
@ -1,24 +1,31 @@
 {
-    volScalarField rAU(1.0/UEqn.A());
-    surfaceScalarField rAUf(fvc::interpolate(rAU));
+    volScalarField rAU("rAU", 1.0/UEqn.A());
+    surfaceScalarField rAUf("Dp", fvc::interpolate(rAU));

-    U = rAU*UEqn.H();
-    surfaceScalarField phiU
+    volVectorField HbyA("HbyA", U);
+    HbyA = rAU*UEqn.H();
+
+    surfaceScalarField phiHbyA
    (
-        "phiU",
-        (fvc::interpolate(U) & mesh.Sf())
+        "phiHbyA",
+        (fvc::interpolate(HbyA) & mesh.Sf())
      + fvc::ddtPhiCorr(rAU, rho, U, phi)
    );

-    adjustPhi(phiU, U, p_rgh);
+    adjustPhi(phiHbyA, U, p_rgh);

-    phi = phiU - ghf*fvc::snGrad(rho)*rAUf*mesh.magSf();
+    surfaceScalarField phig
+    (
+        - ghf*fvc::snGrad(rho)*rAUf*mesh.magSf()
+    );
+
+    phiHbyA += phig;

    while (pimple.correctNonOrthogonal())
    {
        fvScalarMatrix p_rghEqn
        (
-            fvm::laplacian(rAUf, p_rgh) == fvc::div(phi)
+            fvm::laplacian(rAUf, p_rgh) == fvc::div(phiHbyA)
        );

        p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));
@ -27,13 +34,13 @@

        if (pimple.finalNonOrthogonalIter())
        {
-            phi -= p_rghEqn.flux();
+            phi = phiHbyA - p_rghEqn.flux();
+
+            U = HbyA + rAU*fvc::reconstruct((phig - p_rghEqn.flux())/rAUf);
+            U.correctBoundaryConditions();
        }
    }

-    U += rAU*fvc::reconstruct((phi - phiU)/rAUf);
-    U.correctBoundaryConditions();
-
    #include "continuityErrs.H"

    p == p_rgh + rho*gh;