Merge pull request #120 from ParticulateFlow/feature/recurrence_chemistry

Feature/recurrence chemistry

applications/solvers/cfdemSolverPiso/UEqn.H

@@ -1,8 +1,11 @@
+particleCloud.otherForces(fOther);
+
 fvVectorMatrix UEqn
 (
     fvm::ddt(voidfraction,U) - fvm::Sp(fvc::ddt(voidfraction),U)
   + fvm::div(phi,U) - fvm::Sp(fvc::div(phi),U)
   + particleCloud.divVoidfractionTau(U, voidfraction)
+  - fOther/rho
   ==
     fvOptions(U)
   - fvm::Sp(Ksl/rho,U)

applications/solvers/cfdemSolverPiso/createFields.H

@@ -46,6 +46,21 @@
         //dimensionedScalar("0", dimensionSet(1, -3, -1, 0, 0), 1.0)
     );
 
+    Info<< "\nCreating body force field\n" << endl;
+    volVectorField fOther
+    (
+        IOobject
+        (
+            "fOther",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedVector("zero", dimensionSet(1,-2,-2,0,0,0,0), vector::zero)
+    );
+
     Info<< "\nReading voidfraction field voidfraction = (Vgas/Vparticle)\n" << endl;
     volScalarField voidfraction
     (

applications/solvers/cfdemSolverPisoFreeStreaming/Make/files

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

applications/solvers/cfdemSolverPisoFreeStreaming/Make/options

@@ -0,0 +1,28 @@
+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 \
+    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
+    -I$(LIB_SRC)/transportModels \
+    -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
+    -I$(FOAM_SOLVERS)/incompressible/pisoFoam \
+    -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/lnInclude \
+    -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/cfdTools \
+    -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/derived/cfdemCloudRec \
+    -I$(LIB_SRC)/sampling/lnInclude
+
+EXE_LIBS = \
+    -L$(CFDEM_LIB_DIR)\
+    -lturbulenceModels \
+    -lincompressibleTurbulenceModels \
+    -lincompressibleTransportModels \
+    -lfiniteVolume \
+    -lmeshTools \
+    -lfvOptions \
+    -lsampling \
+    -l$(CFDEM_LIB_NAME) \
+     $(CFDEM_ADD_LIB_PATHS) \
+     $(CFDEM_ADD_LIBS)

applications/solvers/cfdemSolverPisoFreeStreaming/cfdemSolverPisoFreeStreaming.C

@@ -0,0 +1,126 @@
+/*---------------------------------------------------------------------------*\
+    CFDEMcoupling - Open Source CFD-DEM coupling
+
+    CFDEMcoupling is part of the CFDEMproject
+    www.cfdem.com
+                                Christoph Goniva, christoph.goniva@cfdem.com
+                                Copyright (C) 1991-2009 OpenCFD Ltd.
+                                Copyright (C) 2009-2012 JKU, Linz
+                                Copyright (C) 2012-     DCS Computing GmbH,Linz
+-------------------------------------------------------------------------------
+License
+    This file is part of CFDEMcoupling.
+
+    CFDEMcoupling is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    CFDEMcoupling is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with CFDEMcoupling.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    cfdemSolverPisoFreeStreaming
+
+Description
+    Transient solver for incompressible flow.
+    Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
+    The code is an evolution of the solver pisoFoam in OpenFOAM(R) 1.6,
+    where additional functionality for CFD-DEM coupling is added.
+    the particles follow the fluid velocity
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "singlePhaseTransportModel.H"
+#include "turbulentTransportModel.H"
+#include "pisoControl.H"
+#include "fvOptions.H"
+
+#include "cfdemCloudRec.H"
+
+#include "cfdemCloud.H"
+#include "implicitCouple.H"
+#include "clockModel.H"
+#include "smoothingModel.H"
+#include "forceModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "postProcess.H"
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createControl.H"
+    #include "createFields.H"
+    #include "createFvOptions.H"
+    #include "initContinuityErrs.H"
+
+    cfdemCloudRec<cfdemCloud> particleCloud(mesh);
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+    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");
+
+        particleCloud.evolve(voidfraction,Us,U);
+
+
+        particleCloud.clockM().stop("Coupling");
+
+        particleCloud.clockM().start(26,"Flow");
+
+        if(particleCloud.solveFlow())
+        {
+            // Pressure-velocity PISO corrector
+            {
+                // Momentum predictor
+                 #include "UEqn.H"
+
+                // --- PISO loop
+
+                while (piso.correct())
+                {
+                    #include "pEqn.H"
+                }
+            }
+
+            laminarTransport.correct();
+            turbulence->correct();
+        }
+        else
+        {
+            Info << "skipping flow solution." << endl;
+        }
+
+        runTime.write();
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+
+        particleCloud.clockM().stop("Flow");
+        particleCloud.clockM().stop("Global");
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/cfdemSolverPisoFreeStreaming/createFields.H

@@ -0,0 +1,110 @@
+    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 voidfraction field voidfraction = (Vgas/Vparticle)\n" << endl;
+    volScalarField voidfraction
+    (
+        IOobject
+        (
+            "voidfraction",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< "\nCreating 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
+    );
+
+//===============================
+
+//#   include "createPhi.H"
+#ifndef createPhi_H
+#define createPhi_H
+Info<< "Reading/calculating face flux field phi\n" << endl;
+surfaceScalarField phi
+(
+    IOobject
+    (
+        "phi",
+        runTime.timeName(),
+         mesh,
+        IOobject::READ_IF_PRESENT,
+        IOobject::AUTO_WRITE
+    ),
+    linearInterpolate(U) & mesh.Sf()
+);
+#endif
+
+
+
+    label pRefCell = 0;
+    scalar pRefValue = 0.0;
+    setRefCell(p, mesh.solutionDict().subDict("PISO"), pRefCell, pRefValue);
+
+
+    singlePhaseTransportModel laminarTransport(U, phi);
+
+    autoPtr<incompressible::turbulenceModel> turbulence
+    (
+        incompressible::turbulenceModel::New(U, phi, laminarTransport)
+    );
+
+#include "createMRF.H"

applications/solvers/cfdemSolverPisoScalar/createFields.H

@@ -46,6 +46,21 @@
         //dimensionedScalar("0", dimensionSet(0, 0, -1, 0, 0), 1.0)
     );
 
+    Info<< "\nCreating body force field\n" << endl;
+    volVectorField fOther
+    (
+        IOobject
+        (
+            "fOther",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedVector("zero", dimensionSet(1,-2,-2,0,0,0,0), vector::zero)
+    );
+
     Info<< "\nReading voidfraction field voidfraction = (Vgas/Vparticle)\n" << endl;
     volScalarField voidfraction
     (

applications/solvers/cfdemSolverRhoPimple/EEqn.H

@@ -20,21 +20,23 @@
 
     Cpv = he.name() == "e" ? thermo.Cv() : thermo.Cp();
 
-    // correct source for the thermodynamic reference temperature
-    dimensionedScalar Tref("Tref", dimTemperature, T[0]-he[0]/(Cpv[0]+SMALL));
-    Qsource += QCoeff*Tref;
+// For implict T terms in the energy/enthalpy transport equation, use
+// (he_n+1 - he_n) / (T_n+1 - T_n) = Cpv to eliminate T_n+1 with he_n+1.
+// This formula is valid for ideal gases with e=e(T) and h=h(T). For
+// incompressible fluids, e=e(T) holds, too, but enthalpy would need correction
+// terms accounting for pressure variations.
 
     fvScalarMatrix EEqn
     (
         fvm::ddt(rhoeps, he) + fvm::div(phi, he)
       + addSource
-     // net heat transfer from particles to fluid
       - Qsource
+      - QCoeff*T
       - fvm::Sp(QCoeff/Cpv, he)
-     // thermal conduction of the fluid with effective conductivity
+      + QCoeff/Cpv*he
+      - fvc::laplacian(voidfraction*thCond,T)
       - fvm::laplacian(voidfraction*thCond/Cpv,he)
-      // + particle-fluid energy transfer due to work
-      // + fluid energy dissipation due to shearing
+      + fvc::laplacian(voidfraction*thCond/Cpv,he)
      ==
         fvOptions(rho, he)
     );

applications/solvers/rcfdemSolverBase/createFields.H

@@ -37,12 +37,20 @@
             "URec",
             runTime.timeName(),
             mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
+            IOobject::READ_IF_PRESENT,
+            IOobject::NO_WRITE
         ),
-        mesh
+        mesh,
+        dimensionedVector("URec", dimensionSet(0, 1, -1, 0, 0), vector::zero)
     );
 
+    bool updateURec = false;
+    if (URec.headerOk())
+    {
+        updateURec = true;
+        URec.writeOpt() = IOobject::AUTO_WRITE;
+    }
+
     volScalarField voidfractionRec
     (
         IOobject
@@ -50,12 +58,20 @@
             "voidfractionRec",
             runTime.timeName(),
             mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
+            IOobject::READ_IF_PRESENT,
+            IOobject::NO_WRITE
         ),
-        mesh
+        mesh,
+        dimensionedScalar("voidfractionRec", dimensionSet(0, 0, 0, 0, 0), 1.0)
     );
 
+    bool updateVoidfractionRec = false;
+    if (voidfractionRec.headerOk())
+    {
+        updateVoidfractionRec = true;
+        voidfractionRec.writeOpt() = IOobject::AUTO_WRITE;
+    }
+
     volVectorField UsRec
     (
         IOobject
@@ -63,12 +79,20 @@
             "UsRec",
             runTime.timeName(),
             mesh,
-            IOobject::MUST_READ,
+            IOobject::READ_IF_PRESENT,
             IOobject::AUTO_WRITE
         ),
-        mesh
+        mesh,
+        dimensionedVector("URec", dimensionSet(0, 1, -1, 0, 0), vector::zero)
     );
 
+    bool updateUsRec = false;
+    if (UsRec.headerOk())
+    {
+        updateUsRec = true;
+        UsRec.writeOpt() = IOobject::AUTO_WRITE;
+    }
+
     // calculated fields
     Info << "\nCreating fields subject to calculation\n" << endl;
     volScalarField voidfraction
@@ -78,7 +102,7 @@
             "voidfraction",
             runTime.timeName(),
             mesh,
-            IOobject::NO_READ,
+            IOobject::READ_IF_PRESENT,
             IOobject::AUTO_WRITE
         ),
         voidfractionRec
@@ -91,7 +115,7 @@
             "Us",
             runTime.timeName(),
             mesh,
-            IOobject::NO_READ,
+            IOobject::READ_IF_PRESENT,
             IOobject::AUTO_WRITE
         ),
         UsRec
@@ -111,11 +135,18 @@
             runTime.timeName(),
             mesh,
             IOobject::READ_IF_PRESENT,
-            IOobject::AUTO_WRITE
+            IOobject::NO_WRITE
         ),
         linearInterpolate(URec*voidfractionRec) & mesh.Sf()
     );
+
+    bool updatePhiRec = false;
+    if (phiRec.headerOk())
+    {
+        updatePhiRec = true;
+        phiRec.writeOpt() = IOobject::AUTO_WRITE;
         phiRec.write();
+    }
 
     singlePhaseTransportModel laminarTransport(URec, phiRec);
 
@@ -123,3 +154,40 @@
     (
         incompressible::turbulenceModel::New(URec, phiRec, laminarTransport)
     );
+
+    IOdictionary recDict
+    (
+        IOobject
+        (
+            "recProperties",
+            runTime.constant(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        )
+    );
+
+    word voidfractionFieldName(recDict.lookupOrDefault<word>("voidfractionFieldName","voidfraction"));
+    word UFieldName(recDict.lookupOrDefault<word>("UFieldName","U"));
+    word UsFieldName(recDict.lookupOrDefault<word>("UsFieldName","Us"));
+    word fluxFieldName(recDict.lookupOrDefault<word>("fluxFieldName","phi"));
+
+
+   // place to put weight functions
+    IOdictionary weightDict
+    (
+        IOobject
+        (
+            "weightDict",
+            runTime.constant(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::NO_WRITE
+        )
+    );
+    if (!weightDict.headerOk())
+    {
+        weightDict.add("weights",scalarList(1,1.0));
+    }
+    scalarList weights(weightDict.lookup("weights"));
+    Info << "database initial weights: " << weights << endl;

applications/solvers/rcfdemSolverBase/rcfdemSolverBase.C

@@ -43,6 +43,7 @@ Rules
 #include "cfdemCloudRec.H"
 #include "recBase.H"
 #include "recModel.H"
+#include "recPath.H"
 
 #include "cfdemCloud.H"
 #include "clockModel.H"
@@ -59,6 +60,8 @@ int main(int argc, char *argv[])
     #include "createFields.H"
     #include "createFvOptions.H"
 
+    #include "readGravitationalAcceleration.H"
+
     cfdemCloudRec<cfdemCloud> particleCloud(mesh);
     recBase recurrenceBase(mesh);
 
@@ -88,8 +91,9 @@ int main(int argc, char *argv[])
 
         if ( runTime.timeOutputValue() - startTime - (recTimeIndex+1)*recTimeStep + 1.0e-5 > 0.0 )
         {
+            Info << "updating recurrence fields at time " << runTime.timeName() << "with recTimeIndex = " << recTimeIndex << nl << endl;
             recurrenceBase.updateRecFields();
-            #include "readFields.H"
+            #include "updateFields.H"
             recTimeIndex++;
         }
 

applications/solvers/rcfdemSolverBase/updateFields.H

@@ -0,0 +1,29 @@
+scalarList wList(weightDict.lookupOrDefault("weights",scalarList(1,1.0)));
+
+recurrenceBase.recP().updateIntervalWeights(wList);
+
+if(recurrenceBase.recM().endOfPath())
+{
+    recurrenceBase.extendPath();
+}
+
+// update fields where necessary
+if (updateVoidfractionRec)
+{
+    recurrenceBase.recM().exportVolScalarField(voidfractionFieldName,voidfractionRec);
+}
+
+if (updateURec)
+{
+    recurrenceBase.recM().exportVolVectorField(UFieldName,URec);
+}
+
+if (updateUsRec)
+{
+    recurrenceBase.recM().exportVolVectorField(UsFieldName,UsRec);
+}
+
+if (updatePhiRec)
+{
+    recurrenceBase.recM().exportSurfaceScalarField(fluxFieldName,phiRec);
+}

applications/solvers/rcfdemSolverCoupledHeattransfer/TEqImp.H

@@ -10,6 +10,7 @@
     // main contribution due to gas expansion, not due to transport of kinetic energy
     // fvc::ddt(rhoeps, K) + fvc::div(phiRec, K)
 
+    // assuming constant Cv such that e = Cv * T
     fvScalarMatrix TEqn =
     (
         fvm::ddt(rhoeps, T)

applications/solvers/rcfdemSolverForcedTracers/Make/files

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

applications/solvers/rcfdemSolverForcedTracers/Make/options

@@ -0,0 +1,27 @@
+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 \
+    -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 \
+    -I$(CFDEM_SRC_DIR)/recurrence/lnInclude \
+    -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/derived/cfdemCloudRec \
+
+EXE_LIBS = \
+    -L$(CFDEM_LIB_DIR)\
+    -lrecurrence \
+    -lturbulenceModels \
+    -lincompressibleTurbulenceModels \
+    -lincompressibleTransportModels \
+    -lfiniteVolume \
+    -lmeshTools \
+    -lfvOptions \
+    -l$(CFDEM_LIB_NAME) \
+     $(CFDEM_ADD_LIB_PATHS) \
+$(CFDEM_ADD_LIBS)

applications/solvers/rcfdemSolverForcedTracers/createFields.H

@@ -0,0 +1,113 @@
+    // dummy fields
+    Info << "\nCreating dummy density field\n" << endl;
+
+    volScalarField rho
+    (
+        IOobject
+        (
+            "rho",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("rho", dimensionSet(1, -3, 0, 0, 0), 1.0)
+    );
+
+    // particle fields
+    Info << "\nCreating voidfraction and particle velocity fields\n" << endl;
+
+    volScalarField voidfraction
+    (
+        IOobject
+        (
+            "voidfraction",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    volVectorField Us
+    (
+        IOobject
+        (
+            "Us",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    // recurrence fields
+    Info << "\nCreating recurrence fields.\n" << endl;
+
+    volScalarField pRec
+    (
+        IOobject
+        (
+            "pRec",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("p", dimensionSet(1, 2, -2, 0, 0), 1.0)
+    );
+
+    volScalarField kRec
+    (
+        IOobject
+        (
+            "kRec",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("k", dimensionSet(0, 2, -2, 0, 0), 0.0)
+    );
+
+    volVectorField URec
+    (
+        IOobject
+        (
+            "URec",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+//===============================
+
+    Info << "Calculating face flux field phi\n" << endl;
+    surfaceScalarField phiRec
+    (
+        IOobject
+        (
+            "phiRec",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        linearInterpolate(URec*voidfraction) & mesh.Sf()
+    );
+    phiRec.write();
+
+    singlePhaseTransportModel laminarTransport(URec, phiRec);
+
+    autoPtr<incompressible::turbulenceModel> turbulence
+    (
+        incompressible::turbulenceModel::New(URec, phiRec, laminarTransport)
+    );

applications/solvers/rcfdemSolverForcedTracers/rcfdemSolverForcedTracers.C

@@ -0,0 +1,112 @@
+/*---------------------------------------------------------------------------*\
+    CFDEMcoupling academic - Open Source CFD-DEM coupling
+
+    Contributing authors:
+    Thomas Lichtenegger
+    Copyright (C) 2015- Johannes Kepler University, Linz
+-------------------------------------------------------------------------------
+License
+    This file is part of CFDEMcoupling academic.
+
+    CFDEMcoupling academic 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 academic 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 academic.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    rcfdemSolverForcedTracers
+
+Description
+    Moves tracers according to the activated force models on pressure and velocity
+    fields provided by a recurrence process
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "singlePhaseTransportModel.H"
+#include "turbulentTransportModel.H"
+#include "fvOptions.H"
+
+#include "recBase.H"
+#include "recModel.H"
+
+#include "cfdemCloud.H"
+#include "clockModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "postProcess.H"
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createControl.H"
+    #include "createFields.H"
+    #include "createFvOptions.H"
+
+    cfdemCloud particleCloud(mesh);
+    recBase recurrenceBase(mesh);
+
+    const IOdictionary& recProps = mesh.lookupObject<IOdictionary>("recProperties");
+    bool useRecP(recProps.lookupOrDefault<bool>("useRecP",false));
+    bool useRecK(recProps.lookupOrDefault<bool>("useRecK",false));
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info << "\nCalculating particle trajectories based on recurrence statistics\n" << endl;
+
+    label recTimeIndex = 0;
+    scalar recTimeStep = recurrenceBase.recM().recTimeStep();
+    scalar startTime = runTime.startTime().value();
+
+    while (runTime.run())
+    {
+        runTime++;
+
+        // do stuff (every lagrangian time step)
+        particleCloud.clockM().start(1,"Global");
+
+        Info << "Time = " << runTime.timeName() << nl << endl;
+
+        particleCloud.clockM().start(2,"Coupling");
+
+        particleCloud.evolve(voidfraction,Us,URec);
+
+        particleCloud.clockM().stop("Coupling");
+
+
+        if ( runTime.timeOutputValue() - startTime - (recTimeIndex+1)*recTimeStep + 1.0e-5 > 0.0 )
+        {
+            recurrenceBase.updateRecFields();
+            #include "updateFields.H"
+            recTimeIndex++;
+        }
+
+        particleCloud.clockM().start(27,"Output");
+        runTime.write();
+        particleCloud.clockM().stop("Output");
+
+        particleCloud.clockM().stop("Global");
+
+        Info << "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+             << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+             << nl << endl;
+
+    }
+
+    Info << "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/rcfdemSolverForcedTracers/updateFields.H

@@ -0,0 +1,13 @@
+recurrenceBase.recM().exportVolVectorField("U",URec);
+
+if (useRecP)
+{
+    recurrenceBase.recM().exportVolScalarField("p",pRec);
+}
+
+if (useRecK)
+{
+    recurrenceBase.recM().exportVolScalarField("k",kRec);
+// in case database contains the velocity variance instead of k, do
+// kRec *= 0.5;
+}

applications/solvers/rcfdemSolverRhoSteadyPimple/EEqn.H

@@ -22,19 +22,33 @@
 
     Cpv = he.name() == "e" ? thermo.Cv() : thermo.Cp();
 
+// For implict T terms in the energy/enthalpy transport equation, use
+// (he_n+1 - he_n) / (T_n+1 - T_n) = Cpv to eliminate T_n+1 with he_n+1.
+// This formula is valid for ideal gases with e=e(T) and h=h(T). For
+// incompressible fluids, e=e(T) holds, too, but enthalpy would need correction
+// terms accounting for pressure variations.
 
     fvScalarMatrix EEqn
     (
         fvm::div(phi, he)
       + addSource
       - Qsource
+      - QCoeff*T
       - fvm::Sp(QCoeff/Cpv, he)
- //     - fvm::laplacian(voidfractionRec*kf/Cpv,he)
+      + QCoeff/Cpv*he
+      - fvc::laplacian(voidfractionRec*thCond,T)
       - fvm::laplacian(voidfractionRec*thCond/Cpv,he)
+      + fvc::laplacian(voidfractionRec*thCond/Cpv,he)
      ==
         fvOptions(rho, he)
     );
 
+    if (transientEEqn)
+    {
+        EEqn += fvm::ddt(rho,voidfractionRec,he);
+    }
+
+
     EEqn.relax();
 
     fvOptions.constrain(EEqn);

applications/solvers/rcfdemSolverRhoSteadyPimple/createFields.H

@@ -168,6 +168,8 @@ Info<< "Reading thermophysical properties\n" << endl;
          linearInterpolate(rho*U*voidfraction) & mesh.Sf()
     );
 
+    bool transientEEqn(pimple.dict().lookupOrDefault<bool>("transientEEqn",false));
+
     dimensionedScalar rhoMax
     (
         dimensionedScalar::lookupOrDefault

applications/solvers/rcfdemSolverRhoSteadyPimpleChem/EEqn.H

@@ -0,0 +1,67 @@
+// contributions to internal energy equation can be found in
+// Crowe et al.: "Multiphase flows with droplets and particles", CRC Press 1998
+{
+    // dim he = J / kg
+    volScalarField& he = thermo.he();
+    particleCloud.energyContributions(Qsource);
+    particleCloud.energyCoefficients(QCoeff);
+
+    addSource =
+        (
+            he.name() == "e"
+          ?
+            fvc::div(phi, K) +
+            fvc::div
+            (
+                fvc::absolute(phi/fvc::interpolate(rho), voidfractionRec*U),
+                p,
+                "div(phiv,p)"
+            )
+          : fvc::div(phi, K)
+        );
+
+    Cpv = he.name() == "e" ? thermo.Cv() : thermo.Cp();
+
+// For implict T terms in the energy/enthalpy transport equation, use
+// (he_n+1 - he_n) / (T_n+1 - T_n) = Cpv to eliminate T_n+1 with he_n+1.
+// This formula is valid for ideal gases with e=e(T) and h=h(T). For
+// incompressible fluids, e=e(T) holds, too, but enthalpy would need correction
+// terms accounting for pressure variations.
+
+    fvScalarMatrix EEqn
+    (
+      fvm::div(phi, he)
+      + addSource
+      - Qsource
+      - QCoeff*T
+      - fvm::Sp(QCoeff/Cpv, he)
+      + QCoeff/Cpv*he
+      - fvc::laplacian(voidfractionRec*thCond,T)
+      - fvm::laplacian(voidfractionRec*thCond/Cpv,he)
+      + fvc::laplacian(voidfractionRec*thCond/Cpv,he)
+     ==
+        fvOptions(rho, he)
+    );
+
+    if (transientEEqn)
+    {
+        EEqn += fvm::ddt(rho,voidfractionRec,he);
+    }
+
+    EEqn.relax();
+
+    fvOptions.constrain(EEqn);
+
+    EEqn.solve();
+
+    fvOptions.correct(he);
+
+    thermo.correct();
+
+    Info<< "T max/min : " << max(T).value() << " " << min(T).value() << endl;
+
+
+    particleCloud.clockM().start(31,"energySolve");
+    particleCloud.solve();
+    particleCloud.clockM().stop("energySolve");
+}

applications/solvers/rcfdemSolverRhoSteadyPimpleChem/Make/files

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

applications/solvers/rcfdemSolverRhoSteadyPimpleChem/Make/options

@@ -0,0 +1,52 @@
+include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
+
+FOAM_VERSION_MAJOR := $(word 1,$(subst ., ,$(WM_PROJECT_VERSION)))
+PFLAGS+= -DOPENFOAM_VERSION_MAJOR=$(FOAM_VERSION_MAJOR)
+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 \
+    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/cfdTools \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(LIB_SRC)/fvOptions/lnInclude \
+    -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/lnInclude \
+    -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/cfdTools \
+    -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \
+    -I$(LIB_SRC)/regionModels/regionModel/lnInclude \
+    -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \
+    -I$(LIB_SRC)/ODE/lnInclude \
+    -I$(LIB_SRC)/combustionModels/lnInclude \
+    -I$(CFDEM_SRC_DIR)/recurrence/lnInclude \
+
+EXE_LIBS = \
+    -L$(CFDEM_LIB_DIR)\
+    -lrecurrence \
+    -lcompressibleTransportModels \
+    -lfluidThermophysicalModels \
+    -lspecie \
+    -lturbulenceModels \
+    -lcompressibleTurbulenceModels \
+    -lfiniteVolume \
+    -lmeshTools \
+    -lsampling \
+    -lfvOptions \
+    -l$(CFDEM_LIB_COMP_NAME) \
+     $(CFDEM_ADD_LIB_PATHS) \
+     $(CFDEM_ADD_LIBS) \
+    -lreactionThermophysicalModels \
+    -lchemistryModel \
+    -lradiationModels \
+    -lregionModels \
+    -lsurfaceFilmModels \
+    -lODE \
+    -lcombustionModels

applications/solvers/rcfdemSolverRhoSteadyPimpleChem/UEqn.H

@@ -0,0 +1,35 @@
+// Solve the Momentum equation
+particleCloud.otherForces(fOther);
+
+fvVectorMatrix UEqn
+(
+    fvm::div(phi, U)
+  + particleCloud.divVoidfractionTau(U, voidfractionRec)
+  + fvm::Sp(Ksl,U)
+  - fOther
+ ==
+    fvOptions(rho, U)
+);
+
+if (stepcounter%nEveryFlow==0)
+{
+    UEqn.relax();
+
+    fvOptions.constrain(UEqn);
+
+    if (modelType=="B" || modelType=="Bfull")
+    {
+        solve(UEqn == -fvc::grad(p)+ Ksl*UsRec);
+    }
+    else
+    {
+        solve(UEqn == -voidfractionRec*fvc::grad(p)+ Ksl*UsRec);
+    }
+
+
+    #include "limitU.H"
+
+    fvOptions.correct(U);
+
+    K = 0.5*magSqr(U);
+}

applications/solvers/rcfdemSolverRhoSteadyPimpleChem/YEqn.H

@@ -0,0 +1,81 @@
+particleCloud.clockM().start(29,"Y");
+
+tmp<fv::convectionScheme<scalar> > mvConvection
+(
+    fv::convectionScheme<scalar>::New
+    (
+        mesh,
+        fields,
+        phi,
+        mesh.divScheme("div(phi,Yi_h)")
+    )
+);
+
+{
+    combustion->correct();
+#if OPENFOAM_VERSION_MAJOR < 5
+    dQ = combustion->dQ();
+#else
+    Qdot = combustion->Qdot();
+#endif
+    label inertIndex = -1;
+    volScalarField Yt(0.0*Y[0]);
+
+    forAll(Y, i)
+    {
+        if (Y[i].name() == inertSpecie) inertIndex = i;
+        if (Y[i].name() != inertSpecie || propagateInertSpecie)
+        {
+            volScalarField& Yi = Y[i];
+
+            fvScalarMatrix YiEqn
+            (
+                mvConvection->fvmDiv(phi, Yi)
+              - fvm::laplacian(voidfractionRec*turbulence->muEff(), Yi)
+              ==
+                combustion->R(Yi)
+              + particleCloud.chemistryM(0).Smi(i)*p/p.prevIter()
+              + fvOptions(rho, Yi)
+            );
+
+            YiEqn.relax();
+
+            fvOptions.constrain(YiEqn);
+
+            YiEqn.solve(mesh.solver("Yi"));
+
+            Yi.relax();
+
+            fvOptions.correct(Yi);
+
+            Yi.max(0.0);
+            if (Y[i].name() != inertSpecie) Yt += Yi;
+        }
+    }
+
+    if (inertIndex!=-1)
+    {
+        Y[inertIndex].max(inertLowerBound);
+        Y[inertIndex].min(inertUpperBound);
+    }
+
+    if (propagateInertSpecie)
+    {
+        if (inertIndex!=-1) Yt /= (1-Y[inertIndex] + ROOTVSMALL);
+        forAll(Y,i)
+        {
+            if (i!=inertIndex)
+            {
+                volScalarField& Yi = Y[i];
+                Yi = Yi/(Yt+ROOTVSMALL);
+            }
+        }
+    }
+    else
+    {
+        Y[inertIndex] = scalar(1) - Yt;
+        Y[inertIndex].max(0.0);
+    }
+}
+
+particleCloud.clockM().stop("Y");

applications/solvers/rcfdemSolverRhoSteadyPimpleChem/createFieldRefs.H

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

applications/solvers/rcfdemSolverRhoSteadyPimpleChem/createFields.H

@@ -0,0 +1,404 @@
+Info<< "Reading thermophysical properties\n" << endl;
+
+#if OPENFOAM_VERSION_MAJOR < 6
+    Info<< "Creating combustion model\n" << endl;
+    autoPtr<combustionModels::rhoCombustionModel> combustion
+    (
+        combustionModels::rhoCombustionModel::New(mesh)
+    );
+    rhoReactionThermo& thermo = combustion->thermo();
+#else
+    Info<< "Reading thermophysical properties\n" << endl;
+    autoPtr<rhoReactionThermo> pThermo(rhoReactionThermo::New(mesh));
+    rhoReactionThermo& thermo = pThermo();
+#endif
+    thermo.validate(args.executable(), "h", "e");
+
+    basicSpecieMixture& composition = thermo.composition();
+    PtrList<volScalarField>& Y = composition.Y();
+
+    // read molecular weight
+#if OPENFOAM_VERSION_MAJOR < 6
+    volScalarField W(composition.W());
+#else
+    volScalarField W(thermo.W());
+#endif
+
+    bool propagateInertSpecie(thermo.lookupOrDefault<bool>("propagateInertSpecie",true));
+
+    const word inertSpecie(thermo.lookupOrDefault<word>("inertSpecie","none"));
+
+    const scalar inertLowerBound(thermo.lookupOrDefault<scalar>("inertLowerBound",0.0));
+
+    const scalar inertUpperBound(thermo.lookupOrDefault<scalar>("inertUpperBound",1.0));
+
+    if (!composition.contains(inertSpecie) && inertSpecie != "none")
+    {
+        FatalErrorIn(args.executable())
+            << "Specified inert specie '" << inertSpecie << "' not found in "
+            << "species list. Available species:" << composition.species()
+            << exit(FatalError);
+    }
+
+    volScalarField& p = thermo.p();
+    const volScalarField& T = thermo.T();
+    const volScalarField& psi = thermo.psi();
+
+    multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
+
+    forAll(Y, i)
+    {
+        fields.add(Y[i]);
+    }
+    fields.add(thermo.he());
+
+    Info<< "Reading field rho\n" << endl;
+    volScalarField rho
+    (
+        IOobject
+        (
+            "rho",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        thermo.rho()
+    );
+
+
+    Info<< "Reading field U\n" << endl;
+    volVectorField U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< "\nReading voidfraction field voidfraction = (Vgas/Vparticle)\n" << endl;
+    volScalarField voidfraction
+    (
+        IOobject
+        (
+            "voidfraction",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    volScalarField voidfractionRec
+    (
+        IOobject
+        (
+            "voidfractionRec",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        voidfraction
+    );
+
+    volScalarField addSource
+    (
+        IOobject
+        (
+           "addSource",
+           runTime.timeName(),
+           mesh,
+           IOobject::READ_IF_PRESENT,
+           IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("zero", dimensionSet(1,-1,-3,0,0,0,0), 0.0)
+    );
+
+    Info<< "\nCreating fluid-particle heat flux field\n" << endl;
+    volScalarField Qsource
+    (
+        IOobject
+        (
+            "Qsource",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("zero", dimensionSet(1,-1,-3,0,0,0,0), 0.0)
+    );
+
+    Info<< "\nCreating fluid-particle heat flux coefficient field\n" << endl;
+    volScalarField QCoeff
+    (
+        IOobject
+        (
+            "QCoeff",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("zero", dimensionSet(1,-1,-3,-1,0,0,0), 0.0)
+    );
+
+    Info<< "\nCreating thermal conductivity field\n" << endl;
+    volScalarField thCond
+    (
+        IOobject
+        (
+            "thCond",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("zero", dimensionSet(1,1,-3,-1,0,0,0), 0.0),
+        "zeroGradient"
+    );
+
+    Info<< "\nCreating heat capacity field\n" << endl;
+    volScalarField Cpv
+    (
+        IOobject
+        (
+            "Cpv",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("zero", dimensionSet(0,2,-2,-1,0,0,0), 0.0)
+    );
+
+    Info<< "\nCreating body force field\n" << endl;
+    volVectorField fOther
+    (
+        IOobject
+        (
+            "fOther",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedVector("zero", dimensionSet(1,-2,-2,0,0,0,0), vector::zero)
+    );
+
+    Info<< "Reading/calculating face flux field phi\n" << endl;
+    surfaceScalarField phi
+    (
+        IOobject
+        (
+            "phi",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+         ),
+         linearInterpolate(rho*U*voidfraction) & mesh.Sf()
+    );
+
+    bool transientEEqn(pimple.dict().lookupOrDefault<bool>("transientEEqn",false));
+
+    dimensionedScalar rhoMax
+    (
+        dimensionedScalar::lookupOrDefault
+        (
+            "rhoMax",
+            pimple.dict(),
+            dimDensity,
+            GREAT
+        )
+    );
+
+    dimensionedScalar rhoMin
+    (
+        dimensionedScalar::lookupOrDefault
+        (
+            "rhoMin",
+            pimple.dict(),
+            dimDensity,
+            0
+        )
+    );
+
+    dimensionedScalar pMax
+    (
+        dimensionedScalar::lookupOrDefault
+        (
+            "pMax",
+            pimple.dict(),
+            dimPressure,
+            GREAT
+        )
+    );
+
+    dimensionedScalar pMin
+    (
+        dimensionedScalar::lookupOrDefault
+        (
+            "pMin",
+            pimple.dict(),
+            dimPressure,
+            -GREAT
+        )
+    );
+
+    dimensionedScalar UMax
+    (
+        dimensionedScalar::lookupOrDefault
+        (
+            "UMax",
+            pimple.dict(),
+            dimVelocity,
+            -1.0
+        )
+    );
+
+    Info<< "Creating turbulence model\n" << endl;
+    autoPtr<compressible::turbulenceModel> turbulence
+    (
+        compressible::turbulenceModel::New
+        (
+            rho,
+            U,
+            phi,
+            thermo
+        )
+    );
+
+#if OPENFOAM_VERSION_MAJOR >= 6
+    Info<< "Creating combustion model\n" << endl;
+    autoPtr<CombustionModel<rhoReactionThermo>> combustion
+    (
+        CombustionModel<rhoReactionThermo>::New(thermo, turbulence())
+    );
+#endif
+
+    label pRefCell = 0;
+    scalar pRefValue = 0.0;
+    setRefCell(p, pimple.dict(), pRefCell, pRefValue);
+
+    mesh.setFluxRequired(p.name());
+
+    Info<< "Creating field dpdt\n" << endl;
+    volScalarField dpdt
+    (
+        IOobject
+        (
+            "dpdt",
+            runTime.timeName(),
+            mesh
+        ),
+        mesh,
+        dimensionedScalar("dpdt", p.dimensions()/dimTime, 0)
+    );
+
+    Info<< "Creating field kinetic energy K\n" << endl;
+    volScalarField K("K", 0.5*magSqr(U));
+
+#if OPENFOAM_VERSION_MAJOR < 5
+    volScalarField dQ
+    (
+        IOobject
+        (
+            "dQ",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
+    );
+#else
+    volScalarField Qdot
+    (
+        IOobject
+        (
+            "Qdot",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("Qdot", dimEnergy/dimVolume/dimTime, 0.0)
+    );
+#endif
+
+    Info<< "\nReading momentum exchange field Ksl\n" << endl;
+    volScalarField Ksl
+    (
+        IOobject
+        (
+            "Ksl",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("0", dimensionSet(1, -3, -1, 0, 0), 0.0)
+    );
+
+
+    Info<< "Reading particle velocity field Us\n" << endl;
+    volVectorField Us
+    (
+        IOobject
+        (
+            "Us",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    volScalarField molarConc
+    (
+        IOobject
+        (
+            "molarConc",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("zero",dimensionSet(0, -3, 0, 0, 1),0)
+    );
+
+    volVectorField UsRec
+    (
+        IOobject
+        (
+            "UsRec",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        Us
+    );
+
+
+    dimensionedScalar kf("0", dimensionSet(1, 1, -3, -1, 0, 0, 0), 0.026);
+
+//===============================

applications/solvers/rcfdemSolverRhoSteadyPimpleChem/limitP.H

@@ -0,0 +1,2 @@
+p = max(p, pMin);
+p = min(p, pMax);

applications/solvers/rcfdemSolverRhoSteadyPimpleChem/limitU.H

@@ -0,0 +1,11 @@
+if (UMax.value() > 0)
+{
+    forAll(U,cellI)
+    {
+        scalar mU(mag(U[cellI]));
+        if (mU > UMax.value())
+        {
+            U[cellI] *= UMax.value() / mU;
+        }
+    }
+}

applications/solvers/rcfdemSolverRhoSteadyPimpleChem/molConc.H

@@ -0,0 +1,12 @@
+{
+    molarConc = 0.0 * molarConc;
+    forAll(Y, i)
+    {
+        volScalarField& Yi = Y[i];
+        dimensionedScalar mi("mi",dimensionSet(1, 0, 0, 0, -1),composition.W(i));
+        mi /= 1000.0; // g to kg
+        molarConc += rho * Yi / mi;
+    }
+}
+
+// ************************************************************************* //

applications/solvers/rcfdemSolverRhoSteadyPimpleChem/monitorMass.H

@@ -0,0 +1,7 @@
+{
+    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;
+}

applications/solvers/rcfdemSolverRhoSteadyPimpleChem/pEqn.H

@@ -0,0 +1,96 @@
+rho = thermo.rho();
+rho = max(rho, rhoMin);
+rho = min(rho, rhoMax);
+rho.relax();
+
+if (stepcounter%nEveryFlow==0)
+{
+
+volScalarField rAU(1.0/UEqn.A());
+surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rhoeps*rAU));
+if (modelType=="A")
+{
+    rhorAUf *= fvc::interpolate(voidfractionRec);
+}
+volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
+
+surfaceScalarField phiUs("phiUs", fvc::interpolate(rhoeps*rAU*Ksl*UsRec)& mesh.Sf());
+
+
+if (pimple.transonic())
+{
+//  transonic version not implemented yet
+}
+else
+{
+    surfaceScalarField phiHbyA
+    (
+        "phiHbyA",
+        (
+            fvc::flux(rhoeps*HbyA)
+        )
+    );
+
+    // flux without pressure gradient contribution
+    phi = phiHbyA + phiUs;
+
+    // Update the pressure BCs to ensure flux consistency
+    constrainPressure(p, rhoeps, U, phi, rhorAUf);
+
+    volScalarField SmbyP(particleCloud.chemistryM(0).Sm() / p);
+
+    while (pimple.correctNonOrthogonal())
+    {
+        // Pressure corrector
+        fvScalarMatrix pEqn
+        (
+            fvc::div(phi)
+          - fvm::laplacian(rhorAUf, p)
+          ==
+            fvm::Sp(SmbyP, p)
+          + fvOptions(psi, p, rho.name())
+        );
+
+        pEqn.setReference(pRefCell, pRefValue);
+
+        pEqn.solve();
+
+        if (pimple.finalNonOrthogonalIter())
+        {
+            phi += pEqn.flux();
+        }
+    }
+}
+
+#include "rhoEqn.H"
+#include "compressibleContinuityErrsPU.H"
+
+// Explicitly relax pressure for momentum corrector
+p.relax();
+
+#include "limitP.H"
+
+// Recalculate density from the relaxed pressure
+rho = thermo.rho();
+rho = max(rho, rhoMin);
+rho = min(rho, rhoMax);
+rho.relax();
+Info<< "rho max/min : " << max(rho).value()
+    << " " << min(rho).value() << endl;
+
+if (modelType=="A")
+{
+    U = HbyA - rAU*(voidfractionRec*fvc::grad(p)-Ksl*UsRec);
+}
+else
+{
+    U = HbyA - rAU*(fvc::grad(p)-Ksl*UsRec);
+}
+
+#include "limitU.H"
+
+U.correctBoundaryConditions();
+fvOptions.correct(U);
+K = 0.5*magSqr(U);
+
+}

applications/solvers/rcfdemSolverRhoSteadyPimpleChem/rcfdemSolverRhoSteadyPimpleChem.C

@@ -0,0 +1,218 @@
+/*---------------------------------------------------------------------------*\
+License
+
+    This is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    This code is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this code.  If not, see <http://www.gnu.org/licenses/>.
+
+    Copyright (C) 2015- Thomas Lichtenegger, JKU Linz, Austria
+
+Application
+    rcfdemSolverRhoSteadyPimpleChem
+
+Description
+    Transient (DEM) + steady-state (CFD) solver for compressible flow using the
+    flexible PIMPLE (PISO-SIMPLE) algorithm. Particle-motion is obtained from
+    a recurrence process.
+    Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
+    The code is an evolution of the solver rhoPimpleFoam in OpenFOAM(R) 4.x,
+    where additional functionality for CFD-DEM coupling is added.
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+//#include "psiThermo.H"
+#include "turbulentFluidThermoModel.H"
+#if OPENFOAM_VERSION_MAJOR < 6
+#include "rhoCombustionModel.H"
+#else
+#include "rhoReactionThermo.H"
+#include "CombustionModel.H"
+#endif
+#include "bound.H"
+#include "pimpleControl.H"
+#include "fvOptions.H"
+#include "localEulerDdtScheme.H"
+#include "fvcSmooth.H"
+
+#include "cfdemCloudRec.H"
+#include "recBase.H"
+#include "recModel.H"
+#include "recPath.H"
+
+#include "cfdemCloudEnergy.H"
+#include "implicitCouple.H"
+#include "clockModel.H"
+#include "smoothingModel.H"
+#include "forceModel.H"
+#include "thermCondModel.H"
+#include "energyModel.H"
+#include "chemistryModel.H"
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "postProcess.H"
+
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createControl.H"
+    #include "createTimeControls.H"
+    #include "createRDeltaT.H"
+
+    #include "initContinuityErrs.H"
+    #include "createFields.H"
+    //#include "createFieldRefs.H"
+    #include "createFvOptions.H"
+
+    // create cfdemCloud
+    //#include "readGravitationalAcceleration.H"
+    cfdemCloudRec<cfdemCloudEnergy> particleCloud(mesh);
+    #include "checkModelType.H"
+    recBase recurrenceBase(mesh);
+    #include "updateFields.H"
+
+    turbulence->validate();
+    //#include "compressibleCourantNo.H"
+    //#include "setInitialDeltaT.H"
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    label recTimeIndex = 0;
+    scalar recTimeStep = recurrenceBase.recM().recTimeStep();
+    scalar startTime = runTime.startTime().value();
+
+    const IOdictionary& couplingProps = particleCloud.couplingProperties();
+    label nEveryFlow(couplingProps.lookupOrDefault<label>("nEveryFlow",1));
+    Info << "Solving flow equations for U and p every " << nEveryFlow << " steps.\n" << endl;
+    label stepcounter = 0;
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    scalar m(0.0);
+    scalar m0(0.0);
+    label counter(0);
+    p.storePrevIter();
+
+    while (runTime.run())
+    {
+        #include "readTimeControls.H"
+        #include "compressibleCourantNo.H"
+        #include "setDeltaT.H"
+
+        runTime++;
+
+        particleCloud.clockM().start(1,"Global");
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        // do particle stuff
+        particleCloud.clockM().start(2,"Coupling");
+        bool hasEvolved = particleCloud.evolve(voidfraction,Us,U);
+
+        //voidfraction = voidfractionRec;
+        //Us = UsRec;
+
+        if(hasEvolved)
+        {
+            particleCloud.smoothingM().smoothen(particleCloud.forceM(0).impParticleForces());
+        }
+
+        Info << "update Ksl.internalField()" << endl;
+        Ksl = particleCloud.momCoupleM(0).impMomSource();
+        Ksl.correctBoundaryConditions();
+
+        //Force Checks
+        vector fTotal(0,0,0);
+        vector fImpTotal = sum(mesh.V()*Ksl.primitiveFieldRef()*(Us.primitiveFieldRef()-U.primitiveFieldRef()));
+        reduce(fImpTotal, sumOp<vector>());
+        Info << "TotalForceExp: " << fTotal << endl;
+        Info << "TotalForceImp: " << fImpTotal << endl;
+
+        #include "solverDebugInfo.H"
+        particleCloud.clockM().stop("Coupling");
+
+        particleCloud.clockM().start(26,"Flow");
+        volScalarField rhoeps("rhoeps",rho*voidfractionRec);
+
+        while (pimple.loop())
+        {
+            // if needed, perform drag update here
+#if OPENFOAM_VERSION_MAJOR < 6
+            if (pimple.nCorrPIMPLE() <= 1)
+#else
+            if (pimple.nCorrPimple() <= 1)
+#endif
+            {
+                #include "rhoEqn.H"
+            }
+
+            // --- Pressure-velocity PIMPLE corrector loop
+
+
+            #include "UEqn.H"
+            #include "EEqn.H"
+
+            // --- Pressure corrector loop
+            while (pimple.correct())
+            {
+                // besides this pEqn, OF offers a "pimple consistent"-option
+                #include "molConc.H"
+                #include "pEqn.H"
+                rhoeps=rho*voidfractionRec;
+            }
+            #include "YEqn.H"
+
+            if (pimple.turbCorr())
+            {
+                turbulence->correct();
+            }
+        }
+
+        #include "monitorMass.H"
+
+        stepcounter++;
+        particleCloud.clockM().stop("Flow");
+
+        particleCloud.clockM().start(31,"postFlow");
+        particleCloud.postFlow();
+        particleCloud.clockM().stop("postFlow");
+
+        particleCloud.clockM().start(32,"ReadFields");
+        if ( runTime.timeOutputValue() - startTime - (recTimeIndex+1)*recTimeStep + 1.0e-5 > 0.0 )
+        {
+            recurrenceBase.updateRecFields();
+            #include "updateFields.H"
+            recTimeIndex++;
+        }
+        particleCloud.clockM().stop("ReadFields");
+
+        runTime.write();
+
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+
+
+        particleCloud.clockM().stop("Global");
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/rcfdemSolverRhoSteadyPimpleChem/rhoEqn.H

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

applications/solvers/rcfdemSolverRhoSteadyPimpleChem/updateFields.H

@@ -1,4 +1,8 @@
+// is it neccessary to extend recurrence path?
+if(recurrenceBase.recM().endOfPath())
+{
+    recurrenceBase.extendPath();
+}
+
 recurrenceBase.recM().exportVolScalarField("voidfraction",voidfractionRec);
-recurrenceBase.recM().exportVolVectorField("U",URec);
 recurrenceBase.recM().exportVolVectorField("Us",UsRec);
-recurrenceBase.recM().exportSurfaceScalarField("phi",phiRec);

applications/solvers/rctfSpeciesTransport/CEq.H

@@ -0,0 +1,19 @@
+
+volScalarField alphaEff("alphaEff", turbulence->nu()/Sc + alphat);
+
+CEqn =
+(
+    fvm::ddt(C)
+    + fvm::div(phiRec, C)
+    - fvm::laplacian(alphaEff, C)
+    ==
+    fvOptions(C)
+);
+
+CEqn.relax(relaxCoeff);
+
+fvOptions.constrain(CEqn);
+
+CEqn.solve();
+
+fvOptions.correct(C);

applications/solvers/rctfSpeciesTransport/Make/files

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

applications/solvers/rctfSpeciesTransport/Make/options

@@ -0,0 +1,27 @@
+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 \
+    -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 \
+    -I$(CFDEM_SRC_DIR)/recurrence/lnInclude \
+    -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/derived/cfdemCloudRec \
+
+EXE_LIBS = \
+    -L$(CFDEM_LIB_DIR)\
+    -lrecurrence \
+    -lturbulenceModels \
+    -lincompressibleTurbulenceModels \
+    -lincompressibleTransportModels \
+    -lfiniteVolume \
+    -lmeshTools \
+    -lfvOptions \
+    -l$(CFDEM_LIB_NAME) \
+     $(CFDEM_ADD_LIB_PATHS) \
+     $(CFDEM_ADD_LIBS)

applications/solvers/rctfSpeciesTransport/continuityErrCalc.H

@@ -0,0 +1,17 @@
+// calculate the continuity error according to phiRec
+
+{
+    volScalarField contErr(fvc::div(phiRec));
+
+    scalar sumLocalContErr = runTime.deltaTValue()*
+        mag(contErr)().weightedAverage(mesh.V()).value();
+
+    scalar globalContErr = runTime.deltaTValue()*
+        contErr.weightedAverage(mesh.V()).value();
+    cumulativeContErr += globalContErr;
+
+    Info<< "time step continuity errors : sum local = " << sumLocalContErr
+        << ", global = " << globalContErr
+        << ", cumulative = " << cumulativeContErr
+        << endl;
+}

applications/solvers/rctfSpeciesTransport/createFields.H

@@ -0,0 +1,281 @@
+//creating the fields according to the recurrence dictionary
+
+IOdictionary recProperties_
+(
+    IOobject
+    (
+        "recProperties0",
+        runTime.constant(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::NO_WRITE
+    )
+);
+
+List<wordList> fieldsDict_(recProperties_.lookup("fieldsPairs"));
+
+wordList fieldNames(fieldsDict_.size());
+
+for(int i = 0; i < fieldsDict_.size(); i++)
+{
+
+    fieldNames[i]= fieldsDict_[i][0];
+}
+
+Info<< "\n list of the fields: \n" << fieldNames << endl;
+
+//reading coherent velocity field name
+label k = findIndex(fieldNames,"coh_velocity");
+
+if (k < 0)
+{
+    FatalError <<"\n No field is defiened for the coherent velocity\n" << abort(FatalError);
+}
+const word Ucoh_pair = fieldsDict_[k][1];
+
+    volVectorField UcohRec
+    (
+        IOobject
+        (
+            Ucoh_pair,
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedVector("zero",dimensionSet(0, 1, -1, 0, 0),vector::zero)
+    );
+
+//reading incoherent velocity field name
+k = findIndex(fieldNames,"inc_velocity");
+
+if (k < 0)
+{
+    FatalError <<"\n No field is defiened for the incoherent velocity\n" << abort(FatalError);
+}
+const word Uinc_pair = fieldsDict_[k][1];
+
+    volVectorField UincRec
+    (
+        IOobject
+        (
+            Uinc_pair,
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedVector("zero",dimensionSet(0, 1, -1, 0, 0),vector::zero)
+    );
+
+//reading coherent turb kinetic energy field name
+k = findIndex(fieldNames,"kSGS_coh");
+if (k < 0)
+{
+    FatalError <<"\n No field is defiened for the coherent subgrid-scale turbulent kinetic energy\n" << abort(FatalError);
+}
+const word kSGScoh_pair = fieldsDict_[k][1];
+
+    volScalarField kcohRec
+    (
+        IOobject
+        (
+            kSGScoh_pair,
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("zero",dimensionSet(0, 2, -2, 0, 0),0.0)
+     );
+
+//reading incoherent turb kinetic energy field name
+k = findIndex(fieldNames,"kSGS_inc");
+if (k < 0)
+{
+    FatalError <<"\n No field is defiened for the coherent subgrid-scale turbulent kinetic energy\n" << abort(FatalError);
+}
+const word kSGSinc_pair = fieldsDict_[k][1];
+
+    volScalarField kincRec
+    (
+        IOobject
+        (
+            kSGSinc_pair,
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("zero",dimensionSet(0, 2, -2, 0, 0),0.0)
+     );
+
+// calculated fields
+Info<< "\nCreating cell volume field\n" << endl;
+
+volScalarField delta
+(
+    IOobject
+    (
+        "delta",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+   mesh,
+   dimensionedScalar("delta", dimLength, 0.0)
+);
+
+delta.primitiveFieldRef()=pow(mesh.V(),1.0/3.0);
+delta.write();
+
+
+volVectorField URec
+(
+    IOobject
+    (
+        "URec",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);
+Info<< "\nCreating turb kinetic energy field\n" << endl;
+
+volScalarField kRec
+(
+    IOobject
+    (
+        "kRec",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    kcohRec+kincRec
+);
+
+// check if there is any negative values
+forAll(kRec, cellI)
+{
+    if (kRec[cellI] < SMALL)
+    {
+        kRec[cellI] = 0.0;
+    }
+
+}
+
+const fvPatchList& patches = mesh.boundary();
+forAll(patches, patchI)
+{
+        kRec.boundaryFieldRef()[patchI] = 0.0;
+}
+
+
+kRec.write();
+
+Info<< "\nCreating turb viscosity field\n" << endl;
+volScalarField nutRec
+(
+    IOobject
+    (
+        "nutRec",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    sqrt(kRec)*delta*0.094
+);
+
+nutRec.write();
+
+
+Info<< "Calculating face flux field phiRec\n" << endl;
+surfaceScalarField phiRec
+ (
+     IOobject
+     (
+        "phiRec",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+     ),
+     linearInterpolate(URec) & mesh.Sf()
+ );
+
+    phiRec.write();
+
+    singlePhaseTransportModel laminarTransport(URec, phiRec);
+
+    autoPtr<incompressible::turbulenceModel> turbulence
+    (
+        incompressible::turbulenceModel::New(URec, phiRec, laminarTransport)
+    );
+  
+    dimensionedScalar Sc("Sc", dimless, laminarTransport);
+    dimensionedScalar Sct("Sct", dimless, laminarTransport);
+  
+volScalarField alphat
+(
+    IOobject
+    (
+        "alphat",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::AUTO_WRITE
+    ),
+    nutRec/Sct
+);
+
+ // create the scalar field
+ Info<< "Creating scalar transport field\n" << endl;
+
+    volScalarField C
+    (
+        IOobject
+        (
+            "C",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+    
+    fvScalarMatrix CEqn(C, dimless*dimVolume/(dimTime));
+
+    scalar relaxCoeff(0.0);
+
+ Info<< "reading clockProperties\n" << endl;
+
+    IOdictionary clockProperties
+    (
+            IOobject
+            (
+                "clockProperties",
+                mesh.time().constant(),
+                mesh,
+                IOobject::MUST_READ,
+                IOobject::NO_WRITE
+            )
+    );
+
+    autoPtr<clockModel> myClock
+    (
+            clockModel::New
+            (
+                clockProperties,
+                mesh.time()
+            )
+    );

applications/solvers/rctfSpeciesTransport/createRecBase.H

@@ -0,0 +1,66 @@
+// check which recProperties dicts are present, read them in and construct a PtrList of recBases
+// names for dicts can be "recProperties" or "recPropertiesN" where N in {0, 1, ...} 
+
+#include "error.H"
+
+word dictName = "recProperties";
+wordList recPropertiesList(0);
+PtrList <recBase> recBases(0);
+label maxDictNumber = 100;
+
+{
+    IOdictionary recPropDict
+    (
+        IOobject
+        (
+            dictName,
+            mesh.time().constant(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::NO_WRITE
+        )
+    );
+
+    if (recPropDict.headerOk())
+    {
+        recPropertiesList.append(dictName);
+    }
+}
+
+
+for (label counter = 0; counter < maxDictNumber; counter++)
+{
+    word dictNameIter = dictName + Foam::name(counter);
+    IOdictionary recPropDict
+    (
+        IOobject
+        (
+            dictNameIter,
+            mesh.time().constant(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::NO_WRITE
+        )
+    );
+
+    if (recPropDict.headerOk())
+    {
+        recPropertiesList.append(dictNameIter);
+    }
+}
+
+if (recPropertiesList.size() == 0)
+{
+    FatalError << "no recProperties dicts found" << endl;
+}
+
+else
+{
+    Info << "found " << recPropertiesList.size() << " dicts with names " << recPropertiesList << endl;
+}
+
+
+for (label counter = 0; counter < recPropertiesList.size(); counter++)
+{
+    recBases.append( new recBase(mesh, recPropertiesList[counter]));
+}

applications/solvers/rctfSpeciesTransport/rctfSpeciesTransport.C

@@ -0,0 +1,131 @@
+/*---------------------------------------------------------------------------*\
+    CFDEMcoupling academic - Open Source CFD-DEM coupling
+    
+    Contributing authors:
+    Thomas Lichtenegger, Gerhard Holzinger, Sanaz Abbasi
+    Copyright (C) 2015- Johannes Kepler University, Linz
+-------------------------------------------------------------------------------
+License
+    This file is part of CFDEMcoupling academic.
+
+    CFDEMcoupling academic 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 academic 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 academic.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    Turbulent Transport Recurrence Solver for modal decomposition
+
+Description
+    Solves a transport equation for a passive scalar on a single-phase solution
+    for a solver based on recurrence statistics
+
+Rules
+	Solution data to compute the recurrence statistics from, needs to 
+        reside in $CASE_ROOT/dataBase(0...N)
+    Time step data in the first dataBase needs to be evenly spaced in time
+    A list of indices for the corresponding incoherent fields to coherent ones
+        should be provided.
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "singlePhaseTransportModel.H"
+#include "turbulentTransportModel.H"
+#include "fvOptions.H"
+
+#include "recBase.H"
+#include "recModel.H"
+
+#include "clockModel.H"
+
+#include "objectRegistry.H"
+#include "VectorSpace.H"
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "postProcess.H"
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createControl.H"
+    #include "createFields.H"
+    #include "createFvOptions.H"
+    scalar cumulativeContErr = 0;
+  
+    //create recBases according to a list of recProperties
+    #include "createRecBase.H"
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+
+    label recTimeIndex(0);
+    label currTimeIndex(0);
+
+    scalar recTimeStep_=recBases[0].recM().recTimeStep();
+    labelPairList incPairTimeIndex_(0);
+
+    IFstream pairFile("incIndexPairList");
+    pairFile >> incPairTimeIndex_;
+
+    while (runTime.run())
+    {
+
+        myClock().start(1,"Global");
+        runTime++;
+
+        myClock().start(11,"Total");
+        
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        myClock().start(2,"fieldUpdate");
+        
+        if ( runTime.timeOutputValue() - (recTimeIndex+1)*recTimeStep_ + 1.0e-5 > 0.0 )
+        {
+	    Info << "Updating fields at run time " << runTime.timeOutputValue()
+	        << " corresponding to recurrence time " << (recTimeIndex+1)*recTimeStep_ << ".\n" << endl;
+            recBases[0].updateRecFields();
+            #include "readFields.H"
+
+            recTimeIndex++;
+        }
+
+        myClock().stop("fieldUpdate");
+
+        #include "continuityErrCalc.H"
+
+        myClock().start(3,"speciesEqn");
+        #include "CEq.H"
+        myClock().stop("speciesEqn");
+
+        myClock().stop("Total");
+
+        runTime.write();
+        
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+        << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+        << nl << endl;
+
+        myClock().stop("Global");
+
+    }
+
+    
+    myClock().evalPar();
+    myClock().normHist();
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/rctfSpeciesTransport/readFields.H

@@ -0,0 +1,40 @@
+
+currTimeIndex = recBases[0].recM().currentTimeIndex();
+
+Info << "current Time Index = " << currTimeIndex << endl;
+
+recBases[0].recM().exportVolVectorField(Ucoh_pair,UcohRec);
+recBases[0].recM().exportVolScalarField(kSGScoh_pair,kcohRec);
+
+
+label incTimeIndex = incPairTimeIndex_[currTimeIndex][1];
+
+Info << "  incoherent pair Time Index = " << incTimeIndex  << endl;
+
+UincRec = recBases[1].recM().exportVolVectorField(Uinc_pair,incTimeIndex);
+kincRec = recBases[1].recM().exportVolScalarField(kSGSinc_pair,incTimeIndex);
+
+kRec = kcohRec+kincRec;
+
+forAll(kRec, cellI)
+{
+    if (kRec[cellI] < SMALL)
+    {
+        kRec[cellI] = 0.0;
+    }
+}
+
+const fvPatchList& patches = mesh.boundary();
+forAll(patches, patchI)
+{
+        kRec.boundaryFieldRef()[patchI] = 0.0;
+}
+
+URec = UcohRec + UincRec;
+phiRec = linearInterpolate(URec) & mesh.Sf();
+
+nutRec = sqrt(kRec)*delta*0.094;
+
+alphat = nutRec/Sct;
+alphat.correctBoundaryConditions();
+

applications/utilities/displacementField/Make/files

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

applications/utilities/displacementField/Make/options

@@ -0,0 +1,7 @@
+EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude
+
+EXE_LIBS = \
+    -lfiniteVolume \
+    -lmeshTools

applications/utilities/displacementField/displacementField.C

@@ -0,0 +1,458 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
+   \\    /   O peration     | Website:  https://openfoam.org
+    \\  /    A nd           | Copyright (C) 2011-2018 OpenFOAM Foundation
+     \\/     M anipulation  |
+-------------------------------------------------------------------------------
+License
+    This file is part of OpenFOAM.
+
+    OpenFOAM 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.
+
+    OpenFOAM 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 OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.
+
+Application
+    displacementField
+
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "vectorList.H"
+#include <sstream>
+#include <string>
+#include <sys/stat.h>
+#include <unistd.h>
+#include <set>
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+void findPairs(labelList &, labelList &, labelPairList &);
+void findPairsUnordered(labelList &, labelList &, labelPairList &);
+void fillEmptyCells(fvMesh &, label , labelList &, volVectorField &, volVectorField &, scalarList &, vector, vector, bool, scalar);
+void nearestNeighborCells(fvMesh &, label, label, labelList &, labelList &);
+void normalizeFields(labelList &, volVectorField &, volVectorField &);
+void readDump(std::string, labelList &, vectorList &);
+scalar weightFun(scalar);
+
+int main(int argc, char *argv[])
+{
+    argList::addOption
+    (
+        "totalProcs",
+        "label",
+        "total number of parallel processes, defaults to 1"
+    );
+    argList::addOption
+    (
+        "thisProc",
+        "label",
+        "number of current process, defaults to 0"
+    );
+
+
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+
+    const label thisProc = args.optionLookupOrDefault("thisProc", 0);
+    const label totalProcs = args.optionLookupOrDefault("totalProcs", 1);
+
+    Info << "This is number " << thisProc << " of " << totalProcs << " processes." << endl;
+
+
+    // user-defined input for each case
+    IOdictionary displacementProperties
+    (
+        IOobject
+        (
+            "displacementProperties",
+            mesh.time().constant(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        )
+    );
+
+    label dumpIndexStart(readLabel(displacementProperties.lookup("dumpIndexStart")));
+    label dumpIndexEnd(readLabel(displacementProperties.lookup("dumpIndexEnd")));
+    label dumpIndexInputIncrement(readLabel(displacementProperties.lookup("dumpIndexInputIncrement")));
+    label dumpIndexDisplacementIncrement(readLabel(displacementProperties.lookup("dumpIndexDisplacementIncrement")));
+    label nNeighMin(readLabel(displacementProperties.lookup("nNeighMin")));
+    scalar timePerInputStep(readScalar(displacementProperties.lookup("timePerInputStep")));
+    scalar timePerDisplacementStep(readScalar(displacementProperties.lookup("timePerDisplacementStep")));
+    scalar startTime(readScalar(displacementProperties.lookup("startTime")));
+    std::string filepath=string(displacementProperties.lookup("filepath"));
+    std::string fileext=string(displacementProperties.lookupOrDefault<string>("fileextension",""));
+    bool interpolate=bool(displacementProperties.lookupOrDefault<bool>("fillEmptyCells",true));
+    bool averageMode=bool(displacementProperties.lookupOrDefault<bool>("averageMode",false));
+    vector defaultUs=vector(displacementProperties.lookupOrDefault<vector>("defaultUs",vector::zero));
+    vector defaultUsDirectedVariance=vector(displacementProperties.lookupOrDefault<vector>("defaultUsDirectedVariance",vector::zero));
+
+    scalar xmin=scalar(displacementProperties.lookupOrDefault<scalar>("xmin",-1e10));
+    scalar xmax=scalar(displacementProperties.lookupOrDefault<scalar>("xmax",1e10));
+    scalar ymin=scalar(displacementProperties.lookupOrDefault<scalar>("ymin",-1e10));
+    scalar ymax=scalar(displacementProperties.lookupOrDefault<scalar>("ymax",1e10));
+    scalar zmin=scalar(displacementProperties.lookupOrDefault<scalar>("zmin",-1e10));
+    scalar zmax=scalar(displacementProperties.lookupOrDefault<scalar>("zmax",1e10));
+    scalarList boundaries(6);
+    boundaries[0]=xmin;
+    boundaries[1]=xmax;
+    boundaries[2]=ymin;
+    boundaries[3]=ymax;
+    boundaries[4]=zmin;
+    boundaries[5]=zmax;
+
+    vectorList probePoints=vectorList(displacementProperties.lookupOrDefault<vectorList>("probePoints",vectorList(0)));
+    bool monitorProbes = false;
+    if (probePoints.size()>0) monitorProbes = true;
+    #include "OFstream.H"
+    OFstream monitoringDataFile("monitoringData.txt");
+    if (monitorProbes)
+    {
+        monitoringDataFile << "# monitoring data file" << endl;
+        monitoringDataFile << "# format: time nPerCell[p1] UDisp[p1] UDispDirectedVariance[p1] nPerCell[p2] ... " << endl;
+        for(label p=0;p<probePoints.size();p++)
+        {
+            vector pos = probePoints[p];
+            monitoringDataFile << "# point[" << p << "] = " << pos << endl;
+        }
+    }
+
+
+    label dumpIndex1 = dumpIndexStart + thisProc * dumpIndexInputIncrement;
+    label dumpIndex2 = dumpIndex1 + dumpIndexDisplacementIncrement;
+
+    volVectorField Us
+    (
+        IOobject
+        (
+            "UDisp",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedVector("zero", dimensionSet(0,1,-1,0,0), vector::zero)
+    );
+
+    volVectorField UsDirectedVariance
+    (
+        IOobject
+        (
+            "UDispDirectedVariance",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedVector("zero", dimensionSet(0,1,-1,0,0), vector::zero)
+    );
+
+    labelList particlesInCell(mesh.nCells(), 0);
+
+    scalar currTime=startTime + thisProc * timePerInputStep;
+    label timeIndex=thisProc;
+
+    while(true)
+    {
+        runTime.setTime(currTime,timeIndex);
+        // read dump files and check which particle indices are present in both
+        labelList indices1, indices2;
+        vectorList positions1, positions2;
+
+        std::stringstream ss;
+        ss << filepath << dumpIndex1 << fileext;
+        std::string filename1 = ss.str();
+        ss.str("");
+        ss << filepath << dumpIndex2 << fileext;
+        std::string filename2 = ss.str();
+
+        if (access( filename1.c_str(), F_OK ) == -1 || access( filename2.c_str(), F_OK ) == -1 || dumpIndex2 > dumpIndexEnd)
+        {
+            if (averageMode)
+            {
+                normalizeFields(particlesInCell, Us, UsDirectedVariance);
+                fillEmptyCells(mesh,nNeighMin,particlesInCell,Us,UsDirectedVariance,boundaries,defaultUs,defaultUsDirectedVariance,interpolate,timePerDisplacementStep);
+
+                Us /= timePerDisplacementStep;
+                UsDirectedVariance /= timePerDisplacementStep;
+                Us.write();
+                UsDirectedVariance.write();
+            }
+            break;
+        }
+
+        Info << "\nReading" << endl;
+        Info << "\t" << filename1 << endl;
+        Info << "\t" << filename2 << endl;
+        Info << "corresponding to time = " << currTime << "." << endl;
+
+        readDump(filename1, indices1, positions1);
+        readDump(filename2, indices2, positions2);
+
+        labelPairList pairs;
+        findPairs(indices1,indices2,pairs);
+
+        // average particle displacements and their variance
+        Info << "Binning particle displacements on mesh." << endl;
+        vector position, displacement;
+        label line1, line2;
+        label cellI;
+
+        if (!averageMode)
+        {
+            Us *= 0.0;
+            UsDirectedVariance *= 0.0;
+            particlesInCell.clear();
+            particlesInCell.setSize(mesh.nCells(), 0);
+        }
+
+        for (label partI = 0; partI < pairs.size(); partI++)
+        {
+            line1 = pairs[partI].first();
+            line2 = pairs[partI].second();
+            position = positions1[line1];
+            displacement = positions2[line2] - positions1[line1];
+            cellI = mesh.findCell(position);
+            if (cellI < 0)  continue;
+            particlesInCell[cellI] += 1;
+            Us[cellI] += displacement;
+
+            for (label comp=0;comp<3;comp++)
+            {
+                UsDirectedVariance[cellI].component(comp) += displacement.component(comp)*displacement.component(comp);
+            }
+        }
+
+        if (!averageMode)
+        {
+            normalizeFields(particlesInCell, Us, UsDirectedVariance);
+            fillEmptyCells(mesh,nNeighMin,particlesInCell,Us,UsDirectedVariance,boundaries,defaultUs,defaultUsDirectedVariance,interpolate,timePerDisplacementStep);
+
+            Us /= timePerDisplacementStep;
+            UsDirectedVariance /= timePerDisplacementStep;
+            Us.write();
+            UsDirectedVariance.write();
+        }
+
+        if (averageMode && monitorProbes)
+        {
+            monitoringDataFile << currTime << " ";
+            for(label p=0;p<probePoints.size();p++)
+            {
+                vector pos = probePoints[p];
+                label cellP = mesh.findCell(pos);
+                monitoringDataFile << " " << particlesInCell[cellP] << " " << Us[cellP]/timePerDisplacementStep << " " << UsDirectedVariance[cellP]/(timePerDisplacementStep*timePerDisplacementStep);
+            }
+            monitoringDataFile << endl;
+        }
+
+        dumpIndex1 += dumpIndexInputIncrement*totalProcs;
+        dumpIndex2 += dumpIndexInputIncrement*totalProcs;
+        currTime += timePerInputStep*totalProcs;
+        timeIndex += totalProcs;
+    }
+    return 0;
+}
+
+void readDump(std::string filename, labelList &indices, vectorList &positions)
+{
+    #include <fstream>
+
+    const label leadingLines = 9;
+    label lineCounter = 0;
+    label partIndex;
+    scalar x, y, z;
+
+    indices.clear();
+    positions.clear();
+
+    std::ifstream file(filename);
+    std::string str;
+    while (std::getline(file, str))
+    {
+        if (lineCounter >= leadingLines)
+        {
+            sscanf(str.c_str(), "%d %lf %lf %lf", &partIndex, &x, &y, &z);
+            indices.append(partIndex);
+            positions.append(vector(x,y,z));
+        }
+        lineCounter++;
+    }
+}
+
+void findPairs(labelList &indices1, labelList &indices2, labelPairList &pairs)
+{
+    // remove all entries from first list if they are not present in second list
+    // this assumes ordered entries
+
+    if (indices2.size() == 0) return;
+
+    for (label i=0;i<indices1.size();i++)
+    {
+        label j1 = -1;
+        label j2 = indices2.size();
+        label jmid = 0;
+        label index1 = indices1[i];
+        while(true)
+        {
+            jmid = (j1+j2)/2;
+            if (indices2[jmid] > index1) j2 = jmid;
+            else if (indices2[jmid] < index1) j1 = jmid;
+            else
+            {
+                pairs.append(labelPair(i,jmid));
+                break;
+            }
+            if (j2-j1 == 1) break;
+        }
+    }
+    Info << "findPairs: " << pairs.size() << " pairs found." << endl;
+}
+
+void findPairsUnordered(labelList &indices1, labelList &indices2, labelPairList &pairs)
+{
+    // remove all entries from first list if they are not present in second list
+
+    for (label i=0;i<indices1.size();i++)
+    {
+        for (label j=0;j<indices2.size();j++)
+        {
+            if (indices1[i] == indices2[j])
+            {
+                pairs.append(labelPair(i,j));
+                break;
+            }
+        }
+    }
+    Info << "findPairs: " << pairs.size() << " pairs found." << endl;
+}
+
+void fillEmptyCells(fvMesh &mesh, label nNeighMin, labelList &particlesInCell, volVectorField &Us, volVectorField& UsDirectedVariance,scalarList& boundaries, vector defaultUs, vector defaultUsDirectedVariance, bool interpolate, scalar dt)
+{
+    labelList neighborsWithValues;
+    scalar neighborSqrDistance;
+    scalar weight;
+    scalar weightSum;
+    scalarList weights;
+
+    Info << "Filling empty cells." << endl;
+    forAll(mesh.C(), cellI)
+    {
+        if (particlesInCell[cellI] > 0) continue;
+
+        vector position = mesh.C()[cellI];
+        label outsideBox = 0;
+        if (position.x() < boundaries[0] || position.x() > boundaries[1]) outsideBox++;
+        if (position.y() < boundaries[2] || position.y() > boundaries[3]) outsideBox++;
+        if (position.z() < boundaries[4] || position.z() > boundaries[5]) outsideBox++;
+
+        if (outsideBox > 0 || !interpolate)
+        {
+            Us[cellI] = defaultUs*dt;
+            UsDirectedVariance[cellI] = defaultUsDirectedVariance*dt;
+            continue;
+        }
+
+        nearestNeighborCells(mesh, cellI, nNeighMin, particlesInCell, neighborsWithValues);
+        weightSum = 0.0;
+        weights.clear();
+        for (label neighI=0; neighI<neighborsWithValues.size(); neighI++)
+        {
+            neighborSqrDistance = magSqr(mesh.C()[cellI] - mesh.C()[neighborsWithValues[neighI]]);
+            weight = weightFun(neighborSqrDistance);
+            weights.append(weight);
+            weightSum += weight;
+        }
+        for (label neighI=0; neighI<neighborsWithValues.size(); neighI++)
+        {
+            weight = weights[neighI]/weightSum;
+            Us[cellI] += weight*Us[neighborsWithValues[neighI]];
+            UsDirectedVariance[cellI] += weight*UsDirectedVariance[neighborsWithValues[neighI]];
+        }
+
+        // make sure no particles are placed outside of domain
+        // TODO: correct following implementation (meshSearch) and test it
+/*
+        vector shiftedPosition = position + dt * Us[cellI];
+        label cellJ = mesh.findCell(shiftedPosition);
+        if (cellJ < 0)
+        {
+            label cellK = mesh.findNearestCellWalk(shiftedPosition,cellI);
+            Us[cellI] = (mesh.C()[cellI] - mesh.C()[cellK]) / dt;
+        }
+*/
+    }
+}
+
+void nearestNeighborCells(fvMesh &mesh, label refCell, label nNeighMin, labelList &particlesInCell, labelList &neighborsWithValues)
+{
+    std::set<label> neighbors;
+    std::set<label> newNeighbors;
+    std::set<label> recentNeighbors;
+
+    neighbors.insert(refCell);
+    recentNeighbors.insert(refCell);
+
+    neighborsWithValues.clear();
+
+    while(neighborsWithValues.size() < nNeighMin)
+    {
+        for (std::set<label>::iterator it=recentNeighbors.begin(); it!=recentNeighbors.end(); ++it)
+        {
+            labelList adjacent = mesh.cellCells()[*it];
+            label adj;
+            for (label j=0; j<adjacent.size(); j++)
+            {
+                adj = adjacent[j];
+                std::set<label>::iterator it2 = neighbors.find(adj);
+                if (it2 == neighbors.end())
+                {
+                    newNeighbors.insert(adj);
+                    neighbors.insert(adj);
+                    if (particlesInCell[adj] > 0) neighborsWithValues.append(adj);
+                }
+            }
+        }
+
+        if (newNeighbors.size() == 0) return;
+        recentNeighbors.clear();
+        recentNeighbors = newNeighbors;
+        newNeighbors.clear();
+    }
+}
+
+void normalizeFields(labelList& particlesInCell, volVectorField& Us, volVectorField & UsDirectedVariance)
+{
+    for (label cellJ = 0; cellJ<particlesInCell.size(); cellJ++)
+    {
+        if (particlesInCell[cellJ] > 0)
+        {
+            Us[cellJ] /= particlesInCell[cellJ];
+            UsDirectedVariance[cellJ] /= particlesInCell[cellJ];
+            for (label comp=0;comp<3;comp++)
+            {
+                UsDirectedVariance[cellJ].component(comp) -= Us[cellJ].component(comp)*Us[cellJ].component(comp);
+                if (UsDirectedVariance[cellJ].component(comp) > 0) UsDirectedVariance[cellJ].component(comp) = Foam::sqrt(UsDirectedVariance[cellJ].component(comp));
+            }
+        }
+    }
+}
+
+scalar weightFun(scalar distSqr)
+{
+    // inverse distance weighting, order 2
+    return 1.0/distSqr;
+}
+// ************************************************************************* //

