driftFluxFoam: Updated in preparation for conversion into a solver module

applications/solvers/multiphase/driftFluxFoam/CourantNo.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
@@ -37,7 +37,7 @@ scalar meanCoNum = 0.0;
     (
         fvc::surfaceSum
         (
-            mag(phi) + mag(fvc::flux(mixture.Udm()))
+            mag(phi) + mag(fvc::flux(relativeVelocity->Udm()))
         )().primitiveField()
     );
 

applications/solvers/multiphase/driftFluxFoam/UEqn.H

@@ -3,8 +3,8 @@
     fvVectorMatrix UEqn
     (
         fvm::ddt(rho, U) + fvm::div(rhoPhi, U)
-      + mixture.MRF().DDt(rho, U)
+      + MRF.DDt(rho, U)
-      + mixture.divTauDm()
+      + relativeVelocity->divDevTau()
       + turbulence->divDevTau(U)
      ==
         fvModels.source(rho, U)

applications/solvers/multiphase/driftFluxFoam/alphaEqnSubCycle.H

@@ -13,7 +13,7 @@
         dimensionedScalar(phi.dimensions(), 0)
     );
 
-    surfaceScalarField phir(fvc::flux(mixture.Udm()));
+    surfaceScalarField phir(fvc::flux(relativeVelocity->Udm()));
 
     if (nAlphaSubCycles > 1)
     {

applications/solvers/multiphase/driftFluxFoam/createFields.H

@@ -35,7 +35,7 @@ volVectorField U
 
 
 Info<< "Reading incompressibleTwoPhaseInteractingMixture\n" << endl;
-incompressibleTwoPhaseInteractingMixture mixture(U, phi, g);
+incompressibleTwoPhaseInteractingMixture mixture(U, phi);
 
 volScalarField& alpha1(mixture.alpha1());
 const volScalarField& rho(mixture.rho());
@@ -55,6 +55,13 @@ surfaceScalarField rhoPhi
 );
 
 
