fixedFluxPressure BC: the snGrad is now pushed into the BC from pEqn.…

…H rather than being evaluated in the BC
OpenCFD · Sep 10, 2013 · f34c618 · f34c618
1 parent 1a41224
commit f34c618
Show file tree

Hide file tree

Showing 77 changed files with 526 additions and 773 deletions.
diff --git a/applications/solvers/DNS/dnsFoam/dnsFoam.C b/applications/solvers/DNS/dnsFoam/dnsFoam.C
@@ -86,20 +86,20 @@ int main(int argc, char *argv[])
         for (int corr=1; corr<=1; corr++)
         {
             volScalarField rAU(1.0/UEqn.A());
-            surfaceScalarField Dp("Dp", fvc::interpolate(rAU));
+            surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU));
             volVectorField HbyA("HbyA", U);
             HbyA = rAU*UEqn.H();
 
             surfaceScalarField phiHbyA
             (
                 "phiHbyA",
                 (fvc::interpolate(HbyA) & mesh.Sf())
-              + Dp*fvc::ddtCorr(U, phi)
+              + rAUf*fvc::ddtCorr(U, phi)
             );
 
             fvScalarMatrix pEqn
             (
-                fvm::laplacian(Dp, p) == fvc::div(phiHbyA)
+                fvm::laplacian(rAUf, p) == fvc::div(phiHbyA)
             );
 
             pEqn.solve();

diff --git a/applications/solvers/combustion/XiFoam/pEqn.H b/applications/solvers/combustion/XiFoam/pEqn.H
@@ -1,7 +1,7 @@
 rho = thermo.rho();
 
 volScalarField rAU(1.0/UEqn.A());
-surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
 
 volVectorField HbyA("HbyA", U);
 HbyA = rAU*UEqn.H();
@@ -14,7 +14,7 @@ if (pimple.transonic())
         fvc::interpolate(psi)
        *(
             (fvc::interpolate(rho*HbyA) & mesh.Sf())
-          + Dp*fvc::ddtCorr(rho, U, phi)
+          + rhorAUf*fvc::ddtCorr(rho, U, phi)
         )/fvc::interpolate(rho)
     );
 
@@ -26,7 +26,7 @@ if (pimple.transonic())
         (
             fvm::ddt(psi, p)
           + fvm::div(phid, p)
-          - fvm::laplacian(Dp, p)
+          - fvm::laplacian(rhorAUf, p)
          ==
             fvOptions(psi, p, rho.name())
         );
@@ -48,7 +48,7 @@ else
         "phiHbyA",
         (
             (fvc::interpolate(rho*HbyA) & mesh.Sf())
-          + Dp*fvc::ddtCorr(rho, U, phi)
+          + rhorAUf*fvc::ddtCorr(rho, U, phi)
         )
     );
 
@@ -60,7 +60,7 @@ else
         (
             fvm::ddt(psi, p)
           + fvc::div(phiHbyA)
-          - fvm::laplacian(Dp, p)
+          - fvm::laplacian(rhorAUf, p)
          ==
             fvOptions(psi, p, rho.name())
         );

diff --git a/applications/solvers/combustion/engineFoam/pEqn.H b/applications/solvers/combustion/engineFoam/pEqn.H
@@ -1,7 +1,7 @@
 rho = thermo.rho();
 
 volScalarField rAU(1.0/UEqn.A());
-surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
 
 volVectorField HbyA("HbyA", U);
 HbyA = rAU*UEqn.H();
@@ -15,7 +15,7 @@ if (pimple.transonic())
        *(
             (
                 (fvc::interpolate(rho*HbyA) & mesh.Sf())
-              //***HGW + Dp*fvc::ddtCorr(rho, U, phi)
+              //***HGW + rhorAUf*fvc::ddtCorr(rho, U, phi)
             )/fvc::interpolate(rho)
           - fvc::meshPhi(rho, U)
         )
