Merge pull request #154 from ParticulateFlow/release

Release 24.01

README.md

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

applications/solvers/cfdemSolverIB/Make/options

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

applications/solvers/cfdemSolverIBContinuousForcing/Make/options

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

applications/solvers/cfdemSolverMultiphase/Make/options

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

applications/solvers/cfdemSolverMultiphaseScalar/Make/options

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

applications/solvers/cfdemSolverPimple/Make/files

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

applications/solvers/cfdemSolverPimple/Make/options

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

applications/solvers/cfdemSolverPimple/UEqn.H

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

applications/solvers/cfdemSolverPimple/cfdemSolverPimple.C

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

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 \

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/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,9 +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;
+    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/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

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/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 \

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

@@ -59,4 +59,6 @@
     thermo.correct();
 
     Info<< "T max/min : " << max(T).value() << " " << min(T).value() << endl;
+
+    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 \

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 \

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

@@ -316,7 +316,7 @@ void readDump(std::string filename, labelList &indices, scalarList &radii, vecto
 
     const label leadingLines = 9;
     label lineCounter = 0;
-    label partIndex;
+    label partIndex = 0;
     scalar r = 1.0, x = 0.0, y = 0.0, z = 0.0;
 
     indices.clear();

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,6 +122,7 @@ isotropicFluctuations,
 particleDeformation,
 "pdCorrelation"_forceModel_pdCorrelation.html,
 potentialRelaxation,
+"staticPressure"_forceModel_staticPressure.html,
 "surfaceTensionForce"_forceModel_surfaceTensionForce.html,
 terminalVelocity,
 "transferFluidProperties"_forceModel_transferFluidProperties.html,
@@ -199,6 +200,7 @@ each DEM particle.
 
 expParticleForces,
 gravity,
+"gravityEff"_otherForceModel_gravityEff.html,
 secondaryPhaseInducedBuoyancy,
 weightSecondaryPhase :tb(c=2,ea=c)
 

doc/CFDEMcoupling_solvers.txt

@@ -13,6 +13,7 @@ This section lists all CFDEMcoupling solvers alphabetically.
 "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/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_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/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
+

etc/OFversion/OFversion.H

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

etc/bashrc

@@ -17,7 +17,7 @@
 #------------------------------------------------------------------------------
 
 export CFDEM_PROJECT=CFDEM
-export CFDEM_VERSION=23.02
+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 23.02
+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/solver-list.txt

@@ -7,6 +7,7 @@ rcfdemSolverBase/dir
 rtfmSolverSpecies/dir
 cfdemSolverPisoMS/dir
 cfdemSolverPiso/dir
+cfdemSolverPimple/dir
 cfdemSolverRhoPimple/dir
 cfdemSolverIB/dir
 cfdemSolverPisoScalar/dir

src/lagrangian/cfdemParticle/Make/files

@@ -109,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
@@ -126,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 \

src/lagrangian/cfdemParticle/cfdTools/versionInfo.H

@@ -34,8 +34,8 @@ Description
 #ifndef versionInfo_H
 #define versionInfo_H
 
-word CFDEMversion="PFM 23.02";
-word compatibleLIGGGHTSversion="PFM 23.02";
+word CFDEMversion="PFM 24.01";
+word compatibleLIGGGHTSversion="PFM 24.01";
 word OFversion="6";
 
 Info << "\nCFDEMcoupling version: " << CFDEMversion <<  endl;

src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.H

@@ -43,7 +43,6 @@ SourceFiles
 #define cfdemCloud_H
 
 // choose version
-#include "OFversion.H"
 #include <vector>
 #include <typeinfo>
 

src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunn/heatTransferGunn.C

@@ -51,6 +51,10 @@ heatTransferGunn::heatTransferGunn
     implicit_(propsDict_.lookupOrDefault<bool>("implicit",true)),
     calcTotalHeatFlux_(propsDict_.lookupOrDefault<bool>("calcTotalHeatFlux",true)),
     initPartTemp_(propsDict_.lookupOrDefault<bool>("initPartTemp",false)),
+    scaleNuCellsName_(propsDict_.lookupOrDefault<word>("scaleNuCellsName","all")),
+    scaleNuCells_(),
+    allScaleNuCells_(false),
+    scaleNu_(false),
     Tmin_(propsDict_.lookupOrDefault<scalar>("Tmin",0.0)),
     Tmax_(propsDict_.lookupOrDefault<scalar>("Tmax",1e6)),
     totalHeatFlux_(0.0),
@@ -165,8 +169,14 @@ heatTransferGunn::heatTransferGunn
 
     if (propsDict_.found("NusseltScalingFactor"))
     {
+        scaleNu_ = true;
         NusseltScalingFactor_=readScalar(propsDict_.lookup ("NusseltScalingFactor"));
-        Info << "NusseltScalingFactor set to: " << NusseltScalingFactor_ << endl;
+        if(scaleNuCellsName_ != "all")
+        {
+            scaleNuCells_.set(new cellSet(particleCloud_.mesh(),scaleNuCellsName_));
+        }
+        else allScaleNuCells_ = true;
+        Info << "NusseltScalingFactor set to: " << NusseltScalingFactor_ << " in cellSet " << scaleNuCellsName_ << endl;
     }
 
     if (propsDict_.found("maxSource"))
@@ -237,6 +247,11 @@ heatTransferGunn::~heatTransferGunn()
 }
 
 // * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
