Merge pull request #15 from ParticulateFlow/release

Release 17.02

.gitignore

@@ -0,0 +1,8 @@
+*.o
+*.d
+*.a
+*.dep
+log_*
+log.*
+*~
+**/linux64GccDPInt32Opt

README

@@ -5,7 +5,8 @@
     www.cfdem.com
                                 Christoph Goniva, christoph.goniva@cfdem.com
                                 Copyright 2009-2012 JKU Linz
-                                Copyright 2012-     DCS Computing GmbH, Linz
+                                Copyright 2012-2015 DCS Computing GmbH, Linz
+                                Copyright 2015-     JKU Linz
 -------------------------------------------------------------------------------
 License
     This file is part of CFDEMcoupling.
@@ -25,16 +26,17 @@ License
     Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
 
 Description
-    This code is designed to realize coupled CFD-DEM simulations using LIGGGHTS
-    and OpenFOAM. Note: this code is not part of OpenFOAM (see DISCLAIMER).
+    This code provides models and solvers to realize coupled CFD-DEM simulations
+    using LIGGGHTS and OpenFOAM.
+    Note: this code is not part of OpenFOAM (see DISCLAIMER).
 \*---------------------------------------------------------------------------*/
 
 
-CFDEM coupling provides an open source parallel coupled CFD-DEM framework 
-combining the strengths of LIGGGHTS DEM code and the Open Source 
-CFD package OpenFOAM(R)(*). The CFDEMcoupling toolbox allows to expand 
+CFDEM(R) coupling provides an open source parallel coupled CFD-DEM framework 
+combining the strengths of the LIGGGHTS(R) DEM code and the Open Source 
+CFD package OpenFOAM(R)(*). The CFDEM(R)coupling toolbox allows to expand 
 standard CFD solvers of OpenFOAM(R)(*) to include a coupling to the DEM 
-code LIGGGHTS. In this toolbox the particle representation within the 
+code LIGGGHTS(R). In this toolbox the particle representation within the 
 CFD solver is organized by "cloud" classes. Key functionalities are organised 
 in sub-models (e.g. force models, data exchange models, etc.) which can easily 
 be selected and combined by dictionary settings.
@@ -45,8 +47,6 @@ Features are:
 
 - its modular approach allows users to easily implement new models
 - its MPI parallelization enables to use it for large scale problems
-- the "forum"_lws on CFD-DEM gives the possibility to exchange with other 
-  users / developers
 - the use of GIT allows to easily update to the latest version
 - basic documentation is provided
 
@@ -54,7 +54,7 @@ The file structure:
 
 - "src" directory including the source files of the coupling toolbox and models
 - "applications" directory including the solver files for coupled CFD-DEM simulations
-- "doc" directory including the documentation of CFDEMcoupling
+- "doc" directory including the documentation of CFDEM(R)coupling
 - "tutorials" directory including basic tutorial cases showing the functionality
 
 
@@ -64,18 +64,17 @@ Details on installation are given on the "www.cfdem.com"
 The functionality of this CFD-DEM framwork is described via "tutorial cases" showing 
 how to use different solvers and models.
 
-CFDEMcoupling stands for Computational Fluid Dynamics (CFD) -
+CFDEM(R)coupling stands for Computational Fluid Dynamics (CFD) -
 Discrete Element Method (DEM) coupling.
 
-CFDEMcoupling is an open-source code, distributed freely under the terms of the 
+CFDEM(R)coupling is an open-source code, distributed freely under the terms of the 
 GNU Public License (GPL).
 
-Core development of CFDEMcoupling is done by 
+Core development of CFDEM(R)coupling is done by 
 Christoph Goniva and Christoph Kloss, both at DCS Computing GmbH, 2012
 
 
-\*---------------------------------------------------------------------------*/
-(*) "OpenFOAM(R)"_of is a registered trade mark of the ESI Group.
-This offering is not affiliated, approved or endorsed by ESI Group, 
-the producer of the OpenFOAM® software and owner of the OpenFOAM® trade mark.
+/*---------------------------------------------------------------------------*\
+(*) "OpenFOAM(R)" is a registered trade mark of OpenCFD Limited, a wholly owned subsidiary of the ESI Group.
+This offering is not approved or endorsed by OpenCFD Limited, the producer of the OpenFOAM software and owner of the OPENFOAM® and OpenCFD® trade marks.
 \*---------------------------------------------------------------------------*/

applications/.gitignore

@@ -0,0 +1,7 @@
+*.o
+*.d
+*.a
+*.dep
+log_*
+log.*
+*~

applications/solvers/cfdemSolverIB/Make/options

@@ -1,10 +1,13 @@
+include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
+
 EXE_INC = \
-    -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
+    -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$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/sampling/lnInclude \
     -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
@@ -13,12 +16,16 @@ EXE_INC = \
     -I$(LIB_SRC)/fvOptions/lnInclude
 
 EXE_LIBS = \
-    -L$(FOAM_USER_LIBBIN)\
-    -lincompressibleRASModels \
-    -lincompressibleLESModels \
+    -L$(CFDEM_LIB_DIR)\
+    -lturbulenceModels \
+    -lincompressibleTurbulenceModels \
     -lincompressibleTransportModels \
     -lfiniteVolume \
+    -lmeshTools \
     -ldynamicFvMesh \
     -ldynamicMesh \
     -lfvOptions \
-    -l$(CFDEM_LIB_NAME)
+    -l$(CFDEM_LIB_NAME) \
+     $(CFDEM_ADD_LIB_PATHS) \
+     $(CFDEM_ADD_LIBS)
+

applications/solvers/cfdemSolverIB/cfdemSolverIB.C

@@ -39,7 +39,8 @@ Contributions
 
 #include "fvCFD.H"
 #include "singlePhaseTransportModel.H"
-#include "turbulenceModel.H"
+#include "turbulentTransportModel.H"
+#include "pisoControl.H"
 
 #include "cfdemCloudIB.H"
 #include "implicitCouple.H"
@@ -52,11 +53,6 @@ Contributions
 
 #include "cellSet.H"
 
-#if defined(version22)
-    #include "meshToMeshNew.H"
-    #include "fvIOoptionList.H"
-#endif
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 int main(int argc, char *argv[])
@@ -67,14 +63,14 @@ int main(int argc, char *argv[])
 
     #include "createDynamicFvMesh.H"
 
+    #include "createControl.H"
+
+    #include "createTimeControls.H"
+
     #include "createFields.H"
 
     #include "initContinuityErrs.H"
 
-    #if defined(version22)
-        #include "createFvOptions.H"
-    #endif
-
     // create cfdemCloud
     #include "readGravitationalAcceleration.H"
     cfdemCloudIB particleCloud(mesh);
@@ -91,8 +87,9 @@ int main(int argc, char *argv[])
         interFace = mag(mesh.lookupObject<volScalarField>("voidfractionNext"));
         mesh.update(); //dyM
 
-        #include "readPISOControls.H"
+        #include "readTimeControls.H"
         #include "CourantNo.H"
+        #include "setDeltaT.H"
 
         // do particle stuff
         Info << "- evolve()" << endl;
@@ -107,43 +104,29 @@ int main(int argc, char *argv[])
                 fvm::ddt(voidfraction,U)
               + fvm::div(phi, U)
               + turbulence->divDevReff(U)
-                #if defined(version22)
-                ==
-                fvOptions(U)
-                #endif
             );
 
             UEqn.relax();
 
-            #if defined(version22)
-            fvOptions.constrain(UEqn);
-            #endif
-
-            if (momentumPredictor)
+            if (piso.momentumPredictor())
             {
                 solve(UEqn == -fvc::grad(p));
             }
 
             // --- PISO loop
-            for (int corr=0; corr<nCorr; corr++)
+            while (piso.correct())
             {
                 volScalarField rUA = 1.0/UEqn.A();
+                surfaceScalarField rUAf(fvc::interpolate(rUA));
 
                 U = rUA*UEqn.H();
-                #ifdef version23
-                phi = (fvc::interpolate(U) & mesh.Sf()); // there is a new version in 23x
-                #else
+
                 phi = (fvc::interpolate(U) & mesh.Sf())
-                    + fvc::ddtPhiCorr(rUA, U, phi);
-                #endif
+                    + rUAf*fvc::ddtCorr(U, phi);
+
                 adjustPhi(phi, U, p);
 
-                #if defined(version22)
-                fvOptions.relativeFlux(phi);
-                #endif
-
-                // Non-orthogonal pressure corrector loop
-                for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+                while (piso.correctNonOrthogonal())
                 {
                     // Pressure corrector
 
@@ -154,20 +137,9 @@ int main(int argc, char *argv[])
 
                     pEqn.setReference(pRefCell, pRefValue);
 
-                    if
-                    (
-                        corr == nCorr-1
-                     && nonOrth == nNonOrthCorr
-                    )
-                    {
-                        pEqn.solve(mesh.solver("pFinal"));
-                    }
-                    else
-                    {
-                        pEqn.solve();
-                    }
+                    pEqn.solve(mesh.solver(p.select(piso.finalInnerIter())));
 
-                    if (nonOrth == nNonOrthCorr)
+                    if (piso.finalNonOrthogonalIter())
                     {
                         phi -= pEqn.flux();
                     }
@@ -186,10 +158,6 @@ int main(int argc, char *argv[])
         volScalarField voidfractionNext=mesh.lookupObject<volScalarField>("voidfractionNext");
         particleCloud.calcVelocityCorrection(p,U,phiIB,voidfractionNext);
 
-        #if defined(version22)
-        fvOptions.correct(U);
-        #endif
-
         runTime.write();
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"

applications/solvers/cfdemSolverPiso/Make/options

@@ -1,15 +1,22 @@
+include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
+
 EXE_INC = \
-    -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
+    -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$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/lnInclude \
     -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/cfdTools \
 
 EXE_LIBS = \
-    -L$(FOAM_USER_LIBBIN)\
-    -lincompressibleRASModels \
-    -lincompressibleLESModels \
+    -L$(CFDEM_LIB_DIR)\
+    -lturbulenceModels \
+    -lincompressibleTurbulenceModels \
     -lincompressibleTransportModels \
     -lfiniteVolume \
-    -l$(CFDEM_LIB_NAME)
+    -lmeshTools \
+    -l$(CFDEM_LIB_NAME) \
+     $(CFDEM_ADD_LIB_PATHS) \
+     $(CFDEM_ADD_LIBS)

applications/solvers/cfdemSolverPiso/UEqn.H

@@ -0,0 +1,23 @@
+fvVectorMatrix UEqn
+(
+    fvm::ddt(voidfraction,U) - fvm::Sp(fvc::ddt(voidfraction),U)
+  + fvm::div(phi,U) - fvm::Sp(fvc::div(phi),U)
+  + particleCloud.divVoidfractionTau(U, voidfraction)
+  ==
+  - fvm::Sp(Ksl/rho,U)
+);
+
+UEqn.relax();
+
+fvOptions.constrain(UEqn);
+
+if (piso.momentumPredictor() && (modelType=="B" || modelType=="Bfull"))
+{
+    solve(UEqn == - fvc::grad(p) + Ksl/rho*Us);
+    fvOptions.correct(U);		    
+}
+else if (piso.momentumPredictor())
+{
+    solve(UEqn == - voidfraction*fvc::grad(p) + Ksl/rho*Us);
+    fvOptions.correct(U);
+}

applications/solvers/cfdemSolverPiso/cfdemSolverPiso.C

@@ -36,12 +36,15 @@ Description
 
 #include "fvCFD.H"
 #include "singlePhaseTransportModel.H"
-#include "turbulenceModel.H"
+#include "turbulentTransportModel.H"
+#include "pisoControl.H"
+#include "fvOptions.H"
 
 #include "cfdemCloud.H"
 #include "implicitCouple.H"
 #include "clockModel.H"
 #include "smoothingModel.H"
+#include "forceModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -50,7 +53,9 @@ 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
@@ -59,26 +64,35 @@ int main(int argc, char *argv[])
     #include "checkModelType.H"
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
+    Info<< "\nStarting time loop\n" << endl;
     while (runTime.loop())
     {
-        Info<< "\nStarting time loop\n" << endl;
-            particleCloud.clockM().start(1,"Global");
+        particleCloud.clockM().start(1,"Global");
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        #include "readPISOControls.H"
         #include "CourantNo.H"
 
         // do particle stuff
         particleCloud.clockM().start(2,"Coupling");
-        particleCloud.evolve(voidfraction,Us,U);
+        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();
-        particleCloud.smoothingM().smoothen(Ksl);
         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");
 
@@ -89,92 +103,19 @@ int main(int argc, char *argv[])
             // Pressure-velocity PISO corrector
             {
                 // Momentum predictor
-                fvVectorMatrix UEqn
-                (
-                    fvm::ddt(voidfraction,U) + fvm::Sp(fvc::ddt(voidfraction),U)
-                  + fvm::div(phi,U) + fvm::Sp(fvc::div(phi),U)
-//                + turbulence->divDevReff(U)
-                  + particleCloud.divVoidfractionTau(U, voidfraction)
-                  ==
-                  - fvm::Sp(Ksl/rho,U)
-                );
-
-                if (modelType=="B")
-                 UEqn == - fvc::grad(p) + Ksl/rho*Us;
-                else
-                    UEqn == - voidfraction*fvc::grad(p) + Ksl/rho*Us;
-
-                UEqn.relax();
-
-                if (momentumPredictor)
-                    solve(UEqn);
+                 #include "UEqn.H"
 
                 // --- PISO loop
 
-                //for (int corr=0; corr<nCorr; corr++)
-                int nCorrSoph = nCorr + 5 * pow((1-particleCloud.dataExchangeM().timeStepFraction()),1);
-
-                for (int corr=0; corr<nCorrSoph; corr++)
+                while (piso.correct())
                 {
-                    volScalarField rUA = 1.0/UEqn.A();
-
-                    surfaceScalarField rUAf("(1|A(U))", fvc::interpolate(rUA));
-                    volScalarField rUAvoidfraction("(voidfraction2|A(U))",rUA*voidfraction);
-
-                    U = rUA*UEqn.H();
-
-                    phi = (fvc::interpolate(U*voidfraction) & mesh.Sf() );
-                     //+ fvc::ddtPhiCorr(rUAvoidfraction, U, phi);
-                    surfaceScalarField phiS(fvc::interpolate(Us*voidfraction) & mesh.Sf());
-                    surfaceScalarField phiGes = phi + rUAf*(fvc::interpolate(Ksl/rho) * phiS);
-
-                    if (modelType=="A")
-                        rUAvoidfraction = volScalarField("(voidfraction2|A(U))",rUA*voidfraction*voidfraction);
-
-                    // Non-orthogonal pressure corrector loop
-                    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
-                    {
-                        // Pressure corrector
-                        fvScalarMatrix pEqn
-                        (
-                            fvm::laplacian(rUAvoidfraction, p) == fvc::div(phiGes) + particleCloud.ddtVoidfraction()
-                        );
-                        pEqn.setReference(pRefCell, pRefValue);
-
-                        if
-                        (
-                            corr == nCorr-1
-                         && nonOrth == nNonOrthCorr
-                        )
-                        {
-                            pEqn.solve(mesh.solver("pFinal"));
-                        }
-                        else
-                        {
-                            pEqn.solve();
-                        }
-
-                        if (nonOrth == nNonOrthCorr)
-                        {
-                            phiGes -= pEqn.flux();
-                        }
-
-                    } // end non-orthogonal corrector loop
-
-                    #include "continuityErrorPhiPU.H"
-
-                    if (modelType=="B")
-                        U -= rUA*fvc::grad(p) - Ksl/rho*Us*rUA;
-                    else
-                        U -= voidfraction*rUA*fvc::grad(p) - Ksl/rho*Us*rUA;
-
-                    U.correctBoundaryConditions();
-
-                } // end piso loop
+                    #include "pEqn.H"
+                }
             }
 
+            laminarTransport.correct();
             turbulence->correct();
-        }// end solveFlow
+        }
         else
         {
             Info << "skipping flow solution." << endl;
@@ -191,7 +132,7 @@ int main(int argc, char *argv[])
     }
 
     Info<< "End\n" << endl;
-    
+
     return 0;
 }
 

applications/solvers/cfdemSolverPiso/createFields.H

@@ -122,3 +122,5 @@ surfaceScalarField phi
     (
         incompressible::turbulenceModel::New(U, phi, laminarTransport)
     );
+
+#include "createMRF.H"

applications/solvers/cfdemSolverPiso/pEqn.H

@@ -0,0 +1,58 @@
+volScalarField rAU = 1.0/UEqn.A();
+
+surfaceScalarField rAUf("(1|A(U))", fvc::interpolate(rAU));
+
+volScalarField rAUvoidfraction("(voidfraction2|A(U))",rAU*voidfraction);
+
+if (modelType=="A")
+    rAUvoidfraction *= voidfraction;
+
+volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
+
+surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(voidfraction)*fvc::flux(HbyA) );
+
+volVectorField Uvoidfraction("(Uvoidfraction)", U * voidfraction);
+
+surfaceScalarField phiS(fvc::interpolate(Us*voidfraction) & mesh.Sf());
+
+phi = phiHbyA + rAUf*(fvc::interpolate(Ksl/rho) * phiS);
+
+// rotating frames of references not tested yet
+MRF.makeRelative(phi);
+
+// adjustment of phi (only in cases w.o. p boundary conditions) not tested yet
+adjustPhi(phi, U, p);
+
+// Update the pressure BCs to ensure flux consistency
+constrainPressure(p, Uvoidfraction, phiHbyA, rAUvoidfraction, MRF);
+
+
+// Non-orthogonal pressure corrector loop
+while (piso.correctNonOrthogonal())
+{
+    // Pressure corrector
+  
+    fvScalarMatrix pEqn
+    (
+        fvm::laplacian(rAUvoidfraction, p) == fvc::div(phi) + particleCloud.ddtVoidfraction()
+    );
+    
+    pEqn.setReference(pRefCell, pRefValue);
+
+    pEqn.solve(mesh.solver(p.select(piso.finalInnerIter())));
+
+    if (piso.finalNonOrthogonalIter())
+    {
+        phi -= pEqn.flux();
+    }
+}
+
+#include "continuityErrorPhiPU.H"
+
+if (modelType=="B" || modelType=="Bfull")
+    U = HbyA - rAU*fvc::grad(p) + Ksl/rho*Us*rAU;
+else
+    U = HbyA - voidfraction*rAU*fvc::grad(p) + Ksl/rho*Us*rAU;
+
+U.correctBoundaryConditions();
+fvOptions.correct(U);

applications/solvers/cfdemSolverPisoMS/Make/options

@@ -1,16 +1,23 @@
+include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
+
 EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/lnInclude \
-    -I ../cfdemSolverPiso \
-    -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
+    -I../cfdemSolverPiso \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/cfdTools \
 
 EXE_LIBS = \
-    -L$(FOAM_USER_LIBBIN)\
-    -lincompressibleRASModels \
-    -lincompressibleLESModels \
+    -L$(CFDEM_LIB_DIR)\
+    -lturbulenceModels \
+    -lincompressibleTurbulenceModels \
     -lincompressibleTransportModels \
     -lfiniteVolume \
-    -l$(CFDEM_LIB_NAME)
+    -lmeshTools \
+    -l$(CFDEM_LIB_NAME) \
+     $(CFDEM_ADD_LIB_PATHS) \
+     $(CFDEM_ADD_LIBS)

applications/solvers/cfdemSolverPisoMS/cfdemSolverPisoMS.C

@@ -25,23 +25,26 @@ License
     along with CFDEMcoupling.  If not, see <http://www.gnu.org/licenses/>.
 
 Application
-    cfdemSolverPisoMS
+    cfdemSolverPiso
 
 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 pisoFoam in OpenFOAM(R) 1.6, 
+    The code is an evolution of the solver pisoFoam in OpenFOAM(R) 1.6,
     where additional functionality for CFD-DEM coupling is added.
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
 #include "singlePhaseTransportModel.H"
-#include "turbulenceModel.H"
+#include "turbulentTransportModel.H"
+#include "pisoControl.H"
+#include "fvOptions.H"
 
 #include "cfdemCloudMS.H"
 #include "implicitCouple.H"
 #include "clockModel.H"
 #include "smoothingModel.H"
+#include "forceModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -50,7 +53,9 @@ 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
@@ -59,118 +64,62 @@ int main(int argc, char *argv[])
     #include "checkModelType.H"
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
+    Info<< "\nStarting time loop\n" << endl;
     while (runTime.loop())
     {
-        Info<< "\nStarting time loop\n" << endl;
-            particleCloud.clockM().start(1,"Global");
+        particleCloud.clockM().start(1,"Global");
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        #include "readPISOControls.H"
         #include "CourantNo.H"
 
         // do particle stuff
         particleCloud.clockM().start(2,"Coupling");
-        particleCloud.evolve(voidfraction,Us,U);
-      
+        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();
-        particleCloud.smoothingM().smoothen(Ksl);
         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");
-        // Pressure-velocity PISO corrector
+
+        if(particleCloud.solveFlow())
         {
-            // Momentum predictor
-            fvVectorMatrix UEqn
-            (
-                fvm::ddt(voidfraction,U) //particleCloud.ddtVoidfractionU(U,voidfraction) //
-              + fvm::div(phi, U)
-//              + turbulence->divDevReff(U)
-              + particleCloud.divVoidfractionTau(U, voidfraction)
-             ==
-              - fvm::Sp(Ksl/rho,U)
-            );
-
-            if (modelType=="B")
-                UEqn == - fvc::grad(p) + Ksl/rho*Us;
-            else
-                UEqn == - voidfraction*fvc::grad(p) + Ksl/rho*Us;
-
-            UEqn.relax();
-
-            if (momentumPredictor)
-                solve(UEqn);
-
-            // --- PISO loop
-
-            //for (int corr=0; corr<nCorr; corr++)
-            int nCorrSoph = nCorr + 5 * pow((1-particleCloud.dataExchangeM().timeStepFraction()),1);
-
-            for (int corr=0; corr<nCorrSoph; corr++)
+            // Pressure-velocity PISO corrector
             {
-                volScalarField rUA = 1.0/UEqn.A();
+                // Momentum predictor
+                 #include "UEqn.H"
 
-                surfaceScalarField rUAf("(1|A(U))", fvc::interpolate(rUA));
-                volScalarField rUAvoidfraction("(voidfraction2|A(U))",rUA*voidfraction);
+                // --- PISO loop
 
-                U = rUA*UEqn.H();
-
-                phi = (fvc::interpolate(U*voidfraction) & mesh.Sf() );
-                     //+ fvc::ddtPhiCorr(rUAvoidfraction, U, phi);
-                surfaceScalarField phiS(fvc::interpolate(Us*voidfraction) & mesh.Sf());
-                surfaceScalarField phiGes = phi + rUAf*(fvc::interpolate(Ksl/rho) * phiS);
-
-                if (modelType=="A")
-                    rUAvoidfraction = volScalarField("(voidfraction2|A(U))",rUA*voidfraction*voidfraction);
-
-                // Non-orthogonal pressure corrector loop
-                for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+                while (piso.correct())
                 {
-                    // Pressure corrector
-                    fvScalarMatrix pEqn
-                    (
-                        fvm::laplacian(rUAvoidfraction, p) == fvc::div(phiGes) + particleCloud.ddtVoidfraction()
-                    );
-                    pEqn.setReference(pRefCell, pRefValue);
+                    #include "pEqn.H"
+                }
+            }
 
-                    if
-                    (
-                        corr == nCorr-1
-                     && nonOrth == nNonOrthCorr
-                    )
-                    {
-                        pEqn.solve(mesh.solver("pFinal"));
-                    }
-                    else
-                    {
-                        pEqn.solve();
-                    }
-
-                    if (nonOrth == nNonOrthCorr)
-                    {
-                        phiGes -= pEqn.flux();
-                    }
-
-                } // end non-orthogonal corrector loop
-
-                #include "continuityErrorPhiPU.H"
-
-                if (modelType=="B")
-                    U -= rUA*fvc::grad(p) - Ksl/rho*Us*rUA;
-                else
-                    U -= voidfraction*rUA*fvc::grad(p) - Ksl/rho*Us*rUA;
-
-                U.correctBoundaryConditions();
-
-            } // end piso loop
+            laminarTransport.correct();
+            turbulence->correct();
+        }
+        else
+        {
+            Info << "skipping flow solution." << endl;
         }