applications/utilities/rBaseMirror/rBaseMirror/Make/files

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

applications/utilities/rBaseMirror/rBaseMirror/createFields.H

@@ -13,5 +13,7 @@
     vector refPoint(mirrorProperties.lookup("refPoint"));
     vector refDirection(mirrorProperties.lookup("refDirection"));
 
-    word fieldName(mirrorProperties.lookup("fieldName"));
     word dataBaseName(mirrorProperties.lookup("dataBaseName"));
+
+    const wordList volScalarFieldNames(mirrorProperties.lookup("volScalarFields"));
+    const wordList volVectorFieldNames(mirrorProperties.lookup("volVectorFields"));

applications/utilities/rBaseMirror/rBaseMirror/rBaseMirror.C

@@ -66,16 +66,8 @@ int main(int argc, char *argv[])
     instantList timeDirs(recTime.times());
     recTime.setTime(timeDirs[0],0);
 
-    #include "readFields.H"
-
-    Info << fieldName << endl;
-
-    volScalarField transformedField = origField;
-
     scalar t;
 
-    label shiftedTimeI = 0;
-
     // check number of time directories
     label shift = 0;
     forAll(timeDirs, timeI)
@@ -95,14 +87,30 @@ int main(int argc, char *argv[])
     label cellI_transformed = -1;
     forAll(timeDirs, timeI)
     {
-
        recTime.setTime(timeDirs[timeI], timeI);
        t = recTime.value();
        if(t < startTime) continue;
        if(t > endTime) continue;
        Info << "time = " << t << ", time index = " << timeI << endl;
 
-       #include "readFields.H"
+       // volScalarFields
+       for (int sf = 0; sf < volScalarFieldNames.size(); sf++)
+       {
+           word fieldName = volScalarFieldNames[sf];
+           volScalarField origField
+           (
+               IOobject
+               (
+                    fieldName,
+                    recTime.timePath(),
+                    mesh,
+                    IOobject::MUST_READ,
+                    IOobject::NO_WRITE
+               ),
+                mesh
+           );
+
+           volScalarField transformedField = origField;
 
            forAll(transformedField, cellI)
            {
@@ -121,13 +129,54 @@ int main(int argc, char *argv[])
                transformedField[cellI] = transformedValue;
            }
 
-       shiftedTimeI = timeI + shift;
-       t = recTime.value() + origTimeRange + dt;
-       runTime.setTime(t, shiftedTimeI);
-       Info << "creating transformed fields for time = " << t << ", time index = " << shiftedTimeI << endl;
+           runTime.setTime(recTime.value() + origTimeRange + dt, timeI + shift);
+           Info << "creating transformed field " << fieldName << " for time = " << recTime.value() + origTimeRange + dt << endl;
            transformedField.write();
        }
 
+
+       // volVectorFields
+       for (int vf = 0; vf < volVectorFieldNames.size(); vf++)
+       {
+           word fieldName = volVectorFieldNames[vf];
+           volVectorField origField
+           (
+               IOobject
+               (
+                    fieldName,
+                    recTime.timePath(),
+                    mesh,
+                    IOobject::MUST_READ,
+                    IOobject::NO_WRITE
+               ),
+                mesh
+           );
+
+           volVectorField transformedField = origField;
+
+           forAll(transformedField, cellI)
+           {
+               vector position = mesh.C()[cellI];
+               vector transformedPosition = 2 * ((refPoint - position) & refDirection) * refDirection / (refDirection & refDirection) + position;
+               cellI_transformed = mesh.findCell(transformedPosition);
+               if(cellI_transformed < 0)
+               {
+                   Info << "Couldn't find transformed cell. Stopping." << endl;
+                   return 0;
+               }
+
+               vector value = origField[cellI_transformed];
+               vector transformedValue = -2 * (value & refDirection) * refDirection / (refDirection & refDirection) + value;
+
+               transformedField[cellI] = transformedValue;
+           }
+
+           runTime.setTime(recTime.value() + origTimeRange + dt, timeI + shift);
+           Info << "creating transformed field " << fieldName << " for time = " << recTime.value() + origTimeRange + dt << endl;
+           transformedField.write();
+        }
+    }
+
     Info << "\nEnd" << endl;
 
     return 0;

