twoPhaseEulerFoam: Interpolate lift, wall-lubrication and turbulent dispersion forces

Reduces or eliminates staggering patterns due to cell-force imbalances
Resolves bug-report http://www.openfoam.org/mantisbt/view.php?id=1363

applications/solvers/multiphase/twoPhaseEulerFoam/UEqns.H

@@ -9,9 +9,6 @@ volScalarField dragCoeff(fluid.dragCoeff());
 
 {
     volScalarField virtualMassCoeff(fluid.virtualMassCoeff());
-    volVectorField liftForce(fluid.liftForce());
-    volVectorField wallLubricationForce(fluid.wallLubricationForce());
-    volVectorField turbulentDispersionForce(fluid.turbulentDispersionForce());
 
     {
         U1Eqn =
@@ -21,9 +18,6 @@ volScalarField dragCoeff(fluid.dragCoeff());
           + mrfZones(alpha1*rho1 + virtualMassCoeff, U1)
           + phase1.turbulence().divDevRhoReff(U1)
          ==
-          - liftForce
-          - wallLubricationForce
-          - turbulentDispersionForce
           - virtualMassCoeff
            *(
                 fvm::ddt(U1)
@@ -46,9 +40,6 @@ volScalarField dragCoeff(fluid.dragCoeff());
           + mrfZones(alpha2*rho2 + virtualMassCoeff, U2)
           + phase2.turbulence().divDevRhoReff(U2)
          ==
-            liftForce
-          + wallLubricationForce
-          + turbulentDispersionForce
           - virtualMassCoeff
            *(
                 fvm::ddt(U2)

applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/Burns/Burns.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
@@ -70,7 +70,7 @@ Foam::turbulentDispersionModels::Burns::~Burns()
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 Foam::tmp<Foam::volScalarField>
-Foam::turbulentDispersionModels::Burns::Fprime() const
+Foam::turbulentDispersionModels::Burns::D() const
 {
     const fvMesh& mesh(pair_.phase1().mesh());
     const dragModel&

applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/Burns/Burns.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
@@ -103,9 +103,9 @@ public:
 
     // Member Functions
 
-        //- Turbulent dispersion force coefficient
+        //- Turbulent diffusivity
         //  multiplying the gradient of the phase-fraction
-        virtual tmp<volScalarField> Fprime() const;
+        virtual tmp<volScalarField> D() const;
 };
 
 

applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/Gosman/Gosman.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
@@ -69,7 +69,7 @@ Foam::turbulentDispersionModels::Gosman::~Gosman()
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 Foam::tmp<Foam::volScalarField>
-Foam::turbulentDispersionModels::Gosman::Fprime() const
+Foam::turbulentDispersionModels::Gosman::D() const
 {
     const fvMesh& mesh(pair_.phase1().mesh());
     const dragModel&

applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/Gosman/Gosman.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
@@ -92,9 +92,9 @@ public:
 
     // Member Functions
 
-        //- Turbulent dispersion force coefficient
+        //- Turbulent diffusivity
         //  multiplying the gradient of the phase-fraction
-        virtual tmp<volScalarField> Fprime() const;
+        virtual tmp<volScalarField> D() const;
 };
 
 

applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/LopezDeBertodano/LopezDeBertodano.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
@@ -67,7 +67,7 @@ Foam::turbulentDispersionModels::LopezDeBertodano::~LopezDeBertodano()
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 Foam::tmp<Foam::volScalarField>
-Foam::turbulentDispersionModels::LopezDeBertodano::Fprime() const
+Foam::turbulentDispersionModels::LopezDeBertodano::D() const
 {
     return
         Ctd_
@@ -75,4 +75,5 @@ Foam::turbulentDispersionModels::LopezDeBertodano::Fprime() const
        *pair_.continuous().turbulence().k();
 }
 
+
 // ************************************************************************* //

applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/LopezDeBertodano/LopezDeBertodano.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,9 +98,9 @@ public:
 
     // Member Functions
 
-        //- Turbulent dispersion force coefficient
+        //- Turbulent diffusivity
         //  multiplying the gradient of the phase-fraction
-        virtual tmp<volScalarField> Fprime() const;
+        virtual tmp<volScalarField> D() const;
 };
 
 

applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/constantTurbulentDispersionCoefficient/constantTurbulentDispersionCoefficient.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
@@ -70,7 +70,7 @@ Foam::turbulentDispersionModels::constantTurbulentDispersionCoefficient::
 
 Foam::tmp<Foam::volScalarField>
 Foam::turbulentDispersionModels::constantTurbulentDispersionCoefficient::
-Fprime() const
+D() const
 {
     return
         Ctd_

applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/constantTurbulentDispersionCoefficient/constantTurbulentDispersionCoefficient.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
@@ -83,9 +83,9 @@ public:
 
     // Member Functions
 
-        //- Turbulent dispersion force coefficient
+        //- Turbulent diffusivity
         //  multiplying the gradient of the phase-fraction
-        virtual tmp<volScalarField> Fprime() const;
+        virtual tmp<volScalarField> D() const;
 };
 
 

applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/noTurbulentDispersion/noTurbulentDispersion.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,6 +65,36 @@ Foam::turbulentDispersionModels::noTurbulentDispersion::
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
+Foam::tmp<Foam::volScalarField>
+Foam::turbulentDispersionModels::noTurbulentDispersion::D() const
+{
+    const fvMesh& mesh(this->pair_.phase1().mesh());
+
+    return tmp<volScalarField>
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "zero",
+                mesh.time().timeName(),
+                mesh,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE,
+                false
+            ),
+            mesh,
+            dimensionedScalar
+            (
+                "zero",
+                dimensionSet(1, -2, 1, 0, 0),
+                0
+            )
+        )
+    );
+}
+
+
 Foam::tmp<Foam::volVectorField>
 Foam::turbulentDispersionModels::noTurbulentDispersion::F() const
 {
@@ -93,32 +123,4 @@ Foam::turbulentDispersionModels::noTurbulentDispersion::F() const
 }
 
 