+// Dispersed phase relative velocity model
+autoPtr<relativeVelocityModel> relativeVelocity
+(
+    relativeVelocityModel::New(mixture, mixture, g)
+);
+
+
 // Construct compressible turbulence model
 autoPtr<compressible::momentumTransportModel> turbulence
 (
@@ -95,6 +102,7 @@ mesh.schemes().setFluxRequired(alpha1.name());
 // MULES Correction
 tmp<surfaceScalarField> talphaPhiCorr0;
 
+#include "createMRF.H"
 #include "createFvModels.H"
 #include "createFvConstraints.H"
 
@@ -102,7 +110,7 @@ tmp<surfaceScalarField> talphaPhiCorr0;
 //  used for ddtCorr with MRF
 autoPtr<surfaceVectorField> Uf;
 
-if (mixture.MRF().size())
+if (MRF.size())
 {
     Info<< "Constructing face momentum Uf" << endl;
 

applications/solvers/multiphase/driftFluxFoam/driftFluxFoam.C

@@ -38,6 +38,7 @@ Description
 #include "CMULES.H"
 #include "subCycle.H"
 #include "incompressibleTwoPhaseInteractingMixture.H"
+#include "relativeVelocityModel.H"
 #include "momentumTransportModel.H"
 #include "compressibleMomentumTransportModels.H"
 #include "pimpleControl.H"
@@ -105,7 +106,7 @@ int main(int argc, char *argv[])
 
             #include "alphaEqnSubCycle.H"
 
-            mixture.correct();
+            relativeVelocity->correct();
 
             if (pimple.predictTransport())
             {

applications/solvers/multiphase/driftFluxFoam/incompressibleTwoPhaseInteractingMixture/incompressibleTwoPhaseInteractingMixture.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
@@ -25,8 +25,6 @@ License
 
 #include "incompressibleTwoPhaseInteractingMixture.H"
 #include "mixtureViscosityModel.H"
-#include "relativeVelocityModel.H"
-#include "fvcDiv.H"
 #include "addToRunTimeSelectionTable.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
@@ -43,25 +41,19 @@ Foam::incompressibleTwoPhaseInteractingMixture::
 incompressibleTwoPhaseInteractingMixture
 (
     volVectorField& U,
-    const surfaceScalarField& phi,
+    const surfaceScalarField& phi
-    const uniformDimensionedVectorField& g
 )
 :
     twoPhaseMixture(U.mesh()),
 
     U_(U),
 
-    muModel_(mixtureViscosityModel::New(*this)),
     nucModel_(viscosityModel::New(U.mesh(), phase2Name())),
+    muModel_(mixtureViscosityModel::New(*this)),
 
-    rhod_("rho", dimDensity, muModel_()),
     rhoc_("rho", dimDensity, nucModel_()),
-    dd_
+    rhod_("rho", dimDensity, muModel_()),
-    (
+
-        "d",
-        dimLength,
-        muModel_->lookupOrDefault("d", 0.0)
-    ),
     alphaMax_(lookupOrDefault("alphaMax", 1.0)),
 
     rho_
@@ -78,22 +70,18 @@ incompressibleTwoPhaseInteractingMixture
         dimensionedScalar("rho", dimDensity, 0)
     ),
 
-    mu_
+    nu_
     (
         IOobject
         (
-            "mu",
+            "nu",
             U_.time().name(),
             U_.db()
         ),
         U_.mesh(),
-        dimensionedScalar(dimensionSet(1, -1, -1, 0, 0), 0),
+        dimensionedScalar(dimViscosity, 0),
         calculatedFvPatchScalarField::typeName
-    ),
+    )
-
-    MRF_(U.mesh()),
-
-    UdmModel_(relativeVelocityModel::New(*this, *this, g))
 {
     correct();
 }
@@ -108,13 +96,6 @@ Foam::incompressibleTwoPhaseInteractingMixture::
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
-const Foam::volScalarField&
-Foam::incompressibleTwoPhaseInteractingMixture::alphad() const
-{
-    return alpha1();
-}
-
-
 const Foam::volScalarField&
 Foam::incompressibleTwoPhaseInteractingMixture::alphac() const
 {
@@ -122,24 +103,10 @@ Foam::incompressibleTwoPhaseInteractingMixture::alphac() const
 }
 
 
-const Foam::mixtureViscosityModel&
+const Foam::volScalarField&
-Foam::incompressibleTwoPhaseInteractingMixture::muModel() const
+Foam::incompressibleTwoPhaseInteractingMixture::alphad() const
 {
-    return muModel_();
+    return alpha1();
-}
-
-
-const Foam::viscosityModel&
-Foam::incompressibleTwoPhaseInteractingMixture::nucModel() const
-{
-    return nucModel_();
-}
-
-
-const Foam::dimensionedScalar&
-Foam::incompressibleTwoPhaseInteractingMixture::rhod() const
-{
-    return rhod_;
 }
 
 
@@ -151,9 +118,9 @@ Foam::incompressibleTwoPhaseInteractingMixture::rhoc() const
 
 
 const Foam::dimensionedScalar&
-Foam::incompressibleTwoPhaseInteractingMixture::dd() const
+Foam::incompressibleTwoPhaseInteractingMixture::rhod() const
 {
-    return dd_;
+    return rhod_;
 }
 
 
@@ -177,60 +144,24 @@ Foam::incompressibleTwoPhaseInteractingMixture::rho() const
 }
 
 
-Foam::tmp<Foam::volScalarField>
-Foam::incompressibleTwoPhaseInteractingMixture::mu() const
-{
-    return mu_;
-}
-
-
-Foam::tmp<Foam::scalarField>
-Foam::incompressibleTwoPhaseInteractingMixture::mu(const label patchi) const
-{
-    return mu_.boundaryField()[patchi];
-}
-
-
 Foam::tmp<Foam::volScalarField>
 Foam::incompressibleTwoPhaseInteractingMixture::nu() const
 {
-    return mu_/rho_;
+    return nu_;
 }
 
 
 Foam::tmp<Foam::scalarField>
 Foam::incompressibleTwoPhaseInteractingMixture::nu(const label patchi) const
 {
-    return mu_.boundaryField()[patchi]/rho_.boundaryField()[patchi];
+    return nu_.boundaryField()[patchi];
-}
-
-
-const Foam::IOMRFZoneList&
-Foam::incompressibleTwoPhaseInteractingMixture::MRF() const
-{
-    return MRF_;
-}
-
-
-const Foam::volVectorField&
-Foam::incompressibleTwoPhaseInteractingMixture::Udm() const
-{
-    return UdmModel_->Udm();
-}
-
-
-Foam::tmp<Foam::volVectorField>
-Foam::incompressibleTwoPhaseInteractingMixture::divTauDm() const
-{
-    return fvc::div(UdmModel_->tauDm());
 }
 
 
 void Foam::incompressibleTwoPhaseInteractingMixture::correct()
 {
     rho_ = alpha1()*rhod_ + alpha2()*rhoc_;
-    mu_ = muModel_->mu(rhoc_*nucModel_->nu(), U_);
+    nu_ = muModel_->mu(rhoc_*nucModel_->nu(), U_)/rho_;
-    UdmModel_->correct();
 }
 
 
