multiphaseEulerFoam, turbulentDispersionModel: Corrected face-force for multiphase case

This is a completion of commit 64da7a2c. The fix has now also been
applied to the face-momentum equation.

applications/solvers/multiphase/multiphaseEulerFoam/functionObjects/phaseForces/phaseForces.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2018-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2018-2022 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -26,6 +26,7 @@ License
 #include "phaseForces.H"
 #include "addToRunTimeSelectionTable.H"
 #include "BlendedInterfacialModel.H"
+#include "fvcGrad.H"
 #include "dragModel.H"
 #include "virtualMassModel.H"
 #include "liftModel.H"
@@ -277,7 +278,12 @@ bool Foam::functionObjects::phaseForces::execute()
             if (fluid_.foundBlendedSubModel<turbulentDispersionModel>(pair))
             {
                 *forceFields_[turbulentDispersionModel::typeName] +=
-                    nonDragForce<turbulentDispersionModel>(pair);
+                    fluid_.lookupBlendedSubModel<turbulentDispersionModel>
+                    (
+                        pair
+                    ).D()
+                   *(&pair.phase1() == &phase_ ? -1 : +1)
+                   *fvc::grad(pair.phase1()/(pair.phase1() + pair.phase2()));
             }
         }
     }

applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/turbulentDispersionModels/noTurbulentDispersion/noTurbulentDispersion.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2014-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2014-2022 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -79,18 +79,4 @@ Foam::turbulentDispersionModels::noTurbulentDispersion::D() const
 }
 
 
-Foam::tmp<Foam::volVectorField>
-Foam::turbulentDispersionModels::noTurbulentDispersion::F() const
-{
-    const fvMesh& mesh(this->pair_.phase1().mesh());
-
-    return volVectorField::New
-    (
-        "zero",
-        mesh,
-        dimensionedVector(dimF, Zero)
-    );
-}
-
-
 // ************************************************************************* //

applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/turbulentDispersionModels/noTurbulentDispersion/noTurbulentDispersion.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2014-2020 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2014-2022 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -79,9 +79,6 @@ public:
         //- Turbulent diffusivity
         //  multiplying the gradient of the phase-fraction
         virtual tmp<volScalarField> D() const;
-
-        //- Turbulent dispersion force
-        virtual tmp<volVectorField> F() const;
 };
 
 

applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/turbulentDispersionModels/turbulentDispersionModel/turbulentDispersionModel.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2014-2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2014-2022 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -25,9 +25,6 @@ License
 
 #include "turbulentDispersionModel.H"
 #include "phasePair.H"
-#include "fvcGrad.H"
-#include "surfaceInterpolate.H"
-#include "fvcSnGrad.H"
 #include "phaseCompressibleMomentumTransportModel.H"
 #include "BlendedInterfacialModel.H"
 
@@ -41,7 +38,6 @@ namespace Foam
 }
 
 const Foam::dimensionSet Foam::turbulentDispersionModel::dimD(1, -1, -2, 0, 0);
-const Foam::dimensionSet Foam::turbulentDispersionModel::dimF(1, -2, -2, 0, 0);
 
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
@@ -82,41 +78,4 @@ Foam::turbulentDispersionModel::continuousTurbulence() const
 }
 
 
-Foam::tmp<Foam::volVectorField>
-Foam::turbulentDispersionModel::F() const
-{
-    return
-        D()
-       *fvc::grad
-        (
-            pair_.dispersed()
-           /max
-            (
-                pair_.dispersed() + pair_.continuous(),
-                pair_.dispersed().residualAlpha()
-            )
-        );
-}
-
-
-Foam::tmp<Foam::surfaceScalarField>
-Foam::turbulentDispersionModel::Ff() const
-{
-    return
-    pair_.phase1().mesh().magSf()
-   *(
-        fvc::interpolate(D())
-       *fvc::snGrad
-        (
-            pair_.dispersed()
-           /max
-            (
-                pair_.dispersed() + pair_.continuous(),
-                pair_.dispersed().residualAlpha()
-            )
-        )
-    );
-}
-
-
 // ************************************************************************* //

applications/solvers/multiphase/multiphaseEulerFoam/interfacialModels/turbulentDispersionModels/turbulentDispersionModel/turbulentDispersionModel.H

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2014-2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2014-2022 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -86,9 +86,6 @@ public:
         //- Diffusivity dimensions
         static const dimensionSet dimD;
 
-        //- Force dimensions
-        static const dimensionSet dimF;
-
 
     // Constructors
 
@@ -122,12 +119,6 @@ public:
         //- 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 on faces
-        virtual tmp<surfaceScalarField> Ff() const;
 };
 
 

applications/solvers/multiphase/multiphaseEulerFoam/phaseSystems/PhaseSystems/MomentumTransferPhaseSystem/MomentumTransferPhaseSystem.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2015-2021 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2015-2022 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -753,24 +753,35 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::phiFfs
         turbulentDispersionModelIter
     )
     {
-        const surfaceScalarField Ff(turbulentDispersionModelIter()->Ff());
         const phasePair&
             pair(this->phasePairs_[turbulentDispersionModelIter.key()]);
 
-        addField
+        const surfaceScalarField Df
         (
-            pair.phase1(),
+            fvc::interpolate(turbulentDispersionModelIter()->D())
-            "phiFf",
-            rAUfs[pair.phase1().index()]*Ff,
-            phiFfs
         );
-        addField
+
+        const surfaceScalarField DByA1f(rAUfs[pair.phase1().index()]*Df);
+        const surfaceScalarField DByA2f(rAUfs[pair.phase2().index()]*Df);
+
+        const volScalarField alpha12(pair.phase1() + pair.phase2());
+        const surfaceScalarField snGradAlpha1By12
         (
-            pair.phase2(),
+            fvc::snGrad
-            "phiFf",
+            (
-           -rAUfs[pair.phase2().index()]*Ff,
+                pair.phase1()/max(alpha12, pair.phase1().residualAlpha())
-            phiFfs
+            )*this->mesh_.magSf()
         );
+        const surfaceScalarField snGradAlpha2By12
+        (
+            fvc::snGrad
+            (
+                pair.phase2()/max(alpha12, pair.phase2().residualAlpha())
+            )*this->mesh_.magSf()
+        );
+
+        addField(pair.phase1(), "phiF", DByA1f*snGradAlpha1By12, phiFfs);
+        addField(pair.phase2(), "phiF", DByA2f*snGradAlpha2By12, phiFfs);
     }
 
     if (this->fillFields_)