compressibleTwoPhaseEulerFoam: Completed update to energy equations and tutorials

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/EEqns.H

@@ -14,14 +14,18 @@
     volScalarField& he1 = thermo1.he();
     volScalarField& he2 = thermo2.he();
 
+    volScalarField Cpv1(thermo1.Cpv());
+    volScalarField Cpv2(thermo2.Cpv());
+
     fvScalarMatrix he1Eqn
     (
         fvm::ddt(alpha1, he1)
       + fvm::div(alphaPhi1, he1)
       - fvm::laplacian(k1, he1)
      ==
-        heatTransferCoeff*(thermo1.he(p, thermo2.T())/thermo1.Cp())/rho1
-      - fvm::Sp(heatTransferCoeff/thermo1.Cp()/rho1, he1)
+        heatTransferCoeff*(thermo2.T() - thermo1.T())/rho1
+      + heatTransferCoeff*he1/Cpv1/rho1
+      - fvm::Sp(heatTransferCoeff/Cpv1/rho1, he1)
       + alpha1*(dpdt/rho1 - (fvc::ddt(K1) + fvc::div(phi1, K1)))
     );
 
@@ -31,8 +35,9 @@
       + fvm::div(alphaPhi2, he2)
       - fvm::laplacian(k2, he2)
      ==
-        heatTransferCoeff*(thermo2.he(p, thermo1.T())/thermo2.Cp())/rho2
-      - fvm::Sp(heatTransferCoeff/thermo2.Cp()/rho2, he2)
+        heatTransferCoeff*(thermo1.T() - thermo2.T())/rho2
+      + heatTransferCoeff*he2/Cpv2/rho2
+      - fvm::Sp(heatTransferCoeff/Cpv2/rho2, he2)
       + alpha2*(dpdt/rho2 - (fvc::ddt(K2) + fvc::div(phi2, K2)))
     );
 

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/alphaEqn.H

@@ -2,8 +2,8 @@ surfaceScalarField alphaPhi1("alphaPhi" + phase1Name, phi1);
 surfaceScalarField alphaPhi2("alphaPhi" + phase2Name, phi2);
 
 {
-    word alphaScheme("div(phi,alpha)");
-    word alpharScheme("div(phir,alpha)");
+    word alphaScheme("div(phi," + alpha1.name() + ')');
+    word alpharScheme("div(phir," + alpha1.name() + ')');
 
     surfaceScalarField phic("phic", phi);
     surfaceScalarField phir("phir", phi1 - phi2);

applications/solvers/multiphase/multiphaseEulerFoam/Make/options

@@ -1,4 +1,4 @@
-EXE_INC = \
+EXE_INC = -g \
     -IphaseModel/lnInclude \
     -ImultiphaseSystem/lnInclude \
     -ImultiphaseFixedFluxPressure \

applications/solvers/multiphase/multiphaseEulerFoam/UEqns.H

@@ -32,9 +32,10 @@ forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
                 "div(Rc)"
             )
          ==
-          - fvm::Sp(fluid.dragCoeff(phase, dragCoeffs())/phase.rho(), U)
+        //- fvm::Sp(fluid.dragCoeff(phase, dragCoeffs())/phase.rho(), U)
         //- (alpha*phase.rho())*fluid.lift(phase)
-          + (alpha/phase.rho())*fluid.Svm(phase)
+          //+
+            (alpha/phase.rho())*fluid.Svm(phase)
           - fvm::Sp
             (
                 slamDampCoeff
@@ -53,7 +54,7 @@ forAllIter(PtrDictionary<phaseModel>, fluid.phases(), iter)
         alpha*(1 + (1/phase.rho())*fluid.Cvm(phase)),
         UEqns[phasei]
     );
-    UEqns[phasei].relax();
+    //UEqns[phasei].relax();
 
     phasei++;
 }

applications/solvers/multiphase/multiphaseEulerFoam/pEqn.H

@@ -67,8 +67,29 @@
         const volScalarField& alpha = phase;
 
         alphafs.set(phasei, fvc::interpolate(alpha).ptr());
-        rAUs.set(phasei, (1.0/UEqns[phasei].A()).ptr());
-        rAlphaAUfs.set(phasei, fvc::interpolate(alpha*rAUs[phasei]).ptr());
+
+        volScalarField dragCoeffi
+        (
+            IOobject
+            (
+                "dragCoeffi",
+                runTime.timeName(),
+                mesh
+            ),
+            fluid.dragCoeff(phase, dragCoeffs())/phase.rho(),
+            zeroGradientFvPatchScalarField::typeName
+        );
+        dragCoeffi.correctBoundaryConditions();
+
+        rAUs.set(phasei, (1.0/(UEqns[phasei].A() + dragCoeffi)).ptr());
+        rAlphaAUfs.set
+        (
+            phasei,
+            (
+                alphafs[phasei]
+               /fvc::interpolate(UEqns[phasei].A() + dragCoeffi)
+            ).ptr()
+        );
 
         HbyAs[phasei] = rAUs[phasei]*UEqns[phasei].H();
 
@@ -115,10 +136,9 @@
             }
 
             phiHbyAs[phasei] +=
-                fvc::interpolate
-                (
-                    (1.0/phase.rho())*rAUs[phasei]*(*dcIter())
-                )*phase2Ptr->phi();
+                fvc::interpolate((*dcIter())/phase.rho())
+               /fvc::interpolate(UEqns[phasei].A() + dragCoeffi)
+               *phase2Ptr->phi();
 
             HbyAs[phasei] +=
                 (1.0/phase.rho())*rAUs[phasei]*(*dcIter())
@@ -240,7 +260,7 @@
                       + rAlphaAUfs[phasei]*mSfGradp/phase.rho()
                     );
 
-                // phase.U() = fvc::reconstruct(phase.phi());
+                //phase.U() = fvc::reconstruct(phase.phi());
                 phase.U().correctBoundaryConditions();
 
                 U += alpha*phase.U();