functionObjects: forces: Added support for Euler-Euler multiphase

src/functionObjects/forces/Make/options

@@ -5,6 +5,8 @@ EXE_INC = \
     -I$(LIB_SRC)/MomentumTransportModels/momentumTransportModels/lnInclude \
     -I$(LIB_SRC)/MomentumTransportModels/incompressible/lnInclude \
     -I$(LIB_SRC)/MomentumTransportModels/compressible/lnInclude \
+    -I$(LIB_SRC)/MomentumTransportModels/phaseIncompressible/lnInclude \
+    -I$(LIB_SRC)/MomentumTransportModels/phaseCompressible/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude
 
@@ -14,6 +16,8 @@ LIB_LIBS = \
     -lmomentumTransportModels \
     -lincompressibleMomentumTransportModels \
     -lcompressibleMomentumTransportModels \
+    -lphaseIncompressibleMomentumTransportModels \
+    -lphaseCompressibleMomentumTransportModels \
     -lspecie \
     -lfileFormats \
     -lfiniteVolume \

src/functionObjects/forces/forces/forces.C

@@ -28,6 +28,8 @@ License
 #include "porosityModel.H"
 #include "kinematicMomentumTransportModel.H"
 #include "dynamicMomentumTransportModel.H"
+#include "phaseKinematicMomentumTransportModel.H"
+#include "phaseDynamicMomentumTransportModel.H"
 #include "fluidThermo.H"
 #include "addToRunTimeSelectionTable.H"
 
