Merge with OF4x; added some heat transfer related models.

applications/solvers/cfdemSolverIB/Make/options

@@ -1,24 +1,27 @@
 include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
 
 EXE_INC = \
-    -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
+    -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$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/lnInclude \
-    -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/sampling/lnInclude \
     -I$(LIB_SRC)/dynamicFvMesh/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/dynamicFvMesh/lnInclude \
     -I$(LIB_SRC)/dynamicMesh/dynamicMesh/lnInclude \
     -I$(LIB_SRC)/fvOptions/lnInclude
+
 EXE_LIBS = \
     -L$(CFDEM_LIB_DIR)\
-    -lincompressibleRASModels \
-    -lincompressibleLESModels \
+    -lturbulenceModels \
+    -lincompressibleTurbulenceModels \
     -lincompressibleTransportModels \
     -lfiniteVolume \
+    -lmeshTools \
     -ldynamicFvMesh \
     -ldynamicMesh \
     -lfvOptions \

applications/solvers/cfdemSolverIB/cfdemSolverIB.C

@@ -39,7 +39,8 @@ Contributions
 
 #include "fvCFD.H"
 #include "singlePhaseTransportModel.H"
-#include "turbulenceModel.H"
+#include "turbulentTransportModel.H"
+#include "pisoControl.H"
 
 #include "cfdemCloudIB.H"
 #include "implicitCouple.H"
@@ -52,10 +53,9 @@ Contributions
 
 #include "cellSet.H"
 
-#if defined(version22)
-    #include "meshToMeshNew.H"
-    #include "fvIOoptionList.H"
-#endif
+#include "meshToMeshNew.H"
+#include "fvIOoptionList.H"
+
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -66,14 +66,14 @@ int main(int argc, char *argv[])
     #include "createTime.H"
 
     #include "createDynamicFvMesh.H"
+  
+    #include "createControl.H"
 
     #include "createFields.H"
 
     #include "initContinuityErrs.H"
 
-    #if defined(version22)
     #include "createFvOptions.H"
-    #endif
 
     // create cfdemCloud
     #include "readGravitationalAcceleration.H"
@@ -91,7 +91,6 @@ int main(int argc, char *argv[])
         interFace = mag(mesh.lookupObject<volScalarField>("voidfractionNext"));
         mesh.update(); //dyM
 
-        #include "readPISOControls.H"
         #include "CourantNo.H"
 
         // do particle stuff
@@ -107,45 +106,33 @@ int main(int argc, char *argv[])
                 fvm::ddt(voidfraction,U)
               + fvm::div(phi, U)
               + turbulence->divDevReff(U)
-                #if defined(version22)
                 ==
                 fvOptions(U)
-                #endif
             );
 
             UEqn.relax();
 
-            #if defined(version22)
             fvOptions.constrain(UEqn);
-            #endif
 
-            if (momentumPredictor)
+            if (piso.momentumPredictor())
             {
                 solve(UEqn == -fvc::grad(p));
             }
 
             // --- PISO loop
-            for (int corr=0; corr<nCorr; corr++)
+           while (piso.correct())
             {
                 volScalarField rUA = 1.0/UEqn.A();
                 surfaceScalarField rUAf(fvc::interpolate(rUA));
 
                 U = rUA*UEqn.H();
-                #ifdef version23
+
                 phi = (fvc::interpolate(U) & mesh.Sf())
                     + rUAf*fvc::ddtCorr(U, phi);
-                #else
-                phi = (fvc::interpolate(U) & mesh.Sf())
-                    + fvc::ddtPhiCorr(rUA, U, phi);
-                #endif
+
                 adjustPhi(phi, U, p);
 
-                #if defined(version22)
-                fvOptions.relativeFlux(phi);
-                #endif
-
-                // Non-orthogonal pressure corrector loop
-                for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+                while (piso.correct())
                 {
                     // Pressure corrector
 
@@ -156,20 +143,9 @@ int main(int argc, char *argv[])
 
                     pEqn.setReference(pRefCell, pRefValue);
 
-                    if
-                    (
-                        corr == nCorr-1
-                     && nonOrth == nNonOrthCorr
-                    )
-                    {
-                        pEqn.solve(mesh.solver("pFinal"));
-                    }
-                    else
-                    {
-                        pEqn.solve();
-                    }
+                    pEqn.solve(mesh.solver(p.select(piso.finalInnerIter())));
 
-                    if (nonOrth == nNonOrthCorr)
+                    if (piso.finalNonOrthogonalIter())
                     {
                         phi -= pEqn.flux();
                     }

applications/solvers/cfdemSolverPiso/Make/options

@@ -1,19 +1,22 @@
 include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
 
 EXE_INC = \
-    -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
+    -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$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/lnInclude \
     -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/cfdTools \
 
 EXE_LIBS = \
     -L$(CFDEM_LIB_DIR)\
-    -lincompressibleRASModels \
-    -lincompressibleLESModels \
+    -lturbulenceModels \
+    -lincompressibleTurbulenceModels \
     -lincompressibleTransportModels \
     -lfiniteVolume \
+    -lmeshTools \
     -l$(CFDEM_LIB_NAME) \
      $(CFDEM_ADD_LIB_PATHS) \
      $(CFDEM_ADD_LIBS)

applications/solvers/cfdemSolverPiso/UEqn.H

@@ -0,0 +1,23 @@
+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)
+  ==
+  - fvm::Sp(Ksl/rho,U)
+);
+
+UEqn.relax();
+
+fvOptions.constrain(UEqn);
+
+if (piso.momentumPredictor() && (modelType=="B" || modelType=="Bfull"))
+{
+    solve(UEqn == - fvc::grad(p) + Ksl/rho*Us);
+    fvOptions.correct(U);		    
+}
+else if (piso.momentumPredictor())
+{
+    solve(UEqn == - voidfraction*fvc::grad(p) + Ksl/rho*Us);
+    fvOptions.correct(U);
+}

applications/solvers/cfdemSolverPiso/cfdemSolverPiso.C

@@ -36,7 +36,9 @@ Description
 
 #include "fvCFD.H"
 #include "singlePhaseTransportModel.H"
-#include "turbulenceModel.H"
+#include "turbulentTransportModel.H"
+#include "pisoControl.H"
+#include "fvOptions.H"
 
 #include "cfdemCloud.H"
 #include "implicitCouple.H"
@@ -51,7 +53,9 @@ int main(int argc, char *argv[])
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
+    #include "createControl.H"
     #include "createFields.H"
+    #include "createFvOptions.H"
     #include "initContinuityErrs.H"
 
     // create cfdemCloud
@@ -67,7 +71,6 @@ int main(int argc, char *argv[])
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        #include "readPISOControls.H"
         #include "CourantNo.H"
 
         // do particle stuff
@@ -85,7 +88,7 @@ int main(int argc, char *argv[])
 
        //Force Checks
        vector fTotal(0,0,0);
-       vector fImpTotal = sum(mesh.V()*Ksl.internalField()*(Us.internalField()-U.internalField()));
+       vector fImpTotal = sum(mesh.V()*Ksl.internalField()*(Us.internalField()-U.internalField())).value();
        reduce(fImpTotal, sumOp<vector>());
        Info << "TotalForceExp: " << fTotal << endl;
        Info << "TotalForceImp: " << fImpTotal << endl;
@@ -100,95 +103,19 @@ int main(int argc, char *argv[])
             // Pressure-velocity PISO corrector
             {
                 // Momentum predictor
-                fvVectorMatrix UEqn
-                (
-                    fvm::ddt(voidfraction,U) - fvm::Sp(fvc::ddt(voidfraction),U)
-                  + fvm::div(phi,U) - fvm::Sp(fvc::div(phi),U)
-//                + turbulence->divDevReff(U)
-                  + particleCloud.divVoidfractionTau(U, voidfraction)
-                  ==
-                  - fvm::Sp(Ksl/rho,U)
-                );
-
-                UEqn.relax();
-                if (momentumPredictor && (modelType=="B" || modelType=="Bfull"))
-                    solve(UEqn == - fvc::grad(p) + Ksl/rho*Us);
-                else if (momentumPredictor)
-                    solve(UEqn == - voidfraction*fvc::grad(p) + Ksl/rho*Us);
+                 #include "UEqn.H"
 
                 // --- PISO loop
 
-                //for (int corr=0; corr<nCorr; corr++)
-                int nCorrSoph = nCorr + 5. * (1. - particleCloud.dataExchangeM().timeStepFraction());
-
-                for (int corr=0; corr<nCorrSoph; corr++)
+                while (piso.correct())
                 {
-                    volScalarField rUA = 1.0/UEqn.A();
-
-                    surfaceScalarField rUAf("(1|A(U))", fvc::interpolate(rUA));
-                    volScalarField rUAvoidfraction("(voidfraction2|A(U))",rUA*voidfraction);
-                    surfaceScalarField rUAfvoidfraction("(voidfraction2|A(U)F)", fvc::interpolate(rUAvoidfraction));
-
-                    U = rUA*UEqn.H();
-
-                    #ifdef version23
-                    phi = ( fvc::interpolate(U*voidfraction) & mesh.Sf() )
-                        + rUAfvoidfraction*fvc::ddtCorr(U, phi);
-                    #else
-                    phi = ( fvc::interpolate(U*voidfraction) & mesh.Sf() )
-                        + fvc::ddtPhiCorr(rUAvoidfraction, U, phi);
-                    #endif
-                    surfaceScalarField phiS(fvc::interpolate(Us*voidfraction) & mesh.Sf());
-                    surfaceScalarField phiGes = phi + rUAf*(fvc::interpolate(Ksl/rho) * phiS);
-
-                    if (modelType=="A")
-                        rUAvoidfraction = volScalarField("(voidfraction2|A(U))",rUA*voidfraction*voidfraction);
-
-                    // Non-orthogonal pressure corrector loop
-                    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
-                    {
-                        // Pressure corrector
-                        fvScalarMatrix pEqn
-                        (
-                            fvm::laplacian(rUAvoidfraction, p) == fvc::div(phiGes) + particleCloud.ddtVoidfraction()
-                        );
-                        pEqn.setReference(pRefCell, pRefValue);
-
-                        if
-                        (
-                            corr == nCorr-1
-                         && nonOrth == nNonOrthCorr
-                        )
-                        {
-                            pEqn.solve(mesh.solver("pFinal"));
-                        }
-                        else
-                        {
-                            pEqn.solve();
-                        }
-
-                        if (nonOrth == nNonOrthCorr)
-                        {
-                            phiGes -= pEqn.flux();
-                            phi = phiGes;
-                        }
-
-                    } // end non-orthogonal corrector loop
-
-                    #include "continuityErrorPhiPU.H"
-
-                    if (modelType=="B" || modelType=="Bfull")
-                        U -= rUA*fvc::grad(p) - Ksl/rho*Us*rUA;
-                    else
-                        U -= voidfraction*rUA*fvc::grad(p) - Ksl/rho*Us*rUA;
-
-                    U.correctBoundaryConditions();
-
-                } // end piso loop
+                    #include "pEqn.H"
+                }
             }
 
+            laminarTransport.correct();
             turbulence->correct();
-        }// end solveFlow
+        }
         else
         {
             Info << "skipping flow solution." << endl;

applications/solvers/cfdemSolverPiso/createFields.H

@@ -122,3 +122,5 @@ surfaceScalarField phi
     (
         incompressible::turbulenceModel::New(U, phi, laminarTransport)
     );
+
+#include "createMRF.H"

applications/solvers/cfdemSolverPiso/pEqn.H

@@ -0,0 +1,58 @@
+volScalarField rAU = 1.0/UEqn.A();
+
+surfaceScalarField rAUf("(1|A(U))", fvc::interpolate(rAU));
+
+volScalarField rAUvoidfraction("(voidfraction2|A(U))",rAU*voidfraction);
+
+if (modelType=="A")
+    rAUvoidfraction *= voidfraction;
+
+volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
+
+surfaceScalarField phiHbyA("phiHbyA", fvc::interpolate(voidfraction)*fvc::flux(HbyA) );
+
+volVectorField Uvoidfraction("(Uvoidfraction)", U * voidfraction);
+
+surfaceScalarField phiS(fvc::interpolate(Us*voidfraction) & mesh.Sf());
+
+phi = phiHbyA + rAUf*(fvc::interpolate(Ksl/rho) * phiS);
+
+// rotating frames of references not tested yet
+MRF.makeRelative(phi);
+
+// adjustment of phi (only in cases w.o. p boundary conditions) not tested yet
+adjustPhi(phi, U, p);
+
+// Update the pressure BCs to ensure flux consistency
+constrainPressure(p, Uvoidfraction, phiHbyA, rAUvoidfraction, MRF);
+
+
+// Non-orthogonal pressure corrector loop
+while (piso.correctNonOrthogonal())
+{
+    // Pressure corrector
+  
+    fvScalarMatrix pEqn
+    (
+        fvm::laplacian(rAUvoidfraction, p) == fvc::div(phi) + particleCloud.ddtVoidfraction()
+    );
+    
+    pEqn.setReference(pRefCell, pRefValue);
+
+    pEqn.solve(mesh.solver(p.select(piso.finalInnerIter())));
+
+    if (piso.finalNonOrthogonalIter())
+    {
+        phi -= pEqn.flux();
+    }
+}
+
+#include "continuityErrorPhiPU.H"
+
+if (modelType=="B" || modelType=="Bfull")
+    U = HbyA - rAU*fvc::grad(p) + Ksl/rho*Us*rAU;
+else
+    U = HbyA - voidfraction*rAU*fvc::grad(p) + Ksl/rho*Us*rAU;
+
+U.correctBoundaryConditions();
+fvOptions.correct(U);

applications/solvers/cfdemSolverPisoMS/Make/options

@@ -1,16 +1,23 @@
+include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
+
 EXE_INC = \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/lnInclude \
     -I ../cfdemSolverPiso \
-    -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude \
     -I$(LIB_SRC)/transportModels \
     -I$(LIB_SRC)/transportModels/incompressible/singlePhaseTransportModel \
-    -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/cfdTools \
 
 EXE_LIBS = \
     -L$(CFDEM_LIB_DIR)\
-    -lincompressibleRASModels \
-    -lincompressibleLESModels \
+    -lturbulenceModels \
+    -lincompressibleTurbulenceModels \
     -lincompressibleTransportModels \
     -lfiniteVolume \
-    -l$(CFDEM_LIB_NAME)
+    -lmeshTools \
+    -l$(CFDEM_LIB_NAME) \
+     $(CFDEM_ADD_LIB_PATHS) \
+     $(CFDEM_ADD_LIBS)

applications/solvers/cfdemSolverPisoMS/cfdemSolverPisoMS.C

@@ -25,23 +25,26 @@ License
     along with CFDEMcoupling.  If not, see <http://www.gnu.org/licenses/>.
 
 Application
-    cfdemSolverPisoMS
+    cfdemSolverPiso
 
 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, 
+    The code is an evolution of the solver pisoFoam in OpenFOAM(R) 1.6,
     where additional functionality for CFD-DEM coupling is added.
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
 #include "singlePhaseTransportModel.H"
-#include "turbulenceModel.H"
+#include "turbulentTransportModel.H"
+#include "pisoControl.H"
+#include "fvOptions.H"
 
 #include "cfdemCloudMS.H"
 #include "implicitCouple.H"
 #include "clockModel.H"
 #include "smoothingModel.H"
+#include "forceModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -50,7 +53,9 @@ int main(int argc, char *argv[])
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
+    #include "createControl.H"
     #include "createFields.H"
+    #include "createFvOptions.H"
     #include "initContinuityErrs.H"
 
     // create cfdemCloud
@@ -59,118 +64,62 @@ int main(int argc, char *argv[])
     #include "checkModelType.H"
 
     // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
+    Info<< "\nStarting time loop\n" << endl;
     while (runTime.loop())
     {
-        Info<< "\nStarting time loop\n" << endl;
-            particleCloud.clockM().start(1,"Global");
+        particleCloud.clockM().start(1,"Global");
 
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        #include "readPISOControls.H"
         #include "CourantNo.H"
 
         // do particle stuff
         particleCloud.clockM().start(2,"Coupling");
-        particleCloud.evolve(voidfraction,Us,U);
-      
+        bool hasEvolved = particleCloud.evolve(voidfraction,Us,U);
+
+        if(hasEvolved)
+        {
+            particleCloud.smoothingM().smoothen(particleCloud.forceM(0).impParticleForces());
+        }
+    
         Info << "update Ksl.internalField()" << endl;
         Ksl = particleCloud.momCoupleM(0).impMomSource();
-        particleCloud.smoothingM().smoothen(Ksl);
         Ksl.correctBoundaryConditions();
 
+       //Force Checks
+       vector fTotal(0,0,0);
+       vector fImpTotal = sum(mesh.V()*Ksl.internalField()*(Us.internalField()-U.internalField())).value();
+       reduce(fImpTotal, sumOp<vector>());
+       Info << "TotalForceExp: " << fTotal << endl;
+       Info << "TotalForceImp: " << fImpTotal << endl;
+
         #include "solverDebugInfo.H"
         particleCloud.clockM().stop("Coupling");
 
         particleCloud.clockM().start(26,"Flow");
-        // Pressure-velocity PISO corrector
+
+        if(particleCloud.solveFlow())
         {
-            // Momentum predictor
-            fvVectorMatrix UEqn
-            (
-                fvm::ddt(voidfraction,U) //particleCloud.ddtVoidfractionU(U,voidfraction) //
-              + fvm::div(phi, U)
-//              + turbulence->divDevReff(U)
-              + particleCloud.divVoidfractionTau(U, voidfraction)
-             ==
-              - fvm::Sp(Ksl/rho,U)
-            );
-
-            if (modelType=="B")
-                UEqn == - fvc::grad(p) + Ksl/rho*Us;
-            else
-                UEqn == - voidfraction*fvc::grad(p) + Ksl/rho*Us;
-
-            UEqn.relax();
-
-            if (momentumPredictor)
-                solve(UEqn);
-
-            // --- PISO loop
-
-            //for (int corr=0; corr<nCorr; corr++)
-            int nCorrSoph = nCorr + 5. * (1. - particleCloud.dataExchangeM().timeStepFraction());
-
-            for (int corr=0; corr<nCorrSoph; corr++)
+            // Pressure-velocity PISO corrector
             {
-                volScalarField rUA = 1.0/UEqn.A();
+                // Momentum predictor
+                 #include "UEqn.H"
 
-                surfaceScalarField rUAf("(1|A(U))", fvc::interpolate(rUA));
-                volScalarField rUAvoidfraction("(voidfraction2|A(U))",rUA*voidfraction);
+                // --- PISO loop
 
-                U = rUA*UEqn.H();
-
-                phi = (fvc::interpolate(U*voidfraction) & mesh.Sf() );
-                     //+ fvc::ddtPhiCorr(rUAvoidfraction, U, phi);
-                surfaceScalarField phiS(fvc::interpolate(Us*voidfraction) & mesh.Sf());
-                surfaceScalarField phiGes = phi + rUAf*(fvc::interpolate(Ksl/rho) * phiS);
-
-                if (modelType=="A")
-                    rUAvoidfraction = volScalarField("(voidfraction2|A(U))",rUA*voidfraction*voidfraction);
-
-                // Non-orthogonal pressure corrector loop
-                for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
+                while (piso.correct())
                 {
-                    // Pressure corrector
-                    fvScalarMatrix pEqn
-                    (
-                        fvm::laplacian(rUAvoidfraction, p) == fvc::div(phiGes) + particleCloud.ddtVoidfraction()
-                    );
-                    pEqn.setReference(pRefCell, pRefValue);
+                    #include "pEqn.H"
+                }
+            }
 
-                    if
-                    (
-                        corr == nCorr-1
-                     && nonOrth == nNonOrthCorr
-                    )
-                    {
-                        pEqn.solve(mesh.solver("pFinal"));
-                    }
-                    else
-                    {
-                        pEqn.solve();
-                    }
-
-                    if (nonOrth == nNonOrthCorr)
-                    {
-                        phiGes -= pEqn.flux();
-                    }
-
-                } // end non-orthogonal corrector loop
-
-                #include "continuityErrorPhiPU.H"
-
-                if (modelType=="B")
-                    U -= rUA*fvc::grad(p) - Ksl/rho*Us*rUA;
-                else
-                    U -= voidfraction*rUA*fvc::grad(p) - Ksl/rho*Us*rUA;
-
-                U.correctBoundaryConditions();
-
-            } // end piso loop
+            laminarTransport.correct();
+            turbulence->correct();
+        }
+        else
+        {
+            Info << "skipping flow solution." << endl;
         }