@@ -51,7 +51,7 @@ else
         "phiHbyA",
         (
             (fvc::interpolate(rho*HbyA) & mesh.Sf())
-          //***HGW + Dp*fvc::ddtCorr(rho, U, phi)
+          //***HGW + rhorAUf*fvc::ddtCorr(rho, U, phi)
         )
       - fvc::interpolate(rho)*fvc::meshPhi(rho, U)
     );

diff --git a/applications/solvers/combustion/fireFoam/fireFoam.C b/applications/solvers/combustion/fireFoam/fireFoam.C
@@ -41,6 +41,7 @@ Description
 #include "psiCombustionModel.H"
 #include "pimpleControl.H"
 #include "fvIOoptionList.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

diff --git a/applications/solvers/combustion/fireFoam/pEqn.H b/applications/solvers/combustion/fireFoam/pEqn.H
@@ -1,27 +1,35 @@
 rho = thermo.rho();
 
 volScalarField rAU(1.0/UEqn.A());
-surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
 
 volVectorField HbyA("HbyA", U);
 HbyA = rAU*UEqn.H();
-phi.boundaryField() =
-    fvc::interpolate(rho.boundaryField()*U.boundaryField())
-  & mesh.Sf().boundaryField();
 
-surfaceScalarField phig(-Dp*ghf*fvc::snGrad(rho)*mesh.magSf());
+surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
 
 surfaceScalarField phiHbyA
 (
     "phiHbyA",
     (
         (fvc::interpolate(rho*HbyA) & mesh.Sf())
-      + Dp*fvc::ddtCorr(rho, U, phi)
+      + rhorAUf*fvc::ddtCorr(rho, U, phi)
     )
   + phig
 );
 
-fvOptions.makeRelative(phiHbyA);
+fvOptions.makeRelative(fvc::interpolate(rho), phiHbyA);
+
+// Update the fixedFluxPressure BCs to ensure flux consistency
+setSnGrad<fixedFluxPressureFvPatchScalarField>
+(
+    p_rgh.boundaryField(),
+    (
+        phiHbyA.boundaryField()
+      - fvOptions.relative(mesh.Sf().boundaryField() & U.boundaryField())
+       *rho.boundaryField()
+    )/(mesh.magSf().boundaryField()*rhorAUf.boundaryField())
+);
 
 while (pimple.correctNonOrthogonal())
 {
@@ -30,7 +38,7 @@ while (pimple.correctNonOrthogonal())
         fvc::ddt(psi, rho)*gh
       + fvc::div(phiHbyA)
       + fvm::ddt(psi, p_rgh)
-      - fvm::laplacian(Dp, p_rgh)
+      - fvm::laplacian(rhorAUf, p_rgh)
      ==
         parcels.Srho()
       + surfaceFilm.Srho()
@@ -44,7 +52,7 @@ while (pimple.correctNonOrthogonal())
     if (pimple.finalNonOrthogonalIter())
     {
         phi = phiHbyA + p_rghEqn.flux();
-        U = HbyA + rAU*fvc::reconstruct((p_rghEqn.flux() + phig)/Dp);
+        U = HbyA + rAU*fvc::reconstruct((p_rghEqn.flux() + phig)/rhorAUf);
         U.correctBoundaryConditions();
         fvOptions.correct(U);
     }

diff --git a/applications/solvers/combustion/reactingFoam/pEqn.H b/applications/solvers/combustion/reactingFoam/pEqn.H
@@ -1,7 +1,7 @@
 rho = thermo.rho();
 
 volScalarField rAU(1.0/UEqn.A());
-surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
 
 volVectorField HbyA("HbyA", U);
 HbyA = rAU*UEqn.H();
@@ -14,7 +14,7 @@ if (pimple.transonic())
         fvc::interpolate(psi)
        *(
             (fvc::interpolate(rho*HbyA) & mesh.Sf())
-          + Dp*fvc::ddtCorr(rho, U, phi)
+          + rhorAUf*fvc::ddtCorr(rho, U, phi)
         )/fvc::interpolate(rho)
     );
 
