Merge pull request #125 from ParticulateFlow/release

Release 21.11

Author: danielque

applications/solvers/cfdemSolverMultiphase/multiphaseMixture/multiphaseMixture.C

@@ -55,14 +55,21 @@ Foam::multiphaseMixture::calcNu() const
 {
     PtrDictionary<phase>::const_iterator iter = phases_.begin();
 
+    // 1/nu
+    tmp<volScalarField> tnuInv = iter()/iter().nu();
+    volScalarField& nuInv = tnuInv.ref();
+
+    // nu
     tmp<volScalarField> tnu = iter()*iter().nu();
     volScalarField& nu = tnu.ref();
 
     for (++iter; iter != phases_.end(); ++iter)
     {
-        nu += iter()*iter().nu();
+        nuInv += iter()/iter().nu();
     }
 
+    nu = 1/nuInv;
+    
     return tnu;
 }
 

applications/solvers/cfdemSolverMultiphaseScalar/Allwclean

@@ -0,0 +1,8 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # Run from this directory
+set -x
+
+wclean libso multiphaseMixture
+wclean
+
+#------------------------------------------------------------------------------

applications/solvers/cfdemSolverMultiphaseScalar/Allwmake

@@ -0,0 +1,12 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # Run from this directory
+
+# Parse arguments for library compilation
+targetType=libso
+. $WM_PROJECT_DIR/wmake/scripts/AllwmakeParseArguments
+set -x
+
+wmake $targetType multiphaseMixture
+wmake
+
+#------------------------------------------------------------------------------

applications/solvers/cfdemSolverMultiphaseScalar/CEqn.H

@@ -0,0 +1,22 @@
+        // get mixture properties
+        Cs = mixture.Cs();
+        diffusionCorrection = mixture.diffusionCorrection();
+        Deff = particleCloud.diffCoeffM().diffCoeff();
+
+        // get scalar source from DEM
+        particleCloud.massContributions(Sm);
+        particleCloud.massCoefficients(Smi);
+
+        fvScalarMatrix CEqn
+        (
+            fvm::ddt(voidfraction,C)
+          + fvm::div(phi,C)
+          - fvm::laplacian(Deff*voidfraction,C)
+          + fvm::div(fvc::interpolate(Deff*voidfraction)*diffusionCorrection*mesh.magSf(), C)
+          ==
+            Sm + fvm::Sp(Smi,C)
+        );
+
+        CEqn.relax();
+        fvOptions.constrain(CEqn);
+        CEqn.solve();

applications/solvers/cfdemSolverMultiphaseScalar/EEqn.H

@@ -0,0 +1,22 @@
+        // get mixture properties
+        Cp = mixture.Cp();
+        kf = mixture.kf();
+
+        // get scalar source from DEM
+        particleCloud.energyContributions(Qsource);
+        particleCloud.energyCoefficients(QCoeff);
+
+        fvScalarMatrix EEqn
+        (
+            rho*Cp*(fvm::ddt(voidfraction,T)
+          + fvm::div(phi,T))
+          - fvm::laplacian(thCond*voidfraction,T)
+          ==
+            Qsource + fvm::Sp(QCoeff,T)
+        );
+
+
+        EEqn.relax();
+        fvOptions.constrain(EEqn);
+        EEqn.solve();
+

applications/solvers/cfdemSolverMultiphaseScalar/Make/files

@@ -0,0 +1,3 @@
+cfdemSolverMultiphaseScalar.C
+
+EXE = $(CFDEM_APP_DIR)/cfdemSolverMultiphaseScalar

applications/solvers/cfdemSolverMultiphaseScalar/Make/options