-Foam::tmp<Foam::volScalarField>
-Foam::turbulentDispersionModels::noTurbulentDispersion::Fprime() const
-{
-    const fvMesh& mesh(this->pair_.phase1().mesh());
-
-    return
-        tmp<volScalarField>
-        (
-            new volScalarField
-            (
-                IOobject
-                (
-                    "zero",
-                    mesh.time().timeName(),
-                    mesh
-                ),
-                mesh,
-                dimensionedScalar
-                (
-                    "zero",
-                    dimensionSet(1, -2, 1, 0, 0),
-                    0
-                )
-            )
-        );
-}
-
-
 // ************************************************************************* //

applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/noTurbulentDispersion/noTurbulentDispersion.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,12 +76,12 @@ public:
 
     // Member Functions
 
+        //- Turbulent diffusivity
+        //  multiplying the gradient of the phase-fraction
+        virtual tmp<volScalarField> D() const;
+
         //- Turbulent dispersion force
         virtual tmp<volVectorField> F() const;
-
-        //- Turbulent dispersion force coefficient
-        //  multiplying the gradient of the phase-fraction
-        virtual tmp<volScalarField> Fprime() const;
 };
 
 

applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/turbulentDispersionModel/turbulentDispersionModel.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
@@ -35,6 +35,7 @@ namespace Foam
     defineRunTimeSelectionTable(turbulentDispersionModel, dictionary);
 }
 
+const Foam::dimensionSet Foam::turbulentDispersionModel::dimD(1, -1, -2, 0, 0);
 const Foam::dimensionSet Foam::turbulentDispersionModel::dimF(1, -2, -2, 0, 0);
 
 
@@ -61,8 +62,7 @@ Foam::turbulentDispersionModel::~turbulentDispersionModel()
 Foam::tmp<Foam::volVectorField>
 Foam::turbulentDispersionModel::F() const
 {
-    return Fprime()*fvc::grad(pair_.dispersed());
+    return D()*fvc::grad(pair_.dispersed());
-
 }
 
 

applications/solvers/multiphase/twoPhaseEulerFoam/interfacialModels/turbulentDispersionModels/turbulentDispersionModel/turbulentDispersionModel.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
@@ -82,6 +82,9 @@ public:
 
     // Static data members
 
+        //- Diffusivity dimensions
+        static const dimensionSet dimD;
+
         //- Force dimensions
         static const dimensionSet dimF;
 
@@ -111,12 +114,12 @@ public:
 
     // Member Functions
 
