twoPhaseSolver: Store the flux of the non-solved-for phase

This flux is needed for boundary conditions, post-processing and
Euler-Euler-like sub-models and functions

applications/modules/compressibleVoF/compressibleInterPhaseTransportModel/compressibleInterPhaseTransportModel.C

@@ -35,6 +35,7 @@ Foam::compressibleInterPhaseTransportModel::compressibleInterPhaseTransportModel
     const surfaceScalarField& phi,
     const surfaceScalarField& rhoPhi,
     const surfaceScalarField& alphaPhi1,
+    const surfaceScalarField& alphaPhi2,
     const surfaceScalarField& alphaRhoPhi1,
     const surfaceScalarField& alphaRhoPhi2,
     const compressibleTwoPhaseVoFMixture& mixture
@@ -44,6 +45,7 @@ Foam::compressibleInterPhaseTransportModel::compressibleInterPhaseTransportModel
     mixture_(mixture),
     phi_(phi),
     alphaPhi1_(alphaPhi1),
+    alphaPhi2_(alphaPhi2),
     alphaRhoPhi1_(alphaRhoPhi1),
     alphaRhoPhi2_(alphaRhoPhi2)
 {

applications/modules/compressibleVoF/compressibleInterPhaseTransportModel/compressibleInterPhaseTransportModel.H

@@ -72,13 +72,16 @@ class compressibleInterPhaseTransportModel
         //- Mixture volumetric flux
         const surfaceScalarField& phi_;
 
-        //- Phase volumetric flux
+        //- Phase-1 volumetric flux
         const surfaceScalarField& alphaPhi1_;
 
-        //- Phase-1 mass-flux (constructed for two-phase transport)
+        //- Phase-2 volumetric flux
+        const surfaceScalarField& alphaPhi2_;
+
+        //- Phase-1 mass flux
         const surfaceScalarField& alphaRhoPhi1_;
 
-        //- Phase-2 mass-flux (constructed for two-phase transport)
+        //- Phase-2 mass flux
         const surfaceScalarField& alphaRhoPhi2_;
 
         //- Mixture transport model (constructed for mixture transport)
@@ -105,6 +108,7 @@ public:
             const surfaceScalarField& phi,
             const surfaceScalarField& rhoPhi,
             const surfaceScalarField& alphaPhi1,
+            const surfaceScalarField& alphaPhi2,
             const surfaceScalarField& alphaRhoPhi1,
             const surfaceScalarField& alphaRhoPhi2,
             const compressibleTwoPhaseVoFMixture& mixture

applications/modules/compressibleVoF/compressibleVoF.C

@@ -88,7 +88,7 @@ Foam::solvers::compressibleVoF::compressibleVoF(fvMesh& mesh)
     alphaRhoPhi2
     (
         IOobject::groupName("alphaRhoPhi", alpha2.group()),
-        fvc::interpolate(mixture_.thermo2().rho())*(phi - alphaPhi1)
+        fvc::interpolate(mixture_.thermo2().rho())*alphaPhi2
     ),
 
     K("K", 0.5*magSqr(U)),
@@ -100,6 +100,7 @@ Foam::solvers::compressibleVoF::compressibleVoF(fvMesh& mesh)
         phi,
         rhoPhi,
         alphaPhi1,
+        alphaPhi2,
         alphaRhoPhi1,
         alphaRhoPhi2,
         mixture_
@@ -144,7 +145,7 @@ void Foam::solvers::compressibleVoF::prePredictor()
     const volScalarField& rho2 = mixture_.thermo2().rho();
 
     alphaRhoPhi1 = fvc::interpolate(rho1)*alphaPhi1;
-    alphaRhoPhi2 = fvc::interpolate(rho2)*(phi - alphaPhi1);
+    alphaRhoPhi2 = fvc::interpolate(rho2)*alphaPhi2;
 
     rhoPhi = alphaRhoPhi1 + alphaRhoPhi2;
 

applications/modules/compressibleVoF/pressureCorrector.C

@@ -54,8 +54,6 @@ void Foam::solvers::compressibleVoF::pressureCorrector()
 
     const surfaceScalarField rAUf("rAUf", fvc::interpolate(rAU()));
 
-    const surfaceScalarField alphaPhi2("alphaPhi2", phi - alphaPhi1);
-
     while (pimple.correct())
     {
         const volVectorField HbyA(constrainHbyA(rAU()*UEqn.H(), U, p_rgh));

applications/modules/incompressibleDriftFlux/incompressibleDriftFlux.C

@@ -193,7 +193,9 @@ void Foam::solvers::incompressibleDriftFlux::prePredictor()
         alpha1Eqn.solve(alpha1.name() + "Diffusion");
 
         alphaPhi1 += alpha1Eqn.flux();
-        alpha2 = 1.0 - alpha1;
+
+        alpha2 = scalar(1) - alpha1;
+        alphaPhi2 = phi - alphaPhi1;
 
         Info<< "Phase-1 volume fraction = "
             << alpha1.weightedAverage(mesh.Vsc()).value()
@@ -207,8 +209,8 @@ void Foam::solvers::incompressibleDriftFlux::prePredictor()
     const dimensionedScalar& rho1 = mixture.rho1();
     const dimensionedScalar& rho2 = mixture.rho2();
 
-    // Calculate the mass-flux from the accumulated alphaPhi1
+    // Calculate the mass-flux
-    rhoPhi = (alphaPhi1*(rho1 - rho2) + phi*rho2);
+    rhoPhi = alphaPhi1*rho1 + alphaPhi2*rho2;
 
     relativeVelocity->correct();
 

applications/modules/incompressibleVoF/incompressibleInterPhaseTransportModel/incompressibleInterPhaseTransportModel.C

@@ -34,13 +34,15 @@ incompressibleInterPhaseTransportModel
     const volVectorField& U,
     const surfaceScalarField& phi,
     const surfaceScalarField& alphaPhi1,
+    const surfaceScalarField& alphaPhi2,
     const incompressibleTwoPhaseVoFMixture& mixture
 )
 :
     twoPhaseTransport_(false),
     mixture_(mixture),
     phi_(phi),
-    alphaPhi1_(alphaPhi1)
+    alphaPhi1_(alphaPhi1),
+    alphaPhi2_(alphaPhi2)
 {
     {
         IOdictionary momentumTransport
@@ -71,15 +73,6 @@ incompressibleInterPhaseTransportModel
         const volScalarField& alpha1(mixture_.alpha1());
         const volScalarField& alpha2(mixture_.alpha2());
 
-        alphaPhi2_ =
-        (
-            new surfaceScalarField
-            (
-                IOobject::groupName("alphaPhi", alpha2.group()),
-                (phi_ - alphaPhi1_)
-            )
-        );
-
         turbulence1_ =
         (
             phaseIncompressible::momentumTransportModel::New
@@ -98,7 +91,7 @@ incompressibleInterPhaseTransportModel
             (
                 alpha2,
                 U,
-                alphaPhi2_(),
+                alphaPhi2_,
                 phi,
                 mixture.nuModel2()
             )
@@ -140,15 +133,6 @@ Foam::incompressibleInterPhaseTransportModel::divDevTau
 }
 
 
-void Foam::incompressibleInterPhaseTransportModel::correctPhasePhi()
-{
-    if (twoPhaseTransport_)
-    {
-        alphaPhi2_.ref() = (phi_ - alphaPhi1_);
-    }
-}
-
-
 void Foam::incompressibleInterPhaseTransportModel::predict()
 {
     if (twoPhaseTransport_)

applications/modules/incompressibleVoF/incompressibleInterPhaseTransportModel/incompressibleInterPhaseTransportModel.H

@@ -69,11 +69,11 @@ class incompressibleInterPhaseTransportModel
         //- Mixture volumetric flux
         const surfaceScalarField& phi_;
 
-        //- Phase volumetric flux
+        //- Phase-1 volumetric flux
         const surfaceScalarField& alphaPhi1_;
 
-        //- Phase-2 mass-flux (constructed from alphaPhi1_ and phi_)
+        //- Phase-2 volumetric flux
-        tmp<surfaceScalarField> alphaPhi2_;
+        const surfaceScalarField& alphaPhi2_;
 
         //- Mixture transport model (constructed for mixture transport)
         autoPtr<incompressible::momentumTransportModel> turbulence_;
@@ -95,6 +95,7 @@ public:
             const volVectorField& U,
             const surfaceScalarField& phi,
             const surfaceScalarField& alphaPhi1,
+            const surfaceScalarField& alphaPhi2,
             const incompressibleTwoPhaseVoFMixture& mixture
         );
 
@@ -114,10 +115,6 @@ public:
             volVectorField& U
         ) const;
 
-        //- Correct the phase mass-fluxes
-        // (required for the two-phase transport option)
-        void correctPhasePhi();
-
         //- Predict the phase or mixture transport models
         void predict();
 

applications/modules/incompressibleVoF/incompressibleVoF.C

@@ -81,6 +81,7 @@ Foam::solvers::incompressibleVoF::incompressibleVoF(fvMesh& mesh)
         U,
         phi,
         alphaPhi1,
+        alphaPhi2,
         mixture
     )
 {
@@ -134,8 +135,8 @@ void Foam::solvers::incompressibleVoF::prePredictor()
     const dimensionedScalar& rho1 = mixture.rho1();
     const dimensionedScalar& rho2 = mixture.rho2();
 
-    // Calculate the mass-flux from the accumulated alphaPhi1
+    // Calculate the mass-flux
-    rhoPhi = (alphaPhi1*(rho1 - rho2) + phi*rho2);
+    rhoPhi = alphaPhi1*rho1 + alphaPhi2*rho2;
 
     if (pimple.predictTransport())
     {

applications/modules/twoPhaseSolver/alphaPredictor.C

@@ -225,7 +225,8 @@ void Foam::solvers::twoPhaseSolver::alphaSolve
         // Cache the upwind-flux
         talphaPhi1Corr0 = talphaPhi1UD;
 
-        alpha2 = 1.0 - alpha1;
+        alpha2 = scalar(1) - alpha1;
+        alphaPhi2 = phi - alphaPhi1;
 
         correctInterface();
     }
@@ -328,7 +329,8 @@ void Foam::solvers::twoPhaseSolver::alphaSolve
             }
         }
 
