Merge branch 'feature/recurrenceLib' of https://github.com/ParticulateFlow/CFDEMcoupling into feature/recurrenceLib

applications/solvers/cfdemSolverRhoSteadyPimple/Make/files

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

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)
     );

applications/solvers/rcfdemSolverRhoSteadyPimple/Make/files

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

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 \

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);

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
+    );
 //===============================

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();

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();
 
 

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);

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

src/lagrangian/cfdemParticle/derived/cfdemCloudRec/cfdemCloudRec.H

@@ -66,6 +66,8 @@ protected:
 
     bool coupleRecForce_;
 
+    bool setForceField_;
+    
 public:
 
     friend class dataExchangeModel;

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()

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