-
-        turbulence->correct();
 
         runTime.write();
 
@@ -183,7 +132,7 @@ int main(int argc, char *argv[])
     }
 
     Info<< "End\n" << endl;
-    
+
     return 0;
 }
 

applications/solvers/cfdemSolverPisoScalar/Make/options

@@ -1,15 +1,23 @@
+include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
+
 EXE_INC = \
-    -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
+    -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$(LIB_SRC)/finiteVolume/lnInclude \
+    -I../cfdemSolverPiso \
     -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/lnInclude \
     -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/cfdTools \
 
 EXE_LIBS = \
-    -L$(FOAM_USER_LIBBIN)\
-    -lincompressibleRASModels \
-    -lincompressibleLESModels \
+    -L$(CFDEM_LIB_DIR)\
+    -lturbulenceModels \
+    -lincompressibleTurbulenceModels \
     -lincompressibleTransportModels \
     -lfiniteVolume \
-    -l$(CFDEM_LIB_NAME)
+    -lmeshTools \
+    -l$(CFDEM_LIB_NAME) \
+     $(CFDEM_ADD_LIB_PATHS) \
+     $(CFDEM_ADD_LIBS)

applications/solvers/cfdemSolverPisoScalar/TEqn.H

@@ -0,0 +1,15 @@
+ // get scalar source from DEM        
+        particleCloud.forceM(1).manipulateScalarField(Tsource);
+        Tsource.correctBoundaryConditions();
+
+        // solve scalar transport equation
+        fvScalarMatrix TEqn
+        (
+           fvm::ddt(voidfraction,T) - fvm::Sp(fvc::ddt(voidfraction),T)
+         + fvm::div(phi, T) - fvm::Sp(fvc::div(phi),T)
+         - fvm::laplacian(DT*voidfraction, T)
+         ==
+           Tsource
+        );
+        TEqn.relax();
+        TEqn.solve();

applications/solvers/cfdemSolverPisoScalar/cfdemSolverPisoScalar.C

@@ -25,174 +25,112 @@ License
     along with CFDEMcoupling.  If not, see <http://www.gnu.org/licenses/>.
 
 Application
-    cfdemSolverPisoScalar
+    cfdemSolverPiso
 
 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 pisoFoam in OpenFOAM(R) 1.6, 
+    The code is an evolution of the solver pisoFoam in OpenFOAM(R) 1.6,
     where additional functionality for CFD-DEM coupling is added.
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
 #include "singlePhaseTransportModel.H"
-#include "turbulenceModel.H"
+#include "turbulentTransportModel.H"
+#include "pisoControl.H"
+#include "fvOptions.H"
 
 #include "cfdemCloud.H"
 #include "implicitCouple.H"
-#include "forceModel.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"
+
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
     Info<< "\nStarting time loop\n" << endl;
-
     while (runTime.loop())
     {
+        particleCloud.clockM().start(1,"Global");
+
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        #include "readPISOControls.H"
         #include "CourantNo.H"
 
         // do particle stuff
-        Info << "- evolve()" << endl;
-        particleCloud.evolve(voidfraction,Us,U);
+        particleCloud.clockM().start(2,"Coupling");
+        bool hasEvolved = particleCloud.evolve(voidfraction,Us,U);
 
-        Ksl.internalField() = particleCloud.momCoupleM(0).impMomSource();
-        particleCloud.smoothingM().smoothen(Ksl);
+        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");
 
-        // get scalar source from DEM        
-        particleCloud.forceM(1).manipulateScalarField(Tsource);
-        Tsource.correctBoundaryConditions();
+        particleCloud.clockM().start(26,"Flow");
+	
+        #include "TEqn.H"
 
-        // solve scalar transport equation
-        phi = fvc::interpolate(U*voidfraction) & mesh.Sf();
-
-        solve
-        (
-            fvm::ddt(voidfraction,T)
-          + fvm::div(phi, T)
-          - fvm::laplacian(DT*voidfraction, T)
-          ==
-            Tsource
-        );
-
-        // Pressure-velocity PISO corrector
+        if(particleCloud.solveFlow())
         {
-            // Momentum predictor
-            fvVectorMatrix UEqn
-            (
-                fvm::ddt(voidfraction,U)
-              + fvm::div(phi, U)
-              + turbulence->divDevReff(U)
-             == 
-              - fvm::Sp(Ksl/rho,U)
-            );
-
-            UEqn.relax();
-
-            if (momentumPredictor)
+            // Pressure-velocity PISO corrector
             {
-                //solve UEqn
-                if (modelType=="B")
-                    solve(UEqn == - fvc::grad(p) + Ksl/rho*Us);
-                else
-                    solve(UEqn == - voidfraction*fvc::grad(p) + Ksl/rho*Us);
+                // Momentum predictor
+                 #include "UEqn.H"
+
+                // --- PISO loop
+
+                while (piso.correct())
+                {
+                    #include "pEqn.H"
+                }
             }
 
-            // --- PISO loop
-
-            //for (int corr=0; corr<nCorr; corr++)
-            int nCorrSoph = nCorr + 5 * pow((1-particleCloud.dataExchangeM().timeStepFraction()),1);
-            for (int corr=0; corr<nCorrSoph; corr++)
-            {
-                volScalarField rUA = 1.0/UEqn.A();
-                surfaceScalarField rUAf("(1|A(U))", fvc::interpolate(rUA));
-                volScalarField rUAvoidfraction("(voidfraction2|A(U))",rUA*voidfraction);
-
-                U = rUA*UEqn.H();
-
-                #ifdef version23
-                phi = ( fvc::interpolate(U*voidfraction) & mesh.Sf() );
-                #else
-                phi = ( fvc::interpolate(U*voidfraction) & mesh.Sf() )
-                      + fvc::ddtPhiCorr(rUAvoidfraction, U, phi);
-                #endif
-                surfaceScalarField phiS(fvc::interpolate(Us*voidfraction) & mesh.Sf());
-                surfaceScalarField phiGes = phi + rUAf*(fvc::interpolate(Ksl/rho) * phiS);
-                if (modelType=="A")
-                    rUAvoidfraction = volScalarField("(voidfraction2|A(U))",rUA*voidfraction*voidfraction);
-
-                // Non-orthogonal pressure corrector loop
-                for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
-                {
-                    // Pressure corrector
-                    fvScalarMatrix pEqn
-                    (
-                        fvm::laplacian(rUAvoidfraction, p) == fvc::div(phiGes) + particleCloud.ddtVoidfraction()
-                    );
-                    pEqn.setReference(pRefCell, pRefValue);
-
-                    if
-                    (
-                        corr == nCorr-1
-                     && nonOrth == nNonOrthCorr
-                    )
-                    {
-                        pEqn.solve(mesh.solver("pFinal"));
-                    }
-                    else
-                    {
-                        pEqn.solve();
-                    }
-
-                    if (nonOrth == nNonOrthCorr)
-                    {
-                        phiGes -= pEqn.flux();
-                    }
-
-                } // end non-orthogonal corrector loop
-
-                #include "continuityErrorPhiPU.H"
-
-                if (modelType=="B")
-                    U -= rUA*fvc::grad(p) - Ksl/rho*Us*rUA;
-                else
-                    U -= voidfraction*rUA*fvc::grad(p) - Ksl/rho*Us*rUA;
-
-                U.correctBoundaryConditions();
-
-            } // end piso loop
+            laminarTransport.correct();
+            turbulence->correct();
+        }
+        else
+        {
+            Info << "skipping flow solution." << endl;
         }
-
-        turbulence->correct();
 
         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;

applications/solvers/cfdemSolverPisoScalar/createFields.H

@@ -1,36 +1,36 @@
-    Info<< "Reading field p\n" << endl;
-    volScalarField p
-    (
-        IOobject
-        (
-            "p",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-    Info<< "Reading physical velocity field U" << endl;
-    Info<< "Note: only if voidfraction at boundary is 1, U is superficial velocity!!!\n" << endl;
-    volVectorField U
-    (
-        IOobject
-        (
-            "U",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-//========================
-// drag law modelling
-//========================
-
+    Info<< "Reading field p\n" << endl;
+    volScalarField p
+    (
+        IOobject
+        (
+            "p",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< "Reading physical velocity field U" << endl;
+    Info<< "Note: only if voidfraction at boundary is 1, U is superficial velocity!!!\n" << endl;
+    volVectorField U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+//========================
+// drag law modelling
+//========================
+
     Info<< "\nReading momentum exchange field Ksl\n" << endl;
     volScalarField Ksl
     (
@@ -44,8 +44,8 @@
         ),
         mesh
         //dimensionedScalar("0", dimensionSet(0, 0, -1, 0, 0), 1.0)
-    );
-
+    );
+
     Info<< "\nReading voidfraction field voidfraction = (Vgas/Vparticle)\n" << endl;
     volScalarField voidfraction
     (
@@ -58,8 +58,8 @@
             IOobject::AUTO_WRITE
         ),
         mesh
-    );
-
+    );
+
     Info<< "\nCreating density field rho\n" << endl;
     volScalarField rho
     (
@@ -71,27 +71,27 @@
             IOobject::READ_IF_PRESENT,
             IOobject::AUTO_WRITE
         ),
-        mesh,
+        mesh,
         dimensionedScalar("0", dimensionSet(1, -3, 0, 0, 0), 1.0)
-    );
-
-    Info<< "Reading particle velocity field Us\n" << endl;
-    volVectorField Us
-    (
-        IOobject
-        (
-            "Us",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-//========================
-// scalar field modelling
-//========================
+    );
+
+    Info<< "Reading particle velocity field Us\n" << endl;
+    volVectorField Us
+    (
+        IOobject
+        (
+            "Us",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+//========================
+// scalar field modelling
+//========================
     Info<< "\nCreating dummy density field rho = 1\n" << endl;
     volScalarField T
     (
@@ -103,10 +103,10 @@
             IOobject::MUST_READ,
             IOobject::AUTO_WRITE
         ),
-        mesh//,
+        mesh//,
         //dimensionedScalar("0", dimensionSet(0, 0, -1, 1, 0), 273.15)
-    );
-
+    );
+
     Info<< "\nCreating fluid-particle heat flux field\n" << endl;
     volScalarField Tsource
     (
@@ -118,57 +118,59 @@
             IOobject::MUST_READ,
             IOobject::AUTO_WRITE
         ),
-        mesh//,
+        mesh//,
         //dimensionedScalar("0", dimensionSet(0, 0, -1, 1, 0), 0.0)
-    );
-
-    IOdictionary transportProperties
-    (
-        IOobject
-        (
-            "transportProperties",
-            runTime.constant(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE
-        )
-    );
-
-    dimensionedScalar DT
-    (
-        transportProperties.lookup("DT")
-    );
-
-//========================
-
-//#   include "createPhi.H"
-#ifndef createPhi_H
-#define createPhi_H
-Info<< "Reading/calculating face flux field phi\n" << endl;
-surfaceScalarField phi
- (
-     IOobject
-     (
-        "phi",
-        runTime.timeName(),
-         mesh,
-        IOobject::READ_IF_PRESENT,
-        IOobject::AUTO_WRITE
-     ),
-     linearInterpolate(U*voidfraction) & mesh.Sf()
- );
-#endif
-
-
-
-    label pRefCell = 0;
-    scalar pRefValue = 0.0;
-    setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);
-
-
-    singlePhaseTransportModel laminarTransport(U, phi);
-
-    autoPtr<incompressible::turbulenceModel> turbulence
-    (
-        incompressible::turbulenceModel::New(U, phi, laminarTransport)
-    );
+    );
+
+    IOdictionary transportProperties
+    (
+        IOobject
+        (
+            "transportProperties",
+            runTime.constant(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        )
+    );
+
+    dimensionedScalar DT
+    (
+        transportProperties.lookup("DT")
+    );
+
+//========================
+
+//#   include "createPhi.H"
+#ifndef createPhi_H
+#define createPhi_H
+Info<< "Reading/calculating face flux field phi\n" << endl;
+surfaceScalarField phi
+ (
+     IOobject
+     (
+        "phi",
+        runTime.timeName(),
+         mesh,
+        IOobject::READ_IF_PRESENT,
+        IOobject::AUTO_WRITE
+     ),
+     linearInterpolate(U*voidfraction) & mesh.Sf()
+ );
+#endif
+
+
+
+    label pRefCell = 0;
+    scalar pRefValue = 0.0;
+    setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);
+
+
+    singlePhaseTransportModel laminarTransport(U, phi);
+
+    autoPtr<incompressible::turbulenceModel> turbulence
+    (
+        incompressible::turbulenceModel::New(U, phi, laminarTransport)
+    );
+
+#include "createMRF.H"

applications/solvers/cfdemSolverRhoPimple/EEqn.H

@@ -0,0 +1,59 @@
+// contributions to internal energy equation can be found in
+// Crowe et al.: "Multiphase flows with droplets and particles", CRC Press 1998
+{
+    // dim he = J / kg
+    volScalarField& he = thermo.he();
+    particleCloud.energyContributions(Qsource);
+    particleCloud.energyCoefficients(QCoeff);
+
+    //thDiff=particleCloud.thermCondM().thermDiff();   
+    thCond=particleCloud.thermCondM().thermCond(); 
+
+    addSource = fvc::ddt(rhoeps, K) + fvc::div(phi, K)
+      + (
+            he.name() == "e"
+          ? fvc::div
+            (
+                fvc::absolute(phi/fvc::interpolate(rho), voidfraction*U),
+                p,
+                "div(phiv,p)"
+            )
+          : -dpdt
+        );
+
+    Cpv = he.name() == "e" ? thermo.Cv() : thermo.Cp();
+
+
+    fvScalarMatrix EEqn
+    (
+        fvm::ddt(rhoeps, he) + fvm::div(phi, he)
+      + addSource
+     // net heat transfer from particles to fluid
+      - Qsource
+      - fvm::Sp(QCoeff/Cpv, he)
+     // thermal conduction of the fluid with effective conductivity
+     // - fvm::laplacian(rhoeps*thDiff,he)
+      - fvm::laplacian(voidfraction*thCond/Cpv,he)
+      // + particle-fluid energy transfer due to work
+      // + fluid energy dissipation due to shearing
+     ==
+        fvOptions(rho, he)
+    );
+    
+
+    EEqn.relax();
+
+    fvOptions.constrain(EEqn);
+
+    EEqn.solve();
+
+    fvOptions.correct(he);
+
+    thermo.correct();
+
+    Info<< "T max/min : " << max(T).value() << " " << min(T).value() << endl;
+
+    particleCloud.clockM().start(31,"postFlow");
+    particleCloud.postFlow();
+    particleCloud.clockM().stop("postFlow");
+}

applications/solvers/cfdemSolverRhoPimple/Make/files

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

applications/solvers/cfdemSolverRhoPimple/Make/options

@@ -0,0 +1,31 @@
+include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
+
+PFLAGS+= -Dcompre
+
+EXE_INC = \
+     $(PFLAGS) \
+    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/cfdTools \
+    -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)\
+    -lcompressibleTransportModels \
+    -lfluidThermophysicalModels \
+    -lspecie \
+    -lturbulenceModels \
+    -lcompressibleTurbulenceModels \
+    -lfiniteVolume \
+    -lmeshTools \
+    -lsampling \
+    -lfvOptions \
+    -l$(CFDEM_LIB_COMP_NAME) \
+     $(CFDEM_ADD_LIB_PATHS) \
+     $(CFDEM_ADD_LIBS)

applications/solvers/cfdemSolverRhoPimple/UEqn.H

@@ -0,0 +1,33 @@
+// Solve the Momentum equation
+particleCloud.otherForces(fOther);
+
+tmp<fvVectorMatrix> tUEqn
+(
+    fvm::ddt(rhoeps, U)
+  + fvm::div(phi, U)
+  + particleCloud.divVoidfractionTau(U, voidfraction)
+  + fvm::Sp(Ksl,U)
+  - fOther
+ ==
+    fvOptions(rho, U)
+);
+fvVectorMatrix& UEqn = tUEqn.ref();
+
+UEqn.relax();
+
+fvOptions.constrain(UEqn);
+
+if (pimple.momentumPredictor() && (modelType=="B" || modelType=="Bfull"))
+{
+    solve(UEqn == -fvc::grad(p)+ Ksl*Us);
+
+    fvOptions.correct(U);
+    K = 0.5*magSqr(U);
+}
+else if (pimple.momentumPredictor())
+{
+    solve(UEqn == -voidfraction*fvc::grad(p)+ Ksl*Us);
+
+    fvOptions.correct(U);
+    K = 0.5*magSqr(U);
+}

applications/solvers/cfdemSolverRhoPimple/cfdemSolverRhoPimple.C

