solvers::multiphaseEuler: New cell momentum/pressure algorithm

The cell-base momentum/pressure algorithm in the multiphaseEuler solver module
has been substantially updated to improve consistency, conservation and reduce
drag generated staggering patterns at sharp interfaces and the boundaries with
stationary phases.  For most if not all cases this new algorithm can be used to
provide well resolved and reliable solutions where the faceMomentum algorithm
would have been chosen previously in order to obtain sufficiently smooth
solutions but at the expense of a noticeable loss in accuracy and resolution.

The first significant change in the momentum/pressure algorithm is in the
interpolation practice used to construct the flux predictor equation from the
cell momentum equation: rather than interpolating the H/A ratio to the faces
i.e. (H/A)_f the terms in the momentum equation are interpolated separately so
that H_f/A_f is used.  The same approach is used for the drag i.e. (D_f/A_f) and
virtual mass contributions.  The advantage of this change is that the phase
forces are now consistent in both the momentum and flux equations, i.e. sum to
zero for each pair of phases.

The second significant change is in the handling of ddtCorr which converts the
old-time time-derivative contributions in H from velocity to flux which is now
consistent due to the change to H/A interpolation and also generalised to use
the fvc::ddtCorr function which has been updated for multiphase.  Additionally
ddtCorr may optionally be applied to the time-derivative in the virtual mass
term in a consistent manner so that the contributions to the flux equation sum
to zero for each pair of phases.

The third significant change is the addition of an optional drag correction term
to the momentum corrector to reduce the staggering patters generated in the
velocity field due to sudden changes in drag force between phase, e.g. at sharp
interfaces between phases or at the boundaries with stationary phases.  This is
particularly beneficial for fluidised bed simulations.  However this correction
is not and cannot be phase consistent, i.e. the correction does not sum to zero
for pairs of phases it is applied to so a small drag error is introduced, but
tests so far have shown that the error is small and outweighed by the benefit in
the reduction in numerical artefacts in the solution.

The final significant change is in the handling of residualAlpha for drag and
virtual mass to provide stable and physical phase velocities in the limit of the
phase-fraction -> 0.  The new approach is phase asymmetric such that the
residual drag is applied only to the phase with a phase-fraction less than
residualAlpha and not to the carrier phase.  This change ensures that the flow
of a pure phase is unaffected by the residualAlpha and residual drag of the
other phases that are stabilised in pure phase region.

There are now four options in the PIMPLE section of the fvSolutions dictionary
relating to the multiphase momentum/pressure algorithm:

PIMPLE
{
    faceMomentum        no;
    VmDdtCorrection     yes;
    dragCorrection      yes;
    partialElimination  no;
}

faceMomentum:
    Switches between the cell and face momentum equation algorithms.
    Provides much smoother and reliable solutions for even the most challenging
    multiphase cases at the expense of a noticeable loss in accuracy and resolution.
    Defaults to 'no'.

VmDdtCorrection:
    Includes the ddtCorr correction term to the time-derivative part of the
    virtual-mass term in the flux equation which ensures consistency between the
    phase virtual mass force on the faces but generates solutions which are
    slightly less smooth and more likely to contain numerical artefacts.
    Defaults to 'no'.

    Testing so far has shown that the loss in smoothness is small and there is
    some noticeable improvement is some cases so in the future the default may
    be changed to 'yes'.

dragCorrection:
    Includes the momentum corrector drag correction term to reduce the
    staggering patters generated in the velocity field due to sudden changes in
    drag force at the expense of a small error in drag consistency.
    Defaults to 'no'

partialElimination:
    Switches the partial-elimination momentum corrector which inverts the drag
    matrix for both the momentum equations and/or flux equations to provide a
    drag implicit correction to the phase velocity and flux fields.  The
    algorithm is the same as previously but updated for the new consistent drag
    interpolation.

All the tutorials/modules/multiphaseEuler tutorial cases have been updated and
tested with the above developments and the four options set appropriately for
each.

applications/solvers/modules/multiphaseEuler/interfacialModels/dragModels/dispersedDragModel/dispersedDragModel.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -44,11 +44,14 @@ SourceFiles
 
 namespace Foam
 {
+
+class swarmCorrection;
+
 namespace dragModels
 {
 
 /*---------------------------------------------------------------------------*\
-                          Class dispersedDragModel Declaration
+                       Class dispersedDragModel Declaration
 \*---------------------------------------------------------------------------*/
 
 class dispersedDragModel

applications/solvers/modules/multiphaseEuler/interfacialModels/dragModels/dragModel/dragModel.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -46,8 +46,6 @@ SourceFiles
 namespace Foam
 {
 
-class swarmCorrection;
-
 /*---------------------------------------------------------------------------*\
                           Class dragModel Declaration
 \*---------------------------------------------------------------------------*/
@@ -83,8 +81,8 @@ public:
         //- Coefficient dimensions
         static const dimensionSet dimK;
 
-        //- Does this model require correcting on fixed flux boundaries?
+        //- This model MUST NOT be set to 0 on fixed flux boundaries
-        static const bool correctFixedFluxBCs = true;
+        static const bool correctFixedFluxBCs = false;
 
 
     // Constructors

applications/solvers/modules/multiphaseEuler/interfacialModels/dragModels/segregated/segregated.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2014-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2014-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -72,29 +72,16 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::segregated::K() const
     const volScalarField& alpha1(interface_.phase1());
     const volScalarField& alpha2(interface_.phase2());
 
-    const volScalarField& rho1(interface_.phase1().rho());
+    const volScalarField::Internal& rho1(interface_.phase1().rho());
-    const volScalarField& rho2(interface_.phase2().rho());
+    const volScalarField::Internal& rho2(interface_.phase2().rho());
 
     tmp<volScalarField> tnu1(interface_.phase1().thermo().nu());
     tmp<volScalarField> tnu2(interface_.phase2().thermo().nu());
 
-    const volScalarField& nu1(tnu1());
+    const volScalarField::Internal& nu1(tnu1());
-    const volScalarField& nu2(tnu2());
+    const volScalarField::Internal& nu2(tnu2());
 
-    volScalarField L
+    const volScalarField::Internal L(cbrt(mesh.V()));
-    (
-        IOobject
-        (
-            "L",
-            mesh.time().name(),
-            mesh
-        ),
-        mesh,
-        dimensionedScalar(dimLength, 0),
-        zeroGradientFvPatchField<scalar>::typeName
-    );
-    L.primitiveFieldRef() = cbrt(mesh.V());
-    L.correctBoundaryConditions();
 
     const dimensionedScalar residualAlpha
     (
@@ -104,44 +91,67 @@ Foam::tmp<Foam::volScalarField> Foam::dragModels::segregated::K() const
         )/2
     );
 
-    const volScalarField I1(alpha1/max(alpha1 + alpha2, residualAlpha));
+    const volScalarField I1
-    const volScalarField I2(alpha2/max(alpha1 + alpha2, residualAlpha));
+    (
-    const volScalarField magGradI
+        alpha1/max(alpha1 + alpha2, residualAlpha)
+    );
+    const volScalarField I2
+    (
+        alpha2/max(alpha1 + alpha2, residualAlpha)
+    );
+    const volScalarField::Internal magGradI
     (
         max
         (
-            (rho2*mag(fvc::grad(I1)) + rho1*mag(fvc::grad(I2)))/(rho1 + rho2),
+            (
+                rho2*mag(fvc::grad(I1)()())
+              + rho1*mag(fvc::grad(I2)()())
+            )/(rho1 + rho2),
             residualAlpha/2/L
         )
     );
 
-    const volScalarField muI(rho1*nu1*rho2*nu2/(rho1*nu1 + rho2*nu2));
+    const volScalarField::Internal muI(rho1*nu1*rho2*nu2/(rho1*nu1 + rho2*nu2));
 
-    const volScalarField limitedAlpha1
+    const volScalarField::Internal limitedAlpha1
     (
         max(alpha1, interface_.phase1().residualAlpha())
     );
 
-    const volScalarField limitedAlpha2
+    const volScalarField::Internal limitedAlpha2
     (
         max(alpha2, interface_.phase2().residualAlpha())
     );
 
-    const volScalarField muAlphaI
+    const volScalarField::Internal muAlphaI
     (
         limitedAlpha1*rho1*nu1*limitedAlpha2*rho2*nu2
        /(limitedAlpha1*rho1*nu1 + limitedAlpha2*rho2*nu2)
     );
 
-    const volScalarField ReI
+    const volScalarField::Internal ReI
     (
-        interface_.rho()*interface_.magUr()
+        (interface_.rho()()()*interface_.magUr()()())
        /(magGradI*limitedAlpha1*limitedAlpha2*muI)
     );
 
-    const volScalarField lambda(m_*ReI + n_*muAlphaI/muI);
+    const volScalarField::Internal lambda(m_*ReI + n_*muAlphaI/muI);
 
-    return lambda*sqr(magGradI)*muI;
+    tmp<volScalarField> tK
+    (
+        volScalarField::New
+        (
+            "K",
+            mesh,
+            dimensionedScalar(dimK, 0),
+            zeroGradientFvPatchField<scalar>::typeName
+        )
+    );
+
+    tK.ref().ref() = lambda*sqr(magGradI)*muI;
+    tK.ref().correctBoundaryConditions();
+
+    return tK;
 }
 
 

applications/solvers/modules/multiphaseEuler/interfacialModels/virtualMassModels/dispersedVirtualMassModel/dispersedVirtualMassModel.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -58,14 +58,12 @@ Foam::virtualMassModels::dispersedVirtualMassModel::Ki() const
 Foam::tmp<Foam::volScalarField>
 Foam::virtualMassModels::dispersedVirtualMassModel::K() const
 {
-    return interface_.dispersed()*Ki();
+    return
-}
+        max
-
+        (
-
+            interface_.dispersed(),
-Foam::tmp<Foam::surfaceScalarField>
+            interface_.dispersed().residualAlpha()
-Foam::virtualMassModels::dispersedVirtualMassModel::Kf() const
+        )*Ki();
-{
-    return fvc::interpolate(interface_.dispersed())*fvc::interpolate(Ki());
 }
 
 

applications/solvers/modules/multiphaseEuler/interfacialModels/virtualMassModels/dispersedVirtualMassModel/dispersedVirtualMassModel.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2011-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2011-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -96,10 +96,6 @@ public:
         //    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;
 };
 
 

