twoPhaseEulerFoam: Added experimental face-based momentum equation formulation

This formulation provides C-grid like pressure-flux staggering on an
unstructured mesh which is hugely beneficial for Euler-Euler multiphase
equations as it allows for all forces to be treated in a consistent
manner on the cell-faces which provides better balance, stability and
accuracy.  However, to achieve face-force consistency the momentum
transport terms must be interpolated to the faces reducing accuracy of
this part of the system but this is offset by the increase in accuracy
of the force-balance.

Currently it is not clear if this face-based momentum equation
formulation is preferable for all Euler-Euler simulations so I have
included it on a switch to allow evaluation and comparison with the
previous cell-based formulation.  To try the new algorithm simply switch
it on, e.g.:

PIMPLE
{
    nOuterCorrectors 3;
    nCorrectors      1;
    nNonOrthogonalCorrectors 0;
    faceMomentum     yes;
}

It is proving particularly good for bubbly flows, eliminating the
staggering patterns often seen in the air velocity field with the
previous algorithm, removing other spurious numerical artifacts in the
velocity fields and improving stability and allowing larger time-steps
For particle-gas flows the advantage is noticeable but not nearly as
pronounced as in the bubbly flow cases.

Please test the new algorithm on your cases and provide feedback.

Henry G. Weller
CFD Direct

applications/solvers/multiphase/twoPhaseEulerFoam/EEqns.H

@@ -5,7 +5,7 @@
     volScalarField Cpv1("Cpv1", thermo1.Cpv());
     volScalarField Cpv2("Cpv2", thermo2.Cpv());
 
-    volScalarField heatTransferCoeff(fluid.heatTransferCoeff());
+    volScalarField Kh(fluid.Kh());
 
     fvScalarMatrix he1Eqn
     (
@@ -32,9 +32,9 @@
 
     he1Eqn -=
     (
-        heatTransferCoeff*(thermo2.T() - thermo1.T())
+        Kh*(thermo2.T() - thermo1.T())
-      + heatTransferCoeff*he1/Cpv1
+      + Kh*he1/Cpv1
-      - fvm::Sp(heatTransferCoeff/Cpv1, he1)
+      - fvm::Sp(Kh/Cpv1, he1)
       + fvOptions(alpha1, rho1, he1)
     );
 
@@ -63,9 +63,9 @@
 
     he2Eqn -=
     (
-        heatTransferCoeff*(thermo1.T() - thermo2.T())
+        Kh*(thermo1.T() - thermo2.T())
-      + heatTransferCoeff*he2/Cpv2
+      + Kh*he2/Cpv2
-      - fvm::Sp(heatTransferCoeff/Cpv2, he2)
+      - fvm::Sp(Kh/Cpv2, he2)
       + fvOptions(alpha2, rho2, he2)
     );
 

applications/solvers/multiphase/twoPhaseEulerFoam/createFields.H

@@ -70,52 +70,6 @@
         fluid.U()
     );
 