+bool heatTransferGunn::scaleNuCell(label cell) const
+{
+    if (allScaleNuCells_) return true;
+    else return scaleNuCells_()[cell];
+}
 
 // * * * * * * * * * * * * * * * * Member Fct  * * * * * * * * * * * * * * * //
 
@@ -359,7 +374,14 @@ void heatTransferGunn::calcEnergyContribution()
                     Pr = max(SMALL, Cp * muf / kf0);
                     Nup = Nusselt(voidfraction, Rep, Pr);
                 }
-                Nup *= NusseltScalingFactor_;
+
+                if (scaleNu_)
+                {
+                    if (scaleNuCell(cellI))
+                    {
+                        Nup *= NusseltScalingFactor_;
+                    }
+                }
 
                 Tsum += partTemp_[index][0];
                 Nsum += 1.0;

src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunn/heatTransferGunn.H

@@ -28,6 +28,8 @@ License
 #include "fvCFD.H"
 #include "cfdemCloudEnergy.H"
 #include "energyModel.H"
+#include "autoPtr.H"
+#include "cellSet.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -60,6 +62,14 @@ protected:
 
     bool initPartTemp_;
 
+    word scaleNuCellsName_;
+
+    autoPtr<cellSet> scaleNuCells_;
+
+    bool allScaleNuCells_;
+
+    bool scaleNu_;
+
     scalar Tmin_;
 
     scalar Tmax_;
@@ -128,6 +138,8 @@ protected:
 
     scalar Nusselt(scalar, scalar, scalar) const;
 
+    bool scaleNuCell(label) const;
+
     virtual void giveData();
 
     virtual void heatFlux(label, scalar, scalar, scalar, scalar cg3 = 1.0);

src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferInterGrain/heatTransferInterGrain.C

@@ -295,6 +295,11 @@ void heatTransferInterGrain::calcEnergyContribution()
 
 void heatTransferInterGrain::calcPartEffThermCond()
 {
+    // clarification wrt to notation: Qian et al. define effective conductivities such that they fit into the
+    //     form \alpha_p * k * \nabla T_p, hence the multiplication with volFrac below.
+    //     Carson et al. define total bed conductivities which would not have an \alpha_p factor in the heat flux
+    //     term, therefore there is no such multiplication below.
+    //     To write down both contributions in one equation, one has to redefine the conductivity of Carson et al.
     calcPartThermCond();
     if (radiativeHeatTransfer_) calcPartThermRad();
 

src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferRanzMarshall/heatTransferRanzMarshall.C

@@ -50,6 +50,10 @@ heatTransferRanzMarshall::heatTransferRanzMarshall
     implicit_(propsDict_.lookupOrDefault<bool>("implicit",true)),
     calcTotalHeatFlux_(propsDict_.lookupOrDefault<bool>("calcTotalHeatFlux",true)),
     initPartTemp_(propsDict_.lookupOrDefault<bool>("initPartTemp",false)),
+    scaleNuCellsName_(propsDict_.lookupOrDefault<word>("scaleNuCellsName","all")),
+    scaleNuCells_(),
+    allScaleNuCells_(false),
+    scaleNu_(false),
     Tmin_(propsDict_.lookupOrDefault<scalar>("Tmin",0.0)),
     Tmax_(propsDict_.lookupOrDefault<scalar>("Tmax",1e6)),
     totalHeatFlux_(0.0),
@@ -165,8 +169,14 @@ heatTransferRanzMarshall::heatTransferRanzMarshall
 
     if (propsDict_.found("NusseltScalingFactor"))
     {
+        scaleNu_ = true;
         NusseltScalingFactor_=readScalar(propsDict_.lookup ("NusseltScalingFactor"));
-        Info << "NusseltScalingFactor set to: " << NusseltScalingFactor_ << endl;
+        if(scaleNuCellsName_ != "all")
+        {
+            scaleNuCells_.set(new cellSet(particleCloud_.mesh(),scaleNuCellsName_));
+        }
+        else allScaleNuCells_ = true;
+        Info << "NusseltScalingFactor set to: " << NusseltScalingFactor_ << " in cellSet " << scaleNuCellsName_ << endl;
     }
 
     if (NusseltConstParameter_ < 0.6 || NusseltConstParameter_ > 1.8)
@@ -246,6 +256,11 @@ heatTransferRanzMarshall::~heatTransferRanzMarshall()
 }
 
 // * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
+bool heatTransferRanzMarshall::scaleNuCell(label cell) const
+{
+    if (allScaleNuCells_) return true;
+    else return scaleNuCells_()[cell];
+}
 
 // * * * * * * * * * * * * * * * * Member Fct  * * * * * * * * * * * * * * * //
 