applications/solvers/modules/multiphaseEuler/interfacialModels/virtualMassModels/noVirtualMass/noVirtualMass.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2014-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2014-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -72,16 +72,4 @@ Foam::virtualMassModels::noVirtualMass::K() const
 }
 
 
-Foam::tmp<Foam::surfaceScalarField>
-Foam::virtualMassModels::noVirtualMass::Kf() const
-{
-    return surfaceScalarField::New
-    (
-        "Kf",
-        interface_.mesh(),
-        dimensionedScalar(dimK, 0)
-    );
-}
-
-
 // ************************************************************************* //

applications/solvers/modules/multiphaseEuler/interfacialModels/virtualMassModels/noVirtualMass/noVirtualMass.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2014-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2014-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -85,10 +85,6 @@ public:
         //- Return the virtual mass coefficient K
         //  used in the momentum equations
         virtual tmp<volScalarField> K() const;
-
-        //- Return the virtual mass coefficient Kf
-        //  used in the face-momentum equations
-        virtual tmp<surfaceScalarField> Kf() const;
 };
 
 

applications/solvers/modules/multiphaseEuler/interfacialModels/virtualMassModels/virtualMassModel/virtualMassModel.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2014-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2014-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -81,10 +81,4 @@ Foam::tmp<Foam::volScalarField> Foam::blendedVirtualMassModel::K() const
 }
 
 
-Foam::tmp<Foam::surfaceScalarField> Foam::blendedVirtualMassModel::Kf() const
-{
-    return evaluate(&virtualMassModel::Kf, "Kf", virtualMassModel::dimK, false);
-}
-
-
 // ************************************************************************* //

applications/solvers/modules/multiphaseEuler/interfacialModels/virtualMassModels/virtualMassModel/virtualMassModel.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2014-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2014-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -81,8 +81,8 @@ public:
         //- Coefficient dimensions
         static const dimensionSet dimK;
 
-        //- Does this model require correcting on fixed flux boundaries?
+        //- This model should not be set to 0 on fixed flux boundaries
-        static const bool correctFixedFluxBCs = true;
+        static const bool correctFixedFluxBCs = false;
 
 
     // Constructors
@@ -119,10 +119,6 @@ public:
         //    ddt(alpha2*rho2*U2) + ... = ... K*(DU1_Dt - DU2_Dt)
         virtual tmp<volScalarField> K() const = 0;
 
-        //- Return the virtual mass coefficient Kf
-        //  used in the face-momentum equations
-        virtual tmp<surfaceScalarField> Kf() const = 0;
-
         // Dummy write for regIOobject
         bool writeData(Ostream& os) const;
 };
@@ -159,9 +155,6 @@ public:
 
         //- Return the lift coefficient K
         tmp<volScalarField> K() const;
-
-        //- Return the lift coefficient Kf
-        tmp<surfaceScalarField> Kf() const;
 };
 
 

applications/solvers/modules/multiphaseEuler/multiphaseEuler/cellPressureCorrector.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2022-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -58,30 +58,39 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector()
     rAUs.clear();
     rAUs.setSize(phases.size());
 
+    PtrList<surfaceScalarField> rAUfs(phases.size());
+
+    PtrList<volScalarField> Kds(fluid.Kds());
+
     forAll(fluid.movingPhases(), movingPhasei)
     {
         phaseModel& phase = fluid.movingPhases()[movingPhasei];
         const volScalarField& alpha = phase;
 
+        const volScalarField AU
+        (
+            UEqns[phase.index()].A() + Kds[phase.index()]
+          + byDt
+            (
+                max(phase.residualAlpha() - alpha, scalar(0))
+               *phase.rho()
+            )
+        );
+
         rAUs.set
         (
             phase.index(),
             new volScalarField
             (
                 IOobject::groupName("rAU", phase.name()),
-                1.0
+                1/AU
-               /(
-                   UEqns[phase.index()].A()
-                 + byDt
-                   (
-                       max(phase.residualAlpha() - alpha, scalar(0))
-                      *phase.rho()
-                   )
-                )
             )
         );
+
+        rAUfs.set(phase.index(), 1/fvc::interpolate(AU));
     }
     fluid.fillFields("rAU", dimTime/dimDensity, rAUs);
+    fluid.fillFields("rAUf", dimTime/dimDensity, rAUfs);
 
     // Phase diagonal coefficients
     PtrList<surfaceScalarField> alpharAUfs(phases.size());
@@ -94,13 +103,14 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector()
         (
             phasei,
             (
-                fvc::interpolate(max(alpha, phase.residualAlpha())*rAUs[phasei])
+                fvc::interpolate(max(alpha, phase.residualAlpha()))
+               *rAUfs[phasei]
             ).ptr()
         );
     }
 
     // Explicit force fluxes
-    PtrList<surfaceScalarField> phiFs(fluid.phiFs(rAUs));
+    PtrList<surfaceScalarField> Fs(fluid.Fs());
 
     // Mass transfer rates
     PtrList<volScalarField> dmdts(fluid.dmdts());
@@ -134,8 +144,7 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector()
                 (
                     phasei,
                     (
-                        alpharAUfs[phasei]
+                        (
-                       *(
                            ghSnGradRho
                          - (fvc::interpolate(phase.rho() - rho))
                           *(buoyancy.g & mesh.Sf())
@@ -143,11 +152,6 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector()
                         )
                     ).ptr()
                 );
-
-                if (phiFs.set(phasei))
-                {
-                    phigFs[phasei] += phiFs[phasei];
-                }
             }
         }
 