@@ -0,0 +1,35 @@
+FOAM_VERSION_MAJOR := $(word 1,$(subst ., ,$(WM_PROJECT_VERSION)))
+PFLAGS+= -DOPENFOAM_VERSION_MAJOR=$(FOAM_VERSION_MAJOR)
+
+include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
+
+EXE_INC = \
+     $(PFLAGS) \
+    -I$(CFDEM_OFVERSION_DIR) \
+    -ImultiphaseMixture/lnInclude \
+    -I$(LIB_SRC)/transportModels \
+    -I$(LIB_SRC)/transportModels/incompressible/lnInclude \
+    -I$(LIB_SRC)/transportModels/interfaceProperties/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/lnInclude \
+    -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/cfdTools \
+    -Wno-deprecated-copy
+
+EXE_LIBS = \
+    -L$(CFDEM_LIB_DIR)\
+    -lcfdemMultiphaseInterFoamScalar \
+    -linterfaceProperties \
+    -lincompressibleTransportModels \
+    -lturbulenceModels \
+    -lincompressibleTurbulenceModels \
+    -lfiniteVolume \
+    -lfvOptions \
+    -lmeshTools \
+    -lsampling \
+    -l$(CFDEM_LIB_NAME) \
+     $(CFDEM_ADD_LIB_PATHS) \
+     $(CFDEM_ADD_LIBS)

applications/solvers/cfdemSolverMultiphaseScalar/UEqn.H

@@ -0,0 +1,61 @@
+const surfaceScalarField& rhoPhi(mixture.rhoPhi());
+
+volScalarField muEff = rho*(turbulence->nu() + turbulence->nut());
+
+if (modelType == "A")
+    muEff *= voidfraction;
+
+fvVectorMatrix UEqn
+(
+    fvm::ddt(rhoEps, U) - fvm::Sp(fvc::ddt(rhoEps),U)
+  + fvm::div(rhoPhi, U) - fvm::Sp(fvc::div(rhoPhi),U)
+  //+ particleCloud.divVoidfractionTau(U, voidfraction)
+  - fvm::laplacian(muEff, U) - fvc::div(muEff*dev2(fvc::grad(U)().T()))
+  ==
+    fvOptions(rho, U)
+  - fvm::Sp(Ksl,U)
+);
+
+UEqn.relax();
+
+fvOptions.constrain(UEqn);
+
+if (pimple.momentumPredictor() && (modelType=="B" || modelType=="Bfull"))
+{
+    solve
+    (
+        UEqn
+     ==
+        fvc::reconstruct
+        (
+            (- ghf*fvc::snGrad(rho) - fvc::snGrad(p_rgh)) * mesh.magSf()
+        )
+     +
+        fvc::reconstruct
+        (
+            mixture.surfaceTensionForce() * mesh.magSf()
+        ) * voidfraction
+     + Ksl*Us
+    );
+
+    fvOptions.correct(U);
+}
+else if (pimple.momentumPredictor())
+{
+    solve
+    (
+        UEqn
+     ==
+        fvc::reconstruct
+        (
+            (
+                mixture.surfaceTensionForce()
+              - ghf*fvc::snGrad(rho)
+              - fvc::snGrad(p_rgh)
+            ) * mesh.magSf()
+        ) * voidfraction
+     + Ksl*Us
+    );
+
+    fvOptions.correct(U);
+}

applications/solvers/cfdemSolverMultiphaseScalar/additionalChecks.H

@@ -0,0 +1,17 @@
+// Additional solver-specific checks
+
+// Useful if one wants to e.g. initialize floating particles using the Archimedes model
+if (particleCloud.couplingProperties().found("unrestrictedForceModelSelection"))
+{
+    Warning << "Using unrestrictedForceModelSelection, results may be incorrect!" << endl;
+} else
+{
+    #include "checkModelType.H"
+}
+
+word modelType = particleCloud.modelType();
+
+if(!particleCloud.couplingProperties().found("useDDTvoidfraction"))
+{
+    Warning << "Suppressing ddt(voidfraction) is not recommended with this solver as it may generate incorrect results!" << endl;
+}

applications/solvers/cfdemSolverMultiphaseScalar/alphaCourantNo.H

@@ -0,0 +1,21 @@
+scalar alphaCoNum = 0.0;
+scalar meanAlphaCoNum = 0.0;
+
+if (mesh.nInternalFaces())
+{
+    scalarField sumPhi
+    (
+        mixture.nearInterface()().primitiveField()
+       *fvc::surfaceSum(mag(phi))().primitiveField()
+    );
+
+    alphaCoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue();
+
+    meanAlphaCoNum =
+        0.5*(gSum(sumPhi)/gSum(mesh.V().field()))*runTime.deltaTValue();
+}
+
+Info<< "Interface Courant Number mean: " << meanAlphaCoNum
+    << " max: " << alphaCoNum << endl;
+
+// ************************************************************************* //