applications/utilities/rBaseMirror/rBaseMirrorScalar/Make/files

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

applications/utilities/rBaseMirror/rBaseMirrorScalar/createFields.H

@@ -1,17 +0,0 @@
-    IOdictionary mirrorProperties
-    (
-        IOobject
-        (
-            "mirrorProperties",
-            mesh.time().constant(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE
-        )
-    );
-
-    vector refPoint(mirrorProperties.lookup("refPoint"));
-    vector refDirection(mirrorProperties.lookup("refDirection"));
-
-    word fieldName(mirrorProperties.lookup("fieldName"));
-    word dataBaseName(mirrorProperties.lookup("dataBaseName"));

applications/utilities/rBaseMirror/rBaseMirrorScalar/readFields.H

@@ -1,13 +0,0 @@
-
-    volScalarField origField
-    (
-        IOobject
-        (
-            fieldName,
-            recTime.timePath(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh
-    );

applications/utilities/rBaseMirror/rBaseMirrorVec/Make/files

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

applications/utilities/rBaseMirror/rBaseMirrorVec/Make/options

@@ -1,17 +0,0 @@
-include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
-
-
-EXE_INC = \
-     $(PFLAGS) \
-    -I$(LIB_SRC)/finiteVolume/cfdTools \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
-    -I$(LIB_SRC)/sampling/lnInclude \
-    -I$(LIB_SRC)/fvOptions/lnInclude \
-    -Wno-deprecated-copy
-
-EXE_LIBS = \
-    -lfiniteVolume \
-    -lmeshTools \
-    -lsampling \
-    -lfvOptions

applications/utilities/rBaseMirror/rBaseMirrorVec/rBaseMirrorVec.C

@@ -1,134 +0,0 @@
-/*---------------------------------------------------------------------------*\
-  =========                 |
-  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
-   \\    /   O peration     |
-    \\  /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
-     \\/     M anipulation  |
--------------------------------------------------------------------------------
-License
-    This file is part of OpenFOAM.
-
-    OpenFOAM 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 2 of the License, or (at your
-    option) any later version.
-
-    OpenFOAM 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 OpenFOAM; if not, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-Application
-    rBaseMirror
-
-Description
-    Read time series and extend it by mirrored fields if geometry possesses
-    the same symmetry
-
-\*---------------------------------------------------------------------------*/
-
-#include "fvCFD.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-//  Main program:
-
-int main(int argc, char *argv[])
-{
-    argList::noParallel();
-    timeSelector::addOptions();
-
-    #include "setRootCase.H"
-    #include "createTime.H"
-
-   // read in start and end time from controlDict
-
-    scalar startTime=runTime.startTime().value();
-    scalar endTime=runTime.endTime().value();
-    scalar origTimeRange = endTime - startTime;
-
-    Info << "start time = " << runTime.startTime() << endl;
-    Info << "end time = " << runTime.endTime() << endl;
-
-    // check which time directories are present
-    //instantList timeDirs = timeSelector::select0(runTime, args);
-    //runTime.setTime(timeDirs[0], 0);
-
-    #include "createMesh.H"
-
-    #include "createFields.H"
-
-    Foam::Time recTime(fileName(dataBaseName), "", "../system", "../constant", false);
-    instantList timeDirs(recTime.times());
-    recTime.setTime(timeDirs[0],0);
-
-    #include "readFields.H"
-
-    Info << fieldName << endl;
-
-    volVectorField transformedField = origField;
-
-    scalar t;
-
-    label shiftedTimeI = 0;
-
-    // check number of time directories
-    label shift = 0;
-    forAll(timeDirs, timeI)
-    {
-
-       if (recTime.timeName() == "constant") continue;
-       recTime.setTime(timeDirs[timeI], timeI);
-       t = recTime.value();
-       if(t < startTime) continue;
-       if(t > endTime) continue;
-       shift++;
-    }
-
-    scalar dt = origTimeRange / (shift - 1.0);
-    recTime.setEndTime(startTime + 2 * origTimeRange + dt);
-
-    label cellI_transformed = -1;
-    forAll(timeDirs, timeI)
-    {
-       recTime.setTime(timeDirs[timeI], timeI);
-       t = recTime.value();
-       if(t < startTime) continue;
-       if(t > endTime) continue;
-       Info << "time = " << t << ", time index = " << timeI << endl;
-
-       #include "readFields.H"
-
-       forAll(transformedField, cellI)
-       {
-           vector position = mesh.C()[cellI];
-           vector transformedPosition = 2 * ((refPoint - position) & refDirection) * refDirection / (refDirection & refDirection) + position;
-           cellI_transformed = mesh.findCell(transformedPosition);
-           if(cellI_transformed < 0)
-           {
-               Info << "Couldn't find transformed cell. Stopping." << endl;
-               return 0;
-           }
-
-           vector value = origField[cellI_transformed];
-           vector transformedValue = -2 * (value & refDirection) * refDirection / (refDirection & refDirection) + value;
-
-           transformedField[cellI] = transformedValue;
-       }
-
-       shiftedTimeI = timeI + shift;
-       t = recTime.value() + origTimeRange + dt;
-       runTime.setTime(t, shiftedTimeI);
-       Info << "creating transformed fields for time = " << t << ", time index = " << shiftedTimeI << endl;
-       transformedField.write();
-    }
-
-    Info << "\nEnd" << endl;
-
-    return 0;
-}
-
-// ************************************************************************* //

applications/utilities/rBaseMirror/rBaseMirrorVec/readFields.H

@@ -1,14 +0,0 @@
-
-volVectorField origField
-    (
-        IOobject
-        (
-            fieldName,
-            recTime.timePath(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE
-        ),
-        mesh
-    );
-

etc/solver-list.txt

@@ -1,5 +1,6 @@
 rcfdemSolverRhoSteadyPimple/dir
 rcfdemSolverCoupledHeattransfer/dir
+rcfdemSolverForcedTracers/dir
 rStatAnalysis/dir
 rcfdemSolverBase/dir
 rtfmSolverSpecies/dir
@@ -10,3 +11,6 @@ cfdemSolverIB/dir
 cfdemSolverPisoScalar/dir
 cfdemSolverRhoPimpleChem/dir
 cfdemSolverMultiphase/dir
+rcfdemSolverRhoSteadyPimpleChem/dir
+rctfSpeciesTransport/dir
+cfdemSolverPisoFreeStreaming/dir

etc/utilities-list.txt

@@ -1,4 +1,4 @@
+displacementField/dir
 cfdemPostproc/dir
 rSmoother/dir
-rBaseMirror/rBaseMirrorScalar/dir
-rBaseMirror/rBaseMirrorVec/dir
+rBaseMirror/rBaseMirror/dir

src/lagrangian/cfdemParticle/Make/files

@@ -37,12 +37,13 @@ $(chemistryModels)/noChemistry/noChemistry.C
 $(chemistryModels)/diffusionCoefficients/diffusionCoefficients.C
 $(chemistryModels)/massTransferCoeff/massTransferCoeff.C
 $(chemistryModels)/reactantPerParticle/reactantPerParticle.C
+$(chemistryModels)/initMultiLayers/initMultiLayers.C
 
 $(energyModels)/energyModel/energyModel.C
 $(energyModels)/energyModel/newEnergyModel.C
 $(energyModels)/heatTransferGunn/heatTransferGunn.C
 $(energyModels)/heatTransferRanzMarshall/heatTransferRanzMarshall.C
-$(energyModels)/heatTransferGranConduction/heatTransferGranConduction.C
+$(energyModels)/heatTransferInterGrain/heatTransferInterGrain.C
 $(energyModels)/reactionHeat/reactionHeat.C
 
 $(thermCondModels)/thermCondModel/thermCondModel.C
@@ -92,6 +93,8 @@ $(forceModels)/surfaceTensionForce/surfaceTensionForce.C
 $(forceModels)/gradPForceSmooth/gradPForceSmooth.C
 $(forceModels)/particleDeformation/particleDeformation.C
 $(forceModels)/turbulentDispersion/turbulentDispersion.C
+$(forceModels)/terminalVelocity/terminalVelocity.C
+$(forceModels)/turbulentVelocityFluctuations/turbulentVelocityFluctuations.C
 
 $(forceModelsMS)/forceModelMS/forceModelMS.C
 $(forceModelsMS)/forceModelMS/newForceModelMS.C
@@ -108,6 +111,7 @@ $(otherForceModels)/otherForceModel/newOtherForceModel.C
 $(otherForceModels)/gravity/gravity.C
 $(otherForceModels)/weightSecondaryPhase/weightSecondaryPhase.C
 $(otherForceModels)/expParticleForces/expParticleForces.C
+$(otherForceModels)/secondaryPhaseInducedBuoyancy/secondaryPhaseInducedBuoyancy.C
 
 $(probeModels)/probeModel/probeModel.C
 $(probeModels)/probeModel/newProbeModel.C

src/lagrangian/cfdemParticle/cfdTools/checkModelType.H

@@ -1,8 +1,10 @@
+// check model type
+// referring to: Zhou et al. 2010,JFM
 
-    // check model type
-    // referring to: Zhou et al. 2010,JFM
-    word modelType = particleCloud.modelType();
+word modelType = particleCloud.modelType();
 
+if (particleCloud.modelCheck())
+{
     //Warning << "model type not being checked" << endl;
     if (modelType=="Bfull"){
         Info << "solving volume averaged Navier Stokes equations of type B\n"<< endl;
@@ -102,3 +104,8 @@
 
     if (particleCloud.smoothingM().type() == "temporalSmoothing")
         FatalError << "the temporalSmoothing model does not support smoothing of the exchange fields, please see documentation!" << endl;
+}
+else
+{
+    Warning << "Model check deactivated." << endl;
+}

src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.C

@@ -83,11 +83,14 @@ cfdemCloud::cfdemCloud
     verbose_(couplingProperties_.found("verbose")),
     ignore_(couplingProperties_.found("ignore")),
     allowCFDsubTimestep_(true),
+    modelCheck_(couplingProperties_.lookupOrDefault<bool>("modelCheck",true)),
     limitDEMForces_(couplingProperties_.found("limitDEMForces")),
     phaseInForces_(couplingProperties_.found("phaseInForcesTime")),
     getParticleDensities_(couplingProperties_.lookupOrDefault<bool>("getParticleDensities",false)),
     getParticleEffVolFactors_(couplingProperties_.lookupOrDefault<bool>("getParticleEffVolFactors",false)),
     getParticleTypes_(couplingProperties_.lookupOrDefault<bool>("getParticleTypes",false)),
+    streamingMode_(couplingProperties_.lookupOrDefault<bool>("streamingMode",false)),
+    streamingFluc_(couplingProperties_.lookupOrDefault<bool>("streamingFluc",false)),
     maxDEMForce_(0.),
     phaseInForcesTime_(couplingProperties_.lookupOrDefault<scalar>("phaseInForcesTime",0.0)),
     modelType_(couplingProperties_.lookup("modelType")),
@@ -438,17 +441,28 @@ void cfdemCloud::getDEMdata()
 
 void cfdemCloud::giveDEMdata()
 {
-    if(forceM(0).coupleForce())
+    if (forceM(0).coupleForce())
+    {
+        if (!streamingMode_)
         {
             dataExchangeM().giveData("dragforce","vector-atom",DEMForces_);
 
-        if(impDEMdrag_)
+            if (impDEMdrag_)
             {
-            if(verbose_) Info << "sending Ksl and uf" << endl;
+                if (verbose_) Info << "sending Ksl and uf" << endl;
                 dataExchangeM().giveData("Ksl","scalar-atom",Cds_);
                 dataExchangeM().giveData("uf","vector-atom",fluidVel_);
             }
         }
+        else
+        {
+            dataExchangeM().giveData("vrec","vector-atom",particleConvVel_);
+            if (streamingFluc_)
+            {
+                dataExchangeM().giveData("vfluc","vector-atom",particleFlucVel_);
+            }
+        }
+    }
     if(verbose_) Info << "giveDEMdata done." << endl;
 }
 

src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.H

@@ -95,6 +95,8 @@ protected:
 
     bool allowCFDsubTimestep_;
 
+    const bool modelCheck_;
+
     const bool limitDEMForces_;
 
     const bool phaseInForces_;
@@ -105,6 +107,10 @@ protected:
 
     const bool getParticleTypes_;
 
+    const bool streamingMode_;
+
+    const bool streamingFluc_;
+
     scalar maxDEMForce_;
 
     scalar phaseInForcesTime_;
@@ -291,6 +297,8 @@ public:
 
        label liggghtsCommandModelIndex(const word&) const;
 
+       inline void scaleForce(int, scalar);
+
        inline void setCG(double);
 
        inline scalar cg() const;
@@ -305,6 +313,8 @@ public:
 
        inline bool ignore() const;
 
+       inline bool modelCheck() const;
+
        inline const fvMesh& mesh() const;
 
        inline bool solveFlow() const;

src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloudI.H

@@ -89,6 +89,15 @@ inline bool cfdemCloud::treatVoidCellsAsExplicitForce() const
     return treatVoidCellsAsExplicitForce_;
 }
 
+inline void cfdemCloud::scaleForce(int index, scalar factor)
+{
+    for (int i=0; i<3; ++i)
+    {
+        expForces_[index][i] *= factor;
+        impForces_[index][i] *= factor;
+    }
+}
+
 inline scalar cfdemCloud::cg() const
 {
     return cg_;
@@ -99,6 +108,11 @@ inline bool cfdemCloud::ignore() const
     return ignore_;
 }
 
+inline bool cfdemCloud::modelCheck() const
+{
+    return modelCheck_;
+}
+
 inline const fvMesh& cfdemCloud::mesh() const
 {
     return mesh_;

src/lagrangian/cfdemParticle/subModels/chemistryModel/initMultiLayers/initMultiLayers.C

@@ -0,0 +1,306 @@
+/*---------------------------------------------------------------------------*\
+License
+
+    This is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    This code is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this code.  If not, see <http://www.gnu.org/licenses/>.
+
+    Copyright (C) 2015- Thomas Lichtenegger, JKU Linz, Austria
+
+\*---------------------------------------------------------------------------*/
+
+#include "error.H"
+#include "initMultiLayers.H"
+#include "addToRunTimeSelectionTable.H"
+#include "dataExchangeModel.H"
+#include <sys/stat.h>
+#include <unistd.h>
+#include <set>
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(initMultiLayers, 0);
+
+addToRunTimeSelectionTable
+(
+        chemistryModel,
+        initMultiLayers,
+        dictionary
+);
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from components
+initMultiLayers::initMultiLayers
+(
+    const dictionary& dict,
+    cfdemCloudEnergy& sm
+)
+:
+    chemistryModel(dict,sm),
+    propsDict_(dict.subDict(typeName + "Props")),
+    interpolation_(propsDict_.lookupOrDefault<bool>("interpolation",false)),
+    mesh_(sm.mesh()),
+    verbose_(propsDict_.lookupOrDefault<bool>("verbose",false)),
+    maxNumLayers_(0),
+    numLayers_(propsDict_.lookupOrDefault<labelList>("numLayers",labelList(1,-1))),
+    partTypes_(propsDict_.lookupOrDefault<labelList>("partTypes",labelList(1,-1))),
+    relRadiiRegName_(typeName + "relRadii"),
+    filepath_(string(propsDict_.lookup("filepath"))),
+    initialized_(false)
+{
+    for (label i=0; i<numLayers_.size(); i++)
+    {
+        if (numLayers_[i] > maxNumLayers_) maxNumLayers_=numLayers_[i];
+    }
+    if (maxNumLayers_ > 4)
+    {
+        FatalError<< "Currently, not more than four layers are supported." << abort(FatalError);
+    }
+    defaultRelRadii_ = vector(0.998, 0.995, 0.98);
+    particleCloud_.checkCG(false);
+    particleCloud_.registerParticleProperty<double**>(relRadiiRegName_,maxNumLayers_);
+
+}
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+initMultiLayers::~initMultiLayers()
+{
+}
+
+// * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
+
+bool initMultiLayers::init()
+{
+    double**& relRadii = particleCloud_.getParticlePropertyRef<double**>(relRadiiRegName_);
+    // for all types of particles
+    for(label types=0; types<partTypes_.size(); types++)
+    {
+        label type = partTypes_[types];
+
+        labelList indices;
+        vectorList positions, relradii;
+
+        std::stringstream ss;
+        ss << filepath_ << type;
+        std::string filename = ss.str();
+
+        if (access( filename.c_str(), F_OK ) == -1)
+        {
+            Info << "initMultiLayers: Data file " << filename << "not found. No initialisation of layer radii possible." << endl;
+            return false;
+        }
+        Info << "\nReading" << endl;
+        Info << "\t" << filename << endl;
+        readDump(filename, type, indices, positions, relradii);
+
+        Info << "Binning particle displacements on mesh for type " << type << endl;
+
+        volScalarField particlesInCell
+        (
+            IOobject
+            (
+                "particlesInCell",
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh_,
+            dimensionedScalar("zero", dimensionSet(0,0,0,0,0), 0.0)
+        );
+
+        vector position, relrad;
+        ss.str("");
+        ss << "relRadiiField" << type;
+        std::string fieldname = ss.str();
+        volVectorField relRadiiField
+        (
+            IOobject
+            (
+                fieldname,
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh_,
+            dimensionedVector("zero", dimensionSet(0,0,0,0,0), vector::zero)
+        );
+
+        label cellI = -1;
+        for (label lineI = 0; lineI < indices.size(); lineI++)
+        {
+            position = positions[lineI];
+            relrad = relradii[lineI];
+            cellI = mesh_.findCell(position);
+            if (cellI < 0)  continue;
+            particlesInCell[cellI] += 1;
+            relRadiiField[cellI] += relrad;
+        }
+
+        forAll (mesh_.C(), cellJ)
+        {
+            if (particlesInCell[cellJ] > 0.5)
+            {
+                relRadiiField[cellJ] /= particlesInCell[cellJ];
+            }
+        }
+
+        // fill field
+
+
+        // fill particle arrays
+        label cellK = -1;
+        label cellKoccupied = -1;
+        label partType = -1;
+        label listIndex = -1;
+        vector relRadiiLoopVar = vector::zero;
+        for(int index = 0;index < particleCloud_.numberOfParticles(); ++index)
+        {
+            cellK = particleCloud_.cellIDs()[index][0];
+            partType = particleCloud_.particleType(index);
+            if(cellK >= 0 && partType == type)
+            {
+                // look for the nearest occupied cell
+                cellKoccupied = findNearestCellWithValue(cellK, particlesInCell);
+                listIndex = getListIndex(partType);
+                if (cellKoccupied >= 0)
+                {
+                    relRadiiLoopVar = relRadiiField[cellKoccupied];
+                }
+                else
+                {
+                    relRadiiLoopVar = defaultRelRadii_;
+                }
+                relRadii[index][0] = 1.0;
+                for (label i=1;i<numLayers_[listIndex];i++)
+                {
+                    relRadii[index][i] = relRadiiLoopVar.component(i-1);
+                }
+            }
+        }
+        if (verbose_) relRadiiField.write();
+    }
+    return true;
+}
+
+// * * * * * * * * * * * * * * * * Member Fct  * * * * * * * * * * * * * * * //
+
+void initMultiLayers::execute()
+{
+    if(!initialized_)
+    {
+        if (init())
+        {
+            double**& relRadii = particleCloud_.getParticlePropertyRef<double**>(relRadiiRegName_);
+            particleCloud_.dataExchangeM().giveData("relRadii","vector-atom", relRadii);
+        }
+        initialized_ = true;
+    }
+}
+
+label initMultiLayers::findNearestCellWithValue(label refCell, volScalarField &particlesInCell) const
+{
+    if (particlesInCell[refCell] > 0.5) return refCell;
+
+    std::set<label> neighbors;
+    std::set<label> newNeighbors;
+    std::set<label> recentNeighbors;
+
+    neighbors.insert(refCell);
+    recentNeighbors.insert(refCell);
+
+    while(true)
+    {
+        for (std::set<label>::iterator it=recentNeighbors.begin(); it!=recentNeighbors.end(); ++it)
+        {
+            labelList adjacent = mesh_.cellCells()[*it];
+            label adj;
+            for (label j=0; j<adjacent.size(); j++)
+            {
+                adj = adjacent[j];
+                std::set<label>::iterator it2 = neighbors.find(adj);
+                if (it2 == neighbors.end())
+                {
+                    newNeighbors.insert(adj);
+                    neighbors.insert(adj);
+                    if (particlesInCell[adj] > 0.5) return adj;
+                }
+            }
+        }
+
+        // if all cells have been searched and none was occupied, return -1; assumes that reasonable default value is available
+        if (newNeighbors.size() == 0) return -1;
+        recentNeighbors.clear();
+        recentNeighbors = newNeighbors;
+        newNeighbors.clear();
+    }
+}
+
+label initMultiLayers::getListIndex(label testElement) const
+{
+    for(label ind = 0; ind<partTypes_.size(); ind++)
+    {
+        if (partTypes_[ind] == testElement) return ind;
+    }
+    // testing
+    Pout << "Cannot find list index for element " << testElement << endl;
+    return -1;
+}
+
+void initMultiLayers::readDump(std::string filename, label type, labelList &indices, vectorList &positions, vectorList &relradii)
+{
+    #include <fstream>
+
+    const label leadingLines = 9;
+    label lineCounter = 0;
+    label partIndex, partType;
+    scalar x, y, z;
+    scalar r0 = 1.0;
+    scalar r1 = 1.0;
+    scalar r2 = 1.0;
+    scalar r3 = 1.0;
+
+    indices.clear();
+    positions.clear();
+    relradii.clear();
+
+    std::ifstream file(filename);
+    std::string str;
+    while (std::getline(file, str))
+    {
+        if (lineCounter >= leadingLines)
+        {
+            sscanf(str.c_str(), "%d %d %lf %lf %lf %lf %lf %lf %lf", &partIndex, &partType, &x, &y, &z, &r0, &r1, &r2, &r3);
+            if (partType != type)
+            {
+                FatalError<< "Particle of type " << partType << " detected in " << filename << abort(FatalError);
+            }
+            indices.append(partIndex);
+            positions.append(vector(x,y,z));
+            relradii.append(vector(r1,r2,r3));
+        }
+        lineCounter++;
+    }
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/lagrangian/cfdemParticle/subModels/chemistryModel/initMultiLayers/initMultiLayers.H

@@ -0,0 +1,115 @@
+/*---------------------------------------------------------------------------*\
+License
+
+    This is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    This code is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this code.  If not, see <http://www.gnu.org/licenses/>.
+
+    Copyright (C) 2015- Thomas Lichtenegger, JKU Linz, Austria
+
+    Description
+    Model to initiate multilayer radii from dump files
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef initMultiLayers_H
+#define initMultiLayers_H
+
+#include <sstream>
+#include <string>
+
+#include "fvCFD.H"
+#include "cfdemCloudEnergy.H"
+#include "chemistryModel.H"
+#include "vectorList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class species Declaration
+\*---------------------------------------------------------------------------*/
+
+// for future use:
+// + communicate every N steps
+
+class initMultiLayers
+:
+    public chemistryModel
+{
+private:
+
+    dictionary propsDict_;
+
+    bool interpolation_;
+
+    const fvMesh& mesh_;
+
+    bool verbose_;
+
+    label maxNumLayers_;
+
+    labelList numLayers_;
+
+    labelList partTypes_;
+
+    vector defaultRelRadii_;
+
+    const word relRadiiRegName_;
+
+    std::string filepath_;
+
+    bool initialized_;
+
+    bool init();
+
+    void readDump(std::string, label, labelList &, vectorList &, vectorList &);
+
+    label findNearestCellWithValue(label, volScalarField &) const;
+
+    label getListIndex(label) const;
+
+public:
+
+
+    //- Runtime type information
+    TypeName("initMultiLayers");
+
+    // Constructors
+
+        //- Construct from components
+        initMultiLayers
+        (
+            const dictionary& dict,
+            cfdemCloudEnergy& sm
+        );
+
+    // Destructor
+
+        virtual ~initMultiLayers();
+
+    // Member Functions
+        void execute();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/cfdemParticle/subModels/chemistryModel/species/species.C

@@ -98,11 +98,17 @@ species::species
     partMolarConcName_(propsDict_.lookupOrDefault<word>("partMoleName","partMolarConc")),
     loopCounter_(-1),
     Nevery_(propsDict_.lookupOrDefault<label>("Nevery",1)),
+    couplingTimestep_(0.0),
     massSourceCurr_(0.0),
     massSourceTot_(0.0),
     initialized_(false)
 {
     particleCloud_.checkCG(false);
+
+    scalar dtDEM = particleCloud_.dataExchangeM().DEMts();
+    scalar dtCFD = mesh_.time().deltaTValue();
+    couplingTimestep_ = max(dtDEM,dtCFD);
+
     particleCloud_.registerParticleProperty<double**>(partTempName_,1);
     particleCloud_.registerParticleProperty<double**>(partRhoName_,1);
     particleCloud_.registerParticleProperty<double**>(partMolarConcName_,1);
@@ -268,7 +274,6 @@ void species::execute()
 
     // pull changeOfSpeciesMass_, transform onto fields changeOfSpeciesMassFields_, add them up on changeOfGasMassField_
     {
-        scalar timestep = mesh_.time().deltaTValue();
         changeOfGasMassField_.primitiveFieldRef() = 0.0;
         changeOfGasMassField_.boundaryFieldRef() = 0.0;
         for (int i=0; i<speciesNames_.size();i++)
@@ -291,17 +296,17 @@ void species::execute()
 
             // take care for implementation in LIGGGHTS: species produced from particles defined positive
             // changeOf...Fields need to be mass per volume per timestep
-            changeOfSpeciesMassFields_[i].primitiveFieldRef() /= (changeOfSpeciesMassFields_[i].mesh().V() * Nevery_ * timestep);
+            changeOfSpeciesMassFields_[i].primitiveFieldRef() /= (changeOfSpeciesMassFields_[i].mesh().V() * Nevery_ * couplingTimestep_);
             changeOfSpeciesMassFields_[i].correctBoundaryConditions();
             changeOfGasMassField_ += changeOfSpeciesMassFields_[i];
 
             if (verbose_)
             {
                 Info << "total conversion of species" << speciesNames_[i] << " = "
-                     << gSum(changeOfSpeciesMassFields_[i]*1.0*changeOfSpeciesMassFields_[i].mesh().V() * Nevery_ * timestep) << endl;
+                     << gSum(changeOfSpeciesMassFields_[i]*1.0*changeOfSpeciesMassFields_[i].mesh().V() * Nevery_ * couplingTimestep_) << endl;
             }
         }
-        massSourceCurr_ = gSum(changeOfGasMassField_*1.0*changeOfGasMassField_.mesh().V() * Nevery_ * timestep);
+        massSourceCurr_ = gSum(changeOfGasMassField_*1.0*changeOfGasMassField_.mesh().V() * Nevery_ * couplingTimestep_);
         massSourceTot_ += massSourceCurr_;
 
         if (verbose_)

src/lagrangian/cfdemParticle/subModels/chemistryModel/species/species.H

@@ -99,6 +99,8 @@ private:
 
     label Nevery_;
 
+    scalar couplingTimestep_;
+
     scalar massSourceCurr_;
 
     scalar massSourceTot_;

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

@@ -1,249 +0,0 @@
-/*---------------------------------------------------------------------------*\
-License
-
-    This is free software: you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by
-    the Free Software Foundation, either version 3 of the License, or
-    (at your option) any later version.
-
-    This code is distributed in the hope that it will be useful, but WITHOUT
-    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-    for more details.
-
-    You should have received a copy of the GNU General Public License
-    along with this code.  If not, see <http://www.gnu.org/licenses/>.
-
-    Copyright (C) 2015- Thomas Lichtenegger, JKU Linz, Austria
-
-\*---------------------------------------------------------------------------*/
-
-#include "error.H"
-#include "heatTransferGranConduction.H"
-#include "addToRunTimeSelectionTable.H"
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-defineTypeNameAndDebug(heatTransferGranConduction, 0);
-
-addToRunTimeSelectionTable(energyModel, heatTransferGranConduction, dictionary);
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-// Construct from components
-heatTransferGranConduction::heatTransferGranConduction
-(
-    const dictionary& dict,
-    cfdemCloudEnergy& sm
-)
-:
-    energyModel(dict,sm),
-    propsDict_(dict.subDict(typeName + "Props")),
-    multiTypes_(false),
-    verbose_(propsDict_.lookupOrDefault<bool>("verbose",false)),
-    calcTotalHeatFlux_(propsDict_.lookupOrDefault<bool>("calcTotalHeatFlux",false)),
-    totalHeatFlux_(0.0),
-    QPartPartName_(propsDict_.lookupOrDefault<word>("QPartPartName","QPartPart")),
-    QPartPart_
-    (
-        IOobject
-        (
-            QPartPartName_,
-            sm.mesh().time().timeName(),
-            sm.mesh(),
-            IOobject::READ_IF_PRESENT,
-            IOobject::NO_WRITE
-        ),
-        sm.mesh(),
-        dimensionedScalar("zero", dimensionSet(1,-1,-3,0,0,0,0), 0.0)
-    ),
-    partEffThermCondField_
-    (
-        IOobject
-        (
-            "partEffThermCondField",
-            sm.mesh().time().timeName(),
-            sm.mesh(),
-            IOobject::READ_IF_PRESENT,
-            IOobject::NO_WRITE
-        ),
-        sm.mesh(),
-        dimensionedScalar("one", dimensionSet(1, 1, -3, -1,0,0,0), 1.0),
-        "zeroGradient"
-    ),
-    partThermCondField_
-    (
-        IOobject
-        (
-            "partThermCondField",
-            sm.mesh().time().timeName(),
-            sm.mesh(),
-            IOobject::READ_IF_PRESENT,
-            IOobject::NO_WRITE
-        ),
-        sm.mesh(),
-        dimensionedScalar("one", dimensionSet(1, 1, -3, -1,0,0,0), 1.0),
-        "zeroGradient"
-    ),
-    partTempField_(sm.mesh().lookupObject<volScalarField>("partTemp")),
-    prescribedVoidfractionFieldName_(propsDict_.lookupOrDefault<word>("prescribedVoidfractionFieldName","voidfraction")),
-    prescribedVoidfraction_(sm.mesh().lookupObject<volScalarField> (prescribedVoidfractionFieldName_)),
-    voidfractionFieldName_(propsDict_.lookupOrDefault<word>("voidfractionFieldName","voidfraction")),
-    voidfraction_(sm.mesh().lookupObject<volScalarField> (voidfractionFieldName_)),
-    partHeatFluxName_(propsDict_.lookupOrDefault<word>("partHeatFluxName","conductiveHeatFlux")),
-    typePartThermCond_(propsDict_.lookupOrDefault<scalarList>("thermalConductivities",scalarList(1,-1.0))),
-    partThermCondRegName_(typeName + "partThermCond")
-{
-    particleCloud_.registerParticleProperty<double**>(partHeatFluxName_,1);
-    particleCloud_.registerParticleProperty<double**>(partThermCondRegName_,1);
-
-    if (typePartThermCond_[0] < 0.0)
-    {
-        FatalError << "heatTransferGranConduction: provide list of thermal conductivities." << abort(FatalError);
-    }
-
-    if (typePartThermCond_.size() > 1)
-    {
-        multiTypes_ = true;
-    }
-
-    if (multiTypes_ && !particleCloud_.getParticleTypes())
-    {
-        FatalError << "heatTransferGranConduction needs data for more than one type, but types are not communicated." << abort(FatalError);
-    }
-
-    if (verbose_)
-    {
-        QPartPart_.writeOpt() = IOobject::AUTO_WRITE;
-        QPartPart_.write();
-    }
-}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-heatTransferGranConduction::~heatTransferGranConduction()
-{
-}
-
-// * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
-// * * * * * * * * * * * * * * * * Member Fct  * * * * * * * * * * * * * * * //
-
-void heatTransferGranConduction::calcEnergyContribution()
-{
-    calcPartEffThermCond();
-
-    QPartPart_ = fvc::laplacian(partEffThermCondField_,partTempField_);
-
-    label cellI=0;
-    scalar partVolume(0);
-    scalar QPartPart(0);
-    scalar voidfraction(1);
-
-    totalHeatFlux_ = 0.0;
-    double**& partHeatFlux_ = particleCloud_.getParticlePropertyRef<double**>(partHeatFluxName_);
-
-    for(int index = 0;index < particleCloud_.numberOfParticles(); ++index)
-    {
-        cellI = particleCloud_.cellIDs()[index][0];
-        if(cellI >= 0)
-        {
-            voidfraction = voidfraction_[cellI];
-            if (voidfraction < 0.01)
-                voidfraction = 0.01;
-
-            partVolume = particleCloud_.particleVolume(index);
-            QPartPart = QPartPart_[cellI];
-
-            heatFlux(index, partVolume, voidfraction, QPartPart);
-            if (calcTotalHeatFlux_) totalHeatFlux_ += partHeatFlux_[index][0];
-        }
-    }
-}
-
-void heatTransferGranConduction::calcPartEffThermCond()
-{
-    calcPartThermCond();
-
-    scalar volFrac = 0.0;
-
-    forAll(partEffThermCondField_, cellI)
-    {
-        volFrac = 1.0 - prescribedVoidfraction_[cellI];
-        if (volFrac < 0.334)
-        {
-            partEffThermCondField_[cellI] = 0.0;
-        }
-        else
-        {
-            partEffThermCondField_[cellI] = 0.5 * (3*volFrac - 1) * partThermCondField_[cellI];
-        }
-    }
-}
-
-void heatTransferGranConduction::calcPartThermCond()
-{
-    label cellI=0;
-    label partType = 1;
-    double**& partThermCond_ = particleCloud_.getParticlePropertyRef<double**>(partThermCondRegName_);
-
-    for(int index = 0;index < particleCloud_.numberOfParticles(); ++index)
-    {
-            cellI = particleCloud_.cellIDs()[index][0];
-            if (cellI >= 0)
-            {
-                if (multiTypes_)
-                {
-                    partType = particleCloud_.particleType(index);
-                }
-                // LIGGGGHTS counts types 1, 2, ..., C++ array starts at 0
-                partThermCond_[index][0] = typePartThermCond_[partType - 1];
-            }
-    }
-
-    partThermCondField_.primitiveFieldRef() = 0.0;
-    particleCloud_.averagingM().resetWeightFields();
-    particleCloud_.averagingM().setScalarAverage
-    (
-        partThermCondField_,
-        partThermCond_,
-        particleCloud_.particleWeights(),
-        particleCloud_.averagingM().UsWeightField(),
-        NULL
-    );
-}
-
-void heatTransferGranConduction::heatFlux(label index, scalar vol, scalar voidfraction, scalar QPartPart)
-{
-    double**& partHeatFlux_ = particleCloud_.getParticlePropertyRef<double**>(partHeatFluxName_);
-    partHeatFlux_[index][0] = vol * QPartPart / (1.0 - voidfraction) ;
-}
-
-void heatTransferGranConduction::giveData()
-{
-    if (calcTotalHeatFlux_)
-    {
-        reduce(totalHeatFlux_, sumOp<scalar>());
-        Info << "total conductive particle-particle heat flux [W] (Eulerian) = " << totalHeatFlux_ << endl;
-    }
-
-    double**& partHeatFlux_ = particleCloud_.getParticlePropertyRef<double**>(partHeatFluxName_);
-    particleCloud_.dataExchangeM().giveData(partHeatFluxName_,"scalar-atom", partHeatFlux_);
-}
-
-void heatTransferGranConduction::postFlow()
-{
-    giveData();
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// ************************************************************************* //
-

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

@@ -0,0 +1,468 @@
+
+/*---------------------------------------------------------------------------*\
+License
+
+    This is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    This code is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this code.  If not, see <http://www.gnu.org/licenses/>.
+
+    Copyright (C) 2015- Thomas Lichtenegger, JKU Linz, Austria
+
+\*---------------------------------------------------------------------------*/
+
+#include "error.H"
+#include "heatTransferInterGrain.H"
+#include "addToRunTimeSelectionTable.H"
+#include "fvCFD.H"
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(heatTransferInterGrain, 0);
+
+addToRunTimeSelectionTable(energyModel, heatTransferInterGrain, dictionary);
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from components
+heatTransferInterGrain::heatTransferInterGrain
+(
+    const dictionary& dict,
+    cfdemCloudEnergy& sm
+)
+:
+    energyModel(dict,sm),
+    propsDict_(dict.subDict(typeName + "Props")),
+    multiTypes_(false),
+    verbose_(propsDict_.lookupOrDefault<bool>("verbose",false)),
+    implicit_(propsDict_.lookupOrDefault<bool>("implicit",false)),
+    calcTotalHeatFlux_(propsDict_.lookupOrDefault<bool>("calcTotalHeatFlux",false)),
+    radiativeHeatTransfer_(propsDict_.lookupOrDefault<bool>("radiativeHeatTransfer",false)),
+    totalHeatFlux_(0.0),
+    partTempName_(""),
+    QPartPartName_(propsDict_.lookupOrDefault<word>("QPartPartName","QPartPart")),
+    QPartPart_
+    (
+        IOobject
+        (
+            QPartPartName_,
+            sm.mesh().time().timeName(),
+            sm.mesh(),
+            IOobject::READ_IF_PRESENT,
+            IOobject::NO_WRITE
+        ),
+        sm.mesh(),
+        dimensionedScalar("zero", dimensionSet(1,-1,-3,0,0,0,0), 0.0)
+    ),
+    partEffThermCondField_
+    (
+        IOobject
+        (
+            "partEffThermCondField",
+            sm.mesh().time().timeName(),
+            sm.mesh(),
+            IOobject::READ_IF_PRESENT,
+            IOobject::NO_WRITE
+        ),
+        sm.mesh(),
+        dimensionedScalar("one", dimensionSet(1, 1, -3, -1,0,0,0), 1.0),
+        "zeroGradient"
+    ),
+    partThermCondField_
+    (
+        IOobject
+        (
+            "partThermCondField",
+            sm.mesh().time().timeName(),
+            sm.mesh(),
+            IOobject::READ_IF_PRESENT,
+            IOobject::NO_WRITE
+        ),
+        sm.mesh(),
+        dimensionedScalar("one", dimensionSet(1, 1, -3, -1,0,0,0), 1.0),
+        "zeroGradient"
+    ),
+    partThermCapField_
+    (
+        IOobject
+        (
+            "partThermCapField",
+            sm.mesh().time().timeName(),
+            sm.mesh(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        sm.mesh(),
+        dimensionedScalar("zero", dimensionSet(1,-1,-2,-1,0,0,0), 0.0),
+        "zeroGradient"
+    ),
+    partThermRadField_
+    (
+        IOobject
+        (
+            "partThermRadField",
+            sm.mesh().time().timeName(),
+            sm.mesh(),
+            IOobject::READ_IF_PRESENT,
+            IOobject::NO_WRITE
+        ),
+        sm.mesh(),
+        dimensionedScalar("zero", dimensionSet(1, 1, -3, -1,0,0,0), 0.0),
+        "zeroGradient"
+    ),
+    partTempField_(sm.mesh().lookupObject<volScalarField>("partTemp")),
+    prescribedVoidfractionFieldName_(propsDict_.lookupOrDefault<word>("prescribedVoidfractionFieldName","voidfraction")),
+    prescribedVoidfraction_(sm.mesh().lookupObject<volScalarField> (prescribedVoidfractionFieldName_)),
+    voidfractionFieldName_(propsDict_.lookupOrDefault<word>("voidfractionFieldName","voidfraction")),
+    voidfraction_(sm.mesh().lookupObject<volScalarField> (voidfractionFieldName_)),
+    typeCG_(propsDict_.lookupOrDefault<scalarList>("coarseGrainingFactors",scalarList(1,1.0))),
+    maxTypeCG_(typeCG_.size()),
+    partHeatFluxName_(propsDict_.lookupOrDefault<word>("partHeatFluxName","conductiveHeatFlux")),
+    typePartThermCond_(propsDict_.lookupOrDefault<scalarList>("thermalConductivities",scalarList(1,-1.0))),
+    partThermCondRegName_(typeName + "partThermCond"),
+    typePartThermCap_(propsDict_.lookupOrDefault<scalarList>("thermalCapacities",scalarList(1,-1.0))),
+    partThermCapRegName_(typeName + "partThermCap"),
+    partThermRadRegName_(typeName + "partThermRad"),
+    StefanBoltzmannConst_(5.67e-8),
+    typePartEmissivity_(propsDict_.lookupOrDefault<scalarList>("thermalEmissivities",scalarList(1,-1.0))),
+    kMax_(propsDict_.lookupOrDefault<scalar>("kMax",-1.0))
+
+{
+    particleCloud_.registerParticleProperty<double**>(partHeatFluxName_,1);
+    particleCloud_.registerParticleProperty<double**>(partThermCondRegName_,1);
+    if (radiativeHeatTransfer_)
+    {
+        particleCloud_.registerParticleProperty<double**>(partThermRadRegName_,1);
+        partTempName_ = word(propsDict_.lookup("partTempName"));
+    }
+
+    if (typePartThermCond_[0] < 0.0)
+    {
+        FatalError << "heatTransferInterGrain: provide list of thermal conductivities." << abort(FatalError);
+    }
+
+    if (typePartEmissivity_[0] < 0.0 && radiativeHeatTransfer_)
+    {
+        FatalError << "heatTransferGranRadiation: provide list of thermal emissivities." << abort(FatalError);
+    }
+
+    if (typePartEmissivity_[0] > 0.0 && !radiativeHeatTransfer_)
+    {
+        FatalError << "heatTransferGranRadiation: thermal emissivities provided but calculation deactivated." << abort(FatalError);
+    }
+
+    if (typePartThermCond_.size() > 1)
+    {
+        multiTypes_ = true;
+    }
+
+    if (multiTypes_ && !particleCloud_.getParticleTypes())
+    {
+        FatalError << "heatTransferInterGrain needs data for more than one type, but types are not communicated." << abort(FatalError);
+    }
+
+    if (verbose_)
+    {
+        QPartPart_.writeOpt() = IOobject::AUTO_WRITE;
+        partEffThermCondField_.writeOpt() = IOobject::AUTO_WRITE;
+        partThermCondField_.writeOpt() = IOobject::AUTO_WRITE;
+        partThermCapField_.writeOpt() = IOobject::AUTO_WRITE;
+        partThermRadField_.writeOpt() = IOobject::AUTO_WRITE;
+        QPartPart_.write();
+        partEffThermCondField_.write();
+        partThermCondField_.write();
+        partThermCapField_.write();
+        partThermRadField_.write();
+    }
+
+    if (implicit_)
+    {
+        particleCloud_.registerParticleProperty<double**>(partThermCapRegName_,1);
+        if (typePartThermCap_[0] < 0.0)
+        {
+            FatalError << "heatTransferInterGrain: provide list of thermal capacities." << abort(FatalError);
+        }
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+heatTransferInterGrain::~heatTransferInterGrain()
+{
+}
+
+// * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
+// * * * * * * * * * * * * * * * * Member Fct  * * * * * * * * * * * * * * * //
+
+void heatTransferInterGrain::calcEnergyContribution()
+{
+    calcPartEffThermCond();
+
+    if (kMax_ > 0.0)
+    {
+        dimensionedScalar kMax("kMax",dimensionSet(1,1,-3,-1,0,0,0),kMax_);
+        partEffThermCondField_ = min(partEffThermCondField_,kMax);
+    }
+
+    if (!implicit_)
+    {
+        QPartPart_ = fvc::laplacian(partEffThermCondField_,partTempField_);
+    }
+    else
+    {
+        double**& partThermCap_ = particleCloud_.getParticlePropertyRef<double**>(partThermCapRegName_);
+        label cellI = -1;
+        label partType = -1;
+        for(int index = 0;index < particleCloud_.numberOfParticles(); ++index)
+        {
+            cellI = particleCloud_.cellIDs()[index][0];
+            if (cellI >= 0)
+            {
+                if (multiTypes_)
+                {
+                    partType = particleCloud_.particleType(index);
+                }
+                // LIGGGGHTS counts types 1, 2, ..., C++ array starts at 0
+                partThermCap_[index][0] = typePartThermCap_[partType - 1] * particleCloud_.particleDensity(index);
+            }
+        }
+
+        partThermCapField_.primitiveFieldRef() = 0.0;
+        particleCloud_.averagingM().resetWeightFields();
+        particleCloud_.averagingM().setScalarAverage
+        (
+            partThermCapField_,
+            partThermCap_,
+            particleCloud_.particleWeights(),
+            particleCloud_.averagingM().UsWeightField(),
+            NULL
+        );
+
+        partThermCapField_ *= (1.0 - voidfraction_);
+        dimensionedScalar CMin(dimensionedScalar("CMin",dimensionSet(1,-1,-2,-1,0,0,0),SMALL));
+        partThermCapField_ = max(partThermCapField_, CMin);
+        partThermCapField_.oldTime() = partThermCapField_;
+
+        volScalarField partTempField(partTempField_);
+        partTempField.oldTime() = partTempField;
+        fvScalarMatrix TpEqn
+        (
+            fvm::ddt(partThermCapField_,partTempField)
+            ==
+            fvm::laplacian(partEffThermCondField_,partTempField)
+        );
+
+        TpEqn.solve();
+        QPartPart_ = fvc::laplacian(partEffThermCondField_,partTempField);
+    }
+
+    label cellI=0;
+    scalar partVolume(0);
+    scalar QPartPart(0);
+    scalar voidfraction(1);
+
+    totalHeatFlux_ = 0.0;
+    double**& partHeatFlux_ = particleCloud_.getParticlePropertyRef<double**>(partHeatFluxName_);
+
+    for(int index = 0;index < particleCloud_.numberOfParticles(); ++index)
+    {
+        cellI = particleCloud_.cellIDs()[index][0];
+        if(cellI >= 0)
+        {
+            voidfraction = voidfraction_[cellI];
+            if (voidfraction < 0.01)
+                voidfraction = 0.01;
+
+            partVolume = particleCloud_.particleVolume(index);
+            QPartPart = QPartPart_[cellI];
+
+            heatFlux(index, partVolume, voidfraction, QPartPart);
+            if (calcTotalHeatFlux_) totalHeatFlux_ += partHeatFlux_[index][0];
+        }
+    }
+}
+
+void heatTransferInterGrain::calcPartEffThermCond()
+{
+    calcPartThermCond();
+    if (radiativeHeatTransfer_) calcPartThermRad();
+
+    scalar volFrac = 0.0;
+
+    if (!radiativeHeatTransfer_)
+    {
+        forAll(partEffThermCondField_, cellI)
+        {
+            volFrac = 1.0 - prescribedVoidfraction_[cellI];
+            if (volFrac < 0.334)
+            {
+                partEffThermCondField_[cellI] = 0.0;
+            }
+            else
+            {
+                partEffThermCondField_[cellI] = 0.5 * (3*volFrac - 1) * partThermCondField_[cellI];
+            }
+        }
+    }
+    else
+    {
+        scalar thRad = 0.0;
+        forAll(partEffThermCondField_, cellI)
+        {
+            volFrac = 1.0 - prescribedVoidfraction_[cellI];
+            thRad = volFrac * partThermRadField_[cellI]; // volume fraction factor for rad. contribution according to Qian et al.
+            if (volFrac < 0.334)
+            {
+                partEffThermCondField_[cellI] = thRad;
+            }
+            else
+            {
+                partEffThermCondField_[cellI] = 0.5 * (3*volFrac - 1) * partThermCondField_[cellI] + thRad;
+            }
+        }
+    }
+}
+
+void heatTransferInterGrain::calcPartThermCond()
+{
+    label cellI=0;
+    label partType = 1;
+    double**& partThermCond_ = particleCloud_.getParticlePropertyRef<double**>(partThermCondRegName_);
+
+    for(int index = 0;index < particleCloud_.numberOfParticles(); ++index)
+    {
+            cellI = particleCloud_.cellIDs()[index][0];
+            if (cellI >= 0)
+            {
+                if (multiTypes_)
+                {
+                    partType = particleCloud_.particleType(index);
+                }
+                // LIGGGGHTS counts types 1, 2, ..., C++ array starts at 0
+                partThermCond_[index][0] = typePartThermCond_[partType - 1];
+            }
+    }
+
+    partThermCondField_.primitiveFieldRef() = 0.0;
+    particleCloud_.averagingM().resetWeightFields();
+    particleCloud_.averagingM().setScalarAverage
+    (
+        partThermCondField_,
+        partThermCond_,
+        particleCloud_.particleWeights(),
+        particleCloud_.averagingM().UsWeightField(),
+        NULL
+    );
+}
+
+void heatTransferInterGrain::calcPartThermRad()
+{
+    label cellI=0;
+    label partType = 1;
+    scalar cg = typeCG_[0];
+    scalar ds = 0.0;
+    scalar L = 0.0;
+    scalar prefac = 0.0;
+    scalar Tp = 0.0;
+    scalar voidfraction = 0.0;
+
+    double**& partTemp_ = particleCloud_.getParticlePropertyRef<double**>(partTempName_);
+    double**& partThermCond_ = particleCloud_.getParticlePropertyRef<double**>(partThermCondRegName_);
+    double**& partThermRad_ = particleCloud_.getParticlePropertyRef<double**>(partThermRadRegName_);
+
+    for(int index = 0;index < particleCloud_.numberOfParticles(); ++index)
+    {
+            cellI = particleCloud_.cellIDs()[index][0];
+            if (cellI >= 0)
+            {
+                if (multiTypes_)
+                {
+                    partType = particleCloud_.particleType(index);
+                    if (partType > maxTypeCG_)
+                    {
+                        FatalError<< "Too few coarse-graining factors provided." << abort(FatalError);
+                    }
+                    cg = typeCG_[partType - 1];
+                }
+                ds = 2.*particleCloud_.radius(index)/cg;
+                // make sure reasonable values are used
+                Tp = partTemp_[index][0];
+                if (Tp < TMin) Tp = TMin;
+                voidfraction = prescribedVoidfraction_[cellI];
+                if (voidfraction < voidfracMin) voidfraction = voidfracMin;
+                else if (voidfraction > voidfracMax) voidfraction = voidfracMax;
+
+                prefac = 4.0*StefanBoltzmannConst_*ds*Tp*Tp*Tp;
+                L = partThermCond_[index][0]/prefac;
+                // LIGGGGHTS counts types 1, 2, ..., C++ array starts at 0
+                partThermRad_[index][0] = prefac*FE(voidfraction,typePartEmissivity_[partType - 1],L);
+            }
+    }
+
+    partThermRadField_.primitiveFieldRef() = 0.0;
+    particleCloud_.averagingM().resetWeightFields();
+    particleCloud_.averagingM().setScalarAverage
+    (
+        partThermRadField_,
+        partThermRad_,
+        particleCloud_.particleWeights(),
+        particleCloud_.averagingM().UsWeightField(),
+        NULL
+    );
+}
+
+scalar heatTransferInterGrain::FE(scalar voidfraction, scalar emissivity, scalar L)
+{
+    if (voidfraction < 1 - SMALL) return 1.0; // limit of voidfraction -> 1.0
+
+    scalar B = 1.25 * pow((1-voidfraction)/voidfraction,1.111);
+    scalar sqrt_aP = sqrt(1-voidfraction);
+    scalar x = 2.0/emissivity - 1.0;
+    scalar y = 1.0 / (x*L);
+    scalar z = sqrt_aP / (2.0/emissivity - 1.0) * (B+1.0)/B * 1.0/(1.0 + y);
+    return (1 - sqrt_aP)*voidfraction + z;
+}
+
+void heatTransferInterGrain::heatFlux(label index, scalar vol, scalar voidfraction, scalar QPartPart)
+{
+    double**& partHeatFlux_ = particleCloud_.getParticlePropertyRef<double**>(partHeatFluxName_);
+    partHeatFlux_[index][0] = vol * QPartPart / (1.0 - voidfraction) ;
+}
+
+void heatTransferInterGrain::giveData()
+{
+    if (calcTotalHeatFlux_)
+    {
+        reduce(totalHeatFlux_, sumOp<scalar>());
+        Info << "total conductive particle-particle heat flux [W] (Eulerian) = " << totalHeatFlux_ << endl;
+    }
+
+    double**& partHeatFlux_ = particleCloud_.getParticlePropertyRef<double**>(partHeatFluxName_);
+    particleCloud_.dataExchangeM().giveData(partHeatFluxName_,"scalar-atom", partHeatFlux_);
+}
+
+void heatTransferInterGrain::postFlow()
+{
+    giveData();
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //
+

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

@@ -19,12 +19,16 @@ License
     Description
     heat transfer from granular conduction according to the
     Effective Medium Theory discussed e.g. by
-    Carson et al. Int. J. Heat Mass Transfer 48 (2005)
+    Carson et al. Int. J. Heat Mass Transfer 48 (2005) 2150–2158
+    optional radiative contribution according to
+    G. Breitbach, H. Barthels. Nucl. Technol. 49 (1980) 392–399
+    as discussed in
+    Qian et al. Int. J. Heat Mass Transf. 127 (2018) 573–584
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef heatTransferGranConduction_H
-#define heatTransferGranConduction_H
+#ifndef heatTransferInterGrain_H
+#define heatTransferInterGrain_H
 
 #include "fvCFD.H"
 #include "cfdemCloudEnergy.H"
@@ -35,10 +39,10 @@ namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class heatTransferGranConduction Declaration
+                           Class heatTransferInterGrain Declaration
 \*---------------------------------------------------------------------------*/
 
-class heatTransferGranConduction
+class heatTransferInterGrain
 :
     public energyModel
 {
@@ -50,10 +54,16 @@ protected:
 
     bool verbose_;
 
+    bool implicit_;
+
     bool calcTotalHeatFlux_;
 
+    bool radiativeHeatTransfer_;
+
     scalar totalHeatFlux_;
 
+    word partTempName_;
+
     const word QPartPartName_;
 
     volScalarField QPartPart_;
@@ -62,6 +72,10 @@ protected:
 
     volScalarField partThermCondField_;
 
+    volScalarField partThermCapField_;
+
+    volScalarField partThermRadField_;
+
     const volScalarField& partTempField_;
 
     const word prescribedVoidfractionFieldName_;
@@ -72,16 +86,42 @@ protected:
 
     const volScalarField& voidfraction_;
 
+    scalarList typeCG_;
+
+    const label maxTypeCG_;
+
     const word partHeatFluxName_;
 
     scalarList typePartThermCond_;
 
     const word partThermCondRegName_;
 
+    scalarList typePartThermCap_;
+
+    const word partThermCapRegName_;
+
+    const word partThermRadRegName_;
+
+    const scalar TMin = 10.0;
+
+    const scalar voidfracMax = 0.95;
+
+    const scalar voidfracMin = 0.05;
+
+    const scalar StefanBoltzmannConst_;
+
+    scalarList typePartEmissivity_;
+
+    scalar kMax_;
+
     void calcPartEffThermCond();
 
     void calcPartThermCond();
 
+    void calcPartThermRad();
+
+    scalar FE(scalar, scalar, scalar);
+
     virtual void giveData();
 
     virtual void heatFlux(label, scalar, scalar, scalar);
@@ -91,12 +131,12 @@ protected:
 public:
 
     //- Runtime type information
-    TypeName("heatTransferGranConduction");
+    TypeName("heatTransferInterGrain");
 
     // Constructors
 
         //- Construct from components
-        heatTransferGranConduction
+        heatTransferInterGrain
         (
             const dictionary& dict,
             cfdemCloudEnergy& sm
@@ -105,7 +145,7 @@ public:
 
     // Destructor
 
-        virtual ~heatTransferGranConduction();
+        virtual ~heatTransferInterGrain();
 
 
     // Member Functions

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

@@ -47,6 +47,7 @@ reactionHeat::reactionHeat
     propsDict_(dict.subDict(typeName + "Props")),
     interpolation_(propsDict_.lookupOrDefault<bool>("interpolation",false)),
     verbose_(propsDict_.lookupOrDefault<bool>("verbose",false)),
+    execution_(true),
     mesh_(sm.mesh()),
     maxSource_(1e30),
     reactionHeatName_(propsDict_.lookupOrDefault<word>("reactionHeatName","reactionHeat")),
@@ -62,8 +63,14 @@ reactionHeat::reactionHeat
         ),
         mesh_,
         dimensionedScalar("zero", dimensionSet(1,-1,-3,0,0,0,0),0.0)
-    )
+    ),
+    Nevery_(propsDict_.lookupOrDefault<label>("Nevery",1)),
+    couplingTimestep_(0.0)
 {
+    scalar dtDEM = particleCloud_.dataExchangeM().DEMts();
+    scalar dtCFD = mesh_.time().deltaTValue();
+    couplingTimestep_ = max(dtDEM,dtCFD);
+
     particleCloud_.registerParticleProperty<double**>(reactionHeatName_,1);
 
     if(propsDict_.found("maxsource"))
@@ -86,6 +93,12 @@ reactionHeat::~reactionHeat()
 
 void reactionHeat::calcEnergyContribution()
 {
+    execution_ = (particleCloud_.dataExchangeM().couplingStep() % Nevery_ == 0);
+    if (!execution_)
+    {
+        return;
+    }
+
     double**& reactionHeat_ = particleCloud_.getParticlePropertyRef<double**>(reactionHeatName_);
 
     particleCloud_.dataExchangeM().getData(reactionHeatName_,"scalar-atom",reactionHeat_);
@@ -109,7 +122,7 @@ void reactionHeat::calcEnergyContribution()
         NULL
     );
 
-    reactionHeatField_.primitiveFieldRef() /= (reactionHeatField_.mesh().V());
+    reactionHeatField_.primitiveFieldRef() /= (reactionHeatField_.mesh().V() * Nevery_ * couplingTimestep_);
 
     forAll(reactionHeatField_,cellI)
     {
@@ -121,6 +134,12 @@ void reactionHeat::calcEnergyContribution()
         }
     }
 
+    if (verbose_)
+    {
+        Info << "reaction heat per unit time = "
+                     << gSum(reactionHeatField_*1.0*reactionHeatField_.mesh().V()) << endl;
+    }
+
     reactionHeatField_.correctBoundaryConditions();
 }
 

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

@@ -50,6 +50,8 @@ protected:
 
     bool verbose_;
 
+    bool execution_;
+
     const fvMesh& mesh_;
 
     scalar maxSource_;
@@ -58,6 +60,10 @@ protected:
 
     volScalarField reactionHeatField_;
 
+    label Nevery_;
+
+    scalar couplingTimestep_;
+
 public:
 
     //- Runtime type information

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

@@ -53,14 +53,30 @@ particleDeformation::particleDeformation
     initialExec_(true),
     refFieldName_(propsDict_.lookup("refFieldName")),
     refField_(),
+    voidfractionFieldName_(propsDict_.lookupOrDefault<word>("voidfractionFieldName","voidfraction")),
+#if OPENFOAM_VERSION_MAJOR < 5
+    voidfraction_(const_cast<volScalarField&>(sm.mesh().lookupObject<volScalarField>(voidfractionFieldName_))),
+#else
+    voidfraction_(sm.mesh().lookupObjectRef<volScalarField> (voidfractionFieldName_)),
+#endif
     defaultDeformCellsName_(propsDict_.lookupOrDefault<word>("defaultDeformCellsName","none")),
     defaultDeformCells_(),
     existDefaultDeformCells_(false),
     defaultDeformation_(propsDict_.lookupOrDefault<scalar>("defaultDeformation",1.0)),
+    existBackgroundVoidage_(false),
+    backgroundVoidage_(propsDict_.lookupOrDefault<scalar>("backgroundVoidage",-1.0)),
+    backgroundRef_(propsDict_.lookupOrDefault<scalar>("backgroundRef",0.35)),
     partTypes_(propsDict_.lookupOrDefault<labelList>("partTypes",labelList(1,-1))),
     lowerBounds_(propsDict_.lookupOrDefault<scalarList>("lowerBounds",scalarList(1,-1.0))),
     upperBounds_(propsDict_.lookupOrDefault<scalarList>("upperBounds",scalarList(1,-1.0))),
-    partDeformationsName_("partDeformations")
+    partDeformationsName_("partDeformations"),
+    controlForceOnDefPart_(propsDict_.lookupOrDefault<bool>("controlForceOnDefPart",false)),
+    controlFieldName_(propsDict_.lookupOrDefault<word>("controlFieldName","")),
+    controlField_(NULL),
+    controlPoint_(propsDict_.lookupOrDefault<vector>("controlPoint",vector::zero)),
+    controlCell_(-1),
+    controlTargetValue_(propsDict_.lookupOrDefault<scalar>("controlTargetValue",0.0)),
+    controlCouplingStrength_(propsDict_.lookupOrDefault<scalar>("controlCouplingStrength",0.0))
 {
     particleCloud_.registerParticleProperty<double**>(partDeformationsName_,1);
 
@@ -89,6 +105,11 @@ particleDeformation::particleDeformation
        }
     }
 
+    if (backgroundVoidage_ <= 1.0 && backgroundVoidage_ >= 0.01)
+    {
+        existBackgroundVoidage_ = true;
+    }
+
     // check if only single value instead of list was provided
     if (propsDict_.found("partType"))
     {
@@ -113,6 +134,11 @@ particleDeformation::particleDeformation
     Info << "partTypes: " << partTypes_ << endl;
     Info << "lowerBounds: " << lowerBounds_ << endl;
     Info << "upperBounds: " << upperBounds_ << endl;
+
+    if (controlForceOnDefPart_)
+    {
+        controlCell_ = sm.mesh().findCell(controlPoint_);
+    }
 }
 
 
@@ -139,6 +165,16 @@ void particleDeformation::setForce() const
         initialExec_ = false;
     }
 
+    scalar controlCurrValue = 0.0;
+    scalar scaleFactor = 0.0;
+    if (controlForceOnDefPart_)
+    {
+        controlCurrValue = (*controlField_)[controlCell_];
+        reduce(controlCurrValue, sumOp<scalar>());
+        scaleFactor = 1.0 + controlCouplingStrength_ * (controlTargetValue_ - controlCurrValue);
+        Info << "particleDeformation: scaleFactor = " << scaleFactor << endl;
+    }
+
     double**& partDeformations_ = particleCloud_.getParticlePropertyRef<double**>(partDeformationsName_);
 
     label cellI = 0;
@@ -187,6 +223,15 @@ void particleDeformation::setForce() const
 
             partDeformations_[index][0] = deformationDegree;
 
+            if (existBackgroundVoidage_ && deformationDegree > 0.1)
+            {
+                voidfraction_[cellI] = backgroundRef_ + deformationDegree * (backgroundVoidage_ - backgroundRef_);
+            }
+
+            if (controlForceOnDefPart_)
+            {
+                particleCloud_.scaleForce(index,scaleFactor);
+            }
 
             if(forceSubM(0).verbose() && index >= 0 && index < 2)
             {
@@ -222,12 +267,18 @@ void particleDeformation::init() const
                     particleCloud_.mesh().time().timeName(tstart),
                     particleCloud_.mesh(),
                     IOobject::MUST_READ,
-                    IOobject::NO_WRITE
+                    IOobject::AUTO_WRITE
                 ),
                 particleCloud_.mesh()
             )
         );
     }
+
+    if (controlForceOnDefPart_)
+    {
+        volScalarField& controlField(const_cast<volScalarField&>(particleCloud_.mesh().lookupObject<volScalarField> (controlFieldName_)));
+        controlField_ = &controlField;
+    }
 }
 
 label particleDeformation::getListIndex(label testElement) const

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

@@ -59,6 +59,10 @@ private:
 
     mutable autoPtr<volScalarField> refField_;
 
+    word voidfractionFieldName_;
+
+    volScalarField& voidfraction_;
+
     // default deformation in region
     const word defaultDeformCellsName_;
 
@@ -68,6 +72,12 @@ private:
 
     scalar defaultDeformation_;
 
+    bool existBackgroundVoidage_;
+
+    scalar backgroundVoidage_;
+
+    scalar backgroundRef_;
+
     labelList partTypes_;
 
     scalarList lowerBounds_;
@@ -76,6 +86,22 @@ private:
 
     const word partDeformationsName_;
 
+    // control functionality
+
+    bool controlForceOnDefPart_;
+
+    word controlFieldName_;
+
+    mutable volScalarField* controlField_;
+
+    vector controlPoint_;
+
+    label controlCell_;
+
+    scalar controlTargetValue_;
+
+    scalar controlCouplingStrength_;
+
     label getListIndex(label) const;
 
     void init() const;

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

@@ -0,0 +1,263 @@
+/*---------------------------------------------------------------------------*\
+    CFDEMcoupling academic - Open Source CFD-DEM coupling
+
+    Contributing authors:
+    Thomas Lichtenegger, Sanaz Abbasi
+    Copyright (C) 2015- Johannes Kepler University, Linz
+-------------------------------------------------------------------------------
+License
+    This file is part of CFDEMcoupling academic.
+
+    CFDEMcoupling academic 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 academic 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 academic.  If not, see <http://www.gnu.org/licenses/>.
+\*---------------------------------------------------------------------------*/
+
+#include "error.H"
+
+#include "terminalVelocity.H"
+#include "addToRunTimeSelectionTable.H"
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(terminalVelocity, 0);
+
+addToRunTimeSelectionTable
+(
+    forceModel,
+    terminalVelocity,
+    dictionary
+);
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from components
+terminalVelocity::terminalVelocity
+(
+    const dictionary& dict,
+    cfdemCloud& sm
+)
+:
+    forceModel(dict,sm),
+    propsDict_(dict.subDict(typeName + "Props")),
+    interpolate_(propsDict_.lookupOrDefault<bool>("interpolation", false)),
+    mesh_(sm.mesh()),
+    ignoreCellsName_(propsDict_.lookupOrDefault<word>("ignoreCellsName","none")),
+    ignoreCells_(),
+    existIgnoreCells_(true),
+    existturbDissipationRateInObjReg_(false),
+    turbulenceCorrection_(propsDict_.lookupOrDefault<bool>("turbulenceCorrection",false)),
+    turbDissipationRate_(NULL),
+    wallIndicatorField_
+    (
+        IOobject
+        (
+            "wallIndicator",
+            sm.mesh().time().timeName(),
+            sm.mesh(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        sm.mesh(),
+        dimensionedScalar("zero", dimensionSet(0,0,0,0,0,0,0), 0.0)
+    ),
+    liquidViscosity_(propsDict_.lookupOrDefault<scalar>("liquidViscosity", 1.0)),
+    dragReductionFactor_(propsDict_.lookupOrDefault<scalar>("dragReductionFactor", 0.0)),
+    gravityFieldName_(propsDict_.lookupOrDefault<word>("gravityFieldName","g")),
+    g_(sm.mesh().lookupObject<uniformDimensionedVectorField> (gravityFieldName_))
+{
+    // init force sub model
+    setForceSubModels(propsDict_);
+    // define switches which can be read from dict
+    forceSubM(0).setSwitchesList(SW_VERBOSE,true); // activate search for verbose switch
+// TODO: remove bool interpolate for this class, let forceSubModel do this
+//    forceSubM(0).setSwitchesList(SW_INTERPOLATION,true); // activate search for interpolate switch
+    forceSubM(0).readSwitches();
+
+    scalar terminalVelMagnitude(propsDict_.lookupOrDefault<scalar>("terminalVelocity", 0.0));
+    terminalVel_ = -terminalVelMagnitude * g_.value() / mag(g_.value());
+
+    if (ignoreCellsName_ != "none")
+    {
+       ignoreCells_.set(new cellSet(particleCloud_.mesh(),ignoreCellsName_));
+       Info<< "terminalVelocity: ignoring rising velocity in cellSet " << ignoreCells_().name()
+           << " with " << ignoreCells_().size() << " cells." << endl;
+    }
+    else
+    {
+        existIgnoreCells_ = false;
+    }
+
+    turbDissipationRate_ = new volScalarField
+    (
+        IOobject
+        (
+            "turbDissipationRate",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh_,
+        dimensionedScalar("zero", dimensionSet(0,2,-3,0,0), 0)
+    );
+
+    // define a field to indicate if a cell is next to boundary
+    label cellI = -1;
+    forAll(mesh_.boundary(), patchI)
+    {
+        word patchName = mesh_.boundary()[patchI].name();
+        if (patchName.rfind("procB",0) == 0) continue;
+
+        forAll(mesh_.boundary()[patchI], faceI)
+        {
+            cellI = mesh_.boundary()[patchI].faceCells()[faceI];
+            wallIndicatorField_[cellI] = 1.0;
+        }
+    }
+
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+terminalVelocity::~terminalVelocity()
+{
+    if (!existturbDissipationRateInObjReg_) delete turbDissipationRate_;
+}
+
+// * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
+
+bool terminalVelocity::ignoreCell(label cell) const
+{
+    if (!existIgnoreCells_) return false;
+    else return ignoreCells_()[cell];
+}
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void terminalVelocity::setForce() const
+{
+    updateEpsilon();
+
+    vector position(0,0,0);
+    label cellI = -1;
+    scalar radius = 0.0;
+    scalar epsilon = 0.0;
+    scalar dLambda = 0.0;
+    scalar velReductionFactor = 0.0;
+    vector Uparticle(0,0,0);
+
+    label patchID = -1;
+    label faceIGlobal = -1;
+    scalar velProjection = 0.0;
+    vector faceINormal = vector::zero;
+    word patchName("");
+
+    interpolationCellPoint<scalar> turbDissipationRateInterpolator_(*turbDissipationRate_);
+
+    for (int index = 0; index < particleCloud_.numberOfParticles(); ++index)
+    {
+        cellI = particleCloud_.cellIDs()[index][0];
+        Uparticle = vector::zero;
+
+        if (cellI > -1 && !ignoreCell(cellI)) // particle found
+        {
+            if (interpolate_)
+            {
+                position = particleCloud_.position(index);
+                epsilon = turbDissipationRateInterpolator_.interpolate(position,cellI);
+            }
+            else
+            {
+                epsilon = (*turbDissipationRate_)[cellI];
+            }
+
+            position = particleCloud_.position(index);
+            radius = particleCloud_.radius(index);
+
+            if (turbulenceCorrection_)
+            {
+                // d * kolmogorov length scale
+                dLambda = 2*radius*pow(epsilon,0.25)/pow(liquidViscosity_,0.75);
+                velReductionFactor = Foam::sqrt(1 + (dragReductionFactor_*pow(dLambda,3)));
+                terminalVel_ =  terminalVel_ / velReductionFactor;
+            }
+
+            // read the new particle velocity
+            for (int j = 0; j < 3; j++)
+            {
+                particleCloud_.particleConvVels()[index][j] += terminalVel_[j];
+                Uparticle[j] = particleCloud_.particleConvVels()[index][j];
+            }
+
+            // prevent particles being pushed through walls
+            // check if cell is adjacent to wall and remove the normal velocity to the wall
+            if (wallIndicatorField_[cellI] > 0.5)
+            {
+                const cell& faces = mesh_.cells()[cellI];
+                forAll(faces, faceI)
+                {
+                    faceIGlobal = faces[faceI];
+                    patchID = mesh_.boundaryMesh().whichPatch(faceIGlobal);
+                    if (patchID < 0) continue;
+                    patchName = mesh_.boundary()[patchID].name();
+
+                    if (patchName.rfind("procB",0) == 0) continue;
+
+                    faceINormal = mesh_.Sf()[faceIGlobal];
+                    faceINormal /= mag(faceINormal);
+                    velProjection = faceINormal&Uparticle;
+                    if (velProjection > 0.0)
+                    {
+                        // removes the value normal to the face
+                        for (int j = 0; j < 3; j++)
+                        {
+                            particleCloud_.particleConvVels()[index][j] -= velProjection*faceINormal[j];
+                        }
+                    }
+                }
+            }
+        }
+
+        if (forceSubM(0).verbose() && index > 0 && index < 2)
+        {
+            Pout<< "cellI = " << cellI << endl;
+            Pout<< "index = " << index << endl;
+            Pout<< "epsilon = " << epsilon << endl;
+            Pout<< "rising velocity = " << terminalVel_ << endl;
+        }
+    }
+}
+
+
+void terminalVelocity::updateEpsilon() const
+{
+    if (!existturbDissipationRateInObjReg_)
+    {
+        Info<< "epsilon is calculated from the turbulence model. " << endl;
+        *turbDissipationRate_ = particleCloud_.turbulence().epsilon()();
+    }
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

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

@@ -0,0 +1,130 @@
+/*---------------------------------------------------------------------------*\
+    CFDEMcoupling academic - Open Source CFD-DEM coupling
+
+    Contributing authors:
+    Thomas Lichtenegger, Sanaz Abbasi
+    Copyright (C) 2015- Johannes Kepler University, Linz
+-------------------------------------------------------------------------------
+License
+    This file is part of CFDEMcoupling academic.
+
+    CFDEMcoupling academic 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 academic 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 academic.  If not, see <http://www.gnu.org/licenses/>.
+\*---------------------------------------------------------------------------*/
+
+#ifndef terminalVelocity_H
+#define terminalVelocity_H
+
+#include "forceModel.H"
+#include "interpolationCellPoint.H"
+#include "cellSet.H"
+//#include "recBase.H"
+//#include "recModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class freeStreaming Declaration
+\*---------------------------------------------------------------------------*/
+
+class terminalVelocity
+:
+    public forceModel
+{
+private:
+    dictionary propsDict_;
+
+    bool interpolate_;
+
+    const fvMesh& mesh_;
+
+    // ignore streaming velocity in region
+    word ignoreCellsName_;
+
+    autoPtr<cellSet> ignoreCells_;
+
+    bool existIgnoreCells_;
+
+//    word turbDissipationRateFieldName_;
+
+//    word turbKineticEnergyFieldName_;
+
+//    word recurrenceBaseName_;
+
+//    volScalarField* turbKineticEnergy_;
+
+    bool existturbDissipationRateInObjReg_;
+
+    bool turbulenceCorrection_;
+
+    volScalarField* turbDissipationRate_;
+
+//    mutable volScalarField delta;
+
+    volScalarField wallIndicatorField_;
+
+    scalar liquidViscosity_;
+
+    // drag reduction factor in turbulent flow
+    scalar dragReductionFactor_;
+
+    // terminal rising velocity of the bubble in the resting liquid
+    mutable vector terminalVel_;
+
+    word gravityFieldName_;
+
+    const uniformDimensionedVectorField& g_;    // ref to gravity
+
+    bool ignoreCell(label) const;
+
+    void updateEpsilon() const;
+
+//    recBase* recurrenceBase_;
+
+public:
+
+    //- Runtime type information
+    TypeName("terminalVelocity");
+
+
+    // Constructors
+
+        //- Construct from components
+        terminalVelocity
+        (
+            const dictionary& dict,
+            cfdemCloud& sm
+        );
+
+    // Destructor
+
+        ~terminalVelocity();
+
+
+    // Member Functions
+        void setForce() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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

@@ -62,6 +62,24 @@ turbulentDispersion::turbulentDispersion
     ignoreCellsName_(propsDict_.lookupOrDefault<word>("ignoreCellsName","none")),
     ignoreCells_(),
     existIgnoreCells_(true),
+    usePreCalcNut_(propsDict_.lookupOrDefault<bool>("usePreCalcNut",false)),
+    nutName_(propsDict_.lookupOrDefault<word>("nutName","nut")),
+    usePreCalcK_(propsDict_.lookupOrDefault<bool>("usePreCalcK",false)),
+    kName_(propsDict_.lookupOrDefault<word>("kName","k")),
+    usePreCalcDispField_(propsDict_.lookupOrDefault<bool>("usePreCalcDispField",false)),
+    dispFieldName_(propsDict_.lookupOrDefault<word>("dispFieldName","dispVarField")),
+    nut_
+    (   IOobject
+        (
+            "nuTurb",
+            sm.mesh().time().timeName(),
+            sm.mesh(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        sm.mesh(),
+        dimensionedScalar("zero", dimensionSet(0,2,-1,0,0,0,0), 0.0)
+    ),
     wallIndicatorField_
     (   IOobject
         (
@@ -74,13 +92,23 @@ turbulentDispersion::turbulentDispersion
         sm.mesh(),
         dimensionedScalar("zero", dimensionSet(0,0,0,0,0,0,0), 0.0)
     ),
-    minTurbKineticEnergy_(propsDict_.lookupOrDefault<scalar>("minTurbKineticEnergy", 0.0)),
-    turbKineticEnergyFieldName_(propsDict_.lookupOrDefault<word>("turbKineticEnergyFieldName","")),
-    turbKineticEnergy_(NULL),
-    existTurbKineticEnergyInObjReg_(false),
+    delta_
+    (   IOobject
+        (
+            "delta",
+            sm.mesh().time().timeName(),
+            sm.mesh(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        sm.mesh(),
+        dimensionedScalar("delta", dimLength, 1.0)
+    ),
+    turbulentSchmidtNumber_(readScalar(propsDict_.lookup("turbulentSchmidtNumber"))),
     voidfractionFieldName_(propsDict_.lookupOrDefault<word>("voidfractionFieldName","voidfraction")),
     voidfraction_(sm.mesh().lookupObject<volScalarField> (voidfractionFieldName_)),
     critVoidfraction_(propsDict_.lookupOrDefault<scalar>("critVoidfraction", 0.9)),
+    Ck_(propsDict_.lookupOrDefault<scalar>("turbKineticEnergyCoeff", 0.094)),
     ranGen_(clock::getTime()+pid())
 {
     if (ignoreCellsName_ != "none")
@@ -91,6 +119,16 @@ turbulentDispersion::turbulentDispersion
     }
     else existIgnoreCells_ = false;
 
+    if ((usePreCalcNut_ && usePreCalcK_) || (usePreCalcNut_ && usePreCalcDispField_) || (usePreCalcDispField_ && usePreCalcK_))
+    {
+        FatalError<< "Cannot use more than one precalculated nut, k and displacement fluctuations at the same time. Choose one." << abort(FatalError);
+    }
+
+    if (usePreCalcK_)
+    {
+        delta_.primitiveFieldRef()=Foam::pow(mesh_.V(),1.0/3.0);
+    }
+
     // define a field to indicate if a cell is next to boundary
      label cellI = -1;
      forAll (mesh_.boundary(),patchI)
@@ -105,45 +143,18 @@ turbulentDispersion::turbulentDispersion
          }
      }
 
-    if (turbKineticEnergyFieldName_ != "")
-    {
-        existTurbKineticEnergyInObjReg_ = true;
-        volScalarField& k(const_cast<volScalarField&>(sm.mesh().lookupObject<volScalarField> (turbKineticEnergyFieldName_)));
-        turbKineticEnergy_ = &k;
-    }
-    else
-    {
-        turbKineticEnergy_ = new volScalarField
-        (
-            IOobject
-            (
-                "turbKinEnergy",
-                mesh_.time().timeName(),
-                mesh_,
-                IOobject::NO_READ,
-                IOobject::NO_WRITE
-            ),
-            mesh_,
-            dimensionedScalar("zero", dimensionSet(0,2,-2,0,0), 0)
-        );
-    }
-
-    // make sure this is the last force model in list so that fluid velocity does not get overwritten
-    label numLastForceModel = sm.nrForceModels();
-    word lastForceModel = sm.forceModels()[numLastForceModel-1];
-    if (lastForceModel != "turbulentDispersion")
-    {
-        FatalError <<"Force model 'turbulentDispersion' needs to be last in list!\n" << abort(FatalError);
-    }
+    scalar dtCFD = voidfraction_.mesh().time().deltaTValue();
+    scalar dtDEM = particleCloud_.dataExchangeM().DEMts();
+    // if CFD step is larger than DEM step, a corresponding number of DEM steps is taken to reach the CFD step;
+    // if DEM step is larger than CFD step, no update occurs so that full DEM step needs to be taken
+    dt_ = max(dtCFD, dtDEM);
 }
 
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
 turbulentDispersion::~turbulentDispersion()
-{
-    if (!existTurbKineticEnergyInObjReg_) delete turbKineticEnergy_;
-}
+{}
 
 // * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
 
@@ -157,45 +168,72 @@ bool turbulentDispersion::ignoreCell(label cell) const
 
 void turbulentDispersion::setForce() const
 {
-    if (!existTurbKineticEnergyInObjReg_)
+    if (!usePreCalcNut_ && !usePreCalcK_ && !usePreCalcDispField_)
     {
-        *turbKineticEnergy_ = particleCloud_.turbulence().k()();
+        nut_ = particleCloud_.turbulence().nut()();
+    }
+    else if (usePreCalcNut_)
+    {
+        volScalarField& nutRef (const_cast<volScalarField&>(mesh_.lookupObject<volScalarField> (nutName_)));
+        nut_ = nutRef;
+    }
+    else if (usePreCalcK_)
+    {
+        volScalarField& kRef (const_cast<volScalarField&>(mesh_.lookupObject<volScalarField> (kName_)));
+        nut_ = Ck_ * delta_ * sqrt(kRef);
+    }
+    else
+    {
+        volVectorField& dispFieldRef (const_cast<volVectorField&>(mesh_.lookupObject<volVectorField> (dispFieldName_)));
+        dispField_ = &dispFieldRef;
     }
 
     label cellI = -1;
     label patchID = -1;
     label faceIGlobal = -1;
     scalar flucProjection = 0.0;
-    scalar k = 0.0;
+    scalar D = 0.0;
+    scalar randScalar = 0.0;
     vector faceINormal = vector::zero;
     vector flucU = vector::zero;
     vector position = vector::zero;
     word patchName("");
 
-    interpolationCellPoint<scalar> turbKineticEnergyInterpolator_(*turbKineticEnergy_);
+    interpolationCellPoint<scalar> nutInterpolator_(nut_);
 
-    for(int index = 0;index <  particleCloud_.numberOfParticles(); ++index)
+    for(int index = 0; index < particleCloud_.numberOfParticles(); ++index)
     {
         cellI = particleCloud_.cellIDs()[index][0];
         if (cellI > -1 && !ignoreCell(cellI))
         {
-            // particles in dilute regions follow fluid without fluctuations
-
-            if (voidfraction_[cellI] < critVoidfraction_)
+            if (usePreCalcDispField_)
+            {
+                for (label comp=0; comp<3; comp++)
+                {
+#if OPENFOAM_VERSION_MAJOR < 6
+                    randScalar = ranGen_.GaussNormal();
+#else
+                    randScalar = ranGen_.scalarNormal();
+#endif
+                    flucU.component(comp) = randScalar * (*dispField_)[cellI].component(comp);
+                }
+            }
+            else
             {
                 if (interpolate_)
                 {
                     position = particleCloud_.position(index);
-                    k = turbKineticEnergyInterpolator_.interpolate(position,cellI);
+                    D = nutInterpolator_.interpolate(position,cellI) / turbulentSchmidtNumber_;
                 }
                 else
                 {
-                    k = (*turbKineticEnergy_)[cellI];
+                    D = nut_[cellI] / turbulentSchmidtNumber_;
                 }
 
-                if (k < minTurbKineticEnergy_) k = minTurbKineticEnergy_;
+                // include concentration dependence on the diffusivity at this point if necessary
 
-                flucU=unitFlucDir()*Foam::sqrt(2.0*k);
+                flucU=unitFlucDir()*Foam::sqrt(6.0*D/dt_);
+            }
 
             // prevent particles being pushed through walls by regulating velocity fluctuations
             // check if cell is adjacent to wall and remove corresponding components
@@ -220,8 +258,7 @@ void turbulentDispersion::setForce() const
 
             for(int j=0;j<3;j++)
             {
-                    particleCloud_.fluidVels()[index][j] += flucU[j];
-                }
+                particleCloud_.particleFlucVels()[index][j] += flucU[j];
             }
         }
     }

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

@@ -58,15 +58,29 @@ protected:
 
     bool existIgnoreCells_;
 
-    mutable volScalarField wallIndicatorField_;
+    bool usePreCalcNut_;
 
-    scalar minTurbKineticEnergy_;
+    word nutName_;
 
-    word turbKineticEnergyFieldName_;
+    bool usePreCalcK_;
 
-    volScalarField* turbKineticEnergy_;
+    word kName_;
 
-    bool existTurbKineticEnergyInObjReg_;
+    bool usePreCalcDispField_;
+
+    word dispFieldName_;
+
+    mutable volVectorField* dispField_;
+
+    mutable volScalarField nut_;
+
+    volScalarField wallIndicatorField_;
+
+    volScalarField delta_;
+
+    scalar dt_;
+
+    scalar turbulentSchmidtNumber_;
 
     word voidfractionFieldName_;
 
@@ -74,10 +88,12 @@ protected:
 
     scalar critVoidfraction_;
 
-    virtual vector unitFlucDir() const;
+    scalar Ck_;
 
     mutable Random ranGen_;
 
+    virtual vector unitFlucDir() const;
+
     bool ignoreCell(label) const;
 
 public:

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

@@ -0,0 +1,256 @@
+/*---------------------------------------------------------------------------*\
+    CFDEMcoupling academic - Open Source CFD-DEM coupling
+
+    Contributing authors:
+    Thomas Lichtenegger
+    Copyright (C) 2015- Johannes Kepler University, Linz
+-------------------------------------------------------------------------------
+License
+    This file is part of CFDEMcoupling academic.
+
+    CFDEMcoupling academic 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 academic 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 academic.  If not, see <http://www.gnu.org/licenses/>.
+\*---------------------------------------------------------------------------*/
+
+#include "error.H"
+
+#include "turbulentVelocityFluctuations.H"
+#include "addToRunTimeSelectionTable.H"
+#include "OFstream.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(turbulentVelocityFluctuations, 0);
+
+addToRunTimeSelectionTable
+(
+    forceModel,
+    turbulentVelocityFluctuations,
+    dictionary
+);
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from components
+turbulentVelocityFluctuations::turbulentVelocityFluctuations
+(
+    const dictionary& dict,
+    cfdemCloud& sm,
+    word type
+)
+:
+    forceModel(dict,sm),
+    propsDict_(dict.subDict(type + "Props")),
+    interpolate_(propsDict_.lookupOrDefault<bool>("interpolation", false)),
+    mesh_(sm.mesh()),
+    ignoreCellsName_(propsDict_.lookupOrDefault<word>("ignoreCellsName","none")),
+    ignoreCells_(),
+    existIgnoreCells_(true),
+    wallIndicatorField_
+    (   IOobject
+        (
+            "wallIndicator",
+            sm.mesh().time().timeName(),
+            sm.mesh(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        sm.mesh(),
+        dimensionedScalar("zero", dimensionSet(0,0,0,0,0,0,0), 0.0)
+    ),
+    minTurbKineticEnergy_(propsDict_.lookupOrDefault<scalar>("minTurbKineticEnergy", 0.0)),
+    turbKineticEnergyFieldName_(propsDict_.lookupOrDefault<word>("turbKineticEnergyFieldName","")),
+    turbKineticEnergy_(NULL),
+    existTurbKineticEnergyInObjReg_(false),
+    voidfractionFieldName_(propsDict_.lookupOrDefault<word>("voidfractionFieldName","voidfraction")),
+    voidfraction_(sm.mesh().lookupObject<volScalarField> (voidfractionFieldName_)),
+    critVoidfraction_(propsDict_.lookupOrDefault<scalar>("critVoidfraction", 0.9)),
+    ranGen_(clock::getTime()+pid())
+{
+    if (ignoreCellsName_ != "none")
+    {
+       ignoreCells_.set(new cellSet(particleCloud_.mesh(),ignoreCellsName_));
+       Info << type << ": ignoring fluctuations in cellSet " << ignoreCells_().name() <<
+        " with " << ignoreCells_().size() << " cells." << endl;
+    }
+    else existIgnoreCells_ = false;
+
+    // define a field to indicate if a cell is next to boundary
+     label cellI = -1;
+     forAll (mesh_.boundary(),patchI)
+     {
+         word patchName = mesh_.boundary()[patchI].name();
+         if (patchName.rfind("procB",0) == 0) continue;
+
+         forAll(mesh_.boundary()[patchI], faceI)
+         {
+             cellI = mesh_.boundary()[patchI].faceCells()[faceI];
+             wallIndicatorField_[cellI] = 1.0;
+         }
+     }
+
+    if (turbKineticEnergyFieldName_ != "")
+    {
+        existTurbKineticEnergyInObjReg_ = true;
+        volScalarField& k(const_cast<volScalarField&>(sm.mesh().lookupObject<volScalarField> (turbKineticEnergyFieldName_)));
+        turbKineticEnergy_ = &k;
+    }
+    else
+    {
+        turbKineticEnergy_ = new volScalarField
+        (
+            IOobject
+            (
+                "turbKinEnergy",
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh_,
+            dimensionedScalar("zero", dimensionSet(0,2,-2,0,0), 0)
+        );
+    }
+
+    // make sure this is the last force model in list so that fluid velocity does not get overwritten
+    label numLastForceModel = sm.nrForceModels();
+    word lastForceModel = sm.forceModels()[numLastForceModel-1];
+    if (lastForceModel != "turbulentVelocityFluctuations")
+    {
+        FatalError <<"Force model 'turbulentVelocityFluctuations' needs to be last in list!\n" << abort(FatalError);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+turbulentVelocityFluctuations::~turbulentVelocityFluctuations()
+{
+    if (!existTurbKineticEnergyInObjReg_) delete turbKineticEnergy_;
+}
+
+// * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
+
+bool turbulentVelocityFluctuations::ignoreCell(label cell) const
+{
+    if (!existIgnoreCells_) return false;
+    else return ignoreCells_()[cell];
+}
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void turbulentVelocityFluctuations::setForce() const
+{
+    if (!existTurbKineticEnergyInObjReg_)
+    {
+        *turbKineticEnergy_ = particleCloud_.turbulence().k()();
+    }
+
+    label cellI = -1;
+    label patchID = -1;
+    label faceIGlobal = -1;
+    scalar flucProjection = 0.0;
+    scalar k = 0.0;
+    vector faceINormal = vector::zero;
+    vector flucU = vector::zero;
+    vector position = vector::zero;
+    word patchName("");
+
+    interpolationCellPoint<scalar> turbKineticEnergyInterpolator_(*turbKineticEnergy_);
+
+    for(int index = 0;index <  particleCloud_.numberOfParticles(); ++index)
+    {
+        cellI = particleCloud_.cellIDs()[index][0];
+        if (cellI > -1 && !ignoreCell(cellI))
+        {
+            // particles in dilute regions follow fluid without fluctuations
+
+            if (voidfraction_[cellI] < critVoidfraction_)
+            {
+                if (interpolate_)
+                {
+                    position = particleCloud_.position(index);
+                    k = turbKineticEnergyInterpolator_.interpolate(position,cellI);
+                }
+                else
+                {
+                    k = (*turbKineticEnergy_)[cellI];
+                }
+
+                if (k < minTurbKineticEnergy_) k = minTurbKineticEnergy_;
+
+                flucU=unitFlucDir()*Foam::sqrt(2.0*k);
+
+                // prevent particles being pushed through walls by regulating velocity fluctuations
+                // check if cell is adjacent to wall and remove corresponding components
+                if (wallIndicatorField_[cellI] > 0.5)
+                {
+                    const cell& faces = mesh_.cells()[cellI];
+                    forAll (faces, faceI)
+                    {
+                        faceIGlobal = faces[faceI];
+                        patchID = mesh_.boundaryMesh().whichPatch(faceIGlobal);
+                        if (patchID < 0) continue;
+                        patchName = mesh_.boundary()[patchID].name();
+
+                        if (patchName.rfind("procB",0) == 0) continue;
+
+                        faceINormal = mesh_.Sf()[faceIGlobal];
+                        faceINormal /= mag(faceINormal);
+                        flucProjection = faceINormal&flucU;
+                        if (flucProjection > 0.0) flucU -= flucProjection*faceINormal;
+                    }
+                }
+
+                for(int j=0;j<3;j++)
+                {
+                    particleCloud_.fluidVels()[index][j] += flucU[j];
+                }
+            }
+        }
+    }
+}
+
+vector turbulentVelocityFluctuations::unitFlucDir() const
+{
+    // unit random vector
+    // algorithm according to:
+    // Marsaglia. "Choosing a point from the surface of a sphere." The Annals of Mathematical Statistics 43.2 (1972): 645-646.
+    scalar v1(0.0);
+    scalar v2(0.0);
+    scalar s(10.0);
+    scalar s2(0.0);
+    vector rvec(0,0,0);
+    while(s>1.0)
+    {
+        v1=2*(ranGen_.scalar01()-0.5);
+        v2=2*(ranGen_.scalar01()-0.5);
+        s=v1*v1+v2*v2;
+    }
+    s2=Foam::sqrt(1-s);
+    rvec[0]=2*v1*s2;
+    rvec[1]=2*v2*s2;
+    rvec[2]=1-2*s;
+    return rvec;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

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

@@ -0,0 +1,117 @@
+/*---------------------------------------------------------------------------*\
+    CFDEMcoupling academic - Open Source CFD-DEM coupling
+
+    Contributing authors:
+    Thomas Lichtenegger
+    Copyright (C) 2015- Johannes Kepler University, Linz
+-------------------------------------------------------------------------------
+License
+    This file is part of CFDEMcoupling academic.
+
+    CFDEMcoupling academic 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 academic 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 academic.  If not, see <http://www.gnu.org/licenses/>.
+\*---------------------------------------------------------------------------*/
+
+#ifndef turbulentVelocityFluctuations_H
+#define turbulentVelocityFluctuations_H
+
+#include "forceModel.H"
+#include "interpolationCellPoint.H"
+#include "Random.H"
+#include "autoPtr.H"
+#include "cellSet.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class turbulentVelocityFluctuations Declaration
+\*---------------------------------------------------------------------------*/
+
+class turbulentVelocityFluctuations
+:
+    public forceModel
+{
+protected:
+    dictionary propsDict_;
+
+    bool interpolate_;
+
+    const fvMesh& mesh_;
+
+    // ignore fluctuations in region
+    word ignoreCellsName_;
+
+    autoPtr<cellSet> ignoreCells_;
+
+    bool existIgnoreCells_;
+
+    volScalarField wallIndicatorField_;
+
+    scalar minTurbKineticEnergy_;
+
+    word turbKineticEnergyFieldName_;
+
+    volScalarField* turbKineticEnergy_;
+
+    bool existTurbKineticEnergyInObjReg_;
+
+    word voidfractionFieldName_;
+
+    const volScalarField& voidfraction_;
+
+    scalar critVoidfraction_;
+
+    virtual vector unitFlucDir() const;
+
+    mutable Random ranGen_;
+
+    bool ignoreCell(label) const;
+
+public:
+
+    //- Runtime type information
+    TypeName("turbulentVelocityFluctuations");
+
+
+    // Constructors
+
+        //- Construct from components
+        turbulentVelocityFluctuations
+        (
+            const dictionary& dict,
+            cfdemCloud& sm,
+            word type = "turbulentVelocityFluctuations"
+        );
+
+    // Destructor
+
+        ~turbulentVelocityFluctuations();
+
+
+    // Member Functions
+        void setForce() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/cfdemParticle/subModels/momCoupleModel/deactivateCouple/deactivateCouple.C

@@ -54,7 +54,34 @@ deactivateCouple::deactivateCouple
     cfdemCloud& sm
 )
 :
-    momCoupleModel(dict,sm)
+    momCoupleModel(dict,sm),
+    fEmpty_
+    (   IOobject
+        (
+            "fEmpty",
+            sm.mesh().time().timeName(),
+            sm.mesh(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        sm.mesh(),
+        dimensionedVector("zero", dimensionSet(1,-2,-2,0,0), vector::zero),
+        "zeroGradient"
+    ),
+    KslEmpty_
+    (
+        IOobject
+        (
+            "KslEmpty",
+            sm.mesh().time().timeName(),
+            sm.mesh(),
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        sm.mesh(),
+        dimensionedScalar("zero", dimensionSet(1,-3,-1,0,0), 0),
+        "zeroGradient"
+    )
 {}
 
 
@@ -65,9 +92,21 @@ deactivateCouple::~deactivateCouple()
 
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-void Foam::deactivateCouple::resetMomSourceField()
+void deactivateCouple::resetMomSourceField()
 {}
 
+tmp<volVectorField> deactivateCouple::expMomSource()
+{
+    tmp<volVectorField> tsource(fEmpty_);
+    return tsource;
+}
+
+tmp<volScalarField> deactivateCouple::impMomSource()
+{
+    tmp<volScalarField> tsource(KslEmpty_);
+    return tsource;
+
+}
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

src/lagrangian/cfdemParticle/subModels/momCoupleModel/deactivateCouple/deactivateCouple.H

@@ -42,6 +42,9 @@ class deactivateCouple
 {
 private:
 
+    volVectorField fEmpty_;
+
+    volScalarField KslEmpty_;
 
 public:
 
@@ -65,6 +68,10 @@ public:
 
     // Member Functions
         void resetMomSourceField();
+
+        tmp<volVectorField> expMomSource();
+
+        tmp<volScalarField> impMomSource();
 };
 
 

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

@@ -47,7 +47,8 @@ gravity::gravity
     voidfraction_(sm.mesh().lookupObject<volScalarField> (voidfractionFieldName_)),
     rhoGFieldName_(propsDict_.lookupOrDefault<word>("rhoGFieldName","rho")),
     rhoG_(sm.mesh().lookupObject<volScalarField> (rhoGFieldName_)),
-    g_("g",dimensionSet(0,1,-2,0,0),vector(0,0,-9.81))
+    gravityFieldName_(propsDict_.lookupOrDefault<word>("gravityFieldName","g")),
+    g_(sm.mesh().lookupObject<uniformDimensionedVectorField> (gravityFieldName_))
 {}
 
 

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

@@ -54,7 +54,9 @@ protected:
 
         const volScalarField& rhoG_;
 
-        const dimensionedVector g_;
+        word gravityFieldName_;
+
+        const uniformDimensionedVectorField& g_;    // ref to gravity
 
 public:
 

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

@@ -0,0 +1,96 @@
+/*---------------------------------------------------------------------------*\
+License
+    This is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+    This code is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+    You should have received a copy of the GNU General Public License
+    along with this code.  If not, see <http://www.gnu.org/licenses/>.
+
+    Copyright (C) 2015- Thomas Lichtenegger, JKU Linz, Austria
+
+\*---------------------------------------------------------------------------*/
+
+#include "error.H"
+#include "secondaryPhaseInducedBuoyancy.H"
+#include "mathExtra.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(secondaryPhaseInducedBuoyancy, 0);
+
+addToRunTimeSelectionTable(otherForceModel, secondaryPhaseInducedBuoyancy, dictionary);
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from components
+secondaryPhaseInducedBuoyancy::secondaryPhaseInducedBuoyancy
+(
+    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_)),
+    diffRho_("diffRho",dimensionSet(1,-3,0,0,0),0.0),
+    gravityFieldName_(propsDict_.lookupOrDefault<word>("gravityFieldName","g")),
+    g_(sm.mesh().lookupObject<uniformDimensionedVectorField> (gravityFieldName_))
+{
+    // density difference of fluid and secondary phase
+    diffRho_.value()=readScalar(propsDict_.lookup ("diffRho"));
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+secondaryPhaseInducedBuoyancy::~secondaryPhaseInducedBuoyancy()
+{}
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+tmp<volVectorField> secondaryPhaseInducedBuoyancy::exportForceField()
+{
+    tmp<volVectorField> tsource
+    (
+        new volVectorField
+        (
+            IOobject
+            (
+                "secondaryPhaseInducedBuoyancy",
+                particleCloud_.mesh().time().timeName(),
+                particleCloud_.mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            particleCloud_.mesh(),
+            dimensionedVector
+            (
+                "zero",dimensionSet(1, -2, -2, 0, 0),vector::zero
+            )
+        )
+    );
+
+    volVectorField& source = tsource.ref();
+
+    source = diffRho_ * (1.0-voidfraction_) * g_;
+
+    return tsource;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

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

@@ -0,0 +1,97 @@
+/*---------------------------------------------------------------------------*\
+License
+    This is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+    This code is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+    You should have received a copy of the GNU General Public License
+    along with this code.  If not, see <http://www.gnu.org/licenses/>.
+
+    Copyright (C) 2015- Thomas Lichtenegger, JKU Linz, Austria
+
+
+Class
+    secondaryPhaseInducedBuoyancy
+
+SourceFiles
+    secondaryPhaseInducedBuoyancy.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef secondaryPhaseInducedBuoyancy_H
+#define secondaryPhaseInducedBuoyancy_H
+
+#include "otherForceModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class secondaryPhaseInducedBuoyancy Declaration
+\*---------------------------------------------------------------------------*/
+
+class secondaryPhaseInducedBuoyancy
+:
+    public otherForceModel
+{
+
+protected:
+
+    // Protected data
+        dictionary propsDict_;
+
+        word voidfractionFieldName_;
+
+        const volScalarField& voidfraction_;
+
+        // density dispersed phase - density continuous phase
+        dimensionedScalar diffRho_;
+
+        word gravityFieldName_;
+
+        const uniformDimensionedVectorField& g_;    // ref to gravity
+
+
+public:
+
+    //- Runtime type information
+    TypeName("secondaryPhaseInducedBuoyancy");
+
+    // Constructors
+
+        //- Construct from components
+        secondaryPhaseInducedBuoyancy
+        (
+            const dictionary& dict,
+            cfdemCloud& sm
+        );
+
+
+    // Destructor
+
+        virtual ~secondaryPhaseInducedBuoyancy();
+
+
+    // Member Functions
+
+
+        tmp<volVectorField> exportForceField();
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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

@@ -46,7 +46,9 @@ weightSecondaryPhase::weightSecondaryPhase
     volfracFieldName_(propsDict_.lookup("volfracFieldName")),
     alpha_(sm.mesh().lookupObject<volScalarField> (volfracFieldName_)),
     rho_("rho",dimensionSet(1,-3,0,0,0),0.0),
-    g_("g",dimensionSet(0,1,-2,0,0),vector(0,0,-9.81))
+    gravityFieldName_(propsDict_.lookupOrDefault<word>("gravityFieldName","g")),
+    g_(sm.mesh().lookupObject<uniformDimensionedVectorField> (gravityFieldName_))
+
 {
     if (propsDict_.found("rho"))
         rho_.value()=readScalar(propsDict_.lookup ("rho"));

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

@@ -52,7 +52,9 @@ protected:
 
         dimensionedScalar rho_;
 
-        const dimensionedVector g_;
+        word gravityFieldName_;
+
+        const uniformDimensionedVectorField& g_;    // ref to gravity
 
 
 public:

src/lagrangian/cfdemParticleComp/Make/files

@@ -35,12 +35,13 @@ $(chemistryModels)/noChemistry/noChemistry.C
 $(chemistryModels)/diffusionCoefficients/diffusionCoefficients.C
 $(chemistryModels)/massTransferCoeff/massTransferCoeff.C
 $(chemistryModels)/reactantPerParticle/reactantPerParticle.C
+$(chemistryModels)/initMultiLayers/initMultiLayers.C
 
 $(energyModels)/energyModel/energyModel.C
 $(energyModels)/energyModel/newEnergyModel.C
 $(energyModels)/heatTransferGunn/heatTransferGunn.C
 $(energyModels)/heatTransferRanzMarshall/heatTransferRanzMarshall.C
-$(energyModels)/heatTransferGranConduction/heatTransferGranConduction.C
+$(energyModels)/heatTransferInterGrain/heatTransferInterGrain.C
 $(energyModels)/reactionHeat/reactionHeat.C
 
 $(thermCondModels)/thermCondModel/thermCondModel.C
@@ -83,6 +84,7 @@ $(forceModels)/Fines/ErgunStatFines.C
 $(forceModels)/granKineticEnergy/granKineticEnergy.C
 $(forceModels)/particleDeformation/particleDeformation.C
 $(forceModels)/turbulentDispersion/turbulentDispersion.C
+$(forceModels)/turbulentVelocityFluctuations/turbulentVelocityFluctuations.C
 
 $(forceSubModels)/forceSubModel/newForceSubModel.C
 $(forceSubModels)/forceSubModel/forceSubModel.C
@@ -95,6 +97,7 @@ $(otherForceModels)/otherForceModel/newOtherForceModel.C
 $(otherForceModels)/gravity/gravity.C
 $(otherForceModels)/weightSecondaryPhase/weightSecondaryPhase.C
 $(otherForceModels)/expParticleForces/expParticleForces.C
+$(otherForceModels)/secondaryPhaseInducedBuoyancy/secondaryPhaseInducedBuoyancy.C
 
 $(probeModels)/probeModel/probeModel.C
 $(probeModels)/probeModel/newProbeModel.C

src/recurrence/recBase/recBase.C

@@ -37,6 +37,18 @@ recBase::recBase
     const fvMesh& mesh
 )
 :
+  regIOobject
+  (
+      IOobject
+      (
+          "recurrenceBase",
+          mesh.time().timeName(),
+          mesh,
+          IOobject::NO_READ,
+          IOobject::NO_WRITE,
+          true
+      )
+  ),
     mesh_(mesh),
     recProperties_
     (
@@ -83,6 +95,7 @@ recBase::recBase
     ),
     couplingSubStep_(recProperties_.lookupOrDefault<label>("couplingSubStep",0))
 {
+
     recModel_        ->  readFieldSeries();
     if (!recStatAnalysis_->suppressMatrixAndPath())
     {
@@ -98,6 +111,80 @@ recBase::recBase
         recModel_        ->  writeRecPath();
     }
 }
+recBase::recBase
+(
+    const fvMesh& mesh,const word recDictName_
+)
+:
+     regIOobject
+  (
+      IOobject
+      (
+          "recurrenceBase",
+          mesh.time().timeName(),
+          mesh,
+          IOobject::NO_READ,
+          IOobject::NO_WRITE,
+          true
+      )
+  ),
+    mesh_(mesh),
+    recProperties_
+    (
+        IOobject
+        (
+            recDictName_,
+            mesh_.time().constant(),
+            mesh_,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        )
+    ),
+    recModel_
+    (
+        recModel::New
+        (
+            recProperties_,
+            *this
+        )
+    ),
+    recNorm_
+    (
+        recNorm::New
+        (
+            recProperties_,
+            *this
+        )
+    ),
+    recPath_
+    (
+        recPath::New
+        (
+            recProperties_,
+            *this
+        )
+    ),
+    recStatAnalysis_
+    (
+        recStatAnalysis::New
+        (
+            recProperties_,
+            *this
+        )
+    ),
+    couplingSubStep_(recProperties_.lookupOrDefault<label>("couplingSubStep",0))
+{
+  recModel_        ->  readFieldSeries();
+  recNorm_         ->  computeRecMatrix();
+  recPath_         ->  getRecPath();
+
+  recModel_        ->  init();
+
+  recModel_        ->  writeRecMatrix();
+  recModel_        ->  writeRecPath();
+}
+
+
 
 // * * * * * * * * * * * * * * * * Destructors  * * * * * * * * * * * * * * //
 recBase::~recBase()

src/recurrence/recBase/recBase.H

@@ -46,6 +46,8 @@ class recStatAnalysis;
 \*---------------------------------------------------------------------------*/
 
 class recBase
+ :
+  public regIOobject
 {
 
 protected:
@@ -65,6 +67,7 @@ protected:
     // in case of substepping: couple passive quantity at this substep
     label couplingSubStep_;
 
+
 public:
 
       // Constructors
@@ -74,6 +77,11 @@ public:
     (
             const fvMesh& mesh
     );
+    //- Construct from mesh and a list of recProperties
+    recBase
+    (
+            const fvMesh& mesh, const word recDictName_
+    );
 
     //- Destructor
     virtual ~recBase();
@@ -94,6 +102,12 @@ public:
 
     label couplingSubStep() const;
 
+   // Dummy function for regIOobject
+    bool writeData(Ostream&) const
+    {
+        return true;
+    }
+
 };
 
 }

src/recurrence/recModel/lruRecModel/lruRecModel.C

@@ -1058,6 +1058,11 @@ void lruRecModel::exportAveragedVolVectorField(volVectorField& smoothfield, word
     smoothfield /= counter;
 }
 
+label lruRecModel::currentTimeIndex() const
+{
+    return virtualTimeIndex;
+}
+
 
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/recurrence/recModel/lruRecModel/lruRecModel.H

@@ -160,6 +160,8 @@ public:
     label numRecFields() const;
     label numDataBaseFields() const;
 
+    label currentTimeIndex() const;
+
     void updateRecFields();
 
     void writeRecMatrix() const;

src/recurrence/recModel/recModel/recModel.H

@@ -224,6 +224,7 @@ public:
     }
 
     virtual void updateRecFields() = 0;
+    virtual label currentTimeIndex() const = 0;
 
     virtual void writeRecMatrix() const = 0;
 

src/recurrence/recModel/standardRecModel/standardRecModel.C

@@ -63,6 +63,7 @@ standardRecModel::standardRecModel
     cumulativeNumRecFields_(),
     totNumRecFields_(0),
     storeAveragedFields_(propsDict_.lookupOrDefault<bool>("storeAveragedFields",false)),
+    checkTimeStep_(propsDict_.lookupOrDefault<bool>("checkTimeStep",true)),
     recurrenceMatrix_(1,scalar(-1.0)),
     timeIndexList_(),
     timeValueList_(),
@@ -127,7 +128,11 @@ standardRecModel::standardRecModel
 
     // check if time steps in databases are consistent
     // if no initial number of time steps has been specified, create path for full runtime immediately
+    if(checkTimeStep_)
+    {
         recTimeStep_ = checkTimeStep();
+    }
+
     if(totRecSteps_ < 0)
     {
         totRecSteps_ = 1 + static_cast<label>( (endTime_-startTime_) / recTimeStep_ );
@@ -196,6 +201,24 @@ void standardRecModel::init()
 {
     recModel::init();
 
+    const objectRegistry& objReg = base_.mesh().thisDb();
+
+    for(int j=0; j<volScalarFieldNames_.size(); j++)
+    {
+        objReg.checkIn(volScalarFieldList_[j][virtualTimeIndex]);
+    }
+
+    for(int j=0; j<volVectorFieldNames_.size(); j++)
+    {
+        objReg.checkIn(volVectorFieldList_[j][virtualTimeIndex]);
+    }
+
+    for(int j=0; j<surfaceScalarFieldNames_.size(); j++)
+    {
+        objReg.checkIn(surfaceScalarFieldList_[j][virtualTimeIndex]);
+    }
+
+
     for(int i = 0; i < numDataBases_; i++)
     {
         if (virtualTimeIndex < cumulativeNumRecFields_[i])
@@ -699,6 +722,27 @@ label standardRecModel::numDataBaseFields() const
 
 void standardRecModel::updateRecFields()
 {
+    // make fields of upcoming virtualTimeIndex available in object registry
+    const objectRegistry& objReg = base_.mesh().thisDb();
+
+    for(int j=0; j<volScalarFieldNames_.size(); j++)
+    {
+        objReg.checkOut(volScalarFieldList_[j][virtualTimeIndex]);
+        objReg.checkIn(volScalarFieldList_[j][virtualTimeIndexNext]);
+    }
+
+    for(int j=0; j<volVectorFieldNames_.size(); j++)
+    {
+        objReg.checkOut(volVectorFieldList_[j][virtualTimeIndex]);
+        objReg.checkIn(volVectorFieldList_[j][virtualTimeIndexNext]);
+    }
+
+    for(int j=0; j<surfaceScalarFieldNames_.size(); j++)
+    {
+        objReg.checkOut(surfaceScalarFieldList_[j][virtualTimeIndex]);
+        objReg.checkIn(surfaceScalarFieldList_[j][virtualTimeIndexNext]);
+    }
+
     virtualTimeIndex = virtualTimeIndexNext;
     virtualTimeIndexNext++;
 
@@ -736,6 +780,10 @@ void standardRecModel::updateRecFields()
     }
 }
 
+label standardRecModel::currentTimeIndex() const
+{
+    return virtualTimeIndex;
+}
 
 void standardRecModel::writeRecMatrix() const
 {

src/recurrence/recModel/standardRecModel/standardRecModel.H

@@ -61,6 +61,7 @@ protected:
     label totNumRecFields_;
     
     bool storeAveragedFields_;
+    bool checkTimeStep_;
 
     // matrix that contains the recurrence ERROR
     SymmetricSquareMatrix<scalar> recurrenceMatrix_;
@@ -151,6 +152,7 @@ public:
     label numDataBaseFields() const;
     
     void updateRecFields();
+    label currentTimeIndex() const;
     
     void writeRecMatrix() const;
 	

src/recurrence/recPath/MarkovPath/MarkovPath.C

@@ -56,6 +56,7 @@ MarkovPath::MarkovPath
 :
     recPath(dict, base),
     propsDict_(dict.subDict(typeName + "Props")),
+    searchMinimum_(propsDict_.lookupOrDefault<bool>("searchMinimum",true)),
     correlationSteps_(readLabel(propsDict_.lookup("correlationSteps"))),
     meanIntervalSteps_(propsDict_.lookupOrDefault<label>("meanIntervalSteps",-1)),
     minIntervalSteps_(propsDict_.lookupOrDefault<label>("minIntervalSteps",0)),
@@ -167,7 +168,7 @@ void MarkovPath::extendPath()
         scalar randInterval = ranGen.scalar01();
 
         label interval = numIntervals_-1;
-        for(int i = numIntervals_-2 ;i >= 0; i--)
+        for(int i = numIntervals_-2; i >= 0; i--)
         {
            if (randInterval < intervalWeightsCumulative_[i]) interval=i;
         }
@@ -176,6 +177,8 @@ void MarkovPath::extendPath()
         if (interval > 0) startLoop = intervalSizesCumulative_[interval-1];
         label endLoop = intervalSizesCumulative_[interval] - meanIntervalSteps_;
 
+        if (searchMinimum_)
+        {
             scalar nextMinimum(GREAT);
             for (label j = startLoop; j <= endLoop; j++)
             {
@@ -186,6 +189,12 @@ void MarkovPath::extendPath()
                     seqStart = j+1;
                 }
             }
+        }
+        else
+        {
+            scalar randIntervalStart = ranGen.scalar01();
+            seqStart = static_cast<label>(startLoop + randIntervalStart*(endLoop - startLoop));
+        }
 
         virtualTimeIndex = seqStart;
     }

src/recurrence/recPath/MarkovPath/MarkovPath.H

@@ -59,6 +59,8 @@ protected:
 
     void weightsNormalization();
 
+    bool searchMinimum_;
+
     // correlation steps determine which temporal neighborhood is excluded from nearest-neighbor search
     label correlationSteps_;
 

tutorials/displacementField/CFD/constant/displacementProperties

@@ -0,0 +1,46 @@
+/*---------------------------------------------------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  1.4                                   |
+|   \\  /    A nd           | Web:      http://www.openfoam.org               |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+
+
+FoamFile
+{
+    version         2.0;
+    format          ascii;
+
+    root            "";
+    case            "";
+    instance        "";
+    local           "";
+
+    class           dictionary;
+    object          displacementProperties;
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+dumpIndexStart 4004000;
+
+dumpIndexEnd 4028000;
+
+dumpIndexDisplacementIncrement 8000;
+
+dumpIndexInputIncrement 4000;
+
+nNeighMin 3;
+
+timePerDisplacementStep 0.1;
+
+timePerInputStep 0.05;
+
+startTime 0.0;
+
+filepath "../DEM/post_disp/dump";
+
+fillEmptyCells false;
+
+// ************************************************************************* //

tutorials/displacementField/CFD/system/blockMeshDict

@@ -0,0 +1,74 @@
+/*--------------------------------*- C++ -*----------------------------------*\
+| =========                 |                                                 |
+| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
+|  \\    /   O peration     | Version:  2.3.0                                 |
+|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
+|    \\/     M anipulation  |                                                 |
+\*---------------------------------------------------------------------------*/
+FoamFile
+{
+    version     2.0;
+    format      ascii;
+    class       dictionary;
+    object      blockMeshDict;
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+convertToMeters 1;
+
+vertices
+(
+    (-0.04 -0.0075 0)
+    (0.04 -0.075 0)
+    (0.04 0.075 0)
+    (-0.04 0.075 0)
+    (-0.04 -0.075 0.25)
+    (0.04 -0.075 0.25)
+    (0.04 0.075 0.25)
+    (-0.04 0.075 0.25)
+);
+
+blocks
+(
+    hex (0 1 2 3 4 5 6 7) (35 6 110) simpleGrading (1 1 1)
+);
+
+edges
+(
+);
+
+boundary
+(
+    walls
+    {
+        type wall;
+        faces
+        (
+            (0 4 7 3)
+            (2 6 5 1)
+            (1 5 4 0)     
+            (3 7 6 2)
+        );
+    }
+    inlet
+    {
+        type patch;
+        faces
+        (
+             (0 3 2 1)
+         );
+    }
+    outlet
+    {
+        type patch;
+        faces
+        (
+            (4 5 6 7)
+         );
+    }
+);
+
+mergePatchPairs
+(
+);
+
+// ************************************************************************* //