@@ -156,21 +160,22 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector()
         PtrList<surfaceScalarField> phiHbyAs(phases.size());
         {
             // Correction force fluxes
-            PtrList<surfaceScalarField> ddtCorrByAs(fluid.ddtCorrByAs(rAUs));
+            PtrList<surfaceScalarField> ddtCorrs(fluid.ddtCorrs());
 
             forAll(fluid.movingPhases(), movingPhasei)
             {
                 phaseModel& phase = fluid.movingPhases()[movingPhasei];
                 const volScalarField& alpha = phase;
+                const label phasei = phase.index();
 
                 HbyAs.set
                 (
-                    phase.index(),
+                    phasei,
                     constrainHbyA
                     (
-                        rAUs[phase.index()]
+                        rAUs[phasei]
                        *(
-                            UEqns[phase.index()].H()
+                            UEqns[phasei].H()
                           + byDt
                             (
                                 max(phase.residualAlpha() - alpha, scalar(0))
@@ -185,13 +190,17 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector()
 
                 phiHbyAs.set
                 (
-                    phase.index(),
+                    phasei,
                     new surfaceScalarField
                     (
                         IOobject::groupName("phiHbyA", phase.name()),
-                        fvc::flux(HbyAs[phase.index()])
+                        rAUfs[phasei]
-                      - phigFs[phase.index()]
+                       *(
-                      - ddtCorrByAs[phase.index()]
+                           fvc::flux(HbyAs[phasei]/rAUs[phasei])
+                         + ddtCorrs[phasei]
+                        )
+                      - alpharAUfs[phasei]*phigFs[phasei]
+                      - rAUfs[phasei]*Fs[phasei]
                     )
                 );
             }
@@ -200,19 +209,19 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector()
         fluid.fillFields("phiHbyA", dimForce/dimDensity/dimVelocity, phiHbyAs);
 
         // Add explicit drag forces and fluxes
-        PtrList<volVectorField> KdUByAs(fluid.KdUByAs(rAUs));
+        PtrList<volVectorField> KdUs(fluid.KdUs());
-        PtrList<surfaceScalarField> phiKdPhis(fluid.phiKdPhis(rAUs));
+        PtrList<surfaceScalarField> KdPhis(fluid.KdPhis());
 
         forAll(phases, phasei)
         {
-            if (KdUByAs.set(phasei))
+            if (KdUs.set(phasei))
             {
-                HbyAs[phasei] -= KdUByAs[phasei];
+                HbyAs[phasei] -= rAUs[phasei]*KdUs[phasei];
             }
 
-            if (phiKdPhis.set(phasei))
+            if (KdPhis.set(phasei))
             {
-                phiHbyAs[phasei] -= phiKdPhis[phasei];
+                phiHbyAs[phasei] -= rAUfs[phasei]*KdPhis[phasei];
             }
         }
 
@@ -255,6 +264,7 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector()
         forAll(phases, phasei)
         {
             rAUf += alphafs[phasei]*alpharAUfs[phasei];
+            // rAUf += alphafs[phasei]*alphafs[phasei]*rAUfs[phasei];
         }
 
         rAUf = mag(rAUf);
@@ -346,22 +356,54 @@ void Foam::solvers::multiphaseEuler::cellPressureCorrector()
 
                 mSfGradp = pEqnIncomp.flux()/rAUf;
 
-                forAll(fluid.movingPhases(), movingPhasei)
+                if (!dragCorrection)
                 {
-                    phaseModel& phase = fluid.movingPhases()[movingPhasei];
+                    forAll(fluid.movingPhases(), movingPhasei)
+                    {
+                        phaseModel& phase = fluid.movingPhases()[movingPhasei];
+                        const label phasei = phase.index();
 
-                    phase.URef() =
+                        phase.URef() =
-                        HbyAs[phase.index()]
+                            HbyAs[phasei]
-                      + fvc::reconstruct
+                          + fvc::reconstruct
-                        (
+                            (
-                            alpharAUfs[phase.index()]*mSfGradp
+                                alpharAUfs[phasei]*(mSfGradp - phigFs[phasei])
-                          - phigFs[phase.index()]
+                              - rAUfs[phasei]*Fs[phasei]
-                        );
+                            );
+                    }
+                }
+                else
+                {
+                    PtrList<volVectorField> dragCorrs(phases.size());
+                    PtrList<surfaceScalarField> dragCorrfs(phases.size());
+                    fluid.dragCorrs(dragCorrs, dragCorrfs);
+
+                    forAll(fluid.movingPhases(), movingPhasei)
+                    {
+                        phaseModel& phase = fluid.movingPhases()[movingPhasei];
+                        const label phasei = phase.index();
+
+                        phase.URef() =
+                            HbyAs[phasei]
+                          + fvc::reconstruct
+                            (
+                                alpharAUfs[phasei]*(mSfGradp - phigFs[phasei])
+                              + rAUfs[phasei]*(dragCorrfs[phasei] - Fs[phasei])
+                            )
+                        - rAUs[phasei]*dragCorrs[phasei];
+                    }
                 }
 
                 if (partialElimination)
                 {
-                    fluid.partialElimination(rAUs, KdUByAs, alphafs, phiKdPhis);
+                    fluid.partialElimination
+                    (
+                        rAUs,
+                        KdUs,
+                        alphafs,
+                        rAUfs,
+                        KdPhis
+                    );
                 }
                 else
                 {

applications/solvers/modules/multiphaseEuler/multiphaseEuler/facePressureCorrector.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2022-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -71,7 +71,7 @@ void Foam::solvers::multiphaseEuler::facePressureCorrector()
     rAUfs.clear();
     rAUfs.setSize(phases.size());
     {
-        PtrList<surfaceScalarField> AFfs(fluid.AFfs());
+        PtrList<surfaceScalarField> KdVmfs(fluid.KdVmfs());
 
         forAll(fluid.movingPhases(), movingPhasei)
         {
@@ -87,7 +87,7 @@ void Foam::solvers::multiphaseEuler::facePressureCorrector()
                    /(
                         byDt(alphaRho0fs[phase.index()])
                       + fvc::interpolate(UEqns[phase.index()].A())
-                      + AFfs[phase.index()]
+                      + KdVmfs[phase.index()]
                     )
                 )
             );
@@ -111,7 +111,7 @@ void Foam::solvers::multiphaseEuler::facePressureCorrector()
     }
 
     // Explicit force fluxes
-    PtrList<surfaceScalarField> phiFfs(fluid.phiFfs(rAUfs));
+    PtrList<surfaceScalarField> Ffs(fluid.Ffs());
 
     // Mass transfer rates
     PtrList<volScalarField> dmdts(fluid.dmdts());
@@ -155,9 +155,9 @@ void Foam::solvers::multiphaseEuler::facePressureCorrector()
                     ).ptr()
                 );
 
-                if (phiFfs.set(phasei))
+                if (Ffs.set(phasei))
                 {
-                    phigFs[phasei] += phiFfs[phasei];
+                    phigFs[phasei] += rAUfs[phasei]*Ffs[phasei];
                 }
             }
         }
@@ -193,13 +193,13 @@ void Foam::solvers::multiphaseEuler::facePressureCorrector()
         fluid.fillFields("phiHbyA", dimForce/dimDensity/dimVelocity, phiHbyAs);
 
         // Add explicit drag forces and fluxes
-        PtrList<surfaceScalarField> phiKdPhifs(fluid.phiKdPhifs(rAUfs));
+        PtrList<surfaceScalarField> KdPhifs(fluid.KdPhifs());
 
         forAll(phases, phasei)
         {
-            if (phiKdPhifs.set(phasei))
+            if (KdPhifs.set(phasei))
             {
-                phiHbyAs[phasei] -= phiKdPhifs[phasei];
+                phiHbyAs[phasei] -= rAUfs[phasei]*KdPhifs[phasei];
             }
         }
 
@@ -330,7 +330,7 @@ void Foam::solvers::multiphaseEuler::facePressureCorrector()
 
                 if (partialElimination)
                 {
-                    fluid.partialEliminationf(rAUfs, alphafs, phiKdPhifs);
+                    fluid.partialEliminationf(rAUfs, alphafs, KdPhifs);
                 }
                 else
                 {

applications/solvers/modules/multiphaseEuler/multiphaseEuler/momentumPredictor.C

@@ -24,6 +24,7 @@ License
 \*---------------------------------------------------------------------------*/
 
 #include "multiphaseEuler.H"
+#include "fvmSup.H"
 
 // * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * //
 
@@ -37,6 +38,8 @@ void Foam::solvers::multiphaseEuler::cellMomentumPredictor()
     phaseSystem::momentumTransferTable&
         momentumTransfer(momentumTransferPtr());
 
+    const PtrList<volScalarField> Kds(fluid.Kds());
+
     forAll(fluid.movingPhases(), movingPhasei)
     {
         phaseModel& phase = fluid.movingPhases()[movingPhasei];
@@ -54,6 +57,7 @@ void Foam::solvers::multiphaseEuler::cellMomentumPredictor()
              ==
                *momentumTransfer[phase.name()]
               + fvModels().source(alpha, rho, U)
+                // - fvm::Sp(Kds[phase.index()], U)
             )
         );
 

applications/solvers/modules/multiphaseEuler/multiphaseEuler/multiphaseEuler.C

@@ -52,6 +52,9 @@ void Foam::solvers::multiphaseEuler::readControls()
     faceMomentum =
         pimple.dict().lookupOrDefault<Switch>("faceMomentum", false);
 
+    dragCorrection =
+        pimple.dict().lookupOrDefault<Switch>("dragCorrection", false);
+
     partialElimination =
         pimple.dict().lookupOrDefault<Switch>("partialElimination", false);
 
@@ -97,6 +100,11 @@ Foam::solvers::multiphaseEuler::multiphaseEuler(fvMesh& mesh)
         pimple.dict().lookupOrDefault<Switch>("faceMomentum", false)
     ),
 