@@ -224,23 +226,45 @@ void Foam::functionObjects::forces::initialise()
 Foam::tmp<Foam::volSymmTensorField>
 Foam::functionObjects::forces::devTau() const
 {
-    typedef compressible::momentumTransportModel cmpModel;
     typedef incompressible::momentumTransportModel icoModel;
+    typedef compressible::momentumTransportModel cmpModel;
+    typedef phaseIncompressible::momentumTransportModel phaseIcoModel;
+    typedef phaseCompressible::momentumTransportModel phaseCmpModel;
 
-    if (obr_.foundObject<cmpModel>(momentumTransportModel::typeName))
-    {
-        const cmpModel& model =
-            obr_.lookupObject<cmpModel>(momentumTransportModel::typeName);
+    const word& modelName = momentumTransportModel::typeName;
+    const word phaseModelName =
+        phaseName_ == word::null
+      ? word::null
+      : IOobject::groupName(momentumTransportModel::typeName, phaseName_);
 
-        return model.devTau();
-    }
-    else if (obr_.foundObject<icoModel>(momentumTransportModel::typeName))
+    if (obr_.foundObject<icoModel>(modelName))
     {
         const incompressible::momentumTransportModel& model =
-            obr_.lookupObject<icoModel>(momentumTransportModel::typeName);
+            obr_.lookupObject<icoModel>(modelName);
+
+        return alpha()*rho()*model.devSigma();
+    }
+    else if (obr_.foundObject<cmpModel>(modelName))
+    {
+        const cmpModel& model =
+            obr_.lookupObject<cmpModel>(modelName);
+
+        return alpha()*model.devTau();
+    }
+    else if (obr_.foundObject<phaseIcoModel>(phaseModelName))
+    {
+        const phaseIcoModel& model =
+            obr_.lookupObject<phaseIcoModel>(phaseModelName);
 
         return rho()*model.devSigma();
     }
+    else if (obr_.foundObject<phaseCmpModel>(phaseModelName))
+    {
+        const phaseCmpModel& model =
+            obr_.lookupObject<phaseCmpModel>(phaseModelName);
+
+        return model.devTau();
+    }
     else if (obr_.foundObject<dictionary>("transportProperties"))
     {
         // Legacy support for icoFoam
@@ -272,25 +296,49 @@ Foam::functionObjects::forces::devTau() const
 
 Foam::tmp<Foam::volScalarField> Foam::functionObjects::forces::mu() const
 {
-    if (obr_.foundObject<fluidThermo>(basicThermo::dictName))
-    {
-        const fluidThermo& thermo =
-             obr_.lookupObject<fluidThermo>(basicThermo::dictName);
+    typedef incompressible::momentumTransportModel icoModel;
+    typedef compressible::momentumTransportModel cmpModel;
+    typedef phaseIncompressible::momentumTransportModel phaseIcoModel;
+    typedef phaseCompressible::momentumTransportModel phaseCmpModel;
 
-        return thermo.mu();
+    const word& modelName = momentumTransportModel::typeName;
+    const word phaseModelName =
+        phaseName_ == word::null
+      ? word::null
+      : IOobject::groupName(momentumTransportModel::typeName, phaseName_);
+
+    if (obr_.foundObject<icoModel>(modelName))
+    {
+        const incompressible::momentumTransportModel& model =
+            obr_.lookupObject<icoModel>(modelName);
+
+        return rho()*model.transport().nu();
     }
-    else if
-    (
-        obr_.foundObject<kinematicTransportModel>("transportProperties")
-    )
+    else if (obr_.foundObject<cmpModel>(modelName))
     {
-        const kinematicTransportModel& laminarT =
-            obr_.lookupObject<kinematicTransportModel>("transportProperties");
+        const cmpModel& model =
+            obr_.lookupObject<cmpModel>(modelName);
 
-        return rho()*laminarT.nu();
+        return model.transport().mu();
+    }
+    else if (obr_.foundObject<phaseIcoModel>(phaseModelName))
+    {
+        const phaseIcoModel& model =
+            obr_.lookupObject<phaseIcoModel>(phaseModelName);
+
+        return rho()*model.transport().nu();
+    }
+    else if (obr_.foundObject<phaseCmpModel>(phaseModelName))
+    {
+        const phaseCmpModel& model =
+            obr_.lookupObject<phaseCmpModel>(phaseModelName);
+
+        return model.transport().mu();
     }
     else if (obr_.foundObject<dictionary>("transportProperties"))
     {
+        // Legacy support for icoFoam
+
         const dictionary& transportProperties =
              obr_.lookupObject<dictionary>("transportProperties");
 
@@ -352,27 +400,45 @@ Foam::scalar Foam::functionObjects::forces::rho(const volScalarField& p) const
 }
 
 
-Foam::tmp<Foam::vectorField> Foam::functionObjects::forces::phaseFilter
-(
-    const tmp<vectorField>& tF,
-    const label patchi
-) const
+Foam::tmp<Foam::volScalarField> Foam::functionObjects::forces::alpha() const
 {
-    if (phaseName_ != word::null)
+    if (phaseName_ == word::null)
     {
-        const volScalarField& alpha
+        return volScalarField::New
         (
-            obr_.lookupObject<volScalarField>
-            (
-                IOobject::groupName("alpha", phaseName_)
-            )
+            "alpha",
+            mesh_,
+            dimensionedScalar(dimless, 1)
         );
-
-        return alpha.boundaryField()[patchi]*tF;
     }
     else
     {
-        return tF;
+        return obr_.lookupObject<volScalarField>
+        (
+            IOobject::groupName("alpha", phaseName_)
+        );
+    }
+}
+
+
+Foam::tmp<Foam::scalarField> Foam::functionObjects::forces::alpha
+(
+    const label patchi
+) const
+{
+    if (phaseName_ == word::null)
+    {
+        return tmp<scalarField>
+        (
+            new scalarField(mesh_.boundary()[patchi].size(), 1)
+        );
+    }
+    else
+    {
+        return obr_.lookupObject<volScalarField>
+        (
+            IOobject::groupName("alpha", phaseName_)
+        ).boundaryField()[patchi];
     }
 }
 
@@ -854,8 +920,8 @@ void Foam::functionObjects::forces::calcForcesMoment()
             mesh_.Sf().boundaryField();
 
         tmp<volSymmTensorField> tdevTau = devTau();
-        const volSymmTensorField::Boundary& devTaub
-            = tdevTau().boundaryField();
+        const volSymmTensorField::Boundary& devTaub =
+            tdevTau().boundaryField();
 
         // Scale pRef by density for incompressible simulations
         const scalar pRef = pRef_/rho(p);
@@ -871,17 +937,13 @@ void Foam::functionObjects::forces::calcForcesMoment()
 
             const vectorField fN
             (
-                phaseFilter
-                (
-                    rho(p)*Sfb[patchi]*(p.boundaryField()[patchi] - pRef),
-                    patchi
-                )
+                alpha(patchi)
+               *rho(p)
+               *Sfb[patchi]
+               *(p.boundaryField()[patchi] - pRef)
             );
 
-            const vectorField fT
-            (
-                phaseFilter(Sfb[patchi] & devTaub[patchi], patchi)
-            );
+            const vectorField fT(Sfb[patchi] & devTaub[patchi]);
 
             const vectorField fP(Md.size(), Zero);
 

src/functionObjects/forces/forces/forces.H

@@ -251,11 +251,11 @@ protected:
         //  otherwise return 1
         scalar rho(const volScalarField& p) const;
 
-        tmp<vectorField> phaseFilter
-        (
-            const tmp<vectorField>& F,
-            const label patchi
-        ) const;
+        //- Get the volume fraction field
+        tmp<volScalarField> alpha() const;
+
+        //- Get the volume fraction field on a patch
+        tmp<scalarField> alpha(const label patchi) const;
 
         //- Accumulate bin data
         void applyBins