diff --git a/applications/solvers/cfdemSolverRhoSteadyPimple/Make/files b/applications/solvers/cfdemSolverRhoSteadyPimple/Make/files
deleted file mode 100644
index 36bb9c95..00000000
--- a/applications/solvers/cfdemSolverRhoSteadyPimple/Make/files
+++ /dev/null
@@ -1,3 +0,0 @@
-cfdemSolverRhoSteadyPimple.C
-
-EXE=$(CFDEM_APP_DIR)/cfdemSolverRhoSteadyPimple
diff --git a/applications/solvers/cfdemSolverRhoSteadyPimple/EEqn.H b/applications/solvers/rcfdemSolverRhoSteadyPimple/EEqn.H
similarity index 94%
rename from applications/solvers/cfdemSolverRhoSteadyPimple/EEqn.H
rename to applications/solvers/rcfdemSolverRhoSteadyPimple/EEqn.H
index 3dbb5c3f..ebff0160 100644
--- a/applications/solvers/cfdemSolverRhoSteadyPimple/EEqn.H
+++ b/applications/solvers/rcfdemSolverRhoSteadyPimple/EEqn.H
@@ -16,7 +16,7 @@
             fvc::div(phi, K) +
             fvc::div
             (
-                fvc::absolute(phi/fvc::interpolate(rho), voidfraction*U),
+                fvc::absolute(phi/fvc::interpolate(rho), voidfractionRec*U),
                 p,
                 "div(phiv,p)"
             )
@@ -32,7 +32,7 @@
       + addSource
       - Qsource
       - fvm::Sp(QCoeff/Cpv, he)
-      - fvm::laplacian(voidfraction*thCond/Cpv,he)
+      - fvm::laplacian(voidfractionRec*thCond/Cpv,he)
      ==
         fvOptions(rho, he)
     );
diff --git a/applications/solvers/rcfdemSolverRhoSteadyPimple/Make/files b/applications/solvers/rcfdemSolverRhoSteadyPimple/Make/files
new file mode 100644
index 00000000..fac2cca3
--- /dev/null
+++ b/applications/solvers/rcfdemSolverRhoSteadyPimple/Make/files
@@ -0,0 +1,3 @@
+rcfdemSolverRhoSteadyPimple.C
+
+EXE=$(CFDEM_APP_DIR)/rcfdemSolverRhoSteadyPimple
diff --git a/applications/solvers/cfdemSolverRhoSteadyPimple/Make/options b/applications/solvers/rcfdemSolverRhoSteadyPimple/Make/options
similarity index 93%
rename from applications/solvers/cfdemSolverRhoSteadyPimple/Make/options
rename to applications/solvers/rcfdemSolverRhoSteadyPimple/Make/options
index 0377ece5..b2b37fe7 100644
--- a/applications/solvers/cfdemSolverRhoSteadyPimple/Make/options
+++ b/applications/solvers/rcfdemSolverRhoSteadyPimple/Make/options
@@ -15,9 +15,11 @@ EXE_INC = \
     -I$(LIB_SRC)/sampling/lnInclude \
     -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/lnInclude \
     -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/cfdTools \
+    -I$(CFDEM_SRC_DIR)/recurrence/lnInclude \
 
 EXE_LIBS = \
     -L$(CFDEM_LIB_DIR)\
+    -lrecurrence \
     -lcompressibleTransportModels \
     -lfluidThermophysicalModels \
     -lspecie \