+    dragCorrection
+    (
+        pimple.dict().lookupOrDefault<Switch>("dragCorrection", false)
+    ),
+
     partialElimination
     (
         pimple.dict().lookupOrDefault<Switch>("partialElimination", false)

applications/solvers/modules/multiphaseEuler/multiphaseEuler/multiphaseEuler.H

@@ -81,6 +81,10 @@ protected:
         //  Defaults to false, i.e. uses the cell momentum equation
         Switch faceMomentum;
 
+        //- Cell/face drag correction for cell momentum corrector
+        //  Defaults to false
+        Switch dragCorrection;
+
         //- Partial elimination drag contribution optimisation
         //  Defaults to false
         Switch partialElimination;

applications/solvers/modules/multiphaseEuler/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2015-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2015-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -63,8 +63,6 @@ class MomentumTransferPhaseSystem
 :
     public BasePhaseSystem
 {
-private:
-
     // Private typedefs
 
         typedef HashPtrTable
@@ -88,13 +86,6 @@ private:
             phaseInterfaceKey::hash
         > VmTable;
 
-        typedef HashPtrTable
-        <
-            surfaceScalarField,
-            phaseInterfaceKey,
-            phaseInterfaceKey::hash
-        > VmfTable;
-
         typedef HashTable
         <
             autoPtr<blendedDragModel>,
@@ -142,9 +133,6 @@ private:
         //- Virtual mass coefficients
         VmTable Vms_;
 
-        //- Face virtual mass coefficients
-        VmfTable Vmfs_;
-
 
         // Sub Models
 
@@ -206,72 +194,33 @@ public:
 
         //- Return implicit force coefficients on the faces, for the face-based
         //  algorithm.
-        virtual PtrList<surfaceScalarField> AFfs() const;
+        virtual PtrList<surfaceScalarField> KdVmfs() const;
 
         //- Return the explicit force fluxes for the cell-based algorithm, that
         //  do not depend on phase mass/volume fluxes, and can therefore be
         //  evaluated outside the corrector loop. This includes things like
         //  lift, turbulent dispersion, and wall lubrication.
-        virtual PtrList<surfaceScalarField> phiFs
+        virtual PtrList<surfaceScalarField> Fs();
-        (
-            const PtrList<volScalarField>& rAUs
-        );
 
-        //- As phiFs, but for the face-based algorithm
+        //- As Fs, but for the face-based algorithm
-        virtual PtrList<surfaceScalarField> phiFfs
+        virtual PtrList<surfaceScalarField> Ffs();
-        (
-            const PtrList<surfaceScalarField>& rAUfs
-        );
 
         //- Return the explicit drag force fluxes for the cell-based algorithm.
         //  These depend on phase mass/volume fluxes, and must therefore be
         //  evaluated inside the corrector loop.
-        virtual PtrList<surfaceScalarField> phiKdPhis
+        virtual PtrList<surfaceScalarField> KdPhis() const;
-        (
-            const PtrList<volScalarField>& rAUs
-        ) const;
 
-        //- As phiKdPhis, but for the face-based algorithm
+        //- As KdPhis, but for the face-based algorithm
-        virtual PtrList<surfaceScalarField> phiKdPhifs
+        virtual PtrList<surfaceScalarField> KdPhifs() const;
-        (
+
-            const PtrList<surfaceScalarField>& rAUfs
+        //- Return the implicit part of the drag force
-        ) const;
+        virtual PtrList<volScalarField> Kds() const;
 
         //- Return the explicit part of the drag force for the cell-based
-        //  algorithm. This is the cell-equivalent of phiKdPhis. These depend on
+        //  algorithm. This is the cell-equivalent of KdPhis. These depend on
         //  phase velocities, and must therefore be evaluated inside the
         //  corrector loop.
-        virtual PtrList<volVectorField> KdUByAs
+        virtual PtrList<volVectorField> KdUs() const;
-        (
-            const PtrList<volScalarField>& rAUs
-        ) const;
-
-        //- Solve the drag system for the velocities and fluxes
-        virtual void partialElimination
-        (
-            const PtrList<volScalarField>& rAUs,
-            const PtrList<volVectorField>& KdUByAs,
-            const PtrList<surfaceScalarField>& alphafs,
-            const PtrList<surfaceScalarField>& phiKdPhis
-        );
-
-        //- As partialElimination, but for the face-based algorithm. Only solves
-        //  for the fluxes.
-        virtual void partialEliminationf
-        (
-            const PtrList<surfaceScalarField>& rAUfs,
-            const PtrList<surfaceScalarField>& alphafs,
-            const PtrList<surfaceScalarField>& phiKdPhifs
-        );
-
-        //- Return the flux corrections for the cell-based algorithm. These
-        //  depend on phase mass/volume fluxes, and must therefore be evaluated
-        //  inside the corrector loop.
-        virtual PtrList<surfaceScalarField> ddtCorrByAs
-        (
-            const PtrList<volScalarField>& rAUs,
-            const bool includeVirtualMass = false
-        ) const;
 
         //- Returns true if the phase pressure is treated implicitly
         //  in the phase fraction equation
@@ -289,6 +238,40 @@ public:
             const PtrList<surfaceScalarField>& rAUfs
         ) const;
 
+        //- Return the flux corrections for the cell-based algorithm. These
+        //  depend on phase mass/volume fluxes, and must therefore be evaluated
+        //  inside the corrector loop.
+        virtual PtrList<surfaceScalarField> ddtCorrs
+        (
+            const bool includeVirtualMass = false
+        ) const;
+
+        //- Set the cell and faces drag correction fields
+        virtual void dragCorrs
+        (
+            PtrList<volVectorField>& dragCorrs,
+            PtrList<surfaceScalarField>& dragCorrf
+        ) const;
+
+        //- Solve the drag system for the velocities and fluxes
+        virtual void partialElimination
+        (
+            const PtrList<volScalarField>& rAUs,
+            const PtrList<volVectorField>& KdUs,
+            const PtrList<surfaceScalarField>& alphafs,
+            const PtrList<surfaceScalarField>& rAUfs,
+            const PtrList<surfaceScalarField>& KdPhis
+        );
+
+        //- As partialElimination, but for the face-based algorithm. Only solves
+        //  for the fluxes.
+        virtual void partialEliminationf
+        (
+            const PtrList<surfaceScalarField>& rAUfs,
+            const PtrList<surfaceScalarField>& alphafs,
+            const PtrList<surfaceScalarField>& KdPhifs
+        );
+
         //- Read base phaseProperties dictionary
         virtual bool read();
 };

