Merge pull request #80 from ParticulateFlow/release

Release 19.02

.gitignore

@@ -7,5 +7,6 @@ log.*
 *~
 
 **/linux*Gcc*/
+**/.vscode
 
 lnInclude

applications/solvers/cfdemSolverMultiphase/Allwclean

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

applications/solvers/cfdemSolverMultiphase/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/cfdemSolverMultiphase/Make/files

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

applications/solvers/cfdemSolverMultiphase/Make/options

@@ -0,0 +1,30 @@
+include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
+
+EXE_INC = \
+    -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 \
+
+EXE_LIBS = \
+    -L$(CFDEM_LIB_DIR)\
+    -lcfdemMultiphaseInterFoam \
+    -linterfaceProperties \
+    -lincompressibleTransportModels \
+    -lturbulenceModels \
+    -lincompressibleTurbulenceModels \
+    -lfiniteVolume \
+    -lfvOptions \
+    -lmeshTools \
+    -lsampling \
+    -l$(CFDEM_LIB_NAME) \
+     $(CFDEM_ADD_LIB_PATHS) \
+     $(CFDEM_ADD_LIBS)

applications/solvers/cfdemSolverMultiphase/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/cfdemSolverMultiphase/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/cfdemSolverMultiphase/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;
+
+// ************************************************************************* //

applications/solvers/cfdemSolverMultiphase/cfdemSolverMultiphase.C

@@ -0,0 +1,148 @@
+/*---------------------------------------------------------------------------*\
+License
+
+    This is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    This code is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this code.  If not, see <http://www.gnu.org/licenses/>.
+
+    Copyright (C) 2018- Mathias Vångö, JKU Linz, Austria
+
+Application
+    cfdemSolverMultiphase
+
+Description
+    CFD-DEM solver for n incompressible fluids which captures the interfaces and
+    includes surface-tension and contact-angle effects for each phase. It is based
+    on the OpenFOAM(R)-4.x solver multiphaseInterFoam but extended to incorporate
+    DEM functionalities from the open-source DEM code LIGGGHTS.
+
+    Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "multiphaseMixture.H"
+#include "turbulentTransportModel.H"
+#include "pimpleControl.H"
+#include "fvOptions.H"
+#include "CorrectPhi.H"
+
+#include "cfdemCloud.H"
+#include "implicitCouple.H"
+#include "clockModel.H"
+#include "smoothingModel.H"
+#include "forceModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "postProcess.H"
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createControl.H"
+    #include "initContinuityErrs.H"
+    #include "createFields.H"
+    #include "createFvOptions.H"
+    #include "correctPhi.H"
+    #include "CourantNo.H"
+
+    turbulence->validate();
+
+    // create cfdemCloud
+    cfdemCloud particleCloud(mesh);
+
+    #include "additionalChecks.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.loop())
+    {
+        #include "CourantNo.H"
+        #include "alphaCourantNo.H"
+
+        particleCloud.clockM().start(1,"Global");
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        particleCloud.clockM().start(2,"Coupling");
+        bool hasEvolved = particleCloud.evolve(voidfraction,Us,U);
+
+        if(hasEvolved)
+        {
+            particleCloud.smoothingM().smoothen(particleCloud.forceM(0).impParticleForces());
+        }
+
+        Info << "update Ksl.internalField()" << endl;
+        Ksl = particleCloud.momCoupleM(0).impMomSource();
+        Ksl.correctBoundaryConditions();
+
+        //Force Checks
+        vector fTotal(0,0,0);
+        vector fImpTotal = sum(mesh.V()*Ksl.internalField()*(Us.internalField()-U.internalField())).value();
+        reduce(fImpTotal, sumOp<vector>());
+        Info << "TotalForceExp: " << fTotal << endl;
+        Info << "TotalForceImp: " << fImpTotal << endl;
+
+        #include "solverDebugInfo.H"
+        particleCloud.clockM().stop("Coupling");
+
+        particleCloud.clockM().start(26,"Flow");
+
+        if(particleCloud.solveFlow())
+        {
+            mixture.solve();
+            rho = mixture.rho();
+            rhoEps = rho * voidfraction;
+
+             // --- Pressure-velocity PIMPLE corrector loop
+            while (pimple.loop())
+            {
+                #include "UEqn.H"
+
+                // --- Pressure corrector loop
+                while (pimple.correct())
+                {
+                    #include "pEqn.H"
+                }
+
+                if (pimple.turbCorr())
+                {
+                    turbulence->correct();
+                }
+            }
+        }
+        else
+        {
+            Info << "skipping flow solution." << endl;
+        }
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+
+        particleCloud.clockM().stop("Flow");
+        particleCloud.clockM().stop("Global");
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/cfdemSolverMultiphase/correctPhi.H

@@ -0,0 +1,11 @@
+CorrectPhi
+(
+    U,
+    phi,
+    p_rgh,
+    dimensionedScalar("rAUf", dimTime/rho.dimensions(), 1),
+    geometricZeroField(),
+    pimple
+);
+
+#include "continuityErrs.H"