@@ -48,7 +48,7 @@ else
         "phiHbyA",
         (
             (fvc::interpolate(rho*HbyA) & mesh.Sf())
-          + Dp*fvc::ddtCorr(rho, U, phi)
+          + rhorAUf*fvc::ddtCorr(rho, U, phi)
         )
     );
 

diff --git a/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/pEqn.H b/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/pEqn.H
@@ -6,25 +6,36 @@
     thermo.rho() -= psi*p;
 
     volScalarField rAU(1.0/UEqn.A());
-    surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+    surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
 
-    surfaceScalarField phig(-Dp*ghf*fvc::snGrad(rho)*mesh.magSf());
+    surfaceScalarField phig(-rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf());
 
     surfaceScalarField phiHbyA
     (
         "phiHbyA",
         (
             (fvc::interpolate(rho*HbyA) & mesh.Sf())
-          + Dp*fvc::ddtCorr(rho, U, phi)
+          + rhorAUf*fvc::ddtCorr(rho, U, phi)
         )
       + phig
     );
 
     fvOptions.makeRelative(fvc::interpolate(rho), phiHbyA);
 
+    // Update the fixedFluxPressure BCs to ensure flux consistency
+    setSnGrad<fixedFluxPressureFvPatchScalarField>
+    (
+        p_rgh.boundaryField(),
+        (
+            phiHbyA.boundaryField()
+          - fvOptions.relative(mesh.Sf().boundaryField() & U.boundaryField())
+           *rho.boundaryField()
+        )/(mesh.magSf().boundaryField()*rhorAUf.boundaryField())
+    );
+
     fvScalarMatrix p_rghDDtEqn
     (
         fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh))
@@ -38,7 +49,7 @@
         fvScalarMatrix p_rghEqn
         (
             p_rghDDtEqn
-          - fvm::laplacian(Dp, p_rgh)
+          - fvm::laplacian(rhorAUf, p_rgh)
         );
 
         fvOptions.constrain(p_rghEqn);
@@ -55,7 +66,7 @@
 
             // Correct the momentum source with the pressure gradient flux
             // calculated from the relaxed pressure
-            U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/Dp);
+            U = HbyA + rAU*fvc::reconstruct((phig + p_rghEqn.flux())/rhorAUf);
             U.correctBoundaryConditions();
             fvOptions.correct(U);
             K = 0.5*magSqr(U);

diff --git a/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/rhoReactingBuoyantFoam.C b/applications/solvers/combustion/reactingFoam/rhoReactingBuoyantFoam/rhoReactingBuoyantFoam.C
@@ -36,6 +36,7 @@ Description
 #include "multivariateScheme.H"
 #include "pimpleControl.H"
 #include "fvIOoptionList.H"
+#include "fixedFluxPressureFvPatchScalarField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 

diff --git a/applications/solvers/combustion/reactingFoam/rhoReactingFoam/pEqn.H b/applications/solvers/combustion/reactingFoam/rhoReactingFoam/pEqn.H
@@ -6,7 +6,7 @@
     thermo.rho() -= psi*p;
 
     volScalarField rAU(1.0/UEqn.A());
-    surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+    surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
 
     volVectorField HbyA("HbyA", U);
     HbyA = rAU*UEqn.H();
@@ -18,7 +18,7 @@
             "phiHbyA",
             (
                 (fvc::interpolate(rho*HbyA) & mesh.Sf())
-              + Dp*fvc::ddtCorr(rho, U, phi)
+              + rhorAUf*fvc::ddtCorr(rho, U, phi)
             )/fvc::interpolate(rho)
         );
 
@@ -39,7 +39,7 @@
             fvScalarMatrix pEqn
             (
                 pDDtEqn
-              - fvm::laplacian(Dp, p)
+              - fvm::laplacian(rhorAUf, p)
              ==
                 fvOptions(psi, p, rho.name())
             );
@@ -61,7 +61,7 @@
             "phiHbyA",
             (
                 (fvc::interpolate(rho*HbyA) & mesh.Sf())
-              + Dp*fvc::ddtCorr(rho, U, phi)
+              + rhorAUf*fvc::ddtCorr(rho, U, phi)
             )
         );
 