applications/solvers/modules/multiphaseEuler/phaseSystems/diameterModels/velocityGroup/velocityGroup.C

@@ -212,7 +212,7 @@ void Foam::diameterModels::velocityGroup::correct()
     const populationBalanceModel& popBal =
         phase().mesh().lookupObject<populationBalanceModel>(popBalName_);
 
-    if (!popBal.solveOnFinalIterOnly() || popBal.pimple().finalIter())
+    if (!popBal.solveOnFinalIterOnly() || popBal.fluid().pimple().finalIter())
     {
         forAll(sizeGroups_, i)
         {

applications/solvers/modules/multiphaseEuler/phaseSystems/phaseModel/MovingPhaseModel/MovingPhaseModel.C

@@ -141,9 +141,7 @@ Foam::MovingPhaseModel<BasePhaseModel>::MovingPhaseModel
         (
             IOobject::groupName("alphaPhi", this->name()),
             fluid.mesh().time().name(),
-            fluid.mesh(),
+            fluid.mesh()
-            IOobject::READ_IF_PRESENT,
-            IOobject::NO_WRITE
         ),
         fluid.mesh(),
         dimensionedScalar(dimensionSet(0, 3, -1, 0, 0), 0)
@@ -154,9 +152,7 @@ Foam::MovingPhaseModel<BasePhaseModel>::MovingPhaseModel
         (
             IOobject::groupName("alphaRhoPhi", this->name()),
             fluid.mesh().time().name(),
-            fluid.mesh(),
+            fluid.mesh()
-            IOobject::READ_IF_PRESENT,
-            IOobject::AUTO_WRITE
         ),
         fluid.mesh(),
         dimensionedScalar(dimensionSet(1, 0, -1, 0, 0), 0)

applications/solvers/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystem.C

@@ -225,6 +225,8 @@ Foam::phaseSystem::phaseSystem
 
     mesh_(mesh),
 
+    pimple_(mesh_.lookupObject<pimpleNoLoopControl>("solutionControl")),
+
     MRF_(mesh_),
 
     referencePhaseName_(lookupOrDefault("referencePhase", word::null)),

applications/solvers/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystem.H

@@ -45,6 +45,8 @@ SourceFiles
 #include "PtrListDictionary.H"
 #include "hashedWordList.H"
 
+#include "pimpleNoLoopControl.H"
+
 #include "IOMRFZoneList.H"
 #include "fvModels.H"
 #include "fvConstraints.H"
@@ -126,6 +128,9 @@ protected:
         //- Reference to the mesh
         const fvMesh& mesh_;
 
+        //- Reference to pimpleNoLoopControl
+        const pimpleNoLoopControl& pimple_;
+
         //- Optional MRF zones
         IOMRFZoneList MRF_;
 
@@ -273,6 +278,9 @@ public:
             //- Return the mesh
             inline const fvMesh& mesh() const;
 
+            //- Return pimpleNoLoopControl
+            inline const pimpleNoLoopControl& pimple() const;
+
             //- Return the phase models
             inline const phaseModelList& phases() const;
 
@@ -515,37 +523,25 @@ public:
 
             //- Return the implicit force coefficients for the face-based
             //  algorithm
-            virtual PtrList<surfaceScalarField> AFfs() const = 0;
+            virtual PtrList<surfaceScalarField> KdVmfs() const = 0;
 
             //- Return the force fluxes for the cell-based algorithm
-            virtual PtrList<surfaceScalarField> phiFs
+            virtual PtrList<surfaceScalarField> Fs() = 0;
-            (
-                const PtrList<volScalarField>& rAUs
-            ) = 0;
 
             //- Return the force fluxes for the face-based algorithm
-            virtual PtrList<surfaceScalarField> phiFfs
+            virtual PtrList<surfaceScalarField> Ffs() = 0;
-            (
-                const PtrList<surfaceScalarField>& rAUfs
-            ) = 0;
 
             //- Return the force fluxes for the cell-based algorithm
-            virtual PtrList<surfaceScalarField> phiKdPhis
+            virtual PtrList<surfaceScalarField> KdPhis() const = 0;
-            (
-                const PtrList<volScalarField>& rAUs
-            ) const = 0;
 
             //- Return the force fluxes for the face-based algorithm
-            virtual PtrList<surfaceScalarField> phiKdPhifs
+            virtual PtrList<surfaceScalarField> KdPhifs() const = 0;
-            (
+
-                const PtrList<surfaceScalarField>& rAUfs
+            //- Return the implicit part of the drag force
-            ) const = 0;
+            virtual PtrList<volScalarField> Kds() const = 0;
 
             //- Return the explicit part of the drag force
-            virtual PtrList<volVectorField> KdUByAs
+            virtual PtrList<volVectorField> KdUs() const = 0;
-            (
-                const PtrList<volScalarField>& rAUs
-            ) const = 0;
 
             //- Returns true if the phase pressure is treated implicitly
             //  in the phase fraction equation
@@ -563,13 +559,27 @@ public:
                 const PtrList<surfaceScalarField>& rAUfs
             ) const = 0;
 
+            //- Return the flux corrections for the cell-based algorithm
+            virtual PtrList<surfaceScalarField> ddtCorrs
+            (
+                const bool includeVirtualMass = false
+            ) const = 0;
+
+            //- Set the cell and faces drag correction fields
+            virtual void dragCorrs
+            (
+                PtrList<volVectorField>& dragCorrs,
+                PtrList<surfaceScalarField>& dragCorrf
+            ) const = 0;
+
             //- Solve the drag system for the new velocities and fluxes
             virtual void partialElimination
             (
                 const PtrList<volScalarField>& rAUs,
-                const PtrList<volVectorField>& KdUByAs,
+                const PtrList<volVectorField>& KdUs,
                 const PtrList<surfaceScalarField>& alphafs,
-                const PtrList<surfaceScalarField>& phiKdPhis
+                const PtrList<surfaceScalarField>& rAUfs,
+                const PtrList<surfaceScalarField>& KdPhis
             ) = 0;
 
             //- Solve the drag system for the new fluxes
@@ -577,7 +587,7 @@ public:
             (
                 const PtrList<surfaceScalarField>& rAUfs,
                 const PtrList<surfaceScalarField>& alphafs,
-                const PtrList<surfaceScalarField>& phiKdPhifs
+                const PtrList<surfaceScalarField>& KdPhifs
             ) = 0;
 
             //- Re-normalise the flux of the phases
@@ -588,13 +598,6 @@ public:
                 const surfaceScalarField& phim
             );
 
-            //- Return the flux corrections for the cell-based algorithm
-            virtual PtrList<surfaceScalarField> ddtCorrByAs
-            (
-                const PtrList<volScalarField>& rAUs,
-                const bool includeVirtualMass = false
-            ) const = 0;
-
             //- Return the heat transfer matrices
             virtual autoPtr<heatTransferTable> heatTransfer() const = 0;
 

applications/solvers/modules/multiphaseEuler/phaseSystems/phaseSystem/phaseSystemI.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2014-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2014-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -31,6 +31,12 @@ inline const Foam::fvMesh& Foam::phaseSystem::mesh() const
 }
 
 