@@ -0,0 +1,159 @@
+/*---------------------------------------------------------------------------*\
+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) 2015- Thomas Lichtenegger, JKU Linz, Austria
+
+Application
+    cfdemSolverRhoPimple
+
+Description
+    Transient solver for compressible flow using the flexible PIMPLE (PISO-SIMPLE)
+    algorithm.
+    Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
+    The code is an evolution of the solver rhoPimpleFoam in OpenFOAM(R) 4.x,
+    where additional functionality for CFD-DEM coupling is added.
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "psiThermo.H"
+#include "turbulentFluidThermoModel.H"
+#include "bound.H"
+#include "pimpleControl.H"
+#include "fvOptions.H"
+#include "localEulerDdtScheme.H"
+#include "fvcSmooth.H"
+
+
+#include "cfdemCloudEnergy.H"
+#include "implicitCouple.H"
+#include "clockModel.H"
+#include "smoothingModel.H"
+#include "forceModel.H"
+#include "thermCondModel.H"
+#include "energyModel.H"
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "postProcess.H"
+
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createControl.H"
+    #include "createTimeControls.H"
+    #include "createRDeltaT.H"
+    #include "initContinuityErrs.H"
+    #include "createFields.H"
+    #include "createFieldRefs.H"
+    #include "createFvOptions.H"
+
+    // create cfdemCloud
+    #include "readGravitationalAcceleration.H"
+    cfdemCloudEnergy particleCloud(mesh);
+    #include "checkModelType.H"
+
+    turbulence->validate();
+  //        #include "compressibleCourantNo.H"
+  //  #include "setInitialDeltaT.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.run())
+    {
+        #include "readTimeControls.H"
+        #include "compressibleCourantNo.H"
+        #include "setDeltaT.H"
+
+        runTime++;
+
+        particleCloud.clockM().start(1,"Global");
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        // 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.primitiveFieldRef()*(Us.primitiveFieldRef()-U.primitiveFieldRef()));
+        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 (pimple.nCorrPIMPLE() <= 1)
+        {
+            #include "rhoEqn.H"
+        }
+
+        volScalarField rhoeps("rhoeps",rho*voidfraction);
+        // --- Pressure-velocity PIMPLE corrector loop
+        while (pimple.loop())
+        {
+            #include "UEqn.H"
+            #include "EEqn.H"
+
+            // --- Pressure corrector loop
+            while (pimple.correct())
+            {
+                // besides this pEqn, OF offers a "pimple consistent"-option
+                #include "pEqn.H"
+                rhoeps=rho*voidfraction;
+            }
+
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
+            }
+        }
+
+        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/cfdemSolverRhoPimple/createFieldRefs.H

@@ -0,0 +1,2 @@
+const volScalarField& T = thermo.T();
+const volScalarField& psi = thermo.psi();

applications/solvers/cfdemSolverRhoPimple/createFields.H

@@ -0,0 +1,253 @@
+Info<< "Reading thermophysical properties\n" << endl;
+
+    autoPtr<psiThermo> pThermo
+    (
+        psiThermo::New(mesh)
+    );
+    psiThermo& thermo = pThermo();
+    thermo.validate(args.executable(), "h", "e");
+    volScalarField& p = thermo.p();
+
+    Info<< "Reading field rho\n" << endl;
+    volScalarField rho
+    (
+        IOobject
+        (
+            "rho",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        thermo.rho()
+    );
+
+
+    Info<< "Reading field U\n" << endl;
+    volVectorField U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< "\nReading voidfraction field voidfraction = (Vgas/Vparticle)\n" << endl;
+    volScalarField voidfraction
+    (
+        IOobject
+        (
+            "voidfraction",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    volScalarField addSource
+    (
+        IOobject
+        (
+           "addSource",
+           runTime.timeName(),
+           mesh,
+           IOobject::MUST_READ,
+           IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< "\nCreating fluid-particle heat flux field\n" << endl;
+    volScalarField Qsource
+    (
+        IOobject
+        (
+            "Qsource",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("zero", dimensionSet(1,-1,-3,0,0,0,0), 0.0)
+    );
+
+    Info<< "\nCreating fluid-particle heat flux coefficient field\n" << endl;
+    volScalarField QCoeff
+    (
+        IOobject
+        (
+            "QCoeff",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("zero", dimensionSet(1,-1,-3,-1,0,0,0), 0.0)
+    );
+
+  /*  Info<< "\nCreating thermal diffusivity field\n" << endl;
+    volScalarField thDiff
+    (
+        IOobject
+        (
+            "thDiff",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("zero", dimensionSet(0,2,-1,0,0,0,0), 0.0)
+    );
+  */
+    Info<< "\nCreating thermal conductivity field\n" << endl;
+    volScalarField thCond
+    (
+        IOobject
+        (
+            "thCond",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("zero", dimensionSet(1,1,-3,-1,0,0,0), 0.0)
+    );
+
+    Info<< "\nCreating heat capacity field\n" << endl;
+    volScalarField Cpv
+    (
+        IOobject
+        (
+            "Cpv",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("zero", dimensionSet(0,2,-2,-1,0,0,0), 0.0)
+    );
+
+    Info<< "\nCreating body force field\n" << endl;
+    volVectorField fOther
+    (
+        IOobject
+        (
+            "fOther",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedVector("zero", dimensionSet(1,-2,-2,0,0,0,0), vector::zero)
+    );
+
+    Info<< "Reading/calculating face flux field phi\n" << endl;
+    surfaceScalarField phi
+    (
+        IOobject
+        (
+            "phi",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+         ),
+         linearInterpolate(rho*U*voidfraction) & mesh.Sf()
+    );
+
+    dimensionedScalar rhoMax
+    (
+        dimensionedScalar::lookupOrDefault
+        (
+            "rhoMax",
+            pimple.dict(),
+            dimDensity,
+            GREAT
+        )
+    );
+
+    dimensionedScalar rhoMin
+    (
+        dimensionedScalar::lookupOrDefault
+        (
+            "rhoMin",
+            pimple.dict(),
+            dimDensity,
+            0
+        )
+    );
+
+    Info<< "Creating turbulence model\n" << endl;
+    autoPtr<compressible::turbulenceModel> turbulence
+    (
+        compressible::turbulenceModel::New
+        (
+            rho,
+            U,
+            phi,
+            thermo
+        )
+    );
+
+    mesh.setFluxRequired(p.name());
+
+    Info<< "Creating field dpdt\n" << endl;
+    volScalarField dpdt
+    (
+        IOobject
+        (
+            "dpdt",
+            runTime.timeName(),
+            mesh
+        ),
+        mesh,
+        dimensionedScalar("dpdt", p.dimensions()/dimTime, 0)
+    );
+
+    Info<< "Creating field kinetic energy K\n" << endl;
+    volScalarField K("K", 0.5*magSqr(U));
+
+    Info<< "\nReading momentum exchange field Ksl\n" << endl;
+    volScalarField Ksl
+    (
+        IOobject
+        (
+            "Ksl",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+        //dimensionedScalar("0", dimensionSet(1, -3, -1, 0, 0), 1.0)
+    );
+
+
+    Info<< "Reading particle velocity field Us\n" << endl;
+    volVectorField Us
+    (
+        IOobject
+        (
+            "Us",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+//===============================

applications/solvers/cfdemSolverRhoPimple/pEqn.H

@@ -0,0 +1,92 @@
+rho = thermo.rho();
+rho = max(rho, rhoMin);
+rho = min(rho, rhoMax);
+rho.relax();
+
+volScalarField rAU(1.0/UEqn.A());
+surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rhoeps*rAU));
+if (modelType=="A")
+{
+    rhorAUf *= fvc::interpolate(voidfraction);
+}
+volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
+
+surfaceScalarField phiUs("phiUs", fvc::interpolate(rhoeps*rAU*Ksl*Us)& mesh.Sf());
+
+if (pimple.nCorrPISO() <= 1)
+{
+    tUEqn.clear();
+}
+
+if (pimple.transonic())
+{
+//  transonic version not implemented yet
+}
+else
+{
+    surfaceScalarField phiHbyA
+    (
+        "phiHbyA",
+        (
+            fvc::flux(rhoeps*HbyA)
+    //      + rhorAUf*fvc::ddtCorr(rho, U, phi)
+        )
+    );
+    
+    // flux without pressure gradient contribution
+    phi = phiHbyA + phiUs;
+
+    // Update the pressure BCs to ensure flux consistency
+    constrainPressure(p, rhoeps, U, phi, rhorAUf);
+
+    while (pimple.correctNonOrthogonal())
+    {
+        // Pressure corrector
+        fvScalarMatrix pEqn
+        (
+            fvm::ddt(psi*voidfraction, p)
+          + fvc::div(phi)
+          - fvm::laplacian(rhorAUf, p)
+          ==
+            fvOptions(psi, p, rho.name())
+        );
+
+        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
+
+        if (pimple.finalNonOrthogonalIter())
+        {
+            phi += pEqn.flux();
+        }
+    }
+}
+
+#include "rhoEqn.H"
+#include "compressibleContinuityErrsPU.H"
+
+// Explicitly relax pressure for momentum corrector
+p.relax();
+
+// Recalculate density from the relaxed pressure
+rho = thermo.rho();
+rho = max(rho, rhoMin);
+rho = min(rho, rhoMax);
+rho.relax();
+Info<< "rho max/min : " << max(rho).value()
+    << " " << min(rho).value() << endl;
+
+if (modelType=="A")
+{
+    U = HbyA - rAU*(voidfraction*fvc::grad(p)-Ksl*Us);
+}
+else
+{
+    U = HbyA - rAU*(fvc::grad(p)-Ksl*Us);
+}
+U.correctBoundaryConditions();
+fvOptions.correct(U);
+K = 0.5*magSqr(U);
+
+if (thermo.dpdt())
+{
+    dpdt = fvc::ddt(voidfraction,p);
+}

applications/solvers/cfdemSolverRhoPimple/rhoEqn.H

@@ -0,0 +1,17 @@
+{
+    fvScalarMatrix rhoEqn
+    (
+        fvm::ddt(voidfraction,rho)
+      + fvc::div(phi)
+      ==
+        fvOptions(rho)
+    );
+
+    fvOptions.constrain(rhoEqn);
+
+    rhoEqn.solve();
+
+    fvOptions.correct(rho);
+}
+
+// ************************************************************************* //

applications/utilities/cfdemPostproc/Make/options

@@ -1,5 +1,8 @@
+include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
+
 EXE_INC = \
-    -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
@@ -8,10 +11,12 @@ EXE_INC = \
 
 
 EXE_LIBS = \
-    -L$(FOAM_USER_LIBBIN)\
-    -lincompressibleRASModels \
-    -lincompressibleLESModels \
+    -L$(CFDEM_LIB_DIR)\
+    -lturbulenceModels \
+    -lincompressibleTurbulenceModels \
     -lincompressibleTransportModels \
     -lfiniteVolume \
     -l$(CFDEM_LIB_NAME) \
+     $(CFDEM_ADD_LIB_PATHS) \
+     $(CFDEM_ADD_LIBS)
 

applications/utilities/cfdemPostproc/cfdemPostproc.C

@@ -33,7 +33,7 @@ Description
 
 #include "fvCFD.H"
 #include "singlePhaseTransportModel.H"
-#include "turbulenceModel.H"
+#include "turbulentTransportModel.H"
 
 #include "cfdemCloud.H"
 #include "dataExchangeModel.H"
@@ -71,8 +71,9 @@ int main(int argc, char *argv[])
     double **voidfractions_;
     double **particleWeights_;
     double **particleVolumes_;
-    double **cellIDs_;
-    
+    double **particleV_;
+    int **cellIDs_;
+
     particleCloud.dataExchangeM().allocateArray(positions_,0.,3);
     particleCloud.dataExchangeM().allocateArray(velocities_,0.,3);
     particleCloud.get_radii(radii_);  // get ref to radii
@@ -80,6 +81,7 @@ int main(int argc, char *argv[])
     particleCloud.dataExchangeM().allocateArray(voidfractions_,0.,1);
     particleCloud.dataExchangeM().allocateArray(particleWeights_,0.,1);
     particleCloud.dataExchangeM().allocateArray(particleVolumes_,0.,1);
+    particleCloud.dataExchangeM().allocateArray(particleV_,0.,1);
     particleCloud.get_cellIDs(cellIDs_);  // get ref to cellIDs
     //particleCloud.dataExchangeM().allocateArray(cellIDs_,0.,1);
     
@@ -96,7 +98,7 @@ int main(int argc, char *argv[])
         particleCloud.averagingM().resetWeightFields();
         particleCloud.momCoupleM(0).resetMomSourceField();
 
-        particleCloud.dataExchangeM().couple();
+        particleCloud.dataExchangeM().couple(0);
 
         particleCloud.dataExchangeM().getData("x","vector-atom",positions_,count);
         particleCloud.dataExchangeM().getData("v","vector-atom",velocities_,count);
@@ -105,9 +107,9 @@ int main(int argc, char *argv[])
         particleCloud.locateM().findCell(NULL,positions_,cellIDs_,particleCloud.numberOfParticles());
         particleCloud.setPos(positions_);
 
-        particleCloud.voidFractionM().setvoidFraction(NULL,voidfractions_,particleWeights_,particleVolumes_);
+        particleCloud.voidFractionM().setvoidFraction(NULL,voidfractions_,particleWeights_,particleVolumes_,particleV_);
 
-        voidfraction.internalField() = particleCloud.voidFractionM().voidFractionInterp();
+        voidfraction.ref() = particleCloud.voidFractionM().voidFractionInterp();
         voidfraction.correctBoundaryConditions();
 
         particleCloud.averagingM().setVectorAverage
@@ -135,6 +137,7 @@ int main(int argc, char *argv[])
     particleCloud.dataExchangeM().destroy(voidfractions_,1);
     particleCloud.dataExchangeM().destroy(particleWeights_,1);
     particleCloud.dataExchangeM().destroy(particleVolumes_,1);
+    particleCloud.dataExchangeM().destroy(particleV_,1);
     //particleCloud.dataExchangeM().destroy(cellIDs_); // destroyed in cloud
 
     Info<< "End\n" << endl;

applications/utilities/cfdemPostproc/createFields.H

@@ -10,7 +10,7 @@
             IOobject::NO_WRITE
         ),
         mesh,
-        dimensionedScalar("0", dimensionSet(0, 2, -2, 0, 0), 1.0)
+        dimensionedScalar("1", dimensionSet(0, 2, -2, 0, 0), 1.0)
     );
 
     Info<< "Reading physical velocity field U" << endl;
@@ -26,15 +26,29 @@
             IOobject::NO_WRITE
         ),
         mesh,
-        dimensionedVector("0", dimensionSet(0, 1, -1, 0, 0), vector::zero)
+        dimensionedVector("0", dimensionSet(0, 1, -1, 0, 0), Foam::vector::zero)
     );
 
 //========================
 // drag law modelling
 //========================
 
+    Info<< "Creating dummy density field rho\n" << endl;
+    volScalarField rho
+    (
+        IOobject
+        (
+            "rho",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("1", dimensionSet(1, -3, 0, 0, 0), 1.0)
+    );
 
-    Info<< "\nReading voidfraction field voidfraction = (Vgas/Vparticle)\n" << endl;
+    Info<< "Reading voidfraction field voidfraction = (Vgas/Vparticle)\n" << endl;
     volScalarField voidfraction
     (
         IOobject
@@ -46,7 +60,7 @@
             IOobject::AUTO_WRITE
         ),
         mesh,
-        dimensionedScalar("0", dimensionSet(0, 0, 0, 0, 0), 1.)
+        dimensionedScalar("1", dimensionSet(0, 0, 0, 0, 0), 1.0)
     );
 
 
@@ -62,7 +76,7 @@
             IOobject::AUTO_WRITE
         ),
         mesh,
-        dimensionedVector("0", dimensionSet(0, 1, -1, 0, 0), vector::zero)
+        dimensionedVector("0", dimensionSet(0, 1, -1, 0, 0), Foam::vector::zero)
     );
 
 //========================

doc/.gitignore

@@ -0,0 +1,7 @@
+*.o
+*.d
+*.a
+*.dep
+log_*
+log.*
+*~

doc/CFDEMcoupling_Manual.html

@@ -164,7 +164,7 @@ In order to get the latest code version, please use the git repository at http:/
 </P>
 <PRE>modelType 
 </PRE>
-<P>"modelType" refers to the formulation of the equations to be solved. Choose "A" or "B", according to Zhou et al. (2010): "Discrete particle simulation of particle-fluid flow: model formulations and their applicability", JFM. "A" requires the use of the force models gradPForce and viscForce, whereas "B" requires the force model "Archimedes".
+<P>"modelType" refers to the formulation of the equations to be solved. Choose "A", "B" or "Bfull", according to Zhou et al. (2010): "Discrete particle simulation of particle-fluid flow: model formulations and their applicability", JFM. "A" requires the use of the force models gradPForce and viscForce, whereas "B" requires the force model "Archimedes". "Bfull" refers to model type I.
 </P>
 <PRE>couplingInterval 
 </PRE>
@@ -212,30 +212,31 @@ listing below of styles within certain commands.
 <TR ALIGN="center"><TD ><A HREF = "clockModel_noClock.html">clockModel_noClock</A></TD><TD ><A HREF = "clockModel_standardClock.html">clockModel_standardClock</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "dataExchangeModel.html">dataExchangeModel</A></TD><TD ><A HREF = "dataExchangeModel_noDataExchange.html">dataExchangeModel_noDataExchange</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "dataExchangeModel_oneWayVTK.html">dataExchangeModel_oneWayVTK</A></TD><TD ><A HREF = "dataExchangeModel_twoWayFiles.html">dataExchangeModel_twoWayFiles</A></TD></TR>
-<TR ALIGN="center"><TD ><A HREF = "dataExchangeModel_twoWayM2M.html">dataExchangeModel_twoWayM2M</A></TD><TD ><A HREF = "dataExchangeModel_twoWayMPI.html">dataExchangeModel_twoWayMPI</A></TD></TR>
+<TR ALIGN="center"><TD ><A HREF = "dataExchangeModel_twoWayMPI.html">dataExchangeModel_twoWayMPI</A></TD><TD ><A HREF = "dataExchangeModel_twoWayMany2Many.html">dataExchangeModel_twoWayMany2Many</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "forceModel.html">forceModel</A></TD><TD ><A HREF = "forceModel_Archimedes.html">forceModel_Archimedes</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "forceModel_ArchimedesIB.html">forceModel_ArchimedesIB</A></TD><TD ><A HREF = "forceModel_DiFeliceDrag.html">forceModel_DiFeliceDrag</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "forceModel_GidaspowDrag.html">forceModel_GidaspowDrag</A></TD><TD ><A HREF = "forceModel_KochHillDrag.html">forceModel_KochHillDrag</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "forceModel_LaEuScalarTemp.html">forceModel_LaEuScalarTemp</A></TD><TD ><A HREF = "forceModel_MeiLift.html">forceModel_MeiLift</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "forceModel_SchillerNaumannDrag.html">forceModel_SchillerNaumannDrag</A></TD><TD ><A HREF = "forceModel_ShirgaonkarIB.html">forceModel_ShirgaonkarIB</A></TD></TR>
-<TR ALIGN="center"><TD ><A HREF = "forceModel_gradPForce.html">forceModel_gradPForce</A></TD><TD ><A HREF = "forceModel_noDrag.html">forceModel_noDrag</A></TD></TR>
-<TR ALIGN="center"><TD ><A HREF = "forceModel_particleCellVolume.html">forceModel_particleCellVolume</A></TD><TD ><A HREF = "forceModel_virtualMassForce.html">forceModel_virtualMassForce</A></TD></TR>
-<TR ALIGN="center"><TD ><A HREF = "forceModel_viscForce.html">forceModel_viscForce</A></TD><TD ><A HREF = "liggghtsCommandModel.html">liggghtsCommandModel</A></TD></TR>
+<TR ALIGN="center"><TD ><A HREF = "forceModel_fieldStore.html">forceModel_fieldStore</A></TD><TD ><A HREF = "forceModel_gradPForce.html">forceModel_gradPForce</A></TD></TR>
+<TR ALIGN="center"><TD ><A HREF = "forceModel_noDrag.html">forceModel_noDrag</A></TD><TD ><A HREF = "forceModel_particleCellVolume.html">forceModel_particleCellVolume</A></TD></TR>
+<TR ALIGN="center"><TD ><A HREF = "forceModel_virtualMassForce.html">forceModel_virtualMassForce</A></TD><TD ><A HREF = "forceModel_viscForce.html">forceModel_viscForce</A></TD></TR>
+<TR ALIGN="center"><TD ><A HREF = "forceSubModel.html">forceSubModel</A></TD><TD ><A HREF = "forceSubModel_ImEx.html">forceSubModel_ImEx</A></TD></TR>
+<TR ALIGN="center"><TD ><A HREF = "forceSubModel_ImExCorr.html">forceSubModel_ImExCorr</A></TD><TD ><A HREF = "liggghtsCommandModel.html">liggghtsCommandModel</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "liggghtsCommandModel_execute.html">liggghtsCommandModel_execute</A></TD><TD ><A HREF = "liggghtsCommandModel_readLiggghtsData.html">liggghtsCommandModel_readLiggghtsData</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "liggghtsCommandModel_runLiggghts.html">liggghtsCommandModel_runLiggghts</A></TD><TD ><A HREF = "liggghtsCommandModel_writeLiggghts.html">liggghtsCommandModel_writeLiggghts</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "locateModel.html">locateModel</A></TD><TD ><A HREF = "locateModel_engineSearch.html">locateModel_engineSearch</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "locateModel_engineSearchIB.html">locateModel_engineSearchIB</A></TD><TD ><A HREF = "locateModel_standardSearch.html">locateModel_standardSearch</A></TD></TR>
-<TR ALIGN="center"><TD ><A HREF = "locateModel_turboEngineSearchM2M.html">locateModel_turboEngineM2MSearch</A></TD><TD ><A HREF = "locateModel_turboEngineSearch.html">locateModel_turboEngineSearch</A></TD></TR>
-<TR ALIGN="center"><TD ><A HREF = "meshMotionModel.html">meshMotionModel</A></TD><TD ><A HREF = "meshMotionModel_noMeshMotion.html">meshMotionModel_noMeshMotion</A></TD></TR>
-<TR ALIGN="center"><TD ><A HREF = "momCoupleModel.html">momCoupleModel</A></TD><TD ><A HREF = "momCoupleModel_explicitCouple.html">momCoupleModel_explicitCouple</A></TD></TR>
-<TR ALIGN="center"><TD ><A HREF = "momCoupleModel_implicitCouple.html">momCoupleModel_implicitCouple</A></TD><TD ><A HREF = "momCoupleModel_noCouple.html">momCoupleModel_noCouple</A></TD></TR>
-<TR ALIGN="center"><TD ><A HREF = "probeModel.html">probeModel</A></TD><TD ><A HREF = "probeModel_noProbe.html">probeModel_noProbe</A></TD></TR>
-<TR ALIGN="center"><TD ><A HREF = "regionModel.html">regionModel</A></TD><TD ><A HREF = "regionModel_allRegion.html">regionModel_allRegion</A></TD></TR>
-<TR ALIGN="center"><TD ><A HREF = "smoothingModel.html">smoothingModel</A></TD><TD ><A HREF = "smoothingModel_constDiffSmoothing.html">smoothingModel_constDiffSmoothing</A></TD></TR>
-<TR ALIGN="center"><TD ><A HREF = "smoothingModel_noSmoothing.html">smoothingModel_noSmoothing</A></TD><TD ><A HREF = "voidFractionModel.html">voidfractionModel</A></TD></TR>
-<TR ALIGN="center"><TD ><A HREF = "voidFractionModel_GaussVoidFraction.html">voidfractionModel_GaussVoidFraction</A></TD><TD ><A HREF = "voidFractionModel_IBVoidFraction.html">voidfractionModel_IBVoidFraction</A></TD></TR>
-<TR ALIGN="center"><TD ><A HREF = "voidFractionModel_bigParticleVoidFraction.html">voidfractionModel_bigParticleVoidFraction</A></TD><TD ><A HREF = "voidFractionModel_centreVoidFraction.html">voidfractionModel_centreVoidFraction</A></TD></TR>
-<TR ALIGN="center"><TD ><A HREF = "voidFractionModel_dividedVoidFraction.html">voidfractionModel_dividedVoidFraction</A> 
+<TR ALIGN="center"><TD ><A HREF = "locateModel_turboEngineSearch.html">locateModel_turboEngineSearch</A></TD><TD ><A HREF = "meshMotionModel.html">meshMotionModel</A></TD></TR>
+<TR ALIGN="center"><TD ><A HREF = "meshMotionModel_noMeshMotion.html">meshMotionModel_noMeshMotion</A></TD><TD ><A HREF = "momCoupleModel.html">momCoupleModel</A></TD></TR>
+<TR ALIGN="center"><TD ><A HREF = "momCoupleModel_explicitCouple.html">momCoupleModel_explicitCouple</A></TD><TD ><A HREF = "momCoupleModel_implicitCouple.html">momCoupleModel_implicitCouple</A></TD></TR>
+<TR ALIGN="center"><TD ><A HREF = "momCoupleModel_noCouple.html">momCoupleModel_noCouple</A></TD><TD ><A HREF = "probeModel.html">probeModel</A></TD></TR>
+<TR ALIGN="center"><TD ><A HREF = "probeModel_noProbe.html">probeModel_noProbe</A></TD><TD ><A HREF = "regionModel.html">regionModel</A></TD></TR>
+<TR ALIGN="center"><TD ><A HREF = "regionModel_allRegion.html">regionModel_allRegion</A></TD><TD ><A HREF = "smoothingModel.html">smoothingModel</A></TD></TR>
+<TR ALIGN="center"><TD ><A HREF = "smoothingModel_constDiffSmoothing.html">smoothingModel_constDiffSmoothing</A></TD><TD ><A HREF = "smoothingModel_noSmoothing.html">smoothingModel_noSmoothing</A></TD></TR>
+<TR ALIGN="center"><TD ><A HREF = "voidFractionModel.html">voidfractionModel</A></TD><TD ><A HREF = "voidFractionModel_GaussVoidFraction.html">voidfractionModel_GaussVoidFraction</A></TD></TR>
+<TR ALIGN="center"><TD ><A HREF = "voidFractionModel_IBVoidFraction.html">voidfractionModel_IBVoidFraction</A></TD><TD ><A HREF = "voidFractionModel_bigParticleVoidFraction.html">voidfractionModel_bigParticleVoidFraction</A></TD></TR>
+<TR ALIGN="center"><TD ><A HREF = "voidFractionModel_centreVoidFraction.html">voidfractionModel_centreVoidFraction</A></TD><TD ><A HREF = "voidFractionModel_dividedVoidFraction.html">voidfractionModel_dividedVoidFraction</A> 
 </TD></TR></TABLE></DIV>
 
 </HTML>

doc/CFDEMcoupling_Manual.txt

@@ -133,7 +133,7 @@ Reasonable example settings for the "couplingProperties" dictionary are given in
 
 modelType :pre
 
-"modelType" refers to the formulation of the equations to be solved. Choose "A" or "B", according to Zhou et al. (2010): "Discrete particle simulation of particle-fluid flow: model formulations and their applicability", JFM. "A" requires the use of the force models gradPForce and viscForce, whereas "B" requires the force model "Archimedes".
+"modelType" refers to the formulation of the equations to be solved. Choose "A", "B" or "Bfull", according to Zhou et al. (2010): "Discrete particle simulation of particle-fluid flow: model formulations and their applicability", JFM. "A" requires the use of the force models gradPForce and viscForce, whereas "B" requires the force model "Archimedes". "Bfull" refers to model type I.
 
 couplingInterval :pre
 
@@ -221,6 +221,10 @@ Reasonable example settings for the "liggghtsCommands" dictionary are given in t
 
 
 
+
+
+
+
 
 
 
@@ -250,8 +254,8 @@ listing below of styles within certain commands.
 "dataExchangeModel_noDataExchange"_dataExchangeModel_noDataExchange.html,
 "dataExchangeModel_oneWayVTK"_dataExchangeModel_oneWayVTK.html,
 "dataExchangeModel_twoWayFiles"_dataExchangeModel_twoWayFiles.html,
-"dataExchangeModel_twoWayM2M"_dataExchangeModel_twoWayM2M.html,
 "dataExchangeModel_twoWayMPI"_dataExchangeModel_twoWayMPI.html,
+"dataExchangeModel_twoWayMany2Many"_dataExchangeModel_twoWayMany2Many.html,
 "forceModel"_forceModel.html,
 "forceModel_Archimedes"_forceModel_Archimedes.html,
 "forceModel_ArchimedesIB"_forceModel_ArchimedesIB.html,
@@ -262,11 +266,15 @@ listing below of styles within certain commands.
 "forceModel_MeiLift"_forceModel_MeiLift.html,
 "forceModel_SchillerNaumannDrag"_forceModel_SchillerNaumannDrag.html,
 "forceModel_ShirgaonkarIB"_forceModel_ShirgaonkarIB.html,
+"forceModel_fieldStore"_forceModel_fieldStore.html,
 "forceModel_gradPForce"_forceModel_gradPForce.html,
 "forceModel_noDrag"_forceModel_noDrag.html,
 "forceModel_particleCellVolume"_forceModel_particleCellVolume.html,
 "forceModel_virtualMassForce"_forceModel_virtualMassForce.html,
 "forceModel_viscForce"_forceModel_viscForce.html,
+"forceSubModel"_forceSubModel.html,
+"forceSubModel_ImEx"_forceSubModel_ImEx.html,
+"forceSubModel_ImExCorr"_forceSubModel_ImExCorr.html,
 "liggghtsCommandModel"_liggghtsCommandModel.html,
 "liggghtsCommandModel_execute"_liggghtsCommandModel_execute.html,
 "liggghtsCommandModel_readLiggghtsData"_liggghtsCommandModel_readLiggghtsData.html,
@@ -276,7 +284,6 @@ listing below of styles within certain commands.
 "locateModel_engineSearch"_locateModel_engineSearch.html,
 "locateModel_engineSearchIB"_locateModel_engineSearchIB.html,
 "locateModel_standardSearch"_locateModel_standardSearch.html,