@@ -80,7 +80,7 @@
             fvScalarMatrix pEqn
             (
                 pDDtEqn
-              - fvm::laplacian(Dp, p)
+              - fvm::laplacian(rhorAUf, p)
             );
 
             fvOptions.constrain(pEqn);

diff --git a/applications/solvers/compressible/rhoPimpleFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/pEqn.H
@@ -4,7 +4,7 @@ rho = min(rho, rhoMax);
 rho.relax();
 
 volScalarField rAU(1.0/UEqn().A());
-surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
 
 volVectorField HbyA("HbyA", U);
 HbyA = rAU*UEqn().H();
@@ -22,7 +22,7 @@ if (pimple.transonic())
         fvc::interpolate(psi)
        *(
             (fvc::interpolate(rho*HbyA) & mesh.Sf())
-          + Dp*fvc::ddtCorr(rho, U, phi)
+          + rhorAUf*fvc::ddtCorr(rho, U, phi)
         )/fvc::interpolate(rho)
     );
 
@@ -34,7 +34,7 @@ if (pimple.transonic())
         (
             fvm::ddt(psi, p)
           + fvm::div(phid, p)
-          - fvm::laplacian(Dp, p)
+          - fvm::laplacian(rhorAUf, p)
           ==
             fvOptions(psi, p, rho.name())
         );
@@ -56,22 +56,20 @@ else
         "phiHbyA",
         (
             (fvc::interpolate(rho*HbyA) & mesh.Sf())
-          + Dp*fvc::ddtCorr(rho, U, phi)
+          + rhorAUf*fvc::ddtCorr(rho, U, phi)
         )
     );
 
     fvOptions.makeRelative(fvc::interpolate(rho), phiHbyA);
 
-    volScalarField Dp("Dp", rho*rAU);
-
     while (pimple.correctNonOrthogonal())
     {
         // Pressure corrector
         fvScalarMatrix pEqn
         (
             fvm::ddt(psi, p)
           + fvc::div(phiHbyA)
-          - fvm::laplacian(Dp, p)
+          - fvm::laplacian(rhorAUf, p)
           ==
             fvOptions(psi, p, rho.name())
         );

diff --git a/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H b/applications/solvers/compressible/rhoPimpleFoam/rhoPimpleDyMFoam/pEqn.H
@@ -4,7 +4,7 @@ rho = min(rho, rhoMax);
 rho.relax();
 
 volScalarField rAU(1.0/UEqn().A());
-surfaceScalarField Dp("Dp", fvc::interpolate(rho*rAU));
+surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU));
 
 volVectorField HbyA("HbyA", U);
 HbyA = rAU*UEqn().H();
@@ -22,7 +22,7 @@ if (pimple.transonic())
         fvc::interpolate(psi)
        *(
             (fvc::interpolate(rho*HbyA) & mesh.Sf())
-          + Dp*fvc::ddtCorr(rho, U, phiAbs)
+          + rhorAUf*fvc::ddtCorr(rho, U, phiAbs)
         )/fvc::interpolate(rho)
     );
 
@@ -34,7 +34,7 @@ if (pimple.transonic())
         (
             fvm::ddt(psi, p)
           + fvm::div(phid, p)
-          - fvm::laplacian(Dp, p)
+          - fvm::laplacian(rhorAUf, p)
           ==
             fvOptions(psi, p, rho.name())
         );
@@ -55,7 +55,7 @@ else
     (
         "phiHbyA",
         (fvc::interpolate(rho*HbyA) & mesh.Sf())
-      + Dp*fvc::ddtCorr(rho, U, phiAbs)
+      + rhorAUf*fvc::ddtCorr(rho, U, phiAbs)
     );
 
     fvOptions.makeRelative(fvc::interpolate(rho), phiHbyA);
@@ -67,7 +67,7 @@ else
         (
             fvm::ddt(psi, p)
           + fvc::div(phiHbyA)
-          - fvm::laplacian(Dp, p)
+          - fvm::laplacian(rhorAUf, p)
           ==
             fvOptions(psi, p, rho.name())
         );