twoPhaseEulerFoam: Improved handling of fluxes and HbyA to support totalPressure BCs etc.

applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H

@@ -11,8 +11,11 @@
     rUaAf = fvc::interpolate(rUaA);
     surfaceScalarField rUbAf(fvc::interpolate(rUbA));
 
-    Ua = rUaA*UaEqn.H();
-    Ub = rUbA*UbEqn.H();
+    volVectorField HbyAa("HbyAa", Ua);
+    HbyAa = rUaA*UaEqn.H();
+
+    volVectorField HbyAb("HbyAb", Ub);
+    HbyAb = rUbA*UbEqn.H();
 
     surfaceScalarField phiDraga
     (
@@ -46,24 +49,19 @@
 
     surfaceScalarField phiHbyAa
     (
-        (fvc::interpolate(Ua) & mesh.Sf())
+        (fvc::interpolate(HbyAa) & mesh.Sf())
       + fvc::ddtPhiCorr(rUaA, Ua, phia)
       + phiDraga
     );
 
     surfaceScalarField phiHbyAb
     (
-        (fvc::interpolate(Ub) & mesh.Sf())
+        (fvc::interpolate(HbyAb) & mesh.Sf())
       + fvc::ddtPhiCorr(rUbA, Ub, phib)
       + phiDragb
     );
 
-    surfaceScalarField phi0("phi0", alphaf*phiHbyAa + betaf*phiHbyAb);
-
-    phia = phiHbyAa;
-    phib = phiHbyAb;
-
-    phi = alphaf*phia + betaf*phib;
+    surfaceScalarField phiHbyA("phiHbyA", alphaf*phiHbyAa + betaf*phiHbyAb);
 
     surfaceScalarField Dp
     (
@@ -75,7 +73,7 @@
     {
         fvScalarMatrix pEqn
         (
-            fvm::laplacian(Dp, p) == fvc::div(phi0)
+            fvm::laplacian(Dp, p) == fvc::div(phiHbyA)
         );
 
         pEqn.setReference(pRefCell, pRefValue);
@@ -86,17 +84,17 @@
         {
             surfaceScalarField SfGradp(pEqn.flux()/Dp);
 
-            phia -= rUaAf*SfGradp/rhoa;
-            phib -= rUbAf*SfGradp/rhob;
+            phia = phiHbyAa - rUaAf*SfGradp/rhoa;
+            phib = phiHbyAb - rUbAf*SfGradp/rhob;
             phi = alphaf*phia + betaf*phib;
 
             p.relax();
             SfGradp = pEqn.flux()/Dp;
 
-            Ua += fvc::reconstruct(phiDraga - rUaAf*SfGradp/rhoa);
+            Ua = HbyAa + fvc::reconstruct(phiDraga - rUaAf*SfGradp/rhoa);
             Ua.correctBoundaryConditions();
 
-            Ub += fvc::reconstruct(phiDragb - rUbAf*SfGradp/rhob);
+            Ub = HbyAb + fvc::reconstruct(phiDragb - rUbAf*SfGradp/rhob);
             Ub.correctBoundaryConditions();
 
             U = alpha*Ua + beta*Ub;

src/finiteVolume/fields/fvPatchFields/derived/multiphaseFixedFluxPressure/multiphaseFixedFluxPressureFvPatchScalarField.C

@@ -39,7 +39,7 @@ multiphaseFixedFluxPressureFvPatchScalarField
 )
 :
     fixedGradientFvPatchScalarField(p, iF),
-    phi0Name_("phi0"),
+    phiHbyAName_("phiHbyA"),
     phiName_("phi"),
     rhoName_("rho")
 {}
@@ -55,7 +55,7 @@ multiphaseFixedFluxPressureFvPatchScalarField
 )
 :
     fixedGradientFvPatchScalarField(ptf, p, iF, mapper),
-    phi0Name_(ptf.phi0Name_),
+    phiHbyAName_(ptf.phiHbyAName_),
     phiName_(ptf.phiName_),
     rhoName_(ptf.rhoName_)
 {}
@@ -70,7 +70,7 @@ multiphaseFixedFluxPressureFvPatchScalarField
 )
 :
     fixedGradientFvPatchScalarField(p, iF),
-    phi0Name_(dict.lookupOrDefault<word>("phi0", "phi0")),
+    phiHbyAName_(dict.lookupOrDefault<word>("phiHbyA", "phiHbyA")),
     phiName_(dict.lookupOrDefault<word>("phi", "phi")),
     rhoName_(dict.lookupOrDefault<word>("rho", "rho"))
 {
@@ -95,7 +95,7 @@ multiphaseFixedFluxPressureFvPatchScalarField
 )
 :
     fixedGradientFvPatchScalarField(wbppsf),
-    phi0Name_(wbppsf.phi0Name_),
+    phiHbyAName_(wbppsf.phiHbyAName_),
     phiName_(wbppsf.phiName_),
     rhoName_(wbppsf.rhoName_)
 {}
@@ -109,7 +109,7 @@ multiphaseFixedFluxPressureFvPatchScalarField
 )
 :
     fixedGradientFvPatchScalarField(wbppsf, iF),
