multiphaseEulerFoam: Added correctPhi support for moving-mesh cases with cell-volume change

applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam/multiphaseEulerFoam.C

@@ -125,6 +125,26 @@ int main(int argc, char *argv[])
             {
                 if (pimple.firstPimpleIter() || moveMeshOuterCorrectors)
                 {
+                    // Store divU from the previous mesh so that it can be
+                    // mapped and used in correctPhi to ensure the corrected phi
+                    // has the same divergence
+                    tmp<volScalarField> divU;
+
+                    if
+                    (
+                        correctPhi
+                    )
+                    {
+                        divU = volScalarField::New
+                        (
+                            "divU0",
+                            fvc::div
+                            (
+                                fvc::absolute(phi, fluid.movingPhases()[0].U())
+                            )
+                        );
+                    }
+
                     mesh.update();
 
                     if (mesh.changing())
@@ -134,6 +154,11 @@ int main(int argc, char *argv[])
 
                         fluid.meshUpdate();
 
+                        if (correctPhi)
+                        {
+                            fluid.correctPhi(p_rgh, divU, pimple);
+                        }
+
                         if (checkMeshCourantNo)
                         {
                             #include "meshCourantNo.H"

applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam/pU/UEqns.H

@@ -31,6 +31,7 @@ PtrList<fvVectorMatrix> UEqns(phases.size());
 
         UEqns[phase.index()].relax();
         fvOptions.constrain(UEqns[phase.index()]);
+        U.correctBoundaryConditions();
         fvOptions.correct(U);
     }
 }

applications/solvers/multiphase/multiphaseEulerFoam/multiphaseEulerFoam/pU/pEqn.H

@@ -183,7 +183,7 @@ while (pimple.correct())
     }
 
     MRF.makeRelative(phiHbyA);
-    fvc::makeRelative(phiHbyA, phases[0].U());
+    fvc::makeRelative(phiHbyA, fluid.movingPhases()[0].U());
 
     // Construct pressure "diffusivity"
     surfaceScalarField rAUf

applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseSystem/phaseSystem.C

@@ -32,11 +32,14 @@ License
 #include "fvcDiv.H"
 #include "fvcGrad.H"
 #include "fvcSnGrad.H"
+#include "CorrectPhi.H"
+#include "fvcMeshPhi.H"
 #include "alphaContactAngleFvPatchScalarField.H"
 #include "unitConversion.H"
 #include "dragModel.H"
 #include "BlendedInterfacialModel.H"
 #include "movingWallVelocityFvPatchVectorField.H"
+#include "pimpleControl.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -887,6 +890,76 @@ void Foam::phaseSystem::correctBoundaryFlux()
 }
 
 
+void Foam::phaseSystem::correctPhi
+(
+    const volScalarField& p_rgh,
+    const tmp<volScalarField>& divU,
+    nonOrthogonalSolutionControl& pimple
+)
+{
+    forAll(movingPhases(), movingPhasei)
+    {
+        phaseModel& phase = movingPhases()[movingPhasei];
+
+        volVectorField::Boundary& Ubf = phase.URef().boundaryFieldRef();
+        surfaceVectorField::Boundary& UfBf = phase.UfRef().boundaryFieldRef();
+
+        forAll(Ubf, patchi)
+        {
+            if (Ubf[patchi].fixesValue())
+            {
+                Ubf[patchi].initEvaluate();
+            }
+        }
+
+        forAll(Ubf, patchi)
+        {
+            if (Ubf[patchi].fixesValue())
+            {
+                Ubf[patchi].evaluate();
+                UfBf[patchi] = Ubf[patchi];
+            }
+        }
+    }
+
+    // Correct fixed-flux BCs to be consistent with the velocity BCs
+    correctBoundaryFlux();
+
+    {
+        phi_ = Zero;
+        PtrList<surfaceScalarField> alphafs(phaseModels_.size());
+        forAll(movingPhases(), movingPhasei)
+        {
+            phaseModel& phase = movingPhases()[movingPhasei];
+            const label phasei = phase.index();
+            const volScalarField& alpha = phase;
+
+            alphafs.set(phasei, fvc::interpolate(alpha).ptr());
+
+            // Calculate absolute flux
+            // from the mapped surface velocity
+            phi_ += alphafs[phasei]*(mesh_.Sf() & phase.Uf());
+        }
+
+        CorrectPhi
+        (
+            phi_,
+            movingPhases()[0].U(),
+            p_rgh,
+            // surfaceScalarField("rAUf", fvc::interpolate(rAU())),
+            dimensionedScalar(dimTime/dimDensity, 1),
+            divU(),
+            pimple
+        );
+
+        // Make the flux relative to the mesh motion
+        fvc::makeRelative(phi_, movingPhases()[0].U());
+
+        setMixturePhi(alphafs, phi_);
+    }
+}
+
+
 bool Foam::phaseSystem::read()
 {
     if (regIOobject::read())

applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/phaseSystem/phaseSystem.H

@@ -60,6 +60,7 @@ class blendingMethod;
 template<class modelType> class BlendedInterfacialModel;
 class surfaceTensionModel;
 class aspectRatioModel;
+class nonOrthogonalSolutionControl;
 
 /*---------------------------------------------------------------------------*\
                          Class phaseSystem Declaration
@@ -274,14 +275,6 @@ protected:
         //  around the specified mixture mean
         void setMixtureU(const volVectorField& Um);
 
-        //- Re-normalise the flux of the phases
-        //  around the specified mixture mean
-        void setMixturePhi
-        (
-            const PtrList<surfaceScalarField>& alphafs,
-            const surfaceScalarField& phim
-        );
-
 
         // Functions required for interface compression
 
@@ -625,6 +618,14 @@ public:
                 const PtrList<surfaceScalarField>& phiKdPhifs
             ) = 0;
 
+            //- Re-normalise the flux of the phases
+            //  around the specified mixture mean
+            void setMixturePhi
+            (
+                const PtrList<surfaceScalarField>& alphafs,
+                const surfaceScalarField& phim
+            );
+
             //- Return the flux corrections for the cell-based algorithm
             virtual PtrList<surfaceScalarField> ddtCorrByAs
             (
@@ -684,6 +685,13 @@ public:
             //- Correct fixed-flux BCs to be consistent with the velocity BCs
             void correctBoundaryFlux();
 
+            void correctPhi
+            (
+                const volScalarField& p_rgh,
+                const tmp<volScalarField>& divU,
+                nonOrthogonalSolutionControl& pimple
+            );
+
 
         // IO
 

tutorials/multiphase/multiphaseEulerFoam/laminar/mixerVesselAMI2D/system/fvSolution

@@ -79,6 +79,7 @@ PIMPLE
     nOuterCorrectors     1;
     nCorrectors     3;
     nNonOrthogonalCorrectors 0;
+    correctPhi      yes;
 
     pRefCell        0;
     pRefValue       0;