@@ -240,15 +171,8 @@ bool Foam::incompressibleTwoPhaseInteractingMixture::read()
     {
         if (muModel_->read() || nucModel_->read())
         {
-            muModel_->lookup("rho") >> rhod_;
             nucModel_->lookup("rho") >> rhoc_;
-
+            muModel_->lookup("rho") >> rhod_;
-            dd_ = dimensionedScalar
-            (
-                "d",
-                dimLength,
-                muModel_->lookupOrDefault("d", 0)
-            );
 
             alphaMax_ = muModel_->lookupOrDefault("alphaMax", 1.0);
 

applications/solvers/multiphase/driftFluxFoam/incompressibleTwoPhaseInteractingMixture/incompressibleTwoPhaseInteractingMixture.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
@@ -39,8 +39,6 @@ SourceFiles
 
 #include "twoPhaseMixture.H"
 #include "viscosityModel.H"
-#include "uniformDimensionedFields.H"
-#include "IOMRFZoneList.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -48,7 +46,6 @@ namespace Foam
 {
 
 class mixtureViscosityModel;
-class relativeVelocityModel;
 
 /*---------------------------------------------------------------------------*\
           Class incompressibleTwoPhaseInteractingMixture Declaration
@@ -63,27 +60,26 @@ class incompressibleTwoPhaseInteractingMixture
 
         volVectorField& U_;
 
-        autoPtr<mixtureViscosityModel> muModel_;
+        //- Continuous phase viscosity model
         autoPtr<viscosityModel> nucModel_;
 
-        dimensionedScalar rhod_;
+        //- Mixture viscosity model
+        autoPtr<mixtureViscosityModel> muModel_;
+
+        //- Continuous phase density
         dimensionedScalar rhoc_;
 
-        //- Optional diameter of the dispersed phase particles
+        //- Dispersed phase density
-        dimensionedScalar dd_;
+        dimensionedScalar rhod_;
 
         //- Optional maximum dispersed phase-fraction (e.g. packing limit)
         scalar alphaMax_;
 
+        //- Mixture density
         volScalarField rho_;
 
-        volScalarField mu_;
+        //- Mixture kinematic viscosity
-
+        volScalarField nu_;
-        //- Optional MRF zones
-        IOMRFZoneList MRF_;
-
-        //- Dispersed phase relative velocity model
-        autoPtr<relativeVelocityModel> UdmModel_;
 
 
 public:
@@ -97,8 +93,7 @@ public:
         incompressibleTwoPhaseInteractingMixture
         (
             volVectorField& U,
-            const surfaceScalarField& phi,
+            const surfaceScalarField& phi
-            const uniformDimensionedVectorField& g
         );
 
 
@@ -114,20 +109,11 @@ public:
         //- Return const-access to the dispersed phase-fraction
         const volScalarField& alphad() const;
 
-        //- Return const-access to the mixture viscosityModel
-        const mixtureViscosityModel& muModel() const;
-
-        //- Return const-access to the continuous-phase viscosityModel
-        const viscosityModel& nucModel() const;
-
-        //- Return const-access to the dispersed-phase density
-        const dimensionedScalar& rhod() const;
-
         //- Return const-access to continuous-phase density
         const dimensionedScalar& rhoc() const;
 
-        //- Return the diameter of the dispersed-phase particles
+        //- Return const-access to the dispersed-phase density
-        const dimensionedScalar& dd() const;
+        const dimensionedScalar& rhod() const;
 
         //- Optional maximum phase-fraction (e.g. packing limit)
         //  Defaults to 1
@@ -139,27 +125,12 @@ public:
         //- Return the mixture density
         virtual const volScalarField& rho() const;
 
-        //- Return the dynamic mixture viscosity
-        tmp<volScalarField> mu() const;
-
-        //- Return the dynamic mixture viscosity for patch
-        virtual tmp<scalarField> mu(const label patchi) const;
-
         //- Return the mixture viscosity
         virtual tmp<volScalarField> nu() const;
 
         //- Return the mixture viscosity for patch
         virtual tmp<scalarField> nu(const label patchi) const;
 
-        //- Return MRF zones
-        const IOMRFZoneList& MRF() const;
-
-        //- Return the diffusion velocity of the dispersed phase
-        const volVectorField& Udm() const;
-
-        //- Return the div stress tensor due to the phase transport
-        tmp<volVectorField> divTauDm() const;
-
         //- Correct the laminar viscosity
         virtual void correct();
 

applications/solvers/multiphase/driftFluxFoam/pEqn.H

@@ -8,7 +8,7 @@
         fvc::flux(HbyA)
       + fvc::interpolate(rho*rAU)*fvc::ddtCorr(U, phi, Uf)
     );
-    mixture.MRF().makeRelative(phiHbyA);
+    MRF.makeRelative(phiHbyA);
     adjustPhi(phiHbyA, U, p_rgh);
 
     surfaceScalarField phig
@@ -21,7 +21,7 @@
     phiHbyA += phig;
 
     // Update the pressure BCs to ensure flux consistency
-    constrainPressure(p_rgh, U, phiHbyA, rAUf, mixture.MRF());
+    constrainPressure(p_rgh, U, phiHbyA, rAUf, MRF);
 
     while (pimple.correctNonOrthogonal())
     {
@@ -67,5 +67,5 @@
     }
 
     // Correct Uf for MRF
-    fvc::correctUf(Uf, U, phi, mixture.MRF());
+    fvc::correctUf(Uf, U, phi, MRF);
 }

applications/solvers/multiphase/driftFluxFoam/relativeVelocityModels/relativeVelocityModel/relativeVelocityModel.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
@@ -28,6 +28,7 @@ License
 #include "slipFvPatchFields.H"
 #include "partialSlipFvPatchFields.H"
 #include "fvcGrad.H"
+#include "fvcDiv.H"
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
@@ -177,4 +178,10 @@ Foam::tmp<Foam::volSymmTensorField> Foam::relativeVelocityModel::tauDm() const
 }
 
 
+Foam::tmp<Foam::volVectorField> Foam::relativeVelocityModel::divDevTau() const
+{
+    return fvc::div(tauDm());
+}
+
+
 // ************************************************************************* //

applications/solvers/multiphase/driftFluxFoam/relativeVelocityModels/relativeVelocityModel/relativeVelocityModel.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
@@ -137,6 +137,9 @@ public:
         //- Return the stress tensor due to the phase transport
         tmp<volSymmTensorField> tauDm() const;
 
+        //- Return the div stress tensor due to the phase transport
+        tmp<volVectorField> divDevTau() const;
+
         //- Update the diffusion velocity
         virtual void correct() = 0;
 

src/finiteVolume/fvMatrices/solvers/MULES/CMULESTemplates.C

@@ -2,7 +2,7 @@
   =========                 |
   \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
    \\    /   O peration     | Website:  https://openfoam.org
-    \\  /    A nd           | Copyright (C) 2013-2022 OpenFOAM Foundation
+    \\  /    A nd           | Copyright (C) 2013-2023 OpenFOAM Foundation
      \\/     M anipulation  |
 -------------------------------------------------------------------------------
 License
@@ -361,8 +361,16 @@ void Foam::MULES::limiterCorr
         }
     }
 
+    if (extremaCoeff > 0)
+    {
         psiMaxn = min(psiMaxn + extremaCoeff*(psiMax - psiMin), psiMax);
         psiMinn = max(psiMinn - extremaCoeff*(psiMax - psiMin), psiMin);
+    }
+    else
+    {
+        psiMaxn = min(psiMaxn, psiMax);
+        psiMinn = max(psiMinn, psiMin);
+    }
 
     if (smoothLimiter > small)
     {

tutorials/multiphase/driftFluxFoam/RAS/dahl/system/fvSolution

@@ -25,8 +25,6 @@ solvers
         nLimiterIter    3;
         alphaApplyPrevCorr  yes;
 
-        boundaryExtremaCoeff 1;
-
         solver          smoothSolver;
         smoother        symGaussSeidel;
         tolerance       1e-6;