diff --git a/applications/solvers/cfdemSolverRhoSteadyPimple/UEqn.H b/applications/solvers/rcfdemSolverRhoSteadyPimple/UEqn.H
similarity index 66%
rename from applications/solvers/cfdemSolverRhoSteadyPimple/UEqn.H
rename to applications/solvers/rcfdemSolverRhoSteadyPimple/UEqn.H
index 8fbd30f2..3b296ac8 100644
--- a/applications/solvers/cfdemSolverRhoSteadyPimple/UEqn.H
+++ b/applications/solvers/rcfdemSolverRhoSteadyPimple/UEqn.H
@@ -4,7 +4,7 @@ particleCloud.otherForces(fOther);
 tmp<fvVectorMatrix> tUEqn
 (
     fvm::div(phi, U)
-  + particleCloud.divVoidfractionTau(U, voidfraction)
+  + particleCloud.divVoidfractionTau(U, voidfractionRec)
   + fvm::Sp(Ksl,U)
   - fOther
  ==
@@ -18,13 +18,15 @@ fvOptions.constrain(UEqn);
 
 if (modelType=="B" || modelType=="Bfull")
 {
-    solve(UEqn == -fvc::grad(p)+ Ksl*Us);
+    solve(UEqn == -fvc::grad(p)+ Ksl*UsRec);
 }
 else
 {
-    solve(UEqn == -voidfraction*fvc::grad(p)+ Ksl*Us);
+    solve(UEqn == -voidfractionRec*fvc::grad(p)+ Ksl*UsRec);
 }
 
+//U.relax();
+
 fvOptions.correct(U);
 
 K = 0.5*magSqr(U);
diff --git a/applications/solvers/cfdemSolverRhoSteadyPimple/createFieldRefs.H b/applications/solvers/rcfdemSolverRhoSteadyPimple/createFieldRefs.H
similarity index 100%
rename from applications/solvers/cfdemSolverRhoSteadyPimple/createFieldRefs.H
rename to applications/solvers/rcfdemSolverRhoSteadyPimple/createFieldRefs.H
diff --git a/applications/solvers/cfdemSolverRhoSteadyPimple/createFields.H b/applications/solvers/rcfdemSolverRhoSteadyPimple/createFields.H
similarity index 87%
rename from applications/solvers/cfdemSolverRhoSteadyPimple/createFields.H
rename to applications/solvers/rcfdemSolverRhoSteadyPimple/createFields.H
index 2219870f..6429ef5f 100644
--- a/applications/solvers/cfdemSolverRhoSteadyPimple/createFields.H
+++ b/applications/solvers/rcfdemSolverRhoSteadyPimple/createFields.H
@@ -51,6 +51,19 @@ Info<< "Reading thermophysical properties\n" << endl;
         mesh
     );
 
+    volScalarField voidfractionRec
+    (
+        IOobject
+        (
+            "voidfractionRec",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        voidfraction
+    );
+
     volScalarField addSource
     (
         IOobject
@@ -58,10 +71,11 @@ Info<< "Reading thermophysical properties\n" << endl;
            "addSource",
            runTime.timeName(),
            mesh,
-           IOobject::MUST_READ,
+           IOobject::READ_IF_PRESENT,
            IOobject::AUTO_WRITE
         ),
-        mesh
+        mesh,
+        dimensionedScalar("zero", dimensionSet(1,-1,-3,0,0,0,0), 0.0)
     );
 
     Info<< "\nCreating fluid-particle heat flux field\n" << endl;
@@ -217,11 +231,11 @@ Info<< "Reading thermophysical properties\n" << endl;
             "Ksl",
             runTime.timeName(),
             mesh,
-            IOobject::MUST_READ,
+            IOobject::READ_IF_PRESENT,
             IOobject::AUTO_WRITE
         ),
-        mesh
-        //dimensionedScalar("0", dimensionSet(1, -3, -1, 0, 0), 1.0)
+        mesh,
+        dimensionedScalar("0", dimensionSet(1, -3, -1, 0, 0), 0.0)
     );
 
 
@@ -239,4 +253,16 @@ Info<< "Reading thermophysical properties\n" << endl;
         mesh
     );
 