-
-        turbulence->correct();
 
         runTime.write();
 
@@ -183,7 +132,7 @@ int main(int argc, char *argv[])
     }
 
     Info<< "End\n" << endl;
-    
+
     return 0;
 }
 

applications/solvers/cfdemSolverPisoScalar/Make/options

@@ -1,19 +1,23 @@
 include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
 
 EXE_INC = \
-    -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
+    -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$(LIB_SRC)/finiteVolume/lnInclude \
+    -I ../cfdemSolverPiso \
     -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/lnInclude \
     -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/cfdTools \
 
 EXE_LIBS = \
     -L$(CFDEM_LIB_DIR)\
-    -lincompressibleRASModels \
-    -lincompressibleLESModels \
+    -lturbulenceModels \
+    -lincompressibleTurbulenceModels \
     -lincompressibleTransportModels \
     -lfiniteVolume \
+    -lmeshTools \
     -l$(CFDEM_LIB_NAME) \
      $(CFDEM_ADD_LIB_PATHS) \
      $(CFDEM_ADD_LIBS)

applications/solvers/cfdemSolverPisoScalar/TEqn.H

@@ -0,0 +1,15 @@
+ // get scalar source from DEM        
+        particleCloud.forceM(1).manipulateScalarField(Tsource);
+        Tsource.correctBoundaryConditions();
+
+        // solve scalar transport equation
+        fvScalarMatrix TEqn
+        (
+           fvm::ddt(voidfraction,T) - fvm::Sp(fvc::ddt(voidfraction),T)
+         + fvm::div(phi, T) - fvm::Sp(fvc::div(phi),T)
+         - fvm::laplacian(DT*voidfraction, T)
+         ==
+           Tsource
+        );
+        TEqn.relax();
+        TEqn.solve();

applications/solvers/cfdemSolverPisoScalar/cfdemSolverPisoScalar.C

@@ -25,7 +25,7 @@ License
     along with CFDEMcoupling.  If not, see <http://www.gnu.org/licenses/>.
 
 Application
-    cfdemSolverPisoScalar
+    cfdemSolverPiso
 
 Description
     Transient solver for incompressible flow.
@@ -36,10 +36,13 @@ Description
 
 #include "fvCFD.H"
 #include "singlePhaseTransportModel.H"
-#include "turbulenceModel.H"
+#include "turbulentTransportModel.H"
+#include "pisoControl.H"
+#include "fvOptions.H"
 
 #include "cfdemCloud.H"
 #include "implicitCouple.H"
+#include "clockModel.H"
 #include "smoothingModel.H"
 #include "forceModel.H"
 
@@ -50,7 +53,9 @@ int main(int argc, char *argv[])
     #include "setRootCase.H"
     #include "createTime.H"
     #include "createMesh.H"
+    #include "createControl.H"
     #include "createFields.H"
+    #include "createFvOptions.H"
     #include "initContinuityErrs.H"
 
     // create cfdemCloud
@@ -62,12 +67,14 @@ int main(int argc, char *argv[])
     Info<< "\nStarting time loop\n" << endl;
     while (runTime.loop())
     {
+        particleCloud.clockM().start(1,"Global");
+
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
-        #include "readPISOControls.H"
         #include "CourantNo.H"
 
         // do particle stuff
+        particleCloud.clockM().start(2,"Coupling");
         bool hasEvolved = particleCloud.evolve(voidfraction,Us,U);
 
         if(hasEvolved)
@@ -79,119 +86,38 @@ int main(int argc, char *argv[])
         Ksl = particleCloud.momCoupleM(0).impMomSource();
         Ksl.correctBoundaryConditions();
 
+       //Force Checks
+       vector fTotal(0,0,0);
+       vector fImpTotal = sum(mesh.V()*Ksl.internalField()*(Us.internalField()-U.internalField())).value();
+       reduce(fImpTotal, sumOp<vector>());
+       Info << "TotalForceExp: " << fTotal << endl;
+       Info << "TotalForceImp: " << fImpTotal << endl;
 
         #include "solverDebugInfo.H"
+        particleCloud.clockM().stop("Coupling");
 
-        // get scalar source from DEM        
-        particleCloud.forceM(1).manipulateScalarField(Tsource);
-        Tsource.correctBoundaryConditions();
-
-        // solve scalar transport equation
-        fvScalarMatrix TEqn
-        (
-           fvm::ddt(voidfraction,T) - fvm::Sp(fvc::ddt(voidfraction),T)
-         + fvm::div(phi, T) - fvm::Sp(fvc::div(phi),T)
-         - fvm::laplacian(DT*voidfraction, T)
-         ==
-           Tsource
-        );
-        TEqn.relax();
-        TEqn.solve();
+        particleCloud.clockM().start(26,"Flow");
+	
+        #include "TEqn.H"
 
         if(particleCloud.solveFlow())
         {
             // Pressure-velocity PISO corrector
             {
                 // Momentum predictor
-                fvVectorMatrix UEqn
-                (
-                    fvm::ddt(voidfraction,U) - fvm::Sp(fvc::ddt(voidfraction),U)
-                  + fvm::div(phi,U) - fvm::Sp(fvc::div(phi),U)
-//                + turbulence->divDevReff(U)
-                  + particleCloud.divVoidfractionTau(U, voidfraction)
-                  ==
-                  - fvm::Sp(Ksl/rho,U)
-                );
-
-                UEqn.relax();
-                if (momentumPredictor && (modelType=="B" || modelType=="Bfull"))
-                    solve(UEqn == - fvc::grad(p) + Ksl/rho*Us);
-                else if (momentumPredictor)
-                    solve(UEqn == - voidfraction*fvc::grad(p) + Ksl/rho*Us);
+                 #include "UEqn.H"
 
                 // --- PISO loop
 
-                //for (int corr=0; corr<nCorr; corr++)
-                int nCorrSoph = nCorr + 5. * (1. - particleCloud.dataExchangeM().timeStepFraction());
-
-                for (int corr=0; corr<nCorrSoph; corr++)
+                while (piso.correct())
                 {
-                    volScalarField rUA = 1.0/UEqn.A();
-
-                    surfaceScalarField rUAf("(1|A(U))", fvc::interpolate(rUA));
-                    volScalarField rUAvoidfraction("(voidfraction2|A(U))",rUA*voidfraction);
-                    surfaceScalarField rUAfvoidfraction("(voidfraction2|A(U)F)", fvc::interpolate(rUAvoidfraction));
-
-                    U = rUA*UEqn.H();
-
-                    #ifdef version23
-                    phi = ( fvc::interpolate(U*voidfraction) & mesh.Sf() )
-                        + rUAfvoidfraction*fvc::ddtCorr(U, phi);
-                    #else
-                    phi = ( fvc::interpolate(U*voidfraction) & mesh.Sf() )
-                        + fvc::ddtPhiCorr(rUAvoidfraction, U, phi);
-                    #endif
-                    surfaceScalarField phiS(fvc::interpolate(Us*voidfraction) & mesh.Sf());
-                    surfaceScalarField phiGes = phi + rUAf*(fvc::interpolate(Ksl/rho) * phiS);
-
-                    if (modelType=="A")
-                        rUAvoidfraction = volScalarField("(voidfraction2|A(U))",rUA*voidfraction*voidfraction);
-
-                    // Non-orthogonal pressure corrector loop
-                    for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
-                    {
-                        // Pressure corrector
-                        fvScalarMatrix pEqn
-                        (
-                            fvm::laplacian(rUAvoidfraction, p) == fvc::div(phiGes) + particleCloud.ddtVoidfraction()
-                        );
-                        pEqn.setReference(pRefCell, pRefValue);
-
-                        if
-                        (
-                            corr == nCorr-1
-                         && nonOrth == nNonOrthCorr
-                        )
-                        {
-                            pEqn.solve(mesh.solver("pFinal"));
-                        }
-                        else
-                        {
-                            pEqn.solve();
-                        }
-
-                        if (nonOrth == nNonOrthCorr)
-                        {
-                            phiGes -= pEqn.flux();
-                            phi = phiGes;
-                        }
-
-                    } // end non-orthogonal corrector loop
-
-                    #include "continuityErrorPhiPU.H"
-
-                    if (modelType=="B" || modelType=="Bfull")
-                        U -= rUA*fvc::grad(p) - Ksl/rho*Us*rUA;
-                    else
-                        U -= voidfraction*rUA*fvc::grad(p) - Ksl/rho*Us*rUA;
-
-                    U.correctBoundaryConditions();
-
-                } // end piso loop
+                    #include "pEqn.H"
+                }
             }
 
+            laminarTransport.correct();
             turbulence->correct();
-        }// end solveFlow
+        }
         else
         {
             Info << "skipping flow solution." << endl;
@@ -202,6 +128,9 @@ int main(int argc, char *argv[])
         Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
             << "  ClockTime = " << runTime.elapsedClockTime() << " s"
             << nl << endl;
+
+        particleCloud.clockM().stop("Flow");
+        particleCloud.clockM().stop("Global");
     }
 
     Info<< "End\n" << endl;

applications/solvers/cfdemSolverPisoScalar/createFields.H

@@ -1,36 +1,36 @@
-    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
-    );
-
-//========================
-// drag law modelling
-//========================
-
+    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
+    );
+
+//========================
+// drag law modelling
+//========================
+
     Info<< "\nReading momentum exchange field Ksl\n" << endl;
     volScalarField Ksl
     (
@@ -44,8 +44,8 @@
         ),
         mesh
         //dimensionedScalar("0", dimensionSet(0, 0, -1, 0, 0), 1.0)
-    );
-
+    );
+
     Info<< "\nReading voidfraction field voidfraction = (Vgas/Vparticle)\n" << endl;
     volScalarField voidfraction
     (
@@ -58,8 +58,8 @@
             IOobject::AUTO_WRITE
         ),
         mesh
-    );
-
+    );
+
     Info<< "\nCreating density field rho\n" << endl;
     volScalarField rho
     (
@@ -71,27 +71,27 @@
             IOobject::READ_IF_PRESENT,
             IOobject::AUTO_WRITE
         ),
-        mesh,
+        mesh,
         dimensionedScalar("0", dimensionSet(1, -3, 0, 0, 0), 1.0)
-    );
-
-    Info<< "Reading particle velocity field Us\n" << endl;
-    volVectorField Us
-    (
-        IOobject
-        (
-            "Us",
-            runTime.timeName(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::AUTO_WRITE
-        ),
-        mesh
-    );
-
-//========================
-// scalar field modelling
-//========================
+    );
+
+    Info<< "Reading particle velocity field Us\n" << endl;
+    volVectorField Us
+    (
+        IOobject
+        (
+            "Us",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+//========================
+// scalar field modelling
+//========================
     Info<< "\nCreating dummy density field rho = 1\n" << endl;
     volScalarField T
     (
@@ -103,10 +103,10 @@
             IOobject::MUST_READ,
             IOobject::AUTO_WRITE
         ),
-        mesh//,
+        mesh//,
         //dimensionedScalar("0", dimensionSet(0, 0, -1, 1, 0), 273.15)
-    );
-
+    );
+
     Info<< "\nCreating fluid-particle heat flux field\n" << endl;
     volScalarField Tsource
     (
@@ -118,57 +118,59 @@
             IOobject::MUST_READ,
             IOobject::AUTO_WRITE
         ),
-        mesh//,
+        mesh//,
         //dimensionedScalar("0", dimensionSet(0, 0, -1, 1, 0), 0.0)
-    );
-
-    IOdictionary transportProperties
-    (
-        IOobject
-        (
-            "transportProperties",
-            runTime.constant(),
-            mesh,
-            IOobject::MUST_READ,
-            IOobject::NO_WRITE
-        )
-    );
-
-    dimensionedScalar DT
-    (
-        transportProperties.lookup("DT")
-    );
-
-//========================
-
-//#   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*voidfraction) & 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)
-    );
+    );
+
+    IOdictionary transportProperties
+    (
+        IOobject
+        (
+            "transportProperties",
+            runTime.constant(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::NO_WRITE
+        )
+    );
+
+    dimensionedScalar DT
+    (
+        transportProperties.lookup("DT")
+    );
+
+//========================
+
+//#   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*voidfraction) & 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"

src/finiteVolume/Make/files

@@ -3,7 +3,6 @@ derivedFvPatchFields = $(fvPatchFields)/derived
 
 general = cfdTools/general
 
-$(derivedFvPatchFields)/uniformFixedValueVoidfraction/uniformFixedValueVoidfractionFvPatchFields.C
 $(derivedFvPatchFields)/uniformFixedValueTube/uniformFixedValueTubeFvPatchFields.C
 
 LIB = $(CFDEM_LIB_DIR)/libfiniteVolumeCFDEM

src/finiteVolume/fields/fvPatchFields/derived/uniformFixedValueTube/uniformFixedValueTubeFvPatchField.C

@@ -88,7 +88,7 @@ uniformFixedValueTubeFvPatchField<Type>::uniformFixedValueTubeFvPatchField
 )
 :
     fixedValueFvPatchField<Type>(p, iF),
-    uniformValue_(DataEntry<Type>::New("uniformValue", dict)),
+    uniformValue_(Function1<Type>::New("uniformValue", dict)),
 	pName_("p"), //JOKER
 	phiName_("phi"), //JOKER
     velocityFieldName_("U"),

src/finiteVolume/fields/fvPatchFields/derived/uniformFixedValueTube/uniformFixedValueTubeFvPatchField.H

@@ -40,7 +40,7 @@ SourceFiles
 
 #include "Random.H"
 #include "fixedValueFvPatchFields.H"
-#include "DataEntry.H"
+#include "Function1.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -58,7 +58,7 @@ class uniformFixedValueTubeFvPatchField
 {
     // Private data
 
-        autoPtr<DataEntry<Type> > uniformValue_;
+        autoPtr<Function1<Type> > uniformValue_;
 
 		word pName_; //JOKER pressure
 		

src/finiteVolume/fields/fvPatchFields/derived/uniformFixedValueVoidfraction/uniformFixedValueVoidfractionFvPatchField.C

@@ -73,7 +73,7 @@ uniformFixedValueVoidfractionFvPatchField<Type>::uniformFixedValueVoidfractionFv
 )
 :
     fixedValueFvPatchField<Type>(p, iF),
