twoPhaseEulerFoam: Add fvOptions support and tutorial

applications/solvers/multiphase/twoPhaseEulerFoam/EEqns.H

@@ -10,10 +10,10 @@
     fvScalarMatrix he1Eqn
     (
         fvm::ddt(alpha1, rho1, he1) + fvm::div(alphaRhoPhi1, he1)
-      - fvm::Sp(fvc::ddt(alpha1, rho1) + fvc::div(alphaRhoPhi1), he1)
+      - fvm::Sp(contErr1, he1)
 
       + fvc::ddt(alpha1, rho1, K1) + fvc::div(alphaRhoPhi1, K1)
-      - (fvc::ddt(alpha1, rho1) + fvc::div(alphaRhoPhi1))*K1
+      - contErr1*K1
 
       + (
             he1.name() == thermo1.phasePropertyName("e")
@@ -33,15 +33,16 @@
         heatTransferCoeff*(thermo2.T() - thermo1.T())
       + heatTransferCoeff*he1/Cpv1
       - fvm::Sp(heatTransferCoeff/Cpv1, he1)
+      + fvOptions(alpha1, rho1, he1)
     );
 
     fvScalarMatrix he2Eqn
     (
         fvm::ddt(alpha2, rho2, he2) + fvm::div(alphaRhoPhi2, he2)
-      - fvm::Sp(fvc::ddt(alpha2, rho2) + fvc::div(alphaRhoPhi2), he2)
+      - fvm::Sp(contErr2, he2)
 
       + fvc::ddt(alpha2, rho2, K2) + fvc::div(alphaRhoPhi2, K2)
-      - (fvc::ddt(alpha2, rho2) + fvc::div(alphaRhoPhi2))*K2
+      - contErr2*K2
 
       + (
             he2.name() == thermo2.phasePropertyName("e")
@@ -61,12 +62,15 @@
         heatTransferCoeff*(thermo1.T() - thermo2.T())
       + heatTransferCoeff*he2/Cpv2
       - fvm::Sp(heatTransferCoeff/Cpv2, he2)
+      + fvOptions(alpha2, rho2, he2)
     );
 
     he1Eqn.relax();
+    fvOptions.constrain(he1Eqn);
     he1Eqn.solve();
 
     he2Eqn.relax();
+    fvOptions.constrain(he2Eqn);
     he2Eqn.solve();
 
     thermo1.correct();

applications/solvers/multiphase/twoPhaseEulerFoam/Make/options

@@ -9,7 +9,10 @@ EXE_INC = \
     -IphaseCompressibleTurbulenceModels/lnInclude \
     -IinterfacialModels/lnInclude \
     -ItwoPhaseSystem/lnInclude \
-    -Iaveraging
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/fvOptions/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude
 
 EXE_LIBS = \
     -lfluidThermophysicalModels \
@@ -21,4 +24,6 @@ EXE_LIBS = \
     -lcompressibleTwoPhaseSystem \
     -lcompressibleEulerianInterfacialModels \
     -lfiniteVolume \
-    -lmeshTools
+    -lfvOptions \
+    -lmeshTools \
+    -lsampling

applications/solvers/multiphase/twoPhaseEulerFoam/UEqns.H

@@ -9,7 +9,6 @@ volScalarField dragCoeff(fluid.dragCoeff());
 
 {
     volScalarField virtualMassCoeff(fluid.virtualMassCoeff());
-
     volVectorField liftForce(fluid.liftForce());
     volVectorField wallLubricationForce(fluid.wallLubricationForce());
     volVectorField turbulentDispersionForce(fluid.turbulentDispersionForce());
@@ -19,15 +18,14 @@ volScalarField dragCoeff(fluid.dragCoeff());
         (
             fvm::ddt(alpha1, rho1, U1)
           + fvm::div(alphaRhoPhi1, U1)
-          - fvm::Sp(fvc::ddt(alpha1, rho1) + fvc::div(alphaRhoPhi1), U1)
+          - fvm::Sp(contErr1, U1)
+          + mrfZones(alpha1*rho1 + virtualMassCoeff, U1)
           + phase1.turbulence().divDevRhoReff(U1)
          ==
           - fvm::Sp(dragCoeff, U1)
-          - (
-                liftForce
-              + wallLubricationForce
-              + turbulentDispersionForce
-            )
+          - liftForce
+          - wallLubricationForce
+          - turbulentDispersionForce
           - virtualMassCoeff
            *(
                 fvm::ddt(U1)
@@ -35,9 +33,10 @@ volScalarField dragCoeff(fluid.dragCoeff());
               - fvm::Sp(fvc::div(phi1), U1)
               - DDtU2
             )
+          + fvOptions(alpha1, rho1, U1)
         );
-        mrfZones.addCoriolis(alpha1*rho1 + virtualMassCoeff, U1Eqn);
         U1Eqn.relax();
+        fvOptions.constrain(U1Eqn);
     }
 
     {
@@ -45,15 +44,14 @@ volScalarField dragCoeff(fluid.dragCoeff());
         (
             fvm::ddt(alpha2, rho2, U2)
           + fvm::div(alphaRhoPhi2, U2)
-          - fvm::Sp(fvc::ddt(alpha2, rho2) + fvc::div(alphaRhoPhi2), U2)
+          - fvm::Sp(contErr2, U2)
+          + mrfZones(alpha2*rho2 + virtualMassCoeff, U2)
           + phase2.turbulence().divDevRhoReff(U2)
          ==
           - fvm::Sp(dragCoeff, U2)
-          + (
-                liftForce
-              + wallLubricationForce
-              + turbulentDispersionForce
-            )
+          + liftForce
+          + wallLubricationForce
+          + turbulentDispersionForce
           - virtualMassCoeff
            *(
                 fvm::ddt(U2)
@@ -61,8 +59,9 @@ volScalarField dragCoeff(fluid.dragCoeff());
               - fvm::Sp(fvc::div(phi2), U2)
               - DDtU1
             )
+          + fvOptions(alpha2, rho2, U2)
         );
-        mrfZones.addCoriolis(alpha2*rho2 + virtualMassCoeff, U2Eqn);
         U2Eqn.relax();
+        fvOptions.constrain(U2Eqn);
     }
 }

applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H

@@ -130,6 +130,7 @@
         pEqnComp1 =
             (
                 fvc::ddt(alpha1, rho1) + fvc::div(alphaRhoPhi1)
+              - (fvOptions(alpha1, rho1)&rho1)
               - fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1)
             )/rho1
           + (alpha1/rho1)*correction