@@ -369,7 +384,14 @@ void heatTransferRanzMarshall::calcEnergyContribution()
                     Pr = max(SMALL, Cp * muf / kf0);
                     Nup = Nusselt(voidfraction, Rep, Pr);
                 }
-                Nup *= NusseltScalingFactor_;
+
+                if (scaleNu_)
+                {
+                    if (scaleNuCell(cellI))
+                    {
+                        Nup *= NusseltScalingFactor_;
+                    }
+                }
 
                 Tsum += partTemp_[index][0];
                 Nsum += 1.0;

src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferRanzMarshall/heatTransferRanzMarshall.H

@@ -28,6 +28,8 @@ License
 #include "fvCFD.H"
 #include "cfdemCloudEnergy.H"
 #include "energyModel.H"
+#include "autoPtr.H"
+#include "cellSet.H"
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -59,6 +61,14 @@ protected:
 
     bool initPartTemp_;
 
+    word scaleNuCellsName_;
+
+    autoPtr<cellSet> scaleNuCells_;
+
+    bool allScaleNuCells_;
+
+    bool scaleNu_;
+
     scalar Tmin_;
 
     scalar Tmax_;
@@ -129,6 +139,8 @@ protected:
 
     scalar Nusselt(scalar, scalar, scalar) const;
 
+    bool scaleNuCell(label) const;
+
     virtual void giveData();
 
     virtual void heatFlux(label, scalar, scalar, scalar, scalar cg3 = 1.0);

src/lagrangian/cfdemParticle/subModels/energyModel/reactionHeat/reactionHeat.C

@@ -22,6 +22,7 @@ License
 #include "reactionHeat.H"
 #include "addToRunTimeSelectionTable.H"
 #include "dataExchangeModel.H"
+#include "smoothingModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -48,6 +49,7 @@ reactionHeat::reactionHeat
     interpolation_(propsDict_.lookupOrDefault<bool>("interpolation",false)),
     verbose_(propsDict_.lookupOrDefault<bool>("verbose",false)),
     execution_(true),
+    smoothen_(propsDict_.lookupOrDefault<bool>("smoothen",false)),
     mesh_(sm.mesh()),
     maxSource_(1e30),
     reactionHeatName_(propsDict_.lookupOrDefault<word>("reactionHeatName","reactionHeat")),
@@ -123,6 +125,10 @@ void reactionHeat::calcEnergyContribution()
     );
 
     reactionHeatField_.primitiveFieldRef() /= (reactionHeatField_.mesh().V() * Nevery_ * couplingTimestep_);
+    if (smoothen_)
+    {
+        particleCloud_.smoothingM().smoothen(reactionHeatField_);
+    }
 
     forAll(reactionHeatField_,cellI)
     {
@@ -134,11 +140,8 @@ void reactionHeat::calcEnergyContribution()
         }
     }
 
-    if (verbose_)
-    {
-        Info << "reaction heat per unit time = "
-                     << gSum(reactionHeatField_*1.0*reactionHeatField_.mesh().V()) << endl;
-    }
+    Info << "reaction heat per unit time = "
+                 << gSum(reactionHeatField_*1.0*reactionHeatField_.mesh().V()) << endl;
 
     reactionHeatField_.correctBoundaryConditions();
 }

src/lagrangian/cfdemParticle/subModels/energyModel/reactionHeat/reactionHeat.H

@@ -52,6 +52,8 @@ protected:
 
     bool execution_;
 
+    bool smoothen_;
+
     const fvMesh& mesh_;
 
     scalar maxSource_;

src/lagrangian/cfdemParticle/subModels/forceModel/staticPressure/staticPressure.C