-    uniformValue_(DataEntry<Type>::New("uniformValue", dict)),
+    uniformValue_(Function1<Type>::New("uniformValue", dict)),
     voidfractionFieldName_("voidfraction")
 {
     const scalar t = this->db().time().timeOutputValue();

src/finiteVolume/fields/fvPatchFields/derived/uniformFixedValueVoidfraction/uniformFixedValueVoidfractionFvPatchField.H

@@ -40,7 +40,7 @@ SourceFiles
 
 #include "Random.H"
 #include "fixedValueFvPatchFields.H"
-#include "DataEntry.H"
+#include "Function1.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -58,7 +58,7 @@ class uniformFixedValueVoidfractionFvPatchField
 {
     // Private data
 
-        autoPtr<DataEntry<Type> > uniformValue_;
+        autoPtr<Function1<Type> > uniformValue_;
 
         word voidfractionFieldName_;
 

src/lagrangian/cfdemParticle/Make/files

@@ -5,6 +5,8 @@ energyModels = subModels/energyModel
 forceModels = subModels/forceModel
 forceSubModels = subModels/forceModel/forceSubModels
 forceModelsMS = subModels/forceModelMS
+thermCondModels = subModels/thermCondModel
+chemistryModels = subModels/chemistryModel
 IOModels = subModels/IOModel
 voidFractionModels = subModels/voidFractionModel
 locateModels = subModels/locateModel
@@ -27,9 +29,15 @@ derived/cfdemCloudEnergy/cfdemCloudEnergy.C
 $(cfdTools)/global.C
 $(cfdTools)/newGlobal.C
 
-$(energyModels)/heatTransferModel/heatTransferModel/heatTransferModel.C
-$(energyModels)/heatTransferModel/heatTransferModel/newHeatTransferModel.C
-$(energyModels)/heatTransferModel/noHeatTransfer/noHeatTransfer.C
+$(energyModels)/energyModel/energyModel.C
+$(energyModels)/energyModel/newEnergyModel.C
+$(energyModels)/heatTransferGunn/heatTransferGunn.C
+$(energyModels)/heatTransferGunnImplicit/heatTransferGunnImplicit.C
+
+$(thermCondModels)/thermCondModel/thermCondModel.C
+$(thermCondModels)/thermCondModel/newThermCondModel.C
+$(thermCondModels)/SyamlalThermCond/SyamlalThermCond.C
+$(thermCondModels)/noTherm/noThermCond.C
 
 $(forceModels)/forceModel/forceModel.C
 $(forceModels)/forceModel/newForceModel.C
@@ -46,7 +54,6 @@ $(forceModels)/ShirgaonkarIB/ShirgaonkarIB.C
 $(forceModels)/KochHillDrag/KochHillDrag.C
 $(forceModels)/KochHillRWDrag/KochHillRWDrag.C
 $(forceModels)/LaEuScalarTemp/LaEuScalarTemp.C
-$(forceModels)/LaEuScalarEnergy/LaEuScalarEnergy.C
 $(forceModels)/virtualMassForce/virtualMassForce.C
 $(forceModels)/gradPForce/gradPForce.C
 $(forceModels)/viscForce/viscForce.C
@@ -57,6 +64,7 @@ $(forceModels)/volWeightedAverage/volWeightedAverage.C
 $(forceModels)/deactivateForce/deactivateForce.C
 $(forceModels)/freeStreaming/freeStreaming.C
 $(forceModels)/isotropicFluctuations/isotropicFluctuations.C
+$(forceModels)/BeetstraDrag/BeetstraDrag.C
 
 $(forceModelsMS)/forceModelMS/forceModelMS.C
 $(forceModelsMS)/forceModelMS/newForceModelMS.C
@@ -68,6 +76,11 @@ $(forceSubModels)/ImEx/ImEx.C
 $(forceSubModels)/recF/recF.C
 $(forceSubModels)/recU/recU.C
 
+$(chemistryModels)/chemistryModel/chemistryModel.C
+$(chemistryModels)/chemistryModel/newChemistryModel.C
+$(chemistryModels)/species/species.C
+$(chemistryModels)/noChemistry/noChemistry.C
+
 $(probeModels)/probeModel/probeModel.C
 $(probeModels)/probeModel/newProbeModel.C
 $(probeModels)/noProbe/noProbe.C
@@ -119,13 +132,11 @@ $(dataExchangeModels)/oneWayVTK/oneWayVTK.C
 $(dataExchangeModels)/twoWayFiles/twoWayFiles.C
 $(dataExchangeModels)/noDataExchange/noDataExchange.C
 $(dataExchangeModels)/twoWayMPI/twoWayMPI.C
-$(dataExchangeModels)/twoWayMany2Many/twoWayMany2Many.C
 
 $(averagingModels)/averagingModel/averagingModel.C
 $(averagingModels)/averagingModel/newAveragingModel.C
 $(averagingModels)/dilute/dilute.C
 $(averagingModels)/dense/dense.C
-$(averagingModels)/noAveraging/noAveraging.C
 
 $(clockModels)/clockModel/clockModel.C
 $(clockModels)/clockModel/newClockModel.C

src/lagrangian/cfdemParticle/Make/options

@@ -12,7 +12,7 @@ EXE_INC = \
     -I ./cfdemParticle \
     -I$(LIB_SRC)/finiteVolume/lnInclude \
     -I$(LIB_SRC)/transportModels \
-    -I$(LIB_SRC)/turbulenceModels/incompressible/turbulenceModel \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
     -I$(LIB_SRC)/meshTools/lnInclude \
     -I$(LIB_SRC)/lagrangian/basic/lnInclude \
     -I$(LIB_SRC)/OpenFOAM/containers/HashTables/labelHashSet \
@@ -24,8 +24,7 @@ LIB_LIBS = \
      $(PLIBS) \
     -L$(CFDEM_LIB_DIR) \
     -lfiniteVolume \
-    -lincompressibleRASModels \
-    -lincompressibleLESModels \
+    -lincompressibleTurbulenceModels \
     -lmeshTools \
     -llagrangian \
     -lmpi_cxx \

src/lagrangian/cfdemParticle/cfdTools/continuityErrorPhiPU.H

@@ -34,7 +34,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 {
-    volScalarField contErr( fvc::div(phiGes) + fvc::ddt(voidfraction) );
+    volScalarField contErr( fvc::div(phi) + fvc::ddt(voidfraction) );
 
     scalar sumLocalContErr = runTime.deltaTValue()*
         mag(contErr)().weightedAverage(mesh.V()).value();

src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.C

@@ -125,15 +125,7 @@ cfdemCloud::cfdemCloud
     ),
     turbulence_
     (
-        #if defined(version21) || defined(version16ext)
-            #ifdef compre
-                mesh.lookupObject<compressible::turbulenceModel>
-            #else
-                mesh.lookupObject<incompressible::turbulenceModel>
-            #endif
-        #elif defined(version15)
-            mesh.lookupObject<incompressible::RASModel>
-        #endif
+        mesh.lookupObject<turbulenceModel>
         (
             turbulenceModelType_
         )
@@ -442,6 +434,11 @@ int cfdemCloud::nrForceModels()
     return forceModels_.size();
 }
 
+int cfdemCloud::nrMomCoupleModels()
+{
+    return momCoupleModels_.size();
+}
+
 scalar cfdemCloud::voidfraction(int index)
 {
     return voidfractions()[index][0];

src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.H

@@ -49,11 +49,7 @@ SourceFiles
 #include "fvCFD.H"
 #include "IFstream.H"
 
-#if defined(version21) || defined(version16ext)
-    #include <turbulenceModel.H>
-#elif defined(version15)
-    #include <RASModel.H>
-#endif
+#include <turbulenceModel.H>
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -157,15 +153,7 @@ protected:
 
     mutable volScalarField ddtVoidfraction_;
 
-    #if defined(version21) || defined(version16ext)
-        #ifdef compre
-            const compressible::turbulenceModel& turbulence_;
-        #else
-            const incompressible::turbulenceModel& turbulence_;
-        #endif
-    #elif defined(version15)
-        const incompressible::RASModel& turbulence_;
-    #endif
+    const turbulenceModel& turbulence_;
 
     autoPtr<forceModel>* forceModel_;
 
@@ -244,6 +232,8 @@ public:
        virtual const forceModel& forceM(int);
 
        virtual int nrForceModels();
+       
+       virtual int nrMomCoupleModels();
 
        scalar voidfraction(int);
 
@@ -364,20 +354,14 @@ public:
 
        inline autoPtr<liggghtsCommandModel>* liggghtsCommand() const;
 
-       #if defined(version21) || defined(version16ext)
-            #ifdef compre
-                inline const compressible::turbulenceModel& turbulence() const;
-            #else
-                inline const incompressible::turbulenceModel& turbulence() const;
-            #endif
-       #elif defined(version15)
-           inline const incompressible::RASModel& turbulence() const;
-       #endif
-
+       inline const turbulenceModel& turbulence() const;
+       
     // Write
 
       // write cfdemCloud internal data
         virtual bool evolve(volScalarField&,volVectorField&,volVectorField&);
+        
+        virtual void postFlow() {}
 
         virtual bool reAllocArrays();
 

src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloudI.H

@@ -302,15 +302,7 @@ inline autoPtr<liggghtsCommandModel>* cfdemCloud::liggghtsCommand() const
     return liggghtsCommand_;
 }
 
-#if defined(version21) || defined(version16ext)
-    #ifdef compre
-        inline const compressible::turbulenceModel& cfdemCloud::turbulence() const
-    #else
-        inline const incompressible::turbulenceModel& cfdemCloud::turbulence() const
-    #endif
-#elif defined(version15)
-    inline const incompressible::RASModel& cfdemCloud::turbulence() const
-#endif
+inline const turbulenceModel& cfdemCloud::turbulence() const
 {
     return turbulence_;
 }

src/lagrangian/cfdemParticle/derived/cfdemCloudEnergy/cfdemCloudEnergy.C

@@ -1,40 +1,27 @@
 /*---------------------------------------------------------------------------*\
-    CFDEMcoupling - Open Source CFD-DEM coupling
-
-    CFDEMcoupling is part of the CFDEMproject
-    www.cfdem.com
-                                Thomas Lichtenegger, thomas.lichtenegger@jku.at
-                                Copyright 2009-2012 JKU Linz
-                                Copyright 2012-2015 DCS Computing GmbH, Linz
-                                Copyright 2015- JKU Linz
--------------------------------------------------------------------------------
 License
-    This file is part of CFDEMcoupling.
 
-    CFDEMcoupling is free software; you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by the
-    Free Software Foundation; either version 3 of the License, or (at your
-    option) any later version.
+    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.
 
-    CFDEMcoupling is distributed in the hope that it will be useful, but WITHOUT
+    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 CFDEMcoupling; if not, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+    along with this code.  If not, see <http://www.gnu.org/licenses/>.
+
+    Copyright (C) 2015- Thomas Lichtenegger, JKU Linz, Austria
 
-Description
-    This code is designed to realize coupled CFD-DEM simulations using LIGGGHTS
-    and OpenFOAM(R). Note: this code is not part of OpenFOAM(R) (see DISCLAIMER).
 \*---------------------------------------------------------------------------*/
 
-#include "fileName.H"
 #include "cfdemCloudEnergy.H"
-#include "heatTransferModel.H"
-#include <mpi.h>
-#include "IOmanip.H"
+#include "energyModel.H"
+#include "thermCondModel.H"
+#include "chemistryModel.H"
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
@@ -52,16 +39,35 @@ cfdemCloudEnergy::cfdemCloudEnergy
 )
 :
     cfdemCloud(mesh),
-    heatTransferModel_
+    energyModels_(couplingProperties_.lookup("energyModels")),
+    implicitEnergyModel_(false),
+    thermCondModel_
     (
-        heatTransferModel::New
+        thermCondModel::New
+        (
+            couplingProperties_,
+            *this
+        )
+    ),
+    chemistryModel_
+    (
+        chemistryModel::New
         (
             couplingProperties_,
             *this
         )
     )
 {
-
+    energyModel_ = new autoPtr<energyModel>[nrEnergyModels()];
+    for (int i=0;i<nrEnergyModels();i++)
+    {
+        energyModel_[i] = energyModel::New
+        (
+            couplingProperties_,
+            *this,
+            energyModels_[i]
+        );
+    }
 }
 
 
@@ -72,11 +78,95 @@ cfdemCloudEnergy::~cfdemCloudEnergy()
 
 }
 
+// * * * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * * * //
+
+void cfdemCloudEnergy::calcEnergyContributions()
+{
+    for (int i=0;i<nrEnergyModels();i++)
+        energyModel_[i]().calcEnergyContribution();
+}
+
+void cfdemCloudEnergy::speciesExecute()
+{
+        chemistryModel_().execute();
+}
 
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
-inline const heatTransferModel& cfdemCloudEnergy::heatTransferM() const
+int cfdemCloudEnergy::nrEnergyModels()
 {
-    return heatTransferModel_;  
+    return energyModels_.size();
+}
+
+bool& cfdemCloudEnergy::implicitEnergyModel()
+{
+    return implicitEnergyModel_;
+}
+
+const energyModel& cfdemCloudEnergy::energyM(int i)
+{
+    return energyModel_[i];
+}
+
+const chemistryModel& cfdemCloudEnergy::chemistryM()
+{
+    return chemistryModel_;
+}
+
+const thermCondModel& cfdemCloudEnergy::thermCondM()
+{
+    return thermCondModel_;
+}
+
+void cfdemCloudEnergy::energyContributions(volScalarField& Qsource)
+{
+    Qsource.primitiveFieldRef()=0.0;
+    Qsource.boundaryFieldRef()=0.0;
+    for (int i=0;i<nrEnergyModels();i++)
+        energyM(i).addEnergyContribution(Qsource);
+}
+
+void cfdemCloudEnergy::energyCoefficients(volScalarField& Qcoeff)
+{
+    Qcoeff.primitiveFieldRef()=0.0;
+    Qcoeff.boundaryFieldRef()=0.0;
+    for (int i=0;i<nrEnergyModels();i++)
+        energyM(i).addEnergyCoefficient(Qcoeff);
+}
+
+bool cfdemCloudEnergy::evolve
+(
+    volScalarField& alpha,
+    volVectorField& Us,
+    volVectorField& U
+)
+{
+    if (cfdemCloud::evolve(alpha, Us, U))
+    {
+        // calc energy contributions
+        // position 26 was already defined as Flow in clockModels and RhoPimpleChem solver.
+        clockM().start(27,"calcEnergyContributions");
+        if(verbose_) Info << "- calcEnergyContributions" << endl;
+        calcEnergyContributions();
+        if(verbose_) Info << "calcEnergyContributions done." << endl;
+        clockM().stop("calcEnergyContributions");
+
+        // execute chemical model species
+        clockM().start(28,"speciesExecute");
+        if(verbose_) Info << "- speciesExecute()" << endl;
+        speciesExecute();
+        if(verbose_) Info << "speciesExecute done" << endl;
+        clockM().stop("speciesExecute");
+
+	return true;
+    }
+    return false;
+}
+
+void cfdemCloudEnergy::postFlow()
+{
+    cfdemCloud::postFlow();
+    for (int i=0;i<nrEnergyModels();i++)
+        energyModel_[i]().postFlow();
 }
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

src/lagrangian/cfdemParticle/derived/cfdemCloudEnergy/cfdemCloudEnergy.H

@@ -1,35 +1,25 @@
 /*---------------------------------------------------------------------------*\
-    CFDEMcoupling - Open Source CFD-DEM coupling
-
-    CFDEMcoupling is part of the CFDEMproject
-    www.cfdem.com
-                                Thomas Lichtenegger, thomas.lichtenegger@jku.at
-                                Copyright 2009-2012 JKU Linz
-                                Copyright 2012-2015 DCS Computing GmbH, Linz
-                                Copyright 2015- JKU Linz
--------------------------------------------------------------------------------
 License
-    This file is part of CFDEMcoupling.
 
-    CFDEMcoupling is free software; you can redistribute it and/or modify it
-    under the terms of the GNU General Public License as published by the
-    Free Software Foundation; either version 3 of the License, or (at your
-    option) any later version.
+    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.
 
-    CFDEMcoupling is distributed in the hope that it will be useful, but WITHOUT
+    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 CFDEMcoupling; if not, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+    along with this code.  If not, see <http://www.gnu.org/licenses/>.
+
+    Copyright (C) 2015- Thomas Lichtenegger, JKU Linz, Austria
+
 
 Description
-    This code is designed to realize coupled CFD-DEM simulations using LIGGGHTS
-    and OpenFOAM(R). Note: this code is not part of OpenFOAM(R) (see DISCLAIMER).
-
-    cloud class managing DEM data for CFD-DEM coupling of compressible flows
+    cloud class managing DEM data for CFD-DEM coupling with energy models,
+    e.g. for compressible, reacting flows
 
 Class
     Foam::cfdemCloudEnergy derived from cfdemCloud
@@ -49,7 +39,9 @@ SourceFiles
 namespace Foam
 {
 // forward declarations
-class heatTransferModel;
+class energyModel;
+class thermCondModel;
+class chemistryModel;
 /*---------------------------------------------------------------------------*\
                            Class cfdemCloudEnergy Declaration
 \*---------------------------------------------------------------------------*/
@@ -60,10 +52,26 @@ class cfdemCloudEnergy
 {
 protected:
 
-     autoPtr<heatTransferModel> heatTransferModel_;
+    const wordList energyModels_;
+    
+    bool implicitEnergyModel_;
+    
+    const wordList chemistryModels_;
 
+    autoPtr<energyModel>* energyModel_;
+    
+    autoPtr<thermCondModel> thermCondModel_;
+    
+    autoPtr<chemistryModel> chemistryModel_;
+    
+    void calcEnergyContributions();
+
+    void speciesExecute();
+    
 public:
 
+    friend class energyModel;
+
     // Constructors
 
         //- Construct from components
@@ -80,8 +88,25 @@ public:
     // Member Functions
 
 
-    // Access
-        inline const heatTransferModel& heatTransferM() const;
+       const energyModel& energyM(int);
+       
+       const thermCondModel& thermCondM();
+       
+       const chemistryModel& chemistryM();
+
+       int nrEnergyModels();
+
+       inline const wordList& energyModels();
+       
+       bool& implicitEnergyModel();
+       
+       void energyContributions(volScalarField&);
+       
+       void energyCoefficients(volScalarField&);
+       
+       bool evolve(volScalarField&,volVectorField&,volVectorField&);
+       
+       void postFlow();
 
 };
 
@@ -92,6 +117,8 @@ public:
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
+#include "cfdemCloudEnergyI.H"
+
 #endif
 
 // ************************************************************************* //

src/lagrangian/cfdemParticle/derived/cfdemCloudEnergy/cfdemCloudEnergyI.H