+        //- Turbulent diffusivity
+        //  multiplying the gradient of the phase-fraction
+        virtual tmp<volScalarField> D() const = 0;
+
         //- Turbulent dispersion force
         virtual tmp<volVectorField> F() const;
-
-        //- Turbulent dispersion force coefficient
-        //  multiplying the gradient of the phase-fraction
-        virtual tmp<volScalarField> Fprime() const = 0;
 };
 
 

applications/solvers/multiphase/twoPhaseEulerFoam/pEqn.H

@@ -1,267 +0,0 @@
-// --- Pressure corrector loop
-while (pimple.correct())
-{
-    surfaceScalarField alpha1f("alpha1f", fvc::interpolate(alpha1));
-    surfaceScalarField alpha2f("alpha2f", scalar(1) - alpha1f);
-
-    rAU1 = 1.0/U1Eqn.A();
-    rAU2 = 1.0/U2Eqn.A();
-
-    surfaceScalarField rAlphaAU1f(fvc::interpolate(alpha1*rho1*rAU1));
-    surfaceScalarField rAlphaAU2f(fvc::interpolate(alpha2*rho2*rAU2));
-
-    volVectorField HbyA1
-    (
-        IOobject::groupName("HbyA", phase1.name()),
-        U1
-    );
-    HbyA1 = rAU1*U1Eqn.H();
-
-    volVectorField HbyA2
-    (
-        IOobject::groupName("HbyA", phase2.name()),
-        U2
-    );
-    HbyA2 = rAU2*U2Eqn.H();
-
-    mrfZones.makeAbsolute(phi1.oldTime());
-    mrfZones.makeAbsolute(phi1);
-    mrfZones.makeAbsolute(phi2.oldTime());
-    mrfZones.makeAbsolute(phi2);
-
-    // Phase-1 pressure flux (e.g. due to particle-particle pressure)
-    surfaceScalarField phiF1
-    (
-        "phiF1",
-        fvc::interpolate(rAU1*phase1.turbulence().pPrime())
-       *fvc::snGrad(alpha1)*mesh.magSf()
-    );
-
-    // Phase-2 pressure flux (e.g. due to particle-particle pressure)
-    surfaceScalarField phiF2
-    (
-        "phiF2",
-        fvc::interpolate(rAU2*phase2.turbulence().pPrime())
-       *fvc::snGrad(alpha2)*mesh.magSf()
-    );
-
-    volScalarField rho("rho", fluid.rho());
-    surfaceScalarField ghSnGradRho(ghf*fvc::snGrad(rho)*mesh.magSf());
-
-    // Add the phase-1 buoyancy force
-    phiF1 +=
-        rAlphaAU1f
-       *(
-            ghSnGradRho
-          - alpha2f*fvc::interpolate(rho1 - rho2)*(g & mesh.Sf())
-        )/fvc::interpolate(rho1);
-
-    // Add the phase-2 buoyancy force
-    phiF2 +=
-        rAlphaAU2f
-       *(
-            ghSnGradRho
-          - alpha1f*fvc::interpolate(rho2 - rho1)*(g & mesh.Sf())
-        )/fvc::interpolate(rho2);
-
-    // Phase-1 predicted flux
-    surfaceScalarField phiHbyA1
-    (
-        IOobject::groupName("phiHbyA", phase1.name()),
-        (fvc::interpolate(HbyA1) & mesh.Sf())
-      + rAlphaAU1f*fvc::ddtCorr(U1, phi1)
-      + fvc::interpolate(rAU1*dragCoeff)*phi2
-      - phiF1
-    );
-
-    // Phase-2 predicted flux
-    surfaceScalarField phiHbyA2
-    (
-        IOobject::groupName("phiHbyA", phase2.name()),
-        (fvc::interpolate(HbyA2) & mesh.Sf())
-      + rAlphaAU2f*fvc::ddtCorr(U2, phi2)
-      + fvc::interpolate(rAU2*dragCoeff)*phi1
-      - phiF2
-    );
-
-    mrfZones.makeRelative(phiHbyA1);
-    mrfZones.makeRelative(phiHbyA2);
-    mrfZones.makeRelative(phi1.oldTime());
-    mrfZones.makeRelative(phi1);
-    mrfZones.makeRelative(phi2.oldTime());
-    mrfZones.makeRelative(phi2);
-
-    surfaceScalarField phiHbyA("phiHbyA", alpha1f*phiHbyA1 + alpha2f*phiHbyA2);
-
-    HbyA1 += rAU1*dragCoeff*U2;
-    HbyA2 += rAU2*dragCoeff*U1;
-
-    surfaceScalarField rAUf
-    (
-        "rAUf",
-        mag
-        (
-            alpha1f*rAlphaAU1f/fvc::interpolate(rho1)
-          + alpha2f*rAlphaAU2f/fvc::interpolate(rho2)
-        )
-    );
-
-    // Update the fixedFluxPressure BCs to ensure flux consistency
-    setSnGrad<fixedFluxPressureFvPatchScalarField>
-    (
-        p_rgh.boundaryField(),
-        (
-            phiHbyA.boundaryField()
-          - mrfZones.relative
-            (
-                alpha1f.boundaryField()
-               *(mesh.Sf().boundaryField() & U1.boundaryField())
-              + alpha2f.boundaryField()
-               *(mesh.Sf().boundaryField() & U2.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);
-
-    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);
-
-            phi1.boundaryField() ==
-                mrfZones.relative
-                (
-                    mesh.Sf().boundaryField() & U1.boundaryField()
-                );
-            phi1 = phiHbyA1 + rAlphaAU1f*mSfGradp/fvc::interpolate(rho1);
-
-            phi2.boundaryField() ==
-                mrfZones.relative
-                (
-                    mesh.Sf().boundaryField() & U2.boundaryField()
-                );
-            phi2 = phiHbyA2 + rAlphaAU2f*mSfGradp/fvc::interpolate(rho2);
-
-            phi = alpha1f*phi1 + alpha2f*phi2;
-
-            fluid.dgdt() =
-            (
-                alpha1*(pEqnComp2 & p_rgh)
-              - alpha2*(pEqnComp1 & p_rgh)
-            );
-
-            p_rgh.relax();
-
-            p = max(p_rgh + rho*gh, pMin);
-            p_rgh = p - rho*gh;
-
-            mSfGradp = pEqnIncomp.flux()/rAUf;
-
-            U1 = HbyA1
-              + fvc::reconstruct
-                (
-                    rAlphaAU1f*mSfGradp/fvc::interpolate(rho1)
-                  - phiF1
-                );
-            U1.correctBoundaryConditions();
-            fvOptions.correct(U1);
-
-            U2 = HbyA2
-              + fvc::reconstruct
-                (
-                    rAlphaAU2f*mSfGradp/fvc::interpolate(rho2)
-                  - phiF2
-                );
-            U2.correctBoundaryConditions();
-            fvOptions.correct(U2);
-
-            U = fluid.U();
-        }
-    }
-
-    // Update densities from change in p
-    rho1 += psi1*(p_rgh - p_rgh_0);
-    rho2 += psi2*(p_rgh - p_rgh_0);
-
-    K1 = 0.5*magSqr(U1);
-    K2 = 0.5*magSqr(U2);
-
-    if (thermo1.dpdt() || thermo2.dpdt())
-    {
-        dpdt = fvc::ddt(p);
-    }
-}