@@ -0,0 +1,100 @@
+/*---------------------------------------------------------------------------*\
+    CFDEMcoupling - Open Source CFD-DEM coupling
+
+    Copyright (C) 2023  Behrad Esgandari, JKU Linz, Austria
+-------------------------------------------------------------------------------
+License
+    This file is part of CFDEMcoupling.
+
+    CFDEMcoupling is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 3 of the License, or (at your
+    option) any later version.
+
+    CFDEMcoupling is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with CFDEMcoupling; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+\*---------------------------------------------------------------------------*/
+
+#include "error.H"
+#include "staticPressure.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(staticPressure, 0);
+
+addToRunTimeSelectionTable
+(
+    forceModel,
+    staticPressure,
+    dictionary
+);
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+// Construct from components
+staticPressure::staticPressure
+(
+    const dictionary& dict,
+    cfdemCloud& sm
+)
+:
+    forceModel(dict,sm),
+    propsDict_(dict.subDict(typeName + "Props")),
+    rhoGas_(readScalar(propsDict_.lookup ("rhoGas"))),
+    g_(propsDict_.lookup ("g_DEM")),
+    voidfractionFieldName_(propsDict_.lookup("voidfractionFieldName")),
+    voidfraction_(sm.mesh().lookupObject<volScalarField> (voidfractionFieldName_)),
+    solidFraction_(readScalar(propsDict_.lookup ("DomainSolidVolumeFraction")))
+{
+    // init force sub model
+    setForceSubModels(propsDict_);
+    // define switches which can be read from dict
+    forceSubM(0).setSwitchesList(SW_TREAT_FORCE_EXPLICIT,true); // activate treatExplicit switch
+    forceSubM(0).setSwitchesList(SW_TREAT_FORCE_DEM,true); // activate treatForceDEM switch
+    forceSubM(0).setSwitchesList(SW_INTERPOLATION,true); // activate search for interpolate switch
+}
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+staticPressure::~staticPressure()
+{}
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+void staticPressure::setForce() const
+{
+    #include "setupProbeModel.H"
+    label cellI;
+    scalar rhoMix_ = solidFraction_*particleCloud_.particleDensity(0) + (1.0-solidFraction_)*rhoGas_;
+    vector force;
+
+    for(int index = 0; index < particleCloud_.numberOfParticles(); ++index)
+    {
+        force = vector::zero;
+        cellI = particleCloud_.cellIDs()[index][0];
+        if (cellI > -1) // particle Found
+        {
+            scalar Vs = particleCloud_.particleVolume(index);
+            // set force on particle
+            force = -Vs * rhoMix_ * g_;
+        }
+        // write particle based data to global array
+        forceSubM(0).partToArray(index,force,vector::zero);
+     }
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/lagrangian/cfdemParticle/subModels/forceModel/staticPressure/staticPressure.H

@@ -0,0 +1,103 @@
+/*---------------------------------------------------------------------------*\
+    CFDEMcoupling - Open Source CFD-DEM coupling
+
+    Copyright (C) 2023  Behrad Esgandari, JKU Linz, Austria
+-------------------------------------------------------------------------------
+License
+    This file is part of CFDEMcoupling.
+
+    CFDEMcoupling is free software; you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by the
+    Free Software Foundation; either version 3 of the License, or (at your
+    option) any later version.
+
+    CFDEMcoupling is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with CFDEMcoupling; if not, write to the Free Software Foundation,
+    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+
+Description
+    A gravitational force applied on the particles just in case of
+    "Periodic Box simulations"
+
+Class
+    staticPressure
+
+SourceFiles
+    staticPressure.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef staticPressure_H
+#define staticPressure_H
+
+#include "forceModel.H"
+#include "interpolationCellPoint.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class staticPressure Declaration
+\*---------------------------------------------------------------------------*/
+
+class staticPressure
+:
+    public forceModel
+{
+private:
+
+    dictionary propsDict_;
+
+    word rhoGFieldName_;
+
+    scalar rhoGas_;
+
+    vector g_;
+
+    word voidfractionFieldName_;
+
+    const volScalarField& voidfraction_;
+
+    scalar solidFraction_;
+
+public:
+
+    //- Runtime type information
+    TypeName("staticPressure");
+
+
+    // Constructors
+
+        //- Construct from components
+        staticPressure
+        (
+            const dictionary& dict,
+            cfdemCloud& sm
+        );
+
+    // Destructor
+
+        ~staticPressure();
+
+
+    // Member Functions
+        void setForce() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/cfdemParticle/subModels/liggghtsCommandModel/writeLiggghts/writeLiggghts.C

@@ -138,6 +138,19 @@ const char* writeLiggghts::command(int commandLine)
 bool writeLiggghts::runCommand(int couplingStep)
 {
     if (!overwrite_) strCommand_ = addTimeStamp(command_);
+
+    // update time settings if restart file is only written at last step and
+    // end time has changed
+    if
+    (
+        runLast_
+        && particleCloud_.mesh().time().runTimeModifiable()
+        && particleCloud_.mesh().time().endTime().value() != endTime_
+    )
+    {
+        checkTimeSettings(dict_);
+    }
+
     return runThisCommand(couplingStep);
 }
 

src/lagrangian/cfdemParticle/subModels/otherForceModel/gravityEff/gravityEff.C

@@ -0,0 +1,106 @@
+/*---------------------------------------------------------------------------*\
+    CFDEMcoupling - Open Source CFD-DEM coupling
+
+    Copyright (C) 2023  Behrad Esgandari, JKU Linz, Austria
+-------------------------------------------------------------------------------
+License
+    This is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    This code is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this code.  If not, see <http://www.gnu.org/licenses/>.
+
+\*---------------------------------------------------------------------------*/
+
+#include "error.H"
+#include "gravityEff.H"
+#include "mathExtra.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+defineTypeNameAndDebug(gravityEff, 0);
+
+addToRunTimeSelectionTable(otherForceModel,gravityEff, dictionary);
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+// Construct from components
+gravityEff::gravityEff
+(
+    const dictionary& dict,
+    cfdemCloud& sm
+)
+:
+    otherForceModel(dict,sm),
+    propsDict_(dict.subDict(typeName + "Props")),
+    voidfractionFieldName_(propsDict_.lookupOrDefault<word>("voidfractionFieldName","voidfraction")),
+    voidfraction_(sm.mesh().lookupObject<volScalarField> (voidfractionFieldName_)),
+    rhoGFieldName_(propsDict_.lookupOrDefault<word>("rhoGFieldName","rho")),
+    rhoG_(sm.mesh().lookupObject<volScalarField> (rhoGFieldName_)),
+    gravityFieldName_(propsDict_.lookupOrDefault<word>("gravityFieldName","g")),
+    g_(sm.mesh().lookupObject<uniformDimensionedVectorField> (gravityFieldName_)),
+    rhoPart_("rhoPart",dimensionSet(1,-3,0,0,0),0.0),
+    unity_(sm.mesh().lookupObject<volScalarField> ("unity"))
+{
+    if (propsDict_.found("rhoPart"))
+        rhoPart_.value()=readScalar(propsDict_.lookup ("rhoPart"));
+    else
+        FatalError <<"Please specify density of solid phase.\n" << abort(FatalError);
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+gravityEff::~gravityEff()
+{}
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+tmp<volVectorField> gravityEff::exportForceField()
+{
+    tmp<volVectorField> tsource
+    (
+        new volVectorField
+        (
+            IOobject
+            (
+                "grav",
+                particleCloud_.mesh().time().timeName(),
+                particleCloud_.mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            particleCloud_.mesh(),
+            dimensionedVector
+            (
+                "zero",dimensionSet(1, -2, -2, 0, 0),vector::zero
+            )
+        )
+    );
+
+    if (propsDict_.found("rhoPart"))
+    {
+        //Mixture (gas and solid) density based on whole domain
+        dimensionedScalar rhoMix_ = fvc::domainIntegrate((1.0-voidfraction_) * rhoPart_ + voidfraction_ * rhoG_) / fvc::domainIntegrate(unity_);
+        volVectorField& source = tsource.ref();
+        source = (rhoG_ - rhoMix_) * voidfraction_ * g_;
+    }
+
+    return tsource;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/lagrangian/cfdemParticle/subModels/otherForceModel/gravityEff/gravityEff.H

@@ -0,0 +1,105 @@
+/*---------------------------------------------------------------------------*\
+    CFDEMcoupling - Open Source CFD-DEM coupling
+
+    Copyright (C) 2023  Behrad Esgandari, JKU Linz, Austria
+-------------------------------------------------------------------------------
+License
+    This is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    This code is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this code.  If not, see <http://www.gnu.org/licenses/>.
+
+Description
+    Gravitaional source term that is necessary to derive the flow just in case
+    of "Periodic Box simulations"
+
+Class
+    gravityEff
+
+SourceFiles
+    gravityEff.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef gravityEff_H
+#define gravityEff_H
+
+#include "otherForceModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                         Class gravityEff Declaration
+\*---------------------------------------------------------------------------*/
+
+class gravityEff
+:
+    public otherForceModel
+{
+
+protected:
+
+    // Protected data
+        dictionary propsDict_;
+
+        word voidfractionFieldName_;
+
+        const volScalarField& voidfraction_;
+
+        word rhoGFieldName_;
+
+        const volScalarField& rhoG_;
+
+        word gravityFieldName_;
+
+        const uniformDimensionedVectorField& g_;// ref to gravity
+
+        dimensionedScalar rhoPart_;
+
+        const volScalarField& unity_;
+public:
+
+    //- Runtime type information
+    TypeName("gravityEff");
+
+    // Constructors
+
+        //- Construct from components
+        gravityEff
+        (
+            const dictionary& dict,
+            cfdemCloud& sm
+        );
+
+
+    // Destructor
+        virtual ~gravityEff();
+
+
+
+    // Member Functions
+        tmp<volVectorField> exportForceField();
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/cfdemParticle/subModels/smoothingModel/smoothingModel/smoothingModel.C

@@ -69,7 +69,8 @@ smoothingModel::smoothingModel
             IOobject::NO_WRITE
         ),
         particleCloud_.mesh(),
-        dimensionedVector("zero", dimensionSet(0,0,0,0,0), vector::zero)
+        dimensionedVector("zero", dimensionSet(0,0,0,0,0), vector::zero),
+        "zeroGradient"
     ),
     sSmoothField_
     (
@@ -82,7 +83,8 @@ smoothingModel::smoothingModel
             IOobject::NO_WRITE
         ),
         particleCloud_.mesh(),
-        dimensionedScalar("zero", dimensionSet(0,0,0,0,0), 0)
+        dimensionedScalar("zero", dimensionSet(0,0,0,0,0), 0),
+        "zeroGradient"
     )
 {}
 

src/lagrangian/cfdemParticle/subModels/voidFractionModel/trilinearVoidFraction/trilinearVoidFraction.C

@@ -76,7 +76,14 @@ trilinearVoidFraction::trilinearVoidFraction
     checkWeightNporosity(propsDict_);
     if (porosity() != 1.) FatalError << "porosity not used in trilinearVoidFraction" << abort(FatalError);
 
-    Warning << "trilinearVoidFraction model is not yet complete and does not work near boundaries" << endl;
+    Warning << "Trilinear void fraction model is restricted to cubic cells and does not work near boundaries" << endl;
+
+    Pout<< "Trilinear void fraction model assumes " << nCellXYZ_ << " cells in x, y, z direction" << endl;
+
+    if (particleCloud_.mesh().nCells() != nCellXYZ_[0]*nCellXYZ_[1]*nCellXYZ_[2])
+    {
+        FatalError<< "Trilinear void fraction model found different number of cells than expected" << exit(FatalError);
+    }
 }
 
 

src/lagrangian/cfdemParticleComp/Make/files

@@ -99,6 +99,7 @@ $(forceModels)/granKineticEnergy/granKineticEnergy.C
 $(forceModels)/particleDeformation/particleDeformation.C
 $(forceModels)/turbulentDispersion/turbulentDispersion.C
 $(forceModels)/turbulentVelocityFluctuations/turbulentVelocityFluctuations.C
+$(forceModels)/staticPressure/staticPressure.C
 
 $(forceSubModels)/forceSubModel/newForceSubModel.C
 $(forceSubModels)/forceSubModel/forceSubModel.C
@@ -112,6 +113,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/cfdemParticleComp/Make/options

@@ -14,7 +14,6 @@ PFLAGS+= -Dcompre
 EXE_INC = \
      $(PFLAGS) \
      $(PINC) \
-    -I$(CFDEM_OFVERSION_DIR) \
     -I../cfdemParticle/lnInclude \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/transportModels \

tutorials/cfdemSolverPimple/PeriodicBox/Allrun.sh

@@ -0,0 +1,31 @@
+#!/bin/bash
+#------------------------------------------------------------------------------
+# allrun script for periodic box test case
+# run periodic box
+# Behrad Esgandari - August 2023
+#------------------------------------------------------------------------------
+
+#- define variables
+casePath="$(dirname "$(readlink -f ${BASH_SOURCE[0]})")"
+
+#- include functions
+source $CFDEM_PROJECT_DIR/etc/functions.sh
+
+# check if mesh was built
+if [ -f "$casePath/CFD/constant/polyMesh/points" ]; then
+    echo "mesh was built before - using old mesh"
+else
+    echo "mesh needs to be built"
+    cd $casePath/CFD
+    blockMesh
+fi
+
+if [ -f "$casePath/DEM/post/restart/liggghts.restart" ];  then
+    echo "LIGGGHTS init was run before - using existing restart file"
+else
+    #- run serial DEM
+    $casePath/parDEMrun.sh
+fi
+
+#- run parallel CFD-DEM
+bash $casePath/parCFDDEMrun.sh

tutorials/cfdemSolverPimple/PeriodicBox/CFD/0/Ksl

@@ -0,0 +1,54 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  6
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      Ksl;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -3 -1 0 0 0 0];
+
+internalField   uniform 0;
+
+boundaryField
+{
+    bottom
+    {
+        type            cyclic;
+    }
+
+    top
+    {
+        type            cyclic;
+    }
+
+    lr1
+    {
+        type            cyclic;
+    }
+
+    lr2
+    {
+        type            cyclic;
+    }
+
+    vh1
+    {
+        type            cyclic;
+    }
+
+    vh2
+    {
+        type            cyclic;
+    }
+}
+
+// ************************************************************************* //

tutorials/cfdemSolverPimple/PeriodicBox/CFD/0/T

@@ -0,0 +1,54 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  6
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      T;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 1 0 0 0];
+
+internalField   uniform 293;
+
+boundaryField
+{
+    bottom
+    {
+        type            cyclic;
+    }
+
+    top
+    {
+        type            cyclic;
+    }
+
+    lr1
+    {
+        type            cyclic;
+    }
+
+    lr2
+    {
+        type            cyclic;
+    }
+
+    vh1
+    {
+        type            cyclic;
+    }
+
+    vh2
+    {
+        type            cyclic;
+    }
+}
+
+// ************************************************************************* //

