Merge pull request #154 from ParticulateFlow/release

Release 24.01

.gitignore

@@ -8,7 +8,7 @@ log.*
 *.swp
 *.swo
 
-**/linux*Gcc*/
+**/linux*cc*/
 **/.vscode
 
 lnInclude

README.md

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

applications/solvers/cfdemSolverIB/Make/options

@@ -1,7 +1,6 @@
 include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
 
 EXE_INC = \
-    -I$(CFDEM_OFVERSION_DIR) \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
@@ -9,12 +8,12 @@ EXE_INC = \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
     -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/lnInclude \
-    -I$(LIB_SRC)/sampling/lnInclude \
     -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/dynamicFvMesh/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/fvOptions/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
     -Wno-deprecated-copy
 
 EXE_LIBS = \
@@ -27,6 +26,7 @@ EXE_LIBS = \
     -ldynamicFvMesh \
     -ldynamicMesh \
     -lfvOptions \
+    -lsampling \
     -l$(CFDEM_LIB_NAME) \
      $(CFDEM_ADD_LIB_PATHS) \
      $(CFDEM_ADD_LIBS)

applications/solvers/cfdemSolverIB/cfdemSolverIB.C

@@ -6,7 +6,8 @@
                                 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
+                                Copyright (C) 2012-2015 DCS Computing GmbH,Linz
+                                Copyright (C) 2015-     JKU, Linz
 -------------------------------------------------------------------------------
 License
     This file is part of CFDEMcoupling.
@@ -29,11 +30,14 @@ Application
 
 Description
     Transient solver for incompressible flow.
-    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 using immersed body
     (fictitious domain) method is added.
 Contributions
     Alice Hager
+    Daniel Queteschiner
+    Thomas Lichtenegger
+    Achuth N. Balachandran Nair
 \*---------------------------------------------------------------------------*/
 
 
@@ -53,23 +57,21 @@ Contributions
 
 #include "cellSet.H"
 
+#include "fvOptions.H"  // added the fvOptions library
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 int main(int argc, char *argv[])
 {
+
     #include "setRootCase.H"
-
     #include "createTime.H"
-
     #include "createDynamicFvMesh.H"
-
     #include "createControl.H"
-
     #include "createTimeControls.H"
-
     #include "createFields.H"
-
     #include "initContinuityErrs.H"
+    #include "createFvOptions.H"
 
     // create cfdemCloud
     #include "readGravitationalAcceleration.H"
@@ -93,24 +95,31 @@ int main(int argc, char *argv[])
 
         // do particle stuff
         Info << "- evolve()" << endl;
-        particleCloud.evolve();
+        particleCloud.evolve(Us);
 
         // Pressure-velocity PISO corrector
         {
+            MRF.correctBoundaryVelocity(U);
+
             // Momentum predictor
 
             fvVectorMatrix UEqn
             (
-                fvm::ddt(voidfraction,U)
+                fvm::ddt(voidfraction,U) + MRF.DDt(U)
               + fvm::div(phi, U)
               + turbulence->divDevReff(U)
+             ==
+                fvOptions(U)
             );
 
             UEqn.relax();
 
+            fvOptions.constrain(UEqn);
+
             if (piso.momentumPredictor())
             {
                 solve(UEqn == -fvc::grad(p));
+                fvOptions.correct(U);
             }
 
             // --- PISO loop
@@ -126,6 +135,7 @@ int main(int argc, char *argv[])
 
                 adjustPhi(phi, U, p);
 
+
                 while (piso.correctNonOrthogonal())
                 {
                     // Pressure corrector
@@ -152,12 +162,15 @@ int main(int argc, char *argv[])
             }
         }
 
+        laminarTransport.correct();
         turbulence->correct();
 
         Info << "particleCloud.calcVelocityCorrection() " << endl;
         volScalarField voidfractionNext=mesh.lookupObject<volScalarField>("voidfractionNext");
         particleCloud.calcVelocityCorrection(p,U,phiIB,voidfractionNext);
 
+        fvOptions.correct(U);
+
         runTime.write();
 
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"

applications/solvers/cfdemSolverIB/createFields.H

@@ -26,21 +26,6 @@
         ),
         mesh
     );
-    //mod by alice
-    Info<< "Reading physical velocity field U" << endl;
-    Info<< "Note: only if voidfraction at boundary is 1, U is superficial velocity!!!\n" << endl;
-    volVectorField Us
-    (
-        IOobject
-        (
-            "Us",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
 
 //========================
 // drag law modelling
@@ -76,9 +61,8 @@
         mesh
     );
 
-
     //mod by alice
-    Info<< "Reading field phiIB\n" << endl;
+    Info<< "Reading field voidfraction\n" << endl;
     volScalarField voidfraction
     (
         IOobject
@@ -91,6 +75,21 @@
         ),
         mesh
     );
