compressibleTwoPhaseEulerFoam: Complete support for sensible internal energy and update tutorials

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/EEqns.H

@@ -19,34 +19,46 @@
 
     fvScalarMatrix he1Eqn
     (
-        fvm::ddt(alpha1, he1)
-      + fvm::div(alphaPhi1, he1)
+        fvm::ddt(alpha1, he1) + fvm::div(alphaPhi1, he1)
+      + fvc::ddt(alpha1, K1) + fvc::div(alphaPhi1, K1)
 
         // Compressibity correction
       - fvm::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), he1)
+      - (fvc::ddt(alpha1) + fvc::div(alphaPhi1))*K1
+
+      + (
+            he1.name() == "e"
+          ? fvc::div(alphaPhi1, p)
+          : -dalpha1pdt
+        )/rho1
 
       - fvm::laplacian(k1, 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)))
     );
 
     fvScalarMatrix he2Eqn
     (
-        fvm::ddt(alpha2, he2)
-      + fvm::div(alphaPhi2, he2)
+        fvm::ddt(alpha2, he2) + fvm::div(alphaPhi2, he2)
+      + fvc::ddt(alpha2, K2) + fvc::div(alphaPhi2, K2)
 
         // Compressibity correction
       - fvm::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), he2)
+      - (fvc::ddt(alpha2) + fvc::div(alphaPhi2))*K2
+
+      + (
+            he2.name() == "e"
+          ? fvc::div(alphaPhi2, p)
+          : -dalpha2pdt
+        )/rho2
 
       - fvm::laplacian(k2, 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)))
     );
 
     he1Eqn.relax();

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/createFields.H

@@ -275,8 +275,8 @@
     );
 
 
-    Info<< "Creating field dpdt\n" << endl;
-    volScalarField dpdt
+    Info<< "Creating field dalpha1pdt\n" << endl;
+    volScalarField dalpha1pdt
     (
         IOobject
         (
@@ -285,7 +285,20 @@
             mesh
         ),
         mesh,
-        dimensionedScalar("dpdt", p.dimensions()/dimTime, 0)
+        dimensionedScalar("dalpha1pdt", p.dimensions()/dimTime, 0)
+    );
+
+    Info<< "Creating field dalpha2pdt\n" << endl;
+    volScalarField dalpha2pdt
+    (
+        IOobject
+        (
+            "dpdt",
+            runTime.timeName(),
+            mesh
+        ),
+        mesh,
+        dimensionedScalar("dalpha2pdt", p.dimensions()/dimTime, 0)
     );
 
     Info<< "Creating field kinetic energy K\n" << endl;

applications/solvers/multiphase/compressibleTwoPhaseEulerFoam/pEqn.H

@@ -192,6 +192,11 @@
 
     if (thermo1.dpdt())
     {
-        dpdt = fvc::ddt(p);
+        dalpha1pdt = fvc::ddt(alpha1, p);
+    }
+
+    if (thermo2.dpdt())
+    {
+        dalpha2pdt = fvc::ddt(alpha2, p);
     }
 }