+    volVectorField UsRec
+    (
+        IOobject
+        (
+            "UsRec",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        Us
+    );
 //===============================
diff --git a/applications/solvers/cfdemSolverRhoSteadyPimple/pEqn.H b/applications/solvers/rcfdemSolverRhoSteadyPimple/pEqn.H
similarity index 88%
rename from applications/solvers/cfdemSolverRhoSteadyPimple/pEqn.H
rename to applications/solvers/rcfdemSolverRhoSteadyPimple/pEqn.H
index 3d38cc82..b1f352f1 100644
--- a/applications/solvers/cfdemSolverRhoSteadyPimple/pEqn.H
+++ b/applications/solvers/rcfdemSolverRhoSteadyPimple/pEqn.H
@@ -7,11 +7,12 @@ volScalarField rAU(1.0/UEqn.A());
 surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rhoeps*rAU));
 if (modelType=="A")
 {
-    rhorAUf *= fvc::interpolate(voidfraction);
+    rhorAUf *= fvc::interpolate(voidfractionRec);
 }
 volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
+//tUEqn.clear();
 
-surfaceScalarField phiUs("phiUs", fvc::interpolate(rhoeps*rAU*Ksl*Us)& mesh.Sf());
+surfaceScalarField phiUs("phiUs", fvc::interpolate(rhoeps*rAU*Ksl*UsRec)& mesh.Sf());
 
 
 if (pimple.transonic())
@@ -69,11 +70,11 @@ Info<< "rho max/min : " << max(rho).value()
 
 if (modelType=="A")
 {
-    U = HbyA - rAU*(voidfraction*fvc::grad(p)-Ksl*Us);
+    U = HbyA - rAU*(voidfractionRec*fvc::grad(p)-Ksl*UsRec);
 }
 else
 {
-    U = HbyA - rAU*(fvc::grad(p)-Ksl*Us);
+    U = HbyA - rAU*(fvc::grad(p)-Ksl*UsRec);
 }
 
 U.correctBoundaryConditions();
diff --git a/applications/solvers/cfdemSolverRhoSteadyPimple/cfdemSolverRhoSteadyPimple.C b/applications/solvers/rcfdemSolverRhoSteadyPimple/rcfdemSolverRhoSteadyPimple.C
similarity index 80%
rename from applications/solvers/cfdemSolverRhoSteadyPimple/cfdemSolverRhoSteadyPimple.C
rename to applications/solvers/rcfdemSolverRhoSteadyPimple/rcfdemSolverRhoSteadyPimple.C
index d86c4706..65a5efb8 100644
--- a/applications/solvers/cfdemSolverRhoSteadyPimple/cfdemSolverRhoSteadyPimple.C
+++ b/applications/solvers/rcfdemSolverRhoSteadyPimple/rcfdemSolverRhoSteadyPimple.C
@@ -17,11 +17,12 @@ License
     Copyright (C) 2015- Thomas Lichtenegger, JKU Linz, Austria
 
 Application
-    cfdemSolverRhoPimple
+    rcfdemSolverRhoSteadyPimple
 
 Description
-    Transient solver for compressible flow using the flexible PIMPLE (PISO-SIMPLE)
-    algorithm.
+    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.
@@ -36,6 +37,10 @@ Description
 #include "localEulerDdtScheme.H"
 #include "fvcSmooth.H"
 
+#include "cfdemCloudRec.H"
+#include "recBase.H"
+#include "recModel.H"
+#include "recPath.H"
 
 #include "cfdemCloudEnergy.H"
 #include "implicitCouple.H"
@@ -64,9 +69,11 @@ int main(int argc, char *argv[])
     #include "createFvOptions.H"
 
     // create cfdemCloud
-    #include "readGravitationalAcceleration.H"
-    cfdemCloudEnergy particleCloud(mesh);
+  //  #include "readGravitationalAcceleration.H"
+    cfdemCloudRec<cfdemCloudEnergy> particleCloud(mesh);
     #include "checkModelType.H"
+    recBase recurrenceBase(mesh);
+    #include "updateFields.H"
 
     turbulence->validate();
   //        #include "compressibleCourantNo.H"