@@ -0,0 +1,38 @@
+/*---------------------------------------------------------------------------*\
+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
+
+\*---------------------------------------------------------------------------*/
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+inline const wordList& cfdemCloudEnergy::energyModels()
+{
+    return energyModels_;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/lagrangian/cfdemParticle/etc/OFversion/OFversion.H

@@ -1,4 +1,5 @@
-#define version23  // currently being used
+#define version4x
+//#define version23  // currently being used
 //#define version22
 //#define version21
 //#define version16ext

src/lagrangian/cfdemParticle/etc/functions.sh

@@ -76,17 +76,17 @@ compileLib()
         i=$((${#str}-4))
         ending=${str:$i:4}
         if [[ $ending == "Comp" ]]; then
-                echo "Compiling a compressible library - so doing an rmdepall of incomp library first."
+                echo "Compiling a compressible library - so doing an wrmdep -a of incomp library first."
                 echo "Please make sure to have the compressible libraries first in the library-list.txt!"
                 cd $CFDEM_SRC_DIR/lagrangian/cfdemParticle
                 echo "changing to $PWD"
-                rmdepall 2>&1 | tee -a $logpath/$logfileName
+                wrmdep -a 2>&1 | tee -a $logpath/$logfileName
                 cd $casePath
                 echo "changing to $PWD"
             else
                 echo "Compiling a incompressible library."
         fi
-        rmdepall 2>&1 | tee -a $logpath/$logfileName
+        wrmdep -a 2>&1 | tee -a $logpath/$logfileName
         wclean 2>&1 | tee -a $logpath/$logfileName
     #fi
     wmake libso 2>&1 | tee -a $logpath/$logfileName
@@ -128,7 +128,7 @@ compileSolver()
 
     #- wclean and wmake
     #if [ $doClean != "noClean" ]; then
-        rmdepall 2>&1 | tee -a $logpath/$logfileName
+        wrmdep -a 2>&1 | tee -a $logpath/$logfileName
         wclean 2>&1 | tee -a $logpath/$logfileName
     #fi
     
@@ -284,7 +284,7 @@ cleanCFDEM()
 
             cd  $path
             echo "cleaning library $PWD"
-            rmdepall
+            wrmdep -a
             wclean    
             rm -r ./Make/linux*
             rm -r ./lnInclude
@@ -338,7 +338,7 @@ cleanCFDEM()
 
             cd  $path            
             echo "cleaning solver $PWD"
-            rmdepall
+            wrmdep -a
             wclean    
     done
 }

src/lagrangian/cfdemParticle/etc/solver-list.txt

@@ -3,5 +3,4 @@ tfmSolverRecurrence/dir
 cfdemSolverRhoPimple/dir
 cfdemSolverPisoMS/dir
 cfdemSolverPiso/dir
-cfdemSolverIB/dir
 cfdemSolverPisoScalar/dir

src/lagrangian/cfdemParticle/subModels/averagingModel/averagingModel/averagingModel.C

@@ -59,7 +59,7 @@ void averagingModel::undoVectorAverage
 {
 // WARNING - not sure if this is valid for dilute model!!!
 
-    if(!single) fieldPrev.internalField() = fieldNext.internalField();
+    if(!single) fieldPrev.ref() = fieldNext.ref();
 
     label cellI;
     vector valueVec;
@@ -135,7 +135,7 @@ void averagingModel::setVectorSum
 (
     volVectorField& field,
     double**& value,
-    double**& weight,
+    double**const& weight,
     double**const& mask
 ) const
 {
@@ -169,7 +169,7 @@ void averagingModel::setVectorSumSimple
 (
     volVectorField& field,
     double**& value,
-    double**& weight,
+    double**const& weight,
     int nP
 ) const
 {
@@ -232,7 +232,7 @@ void averagingModel::setScalarSum
 void averagingModel::setDSauter
 (
     volScalarField& dSauter,
-    double**& weight,
+    double**const& weight,
     volScalarField& weightField,
     label myParticleType
 ) const
@@ -303,13 +303,13 @@ void averagingModel::setDSauter
 
 void averagingModel::resetVectorAverage(volVectorField& prev,volVectorField& next,bool single) const
 {
-    if(!single) prev.internalField() = next.internalField();
-    next.internalField() = vector::zero;
+    if(!single) prev.ref() = next.ref();
+    next.primitiveFieldRef() = vector::zero;
 }
 
 void averagingModel::resetWeightFields() const
 {
-    UsWeightField_.internalField() = 0;
+    UsWeightField_.ref() = 0;
 }
 
 
@@ -352,12 +352,12 @@ tmp<volVectorField> Foam::averagingModel::UsInterp() const
 
     if (particleCloud_.dataExchangeM().couplingStep() > 1)
     {
-        tsource() = (1 - particleCloud_.dataExchangeM().timeStepFraction()) * UsPrev_
+        tsource.ref() = (1 - particleCloud_.dataExchangeM().timeStepFraction()) * UsPrev_
                     + particleCloud_.dataExchangeM().timeStepFraction() * UsNext_;
     }
     else
     {
-        tsource() = UsNext_;
+        tsource.ref() = UsNext_;
     }
 
     return tsource;
@@ -415,55 +415,6 @@ averagingModel::averagingModel
     )
 {}
 
-// Construct from components, specify initial value for fields
-averagingModel::averagingModel
-(
-    const dictionary& dict,
-    cfdemCloud& sm,
-    const vector u0
-)
-:
-    dict_(dict),
-    particleCloud_(sm),
-    UsWeightField_
-    (
-        IOobject
-        (
-            "UsWeightField_",
-            particleCloud_.mesh().time().timeName(),
-            particleCloud_.mesh(),
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        particleCloud_.mesh(),
-        dimensionedScalar("zero", dimensionSet(0,0,0,0,0), 0.0)
-    ),
-    UsPrev_
-    (   IOobject
-        (
-            "UsPrev",
-            sm.mesh().time().timeName(),
-            sm.mesh(),
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        sm.mesh(),
-        dimensionedVector("zero", dimensionSet(0,1,-1,0,0),u0)
-    ),
-    UsNext_
-    (   IOobject
-        (
-            "UsNext",
-            sm.mesh().time().timeName(),
-            sm.mesh(),
-            IOobject::NO_READ,
-            IOobject::NO_WRITE
-        ),
-        sm.mesh(),
-        dimensionedVector("zero", dimensionSet(0,1,-1,0,0),u0)
-    )
-{}
-
 
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 

src/lagrangian/cfdemParticle/subModels/averagingModel/averagingModel/averagingModel.H

@@ -98,14 +98,6 @@ public:
             cfdemCloud& sm
         );
 
-        //- Construct from components, specifiy initial value of fields
-        averagingModel
-        (
-            const dictionary& dict,
-            cfdemCloud& sm,
-            const vector u0
-        );
-
 
     // Destructor
 
@@ -126,7 +118,7 @@ public:
         (
             volScalarField& field,
             double**& value,
-            double**& weight,
+            double**const& weight,
             volScalarField& weightField,
             double**const& mask
         ) const = 0;
@@ -135,7 +127,7 @@ public:
         (
             volVectorField& field,
             double**& value,
-            double**& weight,
+            double**const& weight,
             volScalarField& weightField,
             double**const& mask
         ) const = 0;
@@ -164,7 +156,7 @@ public:
         (
             volVectorField& field,
             double**& value,
-            double**& weight,
+            double**const& weight,
             double**const& mask
         ) const;
 
@@ -172,7 +164,7 @@ public:
         (
             volVectorField& field,
             double**& value,
-            double**& weight,
+            double**const& weight,
             int np
         ) const;
 
@@ -187,7 +179,7 @@ public:
         void setDSauter
         (
             volScalarField& dSauter,
-            double**& weight,
+            double**const& weight,
             volScalarField& weightField,
             label myParticleType = 0  // can be evaluated for a special type
         ) const;

src/lagrangian/cfdemParticle/subModels/averagingModel/dense/dense.C

@@ -78,7 +78,7 @@ void dense::setScalarAverage
 (
     volScalarField& field,
     double**& value,
-    double**& weight,
+    double**const& weight,
     volScalarField& weightField,
     double**const& mask
 ) const
@@ -127,7 +127,7 @@ void dense::setVectorAverage
 (
     volVectorField& field,
     double**& value,
-    double**& weight,
+    double**const& weight,
     volScalarField& weightField,
     double**const& mask
 ) const

src/lagrangian/cfdemParticle/subModels/averagingModel/dense/dense.H

@@ -82,7 +82,7 @@ public:
         (
             volScalarField& field,
             double**& value,
-            double**& weight,
+            double**const& weight,
             volScalarField& weightField,
             double**const& mask
         ) const;
@@ -91,7 +91,7 @@ public:
         (
             volVectorField& field,
             double**& value,
-            double**& weight,
+            double**const& weight,
             volScalarField& weightField,
             double**const& mask
         ) const;

src/lagrangian/cfdemParticle/subModels/averagingModel/dilute/dilute.C

@@ -77,7 +77,7 @@ void dilute::setScalarAverage
 (
     volScalarField& field,
     double**& value,
-    double**& weight,
+    double**const& weight,
     volScalarField& weightField,
     double**const& mask
 ) const
@@ -115,7 +115,7 @@ void dilute::setVectorAverage
 (
     volVectorField& field,
     double**& value,
-    double**& weight,
+    double**const& weight,
     volScalarField& weightField,
     double**const& mask
 ) const

src/lagrangian/cfdemParticle/subModels/averagingModel/dilute/dilute.H

@@ -84,7 +84,7 @@ public:
         (
             volScalarField& field,
             double**& value,
-            double**& weight,
+            double**const& weight,
             volScalarField& weightField,
             double**const& mask
         ) const;
@@ -93,7 +93,7 @@ public:
         (
             volVectorField& field,
             double**& value,
-            double**& weight,
+            double**const& weight,
             volScalarField& weightField,
             double**const& mask
         ) const;

src/lagrangian/cfdemParticle/subModels/averagingModel/noAveraging/noAveraging.C

@@ -55,7 +55,7 @@ noAveraging::noAveraging
     cfdemCloud& sm
 )
 :
-    averagingModel(dict,sm,vector(0,0,0))
+    averagingModel(dict,sm)
 {}
 
 

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

@@ -1,37 +1,25 @@
 /*---------------------------------------------------------------------------*\
-    CFDEMcoupling - Open Source CFD-DEM coupling
-
-    CFDEMcoupling is part of the CFDEMproject
-    www.cfdem.com
-                                Christoph Goniva, christoph.goniva@cfdem.com
-                                Copyright 2009-2012 JKU Linz
-                                Copyright 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.
+    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.
 
-    CFDEMcoupling is distributed in the hope that it will be useful, but WITHOUT
+    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 CFDEMcoupling; if not, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+    along with this code.  If not, see <http://www.gnu.org/licenses/>.
+
+    Copyright (C) 2015- Thomas Lichtenegger, JKU Linz, Austria
 
-Description
-    This code is designed to realize coupled CFD-DEM simulations using LIGGGHTS
-    and OpenFOAM(R). Note: this code is not part of OpenFOAM(R) (see DISCLAIMER).
 \*---------------------------------------------------------------------------*/
 
 #include "error.H"
-#include "heatTransferModel.H"
-#include "mathExtra.H"
+#include "chemistryModel.H"
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -39,16 +27,16 @@ namespace Foam
 
 // * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
 
-defineTypeNameAndDebug(heatTransferModel, 0);
+defineTypeNameAndDebug(chemistryModel, 0);
 
-defineRunTimeSelectionTable(heatTransferModel, dictionary);
+defineRunTimeSelectionTable(chemistryModel, dictionary);
 
 // * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 // * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
 
 // Construct from components
-heatTransferModel::heatTransferModel
+chemistryModel::chemistryModel
 (
     const dictionary& dict,
     cfdemCloudEnergy& sm
@@ -58,16 +46,14 @@ heatTransferModel::heatTransferModel
     particleCloud_(sm)
 {}
 
-
 // * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
 
-heatTransferModel::~heatTransferModel()
+chemistryModel::~chemistryModel()
 {}
 
 // * * * * * * * * * * * * * * * * Member Fct  * * * * * * * * * * * * * * * //
 
 
-
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

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

@@ -1,57 +1,40 @@
 /*---------------------------------------------------------------------------*\
-    CFDEMcoupling - Open Source CFD-DEM coupling
-
-    CFDEMcoupling is part of the CFDEMproject
-    www.cfdem.com
-                                Christoph Goniva, christoph.goniva@cfdem.com
-                                Copyright 2009-2012 JKU Linz
-                                Copyright 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.
+    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.
 
-    CFDEMcoupling is distributed in the hope that it will be useful, but WITHOUT
+    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 CFDEMcoupling; if not, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+    along with this code.  If not, see <http://www.gnu.org/licenses/>.
 
-Description
-    This code is designed to realize coupled CFD-DEM simulations using LIGGGHTS
-    and OpenFOAM(R). Note: this code is not part of OpenFOAM(R) (see DISCLAIMER).
-
-Class
-    heatTransferModel
-
-SourceFiles
-    heatTransferModel.C
+    Copyright (C) 2015- Thomas Lichtenegger, JKU Linz, Austria
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef heatTransferModel_H
-#define heatTransferModel_H
+#ifndef chemistryModel_H
+#define chemistryModel_H
 
 #include "fvCFD.H"
 #include "cfdemCloudEnergy.H"
 #include "probeModel.H"
+#include "interpolationCellPoint.H"
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class heatTransferModel Declaration
+                           Class chemistryModel Declaration
 \*---------------------------------------------------------------------------*/
 
-class heatTransferModel
+class chemistryModel
 {
 
 protected:
@@ -60,19 +43,20 @@ protected:
         const dictionary& dict_;
 
         cfdemCloudEnergy& particleCloud_;
+	
 
        
 public:
 
     //- Runtime type information
-    TypeName("heatTransferModel");
+    TypeName("chemistryModel");
 
     // Declare runtime constructor selection table
 
         declareRunTimeSelectionTable
         (
             autoPtr,
-            heatTransferModel,
+            chemistryModel,
             dictionary,
             (
                 const dictionary& dict,
@@ -85,21 +69,21 @@ public:
     // Constructors
 
         //- Construct from components
-        heatTransferModel
+        chemistryModel
         (
             const dictionary& dict,
-            cfdemCloudEnergy& sm
+            cfdemCloudEnergy& sm            
         );
 
 
     // Destructor
 
-        virtual ~heatTransferModel();
+        virtual ~chemistryModel();
 
 
     // Selector
 
-        static autoPtr<heatTransferModel> New
+        static autoPtr<chemistryModel> New
         (
             const dictionary& dict,
             cfdemCloudEnergy& sm
@@ -108,9 +92,11 @@ public:
 
     // Member Functions
 
-
-        virtual void partFluidHeatTransfer(volScalarField&) const = 0;
-	virtual scalar partThermCond() const = 0;
+       virtual void execute() = 0;
+       
+//        virtual tmp<fvScalarMatrix> Smi (const label i) const = 0;
+//        
+//        virtual tmp<fvScalarMatrix> Sm () const = 0;
 
 
 };

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

@@ -0,0 +1,70 @@
+/*---------------------------------------------------------------------------*\
+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 "chemistryModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+autoPtr<chemistryModel> chemistryModel::New
+(
+    const dictionary& dict,
+    cfdemCloudEnergy& sm
+)
+{
+    word chemistryModelType
+    (
+        dict.lookup("chemistryModel")
+    );
+    Info<< "Selecting chemistryModel "
+         << chemistryModelType << endl;
+
+    dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(chemistryModelType);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalError
+            << "chemistryModel::New(const dictionary&, cfdemCloudEnergy&) : "
+            << endl
+            << "    unknown chemistryModelType type "
+            << chemistryModelType
+            << ", constructor not in hash table" << endl << endl
+            << "    Valid chemistryModel types are :"
+            << endl;
+        Info<< dictionaryConstructorTablePtr_->toc()
+            << abort(FatalError);
+    }
+
+    return autoPtr<chemistryModel>(cstrIter()(dict,sm));
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

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

@@ -0,0 +1,75 @@
+/*---------------------------------------------------------------------------*\
+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 "noChemistry.H"
+#include "addToRunTimeSelectionTable.H"
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(noChemistry, 0);
+
+addToRunTimeSelectionTable(chemistryModel, noChemistry, dictionary);
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from components
+noChemistry::noChemistry
+(
+    const dictionary& dict,
+    cfdemCloudEnergy& sm
+)
+:
+    chemistryModel(dict,sm)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+noChemistry::~noChemistry()
+{}
+
+// * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
+
+// * * * * * * * * * * * * * * * * Member Fct  * * * * * * * * * * * * * * * //
+
+void noChemistry::execute()
+{}
+
+//tmp<Foam::fvScalarMatrix> noChemistry::Smi(const label i) const
+//{
+//    return tmp<fvScalarMatrix>(new fvScalarMatrix(0, dimMass/dimTime)); 
+//}
+
+//tmp<Foam::fvScalarMatrix> noChemistry::Sm() const
+//{
+//    return tmp<fvScalarMatrix>(new fvScalarMatrix(0, dimMass/dimTime)); 
+//}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

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

@@ -0,0 +1,87 @@
+/*---------------------------------------------------------------------------*\
+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
+    Chemistry turned off
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef noChemistry_H
+#define noChemistry_H
+
+#include "fvCFD.H"
+#include "cfdemCloudEnergy.H"
+#include "chemistryModel.H"
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class noChemistry Declaration
+\*---------------------------------------------------------------------------*/
+
+
+class noChemistry
+:
+    public chemistryModel
+{
+private:
+
+       
+public:
+
+    //- Runtime type information
+    TypeName("off");
+
+    // Constructors
+
+        //- Construct from components
+        noChemistry
+        (
+            const dictionary& dict,
+            cfdemCloudEnergy& sm
+        );
+
+
+    // Destructor
+
+        virtual ~noChemistry();
+
+
+    // Member Functions
+
+    void execute();
+    
+//     tmp<fvScalarMatrix> Smi(const label i) const;
+//     
+//     tmp<fvScalarMatrix> Sm() const;
+
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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

@@ -0,0 +1,289 @@
+/*---------------------------------------------------------------------------*\
+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
+                        M.Efe Kinaci, JKU Linz, Austria
+
+\*---------------------------------------------------------------------------*/
+
+#include "error.H"
+#include "species.H"
+#include "addToRunTimeSelectionTable.H"
+
+#include "dataExchangeModel.H"
+#include "IFstream.H"
+
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(species, 0);
+
+addToRunTimeSelectionTable
+(
+        chemistryModel,
+        species,
+        dictionary
+);
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from components
+species::species
+(
+    const dictionary& dict,
+    cfdemCloudEnergy& sm
+)
+:
+    chemistryModel(dict,sm),
+    propsDict_(dict.subDict(typeName + "Props")),
+    interpolation_(propsDict_.lookupOrDefault<bool>("interpolation",false)),
+    mesh_(sm.mesh()),
+    // define a file name in the coupling properties that contains the species
+    specDict_
+    (
+        IFstream
+        (
+            fileName(propsDict_.lookup("ChemistryFile")).expand()
+        )()
+    ),
+    // create a list from the Species table in the specified species dictionary
+    speciesNames_(specDict_.lookup("species")),
+    mod_spec_names_(speciesNames_.size()),
+    Y_(speciesNames_.size()),                           //volumeScalarFields created in the ts folders
+    concentrations_(speciesNames_.size(),NULL),         //the value of species concentration for every species
+    changeOfSpeciesMass_(speciesNames_.size(),NULL),    //the values that are received from DEM with the name of Modified_+species name
+    changeOfSpeciesMassFields_(speciesNames_.size()),   //the scalar fields generated with the values from Modified_+species names
+    changeOfGasMassField_                               //the total change of Gas Mass field (when the Modified species
+
+    (
+        IOobject
+        (
+            "changeOfGasMassField",
+            mesh_.time().timeName(),
+            mesh_,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+         ),
+        mesh_,
+        dimensionedScalar("zero",dimMass/(dimVol*dimTime),0.0)
+    ),
+    tempFieldName_(propsDict_.lookupOrDefault<word>("tempFieldName","T")),
+    tempField_(sm.mesh().lookupObject<volScalarField> (tempFieldName_)),
+    partTempName_(propsDict_.lookupOrDefault<word>("partTempName","partTemp")),
+    partTemp_(NULL),
+    densityFieldName_(propsDict_.lookupOrDefault<word>("densityFieldName","rho")),
+    rho_(sm.mesh().lookupObject<volScalarField> (densityFieldName_)),
+    partRhoName_(propsDict_.lookupOrDefault<word>("partRhoName","partRho")),
+    partRho_(NULL)
+  // voidfraction and velocity fields can be included by wish
+  /*  voidfractionFieldName_(propsDict_.lookup("voidfractionFieldName")),
+      voidfraction_(sm.mesh().lookupObject<volScalarField> (voidfractionFieldName_)),
+      velFieldName_(propsDict_.lookup("velFieldName")),
+      U_(sm.mesh().lookup<volVectorField> (velFieldName_)),*/
+
+{
+    Info << " Read species list from: " << specDict_.name() << endl;
+    Info << " Reading species list: " << speciesNames_ << endl;
+
+    for (int i=0; i<speciesNames_.size(); i++)
+    {
+        // Defining the Species volume scalar fields
+        Info << " Looking up species fields " << speciesNames_[i] << endl;
+        volScalarField& Y = const_cast<volScalarField&>
+                (sm.mesh().lookupObject<volScalarField>(speciesNames_[i]));
+        Y_.set(i, &Y);
+
+        // define the modified species names
+        mod_spec_names_[i] = "Modified_" + speciesNames_[i];
+
+        // Check if mod species are correct
+        Info << "Modified species names are: " << mod_spec_names_[i] << endl;
+
+        // Create new volScalarFields for the changed values of the species mass fields
+        changeOfSpeciesMassFields_.set
+        (
+            i,
+            new volScalarField
+            (
+                IOobject
+                (
+                "ModSpeciesMassField_"+Y_[i].name(),
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::AUTO_WRITE
+                ),
+                mesh_,
+                dimensionedScalar("0",dimMass/(dimVol*dimTime), 0)
+            )
+         );
+
+        Info << "what are the concentration names (Y_i): " << Y_[i].name() << endl;
+
+    }
+    allocateMyArrays();
+}
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+species::~species()
+{
+    delete partTemp_;
+    delete partRho_;
+
+    for (int i=0; i<speciesNames_.size();i++) delete [] concentrations_[i];
+    for (int i=0; i<speciesNames_.size();i++) delete [] changeOfSpeciesMass_[i];
+
+}
+
+// * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
+void species::allocateMyArrays() const
+{
+  // could be implemented similarly as forcemodel LaEuScalarTemp
+
+    // get memory for 2d arrays
+    double initVal=0.0;
+    particleCloud_.dataExchangeM().allocateArray(partRho_,initVal,1);
+    particleCloud_.dataExchangeM().allocateArray(partTemp_,initVal,1);
+
+    for (int i=0; i<speciesNames_.size(); i++)
+    {
+        particleCloud_.dataExchangeM().allocateArray(concentrations_[i],initVal,1);
+        particleCloud_.dataExchangeM().allocateArray(changeOfSpeciesMass_[i],initVal,1);
+    }
+}
+
+// * * * * * * * * * * * * * * * * Member Fct  * * * * * * * * * * * * * * * //
+
+void species::execute()
+{
+  // realloc the arrays
+    allocateMyArrays();
+
+  // get Y_i, T, rho at particle positions, fill arrays with them and push to LIGGGHTS
+
+    label  cellI=0;
+    scalar Tfluid(0);
+    scalar rhofluid(0);
+    List<scalar> Yfluid_;
+    Yfluid_.setSize(speciesNames_.size());
+
+    // defining interpolators for T, rho
+    interpolationCellPoint <scalar> TInterpolator_(tempField_);
+    interpolationCellPoint <scalar> rhoInterpolator_(rho_);
+
+    for (int index=0; index<particleCloud_.numberOfParticles(); index++)
+    {
+        cellI=particleCloud_.cellIDs()[index][0];
+        if (cellI >=0)
+        {
+            if(interpolation_)
+            {
+                vector position     =   particleCloud_.position(index);
+                Tfluid              =   TInterpolator_.interpolate(position,cellI);
+                rhofluid            =   rhoInterpolator_.interpolate(position,cellI);
+            }
+            else
+            {
+                Tfluid = tempField_[cellI];
+                rhofluid=rho_[cellI];
+            }
+
+            //fill arrays
+            partTemp_[index][0]=Tfluid;
+            partRho_[index][0]=rhofluid;
+            for (int i=0; i<speciesNames_.size();i++)
+            {
+	        Yfluid_[i] = Y_[i][cellI];
+            concentrations_[i][index][0]=Yfluid_[i];
+            }
+        }
+
+        if(particleCloud_.verbose() && index >=0 && index < 2)
+        {
+            /*for(int i =0; i<speciesNames_.size();i++)
+            {
+                Info << "Y_i = " << Y_[i].name() << endl;
+                Info << "concentrations = " << concentrations_[i][index][0] << endl;
+                Info << "partRho_[index][0] = " << partRho_[index][0] << endl;
+                Info << "rhofluid =" << rhofluid << endl;
+                Info << "Yfluid = " << Yfluid_[i] << endl;
+                Info << "partTemp_[index][0] = " << partTemp_[index][0] << endl;
+                Info << "Tfluid = " << Tfluid << endl  ;
+            }*/
+        }
+    }
+
+        // give DEM data
+        particleCloud_.dataExchangeM().giveData(partTempName_, "scalar-atom", partTemp_);
+        particleCloud_.dataExchangeM().giveData(partRhoName_, "scalar-atom", partRho_);
+        for (int i=0; i<speciesNames_.size();i++)
+        {
+            particleCloud_.dataExchangeM().giveData(Y_[i].name(),"scalar-atom",concentrations_[i]);
+        };
+
+        Info << "give data done" << endl;
+
+
+
+        // pull changeOfSpeciesMass_, transform onto fields changeOfSpeciesMassFields_, add them up on changeOfGasMassField_
+        changeOfGasMassField_.primitiveFieldRef() = 0.0;
+        changeOfGasMassField_.boundaryFieldRef() = 0.0;
+        for (int i=0; i<speciesNames_.size();i++)
+        {
+            particleCloud_.dataExchangeM().getData(mod_spec_names_[i],"scalar-atom",changeOfSpeciesMass_[i]);
+            changeOfSpeciesMassFields_[i].primitiveFieldRef() = 0.0;
+            changeOfSpeciesMassFields_[i].boundaryFieldRef() = 0.0;
+            particleCloud_.averagingM().setScalarSum
+            (
+                changeOfSpeciesMassFields_[i],
+                changeOfSpeciesMass_[i],
+                particleCloud_.particleWeights(),
+                NULL
+            );
+
+            // take care for implementation in LIGGGHTS: species produced from particles defined positive
+            changeOfSpeciesMassFields_[i].primitiveFieldRef() /= changeOfSpeciesMassFields_[i].mesh().V();
+            changeOfSpeciesMassFields_[i].correctBoundaryConditions();
+            changeOfGasMassField_ += changeOfSpeciesMassFields_[i];
+            Info << "total conversion of species" << speciesNames_[i] << " = " << gSum(changeOfSpeciesMassFields_[i]*1.0*changeOfSpeciesMassFields_[i].mesh().V()) << endl;
+        }
+        Info << "get data done" << endl;
+
+}
+
+/*tmp<Foam::fvScalarMatrix> species::Smi(const label i) const
+{
+    return tmp<fvScalarMatrix>(new fvScalarMatrix(changeOfSpeciesMassFields_[i], dimMass/dimTime));
+}
+
+tmp<Foam::fvScalarMatrix> species::Sm() const
+{
+    return tmp<fvScalarMatrix>(new fvScalarMatrix(changeOfGasMassField_, dimMass/dimTime));
+}*/
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

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

@@ -0,0 +1,144 @@
+/*---------------------------------------------------------------------------*\
+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 communicate species concentration and changes due to reactions
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef species_H
+#define species_H
+
+#include "fvCFD.H"
+#include "cfdemCloudEnergy.H"
+#include "chemistryModel.H"
+
+#include "HashPtrTable.H"
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class species Declaration
+\*---------------------------------------------------------------------------*/
+
+// for future use:
+// + communicate every N steps
+
+class species
+:
+    public chemistryModel
+{
+private:
+
+    dictionary propsDict_;
+
+    bool interpolation_;
+
+    const fvMesh& mesh_;
+
+    // Species dictionary - To be read where the species are from
+    dictionary specDict_;
+
+    wordList speciesNames_;
+    
+    wordList mod_spec_names_;
+
+    PtrList<volScalarField> Y_;
+
+    mutable List<double**> concentrations_;
+    
+    mutable List<double**> changeOfSpeciesMass_;
+
+    PtrList<volScalarField> changeOfSpeciesMassFields_;
+    
+    volScalarField changeOfGasMassField_;
+    
+    word tempFieldName_;
+
+    const volScalarField& tempField_;              // ref to gas temperature field
+
+    word partTempName_;
+
+    mutable double **partTemp_;                    // gas temperature at particle positions
+
+    word densityFieldName_;
+
+    const volScalarField& rho_;
+
+    word partRhoName_;
+    
+    mutable double **partRho_;                     // gas density at particle positions
+    
+    //word velFieldName_;                          // when velocity is to be calculated in the species source code
+
+    //const volVectorField& U_;
+
+    //word voidfractionFieldName_;                 // when voidfraction is to be calculated in the species source code
+
+    //const volScalarField& voidfraction_;
+
+
+    void allocateMyArrays() const;
+       
+public:
+
+    //- Runtime type information
+    TypeName("species");
+
+    // Constructors
+
+        //- Construct from components
+        species
+        (
+            const dictionary& dict,
+            cfdemCloudEnergy& sm
+        );
+
+
+    // Destructor
+
+        virtual ~species();
+
+
+    // Member Functions
+
+    void execute();
+
+    /*tmp<fvScalarMatrix> Smi(const label i) const;
+
+    tmp<fvScalarMatrix> Sm() const;*/
+    
+    
+
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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

@@ -0,0 +1,85 @@
+/*---------------------------------------------------------------------------*\
+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 "energyModel.H"
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(energyModel, 0);
+
+defineRunTimeSelectionTable(energyModel, dictionary);
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from components
+energyModel::energyModel
+(
+    const dictionary& dict,
+    cfdemCloudEnergy& sm
+)
+:
+    dict_(dict),
+    particleCloud_(sm),
+    transportProperties_
+    (
+        IOobject
+        (
+            "transportProperties",
+            sm.mesh().time().constant(),
+            sm.mesh(),
+            IOobject::MUST_READ_IF_MODIFIED,
+            IOobject::NO_WRITE
+        )
+    ),
+    kf0_
+    (
+        dimensionedScalar(transportProperties_.lookup("kf")).value()
+    ),
+    Cp_
+    (
+        dimensionedScalar(transportProperties_.lookup("Cp")).value()
+    )
+{}
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+energyModel::~energyModel()
+{}
+
+// * * * * * * * * * * * * * * * * Member Fct  * * * * * * * * * * * * * * * //
+
+scalar energyModel::Cp() const
+{
+    return Cp_; 
+}
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //
+

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

@@ -0,0 +1,123 @@
+/*---------------------------------------------------------------------------*\
+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
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef energyModel_H
+#define energyModel_H
+
+#include "fvCFD.H"
+#include "cfdemCloudEnergy.H"
+#include "probeModel.H"
+#include "interpolationCellPoint.H"
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class energyModel Declaration
+\*---------------------------------------------------------------------------*/
+
+class energyModel
+{
+
+protected:
+
+    // Protected data
+        const dictionary& dict_;
+
+        cfdemCloudEnergy& particleCloud_;
+	
+	IOdictionary transportProperties_;
+	
+	scalar kf0_;                              // fluid thermal conductivity [W/(m*K)]
+
+        scalar Cp_;                                  // specific heat capacity [W*s/(kg*K)]
+       
+public:
+
+    //- Runtime type information
+    TypeName("energyModel");
+
+    // Declare runtime constructor selection table
+
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            energyModel,
+            dictionary,
+            (
+                const dictionary& dict,
+                cfdemCloudEnergy& sm
+            ),
+            (dict,sm)
+        );
+
+
+    // Constructors
+
+        //- Construct from components
+        energyModel
+        (
+            const dictionary& dict,
+            cfdemCloudEnergy& sm            
+        );
+
+
+    // Destructor
+
+        virtual ~energyModel();
+
+
+    // Selector
+
+        static autoPtr<energyModel> New
+        (
+            const dictionary& dict,
+            cfdemCloudEnergy& sm,
+            word energyType
+        );
+
+
+    // Member Functions
+
+        virtual void addEnergyContribution(volScalarField&) const = 0;
+	
+	virtual void addEnergyCoefficient(volScalarField&) const = 0;
+	
+	virtual void calcEnergyContribution() = 0;
+	
+	virtual void postFlow() {}
+	
+	scalar Cp() const;
+
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
+

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

@@ -0,0 +1,68 @@
+/*---------------------------------------------------------------------------*\
+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 "energyModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+autoPtr<energyModel> energyModel::New
+(
+    const dictionary& dict,
+    cfdemCloudEnergy& sm,
+    word energyType
+)
+{
+    Info<< "Selecting energyModel "
+         << energyType << endl;
+
+    dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(energyType);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalError
+            << "energyModel::New(const dictionary&, const spray&) : "
+            << endl
+            << "    unknown energyModelType type "
+            << energyType
+            << ", constructor not in hash table" << endl << endl
+            << "    Valid energyModel types are :"
+            << endl;
+        Info<< dictionaryConstructorTablePtr_->toc()
+            << abort(FatalError);
+    }
+
+    return autoPtr<energyModel>(cstrIter()(dict,sm));
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //
+

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

@@ -0,0 +1,314 @@
+/*---------------------------------------------------------------------------*\
+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 "heatTransferGunn.H"
+#include "addToRunTimeSelectionTable.H"
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(heatTransferGunn, 0);
+
+addToRunTimeSelectionTable(energyModel, heatTransferGunn, dictionary);
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from components
+heatTransferGunn::heatTransferGunn
+(
+    const dictionary& dict,
+    cfdemCloudEnergy& sm
+)
+:
+    energyModel(dict,sm),
+    propsDict_(dict.subDict(typeName + "Props")),
+    interpolation_(propsDict_.lookupOrDefault<bool>("interpolation",false)),
+    QPartFluidName_(propsDict_.lookupOrDefault<word>("QPartFluidName","QPartFluid")),
+    QPartFluid_
+    (   IOobject
+        (
+            QPartFluidName_,
+            sm.mesh().time().timeName(),
+            sm.mesh(),
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        sm.mesh(),
+        dimensionedScalar("zero", dimensionSet(1,-1,-3,0,0,0,0), 0.0)
+    ),
+    partTempField_
+    (   IOobject
+        (
+            "particleTemp",
+            sm.mesh().time().timeName(),
+            sm.mesh(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        sm.mesh(),
+        dimensionedScalar("zero", dimensionSet(0,0,0,1,0,0,0), 0.0)
+    ),
+    partRelTempField_
+    (   IOobject
+        (
+            "particleRelTemp",
+            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)
+    ),
+    partRefTemp_("partRefTemp", dimensionSet(0,0,0,1,0,0,0), 0.0),
+    calcPartTempField_(propsDict_.lookupOrDefault<bool>("calcPartTempField",false)),
+    tempFieldName_(propsDict_.lookupOrDefault<word>("tempFieldName","T")),
+    tempField_(sm.mesh().lookupObject<volScalarField> (tempFieldName_)),
+    voidfractionFieldName_(propsDict_.lookupOrDefault<word>("voidfractionFieldName","voidfraction")),
+    voidfraction_(sm.mesh().lookupObject<volScalarField> (voidfractionFieldName_)),
+    maxSource_(1e30),
+    velFieldName_(propsDict_.lookupOrDefault<word>("velFieldName","U")),
+    U_(sm.mesh().lookupObject<volVectorField> (velFieldName_)),
+    densityFieldName_(propsDict_.lookupOrDefault<word>("densityFieldName","rho")),
+    rho_(sm.mesh().lookupObject<volScalarField> (densityFieldName_)),
+    partTempName_(propsDict_.lookup("partTempName")),
+    partTemp_(NULL),
+    partHeatFluxName_(propsDict_.lookup("partHeatFluxName")),
+    partHeatFlux_(NULL)
+{
+     allocateMyArrays();
+
+    if (propsDict_.found("maxSource"))
+    {
+        maxSource_=readScalar(propsDict_.lookup ("maxSource"));
+        Info << "limiting eulerian source field to: " << maxSource_ << endl;
+    }
+    if (calcPartTempField_)
+    {
+        if (propsDict_.found("partRefTemp"))
+	    partRefTemp_.value()=readScalar(propsDict_.lookup ("partRefTemp"));
+	partTempField_.writeOpt() = IOobject::AUTO_WRITE;
+	partRelTempField_.writeOpt() = IOobject::AUTO_WRITE;
+	partTempField_.write();
+	partRelTempField_.write();
+	Info <<  "Particle temperature field activated." << endl;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+heatTransferGunn::~heatTransferGunn()
+{
+    delete partTemp_;
+    delete partHeatFlux_;
+}
+
+// * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
+void heatTransferGunn::allocateMyArrays() const
+{
+    // get memory for 2d arrays
+    double initVal=0.0;
+    particleCloud_.dataExchangeM().allocateArray(partTemp_,initVal,1);  // field/initVal/with/lenghtFromLigghts
+    particleCloud_.dataExchangeM().allocateArray(partHeatFlux_,initVal,1);
+}
+
+// * * * * * * * * * * * * * * * * Member Fct  * * * * * * * * * * * * * * * //
+
+void heatTransferGunn::calcEnergyContribution()
+{
+   // realloc the arrays
+    allocateMyArrays();
+
+    // reset Scalar field
+    QPartFluid_.primitiveFieldRef() = 0.0;
+
+    // get DEM data
+    particleCloud_.dataExchangeM().getData(partTempName_,"scalar-atom",partTemp_);
+    
+    if(calcPartTempField_)
+    {       
+	partTempField_.primitiveFieldRef() = 0.0;
+        particleCloud_.averagingM().resetWeightFields();
+        particleCloud_.averagingM().setScalarAverage
+        (
+            partTempField_,
+            partTemp_,
+            particleCloud_.particleWeights(),
+            particleCloud_.averagingM().UsWeightField(),
+            NULL
+        );
+	volScalarField volP (1 - voidfraction_);
+	volScalarField weigthedTp (volP * partTempField_);
+	// average per cell-value, not per volume * cell-value
+	dimensionedScalar aveTemp = weigthedTp.average() / volP.average();
+	partRelTempField_ = (partTempField_ - aveTemp) / (aveTemp - partRefTemp_);
+    }
+
+    #ifdef compre
+       const volScalarField mufField = particleCloud_.turbulence().mu();
+    #else
+       const volScalarField mufField = particleCloud_.turbulence().nu()*rho_;
+    #endif
+    
+
+    // calc La based heat flux
+    scalar voidfraction(1);
+    vector Ufluid(0,0,0);
+    scalar Tfluid(0);
+    label cellI=0;
+    vector Us(0,0,0);
+    scalar ds(0);
+    scalar muf(0);
+    scalar magUr(0);
+    scalar Rep(0);
+    scalar Pr(0);
+    scalar Nup(0);
+
+
+    interpolationCellPoint<scalar> voidfractionInterpolator_(voidfraction_);
+    interpolationCellPoint<vector> UInterpolator_(U_);
+    interpolationCellPoint<scalar> TInterpolator_(tempField_);
+
+    for(int index = 0;index < particleCloud_.numberOfParticles(); ++index)
+    {
+            cellI = particleCloud_.cellIDs()[index][0];
+            if(cellI >= 0)
+            {
+                if(interpolation_)
+                {
+                    vector position = particleCloud_.position(index);
+                    voidfraction = voidfractionInterpolator_.interpolate(position,cellI);
+                    Ufluid = UInterpolator_.interpolate(position,cellI);
+                    Tfluid = TInterpolator_.interpolate(position,cellI);
+                }
+                else
+                {
+                    voidfraction = voidfraction_[cellI];
+                    Ufluid = U_[cellI];
+                    Tfluid = tempField_[cellI];
+                }
+
+                // calc relative velocity
+                Us = particleCloud_.velocity(index);
+                magUr = mag(Ufluid - Us);
+                ds = 2.*particleCloud_.radius(index);
+                muf = mufField[cellI];
+                Rep = ds * magUr * voidfraction * rho_[cellI]/ muf;
+                Pr = max(SMALL, Cp_ * muf / kf0_);
+
+                Nup = Nusselt(voidfraction, Rep, Pr);
+                
+
+                scalar h = kf0_ * Nup / ds;
+                scalar As = ds * ds * M_PI; // surface area of sphere
+
+                // calc convective heat flux [W]
+                heatFlux(index, h, As, Tfluid);
+                heatFluxCoeff(index, h, As);
+		
+		if(particleCloud_.verbose() && index >=0 && index <2)
+                {
+                    Info << "partHeatFlux = " << partHeatFlux_[index][0] << endl;
+                    Info << "magUr = " << magUr << endl;
+                    Info << "As = " << As << endl;
+                    Info << "muf = " << muf << endl;
+                    Info << "Rep = " << Rep << endl;
+                    Info << "Pr = " << Pr << endl;
+                    Info << "Nup = " << Nup << endl;
+                    Info << "voidfraction = " << voidfraction << endl;
+                    Info << "partTemp_[index][0] = " << partTemp_[index][0] << endl;
+                    Info << "Tfluid = " << Tfluid << endl  ;
+                }
+            }
+    }
+
+    particleCloud_.averagingM().setScalarSum
+    (
+        QPartFluid_,
+        partHeatFlux_,
+        particleCloud_.particleWeights(),
+        NULL
+    );
+
+    QPartFluid_.primitiveFieldRef() /= -QPartFluid_.mesh().V();
+
+    // limit source term
+    forAll(QPartFluid_,cellI)
+    {
+        scalar EuFieldInCell = QPartFluid_[cellI];
+
+        if(mag(EuFieldInCell) > maxSource_ )
+        {
+             QPartFluid_[cellI] = sign(EuFieldInCell) * maxSource_;
+        }
+    }
+    
+    QPartFluid_.correctBoundaryConditions();
+    
+    giveData(0);
+
+}
+
+void heatTransferGunn::addEnergyContribution(volScalarField& Qsource) const
+{
+    Qsource += QPartFluid_;
+}
+
+scalar heatTransferGunn::Nusselt(scalar voidfraction, scalar Rep, scalar Pr) const
+{
+    scalar Nup(0.0);
+    Nup = (7 - 10 * voidfraction + 5 * voidfraction * voidfraction) *
+                        (1 + 0.7 * Foam::pow(Rep,0.2) * Foam::pow(Pr,0.33)) +
+			(1.33 - 2.4 * voidfraction + 1.2 * voidfraction * voidfraction) *
+			Foam::pow(Rep,0.7) * Foam::pow(Pr,0.33);
+                        
+    return Nup;
+}
+
+void heatTransferGunn::heatFlux(label index, scalar h, scalar As, scalar Tfluid)
+{
+    partHeatFlux_[index][0] = h * As * (Tfluid - partTemp_[index][0]);
+}
+
+void heatTransferGunn::heatFluxCoeff(label index, scalar h, scalar As)
+{
+   //no heat transfer coefficient in explicit model   
+}
+
+void heatTransferGunn::giveData(int call)
+{
+    if(call == 0)
+    {
+        Info << "total convective particle-fluid heat flux [W] (Eulerian) = " << gSum(QPartFluid_*1.0*QPartFluid_.mesh().V()) << endl;
+
+        particleCloud_.dataExchangeM().giveData(partHeatFluxName_,"scalar-atom", partHeatFlux_);    
+    }
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //
+

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

@@ -0,0 +1,137 @@
+/*---------------------------------------------------------------------------*\
+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
+    Correlation for Nusselt number according to
+    Gunn, D. J. International Journal of Heat and Mass Transfer 21.4 (1978)
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef heatTransferGunn_H
+#define heatTransferGunn_H
+
+#include "fvCFD.H"
+#include "cfdemCloudEnergy.H"
+#include "energyModel.H"
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class heatTransferGunn Declaration
+\*---------------------------------------------------------------------------*/
+
+class heatTransferGunn
+:
+    public energyModel
+{
+protected:
+
+    dictionary propsDict_;
+    
+    bool interpolation_;
+    
+    word QPartFluidName_;
+        
+    volScalarField QPartFluid_;
+    
+    volScalarField partTempField_;
+    
+    volScalarField partRelTempField_;
+    
+    dimensionedScalar partRefTemp_;
+    
+    bool calcPartTempField_;
+
+    word tempFieldName_;
+
+    const volScalarField& tempField_;            // ref to temperature field
+
+    word voidfractionFieldName_;
+
+    const volScalarField& voidfraction_;    // ref to voidfraction field
+
+    scalar maxSource_;                           // max (limited) value of src field
+
+    word velFieldName_;
+
+    const volVectorField& U_;
+    
+    word densityFieldName_;
+    
+    const volScalarField& rho_;
+
+    word partTempName_;
+
+    mutable double **partTemp_;                          // Lagrangian array
+
+    word partHeatFluxName_;
+
+    mutable double **partHeatFlux_;              // Lagrangian array
+
+    void allocateMyArrays() const;
+    
+    scalar Nusselt(scalar, scalar, scalar) const;
+    
+    virtual void giveData(int);
+    
+    virtual void heatFlux(label, scalar, scalar, scalar);
+    
+    virtual void heatFluxCoeff(label, scalar, scalar);
+       
+public:
+
+    //- Runtime type information
+    TypeName("heatTransferGunn");
+
+    // Constructors
+
+        //- Construct from components
+        heatTransferGunn
+        (
+            const dictionary& dict,
+            cfdemCloudEnergy& sm
+        );
+
+
+    // Destructor
+
+        virtual ~heatTransferGunn();
+
+
+    // Member Functions
+
+        void addEnergyContribution(volScalarField&) const;
+	
+	void addEnergyCoefficient(volScalarField&) const {}
+	
+	void calcEnergyContribution();
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
+

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

@@ -0,0 +1,158 @@
+/*---------------------------------------------------------------------------*\
+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 "heatTransferGunnImplicit.H"
+#include "addToRunTimeSelectionTable.H"
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(heatTransferGunnImplicit, 0);
+
+addToRunTimeSelectionTable(energyModel, heatTransferGunnImplicit, dictionary);
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from components
+heatTransferGunnImplicit::heatTransferGunnImplicit
+(
+    const dictionary& dict,
+    cfdemCloudEnergy& sm
+)
+:
+    heatTransferGunn(dict,sm),
+    QPartFluidCoeffName_(propsDict_.lookupOrDefault<word>("QPartFluidCoeffName","QPartFluidCoeff")),
+    QPartFluidCoeff_
+    (   IOobject
+        (
+            QPartFluidCoeffName_,
+            sm.mesh().time().timeName(),
+            sm.mesh(),
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        sm.mesh(),
+        dimensionedScalar("zero", dimensionSet(1,-1,-3,-1,0,0,0), 0.0)
+    ),
+    partHeatFluxCoeff_(NULL)
+{
+    allocateMyArrays();
+    
+    // no limiting necessary for implicit heat transfer
+    maxSource_ = 1e30;
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+heatTransferGunnImplicit::~heatTransferGunnImplicit()
+{
+    delete partHeatFluxCoeff_;
+}
+
+// * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
+void heatTransferGunnImplicit::allocateMyArrays() const
+{
+    heatTransferGunn::allocateMyArrays();
+    double initVal=0.0;
+    particleCloud_.dataExchangeM().allocateArray(partHeatFluxCoeff_,initVal,1);
+}
+// * * * * * * * * * * * * * * * * Member Fct  * * * * * * * * * * * * * * * //
+
+void heatTransferGunnImplicit::calcEnergyContribution()
+{
+    heatTransferGunn::calcEnergyContribution();
+    
+    particleCloud_.averagingM().setScalarSum
+    (
+        QPartFluidCoeff_,
+        partHeatFluxCoeff_,
+        particleCloud_.particleWeights(),
+        NULL
+    );
+
+    QPartFluidCoeff_.primitiveFieldRef() /= -QPartFluidCoeff_.mesh().V();
+    
+    QPartFluidCoeff_.correctBoundaryConditions();
+
+}
+
+void heatTransferGunnImplicit::addEnergyCoefficient(volScalarField& Qsource) const
+{
+    Qsource += QPartFluidCoeff_;
+}
+
+void heatTransferGunnImplicit::heatFlux(label index, scalar h, scalar As, scalar Tfluid)
+{
+    partHeatFlux_[index][0] = -h * As * partTemp_[index][0];
+}
+
+void heatTransferGunnImplicit::heatFluxCoeff(label index, scalar h, scalar As)
+{
+   partHeatFluxCoeff_[index][0] = h * As;
+}
+
+void heatTransferGunnImplicit::giveData(int call)
+{
+    if(call == 1)
+    {
+        //Info << "total convective particle-fluid heat flux [W] (Eulerian) = " << gSum(QPartFluid_*1.0*QPartFluid_.mesh().V()) << endl;
+
+        particleCloud_.dataExchangeM().giveData(partHeatFluxName_,"scalar-atom", partHeatFlux_);    
+    }
+}
+
+void heatTransferGunnImplicit::postFlow()
+{
+    label cellI;
+    scalar Tfluid(0.0);
+    scalar Tpart(0.0);
+    interpolationCellPoint<scalar> TInterpolator_(tempField_);
+        
+    for(int index = 0;index < particleCloud_.numberOfParticles(); ++index)
+    {
+            cellI = particleCloud_.cellIDs()[index][0];
+            if(cellI >= 0)
+            {
+                if(interpolation_)
+                {
+                    vector position = particleCloud_.position(index);
+                    Tfluid = TInterpolator_.interpolate(position,cellI);
+                }
+                else
+                    Tfluid = tempField_[cellI];
+                
+                Tpart = partTemp_[index][0];
+                partHeatFlux_[index][0] = (Tfluid - Tpart) * partHeatFluxCoeff_[index][0];
+            }
+    }
+                
+    giveData(1);
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //
+

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

@@ -0,0 +1,100 @@
+/*---------------------------------------------------------------------------*\
+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
+    Correlation for Nusselt number according to
+    Gunn, D. J. International Journal of Heat and Mass Transfer 21.4 (1978)
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef heatTransferGunnImplicit_H
+#define heatTransferGunnImplicit_H
+
+#include "fvCFD.H"
+#include "cfdemCloudEnergy.H"
+#include "heatTransferGunn.H"
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class heatTransferGunnImplicit Declaration
+\*---------------------------------------------------------------------------*/
+
+class heatTransferGunnImplicit
+:
+    public heatTransferGunn
+{
+private:
+    
+    word QPartFluidCoeffName_;
+        
+    volScalarField QPartFluidCoeff_;
+
+    mutable double **partHeatFluxCoeff_;              // Lagrangian array
+
+    void allocateMyArrays() const;
+    
+    void giveData(int);
+    
+    void heatFlux(label, scalar, scalar, scalar);
+    
+    void heatFluxCoeff(label, scalar, scalar);
+       
+public:
+
+    //- Runtime type information
+    TypeName("heatTransferGunnImplicit");
+
+    // Constructors
+
+        //- Construct from components
+        heatTransferGunnImplicit
+        (
+            const dictionary& dict,
+            cfdemCloudEnergy& sm
+        );
+
+
+    // Destructor
+
+        virtual ~heatTransferGunnImplicit();
+
+
+    // Member Functions
+	
+	void addEnergyCoefficient(volScalarField&) const;
+	
+	void calcEnergyContribution();
+	
+	void postFlow();
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
+

src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferModel/heatTransfer/heatTransfer.C

@@ -1,82 +0,0 @@
-/*---------------------------------------------------------------------------*\
-    CFDEMcoupling - Open Source CFD-DEM coupling
-
-    CFDEMcoupling is part of the CFDEMproject
-    www.cfdem.com
-                                Christoph Goniva, christoph.goniva@cfdem.com
-                                Copyright 2009-2012 JKU Linz
-                                Copyright 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-Description
-    This code is designed to realize coupled CFD-DEM simulations using LIGGGHTS
-    and OpenFOAM(R). Note: this code is not part of OpenFOAM(R) (see DISCLAIMER).
-\*---------------------------------------------------------------------------*/
-
-#include "error.H"
-#include "heatTransfer.H"
-#include "addToRunTimeSelectionTable.H"
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-defineTypeNameAndDebug(heatTransfer, 0);
-
-addToRunTimeSelectionTable(heatTransferModel, heatTransfer, dictionary);
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-// Construct from components
-heatTransfer::heatTransfer
-(
-    const dictionary& dict,
-    cfdemCloudEnergy& sm
-)
-:
-    heatTransferModel(dict,sm),
-    dict_(dict),
-    particleCloud_(sm)
-{}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-heatTransfer::~heatTransfer()
-{}
-
-// * * * * * * * * * * * * * * * * Member Fct  * * * * * * * * * * * * * * * //
-void heatTransfer::partFluidHeatTransfer(volScalarField& EuField) const
-{
-    EuField.internalField() = 0.0;
-}
-
-scalar heatTransfer::partThermCond() const
-{
-    return 0; 
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// ************************************************************************* //

src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferModel/heatTransferModel/newHeatTransferModel.C

@@ -1,81 +0,0 @@
-/*---------------------------------------------------------------------------*\
-    CFDEMcoupling - Open Source CFD-DEM coupling
-
-    CFDEMcoupling is part of the CFDEMproject
-    www.cfdem.com
-                                Christoph Goniva, christoph.goniva@cfdem.com
-                                Copyright 2009-2012 JKU Linz
-                                Copyright 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-Description
-    This code is designed to realize coupled CFD-DEM simulations using LIGGGHTS
-    and OpenFOAM(R). Note: this code is not part of OpenFOAM(R) (see DISCLAIMER).
-\*---------------------------------------------------------------------------*/
-
-#include "error.H"
-
-#include "heatTransferModel.H"
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-autoPtr<heatTransferModel> heatTransferModel::New
-(
-    const dictionary& dict,
-    cfdemCloudEnergy& sm
-)
-{
-    word heatTransferModelType
-    (
-        dict.lookupOrDefault<word>("heatTransferModel","noHeatTransfer")
-    );
-    Info<< "Selecting heatTransferModel "
-         << heatTransferModelType << endl;
-
-    dictionaryConstructorTable::iterator cstrIter =
-        dictionaryConstructorTablePtr_->find(heatTransferModelType);
-
-    if (cstrIter == dictionaryConstructorTablePtr_->end())
-    {
-        FatalError
-            << "heatTransferModel::New(const dictionary&, const spray&) : "
-            << endl
-            << "    unknown heatTransferModelType type "
-            << heatTransferModelType
-            << ", constructor not in hash table" << endl << endl
-            << "    Valid heatTransferModel types are :"
-            << endl;
-        Info<< dictionaryConstructorTablePtr_->toc()
-            << abort(FatalError);
-    }
-
-    return autoPtr<heatTransferModel>(cstrIter()(dict,sm));
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// ************************************************************************* //

src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferModel/noHeatTransfer/noHeatTransfer.C

@@ -1,82 +0,0 @@
-/*---------------------------------------------------------------------------*\
-    CFDEMcoupling - Open Source CFD-DEM coupling
-
-    CFDEMcoupling is part of the CFDEMproject
-    www.cfdem.com
-                                Christoph Goniva, christoph.goniva@cfdem.com
-                                Copyright 2009-2012 JKU Linz
-                                Copyright 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-Description
-    This code is designed to realize coupled CFD-DEM simulations using LIGGGHTS
-    and OpenFOAM(R). Note: this code is not part of OpenFOAM(R) (see DISCLAIMER).
-\*---------------------------------------------------------------------------*/
-
-#include "error.H"
-#include "noHeatTransfer.H"
-#include "addToRunTimeSelectionTable.H"
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
-
-defineTypeNameAndDebug(noHeatTransfer, 0);
-
-addToRunTimeSelectionTable(heatTransferModel, noHeatTransfer, dictionary);
-
-// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
-
-// Construct from components
-noHeatTransfer::noHeatTransfer
-(
-    const dictionary& dict,
-    cfdemCloudEnergy& sm
-)
-:
-    heatTransferModel(dict,sm),
-    dict_(dict),
-    particleCloud_(sm)
-{}
-
-
-// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
-
-noHeatTransfer::~noHeatTransfer()
-{}
-
-// * * * * * * * * * * * * * * * * Member Fct  * * * * * * * * * * * * * * * //
-void noHeatTransfer::partFluidHeatTransfer(volScalarField& EuField) const
-{
-    EuField.internalField() = 0.0;
-}
-
-scalar noHeatTransfer::partThermCond() const
-{
-    return 0; 
-}
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// ************************************************************************* //

src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferModel/noHeatTransfer/noHeatTransfer.H

@@ -1,105 +0,0 @@
-/*---------------------------------------------------------------------------*\
-    CFDEMcoupling - Open Source CFD-DEM coupling
-
-    CFDEMcoupling is part of the CFDEMproject
-    www.cfdem.com
-                                Christoph Goniva, christoph.goniva@cfdem.com
-                                Copyright 2009-2012 JKU Linz
-                                Copyright 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, write to the Free Software Foundation,
-    Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-
-Description
-    This code is designed to realize coupled CFD-DEM simulations using LIGGGHTS
-    and OpenFOAM(R). Note: this code is not part of OpenFOAM(R) (see DISCLAIMER).
-
-Class
-    noHeatTransfer
-
-SourceFiles
-    noHeatTransfer.C
-
-\*---------------------------------------------------------------------------*/
-
-#ifndef noHeatTransfer_H
-#define noHeatTransfer_H
-
-#include "fvCFD.H"
-#include "cfdemCloudEnergy.H"
-#include "heatTransferModel.H"
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-namespace Foam
-{
-
-/*---------------------------------------------------------------------------*\
-                           Class noHeatTransfer Declaration
-\*---------------------------------------------------------------------------*/
-
-class noHeatTransfer
-:
-    public heatTransferModel
-{
-
-protected:
-
-    // Protected data
-        const dictionary& dict_;
-
-        cfdemCloudEnergy& particleCloud_;
-
-       
-public:
-
-    //- Runtime type information
-    TypeName("noHeatTransfer");
-
-    // Constructors
-
-        //- Construct from components
-        noHeatTransfer
-        (
-            const dictionary& dict,
-            cfdemCloudEnergy& sm
-        );
-
-
-    // Destructor
-
-        virtual ~noHeatTransfer();
-
-
-      // Member Functions
-
-
-        void partFluidHeatTransfer(volScalarField&) const;
-	scalar partThermCond() const;
-
-
-};
-
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-} // End namespace Foam
-
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-#endif
-
-// ************************************************************************* //

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

@@ -64,11 +64,7 @@ Archimedes::Archimedes
     propsDict_(dict.subDict(typeName + "Props")),
     twoDimensional_(false),
     gravityFieldName_(propsDict_.lookup("gravityFieldName")),
-    #if defined(version21) || defined(version16ext)
-        g_(sm.mesh().lookupObject<uniformDimensionedVectorField> (gravityFieldName_))
-    #elif defined(version15)
-        g_(dimensionedVector(sm.mesh().lookupObject<IOdictionary>("environmentalProperties").lookup(gravityFieldName_)).value())
-    #endif
+    g_(sm.mesh().lookupObject<uniformDimensionedVectorField> (gravityFieldName_))
 {
 
     //Append the field names to be probed

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

@@ -65,11 +65,7 @@ private:
 
     word gravityFieldName_;
 
-    #ifdef version21
-        const uniformDimensionedVectorField& g_;    // ref to gravity
-    #elif defined(version16ext) || defined(version15)
-        const dimensionedVector& g_;    // ref to gravity
-    #endif
+    const uniformDimensionedVectorField& g_;    // ref to gravity
 
 
 public:

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

@@ -66,11 +66,7 @@ ArchimedesIB::ArchimedesIB
     voidfractionFieldName_(propsDict_.lookup("voidfractionFieldName")), //mod by alice
     voidfractions_(sm.mesh().lookupObject<volScalarField> (voidfractionFieldName_)),//mod by alice
     gravityFieldName_(propsDict_.lookup("gravityFieldName")),
-    #if defined(version21) || defined(version16ext)
-        g_(sm.mesh().lookupObject<uniformDimensionedVectorField> (gravityFieldName_))
-    #elif  defined(version15)
-        g_(dimensionedVector(sm.mesh().lookupObject<IOdictionary>("environmentalProperties").lookup(gravityFieldName_)).value())
-    #endif
+    g_(sm.mesh().lookupObject<uniformDimensionedVectorField> (gravityFieldName_))
 {
     //Append the field names to be probed
     particleCloud_.probeM().initialize(typeName, "archimedesIBF.logDat");

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

@@ -71,12 +71,8 @@ private:
 
     word gravityFieldName_;
 
-    #ifdef version21
-        const uniformDimensionedVectorField& g_;    // ref to gravity
-    #elif defined(version16ext) || defined(version15)
-        const dimensionedVector& g_;    // ref to gravity
-    #endif
-
+    const uniformDimensionedVectorField& g_;    // ref to gravity
+    
 
 public:
 

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

@@ -0,0 +1,228 @@
+/*---------------------------------------------------------------------------*\
+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 "BeetstraDrag.H"
+#include "addToRunTimeSelectionTable.H"
+#include "averagingModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(BeetstraDrag, 0);
+
+addToRunTimeSelectionTable
+(
+    forceModel,
+    BeetstraDrag,
+    dictionary
+);
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from components
+BeetstraDrag::BeetstraDrag
+(
+    const dictionary& dict,
+    cfdemCloud& sm
+)
+:
+    forceModel(dict,sm),
+    propsDict_(dict.subDict(typeName + "Props")),
+    velFieldName_(propsDict_.lookup("velFieldName")),
+    U_(sm.mesh().lookupObject<volVectorField> (velFieldName_)),
+    voidfractionFieldName_(propsDict_.lookup("voidfractionFieldName")),
+    voidfraction_(sm.mesh().lookupObject<volScalarField> (voidfractionFieldName_)),
+    UsFieldName_(propsDict_.lookup("granVelFieldName")),
+    UsField_(sm.mesh().lookupObject<volVectorField> (UsFieldName_)),
+    scaleDia_(1.),
+    scaleDrag_(1.)
+{
+    //Append the field names to be probed
+    particleCloud_.probeM().initialize(typeName, "BeetstraDrag.logDat");
+    particleCloud_.probeM().vectorFields_.append("dragForce"); //first entry must  be the force
+    particleCloud_.probeM().vectorFields_.append("Urel");
+    particleCloud_.probeM().scalarFields_.append("Rep");
+    particleCloud_.probeM().scalarFields_.append("betaP");
+    particleCloud_.probeM().scalarFields_.append("voidfraction");
+    particleCloud_.probeM().writeHeader();
+
+    // init force sub model
+    setForceSubModels(propsDict_);
+    // define switches which can be read from dict
+    forceSubM(0).setSwitchesList(0,true); // activate treatExplicit switch
+    forceSubM(0).setSwitchesList(2,true); // activate implDEM switch
+    forceSubM(0).setSwitchesList(3,true); // activate search for verbose switch
+    forceSubM(0).setSwitchesList(4,true); // activate search for interpolate switch
+    forceSubM(0).setSwitchesList(8,true); // activate scalarViscosity switch
+    forceSubM(0).readSwitches();
+
+    particleCloud_.checkCG(true);
+    if (propsDict_.found("scale"))
+        scaleDia_=scalar(readScalar(propsDict_.lookup("scale")));
+    if (propsDict_.found("scaleDrag"))
+        scaleDrag_=scalar(readScalar(propsDict_.lookup("scaleDrag")));
+
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+BeetstraDrag::~BeetstraDrag()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void BeetstraDrag::setForce() const
+{
+    if (scaleDia_ > 1)
+        Info << "Beetstra using scale = " << scaleDia_ << endl;
+    else if (particleCloud_.cg() > 1){
+        scaleDia_=particleCloud_.cg();
+        Info << "Beetstra using scale from liggghts cg = " << scaleDia_ << endl;
+    }
+
+    const volScalarField& nufField = forceSubM(0).nuField();
+    const volScalarField& rhoField = forceSubM(0).rhoField();
+
+    vector position(0,0,0);
+    scalar voidfraction(1);
+    vector Ufluid(0,0,0);
+    vector drag(0,0,0);
+    label cellI=0;
+
+    vector Us(0,0,0);
+    vector Ur(0,0,0);
+    scalar ds(0);
+    scalar ds_scaled(0);
+    scalar scaleDia3 = scaleDia_*scaleDia_*scaleDia_;
+    scalar nuf(0);
+    scalar rho(0);
+    scalar magUr(0);
+    scalar Rep(0);
+    scalar localPhiP(0);
+
+    vector dragExplicit(0,0,0);
+    scalar dragCoefficient(0);
+    
+    interpolationCellPoint<scalar> voidfractionInterpolator_(voidfraction_);
+    interpolationCellPoint<vector> UInterpolator_(U_);
+
+    #include "setupProbeModel.H"
+
+    for(int index = 0;index <  particleCloud_.numberOfParticles(); ++index)
+    {
+            cellI = particleCloud_.cellIDs()[index][0];
+            drag = vector(0,0,0);
+            dragExplicit = vector(0,0,0);
+            Ufluid =vector(0,0,0);
+            voidfraction=0;
+            dragCoefficient = 0;
+
+            if (cellI > -1) // particle found
+            {
+
+                if( forceSubM(0).interpolation() )
+                {
+                    position     = particleCloud_.position(index);
+                    voidfraction = voidfractionInterpolator_.interpolate(position,cellI);
+                    Ufluid       = UInterpolator_.interpolate(position,cellI);
+                    //Ensure interpolated void fraction to be meaningful
+                    // Info << " --> voidfraction: " << voidfraction << endl;
+                    if(voidfraction>1.00) voidfraction = 1.0;
+                    if(voidfraction<0.10) voidfraction = 0.10;
+                }
+                else
+                {
+                    voidfraction = voidfraction_[cellI];
+                    Ufluid = U_[cellI];
+                }
+
+                Us = particleCloud_.velocity(index);
+                Ur = Ufluid-Us;
+                magUr = mag(Ur);
+                ds = 2*particleCloud_.radius(index);
+		ds_scaled = ds/scaleDia_;
+                rho = rhoField[cellI];
+                nuf = nufField[cellI];
+
+                Rep=0.0;
+                localPhiP = 1.0f-voidfraction+SMALL;
+
+                // calc particle's drag coefficient (i.e., Force per unit slip velocity and Stokes drag)
+               
+                Rep=ds_scaled*voidfraction*magUr/nuf+SMALL;
+                dragCoefficient = 10.0*localPhiP/(voidfraction*voidfraction) +
+                                  voidfraction*voidfraction*(1.0+1.5*Foam::sqrt(localPhiP)) +
+                                  0.413*Rep/(24*voidfraction*voidfraction)*(1.0/voidfraction+3*voidfraction*localPhiP+8.4*Foam::pow(Rep,-0.343))/
+                                  (1+Foam::pow(10,3*localPhiP)*Foam::pow(Rep,-0.5*(1+4*localPhiP)));
+
+                // calc particle's drag
+                dragCoefficient *= 3*M_PI*ds_scaled*nuf*rho*voidfraction*scaleDia3*scaleDrag_;
+                if (modelType_=="B")
+                    dragCoefficient /= voidfraction;
+
+                drag = dragCoefficient * Ur;
+
+                // explicitCorr
+                forceSubM(0).explicitCorr(drag,dragExplicit,dragCoefficient,Ufluid,U_[cellI],Us,UsField_[cellI],forceSubM(0).verbose());
+
+                if(forceSubM(0).verbose() && index >=0 && index <2)
+                {
+                    Pout << "cellI = " << cellI << endl;
+                    Pout << "index = " << index << endl;
+                    Pout << "Us = " << Us << endl;
+                    Pout << "Ur = " << Ur << endl;
+                    Pout << "ds = " << ds << endl;
+                    Pout << "ds/scale = " << ds/scaleDia_ << endl;
+                    Pout << "rho = " << rho << endl;
+                    Pout << "nuf = " << nuf << endl;
+                    Pout << "voidfraction = " << voidfraction << endl;
+                    Pout << "Rep = " << Rep << endl;
+                    Pout << "drag = " << drag << endl;
+                }
+
+                //Set value fields and write the probe
+                if(probeIt_)
+                {
+                    #include "setupProbeModelfields.H"
+                    vValues.append(drag);   //first entry must the be the force
+                    vValues.append(Ur);
+                    sValues.append(Rep);
+                    sValues.append(voidfraction);
+                    particleCloud_.probeM().writeProbe(index, sValues, vValues);
+                }
+            }
+
+            // write particle based data to global array
+            forceSubM(0).partToArray(index,drag,dragExplicit,Ufluid,dragCoefficient);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

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

@@ -0,0 +1,95 @@
+/*---------------------------------------------------------------------------*\
+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
+    drag law for monodisperse systems according to
+    Beetstra et al. AIChE J 53.2 (2007)
+
+SourceFiles
+    BeetstraDrag.C
+\*---------------------------------------------------------------------------*/
+
+#ifndef BeetstraDrag_H
+#define BeetstraDrag_H
+
+#include "forceModel.H"
+#include "interpolationCellPoint.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class BeetstraDrag Declaration
+\*---------------------------------------------------------------------------*/
+
+class BeetstraDrag
+:
+    public forceModel
+{
+private:
+    dictionary propsDict_;
+
+    word velFieldName_;
+
+    const volVectorField& U_;
+
+    word voidfractionFieldName_;
+
+    const volScalarField& voidfraction_;
+
+    word UsFieldName_;   
+
+    const volVectorField& UsField_;
+
+    mutable scalar scaleDia_;
+
+    mutable scalar scaleDrag_;
+
+public:
+
+    //- Runtime type information
+    TypeName("BeetstraDrag");
+
+
+    // Constructors
+
+        //- Construct from components
+        BeetstraDrag
+        (
+            const dictionary& dict,
+            cfdemCloud& sm
+        );
+
+    // Destructor
+
+        ~BeetstraDrag();
+
+
+    // Member Functions
+        void setForce() const;
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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

@@ -130,7 +130,7 @@ void LaEuScalarTemp::manipulateScalarField(volScalarField& EuField) const
     allocateMyArrays();
 
     // reset Scalar field
-    EuField.internalField() = 0.0;
+    EuField.primitiveFieldRef() = 0.0;
 
     // get DEM data
     particleCloud_.dataExchangeM().getData(partTempName_,"scalar-atom",partTemp_);
@@ -232,7 +232,7 @@ void LaEuScalarTemp::manipulateScalarField(volScalarField& EuField) const
     );
 
     // scale with -1/(Vcell*rho*Cp)
-    EuField.internalField() /= -rhoField.internalField()*Cp_*EuField.mesh().V();
+    EuField.primitiveFieldRef() /= -rhoField.internalField()*Cp_*EuField.mesh().V();
 
     // limit source term
     forAll(EuField,cellI)

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

@@ -120,7 +120,7 @@ void particleCellVolume::setForce() const
     {
         if(verbose_) Info << "particleCellVolume.C - setForce()" << endl;
 
-        scalarField_.internalField()=0.;
+        scalarField_.ref()=0.;
 
         // get reference to actual field
         const volScalarField& field = mesh_.lookupObject<volScalarField>(scalarFieldName_);
@@ -143,8 +143,8 @@ void particleCellVolume::setForce() const
                 scalarField2_[cellI] = 0.;
             }
         }
-        scalarField_.internalField() = gSum(scalarField_);
-        scalarField2_.internalField() = gSum(scalarField2_);
+        scalarField_.ref() = gSum(scalarField_);
+        scalarField2_.ref() = gSum(scalarField2_);
 
         if(verbose_)
         {

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

@@ -179,7 +179,7 @@ void volWeightedAverage::setForce() const
             MPI_Allreduce(&totVol, &totVol_all, 3, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
             integralValue = gSum(scalarFields_[i]);
             volWeightedAverage = integralValue / (totVol_all+SMALL);
-            scalarFields_[i].internalField() = volWeightedAverage;
+            scalarFields_[i].ref() = volWeightedAverage;
 
             if(verbose_)
             {
@@ -224,7 +224,7 @@ void volWeightedAverage::setForce() const
 
             MPI_Allreduce(&totVol, &totVol_all, 3, MPI_DOUBLE, MPI_SUM, MPI_COMM_WORLD);
             volWeightedAverage = gSum(vectorFields_[i]) / (totVol_all+SMALL);
-            vectorFields_[i].internalField() = volWeightedAverage;
+            vectorFields_[i].ref() = volWeightedAverage;
 
             if(verbose_)
             {

src/lagrangian/cfdemParticle/subModels/locateModel/engineSearch/engineSearch.C

@@ -61,14 +61,8 @@ engineSearch::engineSearch
 :
     locateModel(dict,sm),
     propsDict_(dict.subDict(typeName + "Props")),
-    //faceDecomp_(propsDict_.lookup("faceDecomp")),
     treeSearch_(propsDict_.lookup("treeSearch")),
-    #ifdef version16ext
-        searchEngine_(particleCloud_.mesh(),false) //(particleCloud_.mesh(),faceDecomp_)
-    #elif defined(version21)
-        searchEngine_(particleCloud_.mesh(),polyMesh::FACEPLANES) // FACEPLANES or FACECENTRETETS; FACEDIAGTETS not stable
-    #endif
-    //searchEngine_(particleCloud_.mesh(),faceDecomp_) // only 2.0.x
+    searchEngine_(particleCloud_.mesh(), polyMesh::FACE_PLANES)
 {}
 
 

src/lagrangian/cfdemParticle/subModels/locateModel/engineSearch/engineSearch.H

@@ -61,8 +61,6 @@ private:
 
     dictionary propsDict_;
 
-    //Switch faceDecomp_;
-
     Switch treeSearch_;
 
 protected:

src/lagrangian/cfdemParticle/subModels/locateModel/standardSearch/standardSearch.C

@@ -89,11 +89,7 @@ label standardSearch::findCell
             for(int i=0;i<3;i++) position[i] = positions[index][i];
 
             // find cell
-            #ifdef version16ext
-                cellIDs[index][0] = particleCloud_.mesh().findCell(position);
-            #elif defined(version21)
-                cellIDs[index][0] = particleCloud_.mesh().findCell(position, polyMesh::FACEPLANES);
-            #endif
+            cellIDs[index][0] = particleCloud_.mesh().findCell(position, polyMesh::FACE_PLANES);
         }
         else cellIDs[index][0] = -1;
     }
@@ -108,11 +104,7 @@ label standardSearch::findSingleCell
 ) const
 {
     // find cell
-    #ifdef version16ext
-        return particleCloud_.mesh().findCell(position);
-    #elif defined(version21)
-        return particleCloud_.mesh().findCell(position, polyMesh::FACEPLANES);
-    #endif
+    return particleCloud_.mesh().findCell(position, polyMesh::FACE_PLANES);
 }
 
 

src/lagrangian/cfdemParticle/subModels/locateModel/turboEngineSearch/turboEngineSearch.C

@@ -64,11 +64,7 @@ turboEngineSearch::turboEngineSearch
     propsDict_(dict.subDict(typeName + "Props")),
     treeSearch_(propsDict_.lookup("treeSearch")),
     bb_(particleCloud_.mesh().points(),false),
-    #ifdef version16ext
-        searchEngine_(particleCloud_.mesh(),false) //(particleCloud_.mesh(),faceDecomp_)
-    #elif defined(version21)
-        searchEngine_(particleCloud_.mesh(),polyMesh::FACEPLANES) // FACEPLANES or FACECENTRETETS; FACEDIAGTETS not stable
-    #endif
+    searchEngine_(particleCloud_.mesh(), polyMesh::FACE_PLANES)
 {}
 
 

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

@@ -145,18 +145,18 @@ tmp<volVectorField> explicitCouple::expMomSource() const
                 if (fNext_[cellI][i] > fLimit_[i]) fNext_[cellI][i] = fLimit_[i];
             }
         }
-        tsource() = fPrev_;
+        tsource.ref() = fPrev_;
     }else
     {
-        tsource() = (1 - tsf) * fPrev_ + tsf * fNext_;
+        tsource.ref() = (1 - tsf) * fPrev_ + tsf * fNext_;
     }
     return tsource;
 }
 
 void Foam::explicitCouple::resetMomSourceField() const
 {
-    fPrev_.internalField() = fNext_.internalField();
-    fNext_.internalField() = vector::zero;
+    fPrev_.ref() = fNext_.ref();
+    fNext_.primitiveFieldRef() = vector::zero;
 }
 
 inline vector Foam::explicitCouple::arrayToField(label cellI) const

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

@@ -166,10 +166,10 @@ tmp<volScalarField> implicitCouple::impMomSource() const
             // limiter
             if (KslNext_[cellI] > KslLimit_) KslNext_[cellI] = KslLimit_;
         }
-        tsource() = KslPrev_;
+        tsource.ref() = KslPrev_;
     }else
     {
-        tsource() = (1 - tsf) * KslPrev_ + tsf * KslNext_;
+        tsource.ref() = (1 - tsf) * KslPrev_ + tsf * KslNext_;
     }
 
     return tsource;
@@ -177,8 +177,8 @@ tmp<volScalarField> implicitCouple::impMomSource() const
 
 void Foam::implicitCouple::resetMomSourceField() const
 {
-    KslPrev_.internalField() = KslNext_.internalField();
-    KslNext_.internalField() = 0;
+    KslPrev_.ref() = KslNext_.ref();
+    KslNext_.primitiveFieldRef() = 0;
 }
 
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

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

@@ -100,7 +100,7 @@ void Foam::constDiffSmoothing::smoothen(volScalarField& fieldSrc) const
     volScalarField sSmoothField = sSmoothField_;
 
     sSmoothField.dimensions().reset(fieldSrc.dimensions());
-    sSmoothField.internalField()=fieldSrc.internalField();
+    sSmoothField.ref()=fieldSrc.internalField();
     sSmoothField.correctBoundaryConditions();
     sSmoothField.oldTime().dimensions().reset(fieldSrc.dimensions());
     sSmoothField.oldTime()=fieldSrc;
@@ -141,7 +141,7 @@ void Foam::constDiffSmoothing::smoothen(volVectorField& fieldSrc) const
     volVectorField vSmoothField = vSmoothField_;
 
     vSmoothField.dimensions().reset(fieldSrc.dimensions());
-    vSmoothField.internalField()=fieldSrc.internalField();
+    vSmoothField.ref()=fieldSrc.internalField();
     vSmoothField.correctBoundaryConditions();
     vSmoothField.oldTime().dimensions().reset(fieldSrc.dimensions());
     vSmoothField.oldTime()=fieldSrc;
@@ -176,7 +176,7 @@ void Foam::constDiffSmoothing::smoothenReferenceField(volVectorField& fieldSrc)
     volVectorField vSmoothField = vSmoothField_;
 
     vSmoothField.dimensions().reset(fieldSrc.dimensions());
-    vSmoothField.internalField()=fieldSrc.internalField();
+    vSmoothField.ref()=fieldSrc.internalField();
     vSmoothField.correctBoundaryConditions();
     vSmoothField.oldTime().dimensions().reset(fieldSrc.dimensions());
     vSmoothField.oldTime()=fieldSrc;
@@ -210,7 +210,7 @@ void Foam::constDiffSmoothing::smoothenReferenceField(volVectorField& fieldSrc)
     forAll(vSmoothField,cellI)
     {
         if ( mag(vSmoothField.oldTime().internalField()[cellI]) > 0.0f)  // have a vector in the OLD vSmoothField, so keep it!
-            NLarge()[cellI] = sourceStrength;
+            NLarge.ref()[cellI] = sourceStrength;
     }
 
     // do the smoothing

src/lagrangian/cfdemParticle/subModels/thermCondModel/SyamlalThermCond/SyamlalThermCond.C

@@ -0,0 +1,129 @@
+/*---------------------------------------------------------------------------*\
+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 "SyamlalThermCond.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(SyamlalThermCond, 0);
+
+addToRunTimeSelectionTable
+(
+    thermCondModel,
+    SyamlalThermCond,
+    dictionary
+);
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from components
+SyamlalThermCond::SyamlalThermCond
+(
+    const dictionary& dict,
+    cfdemCloud& sm
+)
+:
+    thermCondModel(dict,sm),
+    propsDict_(dict.subDict(typeName + "Props")),
+    voidfractionFieldName_(propsDict_.lookupOrDefault<word>("voidfractionFieldName","voidfraction")),
+    voidfraction_(sm.mesh().lookupObject<volScalarField> (voidfractionFieldName_)),
+    rhoFieldName_(propsDict_.lookupOrDefault<word>("rhoFieldName","rho")),
+    rho_(sm.mesh().lookupObject<volScalarField> (rhoFieldName_)),
+    wallQFactorName_(propsDict_.lookupOrDefault<word>("wallQFactorName","wallQFactor")),
+    wallQFactor_
+    (   IOobject
+        (
+            wallQFactorName_,
+            sm.mesh().time().timeName(),
+            sm.mesh(),
+            IOobject::READ_IF_PRESENT,
+            IOobject::NO_WRITE
+        ),
+        sm.mesh(),
+        dimensionedScalar("zero", dimensionSet(0,0,0,0,0,0,0), 1.0)
+    ),
+    hasWallQFactor_(false)
+{
+    if (wallQFactor_.headerOk())
+    {
+        Info << "Found field for scaling wall heat flux.\n" << endl;     
+        hasWallQFactor_ = true;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+SyamlalThermCond::~SyamlalThermCond()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+tmp<volScalarField> SyamlalThermCond::thermCond() const
+{
+    tmp<volScalarField> tvf
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "tmpThCond",
+                voidfraction_.instance(),
+                voidfraction_.mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            voidfraction_.mesh(),
+	    dimensionedScalar("zero", dimensionSet(1,1,-3,-1,0,0,0), 0.0)
+        )
+    );
+    volScalarField& svf = tvf.ref();
+    
+    svf = (1-sqrt(1-voidfraction_)) / (voidfraction_) * kf0_;
+    
+    // if a wallQFactor field is present, use it to scale heat transport through a patch
+    if (hasWallQFactor_)
+        forAll(wallQFactor_.boundaryField(), patchi)
+            svf.boundaryFieldRef()[patchi] *= wallQFactor_.boundaryField()[patchi];
+  
+    return tvf;
+}
+
+tmp<volScalarField> SyamlalThermCond::thermDiff() const
+{
+    return thermCond()/(rho_*Cp_);
+}
+
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/lagrangian/cfdemParticle/subModels/thermCondModel/SyamlalThermCond/SyamlalThermCond.H

@@ -0,0 +1,109 @@
+/*---------------------------------------------------------------------------*\
+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
+    thermal conductivity of fluid phase in presence of particles according to
+    Syamlal, M. and Gidaspow, M. AIChE Journal 31.1 (1985)
+
+Class
+    SyamlalThermCond
+
+SourceFiles
+    SyamlalThermCond.C
+
+\*---------------------------------------------------------------------------*/
+
+
+#ifndef SyamlalThermCond_H
+#define SyamlalThermCond_H
+
+#include "thermCondModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class SyamlalThermCond Declaration
+\*---------------------------------------------------------------------------*/
+
+class SyamlalThermCond
+:
+    public thermCondModel
+{
+  
+private:
+  
+    dictionary propsDict_;
+  
+    word voidfractionFieldName_;
+  
+    const volScalarField& voidfraction_;
+    
+    word rhoFieldName_;
+    
+    const volScalarField& rho_;
+    
+    word wallQFactorName_;
+    
+    // ratio of half-cell-size and near-wall film
+    volScalarField wallQFactor_;
+    
+    bool hasWallQFactor_;
+
+public:
+
+    //- Runtime type information
+    TypeName("SyamlalThermCond");
+
+
+    // Constructors
+
+        //- Construct from components
+        SyamlalThermCond
+        (
+            const dictionary& dict,
+            cfdemCloud& sm
+        );
+
+    // Destructor
+
+        ~SyamlalThermCond();
+
+
+    // Member Functions
+	
+	tmp<volScalarField> thermCond() const;
+	
+	tmp<volScalarField> thermDiff() const;
+       
+       
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

src/lagrangian/cfdemParticle/subModels/thermCondModel/noTherm/noThermCond.C

@@ -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
+
+\*---------------------------------------------------------------------------*/
+
+#include "error.H"
+
+#include "noThermCond.H"
+#include "addToRunTimeSelectionTable.H"
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(noThermCond, 0);
+
+addToRunTimeSelectionTable(thermCondModel, noThermCond, dictionary);
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from components
+noThermCond::noThermCond
+(
+    const dictionary& dict,
+    cfdemCloud& sm
+)
+:
+    thermCondModel(dict,sm)
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+noThermCond::~noThermCond()
+{}
+
+// * * * * * * * * * * * * * * * * Member Fct  * * * * * * * * * * * * * * * //
+
+tmp<volScalarField> noThermCond::thermCond() const
+{
+    tmp<volScalarField> tcond
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "fake1",
+                particleCloud_.mesh().time().timeName(),
+                particleCloud_.mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+             ),
+             particleCloud_.mesh(),
+             dimensionedScalar
+             (
+              "zero",
+              dimensionSet(0,2,-1,0,0,0,0),
+              0.0
+             )
+        )
+    );
+
+    return tcond;
+}
+
+tmp<volScalarField> noThermCond::thermDiff() const
+{
+    tmp<volScalarField> tdif
+    (
+        new volScalarField
+        (
+            IOobject
+            (
+                "fake2",
+                particleCloud_.mesh().time().timeName(),
+                particleCloud_.mesh(),
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+             ),
+             particleCloud_.mesh(),
+             dimensionedScalar
+             (
+              "zero",
+              dimensionSet(0,2,-1,0,0,0,0),
+              0.0
+             )
+        )
+    );
+
+    return tdif;
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/lagrangian/cfdemParticle/subModels/thermCondModel/noTherm/noThermCond.H

@@ -29,68 +29,56 @@ Description
     and OpenFOAM(R). Note: this code is not part of OpenFOAM(R) (see DISCLAIMER).
 
 Class
-    heatTransfer
+    noThermCond
 
 SourceFiles
-    heatTransfer.C
+    noThermCond.C
 
 \*---------------------------------------------------------------------------*/
 
-#ifndef heatTransfer_H
-#define heatTransfer_H
+#ifndef noThermCond_H
+#define noThermCond_H
 
+#include "thermCondModel.H"
+#include "cfdemCloud.H"
 #include "fvCFD.H"
-#include "cfdemCloudEnergy.H"
-#include "heatTransferModel.H"
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
 {
 
 /*---------------------------------------------------------------------------*\
-                           Class heatTransfer Declaration
+                           Class noDrag Declaration
 \*---------------------------------------------------------------------------*/
 
-class heatTransfer
+class noThermCond
 :
-    public heatTransferModel
+    public thermCondModel
 {
 
-protected:
-
-    // Protected data
-        const dictionary& dict_;
-
-        cfdemCloudEnergy& particleCloud_;
-
-       
 public:
 
     //- Runtime type information
-    TypeName("heatTransfer");
+    TypeName("off");
+
 
     // Constructors
 
         //- Construct from components
-        heatTransfer
+        noThermCond
         (
             const dictionary& dict,
-            cfdemCloudEnergy& sm
+            cfdemCloud& sm
         );
 
-
     // Destructor
 
-        virtual ~heatTransfer();
-
-
-      // Member Functions
-
-
-        void partFluidHeatTransfer(volScalarField&) const;
-	scalar partThermCond() const;
-
+        ~noThermCond();
 
+    // Member Functions
+        tmp<volScalarField> thermCond() const;
+        tmp<volScalarField> thermDiff() const;
 };
 
 

src/lagrangian/cfdemParticle/subModels/thermCondModel/thermCondModel/newThermCondModel.C

@@ -0,0 +1,72 @@
+/*---------------------------------------------------------------------------*\
+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 "thermCondModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+autoPtr<thermCondModel> thermCondModel::New
+(
+    const dictionary& dict,
+    cfdemCloud& sm
+)
+{
+    word thermCondModelType
+    (
+        dict.lookup("thermCondModel")
+    );
+
+    Info<< "Selecting thermCondModel "
+         << thermCondModelType << endl;
+
+
+    dictionaryConstructorTable::iterator cstrIter =
+        dictionaryConstructorTablePtr_->find(thermCondModelType);
+
+    if (cstrIter == dictionaryConstructorTablePtr_->end())
+    {
+        FatalError
+            << "thermCondModel::New(const dictionary&, const spray&) : "
+            << endl
+            << "    unknown thermCondModelType type "
+            << thermCondModelType
+            << ", constructor not in hash table" << endl << endl
+            << "    Valid thermCondModel types are :"
+            << endl;
+        Info<< dictionaryConstructorTablePtr_->toc()
+            << abort(FatalError);
+    }
+
+    return autoPtr<thermCondModel>(cstrIter()(dict,sm));
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/lagrangian/cfdemParticle/subModels/thermCondModel/thermCondModel/thermCondModel.C

@@ -0,0 +1,78 @@
+/*---------------------------------------------------------------------------*\
+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 "thermCondModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(thermCondModel, 0);
+
+defineRunTimeSelectionTable(thermCondModel, dictionary);
+
+// * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from components
+thermCondModel::thermCondModel
+(
+    const dictionary& dict,
+    cfdemCloud& sm
+)
+:
+    dict_(dict),
+    particleCloud_(sm),
+    transportProperties_
+    (
+        IOobject
+        (
+            "transportProperties",
+            sm.mesh().time().constant(),
+            sm.mesh(),
+            IOobject::MUST_READ_IF_MODIFIED,
+            IOobject::NO_WRITE
+        )
+    ),
+    kf0_(transportProperties_.lookup("kf")),
+    Cp_(transportProperties_.lookup("Cp"))
+//      kf0_(readScalar(transportProperties_.lookup("kf"))),
+//      Cp_(readScalar(transportProperties_.lookup("Cp")))
+{}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+thermCondModel::~thermCondModel()
+{}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

src/lagrangian/cfdemParticle/subModels/thermCondModel/thermCondModel/thermCondModel.H

@@ -0,0 +1,123 @@
+/*---------------------------------------------------------------------------*\
+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 for thermal conductivity of fluid phase in presence of particles
+
+Class
+    thermCondModel
+
+SourceFiles
+    thermCondModel.C
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef thermCondModel_H
+#define thermCondModel_H
+
+#include "fvCFD.H"
+#include "cfdemCloud.H"
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class thermCondModel Declaration
+\*---------------------------------------------------------------------------*/
+
+class thermCondModel
+{
+
+protected:
+
+    // Protected data
+        const dictionary& dict_;
+
+        cfdemCloud& particleCloud_;
+	
+	IOdictionary transportProperties_;
+	
+	dimensionedScalar kf0_;
+	
+	dimensionedScalar Cp_;
+	
+
+public:
+
+    //- Runtime type information
+    TypeName("thermCondModel");
+
+    // Declare runtime constructor selection table
+
+        declareRunTimeSelectionTable
+        (
+            autoPtr,
+            thermCondModel,
+            dictionary,
+            (
+                const dictionary& dict,
+                cfdemCloud& sm
+            ),
+            (dict,sm)
+        );
+
+
+    // Constructors
+
+        //- Construct from components
+        thermCondModel
+        (
+            const dictionary& dict,
+            cfdemCloud& sm
+        );
+
+
+    // Destructor
+
+        virtual ~thermCondModel();
+
+
+    // Selector
+
+        static autoPtr<thermCondModel> New
+        (
+            const dictionary& dict,
+            cfdemCloud& sm
+        );
+
+
+    // Member Functions
+	
+	virtual tmp<volScalarField> thermCond() const = 0;
+	
+	virtual tmp<volScalarField> thermDiff() const = 0;
+       
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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

@@ -91,7 +91,7 @@ void GaussVoidFraction::setvoidFraction(double** const& mask,double**& voidfract
 {
     reAllocArrays();
 
-    voidfractionNext_.internalField()=1;
+    voidfractionNext_.ref()=1;
 
     scalar radius(-1);
     scalar volume(0);

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

@@ -96,7 +96,7 @@ void IBVoidFraction::setvoidFraction(double** const& mask,double**& voidfraction
 
     reAllocArrays();
 
-    voidfractionNext_.internalField()=1;
+    voidfractionNext_.ref()=1;
 
     for(int index=0; index< particleCloud_.numberOfParticles(); index++)
     {

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

@@ -90,7 +90,7 @@ void bigParticleVoidFraction::setvoidFraction(double** const& mask,double**& voi
 {
     reAllocArrays();
 
-    voidfractionNext_.internalField()=1;
+    voidfractionNext_.ref()=1;
 
     scalar radius(-1);
     scalar volume(0);

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

@@ -169,19 +169,19 @@ tmp<volScalarField> Foam::voidFractionModel::voidFractionInterp() const
     scalar tsf = particleCloud_.dataExchangeM().timeStepFraction();
     if(1-tsf < 1e-4 && particleCloud_.dataExchangeM().couplingStep() > 1) //tsf==1
     {
-        tsource() = voidfractionPrev_;
+        tsource.ref() = voidfractionPrev_;
     }
     else
     {
-        tsource() = (1 - tsf) * voidfractionPrev_ + tsf * voidfractionNext_;
+        tsource.ref() = (1 - tsf) * voidfractionPrev_ + tsf * voidfractionNext_;
     }
     return tsource;
 }
 
 void Foam::voidFractionModel::resetVoidFractions() const
 {
-    voidfractionPrev_.internalField() = voidfractionNext_.internalField();
-    voidfractionNext_.internalField() = 1;
+    voidfractionPrev_.ref() = voidfractionNext_.ref();
+    voidfractionNext_.ref() = 1;
 }
 
 /*void Foam::voidFractionModel::undoVoidFractions(double**const& mask) const

tutorials/cfdemSolverPiso/ErgunTestCG/CFD/0/U

@@ -30,29 +30,7 @@ boundaryField
 
     inlet
     {
-        /*type            flowRateInletVelocity;
-        flowRate        0.001;
-        value           uniform (0 0 0);*/
-
-        /*type            fixedValue;
-        value           uniform (0 0 0.0001);*/
-
-        //type            zeroGradient;
-
-/*        type            groovyBC;
-        variables       "Uend=vector(0,0,0.02);tEnd=0.1;";
-        valueExpression "((time() < tEnd) ? Uend/tEnd*time():Uend)";
-        value           uniform (0 0 0);*/
-
-        /*// 2.0.x, 1.6,ext
-        type            timeVaryingUniformFixedValue;
-        fileName        "steps_0p1s";
-        outOfBounds     clamp;
-        value           uniform (0 0 0);*/
-
-        // 2.1.x
-		//type              uniformFixedValue;
-        type            uniformFixedValueVoidfraction;
+        type            uniformFixedValue;
         uniformValue     table
         (
             (0.000 (0 0 0.002))
@@ -62,9 +40,6 @@ boundaryField
 
     outlet
     {
-        /*type            fluxCorrectedVelocity; //inletOutlet;
-        value           uniform (0 0 0);
-        inletValue      uniform (0 0 0);*/
         type            zeroGradient;
     }
 }

tutorials/cfdemSolverPiso/ErgunTestCG/CFD/constant/RASProperties

@@ -1,25 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  1.6                                   |
-|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       dictionary;
-    location    "constant";
-    object      RASProperties;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-RASModel        laminar;
-
-turbulence      off;
-
-printCoeffs     on;
-
-
-// ************************************************************************* //

tutorials/cfdemSolverPiso/ErgunTestCG/CFD/constant/couplingProperties

@@ -70,7 +70,7 @@ momCoupleModels
     implicitCouple
 );
 
-turbulenceModelType "RASProperties";//"LESProperties";//
+turbulenceModelType "turbulenceProperties";
 
 //===========================================================================//
 // sub-model properties

tutorials/cfdemSolverPiso/ErgunTestCG/CFD/constant/turbulenceProperties

@@ -15,7 +15,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-simulationType  RASModel;
+simulationType  laminar;
 
 
 // ************************************************************************* //

tutorials/cfdemSolverPiso/ErgunTestCG/CFD/steps_0p1s

@@ -1,22 +0,0 @@
-(
-    (0.000 (0 0 0.002))
-    (0.010 (0 0 0.002))
-    (0.011 (0 0 0.004))
-    (0.020 (0 0 0.004))
-    (0.021 (0 0 0.006))
-    (0.030 (0 0 0.006))
-    (0.031 (0 0 0.008))
-    (0.040 (0 0 0.008))
-    (0.041 (0 0 0.010))
-    (0.050 (0 0 0.010))
-    (0.051 (0 0 0.012))
-    (0.060 (0 0 0.012))
-    (0.061 (0 0 0.014))
-    (0.070 (0 0 0.014))
-    (0.071 (0 0 0.016))
-    (0.080 (0 0 0.016))
-    (0.081 (0 0 0.018))
-    (0.090 (0 0 0.018))
-    (0.091 (0 0 0.020))
-    (0.100 (0 0 0.020))
-)

tutorials/cfdemSolverPiso/ErgunTestCG/CFD/system/controlDict

@@ -53,8 +53,7 @@ maxCo 0.1;
 libs ( "libfiniteVolumeCFDEM.so" );
 
 functions
-(
-
+{
     probes
     {
         type        probes;
@@ -94,76 +93,5 @@ functions
         outputControl   timeStep;//outputTime;
         outputInterval  1;
     }
-
-/*
-    probes_inletFlux
-    {
-        type            faceSource;
-        functionObjectLibs ("libfieldFunctionObjects.so");
-        enabled         true;
-        outputControl   timeStep;   //outputTime;
-        log             true;       // log to screen?
-        valueOutput     true;       // Write values at run-time output times?
-        source          patch;   // Type of face source:
-        sourceName      inlet;         // faceZone name, see below
-        operation       sum;
-        fields
-        (
-            phi
-        );
-    }
-
-    probes_outletFlux
-    {
-        type            faceSource;
-        functionObjectLibs ("libfieldFunctionObjects.so");
-        enabled         true;
-        outputControl   timeStep;   //outputTime;
-        log             true;       // log to screen?
-        valueOutput     true;       // Write values at run-time output times?
-        source          patch;   // Type of face source:
-        sourceName      outlet;         // faceZone name, see below
-        operation       sum;
-        fields
-        (
-            phi
-        );
-    }
-
-    probes_wallFlux
-    {
-        type            faceSource;
-        functionObjectLibs ("libfieldFunctionObjects.so");
-        enabled         true;
-        outputControl   timeStep;   //outputTime;
-        log             true;       // log to screen?
-        valueOutput     true;       // Write values at run-time output times?
-        source          patch;   // Type of face source:
-        sourceName      wall;         // faceZone name, see below
-        operation       sum;
-        fields
-        (
-            phi
-        );
-    }*/
-   /*pressureDrop
-   {
-     type patchAverage;
-     functionObjectLibs
-     (
-         "libsimpleFunctionObjects.so"
-     );
-     verbose true;
-     patches
-     (
-         inlet
-         outlet
-     );
-     fields
-     (
-         p
-     );
-     factor  1;
-   }*/
-);
+};
 // ************************************************************************* //

tutorials/cfdemSolverPiso/ErgunTestMPI/CFD/0/U

@@ -30,40 +30,16 @@ boundaryField
 
     inlet
     {
-        /*type            flowRateInletVelocity;
-        flowRate        0.001;
-        value           uniform (0 0 0);*/
-
-        /*type            fixedValue;
-        value           uniform (0 0 0.0001);*/
-
-        //type            zeroGradient;
-
-/*        type            groovyBC;
-        variables       "Uend=vector(0,0,0.02);tEnd=0.1;";
-        valueExpression "((time() < tEnd) ? Uend/tEnd*time():Uend)";
-        value           uniform (0 0 0);*/
-
-        /*// 2.0.x, 1.6,ext
-        type            timeVaryingUniformFixedValue;
-        fileName        "steps_0p1s";
-        outOfBounds     clamp;
-        value           uniform (0 0 0);*/
-
-        // 2.1.x
         type            uniformFixedValue;
         uniformValue     table
         (
             (0.000 (0 0 0.002))
-            (0.100 (0 0 0.020))
+            (0.500 (0 0 0.020))
         );
     }
 
     outlet
     {
-        /*type            fluxCorrectedVelocity; //inletOutlet;
-        value           uniform (0 0 0);
-        inletValue      uniform (0 0 0);*/
         type            zeroGradient;
     }
 }

tutorials/cfdemSolverPiso/ErgunTestMPI/CFD/0/p

@@ -16,7 +16,7 @@ FoamFile
 
 dimensions      [0 2 -2 0 0 0 0];
 
-internalField   uniform 1;
+internalField   uniform 0e0;
 
 boundaryField
 {
@@ -29,13 +29,9 @@ boundaryField
     inlet
     {
         type            zeroGradient;
-        //type            fixedValue;
-        //value           uniform 100000;
     }
     outlet
     {
-        //type            zeroGradient;
-
         type            fixedValue;
         value           $internalField;
     }

tutorials/cfdemSolverPiso/ErgunTestMPI/CFD/constant/RASProperties

@@ -1,25 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  1.6                                   |
-|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       dictionary;
-    location    "constant";
-    object      RASProperties;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-RASModel        laminar;
-
-turbulence      off;
-
-printCoeffs     on;
-
-
-// ************************************************************************* //

tutorials/cfdemSolverPiso/ErgunTestMPI/CFD/constant/couplingProperties

@@ -55,12 +55,12 @@ smoothingModel off;// localPSizeDiffSmoothing;// constDiffSmoothing; //
 
 forceModels
 (
-    //GidaspowDrag
+    GidaspowDrag
     //BeetstraDrag
     //DiFeliceDrag
     gradPForce
     viscForce
-    KochHillDrag
+    //KochHillDrag
     //DEMbasedDrag
     //RongDrag
     //Archimedes
@@ -75,7 +75,7 @@ momCoupleModels
     implicitCouple
 );
 
-turbulenceModelType "RASProperties";//"LESProperties";//
+turbulenceModelType "turbulenceProperties";//"LESProperties";//
 
 //===========================================================================//
 // sub-model properties
@@ -143,6 +143,7 @@ GidaspowDragProps
 {
     verbose true;
     velFieldName "U";
+    granVelFieldName "Us";
     voidfractionFieldName "voidfraction";
     interpolation true;
     phi 1;

tutorials/cfdemSolverPiso/ErgunTestMPI/CFD/constant/turbulenceProperties

@@ -15,7 +15,7 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-simulationType  RASModel;
+simulationType  laminar;
 
 
 // ************************************************************************* //

tutorials/cfdemSolverPiso/ErgunTestMPI/CFD/octave/totalPressureDrop.m

@@ -9,7 +9,7 @@ rhoG = 10			% density in kg/m3
 %path = '../probes/0/p'; % 2.1.x
 path = '../postProcessing/probes/0/p'; % 2.2.x
 columns=22;
-headerlines=4;
+headerlines=23;%4;
 data = loaddata(path,columns,headerlines);
 data=transpose(data);
 [x,y]=size(data)

tutorials/cfdemSolverPiso/ErgunTestMPI/CFD/system/controlDict

@@ -23,7 +23,7 @@ startTime       0;
 
 stopAt          endTime;
 
-endTime         0.1;
+endTime         0.5;
 
 deltaT          0.0005;
 
@@ -49,10 +49,9 @@ adjustTimeStep  no;
 
 maxCo 0.1;
 
-//libs ( "libgroovyBC.so" "libfiniteVolumeCFDEM.so");
 
 functions
-(
+{
 
     probes
     {
@@ -94,24 +93,5 @@ functions
         outputInterval  1;
     }
 
-   /*pressureDrop
-   {
-     type patchAverage;
-     functionObjectLibs
-     (
-         "libsimpleFunctionObjects.so"
-     );
-     verbose true;
-     patches
-     (
-         inlet
-         outlet
-     );
-     fields
-     (
-         p
-     );
-     factor  1;
-   }*/
-);
+};
 // ************************************************************************* //

tutorials/cfdemSolverPiso/ErgunTestMPI_cgs/CFD/0/U

@@ -23,34 +23,11 @@ boundaryField
 
     wall
     {
-        //type            fixedValue;
-        //value           uniform (0 0 0);
         type            slip;
     }
 
     inlet
     {
-        /*type            flowRateInletVelocity;
-        flowRate        0.001;
-        value           uniform (0 0 0);*/
-
-        /*type            fixedValue;
-        value           uniform (0 0 0.0001);*/
-
-        //type            zeroGradient;
-
-        /*type            groovyBC;
-        variables       "Uend=vector(0,0,2);tEnd=0.1;";
-        valueExpression "((time() < tEnd) ? Uend/tEnd*time():Uend)";
-        value           uniform (0 0 0);*/
-
-        /*// 2.0.x, 1.6,ext
-        type            timeVaryingUniformFixedValue;
-        fileName        "steps_0p1s";
-        outOfBounds     clamp;
-        value           uniform (0 0 0);*/
-
-        // 2.1.x
         type            uniformFixedValue;
         uniformValue     table
         (
@@ -79,9 +56,6 @@ boundaryField
 
     outlet
     {
-        /*type            fluxCorrectedVelocity; //inletOutlet;
-        value           uniform (0 0 0);
-        inletValue      uniform (0 0 0);*/
         type            zeroGradient;
     }
 }

tutorials/cfdemSolverPiso/ErgunTestMPI_cgs/CFD/constant/RASProperties

@@ -1,25 +0,0 @@
-/*--------------------------------*- C++ -*----------------------------------*\
-| =========                 |                                                 |
-| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
-|  \\    /   O peration     | Version:  1.6                                   |
-|   \\  /    A nd           | Web:      www.OpenFOAM.org                      |
-|    \\/     M anipulation  |                                                 |
-\*---------------------------------------------------------------------------*/
-FoamFile
-{
-    version     2.0;
-    format      ascii;
-    class       dictionary;
-    location    "constant";
-    object      RASProperties;
-}
-// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
-
-RASModel        laminar;
-
-turbulence      off;
-
-printCoeffs     on;
-
-
-// ************************************************************************* //