applications/solvers/multiphase/twoPhaseEulerFoam/pEqnElim.H

@@ -4,8 +4,8 @@ surfaceScalarField alpha2f("alpha2f", scalar(1) - alpha1f);
 rAU1 = 1.0/U1Eqn.A();
 rAU2 = 1.0/U2Eqn.A();
 
-surfaceScalarField rAlphaAU1f(fvc::interpolate(alpha1*rho1*rAU1));
+surfaceScalarField alpharAU1f(fvc::interpolate(alpha1*rAU1));
-surfaceScalarField rAlphaAU2f(fvc::interpolate(alpha2*rho2*rAU2));
+surfaceScalarField alpharAU2f(fvc::interpolate(alpha2*rAU2));
 
 // --- Pressure corrector loop
 while (pimple.correct())
@@ -24,21 +24,43 @@ while (pimple.correct())
     );
     HbyA2 = rAU2*U2Eqn.H();
 
-    // Phase pressure flux (e.g. due to particle-particle pressure)
+    // Face force fluxes
-    surfaceScalarField phiF1
+    tmp<surfaceScalarField> phiF1;
-    (
+    tmp<surfaceScalarField> phiF2;
-        "phiF1",
-        fvc::interpolate(rAU1*phase1.turbulence().pPrime())
-       *fvc::snGrad(alpha1)*mesh.magSf()
-    );
 