@@ -74,6 +81,10 @@ int main(int argc, char *argv[])
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+    label recTimeIndex = 0;
+    scalar recTimeStep = recurrenceBase.recM().recTimeStep();
+    scalar startTime = runTime.startTime().value();
+
     Info<< "\nStarting time loop\n" << endl;
 
     while (runTime.run())
@@ -92,6 +103,9 @@ int main(int argc, char *argv[])
         particleCloud.clockM().start(2,"Coupling");
         bool hasEvolved = particleCloud.evolve(voidfraction,Us,U);
 
+//voidfraction = voidfractionRec;
+//Us = UsRec;
+
         if(hasEvolved)
         {
             particleCloud.smoothingM().smoothen(particleCloud.forceM(0).impParticleForces());
@@ -112,7 +126,7 @@ int main(int argc, char *argv[])
         particleCloud.clockM().stop("Coupling");
 
         particleCloud.clockM().start(26,"Flow");
-        volScalarField rhoeps("rhoeps",rho*voidfraction);
+        volScalarField rhoeps("rhoeps",rho*voidfractionRec);
         while (pimple.loop())
         {
             // if needed, perform drag update here
@@ -132,7 +146,7 @@ int main(int argc, char *argv[])
             {
                 // besides this pEqn, OF offers a "pimple consistent"-option
                 #include "pEqn.H"
-                rhoeps=rho*voidfraction;
+                rhoeps=rho*voidfractionRec;
             }
 
             if (pimple.turbCorr())
@@ -146,6 +160,15 @@ int main(int argc, char *argv[])
         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();
 
 
diff --git a/applications/solvers/rcfdemSolverRhoSteadyPimple/updateFields.H b/applications/solvers/rcfdemSolverRhoSteadyPimple/updateFields.H
new file mode 100644
index 00000000..194df94c
--- /dev/null
+++ b/applications/solvers/rcfdemSolverRhoSteadyPimple/updateFields.H
@@ -0,0 +1,8 @@
+// is it neccessary to extend recurrence path?
+if(recurrenceBase.recM().endOfPath())
+{
+    recurrenceBase.extendPath();
+}
+
+recurrenceBase.recM().exportVolScalarField("voidfraction",voidfractionRec);
+recurrenceBase.recM().exportVolVectorField("UsMean",UsRec);
diff --git a/src/lagrangian/cfdemParticle/Make/files b/src/lagrangian/cfdemParticle/Make/files
index 1787e6a8..2f572ada 100644
--- a/src/lagrangian/cfdemParticle/Make/files
+++ b/src/lagrangian/cfdemParticle/Make/files
@@ -158,8 +158,6 @@ $(dataExchangeModels)/oneWayVTK/oneWayVTK.C
 $(dataExchangeModels)/twoWayFiles/twoWayFiles.C
 $(dataExchangeModels)/noDataExchange/noDataExchange.C
 $(dataExchangeModels)/twoWayMPI/twoWayMPI.C
-$(dataExchangeModels)/twoWayOne2One/twoWayOne2One.C
-$(dataExchangeModels)/twoWayOne2One/one2one.C
 
 $(averagingModels)/averagingModel/averagingModel.C
 $(averagingModels)/averagingModel/newAveragingModel.C
diff --git a/src/lagrangian/cfdemParticle/derived/cfdemCloudRec/cfdemCloudRec.H b/src/lagrangian/cfdemParticle/derived/cfdemCloudRec/cfdemCloudRec.H
index a861f81a..28ba05f0 100644
--- a/src/lagrangian/cfdemParticle/derived/cfdemCloudRec/cfdemCloudRec.H
+++ b/src/lagrangian/cfdemParticle/derived/cfdemCloudRec/cfdemCloudRec.H
@@ -65,6 +65,8 @@ protected:
     bool coupleRecFluc_;
 
     bool coupleRecForce_;
+
+    bool setForceField_;
     
 public:
 
diff --git a/src/lagrangian/cfdemParticle/derived/cfdemCloudRec/cfdemCloudRecI.H b/src/lagrangian/cfdemParticle/derived/cfdemCloudRec/cfdemCloudRecI.H
index 48679b7b..eda1d29d 100644
--- a/src/lagrangian/cfdemParticle/derived/cfdemCloudRec/cfdemCloudRecI.H
+++ b/src/lagrangian/cfdemParticle/derived/cfdemCloudRec/cfdemCloudRecI.H
@@ -40,7 +40,8 @@ cfdemCloudRec<baseCloud>::cfdemCloudRec
 :
     baseCloud(mesh),
     coupleRecFluc_(true),
-    coupleRecForce_(true)
+    coupleRecForce_(true),
+    setForceField_(false)
 {
     if (baseCloud::couplingProperties().found("coupleRecFluc"))
     {
@@ -51,6 +52,11 @@ cfdemCloudRec<baseCloud>::cfdemCloudRec
     {
         coupleRecForce_ = readBool(baseCloud::couplingProperties().lookup("coupleRecForce"));
     }
+
+    if (baseCloud::couplingProperties().found("setForceField"))
+    {
+        setForceField_ = readBool(baseCloud::couplingProperties().lookup("setForceField"));
+    }
 }
 
 // * * * * * * * * * * * * * * * * Destructors  * * * * * * * * * * * * * * //
@@ -79,13 +85,22 @@ void cfdemCloudRec<baseCloud>::setForces()
 {
     baseCloud::resetArray(baseCloud::particleConvVel_,baseCloud::numberOfParticles(),3);
     baseCloud::resetArray(baseCloud::particleFlucVel_,baseCloud::numberOfParticles(),3);
+
+    baseCloud::resetArray(baseCloud::impForces_,baseCloud::numberOfParticles(),3);
+    baseCloud::resetArray(baseCloud::expForces_,baseCloud::numberOfParticles(),3);
+
     baseCloud::resetArray(baseCloud::DEMForces_,baseCloud::numberOfParticles(),3);
     for (int i=0;i<baseCloud::nrForceModels();i++) baseCloud::forceM(i).setForce();
 }
 
 template <class baseCloud>
 void cfdemCloudRec<baseCloud>::setParticleForceField()
-{}
+{
+    if (setForceField_)
+    {
+        baseCloud::setParticleForceField();
+    }
+}
 
 template <class baseCloud>
 void cfdemCloudRec<baseCloud>::setVectorAverages()
diff --git a/src/lagrangian/cfdemParticleComp/Make/files b/src/lagrangian/cfdemParticleComp/Make/files
index 5bbd2e0a..72e43a2b 100644
--- a/src/lagrangian/cfdemParticleComp/Make/files
+++ b/src/lagrangian/cfdemParticleComp/Make/files
@@ -133,6 +133,7 @@ $(momCoupleModels)/momCoupleModel/newMomCoupleModel.C
 $(momCoupleModels)/explicitCouple/explicitCouple.C
 $(momCoupleModels)/implicitCouple/implicitCouple.C
 $(momCoupleModels)/noCouple/noCouple.C
+$(momCoupleModels)/deactivateCouple/deactivateCouple.C
 
 $(regionModels)/regionModel/regionModel.C
 $(regionModels)/regionModel/newRegionModel.C
@@ -144,8 +145,6 @@ $(dataExchangeModels)/oneWayVTK/oneWayVTK.C
 $(dataExchangeModels)/twoWayFiles/twoWayFiles.C
 $(dataExchangeModels)/noDataExchange/noDataExchange.C
 $(dataExchangeModels)/twoWayMPI/twoWayMPI.C
-$(dataExchangeModels)/twoWayOne2One/twoWayOne2One.C
-$(dataExchangeModels)/twoWayOne2One/one2one.C
 
 $(averagingModels)/averagingModel/averagingModel.C
 $(averagingModels)/averagingModel/newAveragingModel.C