-"locateModel_turboEngineM2MSearch"_locateModel_turboEngineSearchM2M.html,
 "locateModel_turboEngineSearch"_locateModel_turboEngineSearch.html,
 "meshMotionModel"_meshMotionModel.html,
 "meshMotionModel_noMeshMotion"_meshMotionModel_noMeshMotion.html,

doc/cfdemSolverPisoScalar.html

@@ -11,12 +11,14 @@
 </H3>
 <P><B>Description:</B>
 </P>
-<P>"cfdemSolverPisoScalar" is a coupled CFD-DEM solver using CFDEMcoupling, an open source parallel coupled CFD-DEM framework. Based on pisoFoam(R)(*), a finite volume based solver for turbulent Navier-Stokes equations applying PISO algorithm, "cfdemSolverPisoScalar" has additional functionality for a coupling to the DEM code "LIGGGHTS" as well as a scalar transport equation. The volume averaged Navier-Stokes Equations are solved accounting for momentum exchange and volume displacement of discrete particles whose trajectories are calculated in the DEM code LIGGGHTS. The scalar transport equation is coupled to scalar properties of the particle phase, thus convective heat transfer in a fluid granular system can be modeled with "cfdemSolverPisoScalar".
+<P>"cfdemSolverPisoScalar" is a coupled CFD-DEM solver using CFDEMcoupling, an open source parallel coupled CFD-DEM framework. Based on pisoFoam(R)(*), a finite volume based solver for turbulent Navier-Stokes equations applying PISO algorithm, "cfdemSolverPisoScalar" has additional functionality for a coupling to the DEM code "LIGGGHTS" as well as a scalar transport equation. The volume averaged Navier-Stokes Equations are solved accounting for momentum exchange and volume displacement of discrete particles, whose trajectories are calculated in the DEM code LIGGGHTS. The scalar transport equation is coupled to scalar properties of the particle phase, thus convective heat transfer in a fluid granular system can be modeled with "cfdemSolverPisoScalar". 
 </P>
 <P>see:
 </P>
 <P>GONIVA, C., KLOSS, C., HAGER,A. and PIRKER, S. (2010): "An Open Source CFD-DEM Perspective", Proc. of OpenFOAM Workshop, Göteborg, June 22.-24.
 </P>
+<P>The heat transfer equation is implemented according to Nield & Bejan (2013), Convection in Porous Media, DOI 10.1007/978-1-4614-5541-7_2, Springer
+</P>
 <HR>
 
 <P>(*) This offering is not approved or endorsed by OpenCFD Limited, the producer of the OpenFOAM software and owner of the OPENFOAM®  and OpenCFD®  trade marks. OPENFOAM® is a registered trade mark of OpenCFD Limited, a wholly owned subsidiary of the ESI Group.

doc/cfdemSolverPisoScalar.txt

@@ -9,12 +9,15 @@ cfdemSolverPisoScalar command :h3
 
 [Description:]
 
-"cfdemSolverPisoScalar" is a coupled CFD-DEM solver using CFDEMcoupling, an open source parallel coupled CFD-DEM framework. Based on pisoFoam(R)(*), a finite volume based solver for turbulent Navier-Stokes equations applying PISO algorithm, "cfdemSolverPisoScalar" has additional functionality for a coupling to the DEM code "LIGGGHTS" as well as a scalar transport equation. The volume averaged Navier-Stokes Equations are solved accounting for momentum exchange and volume displacement of discrete particles whose trajectories are calculated in the DEM code LIGGGHTS. The scalar transport equation is coupled to scalar properties of the particle phase, thus convective heat transfer in a fluid granular system can be modeled with "cfdemSolverPisoScalar".
+"cfdemSolverPisoScalar" is a coupled CFD-DEM solver using CFDEMcoupling, an open source parallel coupled CFD-DEM framework. Based on pisoFoam(R)(*), a finite volume based solver for turbulent Navier-Stokes equations applying PISO algorithm, "cfdemSolverPisoScalar" has additional functionality for a coupling to the DEM code "LIGGGHTS" as well as a scalar transport equation. The volume averaged Navier-Stokes Equations are solved accounting for momentum exchange and volume displacement of discrete particles, whose trajectories are calculated in the DEM code LIGGGHTS. The scalar transport equation is coupled to scalar properties of the particle phase, thus convective heat transfer in a fluid granular system can be modeled with "cfdemSolverPisoScalar". 
 
 see:
 
 GONIVA, C., KLOSS, C., HAGER,A. and PIRKER, S. (2010): "An Open Source CFD-DEM Perspective", Proc. of OpenFOAM Workshop, Göteborg, June 22.-24.
 
+The heat transfer equation is implemented according to Nield & Bejan (2013), Convection in Porous Media, DOI 10.1007/978-1-4614-5541-7_2, Springer
+
+
 
 :line
 (*) This offering is not approved or endorsed by OpenCFD Limited, the producer of the OpenFOAM software and owner of the OPENFOAM®  and OpenCFD®  trade marks. OPENFOAM® is a registered trade mark of OpenCFD Limited, a wholly owned subsidiary of the ESI Group.

doc/clockModel_standardClock.html

