Merge pull request #154 from ParticulateFlow/release

Release 24.01

README.md

@@ -2,7 +2,7 @@
 
 CFDEM®coupling stands for Computational Fluid Dynamics (CFD) - Discrete Element Method (DEM) coupling. It combines the open source packages OpenFOAM® (CFD) and LIGGGHTS® (DEM) to simulate particle-laden flows. CFDEM®coupling is part of the [CFDEM®project](https://www.cfdem.com).
 
-[![CircleCI](https://circleci.com/gh/ParticulateFlow/CFDEMcoupling.svg?style=shield&circle-token=e4b6af30d3aa7aee109d206116f01600bf9ee9c6)](https://circleci.com/gh/ParticulateFlow/CFDEMcoupling)
+[![CircleCI](https://circleci.com/gh/ParticulateFlow/CFDEMcoupling.svg?style=shield&circle-token=7e8118524babddbefccf4e3608a7545d405acbb4)](https://circleci.com/gh/ParticulateFlow/CFDEMcoupling)
 [![License: GPL v3](https://img.shields.io/badge/License-GPL%20v3-blue.svg)](https://www.gnu.org/licenses/gpl-3.0.html)
 
 ## Disclaimer

applications/solvers/cfdemSolverIB/Make/options

@@ -1,7 +1,6 @@
 include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
 
 EXE_INC = \
-    -I$(CFDEM_OFVERSION_DIR) \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \

applications/solvers/cfdemSolverIBContinuousForcing/Make/options

@@ -1,7 +1,6 @@
 include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
 
 EXE_INC = \
-    -I$(CFDEM_OFVERSION_DIR) \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \

applications/solvers/cfdemSolverMultiphase/Make/options

@@ -5,7 +5,6 @@ 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 \

applications/solvers/cfdemSolverMultiphaseScalar/Make/options

@@ -5,7 +5,6 @@ include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
 
 EXE_INC = \
      $(PFLAGS) \
-    -I$(CFDEM_OFVERSION_DIR) \
     -ImultiphaseMixtureScalar/lnInclude \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/transportModels/incompressible/lnInclude \

applications/solvers/cfdemSolverPimple/Make/files

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

applications/solvers/cfdemSolverPimple/Make/options

@@ -0,0 +1,24 @@
+include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
+
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
+    -I$(LIB_SRC)/transportModels \
+    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
+    -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/lnInclude \
+    -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/cfdTools \
+    -Wno-deprecated-copy
+
+EXE_LIBS = \
+    -L$(CFDEM_LIB_DIR)\
+    -lturbulenceModels \
+    -lincompressibleTurbulenceModels \
+    -lincompressibleTransportModels \
+    -lfiniteVolume \
+    -lmeshTools \
+    -lfvOptions \
+    -l$(CFDEM_LIB_NAME) \
+     $(CFDEM_ADD_LIB_PATHS) \
+     $(CFDEM_ADD_LIBS)

applications/solvers/cfdemSolverPimple/UEqn.H

@@ -0,0 +1,43 @@
+particleCloud.otherForces(fOther);
+
+tmp<fvVectorMatrix> tUEqn
+(
+    fvm::ddt(voidfraction,U) - fvm::Sp(fvc::ddt(voidfraction),U)
+  + fvm::div(phi,U) - fvm::Sp(fvc::div(phi),U)
+    // in case of "periodic box" simulations the viscous term can be commented
+    // during testing the effect of this term on the results was negligible (about 1-2%)
+  + particleCloud.divVoidfractionTau(U, voidfraction)
+  - fOther/rho
+  ==
+    fvOptions(U)
+  - fvm::Sp(Ksl/rho,U)
+);
+
+fvVectorMatrix& UEqn = tUEqn.ref();
+
+UEqn.relax();
+
+fvOptions.constrain(UEqn);
+
+volScalarField rAU = 1.0/UEqn.A();
+
+surfaceScalarField rAUf("(1|A(U))", fvc::interpolate(rAU));
+
+surfaceScalarField voidfractionf = fvc::interpolate(voidfraction);
+
+surfaceScalarField phicForces
+(
+    fvc::interpolate(rAU*(Ksl*Us)/rho) & mesh.Sf()
+);
+
+if (pimple.momentumPredictor() && (modelType=="B" || modelType=="Bfull"))
+{
+    solve(UEqn == fvc::reconstruct(phicForces/rAUf - fvc::snGrad(p)*mesh.magSf()));
+    fvOptions.correct(U);
+}
+else if (pimple.momentumPredictor())
+{
+    solve(UEqn == fvc::reconstruct(phicForces/rAUf - fvc::snGrad(p)*voidfractionf*mesh.magSf()));
+    fvOptions.correct(U);
+}
+

applications/solvers/cfdemSolverPimple/cfdemSolverPimple.C

@@ -0,0 +1,147 @@
+/*---------------------------------------------------------------------------*\
+    Open Source CFD-DEM coupling
+
+    Copyright (C) 2023  Behrad Esgandari, JKU Linz, Austria
+-------------------------------------------------------------------------------
+License
+
+    This program 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 program 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 program.  If not, see <https://www.gnu.org/licenses/>.
+
+Application
+    cfdemSolverPimple
+
+Description
+    Transient solver for incompressible flow.
+    Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
+    The code is an evolution of the solver pimpleFoam in OpenFOAM(R) 6.0,
+    where additional functionality for CFD-DEM coupling is added.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "singlePhaseTransportModel.H"
+#include "turbulentTransportModel.H"
+#include "pimpleControl.H"
+#include "fvOptions.H"
+
+#include "cfdemCloud.H"
+#include "implicitCouple.H"
+#include "clockModel.H"
+#include "smoothingModel.H"
+#include "forceModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createControl.H"
+    #include "createFields.H"
+    #include "createFvOptions.H"
+    #include "initContinuityErrs.H"
+
+    // create cfdemCloud
+    #include "readGravitationalAcceleration.H"
+    cfdemCloud particleCloud(mesh);
+    #include "checkModelType.H"
+
+    // switch for periodic box simulations
+    Switch periodicBoxSwitch
+    (
+        pimple.dict().lookupOrDefault<Switch>("periodicBoxSwitch", false)
+    );
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+    Info<< "\nStarting time loop\n" << endl;
+    while (runTime.loop())
+    {
+        particleCloud.clockM().start(1,"Global");
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        #include "CourantNo.H"
+
+        // do particle stuff
+        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())
+        {
+            // Pressure-velocity PIMPLE corrector
+            while (pimple.loop())
+            {
+                // Momentum predictor
+                #include "UEqn.H"
+
+                // --- Inner PIMPLE loop
+
+                while (pimple.correct())
+                {
+                    #include "pEqn.H"
+                }
+            }
+
+            laminarTransport.correct();
+            turbulence->correct();
+        }
+        else
+        {
+            Info << "skipping flow solution." << endl;
+        }
+
+        if (periodicBoxSwitch)
+        {
+           #include "periodicBoxProperties.H"
+        }
+
+        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;
+}
+
+
+// ************************************************************************* //