reacting*EulerFoam/.../turbulentDispersionModels: Corrected for multiphase case

In two phases, the turbulent dispersion force is modelled for a phase
pair as follows:

    F12 = D12*grad(alpha1)

Where D12 is the turbulent dispersion coefficient between phases 1 and
2. This force is calculated equivalently whichever phase is chosen to be
phase 1 because the volume fractions are related by alpha1 = 1 - alpha2.
This means that F12 == - F21; i.e., the force in one phase equals the
reaction in the other.

In multiple phases, however, a force calculated in this way is no longer
consistent between phases, because the relationship between the volume
fractions no longer applies. The following form has been chosen instead.

    F12 = D12*grad(alpha1/(alpha1 + alpha2))

I.e., rather than using the gradient of a phase directly, we use the
gradient of the phase within the two-phase sub-system associated with
the pair.

This reduces to the two-phase case above, and the models available in
the literature that are explicitly formulated for multiple phases can
also be expressed in this form.

Based on a patch contributed by Institute of Fluid Dynamics, Helmholtz-Zentrum Dresden -
Rossendorf (HZDR) and VTT Technical Research Centre of Finland Ltd.

applications/solvers/multiphase/reactingEulerFoam/interfacialModels/turbulentDispersionModels/Burns/Burns.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2014-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2014-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -56,17 +56,7 @@ Foam::turbulentDispersionModels::Burns::Burns
 )
 :
     turbulentDispersionModel(dict, pair),
-    sigma_("sigma", dimless, dict),
-    residualAlpha_
-    (
-        "residualAlpha",
-        dimless,
-        dict.lookupOrDefault<scalar>
-        (
-            "residualAlpha",
-            pair_.dispersed().residualAlpha().value()
-        )
-    )
+    sigma_("sigma", dimless, dict)
 {}
 
 
@@ -89,19 +79,14 @@ Foam::turbulentDispersionModels::Burns::D() const
         );
 
     return
-        0.75
-       *drag.CdRe()
-       *pair_.continuous().nu()
+        drag.Ki()
        *continuousTurbulence().nut()
-       /(
-            sigma_
-           *sqr(pair_.dispersed().d())
-        )
-       *pair_.continuous().rho()
+       /sigma_
        *pair_.dispersed()
-       *(
-           1.0/max(pair_.dispersed(), residualAlpha_)
-         + 1.0/max(pair_.continuous(), residualAlpha_)
+       *sqr(pair_.dispersed() + pair_.continuous())
+       /(
+            max(pair_.dispersed(), pair_.dispersed().residualAlpha())
+           *max(pair_.continuous(), pair_.continuous().residualAlpha())
         );
 }
 

applications/solvers/multiphase/reactingEulerFoam/interfacialModels/turbulentDispersionModels/Burns/Burns.H

@@ -28,13 +28,6 @@ Description
     Turbulent dispersion model of Burns et al.
 
     References:
-    \verbatim
-        Otromke, M. (2013).
-        Implementation and Comparison of Correlations for interfacial Forces in
-        a Gas-Liquid System within an Euler-Euler Framework.
-        PhD Thesis.
-    \endverbatim
-
     \verbatim
         Burns, A. D., Frank, T., Hamill, I., & Shi, J. M. (2004, May).
         The Favre averaged drag model for turbulent dispersion in Eulerian
@@ -76,9 +69,6 @@ class Burns
         //- Schmidt number
         const dimensionedScalar sigma_;
 
-        //- Residual phase fraction
-        const dimensionedScalar residualAlpha_;
-
 
 public:
 

applications/solvers/multiphase/reactingEulerFoam/interfacialModels/turbulentDispersionModels/Gosman/Gosman.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2014-2018 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2014-2019 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -81,16 +81,10 @@ Foam::turbulentDispersionModels::Gosman::D() const
         );
 
     return
-        0.75
-       *drag.CdRe()
-       *pair_.dispersed()
-       *pair_.continuous().nu()
+        drag.Ki()
        *continuousTurbulence().nut()
-       /(
-            sigma_
-           *sqr(pair_.dispersed().d())
-        )
-       *pair_.continuous().rho();
+       /sigma_
+       *pair_.dispersed();
 }
 
 

applications/solvers/multiphase/reactingEulerFoam/interfacialModels/turbulentDispersionModels/LopezDeBertodano/LopezDeBertodano.H

@@ -28,19 +28,11 @@ Description
     Lopez de Bertodano (1992) turbulent dispersion model.
 
     \verbatim
