Improved SIMPLEC consistency in two-phase flow solver. Only effects flows with gravity and/or surface tension (PIMPLE is still likely better in such cases)

Author: fppimenta

fe40/src/solvers/rheoInterFoam/pEqn.H

@@ -36,7 +36,7 @@
       rAtU().rename("rAtU");
 
       phi += fvc::interpolate(rAtU() - rAU)*fvc::snGrad(pd)*mesh.magSf();
-      U -= (rAU - rAtU())*fvc::grad(pd);
+      U -= (rAU - rAtU()) * fvc::reconstruct( fvc::snGrad(pd) * mesh.magSf() );
     }
 
 
@@ -66,7 +66,8 @@
 
     if (simplec)
      {
-          U += rAU*fvc::reconstruct(phig/rAUf) - rAtU()*fvc::grad(pd);
+          surfaceScalarField rAtUf("rAtUf", fvc::interpolate(rAtU()));
+          U += rAU*fvc::reconstruct(phig/rAUf) - rAtU()*fvc::reconstruct(pdEqn.flux()/rAtUf); 
      }
     else
      {

of60/src/solvers/rheoInterFoam/pEqn.H

@@ -41,7 +41,7 @@
     {
       rAtU = 1.0/(1.0/rAU() - UEqn.H1());
       phiHbyA += fvc::interpolate(rAtU() - rAU())*fvc::snGrad(p_rgh)*mesh.magSf();
-      HbyA -= (rAU() - rAtU())*fvc::grad(p_rgh);
+      HbyA -= (rAU() - rAtU()) * fvc::reconstruct( fvc::snGrad(p_rgh) * mesh.magSf() );
 
       // Update the pressure BCs to ensure flux consistency
       constrainPressure(p_rgh, U, phiHbyA, rAtU(), MRF);
@@ -78,7 +78,8 @@
 
             if (simplec)
             {
-               U = HbyA + rAU()*fvc::reconstruct(phig/rAUf) - rAtU()*fvc::grad(p_rgh);
+               surfaceScalarField rAtUf("rAtUf", fvc::interpolate(rAtU()));
+               U = HbyA + rAU()*fvc::reconstruct(phig/rAUf) - rAtU()*fvc::reconstruct(p_rghEqn.flux()/rAtUf);
             }
             else
             {