multiphaseEuler: Improved mass transfer linearisation w.r.t pressure

Mass transfer rates are now updated following a change in the pressure
if the mass transfer modelling provides a pressure coefficient. This
means that pimple correctors can be used to improve the behaviour of
mass transfer processes that coupled closely to the pressure field.

applications/solvers/modules/fluid/multiphaseEuler/multiphaseEuler/cellPressureCorrector.C

@@ -102,6 +102,10 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector()
     // Explicit force fluxes
     PtrList<surfaceScalarField> phiFs(fluid.phiFs(rAUs));
 
+    // Mass transfer rates
+    PtrList<volScalarField> dmdts(fluid.dmdts());
+    PtrList<volScalarField> d2mdtdps(fluid.d2mdtdps());
+
     // --- Pressure corrector loop
     while (pimple.correct())
     {
@@ -282,7 +286,7 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector()
         }
 
         // Compressible pressure equations
-        PtrList<fvScalarMatrix> pEqnComps(compressibilityEqns());
+        PtrList<fvScalarMatrix> pEqnComps(compressibilityEqns(dmdts, d2mdtdps));
 
         // Cache p prior to solve for density update
         volScalarField p_rgh_0(p_rgh);
@@ -407,6 +411,15 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector()
             phase.thermoRef().rho() += phase.thermo().psi()*(p_rgh - p_rgh_0);
         }
 
+        // Update mass transfer rates for change in p_rgh
+        forAll(phases, phasei)
+        {
+            if (dmdts.set(phasei) && d2mdtdps.set(phasei))
+            {
+                dmdts[phasei] += d2mdtdps[phasei]*(p_rgh - p_rgh_0);
+            }
+        }
+
         // Correct p_rgh for consistency with p and the updated densities
         rho = fluid.rho();
         p_rgh = p - rho*buoyancy.gh;

applications/solvers/modules/fluid/multiphaseEuler/multiphaseEuler/compressibilityEqns.C

@@ -34,13 +34,14 @@ License
 // * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
 
 Foam::PtrList<Foam::fvScalarMatrix>
-Foam::solvers::multiphaseEuler::compressibilityEqns() const
+Foam::solvers::multiphaseEuler::compressibilityEqns
+(
+    const PtrList<volScalarField>& dmdts,
+    const PtrList<volScalarField>& d2mdtdps
+) const
 {
     PtrList<fvScalarMatrix> pEqnComps(phases.size());
 
-    PtrList<volScalarField> dmdts(fluid.dmdts());
-    PtrList<volScalarField> d2mdtdps(fluid.d2mdtdps());
-
     forAll(phases, phasei)
     {
         phaseModel& phase = phases[phasei];

applications/solvers/modules/fluid/multiphaseEuler/multiphaseEuler/facePressureCorrector.C

@@ -113,6 +113,10 @@ void Foam::solvers::multiphaseEuler::facePressureCorrector()
     // Explicit force fluxes
     PtrList<surfaceScalarField> phiFfs(fluid.phiFfs(rAUfs));
 
+    // Mass transfer rates
+    PtrList<volScalarField> dmdts(fluid.dmdts());
+    PtrList<volScalarField> d2mdtdps(fluid.d2mdtdps());
+
     // --- Pressure corrector loop
     while (pimple.correct())
     {
@@ -269,7 +273,7 @@ void Foam::solvers::multiphaseEuler::facePressureCorrector()
         }
 
         // Compressible pressure equations
-        PtrList<fvScalarMatrix> pEqnComps(compressibilityEqns());
+        PtrList<fvScalarMatrix> pEqnComps(compressibilityEqns(dmdts, d2mdtdps));
 
         // Cache p prior to solve for density update
         volScalarField p_rgh_0(p_rgh);
@@ -381,6 +385,15 @@ void Foam::solvers::multiphaseEuler::facePressureCorrector()
             phase.thermoRef().rho() += phase.thermo().psi()*(p_rgh - p_rgh_0);
         }
 
+        // Update mass transfer rates for change in p_rgh
+        forAll(phases, phasei)
+        {
+            if (dmdts.set(phasei) && d2mdtdps.set(phasei))
+            {
+                dmdts[phasei] += d2mdtdps[phasei]*(p_rgh - p_rgh_0);
+            }
+        }
+
         // Correct p_rgh for consistency with p and the updated densities
         rho = fluid.rho();
         p_rgh = p - rho*buoyancy.gh;

applications/solvers/modules/fluid/multiphaseEuler/multiphaseEuler/multiphaseEuler.H

@@ -185,7 +185,11 @@ private:
         void facePressureCorrector();
 
         //- Return the list of pressure equation compressibility contributions
-        PtrList<fvScalarMatrix> compressibilityEqns() const;
+        PtrList<fvScalarMatrix> compressibilityEqns
+        (
+            const PtrList<volScalarField>& dmdts,
+            const PtrList<volScalarField>& d2mdtdps
+        ) const;
 
 
 public:

applications/solvers/modules/fluid/multiphaseEuler/multiphaseEuler/pressureCorrector.C

@@ -44,7 +44,10 @@ void Foam::solvers::multiphaseEuler::pressureCorrector()
     }
     else
     {
-        PtrList<fvScalarMatrix> pEqnComps(compressibilityEqns());
+        PtrList<fvScalarMatrix> pEqnComps
+        (
+            compressibilityEqns(fluid.dmdts(), fluid.d2mdtdps())
+        );
 
         forAll(phases, phasei)
         {