-        Lopez, D. B. M. (1993).
+        Lopez, D. B. M. (1992).
         Turbulent bubbly two-phase flow in a triangular duct.
         PhD Thesis, Rensselaer Polytechnic Institution.
     \endverbatim
 
-    \verbatim
-        Burns, A. D., Frank, T., Hamill, I., & Shi, J. M. (2004, May).
-        The Favre averaged drag model for turbulent dispersion in Eulerian
-        multi-phase flows.
-        In 5th international conference on multiphase flow,
-        ICMF (Vol. 4, pp. 1-17).
-    \endverbatim
-
 SourceFiles
     LopezDeBertodano.C
 

applications/solvers/multiphase/reactingEulerFoam/interfacialModels/turbulentDispersionModels/constantTurbulentDispersionCoefficient/constantTurbulentDispersionCoefficient.H

@@ -48,7 +48,7 @@ namespace turbulentDispersionModels
 {
 
 /*---------------------------------------------------------------------------*\
-              Class constantTurbulentDispersionCoefficient Declaration
+            Class constantTurbulentDispersionCoefficient Declaration
 \*---------------------------------------------------------------------------*/
 
 class constantTurbulentDispersionCoefficient

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

@@ -82,7 +82,17 @@ Foam::turbulentDispersionModel::continuousTurbulence() const
 Foam::tmp<Foam::volVectorField>
 Foam::turbulentDispersionModel::F() const
 {
-    return D()*fvc::grad(pair_.dispersed());
+    return
+        D()
+       *fvc::grad
+        (
+            pair_.dispersed()
+           /max
+            (
+                pair_.dispersed() + pair_.continuous(),
+                pair_.dispersed().residualAlpha()
+            )
+        );
 }
 
 
@@ -90,8 +100,19 @@ Foam::tmp<Foam::surfaceScalarField>
 Foam::turbulentDispersionModel::Ff() const
 {
     return
-        fvc::interpolate(D())*fvc::snGrad(pair_.dispersed())
-       *pair_.phase1().mesh().magSf();
+    pair_.phase1().mesh().magSf()
+   *(
+        fvc::interpolate(D())
+       *fvc::snGrad
+        (
+            pair_.dispersed()
+           /max
+            (
+                pair_.dispersed() + pair_.continuous(),
+                pair_.dispersed().residualAlpha()
+            )
+        )
+    );
 }
 
 

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

@@ -619,17 +619,33 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::phiFs
             fvc::interpolate(rAUs[pair.phase2().index()]*D)
         );
 
-        const surfaceScalarField snGradAlpha1
+        const surfaceScalarField snGradAlpha1By12
         (
-            fvc::snGrad(pair.phase1())*this->mesh_.magSf()
+            fvc::snGrad
+            (
+                pair.phase1()
+               /max
+                (
+                    pair.phase1() + pair.phase2(),
+                    pair.phase1().residualAlpha()
+                )
+            )*this->mesh_.magSf()
         );
-        const surfaceScalarField snGradAlpha2
+        const surfaceScalarField snGradAlpha2By12
         (
-            fvc::snGrad(pair.phase2())*this->mesh_.magSf()
+            fvc::snGrad
+            (
+                pair.phase2()
+               /max
+                (
+                    pair.phase1() + pair.phase2(),
+                    pair.phase2().residualAlpha()
+                )
+            )*this->mesh_.magSf()
         );
 
-        addField(pair.phase1(), "phiF", DByA1f*snGradAlpha1, phiFs);
-        addField(pair.phase2(), "phiF", DByA2f*snGradAlpha2, phiFs);
+        addField(pair.phase1(), "phiF", DByA1f*snGradAlpha1By12, phiFs);
+        addField(pair.phase2(), "phiF", DByA2f*snGradAlpha2By12, phiFs);
     }
 
     if (this->fillFields_)
@@ -745,7 +761,7 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::phiFfs
         addField(phase, "phiFf", pPrimeByAf*snGradAlpha1, phiFfs);
     }
 
-    // Add the turbulent dispersion force and phase pressure
+    // Add the turbulent dispersion force
     forAllConstIter
     (
         turbulentDispersionModelTable,
@@ -757,8 +773,6 @@ Foam::MomentumTransferPhaseSystem<BasePhaseSystem>::phiFfs
         const phasePair&
             pair(this->phasePairs_[turbulentDispersionModelIter.key()]);
 
-        const volScalarField D(turbulentDispersionModelIter()->D());
-
         addField
         (
             pair.phase1(),