-    // Phase-2 pressure flux (e.g. due to particle-particle pressure)
+    // Turbulent diffusion, particle-pressure lift and wall-lubrication fluxes
-    surfaceScalarField phiF2
+    {
-    (
+        volScalarField turbulentDiffusivity(fluid.turbulentDiffusivity());
-        "phiF2",
+        volVectorField liftForce(fluid.liftForce());
-        fvc::interpolate(rAU2*phase2.turbulence().pPrime())
+        volVectorField wallLubricationForce(fluid.wallLubricationForce());
-       *fvc::snGrad(alpha2)*mesh.magSf()
+        surfaceScalarField snGradAlpha1(fvc::snGrad(alpha1)*mesh.magSf());
-    );
+
+        phiF1 =
+        (
+            fvc::interpolate
+            (
+                rAU1*(turbulentDiffusivity + phase1.turbulence().pPrime())
+            )*snGradAlpha1
+
+          + (
+                fvc::interpolate(rAU1*(wallLubricationForce + liftForce))
+              & mesh.Sf()
+            )
+        );
+
+        phiF2 =
+        (
+          - fvc::interpolate
+            (
+                rAU2*(turbulentDiffusivity + phase2.turbulence().pPrime())
+            )*snGradAlpha1
+
+          - (
+                fvc::interpolate(rAU2*(wallLubricationForce + liftForce))
+              & mesh.Sf()
+            )
+        );
+    }
 
     // Mean density for buoyancy force and p_rgh -> p
     volScalarField rho("rho", fluid.rho());
@@ -47,19 +69,19 @@ while (pimple.correct())
     {
         surfaceScalarField ghSnGradRho(ghf*fvc::snGrad(rho)*mesh.magSf());
 
-        phiF1 +=
+        phiF1() +=
-            rAlphaAU1f
+            alpharAU1f
            *(
                 ghSnGradRho
               - alpha2f*fvc::interpolate(rho1 - rho2)*(g & mesh.Sf())
-            )/fvc::interpolate(rho1);
+            );
 
-        phiF2 +=
+        phiF2() +=
-            rAlphaAU2f
+            alpharAU2f
            *(
                 ghSnGradRho
               - alpha1f*fvc::interpolate(rho2 - rho1)*(g & mesh.Sf())
-            )/fvc::interpolate(rho2);
+            );
     }
 
     // ddtPhiCorr filter -- only apply in pure(ish) phases
@@ -86,12 +108,12 @@ while (pimple.correct())
     (
         IOobject::groupName("phiHbyA", phase1.name()),
         (fvc::interpolate(HbyA1) & mesh.Sf())
-      + phiCorrCoeff1*rAlphaAU1f
+      + phiCorrCoeff1*fvc::interpolate(alpha1.oldTime()*rho1.oldTime()*rAU1)
        *(
             mrfZones.absolute(phi1.oldTime())
           - (fvc::interpolate(U1.oldTime()) & mesh.Sf())
         )/runTime.deltaT()
-      - phiF1
+      - phiF1()
     );
 
     // Phase-2 predicted flux
@@ -99,12 +121,12 @@ while (pimple.correct())
     (
         IOobject::groupName("phiHbyA", phase2.name()),
         (fvc::interpolate(HbyA2) & mesh.Sf())
-      + phiCorrCoeff2*rAlphaAU2f
+      + phiCorrCoeff2*fvc::interpolate(alpha2.oldTime()*rho2.oldTime()*rAU2)
        *(
             mrfZones.absolute(phi2.oldTime())
           - (fvc::interpolate(U2.oldTime()) & mesh.Sf())
         )/runTime.deltaT()
-      - phiF2
+      - phiF2()
     );
 
     // Face-drag coefficients
@@ -125,11 +147,7 @@ while (pimple.correct())
     surfaceScalarField rAUf
     (
         "rAUf",
-        mag
+        mag(alpha1f*alpharAU1f + alpha2f*alpharAU2f)
-        (
-            alpha1f*rAlphaAU1f/fvc::interpolate(rho1)
-          + alpha2f*rAlphaAU2f/fvc::interpolate(rho2)
-        )
     );
 
     // Update the fixedFluxPressure BCs to ensure flux consistency
