Merge branch 'master' of /home/dm4/OpenFOAM/OpenFOAM-dev

applications/solvers/compressible/rhoCentralFoam/rhoCentralDyMFoam/rhoCentralDyMFoam.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2012 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2013 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -165,7 +165,6 @@ int main(int argc, char *argv[])
         mesh.movePoints(motionPtr->newPoints());
 
         phi = aphiv_pos*rho_pos + aphiv_neg*rho_neg;
-        Info<< phi.boundaryField()[0] << endl;
 
         surfaceVectorField phiUp
         (

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);
     }
 }