-    Info<< "Calculating field DDtU1 and DDtU2\n" << endl;
-
-    volVectorField DDtU1
-    (
-        "DDtU1",
-        fvc::ddt(U1)
-      + fvc::div(phi1, U1)
-      - fvc::div(phi1)*U1
-    );
-
-    volVectorField DDtU2
-    (
-        "DDtU2",
-        fvc::ddt(U2)
-      + fvc::div(phi2, U2)
-      - fvc::div(phi2)*U2
-    );
-
-    volScalarField rAU1
-    (
-        IOobject
-        (
-            IOobject::groupName("rAU", phase1.name()),
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh,
-        dimensionedScalar("zero", dimensionSet(-1, 3, 1, 0, 0), 0.0)
-    );
-
-    volScalarField rAU2
-    (
-        IOobject
-        (
-            IOobject::groupName("rAU", phase2.name()),
-            runTime.timeName(),
-            mesh,
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh,
-        dimensionedScalar("zero", dimensionSet(-1, 3, 1, 0, 0), 0.0)
-    );
-
     label pRefCell = 0;
     scalar pRefValue = 0.0;
     setRefCell
@@ -144,3 +98,60 @@
     Info<< "Creating field kinetic energy K\n" << endl;
     volScalarField K1(IOobject::groupName("K", phase1.name()), 0.5*magSqr(U1));
     volScalarField K2(IOobject::groupName("K", phase2.name()), 0.5*magSqr(U2));
+
+
+    volScalarField rAU1
+    (
+        IOobject
+        (
+            IOobject::groupName("rAU", phase1.name()),
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("zero", dimensionSet(-1, 3, 1, 0, 0), 0)
+    );
+
+    volScalarField rAU2
+    (
+        IOobject
+        (
+            IOobject::groupName("rAU", phase2.name()),
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("zero", dimensionSet(-1, 3, 1, 0, 0), 0)
+    );
+
+    surfaceScalarField rAU1f
+    (
+        IOobject
+        (
+            IOobject::groupName("rAUf", phase1.name()),
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("zero", dimensionSet(-1, 3, 1, 0, 0), 0)
+    );
+
+    surfaceScalarField rAU2f
+    (
+        IOobject
+        (
+            IOobject::groupName("rAUf", phase2.name()),
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("zero", dimensionSet(-1, 3, 1, 0, 0), 0)
+    );

applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/dragModel/dragModel.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -26,6 +26,7 @@ License
 #include "dragModel.H"
 #include "phasePair.H"
 #include "swarmCorrection.H"
+#include "surfaceInterpolate.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -102,12 +103,11 @@ Foam::dragModel::~dragModel()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-Foam::tmp<Foam::volScalarField> Foam::dragModel::K() const
+Foam::tmp<Foam::volScalarField> Foam::dragModel::Ki() const
 {
     return
         0.75
        *CdRe()
-       *max(pair_.dispersed(), residualAlpha_)
        *swarmCorrection_->Cs()
        *pair_.continuous().rho()
        *pair_.continuous().nu()
@@ -115,6 +115,20 @@ Foam::tmp<Foam::volScalarField> Foam::dragModel::K() const
 }
 
 
+Foam::tmp<Foam::volScalarField> Foam::dragModel::K() const
+{
+    return max(pair_.dispersed(), residualAlpha_)*Ki();
+}
+
+
+Foam::tmp<Foam::surfaceScalarField> Foam::dragModel::Kf() const
+{
+    return
+        max(fvc::interpolate(pair_.dispersed()), residualAlpha_)
+       *fvc::interpolate(Ki());
+}
+
+
 bool Foam::dragModel::writeData(Ostream& os) const
 {
     return os.good();

applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/dragModel/dragModel.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2011-2014 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -130,14 +130,32 @@ public:
 
     // Member Functions
 
+        //- Return the residual phase-fraction
+        //  used to stabilize the phase momentum as the phase-fraction -> 0
+        const dimensionedScalar& residualAlpha() const
+        {
+            return residualAlpha_;
+        }
+
         //- Drag coefficient
         virtual tmp<volScalarField> CdRe() const = 0;
 
-        //- The drag function K used in the momentum equation
+        //- Return the phase-intensive drag coefficient Ki
+        //  used in the momentum equations
+        //    ddt(alpha1*rho1*U1) + ... = ... alphad*K*(U1-U2)
+        //    ddt(alpha2*rho2*U2) + ... = ... alphad*K*(U2-U1)
+        virtual tmp<volScalarField> Ki() const;
+
+        //- Return the drag coefficient K
+        //  used in the momentum equations
         //    ddt(alpha1*rho1*U1) + ... = ... K*(U1-U2)
         //    ddt(alpha2*rho2*U2) + ... = ... K*(U2-U1)
         virtual tmp<volScalarField> K() const;
 
+        //- Return the drag coefficient Kf
+        //  used in the face-momentum equations
+        virtual tmp<surfaceScalarField> Kf() const;
+
         //- Dummy write for regIOobject
         bool writeData(Ostream& os) const;
 };

applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/segregated/segregated.C

@@ -26,6 +26,7 @@ License
 #include "segregated.H"
 #include "phasePair.H"
 #include "fvcGrad.H"
+#include "surfaceInterpolate.H"
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -147,4 +148,10 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::segregated::K() const
 }
 
 
+Foam::tmp<Foam::surfaceScalarField> Foam::dragModels::segregated::Kf() const
+{
+    return fvc::interpolate(K());
+}
+
+
 // ************************************************************************* //

applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/dragModels/segregated/segregated.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2014 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2014-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -99,6 +99,9 @@ public:
 
         //- The drag function used in the momentum equation
         virtual tmp<volScalarField> K() const;
+
+        //- The drag function Kf used in the face-momentum equations
+        virtual tmp<surfaceScalarField> Kf() const;
 };
 
 

applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/liftModels/liftModel/liftModel.C

@@ -26,7 +26,7 @@ License
 #include "liftModel.H"
 #include "phasePair.H"
 #include "fvcCurl.H"
-#include "addToRunTimeSelectionTable.H"
+#include "surfaceInterpolate.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -59,11 +59,10 @@ Foam::liftModel::~liftModel()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-Foam::tmp<Foam::volVectorField> Foam::liftModel::F() const
+Foam::tmp<Foam::volVectorField> Foam::liftModel::Fi() const
 {
     return
         Cl()
-       *pair_.dispersed()
        *pair_.continuous().rho()
        *(
             pair_.Ur() ^ fvc::curl(pair_.continuous().U())
@@ -71,4 +70,20 @@ Foam::tmp<Foam::volVectorField> Foam::liftModel::F() const
 }
 
 
+Foam::tmp<Foam::volVectorField> Foam::liftModel::F() const
+{
+    return pair_.dispersed()*Fi();
+}
+
+
+Foam::tmp<Foam::surfaceScalarField> Foam::liftModel::Ff() const
+{
+    const fvMesh& mesh(this->pair_.phase1().mesh());
+
+    return
+        fvc::interpolate(pair_.dispersed())
+       *(fvc::interpolate(Fi()) & mesh.Sf());
+}
+
+
 // ************************************************************************* //

applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/liftModels/liftModel/liftModel.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2014 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2014-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -112,11 +112,17 @@ public:
 
     // Member Functions
 
-        //- Lift coefficient
+        //- Return lift coefficient
         virtual tmp<volScalarField> Cl() const = 0;
 
-        //- Lift force
+        //- Return phase-intensive lift force
+        virtual tmp<volVectorField> Fi() const;
+
+        //- Return lift force
         virtual tmp<volVectorField> F() const;
+
+        //- Return face lift force
+        virtual tmp<surfaceScalarField> Ff() const;
 };
 
 

applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/liftModels/noLift/noLift.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2014 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2014-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -63,8 +63,7 @@ Foam::tmp<Foam::volScalarField> Foam::liftModels::noLift::Cl() const
 {
     const fvMesh& mesh(this->pair_.phase1().mesh());
 
-    return
+    return tmp<volScalarField>
-        tmp<volScalarField>
     (
         new volScalarField
         (
@@ -72,7 +71,10 @@ Foam::tmp<Foam::volScalarField> Foam::liftModels::noLift::Cl() const
             (
                 "Cl",
                 mesh.time().timeName(),
-                    mesh
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
             ),
             mesh,
             dimensionedScalar("Cl", dimless, 0)
@@ -83,13 +85,48 @@ Foam::tmp<Foam::volScalarField> Foam::liftModels::noLift::Cl() const
 
 Foam::tmp<Foam::volVectorField> Foam::liftModels::noLift::F() const
 {
-    return
+    const fvMesh& mesh(this->pair_.phase1().mesh());
-        Cl()
+
-       *dimensionedVector
+    return tmp<volVectorField>
     (
-            "zero",
+        new volVectorField
-            dimensionSet(1, -2, -2, 0, 0),
+        (
-            vector::zero
+            IOobject
+            (
+                "noLift:F",
+                mesh.time().timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            mesh,
+            dimensionedVector("zero", dimF, vector::zero)
+        )
+    );
+}
+
+
+Foam::tmp<Foam::surfaceScalarField> Foam::liftModels::noLift::Ff() const
+{
+    const fvMesh& mesh(this->pair_.phase1().mesh());
+
+    return tmp<surfaceScalarField>
+    (
+        new surfaceScalarField
+        (
+            IOobject
+            (
+                "noLift:Ff",
+                mesh.time().timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            mesh,
+            dimensionedScalar("zero", dimF*dimArea, 0)
+        )
     );
 }
 

applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/liftModels/noLift/noLift.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2014 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2014-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -81,6 +81,9 @@ public:
 
         //- Lift force
         virtual tmp<volVectorField> F() const;
+
+        //- Lift force on faces
+        virtual tmp<surfaceScalarField> Ff() const;
 };
 
 

applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/virtualMassModels/constantVirtualMassCoefficient/constantVirtualMassCoefficient.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2014 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2014-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -73,16 +73,18 @@ Foam::virtualMassModels::constantVirtualMassCoefficient::Cvm() const
 {
     const fvMesh& mesh(this->pair_.phase1().mesh());
 
-    return
+    return tmp<volScalarField>
-        tmp<volScalarField>
     (
         new volScalarField
         (
             IOobject
             (
-                    "zero",
+                "Cvm",
                 mesh.time().timeName(),
-                    mesh
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
             ),
             mesh,
             Cvm_

applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/virtualMassModels/virtualMassModel/virtualMassModel.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2014 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2014-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,6 +25,7 @@ License
 
 #include "virtualMassModel.H"
 #include "phasePair.H"
+#include "surfaceInterpolate.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -70,9 +71,22 @@ Foam::virtualMassModel::~virtualMassModel()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+Foam::tmp<Foam::volScalarField> Foam::virtualMassModel::Ki() const
+{
+    return Cvm()*pair_.continuous().rho();
+}
+
+
 Foam::tmp<Foam::volScalarField> Foam::virtualMassModel::K() const
 {
-    return Cvm()*pair_.dispersed()*pair_.continuous().rho();
+    return pair_.dispersed()*Ki();
+}
+
+
+Foam::tmp<Foam::surfaceScalarField> Foam::virtualMassModel::Kf() const
+{
+    return
+        fvc::interpolate(pair_.dispersed())*fvc::interpolate(Ki());
 }
 
 

applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/virtualMassModels/virtualMassModel/virtualMassModel.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2014 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2014-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -116,14 +116,25 @@ public:
 
     // Member Functions
 
-        //- Virtual mass coefficient
+        //- Return the virtual mass coefficient
         virtual tmp<volScalarField> Cvm() const = 0;
 
-        //- The virtual mass function K used in the momentum equation
+        //- Return the phase-intensive virtual mass coefficient Ki
+        //  used in the momentum equation
+        //    ddt(alpha1*rho1*U1) + ... = ... alphad*K*(DU1_Dt - DU2_Dt)
+        //    ddt(alpha2*rho2*U2) + ... = ... alphad*K*(DU1_Dt - DU2_Dt)
+        virtual tmp<volScalarField> Ki() const;
+
+        //- Return the virtual mass coefficient K
+        //  used in the momentum equation
         //    ddt(alpha1*rho1*U1) + ... = ... K*(DU1_Dt - DU2_Dt)
         //    ddt(alpha2*rho2*U2) + ... = ... K*(DU1_Dt - DU2_Dt)
         virtual tmp<volScalarField> K() const;
 
+        //- Return the virtual mass coefficient Kf
+        //  used in the face-momentum equations
+        virtual tmp<surfaceScalarField> Kf() const;
+
         // Dummy write for regIOobject
         bool writeData(Ostream& os) const;
 };

applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/Antal/Antal.C

@@ -66,7 +66,7 @@ Foam::wallLubricationModels::Antal::~Antal()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-Foam::tmp<Foam::volVectorField> Foam::wallLubricationModels::Antal::F() const
+Foam::tmp<Foam::volVectorField> Foam::wallLubricationModels::Antal::Fi() const
 {
     volVectorField Ur(pair_.Ur());
 
@@ -78,7 +78,6 @@ Foam::tmp<Foam::volVectorField> Foam::wallLubricationModels::Antal::F() const
             dimensionedScalar("zero", dimless/dimLength, 0),
             Cw1_/pair_.dispersed().d() + Cw2_/yWall()
         )
-       *pair_.dispersed()
        *pair_.continuous().rho()
        *magSqr(Ur - (Ur & n)*n)
        *n;

applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/Antal/Antal.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2014 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2014-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -95,8 +95,8 @@ public:
 
     // Member Functions
 
-        //- Wall lubrication force
+        //- Return phase-intensive wall lubrication force
-        tmp<volVectorField> F() const;
+        tmp<volVectorField> Fi() const;
 };
 
 

applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/Frank/Frank.C

@@ -67,7 +67,7 @@ Foam::wallLubricationModels::Frank::~Frank()
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
-Foam::tmp<Foam::volVectorField> Foam::wallLubricationModels::Frank::F() const
+Foam::tmp<Foam::volVectorField> Foam::wallLubricationModels::Frank::Fi() const
 {
     volVectorField Ur(pair_.Ur());
 
@@ -88,7 +88,6 @@ Foam::tmp<Foam::volVectorField> Foam::wallLubricationModels::Frank::F() const
             dimensionedScalar("zero", dimless/dimLength, 0.0),
             (1.0 - yTilde)/(Cwd_*y*pow(yTilde, p_ - 1.0))
         )
-       *pair_.dispersed()
        *pair_.continuous().rho()
        *magSqr(Ur - (Ur & n)*n)
        *n;

applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/Frank/Frank.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2014 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2014-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -105,8 +105,8 @@ public:
 
     // Member Functions
 
-        //- Wall lubrication force
+        //- Return phase-intensive wall lubrication force
-        tmp<volVectorField> F() const;
+        tmp<volVectorField> Fi() const;
 };
 
 

applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/TomiyamaWallLubrication/TomiyamaWallLubrication.C

@@ -66,7 +66,7 @@ Foam::wallLubricationModels::TomiyamaWallLubrication::~TomiyamaWallLubrication()
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 Foam::tmp<Foam::volVectorField>
-Foam::wallLubricationModels::TomiyamaWallLubrication::F() const
+Foam::wallLubricationModels::TomiyamaWallLubrication::Fi() const
 {
     volVectorField Ur(pair_.Ur());
 
@@ -87,7 +87,6 @@ Foam::wallLubricationModels::TomiyamaWallLubrication::F() const
             1/sqr(y)
           - 1/sqr(D_ - y)
         )
-       *pair_.dispersed()
        *pair_.continuous().rho()
        *magSqr(Ur - (Ur & n)*n)
        *n;

applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/TomiyamaWallLubrication/TomiyamaWallLubrication.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2014 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2014-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -98,8 +98,8 @@ public:
 
     // Member Functions
 
-        //- Wall lubrication force
+        //- Return phase-intensive wall lubrication force
-        tmp<volVectorField> F() const;
+        tmp<volVectorField> Fi() const;
 };
 
 

applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/noWallLubrication/noWallLubrication.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2014 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2014-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -65,28 +65,50 @@ Foam::wallLubricationModels::noWallLubrication::~noWallLubrication()
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 Foam::tmp<Foam::volVectorField>
-Foam::wallLubricationModels::noWallLubrication::F() const
+Foam::wallLubricationModels::noWallLubrication::Fi() const
 {
     const fvMesh& mesh(this->pair_.phase1().mesh());
 
-    return
+    return tmp<volVectorField>
-        tmp<volVectorField>
     (
         new volVectorField
         (
             IOobject
             (
-                    "zero",
+                "noWallLubrication:Fi",
                 mesh.time().timeName(),
-                    mesh
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
             ),
             mesh,
-                dimensionedVector
+            dimensionedVector("zero", dimF, vector::zero)
-                (
-                    "zero",
-                    dimensionSet(1, -2, -2, 0, 0),
-                    vector::zero
         )
+    );
+}
+
+
+Foam::tmp<Foam::volVectorField>
+Foam::wallLubricationModels::noWallLubrication::F() const
+{
+    const fvMesh& mesh(this->pair_.phase1().mesh());
+
+    return tmp<volVectorField>
+    (
+        new volVectorField
+        (
+            IOobject
+            (
+                "noWallLubrication:F",
+                mesh.time().timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            mesh,
+            dimensionedVector("zero", dimF, vector::zero)
         )
     );
 }

applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/noWallLubrication/noWallLubrication.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2014 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2014-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -76,6 +76,9 @@ public:
 
     // Member Functions
 
+        //- Return phase-intensive wall lubrication force
+        tmp<volVectorField> Fi() const;
+
         //- Wall lubrication force
         tmp<volVectorField> F() const;
 };

applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/wallLubricationModel/wallLubricationModel.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2014 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2014-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,8 +25,7 @@ License
 
 #include "wallLubricationModel.H"
 #include "phasePair.H"
-
+#include "surfaceInterpolate.H"
-const Foam::dimensionSet Foam::wallLubricationModel::dimF(1, -2, -2, 0, 0);
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -36,6 +35,8 @@ namespace Foam
     defineRunTimeSelectionTable(wallLubricationModel, dictionary);
 }
 
+const Foam::dimensionSet Foam::wallLubricationModel::dimF(1, -2, -2, 0, 0);
+
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
@@ -56,4 +57,22 @@ Foam::wallLubricationModel::~wallLubricationModel()
 {}
 
 
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+Foam::tmp<Foam::volVectorField> Foam::wallLubricationModel::F() const
+{
+    return pair_.dispersed()*Fi();
+}
+
+
+Foam::tmp<Foam::surfaceScalarField> Foam::wallLubricationModel::Ff() const
+{
+    const fvMesh& mesh(this->pair_.phase1().mesh());
+
+    return
+        fvc::interpolate(pair_.dispersed())
+       *(fvc::interpolate(Fi()) & mesh.Sf());
+}
+
+
 // ************************************************************************* //

applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/wallLubricationModels/wallLubricationModel/wallLubricationModel.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 2014 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2014-2015 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -115,8 +115,14 @@ public:
 
     // Member Functions
 
-        //- Wall lubrication force
+        //- Return phase-intensive wall lubrication force
-        virtual tmp<volVectorField> F() const = 0;
+        virtual tmp<volVectorField> Fi() const = 0;
+
+        //- Return wall lubrication force
+        virtual tmp<volVectorField> F() const;
+
+        //- Return face wall lubrication force
+        virtual tmp<surfaceScalarField> Ff() const;
 };
 
 

applications/solvers/multiphase/twoPhaseEulerFoam/pU/UEqns.H

@@ -5,20 +5,20 @@ mrfZones.correctBoundaryVelocity(U);
 fvVectorMatrix U1Eqn(U1, rho1.dimensions()*U1.dimensions()*dimVol/dimTime);
 fvVectorMatrix U2Eqn(U2, rho2.dimensions()*U2.dimensions()*dimVol/dimTime);
 
-volScalarField dragCoeff(fluid.dragCoeff());
+volScalarField Kd(fluid.Kd());
 
 {
-    volScalarField virtualMassCoeff(fluid.virtualMassCoeff());
+    volScalarField Vm(fluid.Vm());
 
     {
         U1Eqn =
         (
             fvm::ddt(alpha1, rho1, U1) + fvm::div(alphaRhoPhi1, U1)
           - fvm::Sp(contErr1, U1)
-          + mrfZones(alpha1*rho1 + virtualMassCoeff, U1)
+          + mrfZones(alpha1*rho1 + Vm, U1)
           + phase1.turbulence().divDevRhoReff(U1)
          ==
-          - virtualMassCoeff
+          - Vm
            *(
                 fvm::ddt(U1)
               + fvm::div(phi1, U1)
@@ -28,7 +28,7 @@ volScalarField dragCoeff(fluid.dragCoeff());
           + fvOptions(alpha1, rho1, U1)
         );
         U1Eqn.relax();
-        U1Eqn += fvm::Sp(dragCoeff, U1);
+        U1Eqn += fvm::Sp(Kd, U1);
         fvOptions.constrain(U1Eqn);
         U1.correctBoundaryConditions();
         fvOptions.correct(U1);
@@ -39,10 +39,10 @@ volScalarField dragCoeff(fluid.dragCoeff());
         (
             fvm::ddt(alpha2, rho2, U2) + fvm::div(alphaRhoPhi2, U2)
           - fvm::Sp(contErr2, U2)
-          + mrfZones(alpha2*rho2 + virtualMassCoeff, U2)
+          + mrfZones(alpha2*rho2 + Vm, U2)
           + phase2.turbulence().divDevRhoReff(U2)
          ==
-          - virtualMassCoeff
+          - Vm
            *(
                 fvm::ddt(U2)
               + fvm::div(phi2, U2)
@@ -52,7 +52,7 @@ volScalarField dragCoeff(fluid.dragCoeff());
           + fvOptions(alpha2, rho2, U2)
         );
         U2Eqn.relax();
-        U2Eqn += fvm::Sp(dragCoeff, U2);
+        U2Eqn += fvm::Sp(Kd, U2);
         fvOptions.constrain(U2Eqn);
         U2.correctBoundaryConditions();
         fvOptions.correct(U2);

applications/solvers/multiphase/twoPhaseEulerFoam/pU/createDDtU.H

@@ -0,0 +1,17 @@
+    Info<< "Calculating field DDtU1 and DDtU2\n" << endl;
+
+    volVectorField DDtU1
+    (
+        "DDtU1",
+        fvc::ddt(U1)
+      + fvc::div(phi1, U1)
+      - fvc::div(phi1)*U1
+    );
+
+    volVectorField DDtU2
+    (
+        "DDtU2",
+        fvc::ddt(U2)
+      + fvc::div(phi2, U2)
+      - fvc::div(phi2)*U2
+    );

applications/solvers/multiphase/twoPhaseEulerFoam/pU/pEqn.H

@@ -39,35 +39,26 @@ while (pimple.correct())
 
     // Turbulent diffusion, particle-pressure lift and wall-lubrication fluxes
     {
-        volScalarField turbulentDiffusivity(fluid.turbulentDiffusivity());
+        volScalarField D(fluid.D());
-        volVectorField liftForce(fluid.liftForce());
+        volVectorField F(fluid.F());
-        volVectorField wallLubricationForce(fluid.wallLubricationForce());
         surfaceScalarField snGradAlpha1(fvc::snGrad(alpha1)*mesh.magSf());
 
         phiF1 =
         (
             fvc::interpolate
             (
-                rAU1*(turbulentDiffusivity + phase1.turbulence().pPrime())
+                rAU1*(D + phase1.turbulence().pPrime())
             )*snGradAlpha1
-
+          + (fvc::interpolate(rAU1*F) & mesh.Sf())
-          + (
-                fvc::interpolate(rAU1*(wallLubricationForce + liftForce))
-              & mesh.Sf()
-            )
         );
 
         phiF2 =
         (
           - fvc::interpolate
             (
-                rAU2*(turbulentDiffusivity + phase2.turbulence().pPrime())
+                rAU2*(D + phase2.turbulence().pPrime())
             )*snGradAlpha1
-
+          - (fvc::interpolate(rAU2*F) & mesh.Sf())
-          - (
-                fvc::interpolate(rAU2*(wallLubricationForce + liftForce))
-              & mesh.Sf()
-            )
         );
     }
 
@@ -139,8 +130,8 @@ while (pimple.correct())
     );
 
     // Face-drag coefficients
-    surfaceScalarField D1f(fvc::interpolate(rAU1*dragCoeff));
+    surfaceScalarField D1f(fvc::interpolate(rAU1*Kd));
-    surfaceScalarField D2f(fvc::interpolate(rAU2*dragCoeff));
+    surfaceScalarField D2f(fvc::interpolate(rAU2*Kd));
 
     // Construct the mean predicted flux
     // including explicit drag contributions based on absolute fluxes
@@ -306,8 +297,8 @@ while (pimple.correct())
                     HbyA2 + fvc::reconstruct(alpharAU2f*mSfGradp - phiF2)
                 );
 
-                volScalarField D1(rAU1*dragCoeff);
+                volScalarField D1(rAU1*Kd);
-                volScalarField D2(rAU2*dragCoeff);
+                volScalarField D2(rAU2*Kd);
 
                 U = alpha1*(U1s + D1*U2) + alpha2*(U2s + D2*U1);
                 volVectorField Ur(((1 - D2)*U1s - (1 - D1)*U2s)/(1 - D1*D2));

applications/solvers/multiphase/twoPhaseEulerFoam/pUf/DDtU.H

@@ -0,0 +1,9 @@
+    ddtPhi1 =
+    (
+        (phi1 - phi1.oldTime())/runTime.deltaT()
+    );
+
+    ddtPhi2 =
+    (
+        (phi2 - phi2.oldTime())/runTime.deltaT()
+    );

applications/solvers/multiphase/twoPhaseEulerFoam/pUf/UEqns.H

@@ -0,0 +1,50 @@
+mrfZones.correctBoundaryVelocity(U1);
+mrfZones.correctBoundaryVelocity(U2);
+mrfZones.correctBoundaryVelocity(U);
+
+fvVectorMatrix U1Eqn(U1, rho1.dimensions()*U1.dimensions()*dimVol/dimTime);
+fvVectorMatrix U2Eqn(U2, rho2.dimensions()*U2.dimensions()*dimVol/dimTime);
+
+{
+    volScalarField Vm(fluid.Vm());
+
+    fvVectorMatrix UgradU1
+    (
+        fvm::div(phi1, U1) - fvm::Sp(fvc::div(phi1), U1)
+      + mrfZones(U1)
+    );
+
+    fvVectorMatrix UgradU2
+    (
+        fvm::div(phi2, U2) - fvm::Sp(fvc::div(phi2), U2)
+      + mrfZones(U2)
+    );
+
+    {
+        U1Eqn =
+        (
+            fvm::div(alphaRhoPhi1, U1) - fvm::Sp(fvc::div(alphaRhoPhi1), U1)
+          + mrfZones(alpha1*rho1, U1)
+          + phase1.turbulence().divDevRhoReff(U1)
+          + Vm*(UgradU1 - (UgradU2 & U2))
+        );
+        U1Eqn.relax();
+        fvOptions.constrain(U1Eqn);
+        U1.correctBoundaryConditions();
+        fvOptions.correct(U1);
+    }
+
+    {
+        U2Eqn =
+        (
+            fvm::div(alphaRhoPhi2, U2) - fvm::Sp(fvc::div(alphaRhoPhi2), U2)
+          + mrfZones(alpha2*rho2, U2)
+          + phase2.turbulence().divDevRhoReff(U2)
+          + Vm*(UgradU2 - (UgradU1 & U1))
+        );
+        U2Eqn.relax();
+        fvOptions.constrain(U2Eqn);
+        U2.correctBoundaryConditions();
+        fvOptions.correct(U2);
+    }
+}

applications/solvers/multiphase/twoPhaseEulerFoam/pUf/createDDtU.H

@@ -0,0 +1,9 @@
+    surfaceScalarField ddtPhi1
+    (
+        (phi1 - phi1.oldTime())/runTime.deltaT()
+    );
+
+    surfaceScalarField ddtPhi2
+    (
+        (phi2 - phi2.oldTime())/runTime.deltaT()
+    );

applications/solvers/multiphase/twoPhaseEulerFoam/pUf/pEqn.H

@@ -0,0 +1,347 @@
+surfaceScalarField alpha1f("alpha1f", fvc::interpolate(alpha1));
+surfaceScalarField alpha2f("alpha2f", scalar(1) - alpha1f);
+
+surfaceScalarField alphaRho1f0
+(
+    "alphaRho1f0",
+    fvc::interpolate
+    (
+        max(alpha1.oldTime(), fluid.residualAlpha(phase1))
+       *rho1.oldTime()
+    )
+);
+
+surfaceScalarField alphaRho2f0
+(
+    "alphaRho2f0",
+    fvc::interpolate
+    (
+        max(alpha2.oldTime(), fluid.residualAlpha(phase2))
+       *rho2.oldTime()
+    )
+);
+
+// Drag coefficient
+surfaceScalarField Kdf("Kdf", fluid.Kdf());
+
+// Virtual-mass coefficient
+surfaceScalarField Vmf("Vmf", fluid.Vmf());
+
+// Lift and wall-lubrication forces
+surfaceScalarField Ff("Ff", fluid.Ff());
+
+tmp<surfaceScalarField> snGradAlpha1(fvc::snGrad(alpha1)*mesh.magSf());
+tmp<volScalarField> D(fluid.D());
+
+// Turbulent-dispersion force for phase 1
+surfaceScalarField Ftdf1
+(
+    IOobject::groupName("Ftdf", phase1.name()),
+    fvc::interpolate(D() + phase1.turbulence().pPrime())*snGradAlpha1()
+);
+
+// Turbulent-dispersion force for phase 2
+surfaceScalarField Ftdf2
+(
+    IOobject::groupName("Ftdf", phase2.name()),
+    fvc::interpolate(D + phase2.turbulence().pPrime())*snGradAlpha1
+);
+
+
+while (pimple.correct())
+{
+    // Update continuity errors due to temperature changes
+    #include "correctContErrs.H"
+
+    surfaceScalarField rho1f(fvc::interpolate(rho1));
+    surfaceScalarField rho2f(fvc::interpolate(rho2));
+
+    // Correct flux BCs to be consistent with the velocity BCs
+    phi1.boundaryField() ==
+        mrfZones.relative(mesh.Sf().boundaryField() & U1.boundaryField());
+    phi2.boundaryField() ==
+        mrfZones.relative(mesh.Sf().boundaryField() & U2.boundaryField());
+
+    rAU1f =
+    (
+        IOobject::groupName("rAUf", phase1.name()),
+        1.0
+       /(
+           (alphaRho1f0 + Vmf)/runTime.deltaT()
+         + fvc::interpolate(U1Eqn.A())
+         + Kdf
+        )
+    );
+
+    rAU2f =
+    (
+        IOobject::groupName("rAUf", phase2.name()),
+        1.0
+       /(
+           (alphaRho2f0 + Vmf)/runTime.deltaT()
+         + fvc::interpolate(U2Eqn.A())
+         + Kdf
+        )
+    );
+
+    surfaceScalarField rAlphaAU1f
+    (
+        IOobject::groupName("rAlphaAUf", phase1.name()),
+        max(alpha1f, fluid.residualAlpha(phase1))*rAU1f
+    );
+
+    surfaceScalarField rAlphaAU2f
+    (
+        IOobject::groupName("rAlphaAUf", phase2.name()),
+        max(alpha2f, fluid.residualAlpha(phase2))*rAU2f
+    );
+
+
+    volScalarField rho("rho", fluid.rho());
+    surfaceScalarField ghSnGradRho
+    (
+        "ghSnGradRho",
+        ghf*fvc::snGrad(rho)*mesh.magSf()
+    );
+
+    // Add the phase-1 buoyancy force
+    surfaceScalarField phiF1
+    (
+        IOobject::groupName("phiF", phase1.name()),
+        rAlphaAU1f
+       *(
+            ghSnGradRho
+          - alpha2f*(rho1f - rho2f)*(g & mesh.Sf())
+        )
+    );
+
+    // Add the phase-2 buoyancy force
+    surfaceScalarField phiF2
+    (
+        IOobject::groupName("phiF", phase2.name()),
+        rAlphaAU2f
+       *(
+            ghSnGradRho
+          - alpha1f*(rho2f - rho1f)*(g & mesh.Sf())
+        )
+    );
+
+
+    // Phase-1 predicted flux
+    surfaceScalarField phiHbyA1
+    (
+        IOobject::groupName("phiHbyA", phase1.name()),
+        phi1
+    );
+
+    phiHbyA1 =
+        rAU1f
+       *(
+            (alphaRho1f0 + Vmf)
+           *mrfZones.absolute(phi1.oldTime())/runTime.deltaT()
+          + (fvc::interpolate(U1Eqn.H()) & mesh.Sf())
+          + Vmf*ddtPhi2
+          + Kdf*mrfZones.absolute(phi2)
+          - Ff
+          - Ftdf1
+        );
+
+    // Phase-2 predicted flux
+    surfaceScalarField phiHbyA2
+    (
+        IOobject::groupName("phiHbyA", phase2.name()),
+        phi2
+    );
+
+    phiHbyA2 =
+        rAU2f
+       *(
+            (alphaRho2f0 + Vmf)
+           *mrfZones.absolute(phi2.oldTime())/runTime.deltaT()
+          + (fvc::interpolate(U2Eqn.H()) & mesh.Sf())
+          + Vmf*ddtPhi1
+          + Kdf*mrfZones.absolute(phi1)
+          + Ff
+          + Ftdf2
+       );
+
+
+    surfaceScalarField phiHbyA
+    (
+        "phiHbyA",
+        alpha1f*(phiHbyA1 - phiF1) + alpha2f*(phiHbyA2 - phiF2)
+    );
+    mrfZones.makeRelative(phiHbyA);
+
+    phiHbyA1 -= phiF1;
+    phiHbyA2 -= phiF2;
+
+    surfaceScalarField rAUf
+    (
+        "rAUf",
+        mag(alpha1f*rAlphaAU1f + alpha2f*rAlphaAU2f)
+    );
+
+    // Update the fixedFluxPressure BCs to ensure flux consistency
+    setSnGrad<fixedFluxPressureFvPatchScalarField>
+    (
+        p_rgh.boundaryField(),
+        (
+            phiHbyA.boundaryField()
+          - (
+                alpha1f.boundaryField()*phi1.boundaryField()
+              + alpha2f.boundaryField()*phi2.boundaryField()
+            )
+        )/(mesh.magSf().boundaryField()*rAUf.boundaryField())
+    );
+
+    tmp<fvScalarMatrix> pEqnComp1;
+    tmp<fvScalarMatrix> pEqnComp2;
+
+    if (pimple.transonic())
+    {
+        surfaceScalarField phid1
+        (
+            IOobject::groupName("phid", phase1.name()),
+            fvc::interpolate(psi1)*phi1
+        );
+        surfaceScalarField phid2
+        (
+            IOobject::groupName("phid", phase2.name()),
+            fvc::interpolate(psi2)*phi2
+        );
+
+        pEqnComp1 =
+            (
+                contErr1
+              - fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1)
+            )/rho1
+          + (alpha1/rho1)*correction
+            (
+                psi1*fvm::ddt(p_rgh)
+              + fvm::div(phid1, p_rgh) - fvm::Sp(fvc::div(phid1), p_rgh)
+            );
+        deleteDemandDrivenData(pEqnComp1().faceFluxCorrectionPtr());
+        pEqnComp1().relax();
+
+        pEqnComp2 =
+            (
+                contErr2
+              - fvc::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), rho2)
+            )/rho2
+          + (alpha2/rho2)*correction
+            (
+                psi2*fvm::ddt(p_rgh)
+              + fvm::div(phid2, p_rgh) - fvm::Sp(fvc::div(phid2), p_rgh)
+            );
+        deleteDemandDrivenData(pEqnComp2().faceFluxCorrectionPtr());
+        pEqnComp2().relax();
+    }
+    else
+    {
+        pEqnComp1 =
+            (
+                contErr1
+              - fvc::Sp(fvc::ddt(alpha1) + fvc::div(alphaPhi1), rho1)
+            )/rho1
+          + (alpha1*psi1/rho1)*correction(fvm::ddt(p_rgh));
+
+        pEqnComp2 =
+            (
+                contErr2
+              - fvc::Sp(fvc::ddt(alpha2) + fvc::div(alphaPhi2), rho2)
+            )/rho2
+          + (alpha2*psi2/rho2)*correction(fvm::ddt(p_rgh));
+    }
+
+    // Cache p prior to solve for density update
+    volScalarField p_rgh_0("p_rgh_0", p_rgh);
+
+    while (pimple.correctNonOrthogonal())
+    {
+        fvScalarMatrix pEqnIncomp
+        (
+            fvc::div(phiHbyA)
+          - fvm::laplacian(rAUf, p_rgh)
+        );
+
+        solve
+        (
+            pEqnComp1() + pEqnComp2() + pEqnIncomp,
+            mesh.solver(p_rgh.select(pimple.finalInnerIter()))
+        );
+
+        if (pimple.finalNonOrthogonalIter())
+        {
+            surfaceScalarField mSfGradp("mSfGradp", pEqnIncomp.flux()/rAUf);
+
+            phi = phiHbyA + pEqnIncomp.flux();
+
+            surfaceScalarField phi1s
+            (
+                phiHbyA1
+              + rAlphaAU1f*mSfGradp
+              - rAU1f*Kdf*mrfZones.absolute(phi2)
+            );
+
+            surfaceScalarField phi2s
+            (
+                phiHbyA2
+              + rAlphaAU2f*mSfGradp
+              - rAU2f*Kdf*mrfZones.absolute(phi1)
+            );
+
+            surfaceScalarField phir
+            (
+                ((phi2s + rAU2f*Kdf*phi1s) - (phi1s + rAU1f*Kdf*phi2s))
+               /(1.0 - rAU1f*rAU2f*sqr(Kdf))
+            );
+
+            phi1 = phi - alpha2f*phir;
+            phi2 = phi + alpha1f*phir;
+
+            U1 = fvc::reconstruct(mrfZones.absolute(phi1));
+            U1.correctBoundaryConditions();
+            fvOptions.correct(U1);
+
+            U2 = fvc::reconstruct(mrfZones.absolute(phi2));
+            U2.correctBoundaryConditions();
+            fvOptions.correct(U2);
+
+            U = fluid.U();
+
+            fluid.dgdt() =
+            (
+                alpha1*(pEqnComp2 & p_rgh)
+              - alpha2*(pEqnComp1 & p_rgh)
+            );
+        }
+    }
+
+    // Update and limit the static pressure
+    p = max(p_rgh + rho*gh, pMin);
+
+    // Limit p_rgh
+    p_rgh = p - rho*gh;
+
+    // Update densities from change in p_rgh
+    rho1 += psi1*(p_rgh - p_rgh_0);
+    rho2 += psi2*(p_rgh - p_rgh_0);
+
+    // Correct p_rgh for consistency with p and the updated densities
+    rho = fluid.rho();
+    p_rgh = p - rho*gh;
+    p_rgh.correctBoundaryConditions();
+}
+
+// Update the phase kinetic energies
+K1 = 0.5*magSqr(U1);
+K2 = 0.5*magSqr(U2);
+
+// Update the pressure time-derivative if required
+if (thermo1.dpdt() || thermo2.dpdt())
+{
+    dpdt = fvc::ddt(p);
+}
+
+#include "pUf/DDtU.H"

applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseEulerFoam.C

@@ -50,6 +50,11 @@ int main(int argc, char *argv[])
 
     pimpleControl pimple(mesh);
 
+    Switch faceMomentum
+    (
+        pimple.dict().lookupOrDefault<Switch>("faceMomentum", false)
+    );
+
     #include "createFields.H"
     #include "createMRFZones.H"
     #include "createFvOptions.H"
@@ -58,6 +63,9 @@ int main(int argc, char *argv[])
     #include "CourantNos.H"
     #include "setInitialDeltaT.H"
 
+    #include "pUf/createDDtU.H"
+    #include "pU/createDDtU.H"
+
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
@@ -78,10 +86,21 @@ int main(int argc, char *argv[])
             fluid.correct();
 
             #include "contErrs.H"
-            #include "UEqns.H"
             #include "EEqns.H"
-            #include "pEqn.H"
+
-            #include "DDtU.H"
+            if (faceMomentum)
+            {
+                Info<< "Constructing face momentum equations" << endl;
+                #include "pUf/UEqns.H"
+                #include "pUf/pEqn.H"
+            }
+            else
+            {
+                Info<< "Constructing momentum equations" << endl;
+                #include "pU/UEqns.H"
+                #include "pU/pEqn.H"
+                #include "pU/DDtU.H"
+            }
 
             if (pimple.turbCorr())
             {

applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/BlendedInterfacialModel/BlendedInterfacialModel.C

@@ -25,14 +25,15 @@ License
 
 #include "BlendedInterfacialModel.H"
 #include "fixedValueFvsPatchFields.H"
+#include "surfaceInterpolate.H"
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
 template<class modelType>
-template<class Type>
+template<class GeometricField>
 void Foam::BlendedInterfacialModel<modelType>::correctFixedFluxBCs
 (
-    GeometricField<Type, fvPatchField, volMesh>& field
+    GeometricField& field
 ) const
 {
     forAll(pair_.phase1().phi().boundaryField(), patchI)
@@ -45,7 +46,8 @@ void Foam::BlendedInterfacialModel<modelType>::correctFixedFluxBCs
             )
         )
         {
-            field.boundaryField()[patchI] = pTraits<Type>::zero;
+            field.boundaryField()[patchI]
+              = pTraits<typename GeometricField::value_type>::zero;
         }
     }
 }
@@ -139,9 +141,12 @@ Foam::BlendedInterfacialModel<modelType>::K() const
         (
             IOobject
             (
-                modelType::typeName + "Coeff",
+                modelType::typeName + ":K",
                 pair_.phase1().mesh().time().timeName(),
-                pair_.phase1().mesh()
+                pair_.phase1().mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
             ),
             pair_.phase1().mesh(),
             dimensionedScalar("zero", modelType::dimK, 0)
@@ -176,6 +181,74 @@ Foam::BlendedInterfacialModel<modelType>::K() const
 }
 
 
+template<class modelType>
+Foam::tmp<Foam::surfaceScalarField>
+Foam::BlendedInterfacialModel<modelType>::Kf() const
+{
+    tmp<surfaceScalarField> f1, f2;
+
+    if (model_.valid() || model1In2_.valid())
+    {
+        f1 = fvc::interpolate
+        (
+            blending_.f1(pair1In2_.dispersed(), pair2In1_.dispersed())
+        );
+    }
+
+    if (model_.valid() || model2In1_.valid())
+    {
+        f2 = fvc::interpolate
+        (
+            blending_.f2(pair1In2_.dispersed(), pair2In1_.dispersed())
+        );
+    }
+
+    tmp<surfaceScalarField> x
+    (
+        new surfaceScalarField
+        (
+            IOobject
+            (
+                modelType::typeName + ":Kf",
+                pair_.phase1().mesh().time().timeName(),
+                pair_.phase1().mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            pair_.phase1().mesh(),
+            dimensionedScalar("zero", modelType::dimK, 0)
+        )
+    );
+
+    if (model_.valid())
+    {
+        x() += model_->Kf()*(f1() - f2());
+    }
+
+    if (model1In2_.valid())
+    {
+        x() += model1In2_->Kf()*(1 - f1);
+    }
+
+    if (model2In1_.valid())
+    {
+        x() += model2In1_->Kf()*f2;
+    }
+
+    if
+    (
+        correctFixedFluxBCs_
+     && (model_.valid() || model1In2_.valid() || model2In1_.valid())
+    )
+    {
+        correctFixedFluxBCs(x());
+    }
+
+    return x;
+}
+
+
 template<class modelType>
 template<class Type>
 Foam::tmp<Foam::GeometricField<Type, Foam::fvPatchField, Foam::volMesh> >
@@ -199,9 +272,12 @@ Foam::BlendedInterfacialModel<modelType>::F() const
         (
             IOobject
             (
-                modelType::typeName + "Coeff",
+                modelType::typeName + ":F",
                 pair_.phase1().mesh().time().timeName(),
-                pair_.phase1().mesh()
+                pair_.phase1().mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
             ),
             pair_.phase1().mesh(),
             dimensioned<Type>("zero", modelType::dimF, pTraits<Type>::zero)
@@ -236,6 +312,74 @@ Foam::BlendedInterfacialModel<modelType>::F() const
 }
 
 
+template<class modelType>
+Foam::tmp<Foam::surfaceScalarField>
+Foam::BlendedInterfacialModel<modelType>::Ff() const
+{
+    tmp<surfaceScalarField> f1, f2;
+
+    if (model_.valid() || model1In2_.valid())
+    {
+        f1 = fvc::interpolate
+        (
+            blending_.f1(pair1In2_.dispersed(), pair2In1_.dispersed())
+        );
+    }
+
+    if (model_.valid() || model2In1_.valid())
+    {
+        f2 = fvc::interpolate
+        (
+            blending_.f2(pair1In2_.dispersed(), pair2In1_.dispersed())
+        );
+    }
+
+    tmp<surfaceScalarField> x
+    (
+        new surfaceScalarField
+        (
+            IOobject
+            (
+                modelType::typeName + ":Ff",
+                pair_.phase1().mesh().time().timeName(),
+                pair_.phase1().mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            pair_.phase1().mesh(),
+            dimensionedScalar("zero", modelType::dimF*dimArea, 0)
+        )
+    );
+
+    if (model_.valid())
+    {
+        x() += model_->Ff()*(f1() - f2());
+    }
+
+    if (model1In2_.valid())
+    {
+        x() += model1In2_->Ff()*(1 - f1);
+    }
+
+    if (model2In1_.valid())
+    {
+        x() -= model2In1_->Ff()*f2; // note : subtraction
+    }
+
+    if
+    (
+        correctFixedFluxBCs_
+     && (model_.valid() || model1In2_.valid() || model2In1_.valid())
+    )
+    {
+        correctFixedFluxBCs(x());
+    }
+
+    return x;
+}
+
+
 template<class modelType>
 Foam::tmp<Foam::volScalarField>
 Foam::BlendedInterfacialModel<modelType>::D() const
@@ -258,9 +402,12 @@ Foam::BlendedInterfacialModel<modelType>::D() const
         (
             IOobject
             (
-                modelType::typeName + "Coeff",
+                modelType::typeName + ":D",
                 pair_.phase1().mesh().time().timeName(),
-                pair_.phase1().mesh()
+                pair_.phase1().mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
             ),
             pair_.phase1().mesh(),
             dimensionedScalar("zero", modelType::dimD, 0)
@@ -295,6 +442,19 @@ Foam::BlendedInterfacialModel<modelType>::D() const
 }
 
 
+template<class modelType>
+bool Foam::BlendedInterfacialModel<modelType>::hasModel
+(
+    const class phaseModel& phase
+) const
+{
+    return
+       &phase == &(pair_.phase1())
+      ? model1In2_.valid()
+      : model2In1_.valid();
+}
+
+
 template<class modelType>
 const modelType& Foam::BlendedInterfacialModel<modelType>::phaseModel
 (

applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/BlendedInterfacialModel/BlendedInterfacialModel.H

@@ -88,11 +88,8 @@ class BlendedInterfacialModel
         void operator=(const BlendedInterfacialModel<modelType>&);
 
         //- Correct coeff/value on fixed flux boundary conditions
-        template<class Type>
+        template<class GeometricField>
-        void correctFixedFluxBCs
+        void correctFixedFluxBCs(GeometricField& field) const;
-        (
-            GeometricField<Type, fvPatchField, volMesh>& field
-        ) const;
 
 
 public:
@@ -117,18 +114,27 @@ public:
 
     // Member Functions
 
+        //- Return true if a model is specified for the supplied phase
+        bool hasModel(const phaseModel& phase) const;
+
+        //- Return the model for the supplied phase
+        const modelType& phaseModel(const phaseModel& phase) const;
+
         //- Return the blended force coefficient
         tmp<volScalarField> K() const;
 
+        //- Return the face blended force coefficient
+        tmp<surfaceScalarField> Kf() const;
+
         //- Return the blended force
         template<class Type>
         tmp<GeometricField<Type, fvPatchField, volMesh> > F() const;
 
+        //- Return the face blended force
+        tmp<surfaceScalarField> Ff() const;
+
         //- Return the blended diffusivity
         tmp<volScalarField> D() const;
-
-        //- Return the model for the supplied phase
-        const modelType& phaseModel(const phaseModel& phase) const;
 };
 
 

applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.C

@@ -26,7 +26,6 @@ License
 #include "twoPhaseSystem.H"
 #include "PhaseCompressibleTurbulenceModel.H"
 #include "BlendedInterfacialModel.H"
-#include "dragModel.H"
 #include "virtualMassModel.H"
 #include "heatTransferModel.H"
 #include "liftModel.H"
@@ -48,6 +47,15 @@ License
 #include "blendingMethod.H"
 #include "HashPtrTable.H"
 
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+Foam::dimensionedScalar Foam::twoPhaseSystem::zeroResidualAlpha_
+(
+    "zeroResidualAlpha", dimless, 0
+);
+
+
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 Foam::twoPhaseSystem::twoPhaseSystem
@@ -299,46 +307,49 @@ Foam::tmp<Foam::surfaceScalarField> Foam::twoPhaseSystem::calcPhi() const
 }
 
 
-Foam::tmp<Foam::volScalarField> Foam::twoPhaseSystem::dragCoeff() const
+Foam::tmp<Foam::volScalarField> Foam::twoPhaseSystem::Kd() const
 {
     return drag_->K();
 }
 
 
-Foam::tmp<Foam::volScalarField> Foam::twoPhaseSystem::virtualMassCoeff() const
+Foam::tmp<Foam::surfaceScalarField> Foam::twoPhaseSystem::Kdf() const
+{
+    return drag_->Kf();
+}
+
+
+Foam::tmp<Foam::volScalarField> Foam::twoPhaseSystem::Vm() const
 {
     return virtualMass_->K();
 }
 
 
-Foam::tmp<Foam::volScalarField> Foam::twoPhaseSystem::heatTransferCoeff() const
+Foam::tmp<Foam::surfaceScalarField> Foam::twoPhaseSystem::Vmf() const
+{
+    return virtualMass_->Kf();
+}
+
+
+Foam::tmp<Foam::volScalarField> Foam::twoPhaseSystem::Kh() const
 {
     return heatTransfer_->K();
 }
 
 
-Foam::tmp<Foam::volVectorField> Foam::twoPhaseSystem::liftForce() const
+Foam::tmp<Foam::volVectorField> Foam::twoPhaseSystem::F() const
 {
-    return lift_->F<vector>();
+    return lift_->F<vector>() + wallLubrication_->F<vector>();
 }
 
 
-Foam::tmp<Foam::volVectorField>
+Foam::tmp<Foam::surfaceScalarField> Foam::twoPhaseSystem::Ff() const
-Foam::twoPhaseSystem::wallLubricationForce() const
 {
-    return wallLubrication_->F<vector>();
+    return lift_->Ff() + wallLubrication_->Ff();
 }
 
 
-Foam::tmp<Foam::volVectorField>
+Foam::tmp<Foam::volScalarField> Foam::twoPhaseSystem::D() const
-Foam::twoPhaseSystem::turbulentDispersionForce() const
-{
-    return turbulentDispersion_->F<vector>();
-}
-
-
-Foam::tmp<Foam::volScalarField>
-Foam::twoPhaseSystem::turbulentDiffusivity() const
 {
     return turbulentDispersion_->D();
 }
@@ -354,6 +365,13 @@ void Foam::twoPhaseSystem::solve()
     const surfaceScalarField& phi1 = phase1_.phi();
     const surfaceScalarField& phi2 = phase2_.phi();
 
+    Switch faceMomentum
+    (
+        //pimple.dict().lookupOrDefault<Switch>("faceMomentum", false)
+        mesh_.solutionDict().subDict("PIMPLE")
+           .lookupOrDefault<Switch>("faceMomentum", false)
+    );
+
     const dictionary& alphaControls = mesh_.solverDict
     (
         alpha1.name()
@@ -380,6 +398,27 @@ void Foam::twoPhaseSystem::solve()
     tmp<surfaceScalarField> pPrimeByA;
 
     if (implicitPhasePressure)
+    {
+        if (faceMomentum)
+        {
+            const surfaceScalarField& rAU1f =
+                mesh_.lookupObject<surfaceScalarField>
+                (
+                    IOobject::groupName("rAUf", phase1_.name())
+                );
+            const surfaceScalarField& rAU2f =
+                mesh_.lookupObject<surfaceScalarField>
+                (
+                    IOobject::groupName("rAUf", phase2_.name())
+                );
+
+            volScalarField D(this->D());
+
+            pPrimeByA =
+                rAU1f*fvc::interpolate(D + phase1_.turbulence().pPrime())
+              + rAU2f*fvc::interpolate(D + phase2_.turbulence().pPrime());
+        }
+        else
         {
             const volScalarField& rAU1 = mesh_.lookupObject<volScalarField>
             (
@@ -393,6 +432,7 @@ void Foam::twoPhaseSystem::solve()
             pPrimeByA =
                 fvc::interpolate(rAU1*phase1_.turbulence().pPrime())
               + fvc::interpolate(rAU2*phase2_.turbulence().pPrime());
+        }
 
         surfaceScalarField phiP
         (
@@ -596,8 +636,21 @@ bool Foam::twoPhaseSystem::read()
 }
 
 
-const Foam::dragModel&
+const Foam::dimensionedScalar&
-Foam::twoPhaseSystem::drag(const phaseModel& phase) const
+Foam::twoPhaseSystem::residualAlpha(const phaseModel& phase) const
+{
+    if (drag_->hasModel(phase))
+    {
+        return drag_->phaseModel(phase).residualAlpha();
+    }
+    else
+    {
+        return zeroResidualAlpha_;
+    }
+}
+
+
+const Foam::dragModel& Foam::twoPhaseSystem::drag(const phaseModel& phase) const
 {
     return drag_->phaseModel(phase);
 }

applications/solvers/multiphase/twoPhaseEulerFoam/twoPhaseSystem/twoPhaseSystem.H

@@ -40,13 +40,13 @@ SourceFiles
 #include "orderedPhasePair.H"
 #include "volFields.H"
 #include "surfaceFields.H"
+#include "dragModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
-class dragModel;
 class virtualMassModel;
 class heatTransferModel;
 class liftModel;
@@ -113,6 +113,9 @@ class twoPhaseSystem
         autoPtr<BlendedInterfacialModel<turbulentDispersionModel> >
             turbulentDispersion_;
 
+        //-
+        static dimensionedScalar zeroResidualAlpha_;
+
 
     // Private member functions
 
@@ -141,26 +144,29 @@ public:
         tmp<volVectorField> U() const;
 
         //- Return the drag coefficient
-        tmp<volScalarField> dragCoeff() const;
+        tmp<volScalarField> Kd() const;
+
+        //- Return the face drag coefficient
+        tmp<surfaceScalarField> Kdf() const;
 
         //- Return the virtual mass coefficient
-        tmp<volScalarField> virtualMassCoeff() const;
+        tmp<volScalarField> Vm() const;
+
+        //- Return the face virtual mass coefficient
+        tmp<surfaceScalarField> Vmf() const;
 
         //- Return the heat transfer coefficient
-        tmp<volScalarField> heatTransferCoeff() const;
+        tmp<volScalarField> Kh() const;
 
-        //- Return the lift force
+        //- Return the combined force (lift + wall-lubrication)
-        tmp<volVectorField> liftForce() const;
+        tmp<volVectorField> F() const;
 
-        //- Return the wall lubrication force
+        //- Return the combined face-force (lift + wall-lubrication)
-        tmp<volVectorField> wallLubricationForce() const;
+        tmp<surfaceScalarField> Ff() const;
 
         //- Return the turbulent diffusivity
         //  Multiplies the phase-fraction gradient
-        tmp<volScalarField> turbulentDiffusivity() const;
+        tmp<volScalarField> D() const;
-
-        //- Return the turbulent dispersion force
-        tmp<volVectorField> turbulentDispersionForce() const;
 
         //- Solve for the two-phase-fractions
         void solve();
@@ -176,10 +182,17 @@ public:
 
         // Access
 
-            //- Return the drag model for the supplied phase
+            //- Return the residual phase-fraction for given phase
+            //  Used to stabilize the phase momentum as the phase-fraction -> 0
+            const dimensionedScalar& residualAlpha
+            (
+                const phaseModel& phase
+            ) const;
+
+            //- Return the drag model for the given phase
             const dragModel& drag(const phaseModel& phase) const;
 
-            //- Return the virtual mass model for the supplied phase
+            //- Return the virtual mass model for the given phase
             const virtualMassModel& virtualMass(const phaseModel& phase) const;
 
             //- Return the surface tension coefficient