@@ -238,12 +256,12 @@ while (pimple.correct())
             {
                 surfaceScalarField phi1s
                 (
-                    phiHbyA1 + rAlphaAU1f*mSfGradp/fvc::interpolate(rho1)
+                    phiHbyA1 + alpharAU1f*mSfGradp
                 );
 
                 surfaceScalarField phi2s
                 (
-                    phiHbyA2 + rAlphaAU2f*mSfGradp/fvc::interpolate(rho2)
+                    phiHbyA2 + alpharAU2f*mSfGradp
                 );
 
                 surfaceScalarField phir
@@ -286,22 +304,12 @@ while (pimple.correct())
             {
                 volVectorField U1s
                 (
-                    HbyA1
+                    HbyA1 + fvc::reconstruct(alpharAU1f*mSfGradp - phiF1)
-                  + fvc::reconstruct
-                    (
-                        rAlphaAU1f*mSfGradp/fvc::interpolate(rho1)
-                      - phiF1
-                    )
                 );
 
                 volVectorField U2s
                 (
-                    HbyA2
+                    HbyA2 + fvc::reconstruct(alpharAU2f*mSfGradp - phiF2)
-                  + fvc::reconstruct
-                    (
-                        rAlphaAU2f*mSfGradp/fvc::interpolate(rho2)
-                      - phiF2
-                    )
                 );
 
                 volScalarField D1(rAU1*dragCoeff);

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

@@ -226,6 +226,61 @@ Foam::BlendedInterfacialModel<modelType>::F() const
 }
 
 
+template<class modelType>
+Foam::tmp<Foam::volScalarField>
+Foam::BlendedInterfacialModel<modelType>::D() const
+{
+    tmp<volScalarField> f1, f2;
+
+    if (model_.valid() || model1In2_.valid())
+    {
+        f1 = blending_.f1(pair1In2_.dispersed(), pair2In1_.dispersed());
+    }
+
+    if (model_.valid() || model2In1_.valid())
+    {
+        f2 = blending_.f2(pair1In2_.dispersed(), pair2In1_.dispersed());
+    }
+
+    tmp<volScalarField> x
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                modelType::typeName + "Coeff",
+                pair_.phase1().mesh().time().timeName(),
+                pair_.phase1().mesh()
+            ),
+            pair_.phase1().mesh(),
+            dimensionedScalar("zero", modelType::dimD, 0)
+        )
+    );
+
+    if (model_.valid())
+    {
+        x() += model_->D()*(f1() - f2());
+    }
+
+    if (model1In2_.valid())
+    {
+        x() += model1In2_->D()*(1 - f1);
+    }
+
+    if (model2In1_.valid())
+    {
+        x() += model2In1_->D()*f2;
+    }
+
+    if (model_.valid() || model1In2_.valid() || model2In1_.valid())
+    {
+        correctFixedFluxBCs(x());
+    }
+
+    return x;
+}
+
+
 template<class modelType>
 const modelType& Foam::BlendedInterfacialModel<modelType>::phaseModel
 (

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

@@ -113,13 +113,16 @@ public:
 
     // Member Functions
 
-        //- Return the implicit coefficient
+        //- Return the blended force coefficient
         tmp<volScalarField> K() const;
 
-        //- Return the explicit value
+        //- Return the blended force
         template<class Type>
         tmp<GeometricField<Type, fvPatchField, volMesh> > F() 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

@@ -336,6 +336,13 @@ Foam::twoPhaseSystem::turbulentDispersionForce() const
 }
 
 
+Foam::tmp<Foam::volScalarField>
+Foam::twoPhaseSystem::turbulentDiffusivity() const
+{
+    return turbulentDispersion_->D();
+}
+
+
 void Foam::twoPhaseSystem::solve()
 {
     const Time& runTime = mesh_.time();

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

@@ -155,7 +155,11 @@ public:
         //- Return the wall lubrication force
         tmp<volVectorField> wallLubricationForce() const;
 
-        //- Return the wall lubrication force
+        //- Return the turbulent diffusivity
+        //  Multiplies the phase-fraction gradient
+        tmp<volScalarField> turbulentDiffusivity() const;
+
+        //- Return the turbulent dispersion force
         tmp<volVectorField> turbulentDispersionForce() const;
 
         //- Solve for the two-phase-fractions