-    phi0Name_(wbppsf.phi0Name_),
+    phiHbyAName_(wbppsf.phiHbyAName_),
     phiName_(wbppsf.phiName_),
     rhoName_(wbppsf.rhoName_)
 {}
@@ -124,14 +124,14 @@ void Foam::multiphaseFixedFluxPressureFvPatchScalarField::updateCoeffs()
         return;
     }
 
-    const surfaceScalarField& phi0 =
-        db().lookupObject<surfaceScalarField>(phi0Name_);
+    const surfaceScalarField& phiHbyA =
+        db().lookupObject<surfaceScalarField>(phiHbyAName_);
 
     const surfaceScalarField& phi =
         db().lookupObject<surfaceScalarField>(phiName_);
 
-    fvsPatchField<scalar> phi0p =
-        patch().patchField<surfaceScalarField, scalar>(phi0);
+    fvsPatchField<scalar> phiHbyAp =
+        patch().patchField<surfaceScalarField, scalar>(phiHbyA);
 
     fvsPatchField<scalar> phip =
         patch().patchField<surfaceScalarField, scalar>(phi);
@@ -147,7 +147,7 @@ void Foam::multiphaseFixedFluxPressureFvPatchScalarField::updateCoeffs()
     const fvsPatchField<scalar>& Dpp =
         patch().lookupPatchField<surfaceScalarField, scalar>("Dp");
 
-    gradient() = (phi0p - phip)/patch().magSf()/Dpp;
+    gradient() = (phiHbyAp - phip)/patch().magSf()/Dpp;
 
     fixedGradientFvPatchScalarField::updateCoeffs();
 }
@@ -159,7 +159,7 @@ void Foam::multiphaseFixedFluxPressureFvPatchScalarField::write
 ) const
 {
     fvPatchScalarField::write(os);
-    writeEntryIfDifferent<word>(os, "phi0", "phi0", phi0Name_);
+    writeEntryIfDifferent<word>(os, "phiHbyA", "phiHbyA", phiHbyAName_);
     writeEntryIfDifferent<word>(os, "phi", "phi", phiName_);
     writeEntryIfDifferent<word>(os, "rho", "rho", rhoName_);
     gradient().writeEntry("gradient", os);

src/finiteVolume/fields/fvPatchFields/derived/multiphaseFixedFluxPressure/multiphaseFixedFluxPressureFvPatchScalarField.H

@@ -55,7 +55,7 @@ class multiphaseFixedFluxPressureFvPatchScalarField
     // Private data
 
         //- Name of the predicted flux transporting the field
-        word phi0Name_;
+        word phiHbyAName_;
 
         //- Name of the flux transporting the field
         word phiName_;