@@ -143,7 +144,9 @@
         pEqnComp2 =
             (
                 fvc::ddt(alpha2, rho2) + fvc::div(alphaRhoPhi2)
+              - (fvOptions(alpha2, rho2)&rho2)
               - fvc::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), rho2)
+
             )/rho2
           + (alpha2/rho2)*correction
             (
@@ -158,6 +161,7 @@
         pEqnComp1 =
             (
                 fvc::ddt(alpha1, rho1) + fvc::div(alphaRhoPhi1)
+              - (fvOptions(alpha1, rho1)&rho1)
               - fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1)
             )/rho1
           + (alpha1*psi1/rho1)*correction(fvm::ddt(p));
@@ -165,6 +169,7 @@
         pEqnComp2 =
             (
                 fvc::ddt(alpha2, rho2) + fvc::div(alphaRhoPhi2)
+              - (fvOptions(alpha2, rho2)&rho2)
               - fvc::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), rho2)
             )/rho2
           + (alpha2*psi2/rho2)*correction(fvm::ddt(p));
@@ -227,6 +232,7 @@
                   - phiP1
                 );
             U1.correctBoundaryConditions();
+            fvOptions.correct(U1);
 
             U2 = HbyA2
               + fvc::reconstruct
@@ -239,6 +245,7 @@
                   - phiP2
                 );
             U2.correctBoundaryConditions();
+            fvOptions.correct(U2);
 
             U = fluid.U();
         }

applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C

@@ -35,6 +35,7 @@ Description
 #include "PhaseCompressibleTurbulenceModel.H"
 #include "pimpleControl.H"
 #include "IOMRFZoneList.H"
+#include "fvIOoptionList.H"
 #include "fixedFluxPressureFvPatchScalarField.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
@@ -48,6 +49,7 @@ int main(int argc, char *argv[])
     #include "readGravitationalAcceleration.H"
     #include "createFields.H"
     #include "createMRFZones.H"
+    #include "createFvOptions.H"
     #include "initContinuityErrs.H"
     #include "readTimeControls.H"
     #include "CourantNos.H"
@@ -74,8 +76,21 @@ int main(int argc, char *argv[])
             fluid.solve();
             fluid.correct();
 
-            #include "EEqns.H"
+            volScalarField contErr1
+            (
+                fvc::ddt(alpha1, rho1) + fvc::div(alphaRhoPhi1)
+              - (fvOptions(alpha1, rho1)&rho1)
+            );
+
+            volScalarField contErr2
+            (
+                fvc::ddt(alpha2, rho2) + fvc::div(alphaRhoPhi2)
+               - (fvOptions(alpha2, rho2)&rho2)
+            );
+
+
             #include "UEqns.H"
+            #include "EEqns.H"
 
             // --- Pressure corrector loop
             while (pimple.correct())