@@ -23,7 +23,7 @@
 </P>
 <P>The "standardClock" model is a basic clockModel model which measures the run time between every ".start(int arrayPos,string name)" and ".stop(string name)" statement placed in the code. If a ".start(name)" is called more than once (e.g. in a loop) the accumulated times are calculated.  After the simulation has finished, the data is stored in $caseDir/CFD/clockData/$startTime/*.txt .
 Since the measurements are stored in an array, it is necessary to put a variable <I>arrayPos</I> (type integer) at the start command. Those do not need to be in ascending order and positions may be omitted. The standard size of this array is 30 and can be changed at the initialization of the standardClock class. If <I>arrayPos</I> is out of bounds, the array size will be doubled. The stop command does not need <I>arrayPos</I>, since the class remembers the positions. The string name is intended for easier evaluation afterwards an may be omitted like ".start(int arrayPos)" and ".stop()". The command ".stop(string name)" is a safety feature, because if the name is not equal to the started name, output will be produced for information.
-After the case ran you may use the matPlot.py script located in $CFDEM_UT_DIR/vizClock/ to produce a graphical output of your measurements. The usage is like 'python < matPlot.py' and you have to be in the directory of the desired time step, where there is a file called "timeEvalFull.txt", which contains averaged and maximum data with respect to the number of processes.
+After the case ran you may use the matPlot.py script located in $CFDEM_UT_DIR/vizClock/ to produce a graphical output of your measurements. The usage is like 'python < matPlot.py' and you have to be in the directory of the desired time step, where there is a file called "timeEvalFull.txt", which contains averaged and maximum data with respect to the number of processes. There is an alias called "vizClock" to run this python routine for visualizing the data.
 </P>
 <P><B>Restrictions:</B> none.
 </P>

doc/clockModel_standardClock.txt

@@ -21,7 +21,7 @@ clockModel standardClock; :pre
 
 The "standardClock" model is a basic clockModel model which measures the run time between every ".start(int arrayPos,string name)" and ".stop(string name)" statement placed in the code. If a ".start(name)" is called more than once (e.g. in a loop) the accumulated times are calculated.  After the simulation has finished, the data is stored in $caseDir/CFD/clockData/$startTime/*.txt .
 Since the measurements are stored in an array, it is necessary to put a variable {arrayPos} (type integer) at the start command. Those do not need to be in ascending order and positions may be omitted. The standard size of this array is 30 and can be changed at the initialization of the standardClock class. If {arrayPos} is out of bounds, the array size will be doubled. The stop command does not need {arrayPos}, since the class remembers the positions. The string name is intended for easier evaluation afterwards an may be omitted like ".start(int arrayPos)" and ".stop()". The command ".stop(string name)" is a safety feature, because if the name is not equal to the started name, output will be produced for information.
-After the case ran you may use the matPlot.py script located in $CFDEM_UT_DIR/vizClock/ to produce a graphical output of your measurements. The usage is like 'python < matPlot.py' and you have to be in the directory of the desired time step, where there is a file called "timeEvalFull.txt", which contains averaged and maximum data with respect to the number of processes.
+After the case ran you may use the matPlot.py script located in $CFDEM_UT_DIR/vizClock/ to produce a graphical output of your measurements. The usage is like 'python < matPlot.py' and you have to be in the directory of the desired time step, where there is a file called "timeEvalFull.txt", which contains averaged and maximum data with respect to the number of processes. There is an alias called "vizClock" to run this python routine for visualizing the data.
 
 [Restrictions:] none.
 

doc/dataExchangeModel_twoWayM2M.txt

@@ -1,42 +0,0 @@
-"CFDEMproject WWW Site"_lws - "CFDEM Commands"_lc :c
-
-:link(lws,http://www.cfdem.com)
-:link(lc,CFDEMcoupling_Manual.html#comm)
-
-:line
-
-dataExchangeModel_twoWayM2M command :h3
-
-[Syntax:]
-
-Defined in couplingProperties dictionary.
-
-dataExchangeModel twoWayM2M;
-twoWayM2MProps
-\{
-    liggghtsPath "path";
-\}; :pre
-
-{path} = path to the DEM simulation input file :ulb,l
-:ule
-
-[Examples:]
-
-dataExchangeModel twoWayM2M;
-twoWayM2MProps
-\{
-    liggghtsPath "../DEM/in.liggghts_init";
-\} :pre
-
-[Description:]
-
-The data exchange model performs the data exchange between the DEM code and the CFD code. The twoWayM2M model is a model that can exchange particle properties from DEM to CFD and from CFD to DEM. Data is exchanged via MPI technique using the many to many mapping scheme. The DEM run is executed by the coupling model, via a liggghtsCommandModel object.
-
-[Restrictions:]
-
-Should be used in combination with the turboEngineSearchM2M locate model to achieve best performance!
-
-[Related commands:]
-
-"dataExchangeModel"_dataExchangeModel.html
-

doc/dataExchangeModel_twoWayMany2Many.html

@@ -7,14 +7,14 @@
 
 <HR>
 
-<H3>dataExchangeModel_twoWayM2M command 
+<H3>dataExchangeModel_twoWayMany2Many command 
 </H3>
 <P><B>Syntax:</B>
 </P>
 <P>Defined in couplingProperties dictionary.
 </P>
-<PRE>dataExchangeModel twoWayM2M;
-twoWayM2MProps
+<PRE>dataExchangeModel twoWayMany2Many;
+twoWayMany2ManyProps
 {
     liggghtsPath "path";
 }; 
@@ -25,19 +25,19 @@ twoWayM2MProps
 </UL>
 <P><B>Examples:</B>
 </P>
-<PRE>dataExchangeModel twoWayM2M;
-twoWayM2MProps
+<PRE>dataExchangeModel twoWayMany2Many;
+twoWayMany2ManyProps
 {
     liggghtsPath "../DEM/in.liggghts_init";
 } 
 </PRE>
 <P><B>Description:</B>
 </P>
-<P>The data exchange model performs the data exchange between the DEM code and the CFD code. The twoWayM2M model is a model that can exchange particle properties from DEM to CFD and from CFD to DEM. Data is exchanged via MPI technique using the many to many mapping scheme. The DEM run is executed by the coupling model, via a liggghtsCommandModel object.
+<P>The data exchange model performs the data exchange between the DEM code and the CFD code. The twoWayMany2Many model is a model that can exchange particle properties from DEM to CFD and from CFD to DEM. Data is exchanged via MPI technique using the many to many mapping scheme. The DEM run is executed by the coupling model, via a liggghtsCommandModel object.
 </P>
 <P><B>Restrictions:</B>
 </P>
-<P>Should be used in combination with the turboEngineSearchM2M locate model to achieve best performance!
+<P>Must be used in combination with the engineSearchMany2Many locate model!
 </P>
 <P><B>Related commands:</B>
 </P>

doc/dataExchangeModel_twoWayMany2Many.txt

@@ -0,0 +1,42 @@
+"CFDEMproject WWW Site"_lws - "CFDEM Commands"_lc :c
+
+:link(lws,http://www.cfdem.com)
+:link(lc,CFDEMcoupling_Manual.html#comm)
+
+:line
+
+dataExchangeModel_twoWayMany2Many command :h3
+
+[Syntax:]
+
+Defined in couplingProperties dictionary.
+
+dataExchangeModel twoWayMany2Many;
+twoWayMany2ManyProps
+\{
+    liggghtsPath "path";
+\}; :pre
+
+{path} = path to the DEM simulation input file :ulb,l
+:ule
+
+[Examples:]
+
+dataExchangeModel twoWayMany2Many;
+twoWayMany2ManyProps
+\{
+    liggghtsPath "../DEM/in.liggghts_init";
+\} :pre
+
+[Description:]
+
+The data exchange model performs the data exchange between the DEM code and the CFD code. The twoWayMany2Many model is a model that can exchange particle properties from DEM to CFD and from CFD to DEM. Data is exchanged via MPI technique using the many to many mapping scheme. The DEM run is executed by the coupling model, via a liggghtsCommandModel object.
+
+[Restrictions:]
+
+Must be used in combination with the engineSearchMany2Many locate model!
+
+[Related commands:]
+
+"dataExchangeModel"_dataExchangeModel.html
+

doc/forceModel.html

@@ -33,7 +33,7 @@
 </P>
 <P><B>Description:</B>
 </P>
-<P>The force model performs the calculation of forces (e.g. fluid-particle interaction forces) acting on each DEM particle. All force models selected are executed sequentially and the forces on the particles are superposed.
+<P>The force model performs the calculation of forces (e.g. fluid-particle interaction forces) acting on each DEM particle. All force models selected are executed sequentially and the forces on the particles are superposed. If the fluid density field is needed, by default a field named "rho" will be used. Via the forceSubModel an alternative field can be chosen.
 </P>
 <P><B>Restrictions:</B>
 </P>

doc/forceModel.txt

@@ -31,7 +31,7 @@ Note: This examples list might not be complete - please look for other models (f
 
 [Description:]
 
-The force model performs the calculation of forces (e.g. fluid-particle interaction forces) acting on each DEM particle. All force models selected are executed sequentially and the forces on the particles are superposed.
+The force model performs the calculation of forces (e.g. fluid-particle interaction forces) acting on each DEM particle. All force models selected are executed sequentially and the forces on the particles are superposed. If the fluid density field is needed, by default a field named "rho" will be used. Via the forceSubModel an alternative field can be chosen.
 
 [Restrictions:]
 

doc/forceModel_Archimedes.html

@@ -19,13 +19,10 @@
 );
 ArchimedesProps
 {
-    densityFieldName "density";
     gravityFieldName "gravity";
 }; 
 </PRE>
-<UL><LI><I>density</I> = name of the finite volume density field 
-
-<LI><I>gravity</I> = name of the finite volume gravity field 
+<UL><LI><I>gravity</I> = name of the finite volume gravity field 
 
 
 </UL>
@@ -37,7 +34,6 @@ ArchimedesProps
 );
 ArchimedesProps
 {
-    densityFieldName "rho";
     gravityFieldName "g";
 } 
 </PRE>

doc/forceModel_Archimedes.txt

@@ -17,12 +17,10 @@ forceModels
 );
 ArchimedesProps
 \{
-    densityFieldName "density";
     gravityFieldName "gravity";
 \}; :pre
 
-{density} = name of the finite volume density field :ulb,l
-{gravity} = name of the finite volume gravity field :l
+{gravity} = name of the finite volume gravity field :ulb,l
 :ule
 
 [Examples:]
@@ -33,7 +31,6 @@ forceModels
 );
 ArchimedesProps
 \{
-    densityFieldName "rho";
     gravityFieldName "g";
 \} :pre
 

doc/forceModel_ArchimedesIB.html

@@ -19,14 +19,11 @@
 );
 ArchimedesIBProps
 {
-    densityFieldName "density";
     gravityFieldName "gravity";
     voidfractionFieldName "voidfraction";
 }; 
 </PRE>
-<UL><LI><I>density</I> = name of the finite volume density field 
-
-<LI><I>gravity</I> = name of the finite volume gravity field 
+<UL><LI><I>gravity</I> = name of the finite volume gravity field 
 
 <LI><I>voidfraction</I> = name of the finite volume voidfraction field 
 
@@ -40,7 +37,6 @@ ArchimedesIBProps
 );
 ArchimedesIBProps
 {
-    densityFieldName "rho";
     gravityFieldName "g";
     voidfractionFieldName "voidfractionNext";
 } 

doc/forceModel_ArchimedesIB.txt

@@ -17,13 +17,11 @@ forceModels
 );
 ArchimedesIBProps
 \{
-    densityFieldName "density";
     gravityFieldName "gravity";
     voidfractionFieldName "voidfraction";
 \}; :pre
 
-{density} = name of the finite volume density field :ulb,l
-{gravity} = name of the finite volume gravity field :l
+{gravity} = name of the finite volume gravity field :ulb,l
 {voidfraction} = name of the finite volume voidfraction field :l
 :ule
 
@@ -35,7 +33,6 @@ forceModels
 );
 ArchimedesIBProps
 \{
-    densityFieldName "rho";
     gravityFieldName "g";
     voidfractionFieldName "voidfractionNext";
 \} :pre

doc/forceModel_DiFeliceDrag.html

@@ -20,15 +20,12 @@
 DiFeliceDragProps
 {
     velFieldName "U";
-    densityFieldName "density";
-    interpolation;
+    interpolation switch1;
 }; 
 </PRE>
 <UL><LI><I>U</I> = name of the finite volume fluid velocity field 
 
-<LI><I>density</I> = name of the finite volume gravity field 
-
-<LI><I>interpolation</I> = flag to use interpolated voidfraction and velocity values (normally off) 
+<LI><I>switch1</I> = flag to use interpolated voidfraction and velocity values (normally off) 
 
 
 </UL>
@@ -41,8 +38,7 @@ DiFeliceDragProps
 DiFeliceDragProps
 {
     velFieldName "U";
-    densityFieldName "rho";
-    interpolation;
+    interpolation true;
 } 
 </PRE>
 <P><B>Description:</B>

doc/forceModel_DiFeliceDrag.txt

@@ -18,13 +18,11 @@ forceModels
 DiFeliceDragProps
 \{
     velFieldName "U";
-    densityFieldName "density";
-    interpolation;
+    interpolation switch1;
 \}; :pre
 
 {U} = name of the finite volume fluid velocity field :ulb,l
-{density} = name of the finite volume gravity field :l
-{interpolation} = flag to use interpolated voidfraction and velocity values (normally off) :l
+{switch1} = flag to use interpolated voidfraction and velocity values (normally off) :l
 :ule
 
 [Examples:]
@@ -36,8 +34,7 @@ forceModels
 DiFeliceDragProps
 \{
     velFieldName "U";
-    densityFieldName "rho";
-    interpolation;
+    interpolation true;
 \} :pre
 
 [Description:]

doc/forceModel_GidaspowDrag.html

@@ -20,24 +20,24 @@
 GidaspowDragProps
 {
     velFieldName "U";
-    densityFieldName "density";
     voidfractionFieldName "voidfraction";
+    granVelFieldName "Us";
     phi "scalar";
-    interpolation;
-    implDEM;
+    interpolation switch1;
+    implForceDEM switch2;
 }; 
 </PRE>
 <UL><LI><I>U</I> = name of the finite volume fluid velocity field 
 
-<LI><I>density</I> = name of the finite volume gravity field 
-
 <LI><I>voidfraction</I> = name of the finite volume voidfraction field 
 
+<LI><I>Us</I> = name of the finite volume cell averaged particle velocity field 
+
 <LI><I>phi</I> = drag correction factor (in doubt 1) 
 
-<LI><I>interpolation</I> = (optional, normally off) flag to use interpolated voidfraction and fluid velocity values 
+<LI><I>switch1</I> = (optional, normally off) flag to use interpolated voidfraction and fluid velocity values 
 
-<I>implDEM</I> = (optional, normally off) flag to use implicit formulation of drag on DEM side:l
+<I>switch2</I> = (optional, normally off) flag to use implicit formulation of drag on DEM side:l
 
 </UL>
 <P><B>Examples:</B>
@@ -49,8 +49,8 @@ GidaspowDragProps
 GidaspowDragProps
 {
     velFieldName "U";
-    densityFieldName "rho";
     voidfractionFieldName "voidfraction";
+    granVelFieldName "Us";
 } 
 </PRE>
 <P><B>Description:</B>

doc/forceModel_GidaspowDrag.txt

@@ -18,19 +18,19 @@ forceModels
 GidaspowDragProps
 \{
     velFieldName "U";
-    densityFieldName "density";
     voidfractionFieldName "voidfraction";
+    granVelFieldName "Us";
     phi "scalar";
-    interpolation;
-    implDEM;
+    interpolation switch1;
+    implForceDEM switch2;
 \}; :pre
 
 {U} = name of the finite volume fluid velocity field :ulb,l
-{density} = name of the finite volume gravity field :l
 {voidfraction} = name of the finite volume voidfraction field :l
+{Us} = name of the finite volume cell averaged particle velocity field :l
 {phi} = drag correction factor (in doubt 1) :l
-{interpolation} = (optional, normally off) flag to use interpolated voidfraction and fluid velocity values :l
-{implDEM} = (optional, normally off) flag to use implicit formulation of drag on DEM side:l
+{switch1} = (optional, normally off) flag to use interpolated voidfraction and fluid velocity values :l
+{switch2} = (optional, normally off) flag to use implicit formulation of drag on DEM side:l
 :ule
 
 [Examples:]
@@ -42,8 +42,8 @@ forceModels
 GidaspowDragProps
 \{
     velFieldName "U";
-    densityFieldName "rho";
     voidfractionFieldName "voidfraction";
+    granVelFieldName "Us";
 \} :pre
 
 [Description:]

doc/forceModel_KochHillDrag.html

@@ -20,21 +20,18 @@
 KochHillDragProps
 {
     velFieldName "U";
-    densityFieldName "density";
     voidfractionFieldName "voidfraction"; 
-    interpolation;
-    implDEM;
+    interpolation "bool1";
+    implForceDEM "bool2";
 }; 
 </PRE>
 <UL><LI><I>U</I> = name of the finite volume fluid velocity field 
 
-<LI><I>density</I> = name of the finite volume gravity field 
-
 <LI><I>voidfraction</I> = name of the finite volume voidfraction field 
 
-<LI><I>interpolation</I> = (optional, normally off) flag to use interpolated voidfraction and fluid velocity values 
+<LI><I>bool1</I> = (optional, normally off) flag to use interpolated voidfraction and fluid velocity values 
 
-<I>implDEM</I> = (optional, normally off) flag to use implicit formulation of drag on DEM side:l
+<I>bool2</I> = (optional, normally off) flag to use implicit formulation of drag on DEM side:l
 
 </UL>
 <P><B>Examples:</B>
@@ -46,7 +43,6 @@ KochHillDragProps
 KochHillDragProps
 {
     velFieldName "U";
-    densityFieldName "rho";
     voidfractionFieldName "voidfraction"; 
 } 
 </PRE>

doc/forceModel_KochHillDrag.txt

@@ -18,17 +18,15 @@ forceModels
 KochHillDragProps
 \{
     velFieldName "U";
-    densityFieldName "density";
     voidfractionFieldName "voidfraction"; 
-    interpolation;
-    implDEM;
+    interpolation "bool1";
+    implForceDEM "bool2";
 \}; :pre
 
 {U} = name of the finite volume fluid velocity field :ulb,l
-{density} = name of the finite volume gravity field :l
 {voidfraction} = name of the finite volume voidfraction field :l
-{interpolation} = (optional, normally off) flag to use interpolated voidfraction and fluid velocity values :l
-{implDEM} = (optional, normally off) flag to use implicit formulation of drag on DEM side:l
+{bool1} = (optional, normally off) flag to use interpolated voidfraction and fluid velocity values :l
+{bool2} = (optional, normally off) flag to use implicit formulation of drag on DEM side:l
 :ule
 
 [Examples:]
@@ -40,7 +38,6 @@ forceModels
 KochHillDragProps
 \{
     velFieldName "U";
-    densityFieldName "rho";
     voidfractionFieldName "voidfraction"; 
 \} :pre
 

doc/forceModel_LaEuScalarTemp.html

@@ -21,23 +21,19 @@ LaEuScalarTempProps
 {
     velFieldName "U";
     tempFieldName "T";
-    tempSourceFieldName "Tsource";
     voidfractionFieldName "voidfraction";
     partTempName "Temp";
     partHeatFluxName "convectiveHeatFlux";
     lambda value;
     Cp value1;
-    densityFieldName "density";
-    interpolation;
-    verbose;
+    interpolation "switch1";
+    verbose "switch2";
 }; 
 </PRE>
 <UL><LI><I>U</I> = name of the finite volume fluid velocity field 
 
 <LI><I>T</I> = name of the finite volume scalar temperature field 
 
-<LI><I>Tsource</I> = name of the finite volume scalar temperature source field 
-
 <LI><I>voidfraction</I> = name of the finite volume voidfraction field 
 
 <LI><I>Temp</I> = name of the DEM data representing the particles temperature 
@@ -48,11 +44,9 @@ LaEuScalarTempProps
 
 <LI><I>value1</I> = fluid specific heat capacity [W*s/(kg*K)] 
 
-<LI><I>density</I> = name of the finite volume fluid density field 
+<LI><I>switch1</I> = (optional, normally off) flag to use interpolated voidfraction and fluid velocity values 
 
-<LI><I>interpolation</I> = (optional, normally off) flag to use interpolated voidfraction and fluid velocity values 
-
-<LI><I>verbose</I> = (normally off) for verbose run 
+<LI><I>switch2</I> = (normally off) for verbose run 
 
 
 </UL>
@@ -66,13 +60,11 @@ LaEuScalarTempProps
 {
     velFieldName "U";
     tempFieldName "T";
-    tempSourceFieldName "Tsource";
     voidfractionFieldName "voidfraction";
     partTempName "Temp";
     partHeatFluxName "convectiveHeatFlux";
     lambda 0.0256;
     Cp 1007;
-    densityFieldName "rho";
 } 
 </PRE>
 <P><B>Description:</B>
@@ -81,7 +73,7 @@ LaEuScalarTempProps
 </P>
 <P><B>Restrictions:</B>
 </P>
-<P>Goes only with cfdemSolverScalar.
+<P>Goes only with cfdemSolverScalar. The force model has to be the second (!!!) model in the forces list.
 </P>
 <P><B>Related commands:</B>
 </P>

doc/forceModel_LaEuScalarTemp.txt

@@ -19,28 +19,24 @@ LaEuScalarTempProps
 \{
     velFieldName "U";
     tempFieldName "T";
-    tempSourceFieldName "Tsource";
     voidfractionFieldName "voidfraction";
     partTempName "Temp";
     partHeatFluxName "convectiveHeatFlux";
     lambda value;
     Cp value1;
-    densityFieldName "density";
-    interpolation;
-    verbose;
+    interpolation "switch1";
+    verbose "switch2";
 \}; :pre
 
 {U} = name of the finite volume fluid velocity field :ulb,l
 {T} = name of the finite volume scalar temperature field :l
-{Tsource} = name of the finite volume scalar temperature source field :l
 {voidfraction} = name of the finite volume voidfraction field :l
 {Temp} = name of the DEM data representing the particles temperature :l
 {convectiveHeatFlux} = name of the DEM data representing the particle-fluid convective heat flux :l
 {value} = fluid thermal conductivity \[W/(m*K)\] :l
 {value1} = fluid specific heat capacity \[W*s/(kg*K)\] :l
-{density} = name of the finite volume fluid density field :l
-{interpolation} = (optional, normally off) flag to use interpolated voidfraction and fluid velocity values :l
-{verbose} = (normally off) for verbose run :l
+{switch1} = (optional, normally off) flag to use interpolated voidfraction and fluid velocity values :l
+{switch2} = (normally off) for verbose run :l
 :ule
 
 [Examples:]
@@ -53,13 +49,11 @@ LaEuScalarTempProps
 \{
     velFieldName "U";
     tempFieldName "T";
-    tempSourceFieldName "Tsource";
     voidfractionFieldName "voidfraction";
     partTempName "Temp";
     partHeatFluxName "convectiveHeatFlux";
     lambda 0.0256;
     Cp 1007;
-    densityFieldName "rho";
 \} :pre
 
 [Description:]
@@ -68,7 +62,7 @@ This "forceModel" does not influence the particles or the fluid flow! Using the
 
 [Restrictions:]
 
-Goes only with cfdemSolverScalar.
+Goes only with cfdemSolverScalar. The force model has to be the second (!!!) model in the forces list.
 
 [Related commands:]
 

doc/forceModel_MeiLift.html

@@ -20,21 +20,18 @@
 MeiLiftProps
 {
     velFieldName "U";
-    densityFieldName "density";
     useSecondOrderTerms;
-    interpolation;
-    verbose;
+    interpolation "switch1";
+    verbose "switch2";
 }; 
 </PRE>
 <UL><LI><I>U</I> = name of the finite volume fluid velocity field 
 
-<LI><I>density</I> = name of the finite volume fluid density field 
-
 <LI><I>useSecondOrderTerms</I> = switch to activate second order terms in the lift force model 
 
-<LI><I>interpolation</I> = switch to activate tri-linear interpolation of the flow quantities at the particle position 
+<LI><I>switch1</I> = switch to activate tri-linear interpolation of the flow quantities at the particle position 
 
-<LI><I>verbose</I> = switch to activate the report of per-particle quantities to the screen 
+<LI><I>switch2</I> = switch to activate the report of per-particle quantities to the screen 
 
 
 </UL>
@@ -47,10 +44,9 @@ MeiLiftProps
 MeiLiftProps
 {
     velFieldName "U";
-    densityFieldName "rho";
     useSecondOrderTerms;
-    interpolation;
-    verbose;
+    interpolation true;
+    verbose true;
 } 
 </PRE>
 <P><B>Description:</B>

doc/forceModel_MeiLift.txt

@@ -18,17 +18,15 @@ forceModels
 MeiLiftProps
 \{
     velFieldName "U";
-    densityFieldName "density";
     useSecondOrderTerms;
-    interpolation;
-    verbose;
+    interpolation "switch1";
+    verbose "switch2";
 \}; :pre
 
 {U} = name of the finite volume fluid velocity field :ulb,l
-{density} = name of the finite volume fluid density field :l
 {useSecondOrderTerms} = switch to activate second order terms in the lift force model :l
-{interpolation} = switch to activate tri-linear interpolation of the flow quantities at the particle position :l
-{verbose} = switch to activate the report of per-particle quantities to the screen :l
+{switch1} = switch to activate tri-linear interpolation of the flow quantities at the particle position :l
+{switch2} = switch to activate the report of per-particle quantities to the screen :l
 :ule
 
 [Examples:]
@@ -40,10 +38,9 @@ forceModels
 MeiLiftProps
 \{
     velFieldName "U";
-    densityFieldName "rho";
     useSecondOrderTerms;
-    interpolation;
-    verbose;
+    interpolation true;
+    verbose true;
 \} :pre
 
 [Description:]

doc/forceModel_SchillerNaumannDrag.html

@@ -20,13 +20,10 @@
 SchillerNaumannDragProps
 {
     velFieldName "U";
-    densityFieldName "density";
 }; 
 </PRE>
 <UL><LI><I>U</I> = name of the finite volume fluid velocity field 
 
-<LI><I>density</I> = name of the finite volume gravity field 
-
 
 </UL>
 <P><B>Examples:</B>
@@ -38,7 +35,6 @@ SchillerNaumannDragProps
 SchillerNaumannDragProps
 {
     velFieldName "U";
-    densityFieldName "rho";
 } 
 </PRE>
 <P><B>Description:</B>

doc/forceModel_SchillerNaumannDrag.txt

@@ -18,11 +18,9 @@ forceModels
 SchillerNaumannDragProps
 \{
     velFieldName "U";
-    densityFieldName "density";
 \}; :pre
 
 {U} = name of the finite volume fluid velocity field :ulb,l
-{density} = name of the finite volume gravity field :l
 :ule
 
 [Examples:]
@@ -34,7 +32,6 @@ forceModels
 SchillerNaumannDragProps
 \{
     velFieldName "U";
-    densityFieldName "rho";
 \} :pre
 
 [Description:]

doc/forceModel_ShirgaonkarIB.html

@@ -20,14 +20,11 @@
 ShirgaonkarIBProps
 {
     velFieldName "U";
-    densityFieldName "density";
     pressureFieldName "pressure";
 }; 
 </PRE>
 <UL><LI><I>U</I> = name of the finite volume fluid velocity field 
 
-<LI><I>density</I> = name of the finite volume density field 
-
 <LI><I>pressure</I> = name of the finite volume pressure field 
 
 
@@ -41,7 +38,6 @@ ShirgaonkarIBProps
 ShirgaonkarIBProps
 {
     velFieldName "U";
-    densityFieldName "rho";
     pressureFieldName "p";
 } 
 </PRE>

doc/forceModel_ShirgaonkarIB.txt

@@ -18,12 +18,10 @@ forceModels
 ShirgaonkarIBProps
 \{
     velFieldName "U";
-    densityFieldName "density";
     pressureFieldName "pressure";
 \}; :pre
 
 {U} = name of the finite volume fluid velocity field :ulb,l
-{density} = name of the finite volume density field :l
 {pressure} = name of the finite volume pressure field :l
 :ule
 
@@ -36,7 +34,6 @@ forceModels
 ShirgaonkarIBProps
 \{
     velFieldName "U";
-    densityFieldName "rho";
     pressureFieldName "p";
 \} :pre
 

doc/forceModel_fieldStore.html

@@ -0,0 +1,68 @@
+<HTML>
+<CENTER><A HREF = "http://www.cfdem.com">CFDEMproject WWW Site</A> - <A HREF = "CFDEMcoupling_Manual.html#comm">CFDEM Commands</A> 
+</CENTER>
+
+
+
+
+<HR>
+
+<H3>forceModel_fieldStore command 
+</H3>
+<P><B>Syntax:</B>
+</P>
+<P>Defined in couplingProperties dictionary.
+</P>
+<PRE>forceModels
+(
+    fieldStore
+);
+fieldStoreProps
+{
+    scalarFieldNames
+    (
+        "scalarField"
+    );
+    vectorFieldNames
+    (
+        "vectorField"
+    );
+}; 
+</PRE>
+<UL><LI><I>scalarField</I> = names of the finite volume scalar fields to be stored 
+
+<LI><I>vectorField</I> = names of the finite volume vector fields to be stored 
+
+
+</UL>
+<P><B>Examples:</B>
+</P>
+<PRE>forceModels
+(
+    fieldStore
+);
+fieldStoreProps
+{
+    scalarFieldNames
+    (
+        "voidfraction"
+    );
+    vectorFieldNames
+    (
+        "U"
+    );
+} 
+</PRE>
+<P><B>Description:</B>
+</P>
+<P>This "forceModel" does not influence the particles or the flow - it is a tool to store a scalar/vector field! This is especially useful if you use a boundary condition which cannot interpreted correctly in your postporcessor (e.g. paraview).
+</P>
+<P><B>Restrictions:</B>
+</P>
+<P>none.
+</P>
+<P><B>Related commands:</B>
+</P>
+<P><A HREF = "forceModel.html">forceModel</A>
+</P>
+</HTML>

doc/forceModel_fieldStore.txt

@@ -0,0 +1,63 @@
+"CFDEMproject WWW Site"_lws - "CFDEM Commands"_lc :c
+
+:link(lws,http://www.cfdem.com)
+:link(lc,CFDEMcoupling_Manual.html#comm)
+
+:line
+
+forceModel_fieldStore command :h3
+
+[Syntax:]
+
+Defined in couplingProperties dictionary.
+
+forceModels
+(
+    fieldStore
+);
+fieldStoreProps
+\{
+    scalarFieldNames
+    (
+        "scalarField"
+    );
+    vectorFieldNames
+    (
+        "vectorField"
+    );
+\}; :pre
+
+{scalarField} = names of the finite volume scalar fields to be stored :ulb,l
+{vectorField} = names of the finite volume vector fields to be stored :l
+:ule
+
+[Examples:]
+
+forceModels
+(
+    fieldStore
+);
+fieldStoreProps
+\{
+    scalarFieldNames
+    (
+        "voidfraction"
+    );
+    vectorFieldNames
+    (
+        "U"
+    );
+\} :pre
+
+[Description:]
+
+This "forceModel" does not influence the particles or the flow - it is a tool to store a scalar/vector field! This is especially useful if you use a boundary condition which cannot interpreted correctly in your postporcessor (e.g. paraview).
+
+[Restrictions:]
+
+none.
+
+[Related commands:]
+
+"forceModel"_forceModel.html
+

doc/forceModel_gradPForce.html

@@ -20,18 +20,15 @@
 gradPForceProps
 {
     pFieldName "pressure";
-    densityFieldName "density";
     velocityFieldName "U";
-    interpolation;
+    interpolation switch1;
 }; 
 </PRE>
 <UL><LI><I>pressure</I> = name of the finite volume fluid pressure field 
 
-<LI><I>density</I> = name of the finite volume gravity field 
-
 <LI><I>U</I> = name of the finite volume fluid velocity field 
 
-<LI><I>interpolation</I> = flag to use interpolated pressure values (normally off) 
+<LI><I>switch1</I> = flag to use interpolated pressure values (normally off) 
 
 
 </UL>
@@ -44,9 +41,8 @@ gradPForceProps
 gradPForceProps
 {
     pFieldName "p";
-    densityFieldName "rho";
     velocityFieldName "U";
-    interpolation;
+    interpolation true;
 } 
 </PRE>
 <P><B>Description:</B>

doc/forceModel_gradPForce.txt

@@ -18,15 +18,13 @@ forceModels
 gradPForceProps
 \{
     pFieldName "pressure";
-    densityFieldName "density";
     velocityFieldName "U";
-    interpolation;
+    interpolation switch1;
 \}; :pre
 
 {pressure} = name of the finite volume fluid pressure field :ulb,l
-{density} = name of the finite volume gravity field :l
 {U} = name of the finite volume fluid velocity field :l
-{interpolation} = flag to use interpolated pressure values (normally off) :l
+{switch1} = flag to use interpolated pressure values (normally off) :l
 :ule
 
 [Examples:]
@@ -38,9 +36,8 @@ forceModels
 gradPForceProps
 \{
     pFieldName "p";
-    densityFieldName "rho";
     velocityFieldName "U";
-    interpolation;
+    interpolation true;
 \} :pre
 
 [Description:]

doc/forceModel_virtualMassForce.html

@@ -20,13 +20,10 @@
 virtualMassForceProps
 {
     velFieldName "U";
-    densityFieldName "density";
 }; 
 </PRE>
 <UL><LI><I>U</I> = name of the finite volume fluid velocity field 
 
-<LI><I>density</I> = name of the finite volume fluid density field 
-
 
 </UL>
 <P><B>Examples:</B>
@@ -38,7 +35,6 @@ virtualMassForceProps
 virtualMassForceProps
 {
     velFieldName "U";
-    densityFieldName "rho";
 } 
 </PRE>
 <P><B>Description:</B>

doc/forceModel_virtualMassForce.txt

@@ -18,11 +18,9 @@ forceModels
 virtualMassForceProps
 \{
     velFieldName "U";
-    densityFieldName "density";
 \}; :pre
 
 {U} = name of the finite volume fluid velocity field :ulb,l
-{density} = name of the finite volume fluid density field :l
 :ule
 
 [Examples:]
@@ -34,7 +32,6 @@ forceModels
 virtualMassForceProps
 \{
     velFieldName "U";
-    densityFieldName "rho";
 \} :pre
 
 [Description:]

doc/forceModel_viscForce.html

@@ -20,15 +20,12 @@
 viscForceProps
 {
     velocityFieldName "U";
-    densityFieldName "density";
-    interpolation;
+    interpolation "switch";
 }; 
 </PRE>
 <UL><LI><I>U</I> = name of the finite volume fluid velocity field 
 
-<LI><I>density</I> = name of the finite volume gravity field 
-
-<LI><I>interpolation</I> = flag to use interpolated stress values (normally off) 
+<LI><I>switch</I> = flag to use interpolated stress values (normally off) 
 
 
 </UL>
@@ -41,7 +38,6 @@ viscForceProps
 viscForceProps
 {
     velocityFieldName "U";
-    densityFieldName "density";
 } 
 </PRE>
 <P><B>Description:</B>

doc/forceModel_viscForce.txt

@@ -18,13 +18,11 @@ forceModels
 viscForceProps
 \{
     velocityFieldName "U";
-    densityFieldName "density";
-    interpolation;
+    interpolation "switch";
 \}; :pre
 
 {U} = name of the finite volume fluid velocity field :ulb,l
-{density} = name of the finite volume gravity field :l
-{interpolation} = flag to use interpolated stress values (normally off) :l
+{switch} = flag to use interpolated stress values (normally off) :l
 :ule
 
 [Examples:]
@@ -36,7 +34,6 @@ forceModels
 viscForceProps
 \{
     velocityFieldName "U";
-    densityFieldName "density";
 \} :pre
 
 [Description:]

doc/forceSubModel.html

@@ -0,0 +1,49 @@
+<HTML>
+<CENTER><A HREF = "http://www.cfdem.com">CFDEMproject WWW Site</A> - <A HREF = "CFDEMcoupling_Manual.html#comm">CFDEM Commands</A> 
+</CENTER>
+
+
+
+
+<HR>
+
+<H3>forceSubModel command 
+</H3>
+<P><B>Syntax:</B>
+</P>
+<P>Defined in couplingProperties sub-dictionary of the force model in use. If no force sub-model is applied ImEx is used as default. If the keyword "forceSubModels" is provided, a choice of sub model is demanded.
+</P>
+<PRE>forceSubModels
+(
+    model_x
+    model_y
+); 
+</PRE>
+<UL><LI>model = name of force sub-model to be applied 
+</UL>
+<P><B>Examples:</B>
+</P>
+<PRE>forceSubModels
+(
+    ImEx
+); 
+</PRE>
+<P>Note: This examples list might not be complete - please look for other models (forceSubModel_XY) in this documentation.
+</P>
+<P><B>Description:</B>
+</P>
+<P>The force sub model is designed to hold the settings a force model can have. For now it handles the treatExplicit, treatDEM and implDEM option. 
+</P>
+<P><B>Restrictions:</B>
+</P>
+<P>None.
+</P>
+<P><B>Related commands:</B>
+</P>
+<P><A HREF = "forceSubModel_ImEx.html">ImEx</A>
+</P>
+<P>Note: This examples list may be incomplete - please look for other models (forceSubModel_XY) in this documentation.
+</P>
+<P><B>Default:</B> none.
+</P>
+</HTML>

doc/forceSubModel.txt

@@ -0,0 +1,45 @@
+"CFDEMproject WWW Site"_lws - "CFDEM Commands"_lc :c
+
+:link(lws,http://www.cfdem.com)
+:link(lc,CFDEMcoupling_Manual.html#comm)
+
+:line
+
+forceSubModel command :h3
+
+[Syntax:]
+
+Defined in couplingProperties sub-dictionary of the force model in use. If no force sub-model is applied ImEx is used as default. If the keyword "forceSubModels" is provided, a choice of sub model is demanded.
+
+forceSubModels
+(
+    model_x
+    model_y
+); :pre
+
+model = name of force sub-model to be applied :ul
+
+[Examples:]
+
+forceSubModels
+(
+    ImEx
+); :pre
+
+Note: This examples list might not be complete - please look for other models (forceSubModel_XY) in this documentation.
+
+[Description:]
+
+The force sub model is designed to hold the settings a force model can have. For now it handles the treatExplicit, treatDEM and implDEM option. 
+
+[Restrictions:]
+
+None.
+
+[Related commands:]
+
+"ImEx"_forceSubModel_ImEx.html
+
+Note: This examples list may be incomplete - please look for other models (forceSubModel_XY) in this documentation.
+
+[Default:] none.

doc/forceSubModel_ImEx.html

@@ -0,0 +1,45 @@
+<HTML>
+<CENTER><A HREF = "http://www.cfdem.com">CFDEMproject WWW Site</A> - <A HREF = "CFDEMcoupling_Manual.html#comm">CFDEM Commands</A> 
+</CENTER>
+
+
+
+
+<HR>
+
+<H3>forceSubModel_ImEx command 
+</H3>
+<P><B>Syntax:</B>
+</P>
+<P>Defined in couplingProperties sub-dictionary of the force model in use.
+</P>
+<P>forceSubModels
+(
+    ImEx;
+);
+</P>
+<P>treatExplicit true;  // optional for some force models.
+treatDEM true;       // optional for some force models.
+implDEM true;        // optional for some force models.
+</P>
+<P><B>Examples:</B>
+</P>
+<P>forceSubModels
+(
+    ImEx;
+);
+treatExplicit true;  // optional for some force models.
+</P>
+<P><B>Description:</B>
+</P>
+<P> If no force sub-model is applied ImEx is used as default. If the keyword "forceSubModels" is provided, a choice of sub model is demanded. Depending on the force model different keywords are read and can therefrore be set (see the log file). If the keyword is provided, its value is used.
+</P>
+<P><B>Restrictions:</B>
+</P>
+<P>none.
+</P>
+<P><B>Related commands:</B>
+</P>
+<P><A HREF = "forceSubModel.html">forceSubModel</A>
+</P>
+</HTML>

doc/forceSubModel_ImEx.txt

@@ -0,0 +1,42 @@
+"CFDEMproject WWW Site"_lws - "CFDEM Commands"_lc :c
+
+:link(lws,http://www.cfdem.com)
+:link(lc,CFDEMcoupling_Manual.html#comm)
+
+:line
+
+forceSubModel_ImEx command :h3
+
+[Syntax:]
+
+Defined in couplingProperties sub-dictionary of the force model in use.
+
+forceSubModels
+(
+    ImEx;
+);
+
+treatExplicit true;  // optional for some force models.
+treatDEM true;       // optional for some force models.
+implDEM true;        // optional for some force models.
+
+[Examples:]
+
+forceSubModels
+(
+    ImEx;
+);
+treatExplicit true;  // optional for some force models.
+
+[Description:]
+
+ If no force sub-model is applied ImEx is used as default. If the keyword "forceSubModels" is provided, a choice of sub model is demanded. Depending on the force model different keywords are read and can therefrore be set (see the log file). If the keyword is provided, its value is used.
+
+[Restrictions:]
+
+none.
+
+[Related commands:]
+
+"forceSubModel"_forceSubModel.html
+

doc/forceSubModel_ImExCorr.html

@@ -0,0 +1,46 @@
+<HTML>
+<CENTER><A HREF = "http://www.cfdem.com">CFDEMproject WWW Site</A> - <A HREF = "CFDEMcoupling_Manual.html#comm">CFDEM Commands</A> 
+</CENTER>
+
+
+
+
+<HR>
+
+<H3>forceSubModel_ImExCorr command 
+</H3>
+<P><B>Syntax:</B>
+</P>
+<P>Defined in couplingProperties sub-dictionary of the force model in use.
+</P>
+<P>forceSubModels
+(
+    ImExCorr;
+);
+</P>
+<P>treatExplicit true;  // optional for some force models.
+treatDEM true;       // optional for some force models.
+implDEM true;        // optional for some force models.
+explicitInterpCorr true; // optional for some force models.
+</P>
+<P><B>Examples:</B>
+</P>
+<P>forceSubModels
+(
+    ImExCorr;
+);
+treatExplicit true;  // optional for some force models.
+</P>
+<P><B>Description:</B>
+</P>
+<P> Same as ImEx, but it additionally reads "explicitInterpCorr" to correct the error steming from interpolation of Ufluid and averaging of Uparticles.
+</P>
+<P><B>Restrictions:</B>
+</P>
+<P>none.
+</P>
+<P><B>Related commands:</B>
+</P>
+<P><A HREF = "forceSubModel.html">forceSubModel</A>
+</P>
+</HTML>

doc/forceSubModel_ImExCorr.txt

@@ -0,0 +1,43 @@
+"CFDEMproject WWW Site"_lws - "CFDEM Commands"_lc :c
+
+:link(lws,http://www.cfdem.com)
+:link(lc,CFDEMcoupling_Manual.html#comm)
+
+:line
+
+forceSubModel_ImExCorr command :h3
+
+[Syntax:]
+
+Defined in couplingProperties sub-dictionary of the force model in use.
+
+forceSubModels
+(
+    ImExCorr;
+);
+
+treatExplicit true;  // optional for some force models.
+treatDEM true;       // optional for some force models.
+implDEM true;        // optional for some force models.
+explicitInterpCorr true; // optional for some force models.
+
+[Examples:]
+
+forceSubModels
+(
+    ImExCorr;
+);
+treatExplicit true;  // optional for some force models.
+
+[Description:]
+
+ Same as ImEx, but it additionally reads "explicitInterpCorr" to correct the error steming from interpolation of Ufluid and averaging of Uparticles.
+
+[Restrictions:]
+
+none.
+
+[Related commands:]
+
+"forceSubModel"_forceSubModel.html
+

doc/liggghtsCommandModel_writeLiggghts.html

@@ -26,13 +26,12 @@ writeLiggghtsProps
     verbose;
 } 
 </PRE>
-<UL><LI><I>switch1</I> = switch (choose on/off) to select if only last step is stored or every write step. 
+<UL><LI><I>switch1</I> = switch (choose on/off) to select if only last step is stored or every write step (default on). 
 
-<LI><I>name</I> = name of the restart file to be written in /$caseDir/DEM/  default default "liggghts.restartCFDEM" 
+<LI><I>name</I> = name of the restart file to be written in /$caseDir/DEM/  default (default "liggghts.restartCFDEM") 
 
-<LI><I>switch2</I> = switch (choose on/off) to select if only one restart file $name or many files $name_$timeStamp are written 
-
-<LI><I>verbose</I> = (normally off) for verbose run 
+<LI><I>switch2</I> = switch (choose on/off) to select if only one restart file $name or many files $name_$timeStamp are written (default off):l
+<I>verbose</I> = (default off) for verbose run 
 
 
 </UL>

doc/liggghtsCommandModel_writeLiggghts.txt

@@ -24,10 +24,10 @@ writeLiggghtsProps
     verbose;
 \} :pre
 
-{switch1} = switch (choose on/off) to select if only last step is stored or every write step. :ulb,l
-{name} = name of the restart file to be written in /$caseDir/DEM/  default default "liggghts.restartCFDEM" :l
-{switch2} = switch (choose on/off) to select if only one restart file $name or many files $name_$timeStamp are written :l
-{verbose} = (normally off) for verbose run :l
+{switch1} = switch (choose on/off) to select if only last step is stored or every write step (default on). :ulb,l
+{name} = name of the restart file to be written in /$caseDir/DEM/  default (default "liggghts.restartCFDEM") :l
+{switch2} = switch (choose on/off) to select if only one restart file $name or many files $name_$timeStamp are written (default off):l
+{verbose} = (default off) for verbose run :l
 :ule
 
 [Examples:]

doc/locateModel_turboEngineSearchM2M.html

@@ -1,58 +0,0 @@
-<HTML>
-<CENTER><A HREF = "http://www.cfdem.com">CFDEMproject WWW Site</A> - <A HREF = "CFDEMcoupling_Manual.html#comm">CFDEM Commands</A> 
-</CENTER>
-
-
-
-
-<HR>
-
-<H3>locateModel_turboEngineM2MSearch command 
-</H3>
-<P><B>Syntax:</B>
-</P>
-<P>Defined in couplingProperties dictionary.
-</P>
-<PRE>locateModel turboEngineM2M;
-turboEngineM2MProps
-{
-    turboEngineProps
-    {
-        treeSearch switch1;
-    }
-} 
-</PRE>
-<UL><LI><I>switch1</I> = switch to use tree search algorithm 
-
-
-</UL>
-<P><B>Examples:</B>
-</P>
-<PRE>locateModel turboEngineM2M;
-turboEngineM2MProps
-{
-    turboEngineProps
-    {
-        treeSearch true;
-    }
-} 
-</PRE>
-<P><B>Description:</B>
-</P>
-<P>The locateModel "turboEngineM2M" locates the CFD cell and cellID corresponding to a given position. The algorithm is improved compared to engine search to show better parallel performance.
-</P>
-<P>The turboEngineM2MSearch locate Model can be used with different settings to use different algorithms:
-</P>
-<UL><LI>faceDecomp false; treeSearch false;     will execute some geometric (linear) search using the last known cellID (recommended) 
-
-<LI>faceDecomp false; treeSearch true;      will use a recursive tree structure to find the cell. 
-
-
-</UL>
-<P><B>Restrictions:</B> This model can only be used with many to many data exchange model!
-</P>
-<P><B>Related commands:</B>
-</P>
-<P><A HREF = "locateModel.html">locateModel</A>
-</P>
-</HTML>

doc/locateModel_turboEngineSearchM2M.txt

@@ -1,52 +0,0 @@
-"CFDEMproject WWW Site"_lws - "CFDEM Commands"_lc :c
-
-:link(lws,http://www.cfdem.com)
-:link(lc,CFDEMcoupling_Manual.html#comm)
-
-:line
-
-locateModel_turboEngineM2MSearch command :h3
-
-[Syntax:]
-
-Defined in couplingProperties dictionary.
-
-locateModel turboEngineM2M;
-turboEngineM2MProps
-\{
-    turboEngineProps
-    \{
-        treeSearch switch1;
-    \}
-\} :pre
-
-{switch1} = switch to use tree search algorithm :ulb,l
-:ule
-
-[Examples:]
-
-locateModel turboEngineM2M;
-turboEngineM2MProps
-\{
-    turboEngineProps
-    \{
-        treeSearch true;
-    \}
-\} :pre
-
-[Description:]
-
-The locateModel "turboEngineM2M" locates the CFD cell and cellID corresponding to a given position. The algorithm is improved compared to engine search to show better parallel performance.
-
-The turboEngineM2MSearch locate Model can be used with different settings to use different algorithms:
-
-faceDecomp false; treeSearch false;     will execute some geometric (linear) search using the last known cellID (recommended) :ulb,l
-faceDecomp false; treeSearch true;      will use a recursive tree structure to find the cell. :l
-:ule
-
-[Restrictions:] This model can only be used with many to many data exchange model!
-
-[Related commands:]
-
-"locateModel"_locateModel.html
-

doc/momCoupleModel.html

@@ -33,7 +33,8 @@
 </P>
 <P>Note that the variable "imExSplitFactor" can be set in the couplingProperties in order to treat implicitly defined forces (in the implementation of the force model) as explicit ones. "imExSplitFactor 1.0;" is set by default, meaning that all implicit forces will be considered implicitly, whereas "imExSplitFactor 0.0;" would mean that implicitly defined forces will be treated in an explicit fashion.
 </P>
-<P><B>Description:</B>
+<P>Note that the switch "treatVoidCellsAsExplicitForce true;" can be set in the couplingProperties in order to change the treatment of cells which are void of particles. This is only relevant if (i) smoothing is used, and (ii) implicit force coupling is performed. By default, the particle veloctiy field (Us) will be smoothed to obtain a meaningful reference quantity for the implicit force coupling. In case "treatVoidCellsAsExplicitForce true;" is set, however, Us will not be smoothed and implicit forces (after the smoothing has been performed) in cells void of particles be treated as explicit ones. This avoids the problem of defining Us in cells that are void of particles, but for which an implicit coupling force is obtained in the smoothing process.
+<B>Description:</B>
 </P>
 <P>The momCoupleModel is the base class for momentum exchange between DEM and CFD simulation.
 </P>

doc/momCoupleModel.txt

@@ -31,6 +31,7 @@ Forces can be coupled in an implicit way to the fluid solver (i.e., when solving
 
 Note that the variable "imExSplitFactor" can be set in the couplingProperties in order to treat implicitly defined forces (in the implementation of the force model) as explicit ones. "imExSplitFactor 1.0;" is set by default, meaning that all implicit forces will be considered implicitly, whereas "imExSplitFactor 0.0;" would mean that implicitly defined forces will be treated in an explicit fashion.
 
+Note that the switch "treatVoidCellsAsExplicitForce true;" can be set in the couplingProperties in order to change the treatment of cells which are void of particles. This is only relevant if (i) smoothing is used, and (ii) implicit force coupling is performed. By default, the particle veloctiy field (Us) will be smoothed to obtain a meaningful reference quantity for the implicit force coupling. In case "treatVoidCellsAsExplicitForce true;" is set, however, Us will not be smoothed and implicit forces (after the smoothing has been performed) in cells void of particles be treated as explicit ones. This avoids the problem of defining Us in cells that are void of particles, but for which an implicit coupling force is obtained in the smoothing process.
 [Description:]
 
 The momCoupleModel is the base class for momentum exchange between DEM and CFD simulation.

src/.gitignore

@@ -0,0 +1,7 @@
+*.o
+*.d
+*.a
+*.dep
+log_*
+log.*
+*~

src/finiteVolume/Make/files

@@ -1,5 +1,8 @@
 fvPatchFields = fields/fvPatchFields
 derivedFvPatchFields = $(fvPatchFields)/derived
+
+general = cfdTools/general
+
 $(derivedFvPatchFields)/uniformFixedValueVoidfraction/uniformFixedValueVoidfractionFvPatchFields.C
 $(derivedFvPatchFields)/uniformFixedValueTube/uniformFixedValueTubeFvPatchFields.C
 

src/finiteVolume/Make/options

@@ -2,6 +2,9 @@ EXE_INC = \
     -I$(LIB_SRC)/triSurface/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -Wno-old-style-cast \
+    -Wno-unused-result \
+    -Wno-literal-suffix
 
 LIB_LIBS = \
     -lOpenFOAM \

src/finiteVolume/fields/fvPatchFields/derived/uniformFixedValueTube/uniformFixedValueTubeFvPatchField.C

@@ -88,7 +88,7 @@ uniformFixedValueTubeFvPatchField<Type>::uniformFixedValueTubeFvPatchField
 )
 :
     fixedValueFvPatchField<Type>(p, iF),
-    uniformValue_(DataEntry<Type>::New("uniformValue", dict)),
+    uniformValue_(Function1<Type>::New("uniformValue", dict)),
 	pName_("p"), //JOKER
 	phiName_("phi"), //JOKER
     velocityFieldName_("U"),

src/finiteVolume/fields/fvPatchFields/derived/uniformFixedValueTube/uniformFixedValueTubeFvPatchField.H

@@ -40,7 +40,7 @@ SourceFiles
 
 #include "Random.H"
 #include "fixedValueFvPatchFields.H"
-#include "DataEntry.H"
+#include "Function1.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -58,7 +58,7 @@ class uniformFixedValueTubeFvPatchField
 {
     // Private data
 
-        autoPtr<DataEntry<Type> > uniformValue_;
+        autoPtr<Function1<Type> > uniformValue_;
 
 		word pName_; //JOKER pressure
 		

src/finiteVolume/fields/fvPatchFields/derived/uniformFixedValueVoidfraction/uniformFixedValueVoidfractionFvPatchField.C

@@ -41,7 +41,6 @@ uniformFixedValueVoidfractionFvPatchField<Type>::uniformFixedValueVoidfractionFv
 :
     fixedValueFvPatchField<Type>(p, iF),
     uniformValue_(),
-    //voidfractionFieldName_(dict.lookupOrDefault<word>("voidfraction", "voidfraction"))
     voidfractionFieldName_("voidfraction")
 {}
 
@@ -56,7 +55,7 @@ uniformFixedValueVoidfractionFvPatchField<Type>::uniformFixedValueVoidfractionFv
 )
 :
     fixedValueFvPatchField<Type>(p, iF),
-    uniformValue_(ptf.uniformValue_().clone().ptr()),
+    uniformValue_(ptf.uniformValue_, false),
     voidfractionFieldName_("voidfraction")
 {
     const scalar t = this->db().time().timeOutputValue();
@@ -73,7 +72,7 @@ uniformFixedValueVoidfractionFvPatchField<Type>::uniformFixedValueVoidfractionFv
 )
 :
     fixedValueFvPatchField<Type>(p, iF),
-    uniformValue_(DataEntry<Type>::New("uniformValue", dict)),
+    uniformValue_(Function1<Type>::New("uniformValue", dict)),
     voidfractionFieldName_("voidfraction")
 {
     const scalar t = this->db().time().timeOutputValue();
@@ -88,7 +87,7 @@ uniformFixedValueVoidfractionFvPatchField<Type>::uniformFixedValueVoidfractionFv
 )
 :
     fixedValueFvPatchField<Type>(ptf),
-    uniformValue_(ptf.uniformValue_().clone().ptr()),
+    uniformValue_(ptf.uniformValue_, false),
     voidfractionFieldName_("voidfraction")
 {
     const scalar t = this->db().time().timeOutputValue();
@@ -104,7 +103,7 @@ uniformFixedValueVoidfractionFvPatchField<Type>::uniformFixedValueVoidfractionFv
 )
 :
     fixedValueFvPatchField<Type>(ptf, iF),
-    uniformValue_(ptf.uniformValue_().clone().ptr()),
+    uniformValue_(ptf.uniformValue_, false),
     voidfractionFieldName_("voidfraction")
 {
     const scalar t = this->db().time().timeOutputValue();

src/finiteVolume/fields/fvPatchFields/derived/uniformFixedValueVoidfraction/uniformFixedValueVoidfractionFvPatchField.H

@@ -40,7 +40,7 @@ SourceFiles
 
 #include "Random.H"
 #include "fixedValueFvPatchFields.H"
-#include "DataEntry.H"
+#include "Function1.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -58,7 +58,7 @@ class uniformFixedValueVoidfractionFvPatchField
 {
     // Private data
 
-        autoPtr<DataEntry<Type> > uniformValue_;
+        autoPtr<Function1<Type> > uniformValue_;
 
         word voidfractionFieldName_;
 

src/lagrangian/cfdemParticle/.gitignore

@@ -0,0 +1 @@
+lnInclude

src/lagrangian/cfdemParticle/Make/files

@@ -1,81 +1,105 @@
 cfdemCloud = cfdemCloud
+cfdTools = cfdTools
+energyModels = subModels/energyModel
 forceModels = subModels/forceModel
+forceSubModels = subModels/forceModel/forceSubModels
 forceModelsMS = subModels/forceModelMS
+thermCondModels = subModels/thermCondModel
+chemistryModels = subModels/chemistryModel
 IOModels = subModels/IOModel
 voidFractionModels = subModels/voidFractionModel
 locateModels = subModels/locateModel
 meshMotionModels = subModels/meshMotionModel
 momCoupleModels = subModels/momCoupleModel
+regionModels = subModels/regionModel
 dataExchangeModels = subModels/dataExchangeModel
 averagingModels = subModels/averagingModel
 clockModels = subModels/clockModel
 liggghtsCommandModels = subModels/liggghtsCommandModel
 smoothingModels = subModels/smoothingModel
 probeModels = subModels/probeModel
+otherForceModels = subModels/otherForceModel
 
 $(cfdemCloud)/cfdemCloud.C
-derived/cfdemCloudBiDisperse/cfdemCloudBiDisperse.C
+$(cfdemCloud)/cfdemCloudIO.C
 derived/cfdemCloudIB/cfdemCloudIB.C
 derived/cfdemCloudMS/cfdemCloudMS.C
+derived/cfdemCloudEnergy/cfdemCloudEnergy.C
+
+$(cfdTools)/global.C
+$(cfdTools)/newGlobal.C
+
+$(energyModels)/energyModel/energyModel.C
+$(energyModels)/energyModel/newEnergyModel.C
+$(energyModels)/heatTransferGunn/heatTransferGunn.C
+$(energyModels)/heatTransferGunnImplicit/heatTransferGunnImplicit.C
+$(energyModels)/reactionHeat/reactionHeat.C
+
+$(thermCondModels)/thermCondModel/thermCondModel.C
+$(thermCondModels)/thermCondModel/newThermCondModel.C
+$(thermCondModels)/SyamlalThermCond/SyamlalThermCond.C
+$(thermCondModels)/noTherm/noThermCond.C
 
 $(forceModels)/forceModel/forceModel.C
 $(forceModels)/forceModel/newForceModel.C
 $(forceModels)/noDrag/noDrag.C
 $(forceModels)/checkCouplingInterval/checkCouplingInterval.C
 $(forceModels)/DiFeliceDrag/DiFeliceDrag.C
-$(forceModels)/DiFeliceDragNLift/DiFeliceDragNLift.C
+$(forceModels)/fieldStore/fieldStore.C
 $(forceModels)/GidaspowDrag/GidaspowDrag.C
 $(forceModels)/SchillerNaumannDrag/SchillerNaumannDrag.C
 $(forceModels)/Archimedes/Archimedes.C
 $(forceModels)/ArchimedesIB/ArchimedesIB.C
 $(forceModels)/interface/interface.C
 $(forceModels)/ShirgaonkarIB/ShirgaonkarIB.C
-$(forceModels)/interfaceParticleProbe/interfaceParticleProbe.C
-$(forceModels)/fieldStore/fieldStore.C
-$(forceModels)/fieldTimeAverage/fieldTimeAverage.C
-$(forceModels)/fieldBound/fieldBound.C
-$(forceModels)/volWeightedAverage/volWeightedAverage.C
-$(forceModels)/totalMomentumExchange/totalMomentumExchange.C
 $(forceModels)/KochHillDrag/KochHillDrag.C
 $(forceModels)/KochHillRWDrag/KochHillRWDrag.C
-$(forceModels)/BeetstraDrag/multiphaseFlowBasic/multiphaseFlowBasic.C
-$(forceModels)/BeetstraDrag/BeetstraDrag.C
-$(forceModels)/LaEuScalarLiquid/LaEuScalarLiquid.C
 $(forceModels)/LaEuScalarTemp/LaEuScalarTemp.C
-$(forceModels)/LaEuScalarDust/LaEuScalarDust.C
 $(forceModels)/virtualMassForce/virtualMassForce.C
 $(forceModels)/gradPForce/gradPForce.C
-$(forceModels)/gradULiftForce/gradULiftForce.C
-$(forceModels)/HollowayDrag/HollowayDrag.C
 $(forceModels)/viscForce/viscForce.C
 $(forceModels)/MeiLift/MeiLift.C
-$(forceModels)/melting/melting.C
-$(forceModels)/KochHillDragNLift/KochHillDragNLift.C
-$(forceModels)/stokesSpheroidDrag/stokesSpheroidDrag.C
-$(forceModels)/solidsPressureForce/solidsPressureForce.C
-$(forceModels)/periodicPressure/periodicPressure.C
-$(forceModels)/periodicPressureControl/periodicPressureControl.C
-$(forceModels)/averageSlipVel/averageSlipVel.C
 $(forceModels)/particleCellVolume/particleCellVolume.C
+$(forceModels)/fieldTimeAverage/fieldTimeAverage.C
+$(forceModels)/volWeightedAverage/volWeightedAverage.C
+$(forceModels)/BeetstraDrag/BeetstraDrag.C
+$(forceModels)/dSauter/dSauter.C
+$(forceModels)/Fines/Fines.C
+$(forceModels)/Fines/FinesFields.C
+$(forceModels)/Fines/FanningDynFines.C
+$(forceModels)/Fines/ErgunStatFines.C
+$(forceModels)/granKineticEnergy/granKineticEnergy.C
 
 $(forceModelsMS)/forceModelMS/forceModelMS.C
 $(forceModelsMS)/forceModelMS/newForceModelMS.C
 $(forceModelsMS)/DiFeliceDragMS/DiFeliceDragMS.C
-$(forceModelsMS)/GidaspowDragMS/GidaspowDragMS.C
-$(forceModelsMS)/noDragMS/noDragMS.C
+
+$(forceSubModels)/forceSubModel/newForceSubModel.C
+$(forceSubModels)/forceSubModel/forceSubModel.C
+$(forceSubModels)/ImEx/ImEx.C
+$(forceSubModels)/ScaleForce/ScaleForce.C
+$(forceSubModels)/scaleForceBoundary/scaleForceBoundary.C
+
+$(otherForceModels)/otherForceModel/otherForceModel.C
+$(otherForceModels)/otherForceModel/newOtherForceModel.C
+$(otherForceModels)/gravity/gravity.C
+$(otherForceModels)/weightSecondaryPhase/weightSecondaryPhase.C
+$(otherForceModels)/expParticleForces/expParticleForces.C
+
+$(chemistryModels)/chemistryModel/chemistryModel.C
+$(chemistryModels)/chemistryModel/newChemistryModel.C
+$(chemistryModels)/species/species.C
+$(chemistryModels)/noChemistry/noChemistry.C
 
 $(probeModels)/probeModel/probeModel.C
 $(probeModels)/probeModel/newProbeModel.C
 $(probeModels)/noProbe/noProbe.C
 $(probeModels)/particleProbe/particleProbe.C
-/*$(probeModels)/interfaceParticleProbe/interfaceParticleProbe.C*/
 
 $(IOModels)/IOModel/IOModel.C
 $(IOModels)/IOModel/newIOModel.C
 $(IOModels)/noIO/noIO.C
 $(IOModels)/basicIO/basicIO.C
