multiphaseEuler: Update the virtual-mass force implementation

The central coefficient part of the virtual-mass phase acceleration matrix is now included in the phase velocity transport central coefficient + drag matrix so that the all the phase contributions to each phase momentum equation are handled implicitly and consistently without lagging contribution from the other phases in either the pressure equation or phase momentum correctors. This improves the conditioning of the pressure equation and convergence rate of bubbly-flow cases.
OpenFOAM · Sep 1, 2023 · 18200e7 · 18200e7
1 parent a0c391b
commit 18200e7
Show file tree

Hide file tree

Showing 11 changed files with 274 additions and 378 deletions.
diff --git a/applications/modules/multiphaseEuler/functionObjects/phaseForces/phaseForces.C b/applications/modules/multiphaseEuler/functionObjects/phaseForces/phaseForces.C
@@ -241,7 +241,10 @@ bool Foam::functionObjects::phaseForces::execute()
             *forceFields_[virtualMassModel::typeName] +=
                 fluid_.lookupInterfacialModel
                 <blendedVirtualMassModel>(interface).K()
-               *(otherPhase.DUDt() - phase_.DUDt());
+               *(
+                    (otherPhase.DUDt() & otherPhase.U())
+                  - (phase_.DUDt() & phase_.U())
+                );
         }
 
         if (fluid_.foundInterfacialModel<blendedLiftModel>(interface))

diff --git a/applications/modules/multiphaseEuler/multiphaseEuler/cellPressureCorrector.C b/applications/modules/multiphaseEuler/multiphaseEuler/cellPressureCorrector.C
@@ -68,6 +68,7 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector()
         rAs.setSize(phases.size());
     }
 
+    PtrList<volVectorField> HVms(movingPhases.size());
     PtrList<PtrList<volScalarField>> invADs;
     PtrList<PtrList<surfaceScalarField>> invADfs;
     {
@@ -109,7 +110,7 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector()
             }
         }
 
-        fluid.invADs(As, invADs, invADfs);
+        fluid.invADs(As, HVms, invADs, invADfs);
     }
 
     volScalarField rho("rho", fluid.rho());
@@ -168,7 +169,6 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector()
     // --- Optional momentum predictor
     if (predictMomentum)
     {
-        InfoInFunction << endl;
         PtrList<volVectorField> HbyADs;
         {
             PtrList<volVectorField> Hs(movingPhases.size());
@@ -189,6 +189,11 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector()
                     )
                    *phase.U()().oldTime()
                 );
+
+                if (HVms.set(movingPhasei))
+                {
+                    Hs[movingPhasei] += HVms[movingPhasei];
+                }
             }
 
             HbyADs = invADs & Hs;
@@ -253,6 +258,11 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector()
                    *phase.U()().oldTime()
                 );
 
+                if (HVms.set(movingPhasei))
+                {
+                    Hs[movingPhasei] += HVms[movingPhasei];
+                }
+
                 phiHs.set
                 (
                     movingPhasei,

diff --git a/applications/modules/multiphaseEuler/multiphaseEuler/facePressureCorrector.C b/applications/modules/multiphaseEuler/multiphaseEuler/facePressureCorrector.C
@@ -68,9 +68,9 @@ void Foam::solvers::multiphaseEuler::facePressureCorrector()
         rAs.setSize(phases.size());
     }
 
+    PtrList<surfaceScalarField> HVmfs(movingPhases.size());
     PtrList<PtrList<surfaceScalarField>> invADfs;
     {
-        PtrList<surfaceScalarField> Vmfs(fluid.Vmfs());
         PtrList<surfaceScalarField> Afs(movingPhases.size());
 
         forAll(fluid.movingPhases(), movingPhasei)
@@ -110,14 +110,9 @@ void Foam::solvers::multiphaseEuler::facePressureCorrector()
                     fvc::interpolate(A)
                 )
             );
-
-            if (Vmfs.set(phase.index()))
-            {
-                Afs[movingPhasei] += Vmfs[phase.index()];
-            }
         }
 
-        invADfs = fluid.invADfs(Afs);
+        invADfs = fluid.invADfs(Afs, HVmfs);
     }
 
     volScalarField rho("rho", fluid.rho());
@@ -205,6 +200,11 @@ void Foam::solvers::multiphaseEuler::facePressureCorrector()
                      + fvc::flux(UEqns[phase.index()].H())
                     )
                 );
+
+                if (HVmfs.set(movingPhasei))
+                {
+                    phiHs[movingPhasei] += HVmfs[movingPhasei];
+                }
             }
 
             phiHbyADs = invADfs & phiHs;