+inline const Foam::pimpleNoLoopControl& Foam::phaseSystem::pimple() const
+{
+    return pimple_;
+}
+
+
 inline const Foam::phaseSystem::phaseModelList&
 Foam::phaseSystem::phases() const
 {

applications/solvers/modules/multiphaseEuler/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.C

@@ -797,7 +797,6 @@ Foam::diameterModels::populationBalanceModel::populationBalanceModel
     (
         fluid_.subDict("populationBalanceCoeffs").subDict(name_)
     ),
-    pimple_(mesh_.lookupObject<pimpleNoLoopControl>("solutionControl")),
     continuousPhase_
     (
         mesh_.lookupObject<phaseModel>
@@ -1130,7 +1129,7 @@ Foam::diameterModels::populationBalanceModel::continuousTurbulence() const
 
 void Foam::diameterModels::populationBalanceModel::solve()
 {
-    if (!solveOnFinalIterOnly() || pimple_.finalIter())
+    if (!solveOnFinalIterOnly() || fluid_.pimple().finalIter())
     {
         const label nCorr = this->nCorr();
         const scalar tolerance =

applications/solvers/modules/multiphaseEuler/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModel.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2017-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2017-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -206,7 +206,6 @@ SourceFiles
 
 #include "sizeGroup.H"
 #include "phaseSystem.H"
-#include "pimpleNoLoopControl.H"
 #include "phaseCompressibleMomentumTransportModel.H"
 #include "HashPtrTable.H"
 #include "Pair.H"
@@ -232,8 +231,6 @@ class populationBalanceModel
 :
     public regIOobject
 {
-private:
-
     // Private Data
 
         //- Reference to the phaseSystem
@@ -251,9 +248,6 @@ private:
         //- Dictionary
         dictionary dict_;
 
-        //- Reference to pimpleNoLoopControl
-        const pimpleNoLoopControl& pimple_;
-
         //- Continuous phase
         const phaseModel& continuousPhase_;
 
@@ -438,9 +432,6 @@ public:
         //- Return populationBalanceCoeffs dictionary
         inline const dictionary& dict() const;
 
-        //- Return pimpleNoLoopControl
-        const pimpleNoLoopControl& pimple() const;
-
         //- Return the number of corrections
         inline label nCorr() const;
 

applications/solvers/modules/multiphaseEuler/phaseSystems/populationBalanceModel/populationBalanceModel/populationBalanceModelI.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2017-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2017-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -64,13 +64,6 @@ Foam::diameterModels::populationBalanceModel::dict() const
 }
 
 
-inline const Foam::pimpleNoLoopControl&
-Foam::diameterModels::populationBalanceModel::pimple() const
-{
-    return pimple_;
-}
-
-
 inline Foam::label Foam::diameterModels::populationBalanceModel::nCorr() const
 {
     return mesh_.solution().solverDict(name_).lookup<label>("nCorr");

tutorials/modules/multiphaseEuler/Grossetete/system/fvSolution

@@ -74,7 +74,10 @@ PIMPLE
     nNonOrthogonalCorrectors 0;
     nEnergyCorrectors   1;
 
-    faceMomentum        yes;
+    faceMomentum        no;
+    VmDdtCorrection     yes;
+    dragCorrection      yes;
+    partialElimination  no;
 }
 
 relaxationFactors

tutorials/modules/multiphaseEuler/LBend/constant/phaseProperties

@@ -69,13 +69,7 @@ drag
 }
 
 virtualMass
-{
+{}
-    solids_dispersedIn_gas
-    {
-        type            constantCoefficient;
-        Cvm             0;
-    }
-}
 
 heatTransfer
 {

tutorials/modules/multiphaseEuler/LBend/system/fvSolution

@@ -89,8 +89,10 @@ PIMPLE
     nCorrectors      1;
     nNonOrthogonalCorrectors 0;
 
-    faceMomentum     yes;
+    faceMomentum        no;
-
+    VmDdtCorrection     no;
+    dragCorrection      yes;
+    partialElimination  no;
 }
 
 relaxationFactors

tutorials/modules/multiphaseEuler/bubbleColumn/0/p_rgh

@@ -27,7 +27,7 @@ boundaryField
     outlet
     {
         type            prghPressure;
-        p                  $internalField;
+        p               $internalField;
         value           $internalField;
     }
     walls

tutorials/modules/multiphaseEuler/bubbleColumn/system/fvSolution

@@ -69,6 +69,11 @@ PIMPLE
     nOuterCorrectors 3;
     nCorrectors      1;
     nNonOrthogonalCorrectors 0;
+
+    faceMomentum        no;
+    VmDdtCorrection     yes;
+    dragCorrection      yes;
+    partialElimination  no;
 }
 
 relaxationFactors

tutorials/modules/multiphaseEuler/bubbleColumnEvaporating/system/fvSolution

@@ -69,6 +69,11 @@ PIMPLE
     nOuterCorrectors 5;
     nCorrectors      1;
     nNonOrthogonalCorrectors 0;
+
+    faceMomentum        no;
+    VmDdtCorrection     yes;
+    dragCorrection      yes;
+    partialElimination  no;
 }
 
 relaxationFactors

tutorials/modules/multiphaseEuler/bubbleColumnEvaporatingDissolving/system/fvSolution

@@ -66,6 +66,11 @@ PIMPLE
     nOuterCorrectors 5;
     nCorrectors      1;
     nNonOrthogonalCorrectors 0;
+
+    faceMomentum        no;
+    VmDdtCorrection     yes;
+    dragCorrection      yes;
+    partialElimination  no;
 }
 
 relaxationFactors

tutorials/modules/multiphaseEuler/bubbleColumnEvaporatingReacting/system/fvSolution

@@ -78,6 +78,11 @@ PIMPLE
     nOuterCorrectors 5;
     nCorrectors      1;
     nNonOrthogonalCorrectors 0;
+
+    faceMomentum        no;
+    VmDdtCorrection     yes;
+    dragCorrection      yes;
+    partialElimination  no;
 }
 
 relaxationFactors

tutorials/modules/multiphaseEuler/bubbleColumnIATE/system/fvSolution

@@ -60,6 +60,10 @@ PIMPLE
     nCorrectors      1;
     nNonOrthogonalCorrectors 0;
 
+    faceMomentum        no;
+    VmDdtCorrection     yes;
+    dragCorrection      yes;
+    partialElimination  no;
 }
 
 relaxationFactors

tutorials/modules/multiphaseEuler/bubbleColumnLES/system/fvSolution

@@ -69,6 +69,11 @@ PIMPLE
     nOuterCorrectors 3;
     nCorrectors      1;
     nNonOrthogonalCorrectors 0;
+
+    faceMomentum        no;
+    VmDdtCorrection     yes;
+    dragCorrection      yes;
+    partialElimination  no;
 }
 
 relaxationFactors

tutorials/modules/multiphaseEuler/bubbleColumnLaminar/system/fvSolution

@@ -61,6 +61,10 @@ PIMPLE
     nCorrectors      1;
     nNonOrthogonalCorrectors 0;
 
+    faceMomentum        no;
+    VmDdtCorrection     yes;
+    dragCorrection      yes;
+    partialElimination  no;
 }
 
 relaxationFactors

tutorials/modules/multiphaseEuler/bubblePipe/system/fvSolution

@@ -70,8 +70,10 @@ PIMPLE
     nCorrectors         2;
     nNonOrthogonalCorrectors 0;
 
-    faceMomentum        true;
+    faceMomentum        no;
-
+    VmDdtCorrection     yes;
+    dragCorrection      yes;
+    partialElimination  no;
 }
 
 relaxationFactors

tutorials/modules/multiphaseEuler/damBreak4phase/system/controlDict

@@ -34,7 +34,7 @@ writeInterval   0.02;
 
 purgeWrite      0;
 
-writeFormat     binary;
+writeFormat     ascii;
 
 writePrecision  6;
 
@@ -53,5 +53,9 @@ maxAlphaCo      0.5;
 
 maxDeltaT       1;
 
+functions
+{
+    #includeFunc phaseMap
+}
 
 // ************************************************************************* //

tutorials/modules/multiphaseEuler/damBreak4phase/system/fvSolution

@@ -76,6 +76,11 @@ PIMPLE
 {
     nCorrectors         3;
     nNonOrthogonalCorrectors 0;
+
+    faceMomentum        no;
+    VmDdtCorrection     no;
+    dragCorrection      no;
+    partialElimination  no;
 }
 
 relaxationFactors

tutorials/modules/multiphaseEuler/fluidisedBed/system/fvSolution

@@ -88,8 +88,9 @@ PIMPLE
     nCorrectors      2;
     nNonOrthogonalCorrectors 0;
 
-    faceMomentum     no;
+    faceMomentum        no;
-
+    dragCorrection      yes;
+    partialElimination  no;
 }
 
 relaxationFactors

tutorials/modules/multiphaseEuler/fluidisedBedLaminar/constant/phaseProperties

@@ -66,20 +66,14 @@ drag
 }
 
 virtualMass