-$(IOModels)/tempIO/tempIO.C
-$(IOModels)/colorIO/colorIO.C
 $(IOModels)/trackIO/trackIO.C
 $(IOModels)/sophIO/sophIO.C
 
@@ -83,50 +107,44 @@ $(voidFractionModels)/voidFractionModel/voidFractionModel.C
 $(voidFractionModels)/voidFractionModel/newVoidFractionModel.C
 $(voidFractionModels)/centreVoidFraction/centreVoidFraction.C
 $(voidFractionModels)/dividedVoidFraction/dividedVoidFraction.C
-$(voidFractionModels)/dividedVoidFractionBiDi/dividedVoidFractionBiDi.C
 $(voidFractionModels)/dividedVoidFractionMS/dividedVoidFractionMS.C
 $(voidFractionModels)/bigParticleVoidFraction/bigParticleVoidFraction.C
 $(voidFractionModels)/GaussVoidFraction/GaussVoidFraction.C
 $(voidFractionModels)/IBVoidFraction/IBVoidFraction.C
-$(voidFractionModels)/weightedNeigbhorsVoidFraction/weightedNeigbhorsVoidFraction.C
 
 $(locateModels)/locateModel/locateModel.C
 $(locateModels)/locateModel/newLocateModel.C
 $(locateModels)/standardSearch/standardSearch.C
 $(locateModels)/engineSearch/engineSearch.C
-$(locateModels)/engineSearchMany2Many/engineSearchMany2Many.C
 $(locateModels)/turboEngineSearch/turboEngineSearch.C
-$(locateModels)/turboEngineSearchM2M/turboEngineSearchM2M.C
+$(locateModels)/engineSearchMany2Many/engineSearchMany2Many.C
 $(locateModels)/engineSearchIB/engineSearchIB.C
-$(locateModels)/hyperEngineSearch/hyperEngineSearch.C
-$(locateModels)/ijkSearch/ijkSearch.C
 
 $(meshMotionModels)/meshMotionModel/meshMotionModel.C
 $(meshMotionModels)/meshMotionModel/newMeshMotionModel.C
 $(meshMotionModels)/noMeshMotion/noMeshMotion.C