-        alpha2 = 1.0 - alpha1;
+        alpha2 = scalar(1) - alpha1;
+        alphaPhi2 = phi - alphaPhi1;
 
         // Correct only the mixture interface for the interface compression flux
         correctInterface();
@@ -363,6 +365,7 @@ void Foam::solvers::twoPhaseSolver::alphaSolve
             // Calculate the end-of-time-step alpha flux
             alphaPhi1 =
                 (alphaPhi1 - (1.0 - cnCoeff)*alphaPhi1.oldTime())/cnCoeff;
+            alphaPhi2 = phi - alphaPhi1;
         }
     }
 
@@ -420,6 +423,7 @@ void Foam::solvers::twoPhaseSolver::alphaPredictor()
         }
 
         alphaPhi1 = talphaPhi1();
+        alphaPhi2 = phi - talphaPhi1();
     }
     else
     {

applications/modules/twoPhaseSolver/twoPhaseSolver.C

@@ -74,6 +74,16 @@ Foam::solvers::twoPhaseSolver::twoPhaseSolver
             IOobject::AUTO_WRITE
         ),
         phi*fvc::interpolate(alpha1)
+    ),
+    alphaPhi2
+    (
+        IOobject
+        (
+            IOobject::groupName("alphaPhi", alpha2.group()),
+            runTime.name(),
+            mesh
+        ),
+        phi - alphaPhi1
     )
 {
     mesh.schemes().setFluxRequired(alpha1.name());

applications/modules/twoPhaseSolver/twoPhaseSolver.H

@@ -95,6 +95,9 @@ protected:
         // Phase-1 volumetric flux
         surfaceScalarField alphaPhi1;
 
+        // Phase-2 volumetric flux
+        surfaceScalarField alphaPhi2;
+
 
     // Cached temporary fields