tutorials/cfdemSolverPimple/PeriodicBox/CFD/0/U

@@ -0,0 +1,54 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  6
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volVectorField;
+    object      U;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    bottom
+    {
+        type            cyclic;
+    }
+
+    top
+    {
+        type            cyclic;
+    }
+
+    lr1
+    {
+        type            cyclic;
+    }
+
+    lr2
+    {
+        type            cyclic;
+    }
+
+    vh1
+    {
+        type            cyclic;
+    }
+
+    vh2
+    {
+        type            cyclic;
+    }
+}
+
+// ************************************************************************* //

tutorials/cfdemSolverPimple/PeriodicBox/CFD/0/Us

@@ -0,0 +1,54 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  6
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volVectorField;
+    object      Us;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 1 -1 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    bottom
+    {
+        type            cyclic;
+    }
+
+    top
+    {
+        type            cyclic;
+    }
+
+    lr1
+    {
+        type            cyclic;
+    }
+
+    lr2
+    {
+        type            cyclic;
+    }
+
+    vh1
+    {
+        type            cyclic;
+    }
+
+    vh2
+    {
+        type            cyclic;
+    }
+}
+
+// ************************************************************************* //

tutorials/cfdemSolverPimple/PeriodicBox/CFD/0/p

@@ -0,0 +1,54 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  6
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      p;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 2 -2 0 0 0 0];
+
+internalField   uniform 1e5;
+
+boundaryField
+{
+    bottom
+    {
+        type            cyclic;
+    }
+
+    top
+    {
+        type            cyclic;
+    }
+
+    lr1
+    {
+        type            cyclic;
+    }
+
+    lr2
+    {
+        type            cyclic;
+    }
+
+    vh1
+    {
+        type            cyclic;
+    }
+
+    vh2
+    {
+        type            cyclic;
+    }
+}
+
+// ************************************************************************* //