-$(meshMotionModels)/DEMdrivenMeshMotion/DEMdrivenMeshMotion.C
 
 $(momCoupleModels)/momCoupleModel/momCoupleModel.C
 $(momCoupleModels)/momCoupleModel/newMomCoupleModel.C
 $(momCoupleModels)/explicitCouple/explicitCouple.C
-$(momCoupleModels)/explicitCoupleSource/explicitCoupleSource.C
 $(momCoupleModels)/implicitCouple/implicitCouple.C
 $(momCoupleModels)/noCouple/noCouple.C
 
+$(regionModels)/regionModel/regionModel.C
+$(regionModels)/regionModel/newRegionModel.C
+$(regionModels)/allRegion/allRegion.C
+
 $(dataExchangeModels)/dataExchangeModel/dataExchangeModel.C
 $(dataExchangeModels)/dataExchangeModel/newDataExchangeModel.C
 $(dataExchangeModels)/oneWayVTK/oneWayVTK.C
 $(dataExchangeModels)/twoWayFiles/twoWayFiles.C
 $(dataExchangeModels)/noDataExchange/noDataExchange.C
 $(dataExchangeModels)/twoWayMPI/twoWayMPI.C
-$(dataExchangeModels)/twoWayM2M/twoWayM2M.C
-$(dataExchangeModels)/twoWayMany2Many/twoWayMany2Many.C
 
 $(averagingModels)/averagingModel/averagingModel.C
 $(averagingModels)/averagingModel/newAveragingModel.C
 $(averagingModels)/dilute/dilute.C
 $(averagingModels)/dense/dense.C
-$(averagingModels)/denseBiDi/denseBiDi.C
 
 $(clockModels)/clockModel/clockModel.C
 $(clockModels)/clockModel/newClockModel.C
@@ -135,7 +153,6 @@ $(clockModels)/noClock/noClock.C
 
 $(liggghtsCommandModels)/liggghtsCommandModel/liggghtsCommandModel.C
 $(liggghtsCommandModels)/liggghtsCommandModel/newLiggghtsCommandModel.C
-$(liggghtsCommandModels)/colorParticles/colorParticles.C
 $(liggghtsCommandModels)/execute/execute.C
 $(liggghtsCommandModels)/runLiggghts/runLiggghts.C
 $(liggghtsCommandModels)/writeLiggghts/writeLiggghts.C
@@ -145,6 +162,5 @@ $(smoothingModels)/smoothingModel/smoothingModel.C
 $(smoothingModels)/smoothingModel/newSmoothingModel.C
 $(smoothingModels)/noSmoothing/noSmoothing.C
 $(smoothingModels)/constDiffSmoothing/constDiffSmoothing.C
-$(smoothingModels)/localPSizeDiffSmoothing/localPSizeDiffSmoothing.C
 
 LIB = $(CFDEM_LIB_DIR)/lib$(CFDEM_LIB_NAME)

src/lagrangian/cfdemParticle/Make/options

@@ -1,38 +1,38 @@
 sinclude $(GENERAL_RULES)/mplib$(WM_MPLIB)
 sinclude $(RULES)/mplib$(WM_MPLIB)
 
+GIT_VERSION := $(shell git describe --dirty --always --tags)
+PFLAGS+= -DGITVERSION=\"$(GIT_VERSION)\"
+
+include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
 
 EXE_INC = \
      $(PFLAGS) \
      $(PINC) \
-    -I ./cfdemParticle \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/transportModels \
-    -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
     -I$(LIB_SRC)/OpenFOAM/containers/HashTables/labelHashSet \
     -I$(CFDEM_LIGGGHTS_SRC_DIR) \
     -I$(CFDEM_M2MLIB_PATH) \
-    -I$(CFDEM_Many2ManyLIB_PATH) \
     -I$(CFDEM_SRC_DIR)/cfdTools \
+    -Wno-old-style-cast \
+    -Wno-unused-result \
+    -Wno-literal-suffix
+
 
 LIB_LIBS = \
      $(PLIBS) \
-    -L$(FOAM_USER_LIBBIN) \
+    -L$(CFDEM_LIB_DIR) \
     -lfiniteVolume \
-    -lincompressibleRASModels \
-    -lincompressibleLESModels \
+    -lincompressibleTurbulenceModels \
     -lmeshTools \
     -llagrangian \
     -lmpi_cxx \
-    -L$(CFDEM_LIGGGHTS_SRC_DIR) \
-    -l$(CFDEM_LIGGGHTS_LIB_NAME) \
-    -L$(CFDEM_M2MLIB_PATH) \
-    -lcouple \
+    -Wl,-rpath,$(CFDEM_LIGGGHTS_BIN_DIR) \
+    -L$(CFDEM_LIGGGHTS_BIN_DIR) \
+    -lliggghts \
     -L$(CFDEM_Many2ManyLIB_PATH) \
-    -lcoupleMany2Many \
-
-/* add -I$(CFDEM_POEMSLIB_PATH) \ to EXE_INC */
-/*    -L$(CFDEM_POEMSLIB_PATH) \ */
-/*    -lpoems   */
+    -lcoupleMany2Many

src/lagrangian/cfdemParticle/cfdTools/checkModelType.H

@@ -4,8 +4,7 @@
     word modelType = particleCloud.modelType();
 
     //Warning << "model type not being checked" << endl;