-{
+{}
-    particles_dispersedIn_air
-    {
-        type            constantCoefficient;
-        Cvm             0.5;
-    }
-}
 
 heatTransfer
 {
     particles_dispersedIn_air
     {
         type            RanzMarshall;
-        residualAlpha   1e-3;
+        residualAlpha   1e-4;
     }
 }
 

tutorials/modules/multiphaseEuler/fluidisedBedLaminar/system/fvSolution

@@ -84,12 +84,13 @@ solvers
 
 PIMPLE
 {
-    nOuterCorrectors 3;
+    nOuterCorrectors    3;
-    nCorrectors      2;
+    nCorrectors         2;
     nNonOrthogonalCorrectors 0;
 
-    faceMomentum     no;
+    faceMomentum        no;
-
+    dragCorrection      yes;
+    partialElimination  no;
 }
 
 relaxationFactors

tutorials/modules/multiphaseEuler/hydrofoil/system/fvSolution

@@ -81,6 +81,9 @@ PIMPLE
     transportCorrectionFinal yes;
 
     faceMomentum        no;
+    VmDdtCorrection     no;
+    dragCorrection      yes;
+    partialElimination  no;
 }
 
 cache

tutorials/modules/multiphaseEuler/injection/constant/phaseProperties

@@ -65,13 +65,7 @@ blending
 }
 
 surfaceTension
-{
+{}
-    air_water
-    {
-        type            constant;
-        sigma           0.07;
-    }
-}
 
 drag
 {
@@ -103,8 +97,7 @@ virtualMass
 
     water_dispersedIn_air
     {
-        type            constantCoefficient;
+        type            none;
-        Cvm             0.5;
     }
 }
 

tutorials/modules/multiphaseEuler/injection/system/fvSolution

@@ -61,6 +61,10 @@ PIMPLE
     nCorrectors      1;
     nNonOrthogonalCorrectors 0;
 
+    faceMomentum        no;
+    VmDdtCorrection     yes;
+    dragCorrection      yes;
+    partialElimination  no;
 }
 
 relaxationFactors

tutorials/modules/multiphaseEuler/mixerVessel2D/constant/phaseProperties

@@ -33,7 +33,7 @@ water
 oil
 {
     type            pureIsothermalPhaseModel;
-    diameterModel constant;
+    diameterModel   constant;
     constantCoeffs
     {
         d               1e-3;
@@ -45,7 +45,7 @@ oil
 mercury
 {
     type            pureIsothermalPhaseModel;
-    diameterModel constant;
+    diameterModel   constant;
     constantCoeffs
     {
         d               1e-3;

tutorials/modules/multiphaseEuler/mixerVessel2D/constant/physicalProperties.mercury

@@ -22,7 +22,7 @@ thermoType
     thermo          hConst;
     equationOfState rhoConst;
     specie          specie;
-    energy          sensibleEnthalpy;
+    energy          sensibleInternalEnergy;
 }
 
 mixture

tutorials/modules/multiphaseEuler/mixerVessel2D/constant/physicalProperties.oil

@@ -22,7 +22,7 @@ thermoType
     thermo          hConst;
     equationOfState rhoConst;
     specie          specie;
-    energy          sensibleEnthalpy;
+    energy          sensibleInternalEnergy;
 }
 
 mixture

tutorials/modules/multiphaseEuler/mixerVessel2D/constant/physicalProperties.water

@@ -22,7 +22,7 @@ thermoType
     thermo          eConst;
     equationOfState rPolynomial;
     specie          specie;
-    energy          sensibleEnthalpy;
+    energy          sensibleInternalEnergy;
 }
 
 mixture

tutorials/modules/multiphaseEuler/mixerVessel2D/system/fvSolution

@@ -94,6 +94,11 @@ PIMPLE
 
     correctPhi          yes;
     correctMeshPhi      no;
+
+    faceMomentum        no;
+    VmDdtCorrection     yes;
+    dragCorrection      yes;
+    partialElimination  no;
 }
 
 relaxationFactors

tutorials/modules/multiphaseEuler/mixerVessel2DMRF/system/fvSolution

@@ -81,6 +81,11 @@ PIMPLE
 
     pRefCell        0;
     pRefValue       0;
+
+    faceMomentum        no;
+    VmDdtCorrection     yes;
+    dragCorrection      yes;
+    partialElimination  no;
 }
 
 relaxationFactors

tutorials/modules/multiphaseEuler/pipeBend/system/fvSolution

@@ -53,19 +53,19 @@ solvers
 
     "(k|epsilon).*"
     {
-      solver            PBiCGStab;
+        solver            PBiCGStab;
-      preconditioner    DILU;
+        preconditioner    DILU;
-      tolerance         1e-5;
+        tolerance         1e-5;
-      relTol            0;
+        relTol            0;
-      minIter           1;
+        minIter           1;
     }
 
     "f.*"
     {
-      solver            PBiCGStab;
+        solver            PBiCGStab;
-      preconditioner    DILU;
+        preconditioner    DILU;
-      tolerance         1e-6;
+        tolerance         1e-6;
-      relTol            0;
+        relTol            0;
     }
 
     agglomerates
@@ -85,6 +85,11 @@ PIMPLE
     nNonOrthogonalCorrectors 0;
     alphaSpreadMin           5e-4;
     alphaSpreadMax           0.9995;
+
+    faceMomentum        no;
+    VmDdtCorrection     yes;
+    dragCorrection      yes;
+    partialElimination  no;
 }
 
 relaxationFactors

tutorials/modules/multiphaseEuler/steamInjection/0/T.steam

@@ -20,10 +20,12 @@ internalField       uniform 372.76;
 
 boundaryField
 {
-    walls
+    inlet
     {
-        type            zeroGradient;
+        type            fixedValue;
+        value           $internalField;
     }
+
     outlet
     {
         type            inletOutlet;
@@ -31,10 +33,10 @@ boundaryField
         inletValue      $internalField;
         value           $internalField;
     }
-    inlet
+
+    walls
     {
-        type            fixedValue;
+        type            zeroGradient;
-        value           $internalField;
     }
 }
 

tutorials/modules/multiphaseEuler/steamInjection/0/T.water

@@ -20,11 +20,11 @@ internalField       uniform 360;
 
 boundaryField
 {
-    walls
+    inlet
     {
-        type            fixedValue;
+        type             zeroGradient;
-        value           $internalField;
     }
+
     outlet
     {
         type            inletOutlet;
@@ -32,9 +32,11 @@ boundaryField
         inletValue      $internalField;
         value           $internalField;
     }
-    inlet
+
+    walls
     {
-        type             zeroGradient;
+        type            fixedValue;
+        value           $internalField;
     }
 }
 

tutorials/modules/multiphaseEuler/steamInjection/0/U.steam

@@ -16,21 +16,23 @@ FoamFile
 
 dimensions      [0 1 -1 0 0 0 0];
 
-internalField   uniform (0 0.1 0);
+internalField   uniform (0 0 0);
 
 boundaryField
 {
     inlet
     {
         type            fixedValue;
-        value           uniform (0 0.1 0);
+        value           uniform (0 0 0);
     }
+
     outlet
     {
         type            pressureInletOutletVelocity;
         phi             phi.steam;
         value           $internalField;
     }
+
     walls
     {
         type            noSlip;

tutorials/modules/multiphaseEuler/steamInjection/0/U.water

@@ -25,12 +25,14 @@ boundaryField
         type            fixedValue;
         value           uniform (0 0.1 0);
     }
+
     outlet
     {
         type            pressureInletOutletVelocity;
         phi             phi.water;
         value           $internalField;
     }
+
     walls
     {
         type            noSlip;

tutorials/modules/multiphaseEuler/steamInjection/0/alpha.steam

@@ -23,8 +23,9 @@ boundaryField
     inlet
     {
         type            fixedValue;
-        value           uniform 0.;
+        value           uniform 0;
     }
+
     outlet
     {
         type            inletOutlet;
@@ -32,14 +33,11 @@ boundaryField
         inletValue      uniform 0;
         value           uniform 0;
     }
+
     walls
     {
         type            zeroGradient;
     }
-    defaultFaces
-    {
-        type            empty;
-    }
 }
 
 

tutorials/modules/multiphaseEuler/steamInjection/0/alpha.water

@@ -23,23 +23,21 @@ boundaryField
     inlet
     {
         type            fixedValue;
-        value           uniform 1.;
+        value           uniform 1;
     }
+
     outlet
     {
         type            inletOutlet;
         phi             phi.water;
-        inletValue      uniform 1.0;
+        inletValue      uniform 1;
-        value           uniform 1.0;
+        value           uniform 1;
     }
+
     walls
     {
         type            zeroGradient;
     }
-    defaultFaces
-    {
-        type            empty;
-    }
 }
 
 