tutorials/cfdemSolverPimple/PeriodicBox/CFD/0/rho

@@ -0,0 +1,54 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  6
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      rho;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [1 -3 0 0 0 0 0];
+
+internalField   uniform 1.3;
+
+boundaryField
+{
+    bottom
+    {
+        type            cyclic;
+    }
+
+    top
+    {
+        type            cyclic;
+    }
+
+    lr1
+    {
+        type            cyclic;
+    }
+
+    lr2
+    {
+        type            cyclic;
+    }
+
+    vh1
+    {
+        type            cyclic;
+    }
+
+    vh2
+    {
+        type            cyclic;
+    }
+}
+
+// ************************************************************************* //

tutorials/cfdemSolverPimple/PeriodicBox/CFD/0/sSmoothField

@@ -0,0 +1,56 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  6
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    location    "0";
+    object      sSmoothField;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 1;
+
+boundaryField
+{
+    bottom
+    {
+        type            cyclic;
+    }
+
+    top
+    {
+        type            cyclic;
+    }
+
+    lr1
+    {
+        type            cyclic;
+    }
+
+    lr2
+    {
+        type            cyclic;
+    }
+
+    vh1
+    {
+        type            cyclic;
+    }
+
+    vh2
+    {
+        type            cyclic;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/cfdemSolverPimple/PeriodicBox/CFD/0/vSmoothField

@@ -0,0 +1,56 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  6
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volVectorField;
+    location    "0";
+    object      vSmoothField;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform (0 0 0);
+
+boundaryField
+{
+    bottom
+    {
+        type            cyclic;
+    }
+
+    top
+    {
+        type            cyclic;
+    }
+
+    lr1
+    {
+        type            cyclic;
+    }
+
+    lr2
+    {
+        type            cyclic;
+    }
+
+    vh1
+    {
+        type            cyclic;
+    }
+
+    vh2
+    {
+        type            cyclic;
+    }
+}
+
+
+// ************************************************************************* //

tutorials/cfdemSolverPimple/PeriodicBox/CFD/0/voidfraction

@@ -0,0 +1,54 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  6
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       volScalarField;
+    object      voidfraction;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dimensions      [0 0 0 0 0 0 0];
+
+internalField   uniform 1;
+
+boundaryField
+{
+    bottom
+    {
+        type            cyclic;
+    }
+
+    top
+    {
+        type            cyclic;
+    }
+
+    lr1
+    {
+        type            cyclic;
+    }
+
+    lr2
+    {
+        type            cyclic;
+    }
+
+    vh1
+    {
+        type            cyclic;
+    }
+
+    vh2
+    {
+        type            cyclic;
+    }
+}
+
+// ************************************************************************* //

tutorials/cfdemSolverPimple/PeriodicBox/CFD/constant/couplingProperties

@@ -0,0 +1,162 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  6
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      couplingProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+// sub-models & settings
+
+syncMode false;
+
+modelType "A";
+
+couplingInterval 50;
+
+voidFractionModel centre;
+
+locateModel engine;
+//locateModel engineSearchMany2Many;
+
+meshMotionModel noMeshMotion;
+
+regionModel allRegion;
+
+IOModel off;
+
+probeModel off;
+
+dataExchangeModel twoWayMPI;
+//dataExchangeModel twoWayOne2One;
+
+averagingModel dense;
+
+clockModel off;
+
+smoothingModel constDiffSmoothing;
+
+getParticleDensities true;
+
+forceModels
+(
+    BeetstraDrag
+    gradPForce
+    viscForce
+    staticPressure
+);
+
+otherForceModels
+(
+  gravityEff
+);
+
+momCoupleModels
+(
+    implicitCouple
+);
+
+turbulenceModelType "turbulenceProperties";
+
+
+// sub-model properties
+
+
+engineSearchMany2ManyProps
+{
+    engineProps
+    {
+        treeSearch true;
+    }
+}
+
+twoWayOne2OneProps
+{
+    liggghtsPath "../DEM/in.liggghts_run";
+    verbose true;
+}
+
+
+implicitCoupleProps
+{
+    velFieldName "U";
+    granVelFieldName "Us";
+    voidfractionFieldName "voidfraction";
+}
+
+
+gradPForceProps
+{
+    pFieldName "p";
+    voidfractionFieldName "voidfraction";
+    velocityFieldName "U";
+    interpolation true;
+    treatForceExplicit;
+}
+
+viscForceProps
+{
+    velocityFieldName "U";
+    interpolation true;
+    treatForceExplicit;
+}
+
+BeetstraDragProps
+{
+    velFieldName "U";
+    granVelFieldName "Us";
+    rhoParticle 1500.;
+    voidfractionFieldName "voidfraction";
+    interpolation true;
+    implForceDEM;
+}
+
+gravityEffProps
+{
+    voidfractionFieldName "voidfraction";
+    rhoPart 1500;
+}
+
+staticPressureProps
+{
+    rhoGas 1.3; // gas density
+    voidfractionFieldName "voidfraction";
+    DomainSolidVolumeFraction 0.05; // mass loading
+    g_DEM (0 0 -9.81); //gravity vector compatible with the gravity direction in "g"
+    treatForceExplicit;
+}
+
+constDiffSmoothingProps
+{
+    lowerLimit 0.01;
+    upperLimit 1.0;
+    smoothingLength 0.00025;
+}
+
+
+engineProps
+{
+    treeSearch true;
+}
+
+centreProps
+{
+    alphaMin 0.005;
+    scaleUpVol 1.0;
+}
+
+twoWayMPIProps
+{
+    liggghtsPath "../DEM/in.liggghts_run";
+}
+
+// ************************************************************************* //

tutorials/cfdemSolverPimple/PeriodicBox/CFD/constant/g

@@ -1,9 +1,9 @@
 /*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  1.6                                   |
-|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
-|    \\/     M anipulation  |                                                 |
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  6
+     \\/     M anipulation  |
 \*---------------------------------------------------------------------------*/
 FoamFile
 {

tutorials/cfdemSolverPimple/PeriodicBox/CFD/constant/liggghtsCommands

@@ -0,0 +1,31 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  6
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      liggghtsCommands;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+liggghtsCommandModels
+(
+    runLiggghts
+    writeLiggghts
+);
+
+writeLiggghtsProps
+{
+    writeLast off;
+    writeName "post/restart/liggghts.restartCFDEM";
+    overwrite on;
+}
+
+// ************************************************************************* //

tutorials/cfdemSolverPimple/PeriodicBox/CFD/constant/transportProperties

@@ -0,0 +1,23 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  6
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      transportProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+transportModel  Newtonian;
+
+nu              [0 2 -1 0 0 0 0] 1.8e-5;
+
+
+// ************************************************************************* //

tutorials/cfdemSolverPimple/PeriodicBox/CFD/constant/turbulenceProperties

@@ -0,0 +1,20 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  6
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "constant";
+    object      turbulenceProperties;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+simulationType laminar;
+
+// ************************************************************************* //

tutorials/cfdemSolverPimple/PeriodicBox/CFD/system/blockMeshDict

@@ -0,0 +1,85 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  6
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+convertToMeters 1.0;
+
+vertices
+(
+    (-0.004 -0.004 -0.016)
+    ( 0.004 -0.004 -0.016)
+    ( 0.004  0.004 -0.016)
+    (-0.004  0.004 -0.016)
+    (-0.004 -0.004  0.016)
+    ( 0.004 -0.004  0.016)
+    ( 0.004  0.004  0.016)
+    (-0.004  0.004  0.016)
+);
+
+blocks
+(
+    hex (0 1 2 3 4 5 6 7) (32 32 128)
+    simpleGrading (1 1 1)
+);
+
+edges
+(
+);
+
+boundary
+(
+    bottom
+    {
+        type cyclic;
+        neighbourPatch  top;
+        faces ((0 3 2 1));
+    }
+    top
+    {
+        type cyclic;
+        neighbourPatch  bottom;
+        faces ((4 5 6 7));
+    }
+    lr1
+    {
+        type cyclic;
+        neighbourPatch  lr2;
+        faces ((0 3 7 4));
+    }
+    lr2
+    {
+        type cyclic;
+        neighbourPatch  lr1;
+        faces ((1 5 6 2));
+    }
+    vh1
+    {
+        type cyclic;
+        neighbourPatch  vh2;
+        faces ((0 1 4 5));
+    }
+    vh2
+    {
+        type cyclic;
+        neighbourPatch  vh1;
+        faces ((2 6 7 3));
+    }
+);
+
+mergePatchPairs
+(
+);
+
+// ************************************************************************* //

tutorials/cfdemSolverPimple/PeriodicBox/CFD/system/controlDict

@@ -0,0 +1,57 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Version:  6
+     \\/     M anipulation  |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    location    "system";
+    object      controlDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+application     pisoFoam;
+
+startFrom       startTime;
+
+startTime       0;
+
+stopAt          endTime;
+
+endTime         3;
+
+deltaT          0.00005;
+
+writeControl    adjustableRunTime;
+
+writeInterval   0.01;
+
+purgeWrite      0;
+
+writeFormat     ascii;
+
+writePrecision  6;
+
+writeCompression off;
+
+timeFormat      general;
+
+timePrecision   6;
+
+runTimeModifiable true;
+
+adjustTimeStep  no;
+
+maxCo 0.8;
+
+
+functions
+{
+
+}
+// ************************************************************************* //