-
-    if (modelType=="B"){
+    if (modelType=="Bfull"){
         Info << "solving volume averaged Navier Stokes equations of type B\n"<< endl;
 
         // check if Archimedes is used
@@ -18,6 +17,41 @@
         if(!found)
             FatalError <<"Archimedes model not found!\n" << abort(FatalError);
 
+        // check if gradPForce is used
+        found=false;
+        forAll(particleCloud.forceModels(),i)
+        {
+            if(particleCloud.forceModels()[i]=="gradPForce")
+                found=true;
+        }
+        if(!found)
+            FatalError <<"gradPForce model not found!\n" << abort(FatalError);
+
+        // check if viscForce is used
+        found=false;
+        forAll(particleCloud.forceModels(),i)
+        {
+            if(particleCloud.forceModels()[i]=="viscForce")
+                found=true;
+        }
+        if(!found)
+            FatalError <<"viscForce model not found!\n" << abort(FatalError);
+
+    }else if(modelType=="B"){
+        Info << "solving volume averaged Navier Stokes equations of type B\n"<< endl;
+
+        // check if Archimedes is used
+        bool found=false;
+        forAll(particleCloud.forceModels(),i)
+        {
+            if(particleCloud.forceModels()[i]=="Archimedes")
+                found=true;
+        }
+        if(!found)
+            FatalError <<"Archimedes model not found!\n" << abort(FatalError);
+
+
+
         // check if gradP and viscForce are used
         found=false;
         forAll(particleCloud.forceModels(),i)

src/lagrangian/cfdemParticle/cfdTools/compressibleContinuityErrsPU.H

@@ -0,0 +1,53 @@
+/*---------------------------------------------------------------------------*\
+    CFDEMcoupling - Open Source CFD-DEM coupling
+
+    CFDEMcoupling is part of the CFDEMproject
+    www.cfdem.com
+                                Christoph Goniva, christoph.goniva@cfdem.com
+                                Copyright (C) 1991-2009 OpenCFD Ltd.
+                                Copyright (C) 2009-2012 JKU, Linz
+                                Copyright (C) 2012-     DCS Computing GmbH,Linz
+-------------------------------------------------------------------------------
+License
+    This file is part of CFDEMcoupling.
+
+    CFDEMcoupling 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.
+
+    CFDEMcoupling 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 CFDEMcoupling.  If not, see <http://www.gnu.org/licenses/>.
+
+Global
+    continuityErrs
+
+Description
+    Calculates and prints the continuity errors.
+    The code is an evolution of compressibleContinuityErrs.H in OpenFOAM(R) 2.3.x,
+    where additional functionality for CFD-DEM coupling is added.
+\*---------------------------------------------------------------------------*/
+
+{
+    dimensionedScalar totalMass = fvc::domainIntegrate(rho*voidfraction);
+
+    scalar sumLocalContErr =
+        (fvc::domainIntegrate(mag(rho - thermo.rho())*voidfraction)/totalMass).value();
+
+    scalar globalContErr =
+        (fvc::domainIntegrate((rho - thermo.rho())*voidfraction)/totalMass).value();
+
+    cumulativeContErr += globalContErr;
+
+    Info<< "time step continuity errors : sum local = " << sumLocalContErr
+        << ", global = " << globalContErr
+        << ", cumulative = " << cumulativeContErr
+        << endl;
+}
+
+// ************************************************************************* //

src/lagrangian/cfdemParticle/cfdTools/continuityErrorPhiPU.H

@@ -34,7 +34,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 {
-    volScalarField contErr( fvc::div(phiGes) + fvc::ddt(voidfraction) );
+    volScalarField contErr( fvc::div(phi) + fvc::ddt(voidfraction) );
 
     scalar sumLocalContErr = runTime.deltaTValue()*
         mag(contErr)().weightedAverage(mesh.V()).value();

src/lagrangian/cfdemParticle/cfdTools/debugInfo.H

@@ -23,8 +23,8 @@
           meanUs_array[i] += velocities_[index][i];
       }
       meanR_array += radii_[index][0];
-      particleVolume_radius += radii_[index][0]*radii_[index][0]*radii_[index][0]*4./3.*3.1415;
-      
+      particleVolume_radius += radii_[index][0]*radii_[index][0]*radii_[index][0]*(4./3.)*M_PI;
+
       // loop subCells
       for(int subCell=0;subCell<voidFractionM().cellsPerParticle()[index][0];subCell++)
       {
@@ -63,17 +63,17 @@
   if(countCell>0)
   {
       meanAlpha_field /= countCell;
-      meanU_field /= countCell; 
+      meanU_field /= countCell;
       meanUs_field /= countCell;
   }
   else
   {
       meanAlpha_field = 0;
-      meanU_field = vector(0,0,0); 
+      meanU_field = vector(0,0,0);
       meanUs_field = vector(0,0,0);
   }
-  meanUs_array /= numberOfParticles()+SMALL;   
-  meanR_array /= numberOfParticles()+SMALL;   
+  meanUs_array /= numberOfParticles()+SMALL;
+  meanR_array /= numberOfParticles()+SMALL;
 
   Info <<"=============================================================================" << endl;
   Info << "Debug Info, only serial and not tested!" << endl;

src/lagrangian/cfdemParticle/cfdTools/global.C

@@ -0,0 +1,80 @@
+/*---------------------------------------------------------------------------*\
+    CFDEMcoupling - Open Source CFD-DEM coupling
+
+    CFDEMcoupling is part of the CFDEMproject
+    www.cfdem.com
+                                Christoph Goniva, christoph.goniva@cfdem.com
+                                Copyright 2009-2012 JKU Linz
+                                Copyright 2012-     DCS Computing GmbH, Linz
+-------------------------------------------------------------------------------
+License
+    This file is part of CFDEMcoupling.
+
+    CFDEMcoupling 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.
+
+    CFDEMcoupling 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 CFDEMcoupling; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Description
+    This code is designed to realize coupled CFD-DEM simulations using LIGGGHTS
+    and OpenFOAM(R). Note: this code is not part of OpenFOAM(R) (see DISCLAIMER).
+\*---------------------------------------------------------------------------*/
+
+#include "error.H"
+#include "global.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(global, 0);
+
+defineRunTimeSelectionTable(global, dictionary);
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+void global::info()
+{
+    Info << "\nYou are currently using:" << endl;
+    Info << "OpenFOAM version: " << FOAMversion << endl;
+    Info << "OpenFOAM build: " << FOAMbuild << endl;
+    Info << "CFDEM build: " << CFDEMversion << "\n" << endl;
+}
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from components
+global::global
+(
+    const dictionary& dict,
+    cfdemCloud& sm
+)
+:
+    dict_(dict),
+    particleCloud_(sm),
+    CFDEMversion(GITVERSION)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+global::~global()
+{}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/lagrangian/cfdemParticle/cfdTools/global.H

@@ -0,0 +1,128 @@
+/*---------------------------------------------------------------------------*\
+    CFDEMcoupling - Open Source CFD-DEM coupling
+
+    CFDEMcoupling is part of the CFDEMproject
+    www.cfdem.com
+                                Christoph Goniva, christoph.goniva@cfdem.com
+                                Copyright 2009-2012 JKU Linz
+                                Copyright 2012-     DCS Computing GmbH, Linz
+-------------------------------------------------------------------------------
+License
+    This file is part of CFDEMcoupling.
+
+    CFDEMcoupling 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.
+
+    CFDEMcoupling 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 CFDEMcoupling; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Description
+    This code is designed to realize coupled CFD-DEM simulations using LIGGGHTS
+    and OpenFOAM(R). Note: this code is not part of OpenFOAM(R) (see DISCLAIMER).
+
+Class
+    global
+
+SourceFiles
+    global.Cver
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef global_H
+#define global_H
+
+#include "fvCFD.H"
+#include "cfdemCloud.H"
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class global Declaration
+\*---------------------------------------------------------------------------*/
+
+class global
+{
+
+protected:
+
+    // Protected data
+    const dictionary& dict_;
+
+    cfdemCloud& particleCloud_;
+
+    const char* const CFDEMversion;
+
+
+    // Protected member functions
+
+public:
+
+    //- Runtime type information
+    TypeName("global");
+
+    // Declare runtime constructor selection table
+
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            global,
+            dictionary,
+            (
+                const dictionary& dict,
+                cfdemCloud& sm
+            ),
+            (dict,sm)
+        );
+
+
+    // Constructors
+
+        //- Construct from components
+        global
+        (
+            const dictionary& dict,
+            cfdemCloud& sm
+        );
+
+
+    // Destructor
+
+        virtual ~global();
+
+
+    // Selector
+
+        static autoPtr<global> New
+        (
+            const dictionary& dict,
+            cfdemCloud& sm
+        );
+
+    // Member Function
+
+    void info();
+
+    // Access
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/cfdemParticle/cfdTools/mathExtra.H

@@ -0,0 +1,255 @@
+/* ----------------------------------------------------------------------
+    CFDEMcoupling - Open Source CFD-DEM coupling
+
+    CFDEMcoupling is part of the CFDEMproject
+    www.cfdem.com
+                                Christoph Goniva, christoph.goniva@cfdem.com
+                                Copyright 2009-2012 JKU Linz
+                                Copyright 2012-     DCS Computing GmbH, Linz
+-------------------------------------------------------------------------------
+License
+    This file is part of CFDEMcoupling.
+
+    CFDEMcoupling 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.
+
+    CFDEMcoupling 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 CFDEMcoupling; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Description
+    This code is designed to realize coupled CFD-DEM simulations using LIGGGHTS
+    and OpenFOAM(R). Note: this code is not part of OpenFOAM(R) (see DISCLAIMER).
+    
+    Copyright of this contribution:
+    Copyright 2014-     TU Graz, IPPT
+------------------------------------------------------------------------- */
+
+#ifndef CFDEM_MATH_EXTRA_H
+#define CFDEM_MATH_EXTRA_H
+
+#include <vector>
+//#include "math.h"
+#include <stdio.h>
+#include <string.h>
+#include <error.h>
+#include <ctype.h>
+
+#define TOLERANCE_ORTHO 1e-10
+
+namespace MathExtra 
+{
+
+//  inline void     outerProduct(double *vec1, double *vec2, double **m);
+//  inline double  spheroidGeometry(int index, double bi, double ai);
+
+
+
+//--------------------------------------------------------------------
+//   Outer Product of two vectors
+inline void outerProduct(double *vec1, double *vec2, double **m)
+{
+  int i, j;
+  //debug output
+//  for( i = 0; i < 3; ++i )
+//    printf("OUTER PRODUCT: Input: vec1 element %d = %g", i, vec1[i]);
+//  for( i = 0; i < 3; ++i )
+//    printf("OUTER PRODUCT: Input: vec2 element %d=%g", i, vec2[i]);
+
+  //calculation
+  for( i = 0; i < 3; ++i )
+    for( j = 0; j < 3; ++j )
+    {
+      m[i][j] = vec1[i] * vec2[j];
+      printf("OUTER PRODUCT: Result: m[%d][%d]=%g", i, j, m[i][j]);
+    }
+}
+//--------------------------------------------------------------------
+// Compute the major, minor axis and eccentricity parameters of a prolate spheroid
+inline bool spheroidGeometry(double radius, double aspectRatio, //inputs
+                             double& ai,    double& bi,         //outputs
+                             double& ei,    double& Le          //outputs
+                            )
+{
+    //INPUT
+    //  radius      ...volume-equivalent radius of the spheroid
+    //  aspectRatio ...major/minor aspect ratio
+    
+    //OUTPUT
+    //  ai  ... 
+    //  bi  ... 
+    //  ei  ... 
+    //  Le  ... 
+  
+    if(radius<=0.0)     //avoid troubles in case radius is 0 or negative
+        return false;
+
+    ai = radius * std::pow(aspectRatio*aspectRatio,0.33333333333333333333333);
+    bi = ai / aspectRatio;
+    ei = std::sqrt( 1.0 - 1.0 / (aspectRatio*aspectRatio) );
+    Le = std::log( (1.0+ei) / (1.0-ei) );
+
+    return true;
+}
+
+//--------------------------------------------------------------------
+// Compute the major, minor axis and eccentricity parameters of a prolate spheroid
+inline double Pi()
+{
+    return 3.1415926535897932384626433832795;
+}
+
+
+//--------------------------------------------------------------------
+// Compute the major, minor axis and eccentricity parameters of a prolate spheroid
+inline bool spheroidGeometry2(double radius,  double aspectRatio, //inputs
+                                 double& ai,  double& bi,         //outputs
+                                 double& XAe, double& YAe,        //outputs
+                                 double& XCe, double& YCe,        //outputs
+                                 double& YHe
+                                )
+{
+    //INPUT
+    //  radius      ...volume-equivalent radius of the spheroid
+    //  aspectRatio ...major/minor aspect ratio
+    
+    //OUTPUT
+    //  XAe  ...Eccentricity dependet parameter
+    //  YAe  ...Eccentricity dependet parameter
+    //  XCe  ...Eccentricity dependet parameter
+    //  XCe  ...Eccentricity dependet parameter
+    //  YCe  ...Eccentricity dependet parameter
+    //  YHe  ...Eccentricity dependet parameter
+
+    double ei(0.0), Le(0.0);
+    bool   result = spheroidGeometry(radius, aspectRatio, //inputs
+                                     ai,     bi,          //outputs
+                                     ei,     Le           //outputs
+                                    );
+    if(!result)
+        return false;
+
+    XAe= 2.6666666666666666666666667
+         *ei*ei*ei
+         /(-2.0*ei+(1.0+ei*ei)*Le);
+    YAe= 5.333333333333333333333333333
+        *ei*ei*ei
+        /(2.0*ei+(3.0*ei*ei-1.0)*Le);
+    XCe= 1.333333333333333333333333333
+         *ei*ei*ei
+         *(1.0-ei*ei)
+         /(2.0*ei-(1.0-ei*ei)*Le);
+    YCe= 1.3333333333333333333333
+         *ei*ei*ei
+         *(2.0-ei*ei)
+         /(-2.0*ei+(1.0+ei*ei)*Le);
+    YHe= 1.3333333333333333333333
+         *ei*ei*ei*ei*ei
+         /(-2.0*ei+(1.0+ei*ei)*Le);
+
+    return true;
+
+}
+
+//--------------------------------------------------------------------
+// zeroize a 3x3x3 tensor
+inline void zeroize333(double tensor[3][3][3] )
+{
+    for(int iX=0; iX<3; iX++)
+        for(int iY=0; iY<3; iY++)
+            for(int iZ=0; iZ<3; iZ++)
+                tensor[iX][iY][iZ] = 0.0;
+}
+
+//--------------------------------------------------------------------
+// zeroize a 3x3 tensor
+inline void zeroize33(double tensor[3][3] )
+{
+    for(int iX=0; iX<3; iX++)
+        for(int iY=0; iY<3; iY++)
+                tensor[iX][iY] = 0.0;
+}
+
+//--------------------------------------------------------------------
+// multiply a 3x3x3 tensor with a scalar
+inline void multiply333(double scalar, double tensor[3][3][3] )
+{
+    for(int iX=0; iX<3; iX++)
+        for(int iY=0; iY<3; iY++)
+            for(int iZ=0; iZ<3; iZ++)
+                tensor[iX][iY][iZ] *= scalar;
+}
+//--------------------------------------------------------------------
+// Compute a dot and dyadic product of with a vector
+inline void permutationTensor(double tensor[3][3][3] )
+{
+    zeroize333(tensor);
+    tensor[0][1][2] = 1.0;
+    tensor[1][2][0] = 1.0;
+    tensor[2][0][1] = 1.0;
+    tensor[0][2][1] =-1.0;
+    tensor[2][1][0] =-1.0;
+    tensor[1][0][2] =-1.0;
+}
+
+
+
+//--------------------------------------------------------------------
+// Compute a dot product of the permutation tensor and
+// then a dyadic product of with a vector
+inline bool permutationDotDyadic(
+                                 double vector[3],
+                                 double  tensor[3][3][3]
+                                )
+{
+    //Generate permutation tensor
+    double permutationT[3][3][3];
+    permutationTensor(permutationT);
+    
+    //Step 1: compute dot prodcut of permutation tensor
+    double permutationDotProd[3][3];
+    zeroize33(permutationDotProd);
+
+    for(int iX=0; iX<3; iX++)
+        for(int iY=0; iY<3; iY++)
+            for(int iZ=0; iZ<3; iZ++)
+                permutationDotProd[iX][iY] += permutationT[iX][iY][iZ]
+                                            * vector[iZ];
+
+    for(int iX=0; iX<3; iX++)
+        for(int iY=0; iY<3; iY++)
+            for(int iZ=0; iZ<3; iZ++)
+                tensor[iX][iY][iZ] = permutationDotProd[iX][iY]
+                                   * vector[iZ]; 
+    return true;
+
+}
+
+//--------------------------------------------------------------------
+// Compute a dot and dyadic product of with a vector
+inline bool doubleDotTensor333Tensor33(double tensor333[3][3][3],
+                                       double tensor33[3][3],
+                                       double result[3]
+                                      )
+{
+    result[0]=0.0;result[1]=0.0;result[2]=0.0;
+    
+    for(int iX=0; iX<3; iX++)
+        for(int iY=0; iY<3; iY++)
+            for(int iZ=0; iZ<3; iZ++)
+                result[iX] += tensor333[iX][iY][iZ] * tensor33[iY][iZ];
+
+    return true;
+}
+
+
+} //end of namespace
+
+#endif

src/lagrangian/cfdemParticle/cfdTools/newGlobal.C

@@ -0,0 +1,84 @@
+/*---------------------------------------------------------------------------*\
+    CFDEMcoupling - Open Source CFD-DEM coupling
+
+    CFDEMcoupling is part of the CFDEMproject
+    www.cfdem.com
+                                Christoph Goniva, christoph.goniva@cfdem.com
+                                Copyright 2009-2012 JKU Linz
+                                Copyright 2012-     DCS Computing GmbH, Linz
+-------------------------------------------------------------------------------
+License
+    This file is part of CFDEMcoupling.
+
+    CFDEMcoupling 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.
+
+    CFDEMcoupling 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 CFDEMcoupling; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Description
+    This code is designed to realize coupled CFD-DEM simulations using LIGGGHTS
+    and OpenFOAM(R). Note: this code is not part of OpenFOAM(R) (see DISCLAIMER).
+\*---------------------------------------------------------------------------*/
+
+#include "error.H"
+
+#include "global.H"
+#include "dilute.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+autoPtr<global> global::New
+(
+    const dictionary& dict,
+    cfdemCloud& sm
+)
+{
+    word globalType
+    (
+        dict.lookup("global")
+    );
+
+    Info<< "Selecting global "
+         << globalType << endl;
+
+
+    dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(globalType);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalError
+            << "global::New(const dictionary&, const spray&) : "
+            << endl
+            << "    unknown globalType type "
+            << globalType
+            << ", constructor not in hash table" << endl << endl
+            << "    Valid global types are :"
+            << endl;
+        Info<< dictionaryConstructorTablePtr_->toc()
+            << abort(FatalError);
+    }
+
+    return autoPtr<global>(cstrIter()(dict,sm));
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/lagrangian/cfdemParticle/cfdTools/versionInfo.H

@@ -1,9 +1,45 @@
-word CFDEMversion="cfdem-2.6.4";
-word compatibleLIGGGHTSversion="3.0.1";
-word OFversion="2.2.x-commit-61b850bc107bdd60bbf1bf9a6417b9faf701d128";
+/*---------------------------------------------------------------------------*\
+    CFDEMcoupling - Open Source CFD-DEM coupling
 
-Info << "\nCFDEMcoupling version: " << CFDEMversion << "\n" <<  endl;
-Info << "\n, compatible to LIGGGHTS version: " << compatibleLIGGGHTSversion << "\n" <<  endl;
-Info << "\n, compatible to OF version: " << OFversion << "\n" <<  endl;
+    CFDEMcoupling is part of the CFDEMproject
+    www.cfdem.com
+                                Christoph Goniva, christoph.goniva@cfdem.com
+                                Copyright 2009-2012 JKU Linz
+                                Copyright 2012-2015 DCS Computing GmbH, Linz
+                                Copyright 2015-     JKU Linz
+-------------------------------------------------------------------------------
+License
+    This file is part of CFDEMcoupling.
 
+    CFDEMcoupling 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.
 
+    CFDEMcoupling 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 CFDEMcoupling; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Description
+    This code is designed to realize coupled CFD-DEM simulations using LIGGGHTS
+    and OpenFOAM(R). Note: this code is not part of OpenFOAM(R) (see DISCLAIMER).
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef versionInfo_H
+#define versionInfo_H
+
+word CFDEMversion="PFM 17.02";
+word compatibleLIGGGHTSversion="PFM 17.02";
+word OFversion="4.x";
+
+Info << "\nCFDEMcoupling version: " << CFDEMversion <<  endl;
+Info << ", compatible to LIGGGHTS version: " << compatibleLIGGGHTSversion <<  endl;
+Info << ", compatible to OpenFOAM version: " << OFversion <<  endl;
+
+#endif

src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.C

@@ -31,6 +31,7 @@ Description
 
 #include "fileName.H"
 #include "cfdemCloud.H"
+#include "global.H"
 #include "forceModel.H"
 #include "locateModel.H"
 #include "momCoupleModel.H"
@@ -43,9 +44,13 @@ Description
 #include "clockModel.H"
 #include "smoothingModel.H"
 #include "liggghtsCommandModel.H"
+#include "otherForceModel.H"
+
+namespace Foam
+{
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-Foam::cfdemCloud::cfdemCloud
+cfdemCloud::cfdemCloud
 (
     const fvMesh& mesh
 )
@@ -76,10 +81,12 @@ Foam::cfdemCloud::cfdemCloud
     solveFlow_(true),
     verbose_(false),
     ignore_(false),
+    limitDEMForces_(false),
     modelType_(couplingProperties_.lookup("modelType")),
     positions_(NULL),
     velocities_(NULL),
     fluidVel_(NULL),
+    fAcc_(NULL),
     impForces_(NULL),
     expForces_(NULL),
     DEMForces_(NULL),
@@ -89,20 +96,25 @@ Foam::cfdemCloud::cfdemCloud
     cellIDs_(NULL),
     particleWeights_(NULL),
     particleVolumes_(NULL),
+    particleV_(NULL),
     numberOfParticles_(0),
+    d32_(-1),
     numberOfParticlesChanged_(false),
     arraysReallocated_(false),
     forceModels_(couplingProperties_.lookup("forceModels")),
     momCoupleModels_(couplingProperties_.lookup("momCoupleModels")),
     liggghtsCommandModelList_(liggghtsCommandDict_.lookup("liggghtsCommandModels")),
+    otherForceModels_(couplingProperties_.lookupOrDefault<wordList>("otherForceModels",wordList(0))),
     turbulenceModelType_(couplingProperties_.lookup("turbulenceModelType")),
     cg_(1.),
     cgOK_(true),
     impDEMdrag_(false),
+    impDEMdragAcc_(false),
     imExSplitFactor_(1.0),
+    treatVoidCellsAsExplicitForce_(false),
     useDDTvoidfraction_(false),
     ddtVoidfraction_
-    (   
+    (
         IOobject
         (
             "ddtVoidfraction",
@@ -116,15 +128,7 @@ Foam::cfdemCloud::cfdemCloud
     ),
     turbulence_
     (
-        #if defined(version21) || defined(version16ext)
-            #ifdef comp
-                mesh.lookupObject<compressible::turbulenceModel>
-            #else
-                mesh.lookupObject<incompressible::turbulenceModel>
-            #endif
-        #elif defined(version15)
-            mesh.lookupObject<incompressible::RASModel>
-        #endif
+        mesh.lookupObject<turbulenceModel>
         (
             turbulenceModelType_
         )
@@ -213,20 +217,29 @@ Foam::cfdemCloud::cfdemCloud
     )
 {
     #include "versionInfo.H"
+    global buildInfo(couplingProperties_,*this);
+    buildInfo.info();
 
     Info << "If BC are important, please provide volScalarFields -imp/expParticleForces-" << endl;
     if (couplingProperties_.found("solveFlow"))
         solveFlow_=Switch(couplingProperties_.lookup("solveFlow"));
     if (couplingProperties_.found("imExSplitFactor"))
         imExSplitFactor_ = readScalar(couplingProperties_.lookup("imExSplitFactor"));
+    if (couplingProperties_.found("treatVoidCellsAsExplicitForce"))
+        treatVoidCellsAsExplicitForce_ = readBool(couplingProperties_.lookup("treatVoidCellsAsExplicitForce"));
     if (couplingProperties_.found("verbose")) verbose_=true;
     if (couplingProperties_.found("ignore")) ignore_=true;
+    if (couplingProperties_.found("limitDEMForces"))
+    {
+        limitDEMForces_=true;
+        maxDEMForce_ = readScalar(couplingProperties_.lookup("limitDEMForces"));
+    }
     if (turbulenceModelType_=="LESProperties")
         Info << "WARNING - LES functionality not yet tested!" << endl;
 
     if (couplingProperties_.found("useDDTvoidfraction"))
         useDDTvoidfraction_=true;
-    else        
+    else
         Info << "ignoring ddt(voidfraction)" << endl;
 
     forceModel_ = new autoPtr<forceModel>[nrForceModels()];
@@ -264,18 +277,30 @@ Foam::cfdemCloud::cfdemCloud
         );
     }
 
+    otherForceModel_ = new autoPtr<otherForceModel>[otherForceModels_.size()];
+    for (int i=0;i<otherForceModels_.size();i++)
+    {
+        otherForceModel_[i] = otherForceModel::New
+        (
+            couplingProperties_,
+            *this,
+            otherForceModels_[i]
+        );
+    }
+
     dataExchangeM().setCG();
-    if (!cgOK_ && cg_ > 1) FatalError<< "at least one of your models is not fit for cg !!!"<< abort(FatalError); 
+    if (!cgOK_ && cg_ > 1) FatalError<< "at least one of your models is not fit for cg !!!"<< abort(FatalError);
 }
 
 // * * * * * * * * * * * * * * * * Destructors  * * * * * * * * * * * * * * //
-Foam::cfdemCloud::~cfdemCloud()
+cfdemCloud::~cfdemCloud()
 {
     clockM().evalPar();
     clockM().normHist();
     dataExchangeM().destroy(positions_,3);
     dataExchangeM().destroy(velocities_,3);
     dataExchangeM().destroy(fluidVel_,3);
+    dataExchangeM().destroy(fAcc_,3);
     dataExchangeM().destroy(impForces_,3);
     dataExchangeM().destroy(expForces_,3);
     dataExchangeM().destroy(DEMForces_,3);
@@ -285,16 +310,21 @@ Foam::cfdemCloud::~cfdemCloud()
     dataExchangeM().destroy(cellIDs_,1);
     dataExchangeM().destroy(particleWeights_,1);
     dataExchangeM().destroy(particleVolumes_,1);
+    dataExchangeM().destroy(particleV_,1);
 }
+
 // * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
-void Foam::cfdemCloud::getDEMdata()
+void cfdemCloud::getDEMdata()
 {
     dataExchangeM().getData("radius","scalar-atom",radii_);
     dataExchangeM().getData("x","vector-atom",positions_);
     dataExchangeM().getData("v","vector-atom",velocities_);
+
+    if(impDEMdragAcc_)
+        dataExchangeM().getData("dragAcc","vector-atom",fAcc_); // array is used twice - might be necessary to clean it first
 }
 
-void Foam::cfdemCloud::giveDEMdata()
+void cfdemCloud::giveDEMdata()
 {
     if(forceM(0).coupleForce())
     {
@@ -302,6 +332,7 @@ void Foam::cfdemCloud::giveDEMdata()
 
         if(impDEMdrag_)
         {
+            if(verbose_) Info << "sending Ksl and uf" << endl;
             dataExchangeM().giveData("Ksl","scalar-atom",Cds_);
             dataExchangeM().giveData("uf","vector-atom",fluidVel_);
         }
@@ -313,7 +344,7 @@ void Foam::cfdemCloud::giveDEMdata()
 
 // * * * * * * * * * * * * * * * protected Member Functions  * * * * * * * * * * * * * //
 
-void Foam::cfdemCloud::setNumberOfParticles(int nP)
+void cfdemCloud::setNumberOfParticles(int nP)
 {
     if(nP != numberOfParticles())
     {
@@ -322,12 +353,12 @@ void Foam::cfdemCloud::setNumberOfParticles(int nP)
     }
 }
 
-void Foam::cfdemCloud::findCells()
+void cfdemCloud::findCells()
 {
     locateM().findCell(NULL,positions_,cellIDs_,numberOfParticles());
 }
 
-void Foam::cfdemCloud::setForces()
+void cfdemCloud::setForces()
 {
     resetArray(fluidVel_,numberOfParticles(),3);
     resetArray(impForces_,numberOfParticles(),3);
@@ -335,9 +366,22 @@ void Foam::cfdemCloud::setForces()
     resetArray(DEMForces_,numberOfParticles(),3);
     resetArray(Cds_,numberOfParticles(),1);
     for (int i=0;i<cfdemCloud::nrForceModels();i++) cfdemCloud::forceM(i).setForce();
+
+    if (limitDEMForces_)
+    {
+        scalar maxF = 0.0;
+        for (int index = 0;index <  numberOfParticles(); ++index)
+        {
+            scalar F = mag(expForce(index));
+            if (F > maxF) maxF = F;
+            if (F > maxDEMForce_)
+              for(int i=0;i<3;i++) DEMForces_[index][i] *= maxDEMForce_/F;
+        }
+        Info << "largest particle-fluid interaction on particle: " << maxF << endl;
+    }
 }
 
-void Foam::cfdemCloud::setParticleForceField()
+void cfdemCloud::setParticleForceField()
 {
     averagingM().setVectorSum
     (
@@ -355,7 +399,7 @@ void Foam::cfdemCloud::setParticleForceField()
     );
 }
 
-void Foam::cfdemCloud::setVectorAverages()
+void cfdemCloud::setVectorAverages()
 {
     if(verbose_) Info << "- setVectorAverage(Us,velocities_,weights_)" << endl;
     averagingM().setVectorAverage
@@ -369,12 +413,12 @@ void Foam::cfdemCloud::setVectorAverages()
     if(verbose_) Info << "setVectorAverage done." << endl;
 }
 // * * * * * * * * * * * * * * * public Member Functions  * * * * * * * * * * * * * //
-void Foam::cfdemCloud::checkCG(bool ok)
+void cfdemCloud::checkCG(bool ok)
 {
     if(!cgOK_) return;
     if(!ok) cgOK_ = ok;
 }
-void Foam::cfdemCloud::setPos(double**& pos)
+void cfdemCloud::setPos(double**& pos)
 {
     for(int index = 0;index <  numberOfParticles(); ++index)
     {
@@ -385,60 +429,61 @@ void Foam::cfdemCloud::setPos(double**& pos)
 }
 // * * * * * * * * * * * * * * * ACCESS  * * * * * * * * * * * * * //
 
-label Foam::cfdemCloud::particleCell(int index)
+label cfdemCloud::particleCell(int index)
 {
     label cellI = cellIDs()[index][0];
     return cellI;
 }
 
-double Foam::cfdemCloud::d(int index)
-{
-    return 2*radii()[index][0];
-}
-
-vector Foam::cfdemCloud::position(int index)
+vector cfdemCloud::position(int index)
 {
     vector pos;
     for(int i=0;i<3;i++) pos[i] = positions()[index][i];
     return pos;
 }
 
-vector Foam::cfdemCloud::velocity(int index)
+vector cfdemCloud::velocity(int index)
 {
     vector vel;
     for(int i=0;i<3;i++) vel[i] = velocities()[index][i];
     return vel;
 }
 
-vector Foam::cfdemCloud::fluidVel(int index)
+vector cfdemCloud::expForce(int index)
+{
+    vector force;
+    for(int i=0;i<3;i++) force[i] = DEMForces()[index][i];
+    return force;
+}
+
+vector cfdemCloud::fluidVel(int index)
 {
     vector vel;
     for(int i=0;i<3;i++) vel[i] = fluidVels()[index][i];
     return vel;
 }
 
-const forceModel& Foam::cfdemCloud::forceM(int i)
+const forceModel& cfdemCloud::forceM(int i)
 {
     return forceModel_[i];
 }
 
-int Foam::cfdemCloud::nrForceModels()
+int cfdemCloud::nrForceModels()
 {
     return forceModels_.size();
 }
 
-scalar Foam::cfdemCloud::radius(int index)
+int cfdemCloud::nrMomCoupleModels()
 {
-    scalar r = radii()[index][0];
-    return r;
+    return momCoupleModels_.size();
 }
 
-scalar Foam::cfdemCloud::voidfraction(int index)
+scalar cfdemCloud::voidfraction(int index)
 {
     return voidfractions()[index][0];
 }
 
-label Foam::cfdemCloud::liggghtsCommandModelIndex(word name)
+label cfdemCloud::liggghtsCommandModelIndex(word name)
 {
     int index=-1;
     forAll(liggghtsCommandModelList_,i)
@@ -452,11 +497,21 @@ label Foam::cfdemCloud::liggghtsCommandModelIndex(word name)
     return index;
 }
 
+std::vector< std::vector<double*> >* cfdemCloud::getVprobe()
+{
+ return probeModel_->getVprobe();
+}
+
+std::vector< std::vector<double> >* cfdemCloud::getSprobe()
+{
+ return probeModel_->getSprobe();
+}
+
 // * * * * * * * * * * * * * * * WRITE  * * * * * * * * * * * * * //
 
 // * * *   write cfdemCloud internal data   * * * //
 
-bool Foam::cfdemCloud::evolve
+bool cfdemCloud::evolve
 (
     volScalarField& alpha,
     volVectorField& Us,
@@ -469,16 +524,17 @@ bool Foam::cfdemCloud::evolve
 
     if(!ignore())
     {
-        if (dataExchangeM().couple())
+        if (dataExchangeM().doCoupleNow())
         {
             Info << "\n Coupling..." << endl;
+            dataExchangeM().couple(0);
             doCouple=true;
 
             // reset vol Fields
             clockM().start(16,"resetVolFields");
             if(verbose_)
             {
-                Info << "couplingStep:" << dataExchangeM().couplingStep() 
+                Info << "couplingStep:" << dataExchangeM().couplingStep()
                      << "\n- resetVolFields()" << endl;
             }
             averagingM().resetVectorAverage(averagingM().UsPrev(),averagingM().UsNext(),false);
@@ -507,7 +563,7 @@ bool Foam::cfdemCloud::evolve
             // set void fraction field
             clockM().start(19,"setvoidFraction");
             if(verbose_) Info << "- setvoidFraction()" << endl;
-            voidFractionM().setvoidFraction(NULL,voidfractions_,particleWeights_,particleVolumes_);
+            voidFractionM().setvoidFraction(NULL,voidfractions_,particleWeights_,particleVolumes_,particleV_);
             if(verbose_) Info << "setvoidFraction done." << endl;
             clockM().stop("setvoidFraction");
 
@@ -516,16 +572,21 @@ bool Foam::cfdemCloud::evolve
             setVectorAverages();
 
 
-            //Smoothen "next" fields            
+            //Smoothen "next" fields
             smoothingM().dSmoothing();
             smoothingM().smoothen(voidFractionM().voidFractionNext());
-            smoothingM().smoothenReferenceField(averagingM().UsNext());
-            
+
+            //only smoothen if we use implicit force coupling in cells void of particles
+            //because we need unsmoothened Us field to detect cells for explicit
+            //force coupling
+            if(!treatVoidCellsAsExplicitForce())
+                smoothingM().smoothenReferenceField(averagingM().UsNext());
+
             clockM().stop("setVectorAverage");
         }
-        
+
         //============================================
-        //CHECK JUST TIME-INTERPOATE ALREADY SMOOTHENED VOIDFRACTIONNEXT AND UsNEXT FIELD 
+        //CHECK JUST TIME-INTERPOATE ALREADY SMOOTHENED VOIDFRACTIONNEXT AND UsNEXT FIELD
         //      IMPLICIT FORCE CONTRIBUTION AND SOLVER USE EXACTLY THE SAME AVERAGED
         //      QUANTITIES AT THE GRID!
         Info << "\n timeStepFraction() = " << dataExchangeM().timeStepFraction() << endl;
@@ -571,6 +632,8 @@ bool Foam::cfdemCloud::evolve
             clockM().start(23,"giveDEMdata");
             giveDEMdata();
             clockM().stop("giveDEMdata");
+
+            dataExchangeM().couple(1);
         }//end dataExchangeM().couple()
 
 
@@ -587,7 +650,7 @@ bool Foam::cfdemCloud::evolve
     return doCouple;
 }
 
-bool Foam::cfdemCloud::reAllocArrays() const
+bool cfdemCloud::reAllocArrays()
 {
     if(numberOfParticlesChanged_ && !arraysReallocated_)
     {
@@ -595,15 +658,40 @@ bool Foam::cfdemCloud::reAllocArrays() const
         dataExchangeM().allocateArray(positions_,0.,3);
         dataExchangeM().allocateArray(velocities_,0.,3);
         dataExchangeM().allocateArray(fluidVel_,0.,3);
+        dataExchangeM().allocateArray(fAcc_,0.,3);
         dataExchangeM().allocateArray(impForces_,0.,3);
         dataExchangeM().allocateArray(expForces_,0.,3);
         dataExchangeM().allocateArray(DEMForces_,0.,3);
         dataExchangeM().allocateArray(Cds_,0.,1);
         dataExchangeM().allocateArray(radii_,0.,1);
         dataExchangeM().allocateArray(voidfractions_,1.,voidFractionM().maxCellsPerParticle());
-        dataExchangeM().allocateArray(cellIDs_,0.,voidFractionM().maxCellsPerParticle());
+        dataExchangeM().allocateArray(cellIDs_,-1,voidFractionM().maxCellsPerParticle());
         dataExchangeM().allocateArray(particleWeights_,0.,voidFractionM().maxCellsPerParticle());
         dataExchangeM().allocateArray(particleVolumes_,0.,voidFractionM().maxCellsPerParticle());
+        dataExchangeM().allocateArray(particleV_,0.,1);
+        arraysReallocated_ = true;
+        return true;
+    }
+    return false;
+}
+
+bool cfdemCloud::reAllocArrays(int nP, bool forceRealloc)
+{
+    if( (numberOfParticlesChanged_ && !arraysReallocated_) || forceRealloc)
+    {
+        // get arrays of new length
+        dataExchangeM().allocateArray(positions_,0.,3,nP);
+        dataExchangeM().allocateArray(velocities_,0.,3,nP);
+        dataExchangeM().allocateArray(fluidVel_,0.,3,nP);
+        dataExchangeM().allocateArray(impForces_,0.,3,nP);
+        dataExchangeM().allocateArray(expForces_,0.,3,nP);
+        dataExchangeM().allocateArray(DEMForces_,0.,3,nP);
+        dataExchangeM().allocateArray(Cds_,0.,1,nP);
+        dataExchangeM().allocateArray(radii_,0.,1,nP);
+        dataExchangeM().allocateArray(voidfractions_,1.,voidFractionM().maxCellsPerParticle(),nP);
+        dataExchangeM().allocateArray(cellIDs_,-1,voidFractionM().maxCellsPerParticle(),nP);
+        dataExchangeM().allocateArray(particleWeights_,0.,voidFractionM().maxCellsPerParticle(),nP);
+        dataExchangeM().allocateArray(particleVolumes_,0.,voidFractionM().maxCellsPerParticle(),nP);
         arraysReallocated_ = true;
         return true;
     }
@@ -615,7 +703,7 @@ tmp<fvVectorMatrix> cfdemCloud::divVoidfractionTau(volVectorField& U,volScalarFi
     return
     (
       - fvm::laplacian(voidfractionNuEff(voidfraction), U)
-      - fvc::div(voidfractionNuEff(voidfraction)*dev(fvc::grad(U)().T()))
+      - fvc::div(voidfractionNuEff(voidfraction)*dev2(fvc::grad(U)().T()))
     );
 }
 
@@ -642,17 +730,17 @@ void cfdemCloud::calcDdtVoidfraction(volScalarField& voidfraction) const
 /*tmp<fvVectorMatrix> cfdemCloud::ddtVoidfractionU(volVectorField& U,volScalarField& voidfraction) const
 {
     if (dataExchangeM().couplingStep() <= 2) return fvm::ddt(U);
-    
+
     return fvm::ddt(voidfraction,U);
 }*/
 
 tmp<volScalarField> cfdemCloud::voidfractionNuEff(volScalarField& voidfraction) const
 {
-    if (modelType_=="B")
+    if (modelType_=="B" || modelType_=="Bfull")
     {
         return tmp<volScalarField>
         (
-            #ifdef comp
+            #ifdef compre
                 new volScalarField("viscousTerm", (turbulence_.mut() + turbulence_.mu()))
             #else
                 new volScalarField("viscousTerm", (turbulence_.nut() + turbulence_.nu()))
@@ -663,7 +751,7 @@ tmp<volScalarField> cfdemCloud::voidfractionNuEff(volScalarField& voidfraction)
     {
         return tmp<volScalarField>
         (
-            #ifdef comp
+            #ifdef compre
                 new volScalarField("viscousTerm", voidfraction*(turbulence_.mut() + turbulence_.mu()))
             #else
                 new volScalarField("viscousTerm", voidfraction*(turbulence_.nut() + turbulence_.nu()))
@@ -680,8 +768,17 @@ void cfdemCloud::resetArray(double**& array,int length,int width,double resetVal
         }
     }
 }
-// * * * * * * * * * * * * * * * *  IOStream operators * * * * * * * * * * * //
 
-#include "cfdemCloudIO.C"
+void cfdemCloud::otherForces(volVectorField& forcefield)
+{
+  forcefield.primitiveFieldRef() = vector::zero;
+  forcefield.boundaryFieldRef() = vector::zero;
+  for (int i=0;i<otherForceModels_.size();i++)
+      forcefield += otherForceModel_[i]().exportForceField();
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
 
 // ************************************************************************* //