tutorials/modules/multiphaseEuler/steamInjection/0/nut.water

@@ -25,11 +25,13 @@ boundaryField
         type            calculated;
         value           $internalField;
     }
+
     outlet
     {
         type            calculated;
         value           $internalField;
     }
+
     walls
     {
         type            nutkWallFunction;

tutorials/modules/multiphaseEuler/steamInjection/0/p

@@ -24,11 +24,13 @@ boundaryField
         type            calculated;
         value           $internalField;
     }
+
     outlet
     {
         type            calculated;
         value           $internalField;
     }
+
     walls
     {
         type            calculated;

tutorials/modules/multiphaseEuler/steamInjection/0/p_rgh

@@ -24,12 +24,17 @@ boundaryField
         type            fixedFluxPressure;
         value           $internalField;
     }
+
     outlet
     {
-        type            prghPressure;
+        type            prghTotalPressure;
-        p               $internalField;
+        p0              $internalField;
+        U               U.water;
+        phi             phi.water;
+        rho             rho.water;
         value           $internalField;
     }
+
     walls
     {
         type            fixedFluxPressure;

tutorials/modules/multiphaseEuler/steamInjection/0/water.steam

@@ -1,41 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Version:  dev
-     \\/     M anipulation  |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    format      ascii;
-    class       volScalarField;
-    location    "0";
-    object      AIR.steam;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-dimensions      [0 0 0 0 0 0 0];
-
-internalField   uniform 1;
-
-boundaryField
-{
-    inlet
-    {
-        type            fixedValue;
-        value           uniform 1;
-    }
-    outlet
-    {
-        type            inletOutlet;
-        phi             phi.steam;
-        inletValue      $internalField;
-        value           $internalField;
-    }
-    walls
-    {
-        type            zeroGradient;
-    }
-}
-
-// ************************************************************************* //

tutorials/modules/multiphaseEuler/steamInjection/0/water.water

@@ -1,41 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Version:  dev
-     \\/     M anipulation  |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    format      ascii;
-    class       volScalarField;
-    location    "0";
-    object      AIR.water;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-dimensions      [0 0 0 0 0 0 0];
-
-internalField   uniform 1;
-
-boundaryField
-{
-    inlet
-    {
-        type            fixedValue;
-        value           $internalField;
-    }
-    outlet
-    {
-        type            inletOutlet;
-        phi             phi.water;
-        inletValue      $internalField;
-        value           $internalField;
-    }
-    walls
-    {
-        type            zeroGradient;
-    }
-}
-
-// ************************************************************************* //

tutorials/modules/multiphaseEuler/steamInjection/constant/phaseProperties

@@ -51,13 +51,7 @@ blending
 }
 
 surfaceTension
-{
+{}
-    steam_water
-    {
-        type            constant;
-        sigma           0.07;
-    }
-}
 
 saturationTemperature
 {
@@ -90,13 +84,13 @@ heatTransfer
     steam_dispersedIn_water_inThe_steam
     {
         type spherical;
-        residualAlpha 1e-4;
+        residualAlpha 1e-3;
     }
 
     steam_dispersedIn_water_inThe_water
     {
         type RanzMarshall;
-        residualAlpha 1e-4;
+        residualAlpha 1e-3;
     }
 }
 

tutorials/modules/multiphaseEuler/steamInjection/system/fvSolution

@@ -24,10 +24,12 @@ solvers
 
     p_rgh
     {
-        solver          PCG;
+        solver          GAMG;
-        preconditioner  DIC;
+        smoother        DIC;
-        tolerance       1e-10;
+
-        relTol          0.001;
+        tolerance       1e-8;
+        relTol          0.01;
+        maxIter         20;
     }
 
     p_rghFinal
@@ -38,50 +40,53 @@ solvers
 
     "U.*"
     {
-        solver          smoothSolver;
+        solver          PBiCGStab;
-        smoother        symGaussSeidel;
+        preconditioner  DILU;
+
         tolerance       1e-6;
         relTol          0;
-        minIter         1;
     }
 
     "(e|h).*"
     {
-        solver          smoothSolver;
+        solver          PBiCGStab;
-        smoother        symGaussSeidel;
+        preconditioner  DILU;
-        tolerance       1e-6;
+
+        tolerance       1e-8;
         relTol          0;
-        minIter         1;
         maxIter         20;
     }
 
     "(k|epsilon).*"
     {
-        solver          smoothSolver;
+        solver          PBiCGStab;
-        smoother        symGaussSeidel;
+        preconditioner  DILU;
+
         tolerance       1e-8;
         relTol          0;
-        minIter         1;
     }
 
     "Yi.*"
     {
-        solver          smoothSolver;
+        solver          PBiCGStab;
-        smoother        symGaussSeidel;
+        preconditioner  DILU;
-        tolerance       1e-6;
+
+        tolerance       1e-8;
         relTol          0;
-        minIter         1;
     }
 }
 
 PIMPLE
 {
-    nOuterCorrectors 3;
+    nOuterCorrectors    2;
-    nCorrectors      1;
+    nCorrectors         2;
     nNonOrthogonalCorrectors 0;
-    nEnergyCorrectors 2;
+    nEnergyCorrectors   2;
 
-    faceMomentum    yes;
+    faceMomentum        no;
+    VmDdtCorrection     yes;
+    dragCorrection      yes;
+    partialElimination  yes;
 }
 
 relaxationFactors

tutorials/modules/multiphaseEuler/titaniaSynthesis/system/fvSolution

@@ -83,7 +83,10 @@ PIMPLE
     nEnergyCorrectors   1;
     nNonOrthogonalCorrectors 0;
 
-    faceMomentum        yes;
+    faceMomentum        no;
+    VmDdtCorrection     no;
+    dragCorrection      yes;
+    partialElimination  no;
 }
 
 relaxationFactors

tutorials/modules/multiphaseEuler/titaniaSynthesisSurface/system/fvSolution

@@ -83,7 +83,10 @@ PIMPLE
     nEnergyCorrectors   1;
     nNonOrthogonalCorrectors 0;
 
-    faceMomentum        yes;
+    faceMomentum        no;
+    VmDdtCorrection     no;
+    dragCorrection      yes;
+    partialElimination  no;
 }
 
 relaxationFactors

tutorials/modules/multiphaseEuler/wallBoilingIATE/system/fvSolution

@@ -73,8 +73,10 @@ PIMPLE
     nNonOrthogonalCorrectors 0;
     nEnergyCorrectors   1;
 
-    faceMomentum        yes;
+    faceMomentum        no;
-
+    VmDdtCorrection     yes;
+    dragCorrection      yes;
+    partialElimination  no;
 }
 
 relaxationFactors

tutorials/modules/multiphaseEuler/wallBoilingPolydisperse/system/fvSolution

@@ -82,8 +82,10 @@ PIMPLE
     nNonOrthogonalCorrectors 0;
     nEnergyCorrectors   1;
 
-    faceMomentum        yes;
+    faceMomentum        no;
-
+    VmDdtCorrection     yes;
+    dragCorrection      yes;
+    partialElimination  no;
 }
 
 relaxationFactors

tutorials/modules/multiphaseEuler/wallBoilingPolydisperseTwoGroups/system/fvSolution

@@ -82,8 +82,10 @@ PIMPLE
     nNonOrthogonalCorrectors 0;
     nEnergyCorrectors   1;
 
-    faceMomentum        yes;
+    faceMomentum        no;
-
+    VmDdtCorrection     yes;
+    dragCorrection      yes;
+    partialElimination  no;
 }
 
 relaxationFactors