+
+    Info<< "Reading particle velocity field Us\n" << endl;
+    volVectorField Us
+    (
+        IOobject
+        (
+            "Us",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
 //========================
 
 #   include "createPhi.H"
@@ -126,3 +125,5 @@
     );
 
 //===========================
+
+#include "createMRF.H"

applications/solvers/cfdemSolverIBContinuousForcing/Make/files

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

applications/solvers/cfdemSolverIBContinuousForcing/Make/options

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

applications/solvers/cfdemSolverIBContinuousForcing/UEqn.H

@@ -0,0 +1,19 @@
+fvVectorMatrix UEqn
+(
+    fvm::ddt(voidfractionNext,U) + MRF.DDt(U)
+  + fvm::div(phi, U)
+  + turbulence->divDevReff(U)
+ ==
+    fvOptions(U)
+ +  (lambda*(1-voidfractionNext)/U.mesh().time().deltaT())*(fvc::Sp(1,Us)-fvm::Sp(1,U))
+);
+
+UEqn.relax();
+
+fvOptions.constrain(UEqn);
+
+if (piso.momentumPredictor())
+{
+    solve(UEqn == -fvc::grad(p));
+    fvOptions.correct(U);
+}

applications/solvers/cfdemSolverIBContinuousForcing/cfdemSolverIBContinuousForcing.C

@@ -0,0 +1,129 @@
+/*---------------------------------------------------------------------------*\
+    CFDEMcoupling - Open Source CFD-DEM coupling
+
+    CFDEMcoupling is part of the CFDEMproject
+    www.cfdem.com
+                                Copyright (C) 1991-2009 OpenCFD Ltd.
+                                Copyright (C) 2009-2012 JKU, Linz
+                                Copyright (C) 2012-2015 DCS Computing GmbH,Linz
+                                Copyright (C) 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, see <http://www.gnu.org/licenses/>.
+
+Application
+    cfdemSolverIBContinuousForcing
+
+Description
+    Transient solver for incompressible flow.
+    The code is an evolution of the solver pisoFoam in OpenFOAM(R) 1.6,
+    where additional functionality for CFD-DEM coupling using immersed body
+    (fictitious domain) method and a continuous forcing approach is added.
+Contributions
+    Alice Hager
+    Achuth N. Balachandran Nair
+\*---------------------------------------------------------------------------*/
+
+
+#include "fvCFD.H"
+#include "singlePhaseTransportModel.H"
+#include "turbulentTransportModel.H"
+#include "pisoControl.H"
+
+#include "cfdemCloudIB.H"
+#include "implicitCouple.H"
+
+#include "averagingModel.H"
+#include "regionModel.H"
+#include "voidFractionModel.H"
+
+#include "dynamicFvMesh.H"
+
+#include "cellSet.H"
+
+#include "fvOptions.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createDynamicFvMesh.H"
+    #include "createControl.H"
+    #include "createTimeControls.H"
+    #include "createFields.H"
+    #include "initContinuityErrs.H"
+    #include "createFvOptions.H"
+
+    // create cfdemCloud
+    #include "readGravitationalAcceleration.H"
+    cfdemCloudIB particleCloud(mesh);
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (runTime.loop())
+    {
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        //=== dyM ===================
+        interFace = mag(mesh.lookupObject<volScalarField>("voidfractionNext"));
+        mesh.update(); //dyM
+
+        #include "readTimeControls.H"
+        #include "CourantNo.H"
+        #include "setDeltaT.H"
+
+        // do particle stuff
+        Info << "- evolve()" << endl;
+        particleCloud.evolve(Us);
+
+        volScalarField voidfractionNext=mesh.lookupObject<volScalarField>("voidfractionNext");
+
+        // Pressure-velocity PISO corrector
+        {
+            MRF.correctBoundaryVelocity(U);
+
+            // Momentum predictor
+            #include "UEqn.H"
+
+            // --- PISO loop
+            while (piso.correct())
+            {
+                #include "pEqn.H"
+            }
+        }
+
+        laminarTransport.correct();
+        turbulence->correct();
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/cfdemSolverIBContinuousForcing/createFields.H

@@ -0,0 +1,143 @@
+    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
+    );
+
+    Info<< "Reading particle velocity field Us\n" << endl;
+    volVectorField Us
+    (
+        IOobject
+        (
+            "Us",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< "Reading the penalization factor field lambda\n" << endl;
+    volScalarField lambda
+    (
+        IOobject
+        (
+            "lambda",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+//========================
+// drag law modelling
+//========================
+
+    Info<< "\nCreating dummy density field rho = 1\n" << endl;
+    volScalarField rho
+    (
+        IOobject
+        (
+            "rho",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("0", dimensionSet(1, -3, 0, 0, 0), 1.0)
+    );
+
+
+    Info<< "Reading field phiIB\n" << endl;
+    volScalarField phiIB
+    (
+        IOobject
+        (
+            "phiIB",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    //mod by alice
+    Info<< "Reading field voidfraction\n" << endl;
+    volScalarField voidfraction
+    (
+        IOobject
+        (
+            "voidfraction",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+//========================
+
+#   include "createPhi.H"
+
+    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)
+    );
+
+//=== dyM ===================
+
+    Info<< "Reading field interFace\n" << endl;
+    volScalarField interFace
+    (
+        IOobject
+        (
+            "interFace",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        //dimensionedScalar("0", dimensionSet(0, -1, 0, 0, 0), 0.0)
+        dimensionedScalar("0", dimensionSet(0, 0, 0, 0, 0), 0.0)
+    );
+
+//===========================
+
+#include "createMRF.H"

applications/solvers/cfdemSolverIBContinuousForcing/pEqn.H

@@ -0,0 +1,35 @@
+volScalarField rUA = 1.0/UEqn.A();
+surfaceScalarField rUAf(fvc::interpolate(rUA));
+
+U = rUA*UEqn.H();
+
+phi = (fvc::interpolate(U) & mesh.Sf())
+    + rUAf*fvc::ddtCorr(U, phi);            // Is there additional flux term due to the particle presence?
+
+MRF.makeRelative(phi);
+
+adjustPhi(phi, U, p);
+
+while (piso.correctNonOrthogonal())
+{
+    // Pressure corrector
+
+    fvScalarMatrix pEqn
+    (
+        fvm::laplacian(rUA, 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 "continuityErrs.H"
+U -= rUA*fvc::grad(p); // should we add a pressure correction?
+U.correctBoundaryConditions();
+fvOptions.correct(U);

applications/solvers/cfdemSolverMultiphase/Make/options

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

applications/solvers/cfdemSolverMultiphaseScalar/Allwclean

@@ -2,7 +2,7 @@
 cd ${0%/*} || exit 1    # Run from this directory
 set -x
 
-wclean libso multiphaseMixture
+wclean libso multiphaseMixtureScalar
 wclean
 
 #------------------------------------------------------------------------------

applications/solvers/cfdemSolverMultiphaseScalar/Allwmake

@@ -6,7 +6,7 @@ targetType=libso
 . $WM_PROJECT_DIR/wmake/scripts/AllwmakeParseArguments
 set -x
 
-wmake $targetType multiphaseMixture
+wmake $targetType multiphaseMixtureScalar
 wmake
 
 #------------------------------------------------------------------------------

applications/solvers/cfdemSolverMultiphaseScalar/Make/options

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

applications/solvers/cfdemSolverMultiphaseScalar/cfdemSolverMultiphaseScalar.C

@@ -30,7 +30,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
-#include "multiphaseMixture.H"
+#include "multiphaseMixtureScalar.H"
 #include "turbulentTransportModel.H"
 #include "pimpleControl.H"
 #include "fvOptions.H"

applications/solvers/cfdemSolverMultiphaseScalar/createFields.H

@@ -88,7 +88,7 @@ surfaceScalarField phi
      linearInterpolate(U*voidfraction) & mesh.Sf()
  );
 
-multiphaseMixture mixture(U, phi, voidfraction);
+multiphaseMixtureScalar mixture(U, phi, voidfraction);
 
 // Need to store rho for ddt(rho, U)
 volScalarField rho

applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixtureScalar/Make/files

@@ -1,5 +1,5 @@
 phase/phase.C
 alphaContactAngle/alphaContactAngleFvPatchScalarField.C
-multiphaseMixture.C
+multiphaseMixtureScalar.C
 
 LIB = $(CFDEM_LIB_DIR)/libcfdemMultiphaseInterFoamScalar

applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixtureScalar/alphaContactAngle/alphaContactAngleFvPatchScalarField.H

@@ -26,7 +26,7 @@ Class
 
 Description
     Contact-angle boundary condition for multi-phase interface-capturing
-    simulations.  Used in conjuction with multiphaseMixture.
+    simulations.  Used in conjuction with multiphaseMixtureScalar.
 
 SourceFiles
     alphaContactAngleFvPatchScalarField.C
@@ -37,7 +37,7 @@ SourceFiles
 #define alphaContactAngleFvPatchScalarField_H
 
 #include "zeroGradientFvPatchFields.H"
-#include "multiphaseMixture.H"
+#include "multiphaseMixtureScalar.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -117,8 +117,8 @@ public:
     typedef HashTable
     <
         interfaceThetaProps,
-        multiphaseMixture::interfacePair,
-        multiphaseMixture::interfacePair::hash
+        multiphaseMixtureScalar::interfacePair,
+        multiphaseMixtureScalar::interfacePair::hash
     > thetaPropsTable;
 
 

applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixtureScalar/multiphaseMixtureScalar.C

@@ -18,7 +18,7 @@ License
 
 \*---------------------------------------------------------------------------*/
 
-#include "multiphaseMixture.H"
+#include "multiphaseMixtureScalar.H"
 #include "alphaContactAngleFvPatchScalarField.H"
 #include "Time.H"
 #include "subCycle.H"
@@ -31,13 +31,13 @@ License
 
 // * * * * * * * * * * * * * * * Static Member Data  * * * * * * * * * * * * //
 
-const Foam::scalar Foam::multiphaseMixture::convertToRad =
+const Foam::scalar Foam::multiphaseMixtureScalar::convertToRad =
     Foam::constant::mathematical::pi/180.0;
 
 
 // * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * //
 
-void Foam::multiphaseMixture::calcAlphas()
+void Foam::multiphaseMixtureScalar::calcAlphas()
 {
     scalar level = 0.0;
     alphas_ == 0.0;
@@ -51,7 +51,7 @@ void Foam::multiphaseMixture::calcAlphas()
 
 
 Foam::tmp<Foam::volScalarField>
-Foam::multiphaseMixture::calcNu() const
+Foam::multiphaseMixtureScalar::calcNu() const
 {
     PtrDictionary<phase>::const_iterator iter = phases_.begin();
 
@@ -74,7 +74,7 @@ Foam::multiphaseMixture::calcNu() const
 }
 
 Foam::tmp<Foam::surfaceScalarField>
-Foam::multiphaseMixture::calcStf() const
+Foam::multiphaseMixtureScalar::calcStf() const
 {
     tmp<surfaceScalarField> tstf
     (
@@ -134,7 +134,7 @@ Foam::multiphaseMixture::calcStf() const
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
-Foam::multiphaseMixture::multiphaseMixture
+Foam::multiphaseMixtureScalar::multiphaseMixtureScalar
 (
     const volVectorField& U,
     const surfaceScalarField& phi,
@@ -230,7 +230,7 @@ Foam::multiphaseMixture::multiphaseMixture
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * * //
 
 Foam::tmp<Foam::volScalarField>
-Foam::multiphaseMixture::rho() const
+Foam::multiphaseMixtureScalar::rho() const
 {
     PtrDictionary<phase>::const_iterator iter = phases_.begin();
 
@@ -247,7 +247,7 @@ Foam::multiphaseMixture::rho() const
 
 
 Foam::tmp<Foam::scalarField>
-Foam::multiphaseMixture::rho(const label patchi) const
+Foam::multiphaseMixtureScalar::rho(const label patchi) const
 {
     PtrDictionary<phase>::const_iterator iter = phases_.begin();
 
@@ -264,9 +264,9 @@ Foam::multiphaseMixture::rho(const label patchi) const
 
 
 Foam::tmp<Foam::volScalarField>
-Foam::multiphaseMixture::mu() const
+Foam::multiphaseMixtureScalar::mu() const
 {
-    Info << "In multiphasemixture mu()" << endl;
+    Info << "In multiphaseMixtureScalar mu()" << endl;
     return rho()*nu();
 //    PtrDictionary<phase>::const_iterator iter = phases_.begin();
 
@@ -283,7 +283,7 @@ Foam::multiphaseMixture::mu() const
 
 
 Foam::tmp<Foam::scalarField>
-Foam::multiphaseMixture::mu(const label patchi) const
+Foam::multiphaseMixtureScalar::mu(const label patchi) const
 {
     PtrDictionary<phase>::const_iterator iter = phases_.begin();
 
@@ -306,7 +306,7 @@ Foam::multiphaseMixture::mu(const label patchi) const
 
 
 Foam::tmp<Foam::surfaceScalarField>
-Foam::multiphaseMixture::muf() const
+Foam::multiphaseMixtureScalar::muf() const
 {
 
     return nuf()*fvc::interpolate(rho());
@@ -327,13 +327,13 @@ Foam::multiphaseMixture::muf() const
 
 
 Foam::tmp<Foam::volScalarField>
-Foam::multiphaseMixture::nu() const
+Foam::multiphaseMixtureScalar::nu() const
 {
     return nu_;
 }
 
 Foam::tmp<Foam::scalarField>
-Foam::multiphaseMixture::nu(const label patchi) const
+Foam::multiphaseMixtureScalar::nu(const label patchi) const
 {
     //return nu_.boundaryField()[patchi];
     PtrDictionary<phase>::const_iterator iter = phases_.begin();
@@ -355,7 +355,7 @@ Foam::multiphaseMixture::nu(const label patchi) const
 
 
 Foam::tmp<Foam::surfaceScalarField>
-Foam::multiphaseMixture::nuf() const
+Foam::multiphaseMixtureScalar::nuf() const
 {
     //return muf()/fvc::interpolate(rho());
     PtrDictionary<phase>::const_iterator iter = phases_.begin();
@@ -374,7 +374,7 @@ Foam::multiphaseMixture::nuf() const
 }
 
 Foam::tmp<Foam::volScalarField>
-Foam::multiphaseMixture::Cp() const
+Foam::multiphaseMixtureScalar::Cp() const
 {
     PtrDictionary<phase>::const_iterator iter = phases_.begin();
 
@@ -395,7 +395,7 @@ Foam::multiphaseMixture::Cp() const
 }
 
 Foam::tmp<Foam::volScalarField>
-Foam::multiphaseMixture::kf() const
+Foam::multiphaseMixtureScalar::kf() const
 {
     PtrDictionary<phase>::const_iterator iter = phases_.begin();
 
@@ -417,7 +417,7 @@ Foam::multiphaseMixture::kf() const
 }
 
 Foam::tmp<Foam::volScalarField>
-Foam::multiphaseMixture::D() const
+Foam::multiphaseMixtureScalar::D() const
 {
     PtrDictionary<phase>::const_iterator iter = phases_.begin();
 
@@ -439,7 +439,7 @@ Foam::multiphaseMixture::D() const
 }
 
 Foam::tmp<Foam::volScalarField>
-Foam::multiphaseMixture::Cs() const
+Foam::multiphaseMixtureScalar::Cs() const
 {
     PtrDictionary<phase>::const_iterator iter = phases_.begin();
 
@@ -456,7 +456,7 @@ Foam::multiphaseMixture::Cs() const
 }
 
 Foam::tmp<Foam::surfaceScalarField>
-Foam::multiphaseMixture::diffusionCorrection() const
+Foam::multiphaseMixtureScalar::diffusionCorrection() const
 {
 
     surfaceScalarField numerator
@@ -517,7 +517,7 @@ Foam::multiphaseMixture::diffusionCorrection() const
     return correction;
 }
 
-void Foam::multiphaseMixture::solve()
+void Foam::multiphaseMixtureScalar::solve()
 {
     correct();
 
@@ -570,7 +570,7 @@ void Foam::multiphaseMixture::solve()
 }
 
 
-void Foam::multiphaseMixture::correct()
+void Foam::multiphaseMixtureScalar::correct()
 {
     forAllIter(PtrDictionary<phase>, phases_, iter)
     {
@@ -579,7 +579,7 @@ void Foam::multiphaseMixture::correct()
 }
 
 
-Foam::tmp<Foam::surfaceVectorField> Foam::multiphaseMixture::nHatfv
+Foam::tmp<Foam::surfaceVectorField> Foam::multiphaseMixtureScalar::nHatfv
 (
     const volScalarField& alpha1,
     const volScalarField& alpha2
@@ -605,7 +605,7 @@ Foam::tmp<Foam::surfaceVectorField> Foam::multiphaseMixture::nHatfv
 }
 
 
-Foam::tmp<Foam::surfaceScalarField> Foam::multiphaseMixture::nHatf
+Foam::tmp<Foam::surfaceScalarField> Foam::multiphaseMixtureScalar::nHatf
 (
     const volScalarField& alpha1,
     const volScalarField& alpha2
@@ -622,7 +622,7 @@ Foam::tmp<Foam::surfaceScalarField> Foam::multiphaseMixture::nHatf
 // The dynamic contact angle is calculated from the component of the
 // velocity on the direction of the interface, parallel to the wall.
 
-void Foam::multiphaseMixture::correctContactAngle
+void Foam::multiphaseMixtureScalar::correctContactAngle
 (
     const phase& alpha1,
     const phase& alpha2,
@@ -726,7 +726,7 @@ void Foam::multiphaseMixture::correctContactAngle
 }
 
 
-Foam::tmp<Foam::volScalarField> Foam::multiphaseMixture::K
+Foam::tmp<Foam::volScalarField> Foam::multiphaseMixtureScalar::K
 (
     const phase& alpha1,
     const phase& alpha2
@@ -742,7 +742,7 @@ Foam::tmp<Foam::volScalarField> Foam::multiphaseMixture::K
 
 
 Foam::tmp<Foam::volScalarField>
-Foam::multiphaseMixture::nearInterface() const
+Foam::multiphaseMixtureScalar::nearInterface() const
 {
     tmp<volScalarField> tnearInt
     (
@@ -768,7 +768,7 @@ Foam::multiphaseMixture::nearInterface() const
 }
 
 
-void Foam::multiphaseMixture::solveAlphas
+void Foam::multiphaseMixtureScalar::solveAlphas
 (
     const scalar cAlpha
 )
@@ -901,7 +901,7 @@ void Foam::multiphaseMixture::solveAlphas
 }
 
 
-bool Foam::multiphaseMixture::read()
+bool Foam::multiphaseMixtureScalar::read()
 {
     if (transportModel::read())
     {

applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixtureScalar/multiphaseMixtureScalar.H

@@ -17,10 +17,10 @@ License
     Copyright (C) 2018- Mathias Vångö, JKU Linz, Austria
 
 Class
-    multiphaseMixture
+    multiphaseMixtureScalar
 
 Description
-    This class is based on the OpenFOAM(R) Foam::multiphaseMixture class,
+    This class is based on the OpenFOAM(R) Foam::multiphaseMixtureScalar class,
     which is an incompressible multi-phase mixture with built in solution
     for the phase fractions with interface compression for interface-capturing.
     It has been extended to include the void fraction in the volume fraction
@@ -33,11 +33,11 @@ Description
     between each phase-pair.
 
 SourceFiles
-    multiphaseMixture.C
+    multiphaseMixtureScalar.C
 \*---------------------------------------------------------------------------*/
 
-#ifndef multiphaseMixture_H
-#define multiphaseMixture_H
+#ifndef multiphaseMixtureScalar_H
+#define multiphaseMixtureScalar_H
 
 #include "incompressible/transportModel/transportModel.H"
 #include "IOdictionary.H"
@@ -52,10 +52,10 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                      Class multiphaseMixture Declaration
+                      Class multiphaseMixtureScalar Declaration
 \*---------------------------------------------------------------------------*/
 
-class multiphaseMixture
+class multiphaseMixtureScalar
 :
     public IOdictionary,
     public transportModel
@@ -191,7 +191,7 @@ public:
     // Constructors
 
         //- Construct from components
-        multiphaseMixture
+        multiphaseMixtureScalar
         (
             const volVectorField& U,
             const surfaceScalarField& phi,
@@ -200,7 +200,7 @@ public:
 
 
     //- Destructor
-    virtual ~multiphaseMixture()
+    virtual ~multiphaseMixtureScalar()
     {}
 
 

applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixtureScalar/phase/phase.H

@@ -26,7 +26,7 @@ Class
 
 Description
     Single incompressible phase derived from the phase-fraction.
-    Used as part of the multiPhaseMixture for interface-capturing multi-phase
+    Used as part of the multiphaseMixtureScalar for interface-capturing multi-phase
     simulations.
 
 SourceFiles

applications/solvers/cfdemSolverPimple/Make/files

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

applications/solvers/cfdemSolverPimple/Make/options

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

applications/solvers/cfdemSolverPimple/UEqn.H

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

applications/solvers/cfdemSolverPimple/cfdemSolverPimple.C

@@ -0,0 +1,147 @@
+/*---------------------------------------------------------------------------*\
+    Open Source CFD-DEM coupling
+
+    Copyright (C) 2023  Behrad Esgandari, JKU Linz, Austria
+-------------------------------------------------------------------------------
+License
+
+    This program is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this program.  If not, see <https://www.gnu.org/licenses/>.
+
+Application
+    cfdemSolverPimple
+
+Description
+    Transient solver for incompressible flow.
+    Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
+    The code is an evolution of the solver pimpleFoam in OpenFOAM(R) 6.0,
+    where additional functionality for CFD-DEM coupling is added.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "singlePhaseTransportModel.H"
+#include "turbulentTransportModel.H"
+#include "pimpleControl.H"
+#include "fvOptions.H"
+
+#include "cfdemCloud.H"
+#include "implicitCouple.H"
+#include "clockModel.H"
+#include "smoothingModel.H"
+#include "forceModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createControl.H"
+    #include "createFields.H"
+    #include "createFvOptions.H"
+    #include "initContinuityErrs.H"
+
+    // create cfdemCloud
+    #include "readGravitationalAcceleration.H"
+    cfdemCloud particleCloud(mesh);
+    #include "checkModelType.H"
+
+    // switch for periodic box simulations
+    Switch periodicBoxSwitch
+    (
+        pimple.dict().lookupOrDefault<Switch>("periodicBoxSwitch", false)
+    );
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+    Info<< "\nStarting time loop\n" << endl;
+    while (runTime.loop())
+    {
+        particleCloud.clockM().start(1,"Global");
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        #include "CourantNo.H"
+
+        // do particle stuff
+        particleCloud.clockM().start(2,"Coupling");
+        bool hasEvolved = particleCloud.evolve(voidfraction,Us,U);
+
+        if(hasEvolved)
+        {
+            particleCloud.smoothingM().smoothen(particleCloud.forceM(0).impParticleForces());
+        }
+
+        Info << "update Ksl.internalField()" << endl;
+        Ksl = particleCloud.momCoupleM(0).impMomSource();
+        Ksl.correctBoundaryConditions();
+
+        //Force Checks
+        vector fTotal(0,0,0);
+        vector fImpTotal = sum(mesh.V()*Ksl.internalField()*(Us.internalField()-U.internalField())).value();
+        reduce(fImpTotal, sumOp<vector>());
+        Info << "TotalForceExp: " << fTotal << endl;
+        Info << "TotalForceImp: " << fImpTotal << endl;
+
+        #include "solverDebugInfo.H"
+        particleCloud.clockM().stop("Coupling");
+
+        particleCloud.clockM().start(26,"Flow");
+
+        if(particleCloud.solveFlow())
+        {
+            // Pressure-velocity PIMPLE corrector
+            while (pimple.loop())
+            {
+                // Momentum predictor
+                #include "UEqn.H"
+
+                // --- Inner PIMPLE loop
+
+                while (pimple.correct())
+                {
+                    #include "pEqn.H"
+                }
+            }
+
+            laminarTransport.correct();
+            turbulence->correct();
+        }
+        else
+        {
+            Info << "skipping flow solution." << endl;
+        }
+
+        if (periodicBoxSwitch)
+        {
+           #include "periodicBoxProperties.H"
+        }
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+
+        particleCloud.clockM().stop("Flow");
+        particleCloud.clockM().stop("Global");
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/cfdemSolverPimple/createFields.H

@@ -0,0 +1,169 @@
+//===============================
+// Fluid Fields
+//===============================
+
+   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
+    );
+
+//===============================
+// particle interaction modelling
+//===============================
+
+    Info<< "\nReading momentum exchange field Ksl\n" << endl;
+    volScalarField Ksl
+    (
+        IOobject
+        (
+            "Ksl",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    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<< "\nReading voidfraction field voidfraction = (Vgas/Vparticle)\n" << endl;
+    volScalarField voidfraction
+    (
+        IOobject
+        (
+            "voidfraction",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< "\nCreating dummy density field rho\n" << endl;
+    volScalarField rho
+    (
+        IOobject
+        (
+            "rho",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< "Reading particle velocity field Us\n" << endl;
+    volVectorField Us
+    (
+        IOobject
+        (
+            "Us",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+//===============================
+
+#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
+
+    Info<< "Generating interstitial(!) flux field phiByVoidfraction\n" << endl;
+    surfaceScalarField phiByVoidfraction
+    (
+        IOobject
+        (
+            "phiByVoidfraction",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        linearInterpolate(U) & mesh.Sf()
+    );
+
+//Periodic box
+    volScalarField unity
+    (
+        IOobject
+        (
+            "unity",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("unity", dimless, 1.0)
+    );
+
+    label pRefCell = 0;
+    scalar pRefValue = 0.0;
+    setRefCell(p, mesh.solutionDict().subDict("PIMPLE"), pRefCell, pRefValue);
+
+
+    singlePhaseTransportModel laminarTransport(U, phi);
+
+    autoPtr<incompressible::turbulenceModel> turbulence
+    (
+        incompressible::turbulenceModel::New(U, phi, laminarTransport)
+    );
+
+#include "createMRF.H"

applications/solvers/cfdemSolverPimple/pEqn.H

@@ -0,0 +1,59 @@
+volScalarField rAUvoidfraction("(voidfraction2|A(U))",rAU*voidfraction);
+
+surfaceScalarField rAUfvoidfraction("(voidfraction2|A(U)F)", fvc::interpolate(rAUvoidfraction));
+
+volVectorField HbyA("HbyA", U);
+
+HbyA = rAU*UEqn.H();
+
+phi = voidfractionf*phiByVoidfraction;
+
+surfaceScalarField phiHbyA
+(
+    "phiHbyA",
+    (
+        (fvc::interpolate(HbyA) & mesh.Sf() )
+        + phicForces                               //explicit contribution
+        + rAUfvoidfraction*fvc::ddtCorr(U, phiByVoidfraction) //correction
+    )
+);
+
+if (modelType=="A")
+    rAUvoidfraction = volScalarField("(voidfraction2|A(U))",rAU*voidfraction*voidfraction);
+
+// Update the fixedFluxPressure BCs to ensure flux consistency
+if (modelType=="A")
+{
+    volScalarField rUsed = rAU*voidfraction;
+    constrainPressure(p, U, phiHbyA, rUsed,MRF);
+}
+else constrainPressure(p, U, phiHbyA, rAU,MRF);
+
+while (pimple.correctNonOrthogonal())
+{
+    // Pressure corrector
+    fvScalarMatrix pEqn
+    (
+        fvm::laplacian(rAUvoidfraction, p) == fvc::div(voidfractionf*phiHbyA) + particleCloud.ddtVoidfraction()
+    );
+    pEqn.setReference(pRefCell, pRefValue);
+
+    pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
+
+    if (pimple.finalNonOrthogonalIter())
+    {
+        phiByVoidfraction = phiHbyA - pEqn.flux()/voidfractionf;
+        phi = voidfractionf*phiByVoidfraction;
+
+        #include "continuityErrorPhiPU.H"
+
+        // Explicitly relax pressure for momentum corrector
+        p.relax();
+
+        U = fvc::reconstruct(phiHbyA)
+          - rAU*fvc::reconstruct(pEqn.flux()/voidfractionf/rAUf);
+
+        U.correctBoundaryConditions();
+        fvOptions.correct(U);
+    }
+}

applications/solvers/cfdemSolverPimple/periodicBoxProperties.H

@@ -0,0 +1,65 @@
+dimensionedVector gN = mesh.lookupObject<uniformDimensionedVectorField> ("g");
+dimensionedScalar volume = fvc::domainIntegrate(unity);
+
+Info<< "particle_ENSTROPHY: "
+    << (
+            fvc::domainIntegrate( 0.5*magSqr(fvc::curl(Us)))
+            /volume
+       ).value()
+    << endl;
+
+Info<< "air_ENSTROPHY: "
+    << (
+           fvc::domainIntegrate( 0.5*magSqr(fvc::curl(U)))
+           /volume
+       ).value()
+    << endl;
+
+Info<< "slip_velocity: "
+    << - ((
+        fvc::domainIntegrate(voidfraction*(U&gN)).value()
+           /fvc::domainIntegrate(voidfraction*mag(gN)).value()
+         )
+       - (
+        fvc::domainIntegrate((1.0-voidfraction)*(Us&gN)).value()
+           /fvc::domainIntegrate((1.0-voidfraction)*mag(gN)).value()
+         ))
+   << endl;
+
+dimensionedVector alpha1Us = fvc::domainIntegrate((1.0-voidfraction)*(Us))/fvc::domainIntegrate((1.0-voidfraction));
+dimensionedVector alpha2U = fvc::domainIntegrate(voidfraction*(U))/fvc::domainIntegrate(voidfraction);
+dimensionedScalar alpha1M  = fvc::domainIntegrate((1.0-voidfraction))/volume;
+dimensionedScalar alpha2M  = scalar(1.0) - alpha1M;
+
+Info<< "TKE gas: "
+    << 0.5
+      *(
+          fvc::domainIntegrate(voidfraction*(U&U)).value()
+         /fvc::domainIntegrate(voidfraction).value()
+       )
+     - 0.5
+      *(
+          alpha2U.value()
+         &alpha2U.value()
+       )
+    << endl;
+
+Info<< "TKE solid: "
+    << 0.5
+      *(
+          fvc::domainIntegrate((1.0-voidfraction)*(Us&Us)).value()
+         /fvc::domainIntegrate(1.0-voidfraction).value()
+       )
+     - 0.5
+      *(
+          alpha1Us.value()
+         &alpha1Us.value()
+       )
+    << endl;
+
+Info<< "PhiP2: "
+    << fvc::domainIntegrate((1.0-voidfraction)*(1.0-voidfraction)).value()
+       /fvc::domainIntegrate(unity).value()
+       - alpha1M.value()*alpha1M.value()
+    << endl;
+

applications/solvers/cfdemSolverPiso/Make/options

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

applications/solvers/cfdemSolverPisoFreeStreaming/Make/options

@@ -1,7 +1,6 @@
 include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
 
 EXE_INC = \
-    -I$(CFDEM_OFVERSION_DIR) \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
@@ -12,7 +11,8 @@ EXE_INC = \
     -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/lnInclude \
     -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/cfdTools \
     -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/derived/cfdemCloudRec \
-    -I$(LIB_SRC)/sampling/lnInclude
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -Wno-deprecated-copy
 
 EXE_LIBS = \
     -L$(CFDEM_LIB_DIR)\

applications/solvers/cfdemSolverPisoMS/Make/options

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

applications/solvers/cfdemSolverPisoScalar/Make/options

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

applications/solvers/cfdemSolverRhoPimple/EEqn.H

@@ -51,8 +51,4 @@
     thermo.correct();
 
     Info<< "T max/min/ave : " << max(T).value() << " " << min(T).value() << " " << average(T).value() << endl;
-
-    particleCloud.clockM().start(31,"energySolve");
-    particleCloud.solve();
-    particleCloud.clockM().stop("energySolve");
 }

applications/solvers/cfdemSolverRhoPimple/Make/options

@@ -6,7 +6,6 @@ PFLAGS+= -Dcompre
 
 EXE_INC = \
      $(PFLAGS) \
-    -I$(CFDEM_OFVERSION_DIR) \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \

applications/solvers/cfdemSolverRhoPimpleChem/EEqn.H

@@ -50,13 +50,9 @@ fvScalarMatrix EEqn
 
     thermo.correct();
 
-    Info << "Qsource :" << max(Qsource).value() << " " << min(Qsource).value() << endl;
-    Info << "QCoeff :" << max(QCoeff).value() << " " << min(QCoeff).value() << endl;
-    Info << "Cpv :" << max(Cpv).value() << " " << min(Cpv).value() << endl;
-    Info<<  "T max/min : " << max(T).value() << " " << min(T).value() << endl;
-    Info << "he max/min : " << max(he).value() << "  " << min(he).value() << endl;
-
-    particleCloud.clockM().start(31,"energySolve");
-    particleCloud.solve();
-    particleCloud.clockM().stop("energySolve");
+    Info<< "Qsource :"      << max(Qsource).value() << " " << min(Qsource).value()  << endl;
+    Info<< "QCoeff :"       << max(QCoeff).value()  << " " << min(QCoeff).value()   << endl;
+    Info<< "Cpv :"          << max(Cpv).value() << " " << min(Cpv).value()  << endl;
+    Info<< "T max/min/ave : " << max(T).value() << " " << min(T).value() << " " << average(T).value() << endl;
+    Info<< "he max/min : "  << max(he).value()  << " " << min(he).value()   << endl;
 }

applications/solvers/cfdemSolverRhoPimpleChem/Make/options

@@ -6,7 +6,6 @@ PFLAGS+= -Dcompre
 
 EXE_INC = \
      $(PFLAGS) \
-    -I$(CFDEM_OFVERSION_DIR) \
     -I$(LIB_SRC)/finiteVolume/cfdTools \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \

applications/solvers/cfdemSolverRhoPimpleChem/YEqn.H

@@ -60,13 +60,13 @@ tmp<fv::convectionScheme<scalar> > mvConvection
 
     if (propagateInertSpecie)
     {
-        if (inertIndex!=-1) Yt /= (1-Y[inertIndex] + VSMALL);
+        if (inertIndex!=-1) Yt /= (1-Y[inertIndex] + ROOTVSMALL);
         forAll(Y,i)
         {
             if (i!=inertIndex)
             {
                 volScalarField& Yi = Y[i];
-                Yi = Yi/(Yt+VSMALL);
+                Yi = Yi/(Yt+ROOTVSMALL);
             }
         }
     }

applications/solvers/cfdemSolverRhoPimpleChem/cfdemSolverRhoPimpleChem.C

@@ -23,7 +23,7 @@ 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) 2.3,
+    The code is an evolution of the solver rhoPimpleFoam in OpenFOAM(R) 4.x,
     where additional functionality for CFD-DEM coupling is added.
 \*---------------------------------------------------------------------------*/
 
@@ -38,6 +38,9 @@ Description
 #endif
 #include "bound.H"
 #include "pimpleControl.H"
+#if OPENFOAM_VERSION_MAJOR >= 5
+#include "pressureControl.H"
+#endif
 #include "fvOptions.H"
 #include "localEulerDdtScheme.H"
 #include "fvcSmooth.H"
@@ -78,10 +81,15 @@ int main(int argc, char *argv[])
 
     turbulence->validate();
 
+    #include "compressibleCourantNo.H"
+    #include "setInitialDeltaT.H"
+
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
     Info<< "\nStarting time loop\n" << endl;
+    bool firstStep = true;
 
+    // monitor mass variables
     scalar m(0.0);
     scalar m0(0.0);
     label counter(0);
@@ -125,17 +133,28 @@ int main(int argc, char *argv[])
 
 #if OPENFOAM_VERSION_MAJOR < 6
         if (pimple.nCorrPIMPLE() <= 1)
-#else
-        if (pimple.nCorrPimple() <= 1)
-#endif
         {
             #include "rhoEqn.H"
         }
+        rhoeps = rho*voidfraction;
+#endif
 
-        rhoeps = rho * voidfraction;
         // --- Pressure-velocity PIMPLE corrector loop
         while (pimple.loop())
         {
+#if OPENFOAM_VERSION_MAJOR >= 6
+            if (pimple.firstIter())
+            {
+                #include "rhoEqn.H"
+                if (firstStep)
+                {
+                    rhoeps.oldTime() = rho.oldTime()*voidfraction.oldTime();
+                    firstStep = false;
+                }
+                rhoeps = rho*voidfraction;
+            }
+#endif
+
             #include "UEqn.H"
             #include "YEqn.H"
             #include "EEqn.H"

applications/solvers/cfdemSolverRhoPimpleChem/createFields.H

@@ -94,8 +94,48 @@
         mesh
     );
 
-    volScalarField rhoeps ("rhoeps", rho*voidfraction);
+    volScalarField rhoeps("rhoeps", rho*voidfraction);
+    rhoeps.oldTime(); // switch on saving old time
 
+    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()
+    );
+
+#if OPENFOAM_VERSION_MAJOR < 5
+    dimensionedScalar rhoMax
+    (
+        dimensionedScalar::lookupOrDefault
+        (
+            "rhoMax",
+            pimple.dict(),
+            dimDensity,
+            GREAT
+        )
+    );
+
+    dimensionedScalar rhoMin
+    (
+        dimensionedScalar::lookupOrDefault
+        (
+            "rhoMin",
+            pimple.dict(),
+            dimDensity,
+            0
+        )
+    );
+#else
+    pressureControl pressureControl(p, rho, pimple.dict(), false);
+#endif
 
     Info<< "\nCreating fluid-particle heat flux field\n" << endl;
     volScalarField Qsource
@@ -158,41 +198,6 @@
         dimensionedScalar("zero", dimensionSet(0,2,-2,-1,0,0,0), 0.0)
     );
 
-    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
@@ -206,6 +211,8 @@
         )
     );
 
+    mesh.setFluxRequired(p.name());
+
 #if OPENFOAM_VERSION_MAJOR >= 6
     Info<< "Creating combustion model\n" << endl;
     autoPtr<CombustionModel<rhoReactionThermo>> combustion
@@ -303,18 +310,4 @@
         mesh,
         dimensionedScalar("zero",dimensionSet(0, -3, 0, 0, 1),0)
     );
-
-    volScalarField dSauter
-    (
-        IOobject
-        (
-            "dSauter",
-            runTime.timeName(),
-            mesh,
-            IOobject::READ_IF_PRESENT,
-            IOobject::AUTO_WRITE
-        ),
-        mesh,
-        dimensionedScalar("zero",dimensionSet(0, 1, 0, 0, 0,0,0),0)
-    );
 //===============================

applications/solvers/cfdemSolverRhoPimpleChem/pEqn.H

@@ -1,14 +1,19 @@
 rho = thermo.rho();
+#if OPENFOAM_VERSION_MAJOR < 5
 rho = max(rho, rhoMin);
 rho = min(rho, rhoMax);
 rho.relax();
+rhoeps = rho*voidfraction;
+#else
+rhoeps = rho*voidfraction;
+
+// Thermodynamic density needs to be updated by psi*d(p) after the
+// pressure solution
+const volScalarField psip0(psi*p);
+#endif
 
 volScalarField rAU(1.0/UEqn.A());
-surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rhoeps*rAU));
-if (modelType=="A")
-{
-    rhorAUf *= fvc::interpolate(voidfraction);
-}
+surfaceScalarField rhorAUf("rhorAUf", (modelType=="A")?fvc::interpolate(voidfraction*rhoeps*rAU):fvc::interpolate(rhoeps*rAU));
 volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
 
 surfaceScalarField phiUs("phiUs", fvc::interpolate(rhoeps*rAU*Ksl*Us)& mesh.Sf());
@@ -18,21 +23,19 @@ if (pimple.nCorrPISO() <= 1)
     tUEqn.clear();
 }
 
+surfaceScalarField phiHbyA
+(
+    "phiHbyA",
+        fvc::interpolate(rhoeps)*fvc::flux(HbyA)
+      + rhorAUf*fvc::ddtCorr(rhoeps, U, phi)
+);
+
 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;
 
@@ -41,18 +44,34 @@ else
 
     volScalarField SmbyP(particleCloud.chemistryM(0).Sm() / p);
 
+#if OPENFOAM_VERSION_MAJOR >= 5
+    fvScalarMatrix pDDtEqn
+    (
+        fvc::ddt(rhoeps)
+      + psi*voidfraction*correction(fvm::ddt(p))
+      + fvc::div(phi)
+     ==
+        fvm::Sp(SmbyP, p)
+      + fvOptions(psi, p, rho.name())
+    );
+#endif
+
     while (pimple.correctNonOrthogonal())
     {
         // Pressure corrector
+#if OPENFOAM_VERSION_MAJOR < 5
         fvScalarMatrix pEqn
         (
-          fvm::ddt(voidfraction, psi, p)
+            fvm::ddt(voidfraction, psi, p)
           + fvc::div(phi)
           - fvm::laplacian(rhorAUf, p)
           ==
             fvm::Sp(SmbyP, p)
           + fvOptions(psi, p, rho.name())
         );
+#else
+        fvScalarMatrix pEqn(pDDtEqn - fvm::laplacian(rhorAUf, p));
+#endif
 
         pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
 
@@ -63,21 +82,18 @@ else
     }
 }
 
+// Thermodynamic density update
+#if OPENFOAM_VERSION_MAJOR >= 5
+thermo.correctRho(psi*p - psip0);
+#endif
+
 #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;
-
-rhoeps = rho * voidfraction;
+Info<< "p max/min/ave : " << max(p).value()
+    << " " << min(p).value() << " " << average(p).value() << endl;
 
 if (modelType=="A")
 {
@@ -91,6 +107,24 @@ U.correctBoundaryConditions();
 fvOptions.correct(U);
 K = 0.5*magSqr(U);
 
+// Recalculate density from the relaxed pressure
+#if OPENFOAM_VERSION_MAJOR >= 5
+if (pressureControl.limit(p))
+{
+    p.correctBoundaryConditions();
+}
+rho = thermo.rho();
+#else
+rho = thermo.rho();
+rho = max(rho, rhoMin);
+rho = min(rho, rhoMax);
+rho.relax();
+#endif
+rhoeps = rho*voidfraction;
+
+Info<< "rho max/min/ave : " << max(rho).value()
+    << " " << min(rho).value() << " " << average(rho).value() << endl;
+
 if (thermo.dpdt())
 {
     dpdt = fvc::ddt(voidfraction,p);

applications/solvers/rStatAnalysis/Make/options

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

applications/solvers/rcfdemSolverBase/Make/options

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

applications/solvers/rcfdemSolverCoupledHeattransfer/Make/options

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

applications/solvers/rcfdemSolverCoupledHeattransfer/TEqn.H

@@ -1,4 +1,4 @@
-    volScalarField rhoeps = rhoRec*voidfractionRec;
+    rhoeps = rhoRec*voidfractionRec;
 
     particleCloud.energyContributions(Qsource);
 
@@ -23,10 +23,17 @@
         fvOptions(rhoeps, T)    // no fvOptions support yet
     );
 
-    fvOptions.constrain(TEqn); // no fvOptions support yet
+ //   fvOptions.constrain(TEqn); // no fvOptions support yet
+
+    TEqn.relax();
 
     TEqn.solve();
 
+    T = max(T, TMin);
+    T = min(T, TMax);
+
+    Info<< "T max/min/ave : " << max(T).value() << " " << min(T).value() << " " << average(T).value() << endl;
+
     particleCloud.clockM().start(31,"postFlow");
     counter++;
 

applications/solvers/rcfdemSolverCoupledHeattransfer/createFields.H

@@ -73,6 +73,9 @@
         dimensionedVector("zero", dimensionSet(0, 1, -1, 0, 0), vector::zero)
     );
 
+    volScalarField rhoeps("rhoeps", rhoRec*voidfractionRec);
+    rhoeps.oldTime(); // switch on saving old time
+
     // heat transfer fields
     Info << "\nCreating heat transfer fields.\n" << endl;
 
@@ -228,3 +231,25 @@
         )
     );
     weightDict.add("weights",scalarList(1,1.0));
+
+    dimensionedScalar TMax
+    (
+        dimensionedScalar::lookupOrDefault
+        (
+            "TMax",
+            transportProps,
+            dimTemperature,
+            GREAT
+        )
+    );
+
+    dimensionedScalar TMin
+    (
+        dimensionedScalar::lookupOrDefault
+        (
+            "TMin",
+            transportProps,
+            dimTemperature,
+            0.0
+        )
+    );

applications/solvers/rcfdemSolverCoupledHeattransfer/rcfdemSolverCoupledHeattransfer.C

@@ -98,7 +98,7 @@ int main(int argc, char *argv[])
 
         particleCloud.clockM().start(26,"Flow");
         #include "updateRho.H"
-        #include "TEqImp.H"
+        #include "TEqn.H"
         particleCloud.clockM().stop("Flow");
 
         stepCounter++;

applications/solvers/rcfdemSolverCoupledHeattransfer/updateRho.H

@@ -1 +1,2 @@
-rhoRec = pRec / (T * R);
+dimensionedScalar Tave = T.weightedAverage(voidfractionRec);
+rhoRec = pRec / (Tave * R);

applications/solvers/rcfdemSolverForcedTracers/Make/options

@@ -1,7 +1,6 @@
 include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
 
 EXE_INC = \
-    -I$(CFDEM_OFVERSION_DIR) \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
@@ -12,6 +11,7 @@ EXE_INC = \
     -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/cfdTools \
     -I$(CFDEM_SRC_DIR)/recurrence/lnInclude \
     -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/derived/cfdemCloudRec \
+    -Wno-deprecated-copy
 
 EXE_LIBS = \
     -L$(CFDEM_LIB_DIR)\
@@ -24,4 +24,4 @@ EXE_LIBS = \
     -lfvOptions \
     -l$(CFDEM_LIB_NAME) \
      $(CFDEM_ADD_LIB_PATHS) \
-$(CFDEM_ADD_LIBS)
+     $(CFDEM_ADD_LIBS)

applications/solvers/rcfdemSolverRhoSteadyPimple/EEqn.H

@@ -60,9 +60,4 @@
     thermo.correct();
 
     Info<< "T max/min : " << max(T).value() << " " << min(T).value() << endl;
-
-
-    particleCloud.clockM().start(31,"energySolve");
-    particleCloud.solve();
-    particleCloud.clockM().stop("energySolve");
 }

applications/solvers/rcfdemSolverRhoSteadyPimple/Make/options

@@ -6,7 +6,6 @@ PFLAGS+= -Dcompre
 
 EXE_INC = \
      $(PFLAGS) \
-    -I$(CFDEM_OFVERSION_DIR) \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \

applications/solvers/rcfdemSolverRhoSteadyPimpleChem/EEqn.H

@@ -60,8 +60,5 @@
 
     Info<< "T max/min : " << max(T).value() << " " << min(T).value() << endl;
 
-
-    particleCloud.clockM().start(31,"energySolve");
-    particleCloud.solve();
-    particleCloud.clockM().stop("energySolve");
+    QFluidCond = fvc::laplacian(voidfractionRec*thCond,T);
 }

applications/solvers/rcfdemSolverRhoSteadyPimpleChem/Make/options

@@ -6,7 +6,6 @@ PFLAGS+= -Dcompre
 
 EXE_INC = \
      $(PFLAGS) \
-    -I$(CFDEM_OFVERSION_DIR) \
     -I$(LIB_SRC)/transportModels/compressible/lnInclude \
     -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
@@ -27,6 +26,7 @@ EXE_INC = \
     -I$(LIB_SRC)/ODE/lnInclude \
     -I$(LIB_SRC)/combustionModels/lnInclude \
     -I$(CFDEM_SRC_DIR)/recurrence/lnInclude \
+    -Wno-deprecated-copy
 
 EXE_LIBS = \
     -L$(CFDEM_LIB_DIR)\

applications/solvers/rcfdemSolverRhoSteadyPimpleChem/YEqn.H

@@ -20,6 +20,7 @@ tmp<fv::convectionScheme<scalar> > mvConvection
 #endif
     label inertIndex = -1;
     volScalarField Yt(0.0*Y[0]);
+    Sm *= 0.0;
 
     forAll(Y, i)
     {
@@ -28,13 +29,17 @@ tmp<fv::convectionScheme<scalar> > mvConvection
         {
             volScalarField& Yi = Y[i];
 
+            volScalarField sourceField(particleCloud.chemistryM(0).Smi(i));
+            volScalarField Smi0(neg(sourceField)*sourceField/(Yi + Yismall));
+            volScalarField Smi1(pos0(sourceField)*sourceField);
             fvScalarMatrix YiEqn
             (
                 mvConvection->fvmDiv(phi, Yi)
               - fvm::laplacian(voidfractionRec*turbulence->muEff(), Yi)
               ==
                 combustion->R(Yi)
-              + particleCloud.chemistryM(0).Smi(i)*p/p.prevIter()
+                + fvm::Sp(Smi0,Yi)
+                + Smi1
               + fvOptions(rho, Yi)
             );
 
@@ -48,8 +53,11 @@ tmp<fv::convectionScheme<scalar> > mvConvection
 
             fvOptions.correct(Yi);
 
+            #include "monitorMassSinks.H"
             Yi.max(0.0);
             if (Y[i].name() != inertSpecie) Yt += Yi;
+            #include "monitorMassSources.H"
+            Sm += Smi0*Yi+Smi1;
         }
     }
 

applications/solvers/rcfdemSolverRhoSteadyPimpleChem/createFields.H

@@ -42,6 +42,17 @@ Info<< "Reading thermophysical properties\n" << endl;
 
     Info<< "inert will be bounded in [" << inertLowerBound << "," << inertUpperBound << "]" << endl;
 
+    #include "OFstream.H"
+    OFstream Hf("Hf");
+    Hf << "# species  Hf (J/kg)" << endl;
+    Info << "\nspecies-specific heat of formation (J/kg):" << endl;
+    forAll(composition.species(),i)
+    {
+        Info << composition.species()[i] << " " <<  composition.Hc(i) << endl;
+        Hf << composition.species()[i] << " " <<  composition.Hc(i) << endl;
+    }
+    Info << "\n" << endl;
+
     volScalarField& p = thermo.p();
 
     multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
@@ -122,6 +133,20 @@ Info<< "Reading thermophysical properties\n" << endl;
         dimensionedScalar("zero", dimensionSet(1,-1,-3,0,0,0,0), 0.0)
     );
 
+    volScalarField Sm
+    (
+        IOobject
+        (
+            "Sm",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+         ),
+        mesh,
+        dimensionedScalar("zero",dimMass/(dimVol*dimTime),0.0)
+    );
+
     Info<< "\nCreating fluid-particle heat flux field\n" << endl;
     volScalarField Qsource
     (

applications/solvers/rcfdemSolverRhoSteadyPimpleChem/monitorMass.H

@@ -1,9 +0,0 @@
-{
-    m=gSum(rhoeps*1.0*rhoeps.mesh().V());
-    if(counter==0) m0=m;
-    counter++;
-    Info << "\ncurrent gas mass = " << m << "\n" << endl;
-    Info << "\ncurrent added gas mass = " << m-m0 << "\n" << endl;
-
-    QFluidCond = fvc::laplacian(voidfractionRec*thCond,T);
-}

applications/solvers/rcfdemSolverRhoSteadyPimpleChem/monitorMassSinks.H

@@ -0,0 +1,12 @@
+{
+    scalar massSink = 0.0;
+    forAll(Yi,cellI)
+    {
+        if (Yi[cellI] <= 0.0)
+        {
+            massSink += rhoeps[cellI]*Yi[cellI]*Yi.mesh().V()[cellI];
+        }
+    }
+    reduce(massSink, sumOp<scalar>());
+    Info << Y[i].name() << ": mass sink = " << massSink << endl;
+}

applications/solvers/rcfdemSolverRhoSteadyPimpleChem/monitorMassSource.H

@@ -0,0 +1,9 @@
+{
+    scalar sourceStrength = 0.0;
+    forAll(p,cellI)
+    {
+        sourceStrength += (Sm0[cellI]*p[cellI]+Sm1[cellI])*p.mesh().V()[cellI];
+    }
+    reduce(sourceStrength, sumOp<scalar>());
+    Info << "total mass source strength = " << sourceStrength << endl;
+}

applications/solvers/rcfdemSolverRhoSteadyPimpleChem/monitorMassSources.H

@@ -0,0 +1,9 @@
+{
+    scalar sourceStrength = 0.0;
+    forAll(Yi,cellI)
+    {
+        sourceStrength += (Smi0[cellI]*Yi[cellI]+Smi1[cellI])*Yi.mesh().V()[cellI];
+    }
+    reduce(sourceStrength, sumOp<scalar>());
+    Info << Y[i].name() << ": source strength = " << sourceStrength << endl;
+}

applications/solvers/rcfdemSolverRhoSteadyPimpleChem/pEqn.H

@@ -37,7 +37,8 @@ else
     // Update the pressure BCs to ensure flux consistency
     constrainPressure(p, rhoeps, U, phi, rhorAUf);
 
-    volScalarField SmbyP(particleCloud.chemistryM(0).Sm() / p);
+    volScalarField Sm0(neg(Sm)*Sm/(p + psmall));
+    volScalarField Sm1(pos0(Sm)*Sm);
 
     while (pimple.correctNonOrthogonal())
     {
@@ -47,7 +48,8 @@ else
             fvc::div(phi)
           - fvm::laplacian(rhorAUf, p)
           ==
-            fvm::Sp(SmbyP, p)
+            fvm::Sp(Sm0,p)
+          + Sm1
           + fvOptions(psi, p, rho.name())
         );
 
@@ -60,6 +62,7 @@ else
             phi += pEqn.flux();
         }
     }
+    #include "monitorMassSource.H"
 }
 
 #include "rhoEqn.H"

applications/solvers/rcfdemSolverRhoSteadyPimpleChem/rcfdemSolverRhoSteadyPimpleChem.C

@@ -76,16 +76,12 @@ int main(int argc, char *argv[])
     #include "createFieldRefs.H"
     #include "createFvOptions.H"
 
-    // create cfdemCloud
-    //#include "readGravitationalAcceleration.H"
     cfdemCloudRec<cfdemCloudEnergy> particleCloud(mesh);
     #include "checkModelType.H"
     recBase recurrenceBase(mesh);
     #include "updateFields.H"
 
     turbulence->validate();
-    //#include "compressibleCourantNo.H"
-    //#include "setInitialDeltaT.H"
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -104,6 +100,8 @@ int main(int argc, char *argv[])
     scalar m0(0.0);
     label counter(0);
     p.storePrevIter();
+    const dimensionedScalar psmall("psmall", dimPressure, small);
+    const dimensionedScalar Yismall("Yismall", dimless, small);
 
     while (runTime.run())
     {
@@ -121,9 +119,6 @@ int main(int argc, char *argv[])
         particleCloud.clockM().start(2,"Coupling");
         bool hasEvolved = particleCloud.evolve(voidfraction,Us,U);
 
-        //voidfraction = voidfractionRec;
-        //Us = UsRec;
-
         if(hasEvolved)
         {
             particleCloud.smoothingM().smoothen(particleCloud.forceM(0).impParticleForces());
@@ -160,8 +155,8 @@ int main(int argc, char *argv[])
 
             // --- Pressure-velocity PIMPLE corrector loop
 
-
             #include "UEqn.H"
+            #include "YEqn.H"
             #include "EEqn.H"
 
             // --- Pressure corrector loop
@@ -172,7 +167,6 @@ int main(int argc, char *argv[])
                 #include "pEqn.H"
                 rhoeps=rho*voidfractionRec;
             }
-            #include "YEqn.H"
 
             if (pimple.turbCorr())
             {
@@ -180,8 +174,6 @@ int main(int argc, char *argv[])
             }
         }
 
-        #include "monitorMass.H"
-
         totalStepCounter++;
         particleCloud.clockM().stop("Flow");
 

applications/solvers/rctfSpeciesTransport/Make/options

@@ -1,7 +1,6 @@
 include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
 
 EXE_INC = \
-    -I$(CFDEM_OFVERSION_DIR) \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
@@ -12,6 +11,7 @@ EXE_INC = \
     -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/cfdTools \
     -I$(CFDEM_SRC_DIR)/recurrence/lnInclude \
     -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/derived/cfdemCloudRec \
+    -Wno-deprecated-copy
 
 EXE_LIBS = \
     -L$(CFDEM_LIB_DIR)\

applications/solvers/recSolverTurbTransport/Make/options

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

applications/solvers/rtfmSolverSpecies/Make/options

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

applications/utilities/cfdemPostproc/Make/options

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

applications/utilities/cfdemPostproc/cfdemPostproc.C

@@ -65,22 +65,22 @@ int main(int argc, char *argv[])
     int DEM_dump_Interval(particleCloud.couplingProperties().lookupOrDefault<int>("dumpInterval",1000));
     particleCloud.reAllocArrays();
 
-    double **positions_;
-    double **velocities_;
-    double **radii_;
-    double **voidfractions_;
-    double **particleWeights_;
-    double **particleVolumes_;
-    double **particleV_;
-    int **cellIDs_;
+    double **positions_=NULL;
+    double **velocities_=NULL;
+    double **radii_=NULL;
+    double **voidfractions_=NULL;
+    double **particleWeights_=NULL;
+    double **particleVolumes_=NULL;
+    double **particleV_=NULL;
+    int **cellIDs_=NULL;
 
     particleCloud.dataExchangeM().allocateArray(positions_,0.,3);
     particleCloud.dataExchangeM().allocateArray(velocities_,0.,3);
     particleCloud.get_radii(radii_);  // get ref to radii
     //particleCloud.dataExchangeM().allocateArray(radii_,0.,1);
     particleCloud.dataExchangeM().allocateArray(voidfractions_,0.,1);
-    particleCloud.dataExchangeM().allocateArray(particleWeights_,0.,1);
-    particleCloud.dataExchangeM().allocateArray(particleVolumes_,0.,1);
+    particleCloud.dataExchangeM().allocateArray(particleWeights_,0.,particleCloud.voidFractionM().maxCellsPerParticle());
+    particleCloud.dataExchangeM().allocateArray(particleVolumes_,0.,particleCloud.voidFractionM().maxCellsPerParticle());
     particleCloud.dataExchangeM().allocateArray(particleV_,0.,1);
     particleCloud.get_cellIDs(cellIDs_);  // get ref to cellIDs
     //particleCloud.dataExchangeM().allocateArray(cellIDs_,0.,1);

applications/utilities/displacementField/Make/options

@@ -1,6 +1,7 @@
 EXE_INC = \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -Wno-deprecated-copy
 
 EXE_LIBS = \
     -lfiniteVolume \

applications/utilities/displacementField/displacementField.C

@@ -37,11 +37,19 @@ Application
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 void findPairs(labelList &, labelList &, labelPairList &);
 void findPairsUnordered(labelList &, labelList &, labelPairList &);
-void fillEmptyCells(fvMesh &, label, label, labelList &, volVectorField &, volVectorField &, scalarList &, volVectorField &, volVectorField &, bool, scalar);
-void nearestNeighborCells(fvMesh &, label, label, label, labelList &, labelList &);
-void normalizeFields(labelList &, volVectorField &, volVectorField &);
-void readDump(std::string, labelList &, vectorList &);
+void fillEmptyCells(fvMesh &, label, label, scalarList &, volVectorField &, volVectorField &, scalarList &, volVectorField &, volVectorField &, bool, scalar);
+void nearestNeighborCells(fvMesh &, label, label, label, scalarList &, labelList &);
+void normalizeFields(scalarList &, volVectorField &, volVectorField &);
+void readDump(std::string, labelList &, scalarList &, vectorList &);
 scalar weightFun(scalar);
+label maxNumParticles = 1000000;
+scalar minVol = 1e-12;
+scalar Pi43 = 4.1888;
+label posIndex = -1;
+label posRad = -1;
+label posX = -1;
+label posY = -1;
+label posZ = -1;
 
 int main(int argc, char *argv[])
 {
@@ -88,6 +96,11 @@ int main(int argc, char *argv[])
     label dumpIndexDisplacementIncrement(readLabel(displacementProperties.lookup("dumpIndexDisplacementIncrement")));
     label nNeighMin(readLabel(displacementProperties.lookup("nNeighMin")));
     label maxSearchLayers(displacementProperties.lookupOrDefault<label>("maxSearchLayers",0));
+    posIndex = readLabel(displacementProperties.lookup("posIndex"));
+    posRad = readLabel(displacementProperties.lookup("posRad"));
+    posX = readLabel(displacementProperties.lookup("posX"));
+    posY = readLabel(displacementProperties.lookup("posY"));
+    posZ = readLabel(displacementProperties.lookup("posZ"));
     scalar timePerInputStep(readScalar(displacementProperties.lookup("timePerInputStep")));
     scalar timePerDisplacementStep(readScalar(displacementProperties.lookup("timePerDisplacementStep")));
     scalar startTime(readScalar(displacementProperties.lookup("startTime")));
@@ -95,7 +108,7 @@ int main(int argc, char *argv[])
     std::string fileext=string(displacementProperties.lookupOrDefault<string>("fileextension",""));
     bool interpolate=bool(displacementProperties.lookupOrDefault<bool>("fillEmptyCells",true));
     bool averageMode=bool(displacementProperties.lookupOrDefault<bool>("averageMode",false));
-    
+
     volVectorField defaultUs
     (
         IOobject
@@ -110,11 +123,11 @@ int main(int argc, char *argv[])
         dimensionedVector("zero", dimensionSet(0,1,-1,0,0), vector::zero)
     );
 
-    volVectorField defaultUsDirectedVariance
+    volVectorField defaultUsDirectedStdDev
     (
         IOobject
         (
-            "defaultUDispDirectedVariance",
+            "defaultUDispDirectedStdDev",
             runTime.timeName(),
             mesh,
             IOobject::READ_IF_PRESENT,
@@ -172,11 +185,11 @@ int main(int argc, char *argv[])
         dimensionedVector("zero", dimensionSet(0,1,-1,0,0), vector::zero)
     );
 
-    volVectorField UsDirectedVariance
+    volVectorField UsDirectedStdDev
     (
         IOobject
         (
-            "UDispDirectedVariance",
+            "UDispDirectedStdDev",
             runTime.timeName(),
             mesh,
             IOobject::NO_READ,
@@ -186,7 +199,7 @@ int main(int argc, char *argv[])
         dimensionedVector("zero", dimensionSet(0,1,-1,0,0), vector::zero)
     );
 
-    labelList particlesInCell(mesh.nCells(), 0);
+    scalarList particleVolInCell(mesh.nCells(), 0.0);
 
     scalar currTime=startTime + thisProc * timePerInputStep;
     label timeIndex=thisProc;
@@ -196,6 +209,7 @@ int main(int argc, char *argv[])
         runTime.setTime(currTime,timeIndex);
         // read dump files and check which particle indices are present in both
         labelList indices1, indices2;
+        scalarList radii1, radii2;
         vectorList positions1, positions2;
 
         std::stringstream ss;
@@ -209,13 +223,13 @@ int main(int argc, char *argv[])
         {
             if (averageMode)
             {
-                normalizeFields(particlesInCell, Us, UsDirectedVariance);
-                fillEmptyCells(mesh,nNeighMin,maxSearchLayers,particlesInCell,Us,UsDirectedVariance,boundaries,defaultUs,defaultUsDirectedVariance,interpolate,timePerDisplacementStep);
+                normalizeFields(particleVolInCell, Us, UsDirectedStdDev);
+                fillEmptyCells(mesh,nNeighMin,maxSearchLayers,particleVolInCell,Us,UsDirectedStdDev,boundaries,defaultUs,defaultUsDirectedStdDev,interpolate,timePerDisplacementStep);
 
                 Us /= timePerDisplacementStep;
-                UsDirectedVariance /= timePerDisplacementStep;
+                UsDirectedStdDev /= timePerDisplacementStep;
                 Us.write();
-                UsDirectedVariance.write();
+                UsDirectedStdDev.write();
             }
             break;
         }
@@ -225,8 +239,8 @@ int main(int argc, char *argv[])
         Info << "\t" << filename2 << endl;
         Info << "corresponding to time = " << currTime << "." << endl;
 
-        readDump(filename1, indices1, positions1);
-        readDump(filename2, indices2, positions2);
+        readDump(filename1, indices1, radii1, positions1);
+        readDump(filename2, indices2, radii2, positions2);
 
         labelPairList pairs;
         findPairs(indices1,indices2,pairs);
@@ -234,15 +248,16 @@ int main(int argc, char *argv[])
         // average particle displacements and their variance
         Info << "Binning particle displacements on mesh." << endl;
         vector position, displacement;
+        scalar radius, volume;
         label line1, line2;
         label cellI;
 
         if (!averageMode)
         {
             Us *= 0.0;
-            UsDirectedVariance *= 0.0;
-            particlesInCell.clear();
-            particlesInCell.setSize(mesh.nCells(), 0);
+            UsDirectedStdDev *= 0.0;
+            particleVolInCell.clear();
+            particleVolInCell.setSize(mesh.nCells(), 0);
         }
 
         for (label partI = 0; partI < pairs.size(); partI++)
@@ -250,27 +265,29 @@ int main(int argc, char *argv[])
             line1 = pairs[partI].first();
             line2 = pairs[partI].second();
             position = positions1[line1];
-            displacement = positions2[line2] - positions1[line1];
             cellI = mesh.findCell(position);
             if (cellI < 0)  continue;
-            particlesInCell[cellI] += 1;
-            Us[cellI] += displacement;
+            displacement = positions2[line2] - positions1[line1];
+            radius = radii1[line1];
+            volume = Pi43 * radius * radius * radius;
+            particleVolInCell[cellI] += volume;
+            Us[cellI] += displacement*volume;
 
             for (label comp=0;comp<3;comp++)
             {
-                UsDirectedVariance[cellI].component(comp) += displacement.component(comp)*displacement.component(comp);
+                UsDirectedStdDev[cellI].component(comp) += displacement.component(comp)*displacement.component(comp)*volume;
             }
         }
 
         if (!averageMode)
         {
-            normalizeFields(particlesInCell, Us, UsDirectedVariance);
-            fillEmptyCells(mesh,nNeighMin,maxSearchLayers,particlesInCell,Us,UsDirectedVariance,boundaries,defaultUs,defaultUsDirectedVariance,interpolate,timePerDisplacementStep);
+            normalizeFields(particleVolInCell, Us, UsDirectedStdDev);
+            fillEmptyCells(mesh,nNeighMin,maxSearchLayers,particleVolInCell,Us,UsDirectedStdDev,boundaries,defaultUs,defaultUsDirectedStdDev,interpolate,timePerDisplacementStep);
 
             Us /= timePerDisplacementStep;
-            UsDirectedVariance /= timePerDisplacementStep;
+            UsDirectedStdDev /= timePerDisplacementStep;
             Us.write();
-            UsDirectedVariance.write();
+            UsDirectedStdDev.write();
         }
 
         if (averageMode && monitorProbes)
@@ -280,7 +297,7 @@ int main(int argc, char *argv[])
             {
                 vector pos = probePoints[p];
                 label cellP = mesh.findCell(pos);
-                monitoringDataFile << " " << particlesInCell[cellP] << " " << Us[cellP]/timePerDisplacementStep << " " << UsDirectedVariance[cellP]/(timePerDisplacementStep*timePerDisplacementStep);
+                monitoringDataFile << " " << particleVolInCell[cellP] << " " << Us[cellP]/timePerDisplacementStep << " " << UsDirectedStdDev[cellP]/(timePerDisplacementStep*timePerDisplacementStep);
             }
             monitoringDataFile << endl;
         }
@@ -293,30 +310,69 @@ int main(int argc, char *argv[])
     return 0;
 }
 
-void readDump(std::string filename, labelList &indices, vectorList &positions)
+void readDump(std::string filename, labelList &indices, scalarList &radii, vectorList &positions)
 {
     #include <fstream>
 
     const label leadingLines = 9;
     label lineCounter = 0;
-    label partIndex;
-    scalar x, y, z;
+    label partIndex = 0;
+    scalar r = 1.0, x = 0.0, y = 0.0, z = 0.0;
 
     indices.clear();
+    radii.clear();
     positions.clear();
 
+    indices.setSize(maxNumParticles);
+    radii.setSize(maxNumParticles);
+    positions.setSize(maxNumParticles);
+
     std::ifstream file(filename);
     std::string str;
+    std::string word;
+    label wordcounter;
     while (std::getline(file, str))
     {
         if (lineCounter >= leadingLines)
         {
-            sscanf(str.c_str(), "%d %lf %lf %lf", &partIndex, &x, &y, &z);
-            indices.append(partIndex);
-            positions.append(vector(x,y,z));
+            std::istringstream ss(str);
+            wordcounter = 0;
+            while (ss >> word)
+            {
+                if (wordcounter == posIndex)
+                {
+                    partIndex = stoi(word);
+                }
+                else if (wordcounter == posRad)
+                {
+                    r = stod(word);
+                }
+                else if (wordcounter == posX)
+                {
+                    x = stod(word);
+                }
+                else if (wordcounter == posY)
+                {
+                    y = stod(word);
+                }
+                else if (wordcounter == posZ)
+                {
+                    z = stod(word);
+                }
+                wordcounter++;
+            }
+//            sscanf(str.c_str(), "%d %lf %lf %lf", &partIndex, &x, &y, &z);
+            indices[lineCounter-leadingLines] = partIndex;
+            radii[lineCounter-leadingLines] = r;
+            positions[lineCounter-leadingLines] = vector(x,y,z);
         }
         lineCounter++;
     }
+
+    label readLines = lineCounter - leadingLines;
+    indices.resize(readLines);
+    radii.resize(readLines);
+    positions.resize(readLines);
 }
 
 void findPairs(labelList &indices1, labelList &indices2, labelPairList &pairs)
@@ -324,6 +380,10 @@ void findPairs(labelList &indices1, labelList &indices2, labelPairList &pairs)
     // remove all entries from first list if they are not present in second list
     // this assumes ordered entries
 
+    pairs.clear();
+    pairs.setSize(maxNumParticles);
+    label pairCounter = 0;
+
     if (indices2.size() == 0) return;
 
     for (label i=0;i<indices1.size();i++)
@@ -339,18 +399,23 @@ void findPairs(labelList &indices1, labelList &indices2, labelPairList &pairs)
             else if (indices2[jmid] < index1) j1 = jmid;
             else
             {
-                pairs.append(labelPair(i,jmid));
+                pairs[pairCounter]=labelPair(i,jmid);
+                pairCounter++;
                 break;
             }
             if (j2-j1 == 1) break;
         }
     }
+    pairs.resize(pairCounter);
     Info << "findPairs: " << pairs.size() << " pairs found." << endl;
 }
 
 void findPairsUnordered(labelList &indices1, labelList &indices2, labelPairList &pairs)
 {
     // remove all entries from first list if they are not present in second list
+    pairs.clear();
+    pairs.setSize(maxNumParticles);
+    label pairCounter = 0;
 
     for (label i=0;i<indices1.size();i++)
     {
@@ -358,15 +423,17 @@ void findPairsUnordered(labelList &indices1, labelList &indices2, labelPairList
         {
             if (indices1[i] == indices2[j])
             {
-                pairs.append(labelPair(i,j));
+                pairs[pairCounter]=labelPair(i,j);
+                pairCounter++;
                 break;
             }
         }
     }
+    pairs.resize(pairCounter);
     Info << "findPairs: " << pairs.size() << " pairs found." << endl;
 }
 
-void fillEmptyCells(fvMesh &mesh, label nNeighMin, label maxSearchLayers, labelList &particlesInCell, volVectorField &Us, volVectorField& UsDirectedVariance,scalarList& boundaries, volVectorField &defaultUs, volVectorField &defaultUsDirectedVariance, bool interpolate, scalar dt)
+void fillEmptyCells(fvMesh &mesh, label nNeighMin, label maxSearchLayers, scalarList &particleVolInCell, volVectorField &Us, volVectorField& UsDirectedStdDev,scalarList& boundaries, volVectorField &defaultUs, volVectorField &defaultUsDirectedStdDev, bool interpolate, scalar dt)
 {
     labelList neighborsWithValues;
     scalar neighborSqrDistance;
@@ -377,7 +444,7 @@ void fillEmptyCells(fvMesh &mesh, label nNeighMin, label maxSearchLayers, labelL
     Info << "Filling empty cells." << endl;
     forAll(mesh.C(), cellI)
     {
-        if (particlesInCell[cellI] > 0) continue;
+        if (particleVolInCell[cellI] > minVol) continue;
 
         vector position = mesh.C()[cellI];
         label outsideBox = 0;
@@ -388,11 +455,11 @@ void fillEmptyCells(fvMesh &mesh, label nNeighMin, label maxSearchLayers, labelL
         if (outsideBox > 0 || !interpolate)
         {
             Us[cellI] = defaultUs[cellI]*dt;
-            UsDirectedVariance[cellI] = defaultUsDirectedVariance[cellI]*dt;
+            UsDirectedStdDev[cellI] = defaultUsDirectedStdDev[cellI]*dt;
             continue;
         }
 
-        nearestNeighborCells(mesh, cellI, nNeighMin, maxSearchLayers, particlesInCell, neighborsWithValues);
+        nearestNeighborCells(mesh, cellI, nNeighMin, maxSearchLayers, particleVolInCell, neighborsWithValues);
         weightSum = 0.0;
         weights.clear();
         for (label neighI=0; neighI<neighborsWithValues.size(); neighI++)
@@ -406,13 +473,13 @@ void fillEmptyCells(fvMesh &mesh, label nNeighMin, label maxSearchLayers, labelL
         {
             weight = weights[neighI]/weightSum;
             Us[cellI] += weight*Us[neighborsWithValues[neighI]];
-            UsDirectedVariance[cellI] += weight*UsDirectedVariance[neighborsWithValues[neighI]];
+            UsDirectedStdDev[cellI] += weight*UsDirectedStdDev[neighborsWithValues[neighI]];
         }
-        
+
         if (neighborsWithValues.size() == 0)
         {
             Us[cellI] = defaultUs[cellI]*dt;
-            UsDirectedVariance[cellI] = defaultUsDirectedVariance[cellI]*dt;
+            UsDirectedStdDev[cellI] = defaultUsDirectedStdDev[cellI]*dt;
         }
 
         // make sure no particles are placed outside of domain
@@ -429,7 +496,7 @@ void fillEmptyCells(fvMesh &mesh, label nNeighMin, label maxSearchLayers, labelL
     }
 }
 
-void nearestNeighborCells(fvMesh &mesh, label refCell, label nNeighMin, label maxSearchLayers, labelList &particlesInCell, labelList &neighborsWithValues)
+void nearestNeighborCells(fvMesh &mesh, label refCell, label nNeighMin, label maxSearchLayers, scalarList &particleVolInCell, labelList &neighborsWithValues)
 {
     label numSearchLayers = 0;
     std::set<label> neighbors;
@@ -455,11 +522,11 @@ void nearestNeighborCells(fvMesh &mesh, label refCell, label nNeighMin, label ma
                 {
                     newNeighbors.insert(adj);
                     neighbors.insert(adj);
-                    if (particlesInCell[adj] > 0) neighborsWithValues.append(adj);
+                    if (particleVolInCell[adj] > minVol) neighborsWithValues.append(adj);
                 }
             }
         }
-        
+
         numSearchLayers++;
         if (numSearchLayers > maxSearchLayers && maxSearchLayers > 0) return;
 
@@ -470,18 +537,18 @@ void nearestNeighborCells(fvMesh &mesh, label refCell, label nNeighMin, label ma
     }
 }
 
-void normalizeFields(labelList& particlesInCell, volVectorField& Us, volVectorField & UsDirectedVariance)
+void normalizeFields(scalarList& particleVolInCell, volVectorField& Us, volVectorField & UsDirectedStdDev)
 {
-    for (label cellJ = 0; cellJ<particlesInCell.size(); cellJ++)
+    for (label cellJ = 0; cellJ<particleVolInCell.size(); cellJ++)
     {
-        if (particlesInCell[cellJ] > 0)
+        if (particleVolInCell[cellJ] > minVol)
         {
-            Us[cellJ] /= particlesInCell[cellJ];
-            UsDirectedVariance[cellJ] /= particlesInCell[cellJ];
+            Us[cellJ] /= particleVolInCell[cellJ];
+            UsDirectedStdDev[cellJ] /= particleVolInCell[cellJ];
             for (label comp=0;comp<3;comp++)
             {
-                UsDirectedVariance[cellJ].component(comp) -= Us[cellJ].component(comp)*Us[cellJ].component(comp);
-                if (UsDirectedVariance[cellJ].component(comp) > 0) UsDirectedVariance[cellJ].component(comp) = Foam::sqrt(UsDirectedVariance[cellJ].component(comp));
+                UsDirectedStdDev[cellJ].component(comp) -= Us[cellJ].component(comp)*Us[cellJ].component(comp);
+                if (UsDirectedStdDev[cellJ].component(comp) > 0) UsDirectedStdDev[cellJ].component(comp) = Foam::sqrt(UsDirectedStdDev[cellJ].component(comp));
             }
         }
     }

applications/utilities/rSmoother/Make/options

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

doc/CFDEMcoupling_about.txt

@@ -44,7 +44,7 @@ the use of Git allows to easily update to the latest version :l
 basic documentation is provided :l
 :ule
 
-The file structure:
+The project structure:
 
 {src} directory including the source files of the coupling toolbox and models :ulb,l
 {applications} directory including the solver files for coupled CFD-DEM simulations :l
@@ -52,7 +52,8 @@ The file structure:
 {tutorials} directory including basic tutorial cases showing the functionality  :l
 :ule
 
-Details on installation are given on the "CFDEMproject Website"_lig .
+Details on installation are given in the "Installation"_CFDEMcoupling_install.html
+section.
 The functionality of this CFD-DEM framework is described via "tutorial
 cases"_CFDEMcoupling_tutorials.html showing how to use different solvers and
 models.

doc/CFDEMcoupling_install.txt

@@ -45,6 +45,11 @@ following commands:
 apt-get install build-essential flex bison git-core cmake zlib1g-dev libboost-system-dev libboost-thread-dev libopenmpi-dev openmpi-bin gnuplot libreadline-dev libncurses-dev libxt-dev
 apt-get install libqt5x11extras5-dev libxt-dev qt5-default qttools5-dev curl :pre
 
+NOTE: Ubuntu 21.04 and newer dropped the {qt5-default} package. Instead use:
+
+apt-get install build-essential flex bison git-core cmake zlib1g-dev libboost-system-dev libboost-thread-dev libopenmpi-dev openmpi-bin gnuplot libreadline-dev libncurses-dev libxt-dev
+apt-get install libqt5x11extras5-dev libxt-dev qttools5-dev qttools5-dev-tools qtwebengine5-dev libqt5svg5-dev libqt5websockets5-dev libqt5xmlpatterns5 qtxmlpatterns5-dev-tools curl :pre
+
 2.1.2 Setup the environment :h5
 
 Open your bash startup file

doc/CFDEMcoupling_models.txt

@@ -122,8 +122,10 @@ isotropicFluctuations,
 particleDeformation,
 "pdCorrelation"_forceModel_pdCorrelation.html,
 potentialRelaxation,
+"staticPressure"_forceModel_staticPressure.html,
 "surfaceTensionForce"_forceModel_surfaceTensionForce.html,
 terminalVelocity,
+"transferFluidProperties"_forceModel_transferFluidProperties.html,
 turbulentDispersion,
 turbulentVelocityFluctuations,
 "virtualMassForce"_forceModel_virtualMassForce.html,
@@ -198,6 +200,7 @@ each DEM particle.
 
 expParticleForces,
 gravity,
+"gravityEff"_otherForceModel_gravityEff.html,
 secondaryPhaseInducedBuoyancy,
 weightSecondaryPhase :tb(c=2,ea=c)
 

doc/CFDEMcoupling_solvers.txt

@@ -10,8 +10,10 @@
 This section lists all CFDEMcoupling solvers alphabetically.
 
 "cfdemSolverIB"_cfdemSolverIB.html,
+"cfdemSolverIBContinuousForcing"_cfdemSolverIBContinuousForcing.html,
 "cfdemSolverMultiphase"_cfdemSolverMultiphase.html,
 "cfdemSolverMultiphaseScalar"_cfdemSolverMultiphaseScalar.html,
+"cfdemSolverPimple"_cfdemSolverPimple.html,
 "cfdemSolverPiso"_cfdemSolverPiso.html,
 "cfdemSolverPisoScalar"_cfdemSolverPisoScalar.html,
 "cfdemSolverRhoPimple"_cfdemSolverRhoPimple.html,

doc/CFDEMcoupling_tutorials.txt

@@ -62,7 +62,7 @@ OpenFOAM\ |reg|\ (*) documentations.
 END_RST -->.
 
 Settings of the coupling routines are defined in
-$caseDir/CFD/constant/"couplingProperies"_CFDEMcoupling_dicts.html#couplingProperties
+$caseDir/CFD/constant/"couplingProperties"_CFDEMcoupling_dicts.html#couplingProperties
 (e.g. force models, data exchange model, etc.) and
 $caseDir/CFD/constant/"liggghtsCommands"_CFDEMcoupling_dicts.html#liggghtsCommands
 (allows to execute a LIGGGHTS command during a coupled simulation).

doc/IOModel_basicIO.txt

@@ -17,7 +17,8 @@ IOModel basicIO; :pre
 [Examples:]
 
 IOModel basicIO;
-serialOutput; :pre
+serialOutput;
+cartesianOutput; :pre
 
 [Description:]
 
@@ -29,6 +30,11 @@ Using the keyword {serialOutput;} in the
 the IO is serial to the directory {$casePath/CFD/lagrangian}. In this case
 only the data on processor 0 is written!
 
+Using the keyword {cartesianOutput;} in the
+"couplingProperties"_CFDEMcoupling_dicts.html#couplingProperties dictionary,
+particle positions will be written as cartesian coordinates instead of
+barycentric coordinates.
+
 Data is written every write time of the CFD simulation.
 
 [Restrictions:]

doc/cfdemSolverIBContinuousForcing.txt

@@ -0,0 +1,82 @@
+"CFDEMproject Website"_lws - "Main Page"_main :c
+
+:link(lws,http://www.cfdem.com)
+:link(main,CFDEMcoupling_Manual.html)
+
+:line
+
+cfdemSolverIBContinuousForcing command :h3
+
+[Description:]
+
+"cfdemSolverIBContinuousForcing" is a coupled CFD-DEM solver using CFDEMcoupling,
+an open-source parallel coupled CFD-DEM framework, for calculating the dynamics
+between immersed bodies and the surrounding fluid. Being an implementation of a
+continuous forcing immersed boundary method it allows tackling problems where
+the body diameter exceeds the maximal size of a fluid cell.
+
+<!-- HTML_ONLY -->
+Using the toolbox of OpenFOAM&reg;(*) the governing equations of the fluid are
+computed and the corrections of velocity and pressure field with respect to the
+body-movement information, gained from LIGGGHTS, are incorporated.
+<!-- END_HTML_ONLY -->
+
+<!-- RST
+
+Using the toolbox of OpenFOAM\ |reg|\ (*) the governing equations of the fluid are
+computed and the corrections of velocity and pressure field with respect to the
+body-movement information, gained from LIGGGHTS, are incorporated.
+
+.. |reg|    unicode:: U+000AE .. REGISTERED SIGN
+
+END_RST -->
+
+The code of this solver was contributed by A.N. Balachandran Nair, JKU. For more
+details, see "Balachandran Nair et al. (2021)"_#BalachandranNair2021
+
+[Use:]
+
+The solver is realized within the open-source framework CFDEMcoupling. Just as
+for the unresolved CFD-DEM solver cfdemSolverPiso the file
+CFD/constant/couplingProperties contains information about the settings for the
+different models. While IOmodel, DataExchangeModel etc. are applicable for all
+CFDEMcoupling-solvers, special locate-, force- and void fraction models were
+designed for this solver:
+
+"engineSearchIB"_locateModel_engineSearchIB.html,
+"ArchimedesIB"_forceModel_ArchimedesIB.html,
+"ShirgaonkarIB"_forceModel_ShirgaonkarIB.html,
+"IBVoidfraction"_voidFractionModel_IBVoidFraction.html
+
+:line
+
+:link(BalachandranNair2021)
+[(Balachandran Nair, 2021)] Balachandran Nair, A.N., Pirker, S. and Saeedipour, M.,
+"Resolved CFD-DEM simulation of blood flow with a reduced-order RBC model",
+Comp. Part. Mech. (2021)
+
+:line
+
+<!-- HTML_ONLY -->
+NOTE:
+(*) This offering is not approved or endorsed by OpenCFD Limited, producer and
+distributor of the OpenFOAM software via www.openfoam.com, and owner of the
+OPENFOAM&reg; and OpenCFD&reg; trade marks.
+OPENFOAM&reg; is a registered trade mark of OpenCFD Limited, producer and
+distributor of the OpenFOAM software via www.openfoam.com.
+<!-- END_HTML_ONLY -->
+
+<!-- RST
+
+.. note::
+
+   (*) This offering is not approved or endorsed by OpenCFD Limited, producer
+   and distributor of the OpenFOAM software via www.openfoam.com, and owner of
+   the OPENFOAM\ |reg| and OpenCFD\ |reg| trade marks.
+   OPENFOAM\ |reg| is a registered trade mark of OpenCFD Limited, producer and
+   distributor of the OpenFOAM software via www.openfoam.com.
+
+.. |reg|    unicode:: U+000AE .. REGISTERED SIGN
+
+END_RST -->
+

doc/cfdemSolverPimple.txt

@@ -0,0 +1,67 @@
+<!-- HTML_ONLY -->
+<HEAD>
+  <META CHARSET="utf-8">
+</HEAD>
+<!-- END_HTML_ONLY -->
+"CFDEMproject Website"_lws - "Main Page"_main :c
+
+:link(lws,http://www.cfdem.com)
+:link(main,CFDEMcoupling_Manual.html)
+
+:line
+
+cfdemSolverPimple command :h3
+
+[Description:]
+
+<!-- HTML_ONLY -->
+"cfdemSolverPimple" is a coupled CFD-DEM solver using CFDEMcoupling, an open-\
+source parallel coupled CFD-DEM framework. Based on DPMFoam&reg;(*), a finite
+volume based solver for turbulent Navier-Stokes equations applying the PIMPLE
+algorithm, "cfdemSolverPimple" has additional functionality for a coupling to the
+DEM code "LIGGGHTS".
+<!-- END_HTML_ONLY -->
+
+<!-- RST
+
+"cfdemSolverPimple" is a coupled CFD-DEM solver using CFDEMcoupling, an open-\
+source parallel coupled CFD-DEM framework. Based on DPMFoam\ |reg|\ (*), a finite
+volume based solver for turbulent Navier-Stokes equations applying the PIMPLE
+algorithm, "cfdemSolverPimple" has additional functionality for a coupling to the
+DEM code "LIGGGHTS".
+
+.. |reg|    unicode:: U+000AE .. REGISTERED SIGN
+
+END_RST -->
+
+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.
+
+
+
+
+:line
+
+<!-- HTML_ONLY -->
+NOTE:
+(*) This offering is not approved or endorsed by OpenCFD Limited, producer and
+distributor of the OpenFOAM software via www.openfoam.com, and owner of the
+OPENFOAM&reg; and OpenCFD&reg; trade marks.
+OPENFOAM&reg; is a registered trade mark of OpenCFD Limited, producer and
+distributor of the OpenFOAM software via www.openfoam.com.
+<!-- END_HTML_ONLY -->
+
+<!-- RST
+
+.. note::
+
+   (*) This offering is not approved or endorsed by OpenCFD Limited, producer
+   and distributor of the OpenFOAM software via www.openfoam.com, and owner of
+   the OPENFOAM\ |reg| and OpenCFD\ |reg| trade marks.
+   OPENFOAM\ |reg| is a registered trade mark of OpenCFD Limited, producer and
+   distributor of the OpenFOAM software via www.openfoam.com.
+
+.. |reg|    unicode:: U+000AE .. REGISTERED SIGN
+
+END_RST -->

doc/forceModel_ParmarBassetForce.txt

@@ -23,17 +23,19 @@ ParmarBassetForceProps
     nIntegral           scalar1;
     discretisationOrder scalar2;
     treatForceExplicit  switch1;
-    interpolation       switch2;
+    verbose             switch2;
+    interpolation       switch3;
     smoothingModel      "smoothingModel";
 \} :pre
 
 {U} = name of the finite volume fluid velocity field :ulb,l
 {Us} = name of the finite volume cell averaged particle velocity field :l
-{scalar1} = number of timesteps considered in the near history :l
-{scalar2} = (1 or 2) discretisation order of the far history differential equations :l
+{scalar1} = number of timesteps considered in the near history (integer > 1):l
+{scalar2} = (optional, default 1) discretisation order of the far history differential equations (1 or 2) :l
 {switch1} = (optional, default true) sub model switch, see "forceSubModel"_forceSubModel.html for details :l
 {switch2} = (optional, default false) sub model switch, see "forceSubModel"_forceSubModel.html for details :l
-{smoothingModel} = name of smoothing model for the dU/dt field :l
+{switch3} = (optional, default false) sub model switch, see "forceSubModel"_forceSubModel.html for details :l
+{smoothingModel} = name of smoothing model for the dUrel/dt field :l
 :ule
 
 [Examples:]

doc/forceModel_ShirgaonkarIB.txt

@@ -21,12 +21,14 @@ ShirgaonkarIBProps
     velFieldName      "U";
     pressureFieldName "pressure";
     twoDimensional;
+    useTorque;
     treatForceExplicit switch1;
 \} :pre
 
 {U} = name of the finite volume fluid velocity field :ulb,l
 {pressure} = name of the finite volume pressure field :l
 {twoDimensional} = optional keyword for conducting a two dimensional calculation :l
+{useTorque} =  optional keyword for calculating particle torque :l
 {switch1} = (optional, default false) sub model switch, see "forceSubModel"_forceSubModel.html for details :l
 :ule
 

doc/forceModel_staticPressure.txt

@@ -0,0 +1,70 @@
+"CFDEMproject Website"_lws - "Main Page"_main :c
+
+:link(lws,http://www.cfdem.com)
+:link(main,CFDEMcoupling_Manual.html)
+
+:line
+
+forceModel staticPressure command :h3
+
+[Syntax:]
+
+Defined in "couplingProperties"_CFDEMcoupling_dicts.html#couplingProperties
+dictionary.
+
+forceModels
+(
+   staticPressure
+);
+staticPressureProps
+\{
+   rhoGas number1; 
+   voidfractionFieldName "voidfraction";
+   rhoPart number2; 
+   DomainSolidVolumeFraction number3; 
+   g_DEM gravityVector; 
+\} :pre
+
+{number1} = density of the gas phase :ulb,l
+{voidfraction} = name of the finite volume void fraction field :l
+{number2} = density of the particles :l
+{number3} = domain-average particle volume fraction (mass loading) :l
+{gravityVector} = gravity vector :l
+:ule
+
+[Examples:]
+
+forceModels
+(
+   staticPressure
+);
+staticPressureProps
+\{
+   rhoGas 1.2; 
+   voidfractionFieldName "voidfraction";
+   rhoPart 1500; 
+   DomainSolidVolumeFraction 0.05; 
+   g_DEM (0 0 -9.81); 
+\} :pre
+
+[Description:]
+
+As discussed in "gravityEff"_otherForceModel_gravityEff.html, in case of
+simulating periodic boxes, it is necessary to have a gravitational source term.
+Since that force is added to the gas momentum equation, we need to add the same
+force excluding the effect of the gravitational force on the particles
+(the gravitational force on the particles is already available on the LIGGGHTS side)
+to drive the particles.
+This force model performs the calculation of the following force acting on each
+DEM particle,
+
+:c,image(Eqs/forceModels_staticPressure_pic1.jpg)
+
+[Restrictions:]
+
+none
+
+[Related commands:]
+
+"forceModel"_forceModel.html
+

doc/forceModel_transferFluidProperties.txt

@@ -0,0 +1,42 @@
+"CFDEMproject Website"_lws - "Main Page"_main :c
+
+:link(lws,http://www.cfdem.com)
+:link(main,CFDEMcoupling_Manual.html)
+
+:line
+
+forceModel transferFluidProperties command :h3
+
+[Syntax:]
+
+Defined in "couplingProperties"_CFDEMcoupling_dicts.html#couplingProperties
+dictionary.
+
+forceModels
+(
+    transferFluidProperties
+);
+transferFluidPropertiesProps
+\{
+    verbose       switch1;
+    interpolation switch2;
+\} :pre
+
+{switch1} = (optional, default false) sub model switch, see "forceSubModel"_forceSubModel.html for details :ulb,l
+{switch2} = (optional, default false) sub model switch, see "forceSubModel"_forceSubModel.html for details :l
+:ule
+
+[Description:]
+
+This "force model" does not influence the particles or the flow - it transfer to fluid density and (dynamic)
+viscosity from OpenFOAM to LIGGGHTS.
+
+
+[Restrictions:]
+
+This model requires {fix cfd/coupling/fluidproperties} to work.
+
+[Related commands:]
+
+"forceModel"_forceModel.html
+

doc/otherForceModel_gravityEff.txt

@@ -0,0 +1,80 @@
+"CFDEMproject Website"_lws - "Main Page"_main :c
+
+:link(lws,http://www.cfdem.com)
+:link(main,CFDEMcoupling_Manual.html)
+
+:line
+
+otherForceModel gravityEff command :h3
+
+[Syntax:]
+
+Defined in "couplingProperties"_CFDEMcoupling_dicts.html#couplingProperties
+dictionary.
+
+otherForceModels
+(
+  gravityEff
+);
+gravityEffProps
+\{
+    voidfractionFieldName "voidfraction";
+    rhoPart number1;
+\} :pre
+
+{voidfraction} = name of the finite volume void fraction field :l
+{number1} = particle density :l
+
+:ule
+
+[Examples:]
+
+otherForceModels
+(
+  gravityEff
+);
+gravityEffProps
+\{
+  voidfractionFieldName "voidfraction";
+  rhoPart 1500;
+\} :pre
+
+[Description:]
+
+An effective gravitational source term, {gravityEff}, that is necessery to drive
+the flow only incase of periodic box simulations.
+This source term is added to the momentum equation.
+In case of simulating periodic boxes and assuming "homogeneous no-flow condition"
+we will have:
+
+:c,image(Eqs/otherForceModels_gravityEff_pic1.jpg)
+
+This equation depicts that in "no-flow condition" static pressure gradient balances the mixture weight.
+Therefore, to derive the flow in periodic box simulations, a dynamic pressure gradient is defined as
+
+:c,image(Eqs/otherForceModels_gravityEff_pic2.jpg)
+
+If one replaces the static pressure gradient of the fluid momentum equation with
+
+:c,image(Eqs/otherForceModels_gravityEff_pic3.jpg)
+
+the fluid momentum equation solved in case of periodic box simulations will be
+
+:c,image(Eqs/otherForceModels_gravityEff_pic4.jpg)
+
+The {gravityEff} source term is the third term on the RHS of the momentum equation.
+Note that the gravitational body force of the fluid momentum equation is considered
+with this source term.
+
+In the basic implementation of {gravityEff}, first we calculate the density of
+the mixture and then the source term is calculated based on the third term on the
+RHS of the momentum equation.
+
+[Restrictions:]
+
+none
+
+[Related commands:]
+
+"otherForceModels"_otherForceModels.html
+

doc/smoothingModel_constDiffSmoothing.txt

@@ -15,17 +15,22 @@ dictionary.
 smoothingModel constDiffSmoothing;
 constDiffSmoothingProps
 \{
-    lowerLimit                    number1;
-    upperLimit                    number2;
-    smoothingLength               lengthScale;
-    smoothingLengthReferenceField lengthScaleRefField;
+    lowerLimit                        number1;
+    upperLimit                        number2;
+    smoothingLength                   lengthScale;
+    smoothingLengthReference          lengthScaleRef;
+    smoothingLengthFieldName          fieldName1;
+    smoothingLengthReferenceFieldName fieldName2;
     verbose;
 \} :pre
 
 {number1} = scalar fields will be bound to this lower value :ulb,l
 {number2} = scalar fields will be bound to this upper value :l
 {lengthScale} = length scale over which the exchange fields will be smoothed out :l
-{lengthScaleRefField} = length scale over which reference fields (e.g., the average particle velocity) will be smoothed out. Should be always larger than lengthScale. If not specified, will be equal to lengthScale. :l
+{lengthScaleRef} = (optional) length scale over which reference fields (e.g., the average particle velocity) will be smoothed out. Should be always larger than lengthScale. If not specified, will be equal to lengthScale. :l
+{fieldName1} = (optional) name of scalar field to be used as local smoothing length. :l
+{fieldName2} = (optional) name of scalar field to be used as local smoothing length for reference fields. :l
+
 {verbose} = (optional, default false) flag for debugging output :l
 :ule
 
@@ -51,6 +56,11 @@ which these reference fields are not specified. Values calculated in the cells
 e.g. the average particle velocity, which are not specified in all cells in case
 the flow is rather dilute.
 
+Alternative to {smoothingLength} and {smoothingLengthReference}, 
+{smoothingLengthFieldName} and/or {smoothingLengthReferenceFieldName} can be used
+to define spatial variation of the smoothing lengths. Either the scalar or field 
+options must be used, giving both will result in errors.
+
 [Restrictions:]
 
 This model is tested in a limited number of flow situations.

etc/OFversion/OFversion.H

@@ -1,2 +0,0 @@
-#define version4x
-

etc/bashrc

@@ -17,7 +17,7 @@
 #------------------------------------------------------------------------------
 
 export CFDEM_PROJECT=CFDEM
-export CFDEM_VERSION=21.11
+export CFDEM_VERSION=24.01
 
 ################################################################################
 # USER EDITABLE PART: Changes made here may be lost with the next upgrade
@@ -106,9 +106,6 @@ export CFDEM_APP_DIR=$CFDEM_PROJECT_DIR/platforms/$WM_OPTIONS/bin
 mkdir -p $CFDEM_LIB_DIR
 mkdir -p $CFDEM_APP_DIR
 
-#- path to OF version flag file
-export CFDEM_OFVERSION_DIR=$CFDEM_PROJECT_DIR/etc/OFversion
-
 #------------------------------------------------------------------------------
 # Clean standard environment variables (PATH, LD_LIBRARY_PATH)
 #------------------------------------------------------------------------------

etc/cshrc

@@ -15,7 +15,7 @@
 #------------------------------------------------------------------------------
 
 setenv CFDEM_PROJECT CFDEM
-setenv CFDEM_VERSION 21.11
+setenv CFDEM_VERSION 24.01
 
 ################################################################################
 # USER EDITABLE PART: Changes made here may be lost with the next upgrade
@@ -139,9 +139,6 @@ setenv CFDEM_APP_DIR $CFDEM_PROJECT_DIR/platforms/$WM_OPTIONS/bin
 mkdir -p $CFDEM_LIB_DIR
 mkdir -p $CFDEM_APP_DIR
 
-#- path to OF version flag file
-setenv CFDEM_OFVERSION_DIR $CFDEM_PROJECT_DIR/etc/OFversion
-
 #------------------------------------------------------------------------------
 # Clean standard environment variables (PATH, LD_LIBRARY_PATH)
 #------------------------------------------------------------------------------

etc/functions.sh

@@ -359,7 +359,7 @@ cleanCFDEMcase()
     cd $casePath/CFD
     cleanCase
     rm -r $casePath/DEM/post/*
-    echo "dummyfile" >> $casePath/DEM/post/dummy
+    touch $casePath/DEM/post/.gitignore
     cd $casePath
     echo "done"
 }

etc/library-list.txt

@@ -3,4 +3,4 @@ lagrangian/cfdemParticleComp/dir
 recurrence/dir
 finiteVolume/dir
 ../applications/solvers/cfdemSolverMultiphase/multiphaseMixture/dir
-../applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixture/dir
+../applications/solvers/cfdemSolverMultiphaseScalar/multiphaseMixtureScalar/dir

etc/solver-list.txt

@@ -7,9 +7,11 @@ rcfdemSolverBase/dir
 rtfmSolverSpecies/dir
 cfdemSolverPisoMS/dir
 cfdemSolverPiso/dir
+cfdemSolverPimple/dir
 cfdemSolverRhoPimple/dir
 cfdemSolverIB/dir
 cfdemSolverPisoScalar/dir
+cfdemSolverIBContinuousForcing/dir
 cfdemSolverRhoPimpleChem/dir
 cfdemSolverMultiphase/dir
 cfdemSolverMultiphaseScalar/dir

etc/tutorial-list.txt

@@ -19,3 +19,5 @@ cfdemSolverIB/twoSpheresGlowinskiMPI/dir
 cfdemSolverPisoScalar/packedBedTemp/dir
 
 cfdemSolverPiso/ErgunTestCG/dir
+
+cfdemSolverIB/falling_sphere_two_way_coupling/dir

src/lagrangian/cfdemParticle/Make/files

@@ -96,6 +96,7 @@ $(forceModels)/potentialRelaxation/potentialRelaxation.C
 $(forceModels)/BeetstraDrag/BeetstraDrag.C
 $(forceModels)/BeetstraDragPoly/BeetstraDragPoly.C
 $(forceModels)/dSauter/dSauter.C
+$(forceModels)/transferFluidProperties/transferFluidProperties.C
 $(forceModels)/Fines/Fines.C
 $(forceModels)/Fines/FinesFields.C
 $(forceModels)/Fines/FanningDynFines.C
@@ -108,6 +109,7 @@ $(forceModels)/particleDeformation/particleDeformation.C
 $(forceModels)/turbulentDispersion/turbulentDispersion.C
 $(forceModels)/terminalVelocity/terminalVelocity.C
 $(forceModels)/turbulentVelocityFluctuations/turbulentVelocityFluctuations.C
+$(forceModels)/staticPressure/staticPressure.C
 
 $(forceModelsMS)/forceModelMS/forceModelMS.C
 $(forceModelsMS)/forceModelMS/newForceModelMS.C
@@ -125,6 +127,7 @@ $(otherForceModels)/gravity/gravity.C
 $(otherForceModels)/weightSecondaryPhase/weightSecondaryPhase.C
 $(otherForceModels)/expParticleForces/expParticleForces.C
 $(otherForceModels)/secondaryPhaseInducedBuoyancy/secondaryPhaseInducedBuoyancy.C
+$(otherForceModels)/gravityEff/gravityEff.C
 
 $(probeModels)/probeModel/probeModel.C
 $(probeModels)/probeModel/newProbeModel.C

src/lagrangian/cfdemParticle/Make/options

@@ -11,7 +11,6 @@ include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
 EXE_INC = \
      $(PFLAGS) \
      $(PINC) \
-    -I$(CFDEM_OFVERSION_DIR) \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \