Merge with develop.

applications/solvers/cfdemSolverRhoPimple/EEqn.H

@@ -53,7 +53,8 @@
 
     Info<< "T max/min : " << max(T).value() << " " << min(T).value() << endl;
 
-    particleCloud.clockM().start(31,"postFlow");
+    particleCloud.clockM().start(31,"energySolve");
-    particleCloud.postFlow();
+    particleCloud.solve();
-    particleCloud.clockM().stop("postFlow");
+    particleCloud.clockM().stop("energySolve");
+
 }

applications/solvers/cfdemSolverRhoPimple/cfdemSolverRhoPimple.C

@@ -139,6 +139,10 @@ int main(int argc, char *argv[])
             }
         }
 
+        particleCloud.clockM().start(31,"postFlow");
+        particleCloud.postFlow();
+        particleCloud.clockM().stop("postFlow");
+
         runTime.write();
 
 

applications/solvers/cfdemSolverRhoSimple/EEqn.H

@@ -0,0 +1,57 @@
+// contributions to internal energy equation can be found in
+// Crowe et al.: "Multiphase flows with droplets and particles", CRC Press 1998
+{
+    // dim he = J / kg
+    volScalarField& he = thermo.he();
+    particleCloud.energyContributions(Qsource);
+    particleCloud.energyCoefficients(QCoeff);
+
+    //thDiff=particleCloud.thermCondM().thermDiff();   
+    thCond=particleCloud.thermCondM().thermCond(); 
+
+    addSource = 
+        (
+            he.name() == "e"
+          ? 
+            fvc::div(phi, K) +
+            fvc::div
+            (
+                fvc::absolute(phi/fvc::interpolate(rho), voidfraction*U),
+                p,
+                "div(phiv,p)"
+            )
+          : fvc::div(phi, K)
+        );
+
+    Cpv = he.name() == "e" ? thermo.Cv() : thermo.Cp();
+
+
+    fvScalarMatrix EEqn
+    (
+        fvm::div(phi, he)
+      + addSource
+      - Qsource
+      - fvm::Sp(QCoeff/Cpv, he)
+      - fvm::laplacian(voidfraction*thCond/Cpv,he)
+     ==
+        fvOptions(rho, he)
+    );
+
+
+    EEqn.relax();
+
+    fvOptions.constrain(EEqn);
+
+    EEqn.solve();
+
+    fvOptions.correct(he);
+
+    thermo.correct();
+
+    Info<< "T max/min : " << max(T).value() << " " << min(T).value() << endl;
+
+
+    particleCloud.clockM().start(31,"energySolve");
+    particleCloud.solve();
+    particleCloud.clockM().stop("energySolve");
+}

applications/solvers/cfdemSolverRhoSimple/Make/files

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

applications/solvers/cfdemSolverRhoSimple/Make/options

@@ -0,0 +1,32 @@
+include $(CFDEM_ADD_LIBS_DIR)/additionalLibs
+
+PFLAGS+= -Dcompre
+
+EXE_INC = \
+     $(PFLAGS) \
+    -I$(CFDEM_OFVERSION_DIR) \
+    -I$(LIB_SRC)/transportModels/compressible/lnInclude \
+    -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
+    -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \
+    -I$(LIB_SRC)/finiteVolume/cfdTools \
+    -I$(LIB_SRC)/finiteVolume/lnInclude \
+    -I$(LIB_SRC)/meshTools/lnInclude \
+    -I$(LIB_SRC)/sampling/lnInclude \
+    -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/lnInclude \
+    -I$(CFDEM_SRC_DIR)/lagrangian/cfdemParticle/cfdTools \
+
+EXE_LIBS = \
+    -L$(CFDEM_LIB_DIR)\
+    -lcompressibleTransportModels \
+    -lfluidThermophysicalModels \
+    -lspecie \
+    -lturbulenceModels \
+    -lcompressibleTurbulenceModels \
+    -lfiniteVolume \
+    -lmeshTools \
+    -lsampling \
+    -lfvOptions \
+    -l$(CFDEM_LIB_COMP_NAME) \
+     $(CFDEM_ADD_LIB_PATHS) \
+     $(CFDEM_ADD_LIBS)

applications/solvers/cfdemSolverRhoSimple/UEqn.H

@@ -0,0 +1,30 @@
+// Solve the Momentum equation
+particleCloud.otherForces(fOther);
+
+tmp<fvVectorMatrix> tUEqn
+(
+    fvm::div(phi, U)
+  + particleCloud.divVoidfractionTau(U, voidfraction)
+  + fvm::Sp(Ksl,U)
+  - fOther
+ ==
+    fvOptions(rho, U)
+);
+fvVectorMatrix& UEqn = tUEqn.ref();
+
+UEqn.relax();
+
+fvOptions.constrain(UEqn);
+
+if (modelType=="B" || modelType=="Bfull")
+{
+    solve(UEqn == -fvc::grad(p)+ Ksl*Us);
+}
+else
+{
+    solve(UEqn == -voidfraction*fvc::grad(p)+ Ksl*Us);
+}
+
+fvOptions.correct(U);
+
+K = 0.5*magSqr(U);

applications/solvers/cfdemSolverRhoSimple/cfdemSolverRhoSimple.C

@@ -0,0 +1,142 @@
+/*---------------------------------------------------------------------------*\
+License
+
+    This is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    This code is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this code.  If not, see <http://www.gnu.org/licenses/>.
+
+    Copyright (C) 2015- Thomas Lichtenegger, JKU Linz, Austria
+
+Application
+    cfdemSolverRhoPimple
+
+Description
+    Transient solver for compressible flow using the flexible PIMPLE (PISO-SIMPLE)
+    algorithm.
+    Turbulence modelling is generic, i.e. laminar, RAS or LES may be selected.
+    The code is an evolution of the solver rhoPimpleFoam in OpenFOAM(R) 4.x,
+    where additional functionality for CFD-DEM coupling is added.
+\*---------------------------------------------------------------------------*/
+
+#include "fvCFD.H"
+#include "psiThermo.H"
+#include "turbulentFluidThermoModel.H"
+#include "bound.H"
+#include "simpleControl.H"
+#include "fvOptions.H"
+#include "localEulerDdtScheme.H"
+#include "fvcSmooth.H"
+
+#include "cfdemCloudEnergy.H"
+#include "implicitCouple.H"
+#include "clockModel.H"
+#include "smoothingModel.H"
+#include "forceModel.H"
+#include "thermCondModel.H"
+#include "energyModel.H"
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+int main(int argc, char *argv[])
+{
+    #include "postProcess.H"
+
+    #include "setRootCase.H"
+    #include "createTime.H"
+    #include "createMesh.H"
+    #include "createControl.H"
+    #include "createTimeControls.H"
+    #include "createRDeltaT.H"
+    #include "initContinuityErrs.H"
+    #include "createFields.H"
+    #include "createFieldRefs.H"
+    #include "createFvOptions.H"
+
+    // create cfdemCloud
+    #include "readGravitationalAcceleration.H"
+    cfdemCloudEnergy particleCloud(mesh);
+    #include "checkModelType.H"
+
+    turbulence->validate();
+
+    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+    Info<< "\nStarting time loop\n" << endl;
+
+    while (simple.loop())
+    {
+        particleCloud.clockM().start(1,"Global");
+
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        // do particle stuff
+        particleCloud.clockM().start(2,"Coupling");
+        bool hasEvolved = particleCloud.evolve(voidfraction,Us,U);
+
+        if(hasEvolved)
+        {
+            particleCloud.smoothingM().smoothen(particleCloud.forceM(0).impParticleForces());
+        }
+
+        Info << "update Ksl.internalField()" << endl;
+        Ksl = particleCloud.momCoupleM(0).impMomSource();
+        Ksl.correctBoundaryConditions();
+
+        //Force Checks
+        vector fTotal(0,0,0);
+        vector fImpTotal = sum(mesh.V()*Ksl.primitiveFieldRef()*(Us.primitiveFieldRef()-U.primitiveFieldRef()));
+        reduce(fImpTotal, sumOp<vector>());
+        Info << "TotalForceExp: " << fTotal << endl;
+        Info << "TotalForceImp: " << fImpTotal << endl;
+
+        #include "solverDebugInfo.H"
+        particleCloud.clockM().stop("Coupling");
+
+        particleCloud.clockM().start(26,"Flow");
+
+        volScalarField rhoeps("rhoeps",rho*voidfraction);
+        // Pressure-velocity SIMPLE corrector
+
+        #include "UEqn.H"
+
+
+        // besides this pEqn, OF offers a "simple consistent"-option
+        #include "pEqn.H"
+        rhoeps=rho*voidfraction;
+
+        #include "EEqn.H"
+
+        turbulence->correct();
+
+        particleCloud.clockM().start(32,"postFlow");
+        if(hasEvolved) particleCloud.postFlow();
+        particleCloud.clockM().stop("postFlow");
+
+        runTime.write();
+
+
+        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
+            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
+            << nl << endl;
+
+        particleCloud.clockM().stop("Flow");
+        particleCloud.clockM().stop("Global");
+    }
+
+    Info<< "End\n" << endl;
+
+    return 0;
+}
+
+
+// ************************************************************************* //

applications/solvers/cfdemSolverRhoSimple/createFieldRefs.H

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

applications/solvers/cfdemSolverRhoSimple/createFields.H

@@ -0,0 +1,242 @@
+Info<< "Reading thermophysical properties\n" << endl;
+
+    autoPtr<psiThermo> pThermo
+    (
+        psiThermo::New(mesh)
+    );
+    psiThermo& thermo = pThermo();
+    thermo.validate(args.executable(), "h", "e");
+    volScalarField& p = thermo.p();
+
+    Info<< "Reading field rho\n" << endl;
+    volScalarField rho
+    (
+        IOobject
+        (
+            "rho",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        thermo.rho()
+    );
+
+
+    Info<< "Reading field U\n" << endl;
+    volVectorField U
+    (
+        IOobject
+        (
+            "U",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< "\nReading voidfraction field voidfraction = (Vgas/Vparticle)\n" << endl;
+    volScalarField voidfraction
+    (
+        IOobject
+        (
+            "voidfraction",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    volScalarField addSource
+    (
+        IOobject
+        (
+           "addSource",
+           runTime.timeName(),
+           mesh,
+           IOobject::MUST_READ,
+           IOobject::AUTO_WRITE
+        ),
+        mesh
+    );
+
+    Info<< "\nCreating fluid-particle heat flux field\n" << endl;
+    volScalarField Qsource
+    (
+        IOobject
+        (
+            "Qsource",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("zero", dimensionSet(1,-1,-3,0,0,0,0), 0.0)
+    );
+
+    Info<< "\nCreating fluid-particle heat flux coefficient field\n" << endl;
+    volScalarField QCoeff
+    (
+        IOobject
+        (
+            "QCoeff",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("zero", dimensionSet(1,-1,-3,-1,0,0,0), 0.0)
+    );
+
+    Info<< "\nCreating thermal conductivity field\n" << endl;
+    volScalarField thCond
+    (
+        IOobject
+        (
+            "thCond",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("zero", dimensionSet(1,1,-3,-1,0,0,0), 0.0)
+    );
+
+    Info<< "\nCreating heat capacity field\n" << endl;
+    volScalarField Cpv
+    (
+        IOobject
+        (
+            "Cpv",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("zero", dimensionSet(0,2,-2,-1,0,0,0), 0.0)
+    );
+
+    Info<< "\nCreating body force field\n" << endl;
+    volVectorField fOther
+    (
+        IOobject
+        (
+            "fOther",
+            runTime.timeName(),
+            mesh,
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedVector("zero", dimensionSet(1,-2,-2,0,0,0,0), vector::zero)
+    );
+
+    Info<< "Reading/calculating face flux field phi\n" << endl;
+    surfaceScalarField phi
+    (
+        IOobject
+        (
+            "phi",
+            runTime.timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+         ),
+         linearInterpolate(rho*U*voidfraction) & mesh.Sf()
+    );
+
+    dimensionedScalar rhoMax
+    (
+        dimensionedScalar::lookupOrDefault
+        (
+            "rhoMax",
+            simple.dict(),
+            dimDensity,
+            GREAT
+        )
+    );
+
+    dimensionedScalar rhoMin
+    (
+        dimensionedScalar::lookupOrDefault
+        (
+            "rhoMin",
+            simple.dict(),
+            dimDensity,
+            0
+        )
+    );
+
+    Info<< "Creating turbulence model\n" << endl;
+    autoPtr<compressible::turbulenceModel> turbulence
+    (
+        compressible::turbulenceModel::New
+        (
+            rho,
+            U,
+            phi,
+            thermo
+        )
+    );
+
+    label pRefCell = 0;
+    scalar pRefValue = 0.0;
+    setRefCell(p, simple.dict(), pRefCell, pRefValue);
+
+    mesh.setFluxRequired(p.name());
+
+    Info<< "Creating field dpdt\n" << endl;
+    volScalarField dpdt
+    (
+        IOobject
+        (
+            "dpdt",
+            runTime.timeName(),
+            mesh
+        ),
+        mesh,
+        dimensionedScalar("dpdt", p.dimensions()/dimTime, 0)
+    );
+
+    Info<< "Creating field kinetic energy K\n" << endl;
+    volScalarField K("K", 0.5*magSqr(U));
+
+    Info<< "\nReading momentum exchange field Ksl\n" << endl;
+    volScalarField Ksl
+    (
+        IOobject
+        (
+            "Ksl",
+            runTime.timeName(),
+            mesh,
+            IOobject::MUST_READ,
+            IOobject::AUTO_WRITE
+        ),
+        mesh
+        //dimensionedScalar("0", dimensionSet(1, -3, -1, 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
+    );
+
+//===============================

applications/solvers/cfdemSolverRhoSimple/pEqn.H

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

doc/CFDEMcoupling_Manual.html

@@ -215,10 +215,11 @@ listing below of styles within certain commands.
 <TR ALIGN="center"><TD ><A HREF = "dataExchangeModel_twoWayMPI.html">dataExchangeModel_twoWayMPI</A></TD><TD ><A HREF = "dataExchangeModel_twoWayMany2Many.html">dataExchangeModel_twoWayMany2Many</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "forceModel.html">forceModel</A></TD><TD ><A HREF = "forceModel_Archimedes.html">forceModel_Archimedes</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "forceModel_ArchimedesIB.html">forceModel_ArchimedesIB</A></TD><TD ><A HREF = "forceModel_DiFeliceDrag.html">forceModel_DiFeliceDrag</A></TD></TR>
-<TR ALIGN="center"><TD ><A HREF = "forceModel_GidaspowDrag.html">forceModel_GidaspowDrag</A></TD><TD ><A HREF = "forceModel_KochHillDrag.html">forceModel_KochHillDrag</A></TD></TR>
+<TR ALIGN="center"><TD ><A HREF = "forceModel_dSauter.html">forceModel_dSauter</A></TD><TD ><A HREF = "forceModel_GidaspowDrag.html">forceModel_GidaspowDrag</A></TD></TR>
-<TR ALIGN="center"><TD ><A HREF = "forceModel_LaEuScalarTemp.html">forceModel_LaEuScalarTemp</A></TD><TD ><A HREF = "forceModel_MeiLift.html">forceModel_MeiLift</A></TD></TR>
+<TR ALIGN="center"><TD ><A HREF = "forceModel_KochHillDrag.html">forceModel_KochHillDrag</A></TD><TD ><A HREF = "forceModel_LaEuScalarTemp.html">forceModel_LaEuScalarTemp</A></TD></TR>
-<TR ALIGN="center"><TD ><A HREF = "forceModel_SchillerNaumannDrag.html">forceModel_SchillerNaumannDrag</A></TD><TD ><A HREF = "forceModel_ShirgaonkarIB.html">forceModel_ShirgaonkarIB</A></TD></TR>
+<TR ALIGN="center"><TD ><A HREF = "forceModel_MeiLift.html">forceModel_MeiLift</A></TD><TD ><A HREF = "forceModel_SchillerNaumannDrag.html">forceModel_SchillerNaumannDrag</A></TD></TR>
-<TR ALIGN="center"><TD ><A HREF = "forceModel_fieldStore.html">forceModel_fieldStore</A></TD><TD ><A HREF = "forceModel_gradPForce.html">forceModel_gradPForce</A></TD></TR>
+<TR ALIGN="center"><TD ><A HREF = "forceModel_ShirgaonkarIB.html">forceModel_ShirgaonkarIB</A></TD><TD ><A HREF = "forceModel_fieldStore.html">forceModel_fieldStore</A></TD></TR>
+<TR ALIGN="center"><TD ><A HREF = "forceModel_pdCorrelation.html">forceModel_pdCorrelation</A></TD><TD ><A HREF = "forceModel_gradPForce.html">forceModel_gradPForce</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "forceModel_noDrag.html">forceModel_noDrag</A></TD><TD ><A HREF = "forceModel_particleCellVolume.html">forceModel_particleCellVolume</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "forceModel_virtualMassForce.html">forceModel_virtualMassForce</A></TD><TD ><A HREF = "forceModel_viscForce.html">forceModel_viscForce</A></TD></TR>
 <TR ALIGN="center"><TD ><A HREF = "forceSubModel.html">forceSubModel</A></TD><TD ><A HREF = "forceSubModel_ImEx.html">forceSubModel_ImEx</A></TD></TR>

doc/CFDEMcoupling_Manual.txt

@@ -260,6 +260,7 @@ listing below of styles within certain commands.
 "forceModel_Archimedes"_forceModel_Archimedes.html,
 "forceModel_ArchimedesIB"_forceModel_ArchimedesIB.html,
 "forceModel_DiFeliceDrag"_forceModel_DiFeliceDrag.html,
+"forceModel_dSauter"_forceModel_dSauter.html,
 "forceModel_GidaspowDrag"_forceModel_GidaspowDrag.html,
 "forceModel_KochHillDrag"_forceModel_KochHillDrag.html,
 "forceModel_LaEuScalarTemp"_forceModel_LaEuScalarTemp.html,
@@ -267,6 +268,7 @@ listing below of styles within certain commands.
 "forceModel_SchillerNaumannDrag"_forceModel_SchillerNaumannDrag.html,
 "forceModel_ShirgaonkarIB"_forceModel_ShirgaonkarIB.html,
 "forceModel_fieldStore"_forceModel_fieldStore.html,
+"forceModel_pdCorrelation"_forceModel_pdCorrelation.html,
 "forceModel_gradPForce"_forceModel_gradPForce.html,
 "forceModel_noDrag"_forceModel_noDrag.html,
 "forceModel_particleCellVolume"_forceModel_particleCellVolume.html,

doc/forceModel.html

@@ -41,7 +41,7 @@
 </P>
 <P><B>Related commands:</B>
 </P>
-<P><A HREF = "forceModel_Archimedes.html">Archimedes</A>, <A HREF = "forceModel_DiFeliceDrag.html">DiFeliceDrag</A>, <A HREF = "forceModel_gradPForce.html">gradPForce</A>, <A HREF = "forceModel_viscForce.html">viscForce</A>
+<P><A HREF = "forceModel_Archimedes.html">Archimedes</A>, <A HREF = "forceModel_DiFeliceDrag.html">DiFeliceDrag</A>, <A HREF = "forceModel_gradPForce.html">gradPForce</A>, <A HREF = "forceModel_viscForce.html,<A HREF = "forceModel_dSauter.html">viscForce</A>">dSauter</A>
 </P>
 <P>Note: This examples list may be incomplete - please look for other models (forceModel_XY) in this documentation.
 </P>

doc/forceModel.txt

@@ -39,7 +39,7 @@ None.
 
 [Related commands:]
 
-"Archimedes"_forceModel_Archimedes.html, "DiFeliceDrag"_forceModel_DiFeliceDrag.html, "gradPForce"_forceModel_gradPForce.html, "viscForce"_forceModel_viscForce.html
+"Archimedes"_forceModel_Archimedes.html, "DiFeliceDrag"_forceModel_DiFeliceDrag.html, "gradPForce"_forceModel_gradPForce.html, "viscForce"_forceModel_viscForce.html,"dSauter"_forceModel_dSauter.html
 
 Note: This examples list may be incomplete - please look for other models (forceModel_XY) in this documentation.
 

doc/forceModel_dSauter.html

@@ -0,0 +1,47 @@
+<HTML>
+<CENTER><A HREF = "http://www.cfdem.com">CFDEMproject WWW Site</A> - <A HREF = "CFDEMcoupling_Manual.html#comm">CFDEM Commands</A> 
+</CENTER>
+
+
+
+
+<HR>
+
+<H3>forceModel_dSauter command 
+</H3>
+<P><B>Syntax:</B>
+</P>
+<P>Defined in couplingProperties dictionary.
+</P>
+<PRE>forceModels
+(
+    dSauter
+);
+dSauterProps
+{
+    coarseGrainingFactors
+    (
+        X Y Z
+    );
+}; 
+</PRE>
+<UL><LI><I>coarseGrainingFactors</I> = list of coarse graining factors by type, separated by whitespace, optional 
+
+
+</UL>
+<P><B>Description:</B>
+</P>
+<P>This "forceModel" does not influence the particles or the flow - it calculates the Sauter diameter
+</P>
+<CENTER><IMG SRC = "Eqs/d32.png">
+</CENTER>
+<P>.
+</P>
+<P><B>Restrictions:</B>
+</P>
+<P>none.
+</P>
+<P><B>Related commands:</B>
+</P>
+<P><A HREF = "forceModel.html">forceModel</A></P>
+</HTML>

doc/forceModel_dSauter.txt

@@ -0,0 +1,44 @@
+"CFDEMproject WWW Site"_lws - "CFDEM Commands"_lc :c
+
+:link(lws,http://www.cfdem.com)
+:link(lc,CFDEMcoupling_Manual.html#comm)
+
+:line
+
+forceModel_dSauter command :h3
+
+[Syntax:]
+
+Defined in couplingProperties dictionary.
+
+forceModels
+(
+    dSauter
+);
+dSauterProps
+\{
+    coarseGrainingFactors
+    (
+        X Y Z
+    );
+\}; :pre
+
+{coarseGrainingFactors} = list of coarse graining factors by type, separated by whitespace, optional :ulb,l
+
+:ule
+
+[Description:]
+
+This "forceModel" does not influence the particles or the flow - it calculates the Sauter diameter
+
+:c,image(Eqs/d32.png)
+.
+
+
+[Restrictions:]
+
+none.
+
+[Related commands:]
+
+"forceModel"_forceModel.html

doc/forceModel_pdCorrelation.html

@@ -0,0 +1,58 @@
+<HTML>
+<CENTER><A HREF = "http://www.cfdem.com">CFDEMproject WWW Site</A> - <A HREF = "CFDEMcoupling_Manual.html#comm">CFDEM Commands</A> 
+</CENTER>
+
+
+
+
+<HR>
+
+<H3>forceModel_pdCorrelation 
+</H3>
+<P><B>Syntax:</B>
+</P>
+<P>Defined in couplingProperties dictionary.
+</P>
+<PRE>forceModels
+(
+    pdCorrelation
+);
+pdCorrelationProps
+{
+    coarseGrainingFactors
+    (
+        X Y Z
+    );
+    particleDensities
+    (
+        A B C
+    );
+    runOnWriteOnly true;
+}; 
+</PRE>
+<UL><LI><I>coarseGrainingFactors</I> = list of coarse graining factors by type, separated by whitespace, optional 
+
+<LI><I>particleDensities</I> = list of particle densities by type, separated by whitespace, optional 
+
+<LI><I>runOnWriteOnly</I> = switch if this should be executed on write, optional (default: false - execute every coupling step). 
+
+
+</UL>
+<P><B>Description:</B>
+</P>
+<P>This "forceModel" does not influence the particles or the flow - it calculates the particle momentum-diameter correlation 
+</P>
+<CENTER><IMG SRC = "Eqs/pdCorrelation.png">
+</CENTER>
+<P>where delta is the type-specific coarsegraining factor. 
+</P>
+<P>This model is sensitive to additionally pulled particle type info, and can either use the type-specific densities from the dict or those pulled from LIGGGHTS.
+</P>
+<P><B>Restrictions:</B>
+</P>
+<P>none.
+</P>
+<P><B>Related commands:</B>
+</P>
+<P><A HREF = "forceModel.html">forceModel</A></P>
+</HTML>

doc/forceModel_pdCorrelation.txt

@@ -0,0 +1,53 @@
+"CFDEMproject WWW Site"_lws - "CFDEM Commands"_lc :c
+
+:link(lws,http://www.cfdem.com)
+:link(lc,CFDEMcoupling_Manual.html#comm)
+
+:line
+
+forceModel_pdCorrelation :h3
+
+[Syntax:]
+
+Defined in couplingProperties dictionary.
+
+forceModels
+(
+    pdCorrelation
+);
+pdCorrelationProps
+\{
+    coarseGrainingFactors
+    (
+        X Y Z
+    );
+    particleDensities
+    (
+        A B C
+    );
+    runOnWriteOnly true;
+\}; :pre
+
+{coarseGrainingFactors} = list of coarse graining factors by type, separated by whitespace, optional :ulb,l
+{particleDensities} = list of particle densities by type, separated by whitespace, optional :l
+{runOnWriteOnly} = switch if this should be executed on write, optional (default: false - execute every coupling step). :l
+
+
+:ule
+
+[Description:]
+
+This "forceModel" does not influence the particles or the flow - it calculates the particle momentum-diameter correlation 
+
+:c,image(Eqs/pdCorrelation.png)
+where delta is the type-specific coarsegraining factor. 
+
+This model is sensitive to additionally pulled particle type info, and can either use the type-specific densities from the dict or those pulled from LIGGGHTS.
+
+[Restrictions:]
+
+none.
+
+[Related commands:]
+
+"forceModel"_forceModel.html

etc/bashrc

@@ -17,7 +17,7 @@
 #------------------------------------------------------------------------------
 
 export CFDEM_PROJECT=CFDEM
-export CFDEM_VERSION=17.08
+export CFDEM_VERSION=18.03
 
 ################################################################################
 # USER EDITABLE PART: Changes made here may be lost with the next upgrade

etc/cshrc

@@ -15,7 +15,7 @@
 #------------------------------------------------------------------------------
 
 setenv CFDEM_PROJECT CFDEM
-setenv CFDEM_VERSION 17.08
+setenv CFDEM_VERSION 18.03
 
 ################################################################################
 # USER EDITABLE PART: Changes made here may be lost with the next upgrade

etc/solver-list.txt

@@ -6,5 +6,7 @@ rtfmSolverSpecies/dir
 cfdemSolverRhoPimple/dir
 cfdemSolverPisoMS/dir
 cfdemSolverPiso/dir
+cfdemSolverRhoPimple/dir
+cfdemSolverRhoSimple/dir
 cfdemSolverIB/dir
 cfdemSolverPisoScalar/dir

src/lagrangian/cfdemParticle/Make/files

@@ -33,12 +33,13 @@ $(cfdTools)/newGlobal.C
 $(energyModels)/energyModel/energyModel.C
 $(energyModels)/energyModel/newEnergyModel.C
 $(energyModels)/heatTransferGunn/heatTransferGunn.C
-$(energyModels)/heatTransferGunnImplicit/heatTransferGunnImplicit.C
+$(energyModels)/heatTransferGunnPartField/heatTransferGunnPartField.C
 $(energyModels)/reactionHeat/reactionHeat.C
 
 $(thermCondModels)/thermCondModel/thermCondModel.C
 $(thermCondModels)/thermCondModel/newThermCondModel.C
 $(thermCondModels)/SyamlalThermCond/SyamlalThermCond.C
+$(thermCondModels)/ZehnerSchluenderThermCond/ZehnerSchluenderThermCond.C
 $(thermCondModels)/noTherm/noThermCond.C
 
 $(forceModels)/forceModel/forceModel.C
@@ -70,12 +71,14 @@ $(forceModels)/evaluateFluctuations/evaluateFluctuations.C
 $(forceModels)/directedDiffusiveRelaxation/directedDiffusiveRelaxation.C
 $(forceModels)/potentialRelaxation/potentialRelaxation.C
 $(forceModels)/BeetstraDrag/BeetstraDrag.C
+$(forceModels)/BeetstraDragPoly/BeetstraDragPoly.C
 $(forceModels)/dSauter/dSauter.C
 $(forceModels)/Fines/Fines.C
 $(forceModels)/Fines/FinesFields.C
 $(forceModels)/Fines/FanningDynFines.C
 $(forceModels)/Fines/ErgunStatFines.C
 $(forceModels)/granKineticEnergy/granKineticEnergy.C
+$(forceModels)/pdCorrelation/pdCorrelation.C
 
 $(forceModelsMS)/forceModelMS/forceModelMS.C
 $(forceModelsMS)/forceModelMS/newForceModelMS.C

src/lagrangian/cfdemParticle/cfdTools/debugInfo.H

@@ -69,8 +69,8 @@
   else
   {
       meanAlpha_field = 0;
-      meanU_field = vector(0,0,0);
+      meanU_field = vector::zero;
-      meanUs_field = vector(0,0,0);
+      meanUs_field = vector::zero;
   }
   meanUs_array /= numberOfParticles()+SMALL;
   meanR_array /= numberOfParticles()+SMALL;

src/lagrangian/cfdemParticle/cfdTools/versionInfo.H

@@ -34,8 +34,8 @@ Description
 #ifndef versionInfo_H
 #define versionInfo_H
 
-word CFDEMversion="PFM 17.08";
+word CFDEMversion="PFM 18.03";
-word compatibleLIGGGHTSversion="PFM 17.08";
+word compatibleLIGGGHTSversion="PFM 18.03";
 word OFversion="4.x";
 
 Info << "\nCFDEMcoupling version: " << CFDEMversion <<  endl;

src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.C

@@ -83,6 +83,9 @@ cfdemCloud::cfdemCloud
     ignore_(false),
     allowCFDsubTimestep_(true),
     limitDEMForces_(false),
+    getParticleDensities_(couplingProperties_.lookupOrDefault<bool>("getParticleDensities",false)),
+    getParticleEffVolFactors_(couplingProperties_.lookupOrDefault<bool>("getParticleEffVolFactors",false)),
+    getParticleTypes_(couplingProperties_.lookupOrDefault<bool>("getParticleTypes",false)),
     modelType_(couplingProperties_.lookup("modelType")),
     positions_(NULL),
     velocities_(NULL),
@@ -95,6 +98,9 @@ cfdemCloud::cfdemCloud
     radii_(NULL),
     voidfractions_(NULL),
     cellIDs_(NULL),
+    particleDensities_(NULL),
+    particleEffVolFactors_(NULL),
+    particleTypes_(NULL),
     particleWeights_(NULL),
     particleVolumes_(NULL),
     particleV_(NULL),
@@ -129,6 +135,19 @@ cfdemCloud::cfdemCloud
         mesh,
         dimensionedScalar("zero", dimensionSet(0,0,-1,0,0), 0)  // 1/s
     ),
+    particleDensityField_
+    (
+        IOobject
+        (
+            "partRho",
+            mesh.time().timeName(),
+            mesh,
+            IOobject::READ_IF_PRESENT,
+            IOobject::NO_WRITE
+        ),
+        mesh,
+        dimensionedScalar("zero", dimensionSet(1,-3,0,0,0), 0.0)
+    ),
     checkPeriodicCells_(false),
     turbulence_
     (
@@ -362,6 +381,9 @@ cfdemCloud::~cfdemCloud()
     dataExchangeM().destroy(particleV_,1);
     dataExchangeM().destroy(particleConvVel_,3);
     dataExchangeM().destroy(particleFlucVel_,3);
+    if(getParticleDensities_) dataExchangeM().destroy(particleDensities_,1);
+    if(getParticleEffVolFactors_) dataExchangeM().destroy(particleEffVolFactors_,1);
+    if(getParticleTypes_) dataExchangeM().destroy(particleTypes_,1);
 }
 
 // * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
@@ -373,6 +395,10 @@ void cfdemCloud::getDEMdata()
 
     if(impDEMdragAcc_)
         dataExchangeM().getData("dragAcc","vector-atom",fAcc_); // array is used twice - might be necessary to clean it first
+
+    if(getParticleDensities_) dataExchangeM().getData("density","scalar-atom",particleDensities_);
+    if(getParticleEffVolFactors_) dataExchangeM().getData("effvolfactor","scalar-atom",particleEffVolFactors_);
+    if(getParticleTypes_) dataExchangeM().getData("type","scalar-atom",particleTypes_);
 }
 
 void cfdemCloud::giveDEMdata()
@@ -452,6 +478,25 @@ void cfdemCloud::setParticleForceField()
     );
 }
 
+void cfdemCloud::setScalarAverages()
+{
+    if(!getParticleDensities_) return;
+    if(verbose_) Info << "- setScalarAverage" << endl;
+
+    particleDensityField_.primitiveFieldRef() = 0.0;
+    averagingM().resetWeightFields();
+    averagingM().setScalarAverage
+    (
+        particleDensityField_,
+        particleDensities_,
+        particleWeights_,
+        averagingM().UsWeightField(),
+        NULL
+    );
+
+    if(verbose_) Info << "setScalarAverage done." << endl;
+}
+
 void cfdemCloud::setVectorAverages()
 {
     if(verbose_) Info << "- setVectorAverage(Us,velocities_,weights_)" << endl;
@@ -604,7 +649,8 @@ bool cfdemCloud::evolve
             clockM().stop("setvoidFraction");
 
             // set average particles velocity field
-            clockM().start(20,"setVectorAverage");
+            clockM().start(20,"setAverages");
+            setScalarAverages();
             setVectorAverages();
 
 
@@ -618,7 +664,7 @@ bool cfdemCloud::evolve
             if(!treatVoidCellsAsExplicitForce())
                 smoothingM().smoothenReferenceField(averagingM().UsNext());
 
-            clockM().stop("setVectorAverage");
+            clockM().stop("setAverages");
         }
 
         //============================================
@@ -711,6 +757,9 @@ bool cfdemCloud::reAllocArrays()
         dataExchangeM().allocateArray(particleWeights_,0.,voidFractionM().maxCellsPerParticle());
         dataExchangeM().allocateArray(particleVolumes_,0.,voidFractionM().maxCellsPerParticle());
         dataExchangeM().allocateArray(particleV_,0.,1);
+        if(getParticleDensities_) dataExchangeM().allocateArray(particleDensities_,0.,1);
+        if(getParticleEffVolFactors_) dataExchangeM().allocateArray(particleEffVolFactors_,0.,1);
+        if(getParticleTypes_) dataExchangeM().allocateArray(particleTypes_,0,1);
         arraysReallocated_ = true;
         return true;
     }

src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloud.H

@@ -96,6 +96,12 @@ protected:
 
     bool limitDEMForces_;
 
+    bool getParticleDensities_;
+
+    bool getParticleEffVolFactors_;
+
+    bool getParticleTypes_;
+
     scalar maxDEMForce_;
 
     const word modelType_;
@@ -122,6 +128,12 @@ protected:
 
     mutable int **cellIDs_;
 
+    mutable double **particleDensities_;
+
+    mutable double **particleEffVolFactors_;
+
+    mutable int **particleTypes_;
+
     mutable double **particleWeights_;
 
     mutable double **particleVolumes_;
@@ -166,6 +178,8 @@ protected:
 
     mutable volScalarField ddtVoidfraction_;
 
+    mutable volScalarField particleDensityField_;
+
     mutable Switch checkPeriodicCells_;
 
     const turbulenceModel& turbulence_;
@@ -209,6 +223,8 @@ protected:
 
     virtual void setParticleForceField();
 
+    virtual void setScalarAverages();
+
     virtual void setVectorAverages();
 
 public:
@@ -331,18 +347,26 @@ public:
        virtual inline int maxType() {return -1;}
        virtual inline bool multipleTypesDMax() {return false;}
        virtual inline bool multipleTypesDMin() {return false;}
-       virtual inline double ** particleDensity() const {return NULL;}
+
-       virtual inline int ** particleTypes() const {return NULL;}
+       inline bool getParticleDensities() const;
-       virtual label particleType(label index) const {return -1;}
+       virtual inline double ** particleDensities() const;
+       virtual inline double particleDensity(label index) const;
+
+       inline bool getParticleEffVolFactors() const;
+       virtual inline double particleEffVolFactor(label index) const;
+
+       inline bool getParticleTypes() const;
+       virtual inline int ** particleTypes() const;
+       virtual inline label particleType(label index) const;
 
        //access to the particle's rotation and torque data
        virtual inline double ** DEMTorques() const {return NULL;}
        virtual inline double ** omegaArray() const {return NULL;}
-       virtual vector omega(int) const {return vector(0,0,0);}
+       virtual vector omega(int) const {return vector::zero;}
 
        //access to the particles' orientation information
        virtual inline double ** exArray() const {return NULL;}
-       virtual vector ex(int) const {return vector(0,0,0);}
+       virtual vector ex(int) const {return vector::zero;}
 
        //Detector if SRF module is enable or not
        virtual inline bool SRFOn(){return false;}

src/lagrangian/cfdemParticle/cfdemCloud/cfdemCloudI.H

@@ -226,6 +226,58 @@ inline double cfdemCloud::d32(bool recalc)
     return d32_;
 }
 
+inline bool cfdemCloud::getParticleDensities() const
+{
+    return getParticleDensities_;
+}
+
+inline double ** cfdemCloud::particleDensities() const
+{
+    return particleDensities_;
+}
+
+inline double cfdemCloud::particleDensity(label index) const
+{
+    if(!getParticleDensities_) return -1.0;
+    else
+    {
+        return particleDensities_[index][0];
+    }
+}
+
+inline bool cfdemCloud::getParticleEffVolFactors() const
+{
+    return getParticleEffVolFactors_;
+}
+
+inline double cfdemCloud::particleEffVolFactor(label index) const
+{
+    if(!getParticleEffVolFactors_) return -1.0;
+    else
+    {
+        return particleEffVolFactors_[index][0];
+    }
+}
+
+inline bool cfdemCloud::getParticleTypes() const
+{
+    return getParticleTypes_;
+}
+
+inline int ** cfdemCloud::particleTypes() const
+{
+    return particleTypes_;
+}
+
+inline label cfdemCloud::particleType(label index) const
+{
+    if(!getParticleDensities_) return -1;
+    else
+    {
+        return particleTypes_[index][0];
+    }
+}
+
 inline int cfdemCloud::numberOfParticles() const
 {
     return numberOfParticles_;

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

@@ -169,6 +169,12 @@ void cfdemCloudEnergy::postFlow()
         energyModel_[i]().postFlow();
 }
 
+void cfdemCloudEnergy::solve()
+{
+    for (int i=0;i<nrEnergyModels();i++)
+        energyModel_[i]().solve();
+}
+
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

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

@@ -108,6 +108,8 @@ public:
 
        void postFlow();
 
+       void solve();
+
 };
 
 

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

@@ -81,7 +81,7 @@ void averagingModel::undoVectorAverage
 
                     if(weightField[cellI] == weightP)
                     {
-                        fieldNext[cellI] = vector(0,0,0);
+                        fieldNext[cellI] = vector::zero;
                     }else
                     {
                         fieldNext[cellI] = (fieldNext[cellI]*weightField[cellI]-valueVec*weightP)/(weightField[cellI]-weightP);

src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMany2Many/library/Makefile.fedora_fpic

@@ -27,7 +27,7 @@ lib:	$(OBJ)
 	$(ARCHIVE) $(ARFLAGS) $(LIB) $(OBJ)
 
 clean:
-	rm $(LIB) *.o *.d
+	rm -f $(LIB) *.o *.d
 
 # Compilation rules
 

src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMany2Many/library/Makefile.fedora_fpic_debug

@@ -27,7 +27,7 @@ lib:	$(OBJ)
 	$(ARCHIVE) $(ARFLAGS) $(LIB) $(OBJ)
 
 clean:
-	rm $(LIB) *.o *.d
+	rm -f $(LIB) *.o *.d
 
 # Compilation rules
 

src/lagrangian/cfdemParticle/subModels/dataExchangeModel/twoWayMany2Many/library/Makefile.g++

@@ -26,7 +26,7 @@ lib:	$(OBJ)
 	$(ARCHIVE) $(ARFLAGS) $(LIB) $(OBJ)
 
 clean:
-	rm $(LIB) *.o *.d
+	rm -f $(LIB) *.o *.d
 
 # Compilation rules
 

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

@@ -105,6 +105,8 @@ public:
 	
 	virtual void postFlow() {}
 
+	virtual void solve() {}
+	
 	scalar Cp() const;
 
         virtual scalar aveTpart() const

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

@@ -46,7 +46,7 @@ heatTransferGunn::heatTransferGunn
     expNusselt_(propsDict_.lookupOrDefault<bool>("expNusselt",false)),
     interpolation_(propsDict_.lookupOrDefault<bool>("interpolation",false)),
     verbose_(propsDict_.lookupOrDefault<bool>("verbose",false)),
-    NusseltFieldName_(propsDict_.lookupOrDefault<word>("NusseltFieldName","NuField")),
+    implicit_(propsDict_.lookupOrDefault<bool>("implicit",true)),
     QPartFluidName_(propsDict_.lookupOrDefault<word>("QPartFluidName","QPartFluid")),
     QPartFluid_
     (   IOobject
@@ -61,13 +61,26 @@ heatTransferGunn::heatTransferGunn
         dimensionedScalar("zero", dimensionSet(1,-1,-3,0,0,0,0), 0.0),
 	"zeroGradient"
     ),
+    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)
+    ),
     partTempField_
     (   IOobject
         (
-            "particleTemp",
+            "partTemp",
             sm.mesh().time().timeName(),
             sm.mesh(),
-            IOobject::NO_READ,
+            IOobject::READ_IF_PRESENT,
             IOobject::NO_WRITE
         ),
         sm.mesh(),
@@ -100,6 +113,7 @@ heatTransferGunn::heatTransferGunn
         dimensionedScalar("zero", dimensionSet(0,0,0,0,0,0,0), 0.0),
 	"zeroGradient"
     ),
+    NusseltFieldName_(propsDict_.lookupOrDefault<word>("NusseltFieldName","NuField")),
     NuField_
     (   IOobject
         (
@@ -130,6 +144,7 @@ heatTransferGunn::heatTransferGunn
     partTemp_(NULL),
     partHeatFluxName_(propsDict_.lookup("partHeatFluxName")),
     partHeatFlux_(NULL),
+    partHeatFluxCoeff_(NULL),
     partRe_(NULL),
     partNu_(NULL)
 {
@@ -155,6 +170,11 @@ heatTransferGunn::heatTransferGunn
         Info <<  "Particle temperature field activated." << endl;
     }
 
+    if (!implicit_)
+    {
+        QPartFluidCoeff_.writeOpt() = IOobject::NO_WRITE;
+    }
+
     if (verbose_)
     {
         ReField_.writeOpt() = IOobject::AUTO_WRITE;
@@ -184,6 +204,10 @@ heatTransferGunn::~heatTransferGunn()
     particleCloud_.dataExchangeM().destroy(partHeatFlux_,1);
     particleCloud_.dataExchangeM().destroy(partRe_,1);
     particleCloud_.dataExchangeM().destroy(partNu_,1);
+    if (implicit_)
+    {
+        particleCloud_.dataExchangeM().destroy(partHeatFluxCoeff_,1);
+    }
 }
 
 // * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
@@ -193,6 +217,10 @@ void heatTransferGunn::allocateMyArrays() const
     double initVal=0.0;
     particleCloud_.dataExchangeM().allocateArray(partTemp_,initVal,1);  // field/initVal/with/lenghtFromLigghts
     particleCloud_.dataExchangeM().allocateArray(partHeatFlux_,initVal,1);
+    if(implicit_)
+    {
+        particleCloud_.dataExchangeM().allocateArray(partHeatFluxCoeff_,initVal,1);
+    }
 
     if(verbose_)
     {
@@ -201,7 +229,68 @@ void heatTransferGunn::allocateMyArrays() const
     }
 }
 
+
+void heatTransferGunn::partTempField()
+{
+    partTempField_.primitiveFieldRef() = 0.0;
+    particleCloud_.averagingM().resetWeightFields();
+    particleCloud_.averagingM().setScalarAverage
+    (
+        partTempField_,
+        partTemp_,
+        particleCloud_.particleWeights(),
+        particleCloud_.averagingM().UsWeightField(),
+        NULL
+    );
+
+    dimensionedScalar aveTemp("aveTemp",dimensionSet(0,0,0,1,0,0,0), partTempAve_);
+    partRelTempField_ = (partTempField_ - aveTemp) / (aveTemp - partRefTemp_);
+}
+
+
+scalar heatTransferGunn::Nusselt(scalar voidfraction, scalar Rep, scalar Pr) const
+{
+    return (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);
+}
+
+void heatTransferGunn::giveData()
+{
+    Info << "total convective particle-fluid heat flux [W] (Eulerian) = " << gSum(QPartFluid_*1.0*QPartFluid_.mesh().V()) << endl;
+    particleCloud_.clockM().start(30,"giveDEM_Tdata");
+    particleCloud_.dataExchangeM().giveData(partHeatFluxName_,"scalar-atom", partHeatFlux_);
+    particleCloud_.clockM().stop("giveDEM_Tdata");
+}
+
+void heatTransferGunn::heatFlux(label index, scalar h, scalar As, scalar Tfluid)
+{
+    scalar hAs = h * As;
+    partHeatFlux_[index][0] = - hAs * partTemp_[index][0];
+    if(!implicit_)
+    {
+        partHeatFlux_[index][0] += hAs * Tfluid;
+    }
+    else
+    {
+        partHeatFluxCoeff_[index][0] = hAs;
+    }
+}
+
 // * * * * * * * * * * * * * * * * Member Fct  * * * * * * * * * * * * * * * //
+void heatTransferGunn::addEnergyContribution(volScalarField& Qsource) const
+{
+    Qsource += QPartFluid_;
+}
+
+void heatTransferGunn::addEnergyCoefficient(volScalarField& Qsource) const
+{
+    if(implicit_)
+    {
+        Qsource += QPartFluidCoeff_;
+    }
+}
 
 void heatTransferGunn::calcEnergyContribution()
 {
@@ -292,7 +381,6 @@ void heatTransferGunn::calcEnergyContribution()
 
                 // calc convective heat flux [W]
                 heatFlux(index, h, As, Tfluid);
-                heatFluxCoeff(index, h, As);
 
                 if(verbose_)
                 {
@@ -326,6 +414,16 @@ void heatTransferGunn::calcEnergyContribution()
     
     if(calcPartTempField_) partTempField();
 
+    // gather particle temperature sums and obtain average
+    if(calcPartTempAve_)
+    {
+        reduce(Tsum, sumOp<scalar>());
+        partTempAve_ = Tsum / particleCloud_.numberOfParticles();
+        Info << "mean particle temperature = " << partTempAve_ << endl;
+    }
+
+    if(calcPartTempField_) partTempField();
+
     particleCloud_.averagingM().setScalarSum
     (
         QPartFluid_,
@@ -336,6 +434,21 @@ void heatTransferGunn::calcEnergyContribution()
 
     QPartFluid_.primitiveFieldRef() /= -QPartFluid_.mesh().V();
 
+    if(implicit_)
+    {
+        QPartFluidCoeff_.primitiveFieldRef() = 0.0;
+
+        particleCloud_.averagingM().setScalarSum
+        (
+            QPartFluidCoeff_,
+            partHeatFluxCoeff_,
+            particleCloud_.particleWeights(),
+            NULL
+        );
+
+        QPartFluidCoeff_.primitiveFieldRef() /= -QPartFluidCoeff_.mesh().V();
+    }
+
     if(verbose_)
     {
         ReField_.primitiveFieldRef() = 0.0;
@@ -360,7 +473,9 @@ void heatTransferGunn::calcEnergyContribution()
         );
     }
 
-    // limit source term
+    // limit source term in explicit treatment
+    if(!implicit_)
+    {
         forAll(QPartFluid_,cellI)
         {
             scalar EuFieldInCell = QPartFluid_[cellI];
@@ -371,71 +486,47 @@ void heatTransferGunn::calcEnergyContribution()
                  QPartFluid_[cellI] = sign(EuFieldInCell) * maxSource_;
             }
         }
+    }
 
     QPartFluid_.correctBoundaryConditions();
-
-    giveData(0);
-
 }
 
-void heatTransferGunn::addEnergyContribution(volScalarField& Qsource) const
+void heatTransferGunn::postFlow()
 {
-    Qsource += QPartFluid_;
+    if(implicit_)
-}
-
-scalar heatTransferGunn::Nusselt(scalar voidfraction, scalar Rep, scalar Pr) const
-{
-    return (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);
-}
-
-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;
+        label cellI;
+        scalar Tfluid(0.0);
+        scalar Tpart(0.0);
+        interpolationCellPoint<scalar> TInterpolator_(tempField_);
 
-        particleCloud_.clockM().start(30,"giveDEM_Tdata");
+        for(int index = 0;index < particleCloud_.numberOfParticles(); ++index)
-        particleCloud_.dataExchangeM().giveData(partHeatFluxName_,"scalar-atom", partHeatFlux_);
+        {
-        particleCloud_.clockM().stop("giveDEM_Tdata");
+                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();
 }
 
 scalar heatTransferGunn::aveTpart() const
 {
     return partTempAve_;
 }
-
-void heatTransferGunn::partTempField()
-{
-    partTempField_.primitiveFieldRef() = 0.0;
-    particleCloud_.averagingM().resetWeightFields();
-    particleCloud_.averagingM().setScalarAverage
-    (
-        partTempField_,
-        partTemp_,
-        particleCloud_.particleWeights(),
-        particleCloud_.averagingM().UsWeightField(),
-        NULL
-    );
-
-    //volScalarField sumTp (particleCloud_.averagingM().UsWeightField() * partTempField_);
-    // dimensionedScalar aveTemp("aveTemp",dimensionSet(0,0,0,1,0,0,0), gSum(sumTp) / particleCloud_.numberOfParticles());
-    dimensionedScalar aveTemp("aveTemp",dimensionSet(0,0,0,1,0,0,0), partTempAve_);
-    partRelTempField_ = (partTempField_ - aveTemp) / (aveTemp - partRefTemp_);
-}
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

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

@@ -51,18 +51,24 @@ protected:
 
     bool verbose_;
 
-    word NusseltFieldName_;
+    bool implicit_;
 
     word QPartFluidName_;
 
     volScalarField QPartFluid_;
 
-    volScalarField partTempField_;
+    word QPartFluidCoeffName_;
+
+    volScalarField QPartFluidCoeff_;
+
+    mutable volScalarField partTempField_;
 
     volScalarField partRelTempField_;
 
     volScalarField ReField_;
 
+    word NusseltFieldName_;
+
     volScalarField NuField_;
 
     dimensionedScalar partRefTemp_;
@@ -93,11 +99,13 @@ protected:
 
     word partTempName_;
 
-    mutable double **partTemp_;                          // Lagrangian array
+    mutable double **partTemp_;
 
     word partHeatFluxName_;
 
-    mutable double **partHeatFlux_;              // Lagrangian array
+    mutable double **partHeatFlux_;
+
+    mutable double **partHeatFluxCoeff_;
 
     mutable double **partRe_;
 
@@ -109,12 +117,10 @@ protected:
 
     scalar Nusselt(scalar, scalar, scalar) const;
 
-    virtual void giveData(int);
+    virtual void giveData();
 
     virtual void heatFlux(label, scalar, scalar, scalar);
 
-    virtual void heatFluxCoeff(label, scalar, scalar);
-       
 public:
 
     //- Runtime type information
@@ -139,10 +145,12 @@ public:
 
         void addEnergyContribution(volScalarField&) const;
 
-	void addEnergyCoefficient(volScalarField&) const {}
+	void addEnergyCoefficient(volScalarField&) const;
 
 	void calcEnergyContribution();
 
+	void postFlow();
+
         scalar aveTpart() const;
 };
 

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

@@ -0,0 +1,244 @@
+/*---------------------------------------------------------------------------*\
+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 "heatTransferGunnPartField.H"
+#include "addToRunTimeSelectionTable.H"
+#include "thermCondModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(heatTransferGunnPartField, 0);
+
+addToRunTimeSelectionTable(energyModel, heatTransferGunnPartField, dictionary);
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from components
+heatTransferGunnPartField::heatTransferGunnPartField
+(
+    const dictionary& dict,
+    cfdemCloudEnergy& sm
+)
+:
+    heatTransferGunn(dict,sm),
+    partCpField_
+    (   
+        IOobject
+        (
+            "partCp",
+            sm.mesh().time().timeName(),
+            sm.mesh(),
+            IOobject::READ_IF_PRESENT,
+            IOobject::NO_WRITE
+        ),
+        sm.mesh(),
+        dimensionedScalar("zero", dimensionSet(0,2,-2,-1,0,0,0), 0.0),
+        "zeroGradient"
+    ),
+    partRhoField_(sm.mesh().lookupObject<volScalarField>("partRho")),
+    typeCp_(propsDict_.lookupOrDefault<scalarList>("specificHeatCapacities",scalarList(1,-1.0))),
+    partCp_(NULL),
+    pTMax_(dimensionedScalar("pTMax",dimensionSet(0,0,0,1,0,0,0), -1.0)),
+    pTMin_(dimensionedScalar("pTMin",dimensionSet(0,0,0,1,0,0,0), -1.0)),
+    thermCondModel_
+    (
+        thermCondModel::New
+        (
+            propsDict_,
+            sm
+        )
+    ),
+    fvOptions(fv::options::New(sm.mesh()))
+{
+    if (!implicit_)
+    {
+        FatalError << "heatTransferGunnPartField requires implicit heat transfer treatment." << abort(FatalError);
+    }
+
+    if (typeCp_[0] < 0.0)
+    {
+        FatalError << "heatTransferGunnPartField: provide list of specific heat capacities." << abort(FatalError);
+    }
+    
+    if (propsDict_.found("pTMax"))
+    {
+        pTMax_.value()=scalar(readScalar(propsDict_.lookup("pTMax"))); 
+    }
+    
+    if (propsDict_.found("pTMin"))
+    {
+        pTMin_.value()=scalar(readScalar(propsDict_.lookup("pTMin"))); 
+    }
+
+    partTempField_.writeOpt() = IOobject::AUTO_WRITE;
+
+    allocateMyArrays();
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+heatTransferGunnPartField::~heatTransferGunnPartField()
+{
+    particleCloud_.dataExchangeM().destroy(partCp_,1);
+}
+
+// * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
+
+void heatTransferGunnPartField::allocateMyArrays() const
+{
+    double initVal=0.0;
+    particleCloud_.dataExchangeM().allocateArray(partCp_,initVal,1);
+}
+// * * * * * * * * * * * * * * * * Member Fct  * * * * * * * * * * * * * * * //
+void heatTransferGunnPartField::calcEnergyContribution()
+{
+    allocateMyArrays();
+    heatTransferGunn::calcEnergyContribution();
+
+    // if heat sources in particles present, pull them here
+
+    // loop over all particles to fill partCp_ based on type
+    label cellI=0;
+    label partType = 0;
+    for(int index = 0;index < particleCloud_.numberOfParticles(); ++index)
+    {
+            cellI = particleCloud_.cellIDs()[index][0];
+            if(cellI >= 0)
+            {
+                partType = particleCloud_.particleType(index);
+                // LIGGGGHTS counts types 1, 2, ..., C++ array starts at 0
+                partCp_[index][0] = typeCp_[partType - 1];
+            }
+    }
+
+    partCpField_.primitiveFieldRef() = 0.0;
+    particleCloud_.averagingM().resetWeightFields();
+    particleCloud_.averagingM().setScalarAverage
+    (
+        partCpField_,
+        partCp_,
+        particleCloud_.particleWeights(),
+        particleCloud_.averagingM().UsWeightField(),
+        NULL
+    );
+
+}
+
+void heatTransferGunnPartField::addEnergyContribution(volScalarField& Qsource) const
+{
+    Qsource -= QPartFluidCoeff_*partTempField_;
+}
+
+void heatTransferGunnPartField::giveData()
+{
+    particleCloud_.dataExchangeM().giveData(partTempName_,"scalar-atom", partTemp_);
+}
+
+void heatTransferGunnPartField::postFlow()
+{
+    label cellI;
+    scalar Tpart(0.0);
+    interpolationCellPoint<scalar> partTInterpolator_(partTempField_);
+
+    particleCloud_.dataExchangeM().allocateArray(partTemp_,0.0,1);
+
+    for(int index = 0;index < particleCloud_.numberOfParticles(); ++index)
+    {
+            cellI = particleCloud_.cellIDs()[index][0];
+            if(cellI >= 0)
+            {
+                if(interpolation_)
+                {
+                    vector position = particleCloud_.position(index);
+                    Tpart = partTInterpolator_.interpolate(position,cellI);
+                }
+                else
+                {
+                    Tpart = partTempField_[cellI];
+                }
+                partTemp_[index][0] = Tpart;
+            }
+    }
+
+    giveData();
+}
+
+void heatTransferGunnPartField::solve()
+{
+    // for some weird reason, the particle-fluid heat transfer fields were defined with a negative sign
+
+    volScalarField alphaP = 1.0 - voidfraction_;
+    volScalarField correctedQPartFluidCoeff(QPartFluidCoeff_);
+    // no heattransfer in empty cells -- for numerical stability, add small amount
+    forAll(correctedQPartFluidCoeff,cellI)
+    {
+        if (-correctedQPartFluidCoeff[cellI] < SMALL) correctedQPartFluidCoeff[cellI] = -SMALL;
+    }
+
+    volScalarField Qsource = correctedQPartFluidCoeff*tempField_;
+    volScalarField partCpEff = alphaP*partRhoField_*partCpField_;
+    volScalarField thCondEff = alphaP*thermCondModel_().thermCond();
+//    thCondEff.correctBoundaryConditions();
+
+
+
+    fvScalarMatrix partTEqn
+    (
+     //   fvm::ddt(partCpEff, partTempField_)
+     // + Qsource
+        Qsource
+      - fvm::Sp(correctedQPartFluidCoeff, partTempField_)
+      - fvm::laplacian(thCondEff,partTempField_)
+     ==
+        fvOptions(partCpEff, partTempField_)
+    );
+  // if transient add time derivative - need particle density and specific heat fields
+  // if sources activated add sources
+  // if convection activated add convection
+
+    partTEqn.relax();
+
+    fvOptions.constrain(partTEqn);
+
+    partTEqn.solve();
+
+    partTempField_.relax();
+
+    fvOptions.correct(partTempField_);
+
+    if (pTMax_.value() > 0.0) partTempField_ = min(partTempField_, pTMax_);
+    if (pTMin_.value() > 0.0) partTempField_ = max(partTempField_, pTMin_);
+
+
+    Info<< "partT max/min : " << max(partTempField_).value() << " " << min(partTempField_).value() << endl;
+
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

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

@@ -0,0 +1,108 @@
+/*---------------------------------------------------------------------------*\
+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 heatTransferGunnPartField_H
+#define heatTransferGunnPartField_H
+
+#include "fvCFD.H"
+#include "cfdemCloudEnergy.H"
+#include "heatTransferGunn.H"
+#include "fvOptions.H"
+#include "scalarList.H"
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+class thermCondModel;
+/*---------------------------------------------------------------------------*\
+                           Class heatTransferGunnPartField Declaration
+\*---------------------------------------------------------------------------*/
+
+class heatTransferGunnPartField
+:
+    public heatTransferGunn
+{
+private:
+
+    volScalarField partCpField_;
+
+    const volScalarField& partRhoField_;
+
+    scalarList typeCp_;
+
+    mutable double **partCp_;
+
+    dimensionedScalar pTMax_;
+
+    dimensionedScalar pTMin_;
+
+    autoPtr<thermCondModel> thermCondModel_;
+
+    fv::options& fvOptions;
+
+    void allocateMyArrays() const;
+
+    void giveData();
+
+public:
+
+    //- Runtime type information
+    TypeName("heatTransferGunnPartField");
+
+    // Constructors
+
+        //- Construct from components
+        heatTransferGunnPartField
+        (
+            const dictionary& dict,
+            cfdemCloudEnergy& sm
+        );
+
+
+    // Destructor
+
+        virtual ~heatTransferGunnPartField();
+
+
+    // Member Functions
+        void addEnergyContribution(volScalarField&) const;
+
+	void calcEnergyContribution();
+
+	void postFlow();
+
+        void solve();
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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

@@ -12,6 +12,7 @@ License
     along with this code.  If not, see <http://www.gnu.org/licenses/>.
 
     Copyright (C) 2015- Thomas Lichtenegger, JKU Linz, Austria
+    Contributing author: Paul Kieckhefen, TU Hamburg, Germany
 
 \*---------------------------------------------------------------------------*/
 
@@ -49,14 +50,26 @@ BeetstraDrag::BeetstraDrag
 :
     forceModel(dict,sm),
     propsDict_(dict.subDict(typeName + "Props")),
+    multiTypes_(false),
     velFieldName_(propsDict_.lookup("velFieldName")),
     U_(sm.mesh().lookupObject<volVectorField> (velFieldName_)),
     voidfractionFieldName_(propsDict_.lookup("voidfractionFieldName")),
     voidfraction_(sm.mesh().lookupObject<volScalarField> (voidfractionFieldName_)),
+    minVoidfraction_(propsDict_.lookupOrDefault<scalar>("minVoidfraction",0.1)),
     UsFieldName_(propsDict_.lookup("granVelFieldName")),
     UsField_(sm.mesh().lookupObject<volVectorField> (UsFieldName_)),
     scaleDia_(1.),
-    scaleDrag_(1.)
+    typeCG_(propsDict_.lookupOrDefault<scalarList>("coarseGrainingFactors",scalarList(1,1.0))),
+    scaleDrag_(1.),
+    rhoP_(0.),
+    rho_(0.),
+    Lc2_(0.),
+    dPrim_(0.),
+    nuf_(0.),
+    g_(9.81),
+    k_(0.05),
+    useGC_(false),
+    usePC_(false)
 {
     //Append the field names to be probed
     particleCloud_.probeM().initialize(typeName, typeName+".logDat");
@@ -78,10 +91,42 @@ BeetstraDrag::BeetstraDrag
     forceSubM(0).readSwitches();
 
     particleCloud_.checkCG(true);
-    if (propsDict_.found("scale"))
+    if (propsDict_.found("scale") && typeCG_.size()==1)
+    {
+        // if "scale" is specified and there's only one single type, use "scale"
         scaleDia_=scalar(readScalar(propsDict_.lookup("scale")));
+        typeCG_[0] = scaleDia_;
+    }
     if (propsDict_.found("scaleDrag"))
+    {
         scaleDrag_=scalar(readScalar(propsDict_.lookup("scaleDrag")));
+    }
+
+    if (typeCG_.size()>1) multiTypes_ = true;
+
+    if (propsDict_.found("useFilteredDragModel"))
+    {
+        useGC_ = true;
+        g_=propsDict_.lookupOrDefault<scalar>("g", 9.81);
+        dPrim_=scalar(readScalar(propsDict_.lookup("dPrim")));
+        rhoP_=scalar(readScalar(propsDict_.lookup("rhoP")));
+        rho_=scalar(readScalar(propsDict_.lookup("rho")));
+        nuf_=scalar(readScalar(propsDict_.lookup("nuf")));
+	    scalar ut = terminalVelocity(1., dPrim_, nuf_, rho_, rhoP_, g_);
+	    scalar Frp = ut*ut/g_/dPrim_;
+        Lc2_ = ut*ut/g_*pow(Frp, -.6666667); // n is hardcoded as -2/3
+        Info << "using grid coarsening correction with Lc2 = " << Lc2_ << " and ut = " << ut << " and Frp = " << Frp<< endl;
+
+        if (propsDict_.found("useParcelSizeDependentFilteredDrag"))
+        {
+            usePC_ = true;
+            if (propsDict_.found("k"))
+                k_=scalar(readScalar(propsDict_.lookup("k")));
+            Info << "using particle coarsening correction with k = " << k_ << endl;
+        }
+    }
+
+
 
 }
 
@@ -96,9 +141,9 @@ BeetstraDrag::~BeetstraDrag()
 
 void BeetstraDrag::setForce() const
 {
-    if (scaleDia_ > 1)
+    if (typeCG_.size()>1 || typeCG_[0] > 1)
     {
-        Info << "Beetstra using scale = " << scaleDia_ << endl;
+        Info << "Beetstra using scale = " << typeCG_ << endl;
     }
     else if (particleCloud_.cg() > 1)
     {
@@ -119,12 +164,18 @@ void BeetstraDrag::setForce() const
     vector Ur(0,0,0);
     scalar ds(0);
     scalar ds_scaled(0);
-    scalar scaleDia3 = scaleDia_*scaleDia_*scaleDia_;
+    scalar dSauter(0);
+    scalar scaleDia3 = typeCG_[0]*typeCG_[0]*typeCG_[0];
     scalar nuf(0);
     scalar rho(0);
     scalar magUr(0);
     scalar Rep(0);
     scalar localPhiP(0);
+    scalar GCcorr(1.);
+    scalar PCcorr(1.);
+
+    scalar cg = typeCG_[0];
+    label partType = 1;
 
     vector dragExplicit(0,0,0);
     scalar dragCoefficient(0);
@@ -137,52 +188,79 @@ void BeetstraDrag::setForce() const
     for(int index = 0; index < particleCloud_.numberOfParticles(); ++index)
     {
             cellI = particleCloud_.cellIDs()[index][0];
-            drag = vector(0,0,0);
+            drag = vector::zero;
-            dragExplicit = vector(0,0,0);
+            dragExplicit = vector::zero;
-            Ufluid =vector(0,0,0);
+            Ufluid = vector::zero;
             voidfraction=0;
             dragCoefficient = 0;
 
             if (cellI > -1) // particle found
             {
 
-                if( forceSubM(0).interpolation() )
+                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 > 1.00) voidfraction = 1.0;
-                    if(voidfraction<0.10) voidfraction = 0.10;
+                    if (voidfraction < minVoidfraction_) voidfraction = minVoidfraction_;
                 }
                 else
                 {
                     voidfraction = voidfraction_[cellI];
                     Ufluid = U_[cellI];
                 }
+                // in case a fines phase is present, void fraction needs to be adapted
+                adaptVoidfraction(voidfraction, cellI);
+
+                if (multiTypes_)
+                {
+                    partType = particleCloud_.particleType(index);
+                    cg = typeCG_[partType - 1];
+                    scaleDia3 = cg*cg*cg;
+                }
 
                 Us = particleCloud_.velocity(index);
                 Ur = Ufluid-Us;
                 magUr = mag(Ur);
                 ds = 2*particleCloud_.radius(index);
-                ds_scaled = ds/scaleDia_;
+                ds_scaled = ds/cg;
+                dSauter = meanSauterDiameter(ds_scaled, cellI);
+
                 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;
+                Rep=dSauter*voidfraction*magUr/nuf + SMALL;
-                dragCoefficient = 10.0*localPhiP/(voidfraction*voidfraction) +
+
-                                  voidfraction*voidfraction*(1.0+1.5*Foam::sqrt(localPhiP)) +
+                dragCoefficient = F(voidfraction, Rep)
-                                  0.413*Rep/(24*voidfraction*voidfraction)*(1.0/voidfraction+3*voidfraction*localPhiP+8.4*Foam::pow(Rep,-0.343))/
+                                   *3*M_PI*nuf*rho*voidfraction
-                                  (1+Foam::pow(10,3*localPhiP)*Foam::pow(Rep,-0.5*(1+4*localPhiP)));
+                                   *effDiameter(ds_scaled, cellI, index)
+                                   *scaleDia3*scaleDrag_;
+
+                // calculate filtering corrections
+                if (useGC_)
+                {
+                    GCcorr = 1.-h(localPhiP)
+                                /(  a(localPhiP)
+                                    *Lc2_
+                                    /std::cbrt(U_.mesh().V()[cellI])
+                                  + 1.
+                                 );
+                    if (usePC_)
+                    {
+                        PCcorr = Foam::exp(k_*(1.-cg));
+                    }
+                }
+
+                // apply filtering corrections
+                dragCoefficient *= GCcorr*PCcorr;
 
-                // calc particle's drag
-                dragCoefficient *= 3*M_PI*ds_scaled*nuf*rho*voidfraction*scaleDia3*scaleDrag_;
                 if (modelType_=="B")
                     dragCoefficient /= voidfraction;
 
@@ -198,11 +276,13 @@ void BeetstraDrag::setForce() const
                     Pout << "Us = " << Us << endl;
                     Pout << "Ur = " << Ur << endl;
                     Pout << "ds = " << ds << endl;
-                    Pout << "ds/scale = " << ds/scaleDia_ << endl;
+                    Pout << "ds/scale = " << ds/cg << endl;
                     Pout << "rho = " << rho << endl;
                     Pout << "nuf = " << nuf << endl;
                     Pout << "voidfraction = " << voidfraction << endl;
                     Pout << "Rep = " << Rep << endl;
+                    Pout << "GCcorr = " << GCcorr << endl;
+		            Pout << "PCcorr = " << PCcorr << endl;
                     Pout << "drag = " << drag << endl;
                 }
 
@@ -223,7 +303,149 @@ void BeetstraDrag::setForce() const
     }
 }
 
+/*********************************************************
+ * "Drag Force of Intermediate Reynolds Number Flow Past *
+ *  Mono- and Bidisperse Arrays of Spheres", eq. 16      *
+ *  R Beetstra, M. A. van der Hoef, JAM Kuipers          *
+ *  AIChE Journal 53(2) (2007)                           *
+ *********************************************************/
+double BeetstraDrag::F(double voidfraction, double Rep) const
+{
+    double localPhiP = max(SMALL,min(1.-SMALL,1.-voidfraction));
+    return   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))
+                         );
 
+}
+ 
+
+/*********************************************************
+ * "A drag model for filtered Euler-Lagange simulations  *
+ *  of clustered gas-particle suspension",               *
+ *  S. Radl, S. Sundaresan,                              *
+ *  Chemical Engineering Science 117 (2014).             *
+ *********************************************************/
+double BeetstraDrag::a(double phiP) const
+{
+    double a0m = 0.;
+    double a1m = 0.;
+    double a2m = 0.;
+    double a3m = 0.;
+    double phipam = 0.;
+    if      (phiP < 0.016)
+    {
+      a0m = 21.51;
+    }
+    else if (phiP < 0.100)
+    {
+      a0m =  1.96; a1m =  29.40; a2m =  164.91; a3m = -1923.;
+    }
+    else if (phiP < 0.180)
+    {
+      a0m =  4.63; a1m =   4.68; a2m = -412.04; a3m =  2254.; phipam = 0.10;
+    }
+    else if (phiP < 0.250)
+    {
+      a0m =  3.52; a1m = -17.99; a2m =  128.80; a3m =  -603.; phipam = 0.18;
+    }
+    else if (phiP < 0.400)
+    {
+      a0m =  2.68; a1m =  -8.20; a2m =    2.18; a3m = 112.33; phipam = 0.25;
+    }
+    else
+    {
+      a0m =  1.79;
+    }
+    const double phiP_minus_phipam = phiP-phipam;
+    return a0m + a1m*phiP_minus_phipam + a2m*phiP_minus_phipam*phiP_minus_phipam
+           + a3m*phiP_minus_phipam*phiP_minus_phipam*phiP_minus_phipam;
+}
+
+double BeetstraDrag::h(double phiP) const
+{
+  double h0m = 0.;
+  double h1m = 0.;
+  double h2m = 0.;
+  double h3m = 0.;
+  double phiphm = 0.;
+
+  if      (phiP < 0.03)
+  {
+                 h1m = 7.97;
+  }
+  else if (phiP < 0.08)
+  {
+    h0m = 0.239; h1m =  4.640; h2m =   -4.410; h3m =  253.630; phiphm = 0.03;
+  }
+  else if (phiP < 0.12)
+  {
+    h0m = 0.492; h1m =  6.100; h2m =   33.630; h3m = -789.600; phiphm = 0.08;
+  }
+  else if (phiP < 0.18)
+  {
+    h0m = 0.739; h1m =  5.010; h2m =  -61.100; h3m =  310.800; phiphm = 0.12;
+  }
+  else if (phiP < 0.34)
+  {
+    h0m = 0.887; h1m =  1.030; h2m =   -5.170; h3m =    5.990; phiphm = 0.18;
+  }
+  else if (phiP < 0.48)
+  {
+    h0m = 0.943; h1m = -0.170; h2m =   -2.290; h3m =   -9.120; phiphm = 0.34;
+  }
+  else if (phiP < 0.55)
+  {
+    h0m = 0.850; h1m = -1.350; h2m =   -6.130; h3m = -132.600; phiphm = 0.48;
+  }
+  else
+  {
+    h0m = 0.680; h1m = -2.340; h2m = -225.200;                 phiphm = 0.55;
+  }
+  const double phiP_minus_phiphm = phiP-phiphm;
+  return h0m + h1m*phiP_minus_phiphm + h2m*phiP_minus_phiphm*phiP_minus_phiphm
+         + h3m*phiP_minus_phiphm*phiP_minus_phiphm*phiP_minus_phiphm;
+}
+
+double BeetstraDrag::terminalVelocity(double voidfraction, double dp, double nuf, double rhof, double rhop, double g) const
+{
+    scalar u0(dp*dp*fabs(rhof-rhop)*g/18./rhof/nuf);
+    scalar Re(u0*dp/nuf);
+    scalar res(1.);
+    scalar u(u0);
+    scalar Fi(0);
+    scalar CdSt(0);
+    Info << "uo: " << u0<<endl;
+    int i = 0;
+    while ((res > 1.e-6) && (i<100))
+    {
+        Info << "Iteration " << i;
+        u0 = u;
+        Info << ", u0 = " << u0;
+        CdSt = 24/Re;
+        Info << ", CdSt = " << CdSt;
+        Fi = F(voidfraction, Re);
+        Info << ", F = ";
+        u = sqrt(1.333333333*fabs(rhof-rhop)*g*dp
+                  /(CdSt*voidfraction*Fi*rhof)
+            );
+        Info << ", u = " << u;
+        Re = fabs(u)*dp/nuf*voidfraction;
+        res = fabs((u-u0)/u);
+        Info << "Res: " << res << endl;
+        i++;
+    }
+    if (res >1.e-6)
+        FatalError << "Terminal velocity calculation diverged!" << endl;
+
+    return u;
+}
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam

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

@@ -40,9 +40,11 @@ class BeetstraDrag
 :
     public forceModel
 {
-private:
+protected:
     dictionary propsDict_;
 
+    bool multiTypes_;
+
     word velFieldName_;
 
     const volVectorField& U_;
@@ -51,14 +53,52 @@ private:
 
     const volScalarField& voidfraction_;
 
+    const scalar minVoidfraction_;
+
     word UsFieldName_;
 
     const volVectorField& UsField_;
 
     mutable scalar scaleDia_;
 
+    scalarList typeCG_;
+
     mutable scalar scaleDrag_;
 
+    mutable scalar rhoP_;
+
+    mutable scalar rho_;
+
+    mutable scalar Lc2_;
+
+    mutable scalar dPrim_;
+
+    mutable scalar nuf_;
+
+    mutable scalar g_;
+
+    mutable scalar k_;
+
+    bool useGC_;
+
+    bool usePC_;
+
+    virtual void adaptVoidfraction(double&, label) const {}
+
+    virtual scalar effDiameter(double d, label cellI, label index) const {return d;}
+    
+    virtual scalar meanSauterDiameter(double d, label cellI) const {return d;}
+
+    double F(double, double) const;
+
+    double terminalVelocity(double, double, double, double, double, double) const; 
+
+    double a(double) const;
+
+    double h(double) const;
+
+
+
 public:
 
     //- Runtime type information

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

@@ -0,0 +1,116 @@
+/*---------------------------------------------------------------------------*\
+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 "BeetstraDragPoly.H"
+#include "addToRunTimeSelectionTable.H"
+#include "averagingModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(BeetstraDragPoly, 0);
+
+addToRunTimeSelectionTable
+(
+    forceModel,
+    BeetstraDragPoly,
+    dictionary
+);
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from components
+BeetstraDragPoly::BeetstraDragPoly
+(
+    const dictionary& dict,
+    cfdemCloud& sm
+)
+:
+    BeetstraDrag(dict,sm),
+    fines_(propsDict_.lookupOrDefault<bool>("fines",false)),
+    dFine_(1.0)
+{
+    // if fines are present, take mixture dSauter, otherwise normal dSauter
+    if (fines_)
+    {
+        dFine_ = readScalar(propsDict_.lookup("dFine"));
+        volScalarField& alphaP(const_cast<volScalarField&>(sm.mesh().lookupObject<volScalarField> ("alphaP")));
+        alphaP_.set(&alphaP);
+        volScalarField& alphaSt(const_cast<volScalarField&>(sm.mesh().lookupObject<volScalarField> ("alphaSt")));
+        alphaSt_.set(&alphaSt);
+        volScalarField& dSauter(const_cast<volScalarField&>(sm.mesh().lookupObject<volScalarField> ("dSauterMix")));
+        dSauter_.set(&dSauter);
+    }
+    else
+    {
+        volScalarField& dSauter(const_cast<volScalarField&>(sm.mesh().lookupObject<volScalarField> ("dSauter")));
+        dSauter_.set(&dSauter);
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+BeetstraDragPoly::~BeetstraDragPoly()
+{}
+
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+void BeetstraDragPoly::adaptVoidfraction(double& voidfraction, label cellI) const
+{
+    if (fines_) voidfraction -= alphaSt_()[cellI];
+    if (voidfraction < minVoidfraction_) voidfraction = minVoidfraction_;
+}
+
+scalar BeetstraDragPoly::effDiameter(double d, label cellI, label index) const
+{
+    scalar dS = dSauter_()[cellI];
+    scalar effD = d*d / dS + 0.064*d*d*d*d / (dS*dS*dS);
+    
+    if (fines_)
+    {
+        scalar fineCorr = dFine_*dFine_ / dS + 0.064*dFine_*dFine_*dFine_*dFine_ / (dS*dS*dS);
+        fineCorr *= d*d*d / (dFine_*dFine_*dFine_) * alphaSt_()[cellI] / alphaP_()[cellI];
+        effD += fineCorr;
+    }
+    
+    if (particleCloud_.getParticleEffVolFactors()) 
+    {
+        scalar effVolFac = particleCloud_.particleEffVolFactor(index);
+        effD *= effVolFac;
+    }
+
+    return effD;
+}
+
+scalar BeetstraDragPoly::meanSauterDiameter(double d, label cellI) const
+{
+    return dSauter_()[cellI];
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

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

@@ -0,0 +1,93 @@
+/*---------------------------------------------------------------------------*\
+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
+    BeetstraDragPoly.C
+\*---------------------------------------------------------------------------*/
+
+#ifndef BeetstraDragPoly_H
+#define BeetstraDragPoly_H
+
+#include "BeetstraDrag.H"
+#include "interpolationCellPoint.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class BeetstraDragPoly Declaration
+\*---------------------------------------------------------------------------*/
+
+class BeetstraDragPoly
+:
+    public BeetstraDrag
+{
+protected:
+
+    const bool fines_;
+    
+    scalar dFine_;
+
+    autoPtr<volScalarField> alphaP_;
+    
+    autoPtr<volScalarField> alphaSt_;
+
+    autoPtr<volScalarField> dSauter_;
+
+    void adaptVoidfraction(double&, label) const;
+
+    scalar effDiameter(double, label, label) const;
+    
+    scalar meanSauterDiameter(double, label) const;
+
+public:
+
+    //- Runtime type information
+    TypeName("BeetstraDragPoly");
+
+
+    // Constructors
+
+        //- Construct from components
+        BeetstraDragPoly
+        (
+            const dictionary& dict,
+            cfdemCloud& sm
+        );
+
+    // Destructor
+
+        ~BeetstraDragPoly();
+
+
+    // Member Functions
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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

@@ -148,10 +148,10 @@ void DiFeliceDrag::setForce() const
     for(int index = 0; index < particleCloud_.numberOfParticles(); ++index)
     {
             cellI = particleCloud_.cellIDs()[index][0];
-            drag = vector(0,0,0);
+            drag = vector::zero;
-            dragExplicit = vector(0,0,0);
+            dragExplicit = vector::zero;
             dragCoefficient=0;
-            Ufluid =vector(0,0,0);
+            Ufluid = vector::zero;
 
             if (cellI > -1) // particle Found
             {

src/lagrangian/cfdemParticle/subModels/forceModel/Fines/ErgunStatFines.C

@@ -164,10 +164,10 @@ void ErgunStatFines::setForce() const
     for(int index = 0; index < particleCloud_.numberOfParticles(); ++index)
     {
             cellI = particleCloud_.cellIDs()[index][0];
-            drag = vector(0,0,0);
+            drag = vector::zero;
-            dragExplicit = vector(0,0,0);
+            dragExplicit = vector::zero;
             betaP = 0;
-            Ufluid = vector(0,0,0);
+            Ufluid = vector::zero;
             voidfraction = 0;
             dragCoefficient = 0;
 

src/lagrangian/cfdemParticle/subModels/forceModel/Fines/FanningDynFines.C

@@ -120,8 +120,8 @@ void FanningDynFines::setForce() const
     for(int index = 0; index < particleCloud_.numberOfParticles(); ++index)
     {
             cellI = particleCloud_.cellIDs()[index][0];
-            drag = vector(0,0,0);
+            drag = vector::zero;
-            UDyn = vector(0,0,0);
+            UDyn = vector::zero;
             dragCoefficient = 0;
 
             if (cellI > -1) // particle found

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

@@ -159,11 +159,11 @@ void GidaspowDrag::setForce() const
         //if(mask[index][0])
         //{
             cellI = particleCloud_.cellIDs()[index][0];
-            drag = vector(0,0,0);
+            drag = vector::zero;
-            dragExplicit = vector(0,0,0);
+            dragExplicit = vector::zero;
             betaP = 0;
             Vs = 0;
-            Ufluid = vector(0,0,0);
+            Ufluid = vector::zero;
             voidfraction = 0;
             dragCoefficient = 0;
 

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

@@ -157,12 +157,12 @@ void KochHillDrag::setForce() const
     for (int index = 0; index<particleCloud_.numberOfParticles(); ++index)
     {
             cellI = particleCloud_.cellIDs()[index][0];
-            drag = vector(0,0,0);
+            drag = vector::zero;
-            dragExplicit = vector(0,0,0);
+            dragExplicit = vector::zero;
             dragCoefficient = 0;
             betaP = 0;
             Vs = 0;
-            Ufluid = vector(0,0,0);
+            Ufluid = vector::zero;
             voidfraction = 0;
 
             if (cellI > -1) // particle Found

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

@@ -200,12 +200,12 @@ void KochHillRWDrag::setForce() const
         //if (mask[index][0])
         //{
             cellI = particleCloud_.cellIDs()[index][0];
-            drag = vector(0,0,0);
+            drag = vector::zero;
-            dragExplicit = vector(0,0,0);
+            dragExplicit = vector::zero;
             dragCoefficient = 0;
             betaP = 0;
             Vs = 0;
-            Ufluid = vector(0,0,0);
+            Ufluid = vector::zero;
 
             // Pout << "RW-TEST: cellI = " << cellI << endl; // TEST-Output
             if (cellI > -1) // particle Found

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

@@ -53,10 +53,10 @@ dSauter::dSauter
 :
     forceModel(dict,sm),
     propsDict_(dict.subDict(typeName + "Props")),
+    multiTypes_(false),
     d2_(NULL),
     d3_(NULL),
-    scaleDia_(1.0),
+    typeCG_(propsDict_.lookupOrDefault<scalarList>("coarseGrainingFactors",scalarList(1,1.0))),
-    scaleDiaDist_(1.0),
     d2Field_
     (   IOobject
         (
@@ -95,17 +95,13 @@ dSauter::dSauter
         "zeroGradient"
     )
 {
+    if (typeCG_.size()>1) multiTypes_ = true;
     allocateMyArrays();
     dSauter_.write();
 
 
     // init force sub model
     setForceSubModels(propsDict_);
-
-    if (propsDict_.found("scaleCG"))
-        scaleDia_ = scalar(readScalar(propsDict_.lookup("scaleCG")));
-    if (propsDict_.found("scaleDist"))
-        scaleDiaDist_ = scalar(readScalar(propsDict_.lookup("scaleDist")));
 }
 
 
@@ -131,30 +127,37 @@ void dSauter::allocateMyArrays() const
 
 void dSauter::setForce() const
 {
-    if (scaleDia_ > 1)
+    if (typeCG_.size()>1 || typeCG_[0] > 1.0)
     {
-        Info << "dSauter using scaleCG = " << scaleDia_ << endl;
+        Info << "dSauter using CG factor(s) = " << typeCG_ << endl;
-    }
-    else if (particleCloud_.cg() > 1)
-    {
-        scaleDia_ = particleCloud_.cg();
-        Info << "dSauter using scaleCG from liggghts cg = " << scaleDia_ << endl;
     }
 
     allocateMyArrays();
 
-    label cellI=0;
+    label cellI = 0;
-    scalar ds(0);
+    label partType = 1;
-    scalar scale = scaleDiaDist_/scaleDia_;
+    scalar cg = typeCG_[0];
+    scalar ds = 0.0;
+    scalar effVolFac = 1.0;
+
 
     for(int index = 0; index < particleCloud_.numberOfParticles(); ++index)
     {
         cellI = particleCloud_.cellIDs()[index][0];
         if (cellI >= 0)
         {
+            if (particleCloud_.getParticleEffVolFactors()) 
+            {
+                effVolFac = particleCloud_.particleEffVolFactor(index);
+            }
+            if (multiTypes_) 
+            {
+                partType = particleCloud_.particleType(index);
+                cg = typeCG_[partType - 1];
+            }
             ds = particleCloud_.d(index);
-            d2_[index][0] = ds*ds;
+            d2_[index][0] = ds*ds*effVolFac*cg;
-            d3_[index][0] = ds*ds*ds;
+            d3_[index][0] = ds*ds*ds*effVolFac;
         }
     }
 
@@ -181,7 +184,7 @@ void dSauter::setForce() const
     {
         if (d2Field_[cellI] > ROOTVSMALL)
         {
-            dSauter_[cellI] = d3Field_[cellI] / d2Field_[cellI] * scale;
+            dSauter_[cellI] = d3Field_[cellI] / d2Field_[cellI];
         }
         else
         {

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

@@ -43,13 +43,13 @@ private:
 
     dictionary propsDict_;
 
+    bool multiTypes_;
+
     mutable double **d2_;
     
     mutable double **d3_;
 
-    mutable scalar scaleDia_;      // scaling due to coarse graining
+    scalarList typeCG_;
-
-    mutable scalar scaleDiaDist_;  // scaling due to modified particle size distribution
     
     mutable volScalarField d2Field_;
     

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

@@ -69,7 +69,7 @@ forceModel::forceModel
             IOobject::AUTO_WRITE
         ),
         sm.mesh(),
-        dimensionedVector("zero", dimensionSet(1,1,-2,0,0), vector(0,0,0)) // N
+        dimensionedVector("zero", dimensionSet(1,1,-2,0,0), vector::zero) // N
     ),
     expParticleForces_
     (   IOobject
@@ -81,7 +81,7 @@ forceModel::forceModel
             IOobject::AUTO_WRITE
         ),
         sm.mesh(),
-        dimensionedVector("zero", dimensionSet(1,1,-2,0,0), vector(0,0,0)) // N
+        dimensionedVector("zero", dimensionSet(1,1,-2,0,0), vector::zero) // N
     ),
     coupleForce_(true),
     modelType_(sm.modelType()),

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

@@ -160,7 +160,7 @@ void gradPForce::setForce() const
     {
         //if(mask[index][0])
         //{
-            force = vector(0,0,0);
+            force = vector::zero;
             cellI = particleCloud_.cellIDs()[index][0];
 
             if (cellI > -1) // particle Found

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

@@ -142,7 +142,7 @@ Info << "interface::setForce" << endl;
                 }
 
                 // Initialize an interfaceForce vector
-                vector interfaceForce = Foam::vector(0,0,0);
+                vector interfaceForce = Foam::vector::zero;
 
                 // Calculate the interfaceForce (range of alphap needed for stability)
 

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

@@ -0,0 +1,317 @@
+/*---------------------------------------------------------------------------*\
+License
+    This is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+    This code is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+    You should have received a copy of the GNU General Public License
+    along with this code.  If not, see <http://www.gnu.org/licenses/>.
+
+Description
+    calculates the correlation between particle momentum and diameter
+
+SourceFiles
+    pdCorrelation.C
+
+Contributing Author
+    2018    Paul Kieckhefen, TUHH 
+\*---------------------------------------------------------------------------*/
+
+#include "error.H"
+
+#include "pdCorrelation.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(pdCorrelation, 0);
+
+addToRunTimeSelectionTable
+(
+    forceModel,
+    pdCorrelation,
+    dictionary
+);
+
+const dimensionSet dimMomentum(dimForce*dimTime);
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from components
+pdCorrelation::pdCorrelation
+(
+    const dictionary& dict,
+    cfdemCloud& sm
+)
+:
+    forceModel(dict,sm),
+    propsDict_(dict.subDict(typeName + "Props")),
+    d_(nullptr),
+    p_(nullptr),
+    d2_(nullptr),
+    pd_(nullptr),
+    cg3_(nullptr),
+    dField_
+    (   IOobject
+        (
+            "dField",
+            sm.mesh().time().timeName(),
+            sm.mesh(),
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        sm.mesh(),
+        dimensionedScalar("zero", dimLength, 0),
+        "zeroGradient"
+    ),
+    pField_
+    (   IOobject
+        (
+            "pField",
+            sm.mesh().time().timeName(),
+            sm.mesh(),
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        sm.mesh(),
+        dimensionedVector("zero", dimMomentum, vector::zero)
+    ),
+    d2Field_
+    (   IOobject
+        (
+            "d2Field",
+            sm.mesh().time().timeName(),
+            sm.mesh(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        sm.mesh(),
+        dimensionedScalar("zero", dimLength*dimLength, 0)
+    ),
+    pdField_
+    (   IOobject
+        (
+            "pdCorrelation",
+            sm.mesh().time().timeName(),
+            sm.mesh(),
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        sm.mesh(),
+        dimensionedVector("zero", dimMomentum*dimLength, vector::zero),
+        "zeroGradient"
+    ),
+    cg3Field_
+    (   IOobject
+        (
+            "cg3Field",
+            sm.mesh().time().timeName(),
+            sm.mesh(),
+            IOobject::NO_READ,
+            IOobject::NO_WRITE
+        ),
+        sm.mesh(),
+        dimensionedScalar("zero", dimless, 0)
+    ),
+    typeCG_
+    (
+        propsDict_.lookupOrDefault
+        (
+            "coarseGrainingFactors",
+            scalarList(1, sm.cg())
+        )
+    ),
+    particleDensities_
+    (
+        propsDict_.lookupOrDefault
+        (
+            "particleDensities",
+            scalarList(1, -1.)
+        )
+    ),
+    constantCG_(typeCG_.size() < 2),
+    CG_(!constantCG_ || typeCG_[0] > 1. + SMALL),
+    runOnWriteOnly_(propsDict_.lookupOrDefault("runOnWriteOnly", false))
+{
+    if ((particleDensities_[0] < 0) && !particleCloud_.getParticleDensities())
+    {
+        FatalError<< "Please set the particle density either in LIGGGHTS"
+                  << "or the pdCorrelationPropsDict."
+                  << "In the first case, set getParticleDensities in the couplingProperties."
+                  << abort(FatalError);
+    }
+
+    allocateMyArrays();
+
+    dField_.write();
+    pdField_.write();
+
+
+    // init force sub model
+    setForceSubModels(propsDict_);
+}
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+pdCorrelation::~pdCorrelation()
+{
+    particleCloud_.dataExchangeM().destroy(cg3_, 1);
+    particleCloud_.dataExchangeM().destroy(d_,  1);
+    particleCloud_.dataExchangeM().destroy(p_,  3);
+    particleCloud_.dataExchangeM().destroy(d2_, 1);
+    particleCloud_.dataExchangeM().destroy(pd_, 3);
+}
+
+// * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
+void pdCorrelation::allocateMyArrays() const
+{
+    // get memory for 2d arrays
+    double initVal = 0.0;
+    particleCloud_.dataExchangeM().allocateArray(d_,  initVal, 1);
+    particleCloud_.dataExchangeM().allocateArray(p_,  initVal, 3);
+    particleCloud_.dataExchangeM().allocateArray(d2_, initVal, 1);
+    particleCloud_.dataExchangeM().allocateArray(pd_, initVal, 3);
+    particleCloud_.dataExchangeM().allocateArray(cg3_, initVal, 1);
+}
+// * * * * * * * * * * * * * * * public Member Functions  * * * * * * * * * * * * * //
+
+void pdCorrelation::setForce() const
+{
+    const fvMesh& mesh = particleCloud_.mesh();
+
+    if (runOnWriteOnly_ && !mesh.write()) return; // skip if it's not write time
+
+    allocateMyArrays();
+
+    const Switch densityFromList
+    (
+        particleCloud_.getParticleTypes()
+        && !particleCloud_.getParticleDensities()
+    );
+    const scalar piOverSix = constant::mathematical::pi / 6.;
+
+    dField_.primitiveFieldRef()   = 0.0;
+    pField_.primitiveFieldRef()   = vector::zero;
+    d2Field_.primitiveFieldRef()  = 0.0;
+    pdField_.primitiveFieldRef()  = vector::zero;
+    cg3Field_.primitiveFieldRef() = 0.;
+
+    label  celli(0);
+    scalar dp(0.);
+    scalar cg(typeCG_[0]);
+    scalar rhop(particleDensities_[0]);
+    label  typei(0);
+    scalar particleV(0.);
+    for(int pi = 0; pi < particleCloud_.numberOfParticles(); ++pi)
+    {
+        celli = particleCloud_.cellIDs()[pi][0];
+        if (celli < 0) continue;
+
+        dp = particleCloud_.d(pi);
+        if (particleCloud_.getParticleTypes())
+        {
+            typei = particleCloud_.particleTypes()[pi][0] - 1;
+            rhop  = densityFromList ? particleDensities_[typei]
+                                    : particleCloud_.particleDensity(pi);
+            cg    = constantCG_     ? typeCG_[0] : typeCG_[typei];
+        }
+
+        particleV  = piOverSix * rhop * dp * dp * dp;
+        d_ [pi][0] = cg * cg * dp;
+        d2_[pi][0] = cg * dp * dp;
+        p_ [pi][0] = particleV * particleCloud_.velocities()[pi][0];
+        p_ [pi][1] = particleV * particleCloud_.velocities()[pi][1];
+        p_ [pi][2] = particleV * particleCloud_.velocities()[pi][2];
+        pd_[pi][0] = p_[pi][0] * dp / cg;
+        pd_[pi][1] = p_[pi][1] * dp / cg;
+        pd_[pi][2] = p_[pi][2] * dp / cg;
+
+        cg3_[pi][0] = CG_ ? cg * cg * cg : 1.;
+    }
+
+
+    particleCloud_.averagingM().setScalarSum
+    (
+        dField_,
+        d_,
+        particleCloud_.particleWeights(),
+        nullptr
+    );
+    particleCloud_.averagingM().setVectorSum
+    (
+        pField_,
+        p_,
+        particleCloud_.particleWeights(),
+        nullptr
+    );
+    particleCloud_.averagingM().setScalarSum
+    (
+        d2Field_,
+        d2_,
+        particleCloud_.particleWeights(),
+        nullptr
+    );
+    particleCloud_.averagingM().setVectorSum
+    (
+        pdField_,
+        pd_,
+        particleCloud_.particleWeights(),
+        nullptr
+    );
+
+    particleCloud_.averagingM().setScalarSum
+    (
+        cg3Field_,
+        cg3_,
+        particleCloud_.particleWeights(),
+        nullptr
+    );
+
+    scalar oneOverCG3(1.);
+    const scalar oneMinusSmall(1.-SMALL);
+    forAll(cg3Field_, celli)
+    {
+        if (cg3Field_[celli] < oneMinusSmall) continue;
+        oneOverCG3 = 1./cg3Field_[celli];
+
+        dField_ [celli] *= oneOverCG3;
+        pField_ [celli] *= oneOverCG3;
+        d2Field_[celli] *= oneOverCG3;
+        pdField_[celli] *= oneOverCG3;
+    }
+
+    scalar denom(0.);
+    forAll(dField_, celli)
+    {
+        if (dField_[celli] < SMALL)
+        {
+            pdField_[celli] = vector::zero;
+            continue;
+        }
+        denom = d2Field_[celli] - dField_[celli] * dField_[celli];
+        if (denom < SMALL)
+        {
+            pdField_[celli] = vector::zero;
+            continue;
+        }
+
+        pdField_[celli] -= dField_[celli] * pField_[celli];
+        pdField_[celli] /= denom;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

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

@@ -0,0 +1,102 @@
+/*---------------------------------------------------------------------------*\
+License
+    This is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+    This code is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+    You should have received a copy of the GNU General Public License
+    along with this code.  If not, see <http://www.gnu.org/licenses/>.
+
+
+Description
+    calculates the correlation between particle diameter and momentum
+
+SourceFiles
+    pdCorrelation.C
+
+Contributing Author
+    2018    Paul Kieckhefen, TUHH
+\*---------------------------------------------------------------------------*/
+
+#ifndef pdCorrelation_H
+#define pdCorrelation_H
+
+#include "forceModel.H"
+#include "averagingModel.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class pdCorrelation Declaration
+\*---------------------------------------------------------------------------*/
+
+class pdCorrelation
+:
+    public forceModel
+{
+private:
+
+    dictionary propsDict_;
+
+    mutable double **d_;
+    mutable double **p_;
+    mutable double **d2_;
+    mutable double **pd_;
+    mutable double **cg3_;
+
+    mutable volScalarField dField_;
+    mutable volVectorField pField_;
+    mutable volScalarField d2Field_;
+    mutable volVectorField pdField_;
+
+    mutable volScalarField cg3Field_;
+
+    const scalarList typeCG_;
+    const scalarList particleDensities_;
+    const Switch     constantCG_;
+    const Switch     CG_;
+    const Switch     runOnWriteOnly_;
+
+    void allocateMyArrays() const;
+
+public:
+
+    //- Runtime type information
+    TypeName("pdCorrelation");
+
+    // Constructors
+
+        //- Construct from components
+        pdCorrelation
+        (
+            const dictionary& dict,
+            cfdemCloud& sm
+        );
+
+    // Destructor
+
+        ~pdCorrelation();
+
+
+    // Member Functions
+        void setForce() const;
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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

@@ -142,7 +142,7 @@ void viscForce::setForce() const
     {
         //if(mask[index][0])
         //{
-            force=vector(0,0,0);
+            force = vector::zero;
             cellI = particleCloud_.cellIDs()[index][0];
 
             if (cellI > -1) // particle Found

src/lagrangian/cfdemParticle/subModels/forceModelMS/DiFeliceDragMS/DiFeliceDragMS.C

@@ -156,7 +156,7 @@ void DiFeliceDragMS::setForce() const
         //{
 
             cellI = cloudRefMS().cellIDCM(index);
-            drag = vector(0,0,0);
+            drag = vector::zero;
 
             if (cellI > -1) // particle Found
             {

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

@@ -73,7 +73,7 @@ explicitCouple::explicitCouple
             IOobject::NO_WRITE
         ),
         sm.mesh(),
-        dimensionedVector("zero", dimensionSet(1,-2,-2,0,0), vector(0,0,0)), // N/m3
+        dimensionedVector("zero", dimensionSet(1,-2,-2,0,0), vector::zero), // N/m3
         "zeroGradient"
     ),
     fNext_
@@ -86,7 +86,7 @@ explicitCouple::explicitCouple
             IOobject::NO_WRITE
         ),
         sm.mesh(),
-        dimensionedVector("zero", dimensionSet(1,-2,-2,0,0), vector(0,0,0)), // N/m3
+        dimensionedVector("zero", dimensionSet(1,-2,-2,0,0), vector::zero), // N/m3
         "zeroGradient"
     ),
     fLimit_(1e10,1e10,1e10)

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

@@ -69,7 +69,7 @@ tmp<volVectorField> expParticleForces::exportForceField()
             particleCloud_.mesh(),
             dimensionedVector
             (
-                "zero",dimensionSet(1, -2, -2, 0, 0),vector(0,0,0)
+                "zero",dimensionSet(1, -2, -2, 0, 0),vector::zero
             )
         )
     );

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

@@ -45,7 +45,6 @@ protected:
     // Protected data
     //    dictionary propsDict_;
 
-    
 public:
 
     //- Runtime type information

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

@@ -74,7 +74,7 @@ tmp<volVectorField> gravity::exportForceField()
             particleCloud_.mesh(),
             dimensionedVector
             (
-                "zero",dimensionSet(1, -2, -2, 0, 0),vector(0,0,0)
+                "zero",dimensionSet(1, -2, -2, 0, 0),vector::zero
             )
         )
     );

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

@@ -55,7 +55,6 @@ protected:
 
         const dimensionedVector g_;
 
-    
 public:
 
     //- Runtime type information

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

@@ -47,7 +47,6 @@ protected:
 
         cfdemCloud& particleCloud_;
 
-        
 public:
 
     //- Runtime type information
@@ -95,7 +94,6 @@ public:
 
     // Member Functions
 
-
         virtual tmp<volVectorField> exportForceField() = 0;
 
 };

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

@@ -78,7 +78,7 @@ tmp<volVectorField> weightSecondaryPhase::exportForceField()
             particleCloud_.mesh(),
             dimensionedVector
             (
-                "zero",dimensionSet(1, -2, -2, 0, 0),vector(0,0,0)
+                "zero",dimensionSet(1, -2, -2, 0, 0),vector::zero
             )
         )
     );

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

@@ -0,0 +1,225 @@
+/*---------------------------------------------------------------------------*\
+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 "ZehnerSchluenderThermCond.H"
+#include "addToRunTimeSelectionTable.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(ZehnerSchluenderThermCond, 0);
+
+addToRunTimeSelectionTable
+(
+    thermCondModel,
+    ZehnerSchluenderThermCond,
+    dictionary
+);
+
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from components
+ZehnerSchluenderThermCond::ZehnerSchluenderThermCond
+(
+    const dictionary& dict,
+    cfdemCloud& sm
+)
+:
+    thermCondModel(dict,sm),
+    propsDict_(dict.subDict(typeName + "Props")),
+    partKsField_
+    (
+        IOobject
+        (
+            "partKs",
+            sm.mesh().time().timeName(),
+            sm.mesh(),
+            IOobject::READ_IF_PRESENT,
+            IOobject::NO_WRITE
+        ),
+        sm.mesh(),
+        dimensionedScalar("one", dimensionSet(1, 1, -3, -1,0,0,0), 1.0),
+	"zeroGradient"
+    ),
+    voidfractionFieldName_(propsDict_.lookupOrDefault<word>("voidfractionFieldName","voidfraction")),
+    voidfraction_(sm.mesh().lookupObject<volScalarField> (voidfractionFieldName_)),
+    typeKs_(propsDict_.lookupOrDefault<scalarList>("thermalConductivities",scalarList(1,-1.0))),
+    partKs_(NULL),
+    wallQFactorName_(propsDict_.lookupOrDefault<word>("wallQFactorName","wallQFactor")),
+    wallQFactor_
+    (   IOobject
+        (
+            wallQFactorName_,
+            sm.mesh().time().timeName(),
+            sm.mesh(),
+            IOobject::READ_IF_PRESENT,
+            IOobject::AUTO_WRITE
+        ),
+        sm.mesh(),
+        dimensionedScalar("zero", dimensionSet(0,0,0,0,0,0,0), 1.0)
+    ),
+    hasWallQFactor_(false)
+{
+    if (typeKs_[0] < 0.0)
+    {
+        FatalError << "ZehnerSchluenderThermCond: provide list of thermal conductivities." << abort(FatalError);
+    }
+
+    if (typeKs_.size() > 1) allocateMyArrays();
+    else
+    {
+        partKsField_ *= typeKs_[0];
+    }
+
+    if (wallQFactor_.headerOk())
+    {
+        Info << "Found field for scaling wall heat flux.\n" << endl;
+        hasWallQFactor_ = true;
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+ZehnerSchluenderThermCond::~ZehnerSchluenderThermCond()
+{
+    if (typeKs_.size() > 1) particleCloud_.dataExchangeM().destroy(partKs_,1);
+}
+
+
+void ZehnerSchluenderThermCond::allocateMyArrays() const
+{
+    double initVal=0.0;
+    particleCloud_.dataExchangeM().allocateArray(partKs_,initVal,1);
+}
+
+// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
+
+tmp<volScalarField> ZehnerSchluenderThermCond::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),
+            "zeroGradient"
+        )
+    );
+
+    volScalarField& svf = tvf.ref();
+
+    calcPartKsField();
+    scalar A = 0.0;
+    scalar B = 0.0;
+    scalar C = 0.0;
+    scalar k = 0.0;
+    scalar InvOnemBoA = 0.0;
+    scalar voidfraction = 0.0;
+    scalar w = 7.26e-3;
+
+    forAll(svf, cellI)
+    {
+        voidfraction = voidfraction_[cellI];
+        if(voidfraction > 1.0 - SMALL || partKsField_[cellI] < SMALL) svf[cellI] = 0.0;
+        else
+        {
+            A = partKsField_[cellI]/kf0_.value();
+            B = 1.25 * Foam::pow((1 - voidfraction) / voidfraction, 1.11);
+            InvOnemBoA = 1.0/(1.0 - B/A);
+            C = (A - 1) * InvOnemBoA * InvOnemBoA * B/A * log(A/B) - (B - 1) * InvOnemBoA - 0.5 * (B + 1);
+            C *= 2.0 * InvOnemBoA;
+            k = Foam::sqrt(1 - voidfraction) * (w * A + (1 - w) * C) * kf0_.value();
+            svf[cellI] = k / (1 - voidfraction);
+        }
+    }
+
+    svf.correctBoundaryConditions();
+
+    // if a wallQFactor field is present, use it to scale heat transport through a patch
+    if (hasWallQFactor_)
+    {
+        wallQFactor_.correctBoundaryConditions();
+        forAll(wallQFactor_.boundaryField(), patchi)
+            svf.boundaryFieldRef()[patchi] *= wallQFactor_.boundaryField()[patchi];
+    }
+
+    return tvf;
+}
+
+tmp<volScalarField> ZehnerSchluenderThermCond::thermDiff() const
+{
+    FatalError << "ZehnerSchluenderThermCond does not provide thermal diffusivity." << abort(FatalError);
+    return tmp<volScalarField>(NULL);
+}
+
+void ZehnerSchluenderThermCond::calcPartKsField() const
+{
+    if (typeKs_.size() <= 1) return;
+
+    if (!particleCloud_.getParticleTypes())
+    {
+        FatalError << "ZehnerSchluenderThermCond needs data for more than one type, but types are not communicated." << abort(FatalError);
+    }
+
+    allocateMyArrays();
+    label cellI=0;
+    label partType = 0;
+    for(int index = 0;index < particleCloud_.numberOfParticles(); ++index)
+    {
+            cellI = particleCloud_.cellIDs()[index][0];
+            if(cellI >= 0)
+            {
+                partType = particleCloud_.particleType(index);
+                // LIGGGGHTS counts types 1, 2, ..., C++ array starts at 0
+                partKs_[index][0] = typeKs_[partType - 1];
+            }
+    }
+
+    partKsField_.primitiveFieldRef() = 0.0;
+    particleCloud_.averagingM().resetWeightFields();
+    particleCloud_.averagingM().setScalarAverage
+    (
+        partKsField_,
+        partKs_,
+        particleCloud_.particleWeights(),
+        particleCloud_.averagingM().UsWeightField(),
+        NULL
+    );
+}
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //

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

@@ -0,0 +1,116 @@
+/*---------------------------------------------------------------------------*\
+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 PARTICLE phase in presence of a fluid according to
+    Zehner and Schluender (1970)
+
+Class
+    ZehnerSchluenderThermCond
+
+SourceFiles
+    ZehnerSchluenderThermCond.C
+
+\*---------------------------------------------------------------------------*/
+
+
+#ifndef ZehnerSchluenderThermCond_H
+#define ZehnerSchluenderThermCond_H
+
+#include "thermCondModel.H"
+#include "scalarList.H"
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class ZehnerSchluenderThermCond Declaration
+\*---------------------------------------------------------------------------*/
+
+class ZehnerSchluenderThermCond
+:
+    public thermCondModel
+{
+  
+private:
+  
+    dictionary propsDict_;
+
+    mutable volScalarField partKsField_;
+  
+    word voidfractionFieldName_;
+  
+    const volScalarField& voidfraction_;
+    
+    scalarList typeKs_;
+
+    mutable double **partKs_;
+    
+    word wallQFactorName_;
+    
+    // ratio of half-cell-size and near-wall film
+    mutable volScalarField wallQFactor_;
+    
+    bool hasWallQFactor_;
+
+    void allocateMyArrays() const;
+
+    void calcPartKsField() const;
+
+public:
+
+    //- Runtime type information
+    TypeName("ZehnerSchluenderThermCond");
+
+
+    // Constructors
+
+        //- Construct from components
+        ZehnerSchluenderThermCond
+        (
+            const dictionary& dict,
+            cfdemCloud& sm
+        );
+
+    // Destructor
+
+        ~ZehnerSchluenderThermCond();
+
+
+    // Member Functions
+	
+	tmp<volScalarField> thermCond() const;
+	
+	tmp<volScalarField> thermDiff() const;
+       
+       
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //

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

@@ -95,7 +95,7 @@ void GaussVoidFraction::setvoidFraction(double** const& mask,double**& voidfract
 
     scalar radius(-1);
     scalar volume(0);
-    scalar scaleVol= weight();
+    scalar scaleVol = weight();
     scalar scaleRadius = pow(porosity(),1./3.);
 
     for(int index=0; index< particleCloud_.numberOfParticles(); index++)
@@ -115,6 +115,7 @@ void GaussVoidFraction::setvoidFraction(double** const& mask,double**& voidfract
             label particleCenterCellID=particleCloud_.cellIDs()[index][0];
 
             radius = particleCloud_.radius(index);
+            if (multiWeights_) scaleVol = weight(index);
             volume = constant::mathematical::fourPiByThree*radius*radius*radius*scaleVol;
             radius *= scaleRadius;
 

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

@@ -94,7 +94,7 @@ void bigParticleVoidFraction::setvoidFraction(double** const& mask,double**& voi
 
     scalar radius(-1);
     scalar volume(0);
-    scalar scaleVol= weight();
+    scalar scaleVol = weight();
     scalar scaleRadius = pow(porosity(),1./3.);
 
     for(int index=0; index< particleCloud_.numberOfParticles(); index++)
@@ -113,6 +113,7 @@ void bigParticleVoidFraction::setvoidFraction(double** const& mask,double**& voi
             //collecting data
             label particleCenterCellID=particleCloud_.cellIDs()[index][0];
             radius =  particleCloud_.radius(index);
+            if (multiWeights_) scaleVol = weight(index);
             volume = constant::mathematical::fourPiByThree*radius*radius*radius*scaleVol;
             radius *= scaleRadius;
             vector positionCenter=particleCloud_.position(index);

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

@@ -85,7 +85,7 @@ void centreVoidFraction::setvoidFraction(double** const& mask,double**& voidfrac
     scalar radius(-1);
     scalar volume(0);
     scalar cellVol(0);
-    scalar scaleVol= weight();
+    scalar scaleVol = weight();
 
     for(int index=0; index< particleCloud_.numberOfParticles(); index++)
     {
@@ -99,6 +99,7 @@ void centreVoidFraction::setvoidFraction(double** const& mask,double**& voidfrac
 
             if (cellI >= 0)  // particel centre is in domain
             {
+                if (multiWeights_) scaleVol = weight(index);
                 cellVol = voidfractionNext_.mesh().V()[cellI];
                 radius = particleCloud_.radius(index);
                 volume = constant::mathematical::fourPiByThree*radius*radius*radius*scaleVol;

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

@@ -200,6 +200,7 @@ void dividedVoidFraction::setvoidFraction(double** const& mask,double**& voidfra
             position = particleCloud_.position(index);
             cellID = particleCloud_.cellIDs()[index][0];
             radius = particleCloud_.radius(index);
+            if (multiWeights_) scaleVol = weight(index);
             volume = Vp(index,radius,scaleVol);
             radius *= scaleRadius;
             cellVol = 0;

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

@@ -94,7 +94,7 @@ void trilinearVoidFraction::setvoidFraction(double** const& mask,double**& voidf
 
     scalar radius(-1.);
     scalar volume(0.);
-    const scalar scaleVol = weight();
+    scalar scaleVol = weight();
 
     vector partPos(0.,0.,0.);
     vector pt(0.,0.,0.);
@@ -140,6 +140,7 @@ void trilinearVoidFraction::setvoidFraction(double** const& mask,double**& voidf
         if (cellI >= 0)  // particel centre is in domain
         {
             radius = particleCloud_.radius(index);
+            if (multiWeights_) scaleVol = weight(index);
             volume = constant::mathematical::fourPiByThree * radius * radius * radius * scaleVol;
 
             // store volume for each particle

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

@@ -57,6 +57,7 @@ voidFractionModel::voidFractionModel
 :
     dict_(dict),
     particleCloud_(sm),
+    multiWeights_(false),
     voidfractionPrev_
     (   IOobject
         (
@@ -89,6 +90,7 @@ voidFractionModel::voidFractionModel
     porosity_(1.)
 {
     particleCloud_.dataExchangeM().allocateArray(cellsPerParticle_,1,1);
+    if (particleCloud_.getParticleEffVolFactors()) multiWeights_ = true;
 }
 
 

src/lagrangian/cfdemParticle/subModels/voidFractionModel/voidFractionModel/voidFractionModel.H

@@ -63,6 +63,8 @@ protected:
 
         cfdemCloud& particleCloud_;
 
+        bool multiWeights_;
+
         mutable volScalarField voidfractionPrev_;
 
         mutable volScalarField voidfractionNext_;
@@ -140,6 +142,11 @@ public:
 
         inline scalar weight()const { return weight_; }
 
+        inline scalar weight(label index) const
+        {
+            return particleCloud_.particleEffVolFactor(index);
+        }
+
         inline scalar porosity()const { return porosity_; }
 
         inline void checkWeightNporosity(dictionary& propsDict) const

src/lagrangian/cfdemParticleComp/Make/files

@@ -31,12 +31,13 @@ $(cfdTools)/newGlobal.C
 $(energyModels)/energyModel/energyModel.C
 $(energyModels)/energyModel/newEnergyModel.C
 $(energyModels)/heatTransferGunn/heatTransferGunn.C
-$(energyModels)/heatTransferGunnImplicit/heatTransferGunnImplicit.C
+$(energyModels)/heatTransferGunnPartField/heatTransferGunnPartField.C
 $(energyModels)/reactionHeat/reactionHeat.C
 
 $(thermCondModels)/thermCondModel/thermCondModel.C
 $(thermCondModels)/thermCondModel/newThermCondModel.C
 $(thermCondModels)/SyamlalThermCond/SyamlalThermCond.C
+$(thermCondModels)/ZehnerSchluenderThermCond/ZehnerSchluenderThermCond.C
 $(thermCondModels)/noTherm/noThermCond.C
 
 $(forceModels)/forceModel/forceModel.C
@@ -64,6 +65,7 @@ $(forceModels)/isotropicFluctuations/isotropicFluctuations.C
 $(forceModels)/evaluateFluctuations/evaluateFluctuations.C
 $(forceModels)/directedDiffusiveRelaxation/directedDiffusiveRelaxation.C
 $(forceModels)/BeetstraDrag/BeetstraDrag.C
+$(forceModels)/BeetstraDragPoly/BeetstraDragPoly.C
 $(forceModels)/dSauter/dSauter.C
 $(forceModels)/Fines/Fines.C
 $(forceModels)/Fines/FinesFields.C

tutorials/.gitignore

@@ -6,3 +6,16 @@
 log_*
 log.*
 *~
+**/DEM/post/*
+**/CFD/processor*/
+**/CFD/0*/
+**/CFD/1*/
+**/CFD/2*/
+**/CFD/3*/
+**/CFD/4*/
+**/CFD/5*/
+**/CFD/6*/
+**/CFD/7*/
+**/CFD/8*/
+**/CFD/9*/
+!**/CFD/0/

tutorials/cfdemSolverIB/twoSpheresGlowinskiMPI/parCFDDEMrun.sh

@@ -22,6 +22,7 @@ solverName="cfdemSolverIB"
 nrProcs="4"
 machineFileName="none"   # yourMachinefileName | none
 debugMode="off"          # on | off| strict
+reconstructCase="false"  # true | false
 testHarnessPath="$CFDEM_TEST_HARNESS_PATH"
 runOctave="true"
 postproc="false"
@@ -30,6 +31,14 @@ postproc="false"
 #- call function to run a parallel CFD-DEM case
 parCFDDEMrun $logpath $logfileName $casePath $headerText $solverName $nrProcs $machineFileName $debugMode
 
+#- case needs special reconstruction
+if [ $reconstructCase == "true" ]
+  then
+    cd $casePath/CFD
+    reconstructParMesh -mergeTol 1e-06
+    reconstructPar -noLagrangian
+fi
+
 if [ $runOctave == "true" ]
   then
 
@@ -60,9 +69,6 @@ if [ $postproc == "true" ]
     read
 fi
 
-#- copy log file to test harness
-cp ../../$logfileName $testHarnessPath
-
 #- clean up case
 echo "deleting data at: $casePath"
 source $WM_PROJECT_DIR/bin/tools/CleanFunctions

tutorials/cfdemSolverPiso/ErgunTestCG/CFD/constant/couplingProperties

@@ -33,7 +33,7 @@ couplingInterval 100;
 
 voidFractionModel divided;//centre;//
 
-locateModel engine;//turboEngine;//
+locateModel engine;//turboEngineM2M;//
 
 meshMotionModel noMeshMotion;
 
@@ -49,7 +49,7 @@ averagingModel dense;//dilute;//
 
 clockModel off;//standardClock;//
 
-smoothingModel off;// constDiffSmoothing; //
+smoothingModel off;// localPSizeDiffSmoothing;// constDiffSmoothing; //
 
 forceModels
 (
@@ -61,6 +61,7 @@ forceModels
     GidaspowDrag
     //Archimedes
     //volWeightedAverage
+    //totalMomentumExchange
     particleCellVolume
 );
 
@@ -74,6 +75,14 @@ turbulenceModelType "turbulenceProperties";
 //===========================================================================//
 // sub-model properties
 
+localPSizeDiffSmoothingProps
+{
+    lowerLimit 0.1;
+    upperLimit 1e10;
+    dSmoothingLength 1.5e-3;
+    Csmoothing    1.0; 
+}
+
 constDiffSmoothingProps
 {
     lowerLimit 0.1;
@@ -123,7 +132,13 @@ volWeightedAverageProps
     lowerThreshold 0;
     verbose true;
 }
-
+totalMomentumExchangeProps
+{
+    implicitMomExFieldName "Ksl";
+    explicitMomExFieldName "none";
+    fluidVelFieldName "U";
+    granVelFieldName "Us";
+}
 GidaspowDragProps
 {
     verbose true;
@@ -147,12 +162,15 @@ KochHillDragProps
 BeetstraDragProps
 {
     velFieldName "U";
-    gravityFieldName "g";
-    rhoParticle     2000.;
     voidfractionFieldName "voidfraction";
+    granVelFieldName "Us";
     interpolation true;
-    useFilteredDragModel ;
+//    useFilteredDragModel;
-    useParcelSizeDependentFilteredDrag ;
+//    useParcelSizeDependentFilteredDrag;
+    g       9.81;
+    rhoP    7000.;
+    rho     10.;
+    nuf     1.5e-4;
     k       0.05;
     aLimit  0.0;
 //    verbose true;

tutorials/cfdemSolverPiso/ErgunTestMPI_restart/parCFDDEMrun.sh

@@ -22,10 +22,10 @@ solverName="cfdemSolverPiso"
 nrProcs="4"
 machineFileName="none"   # yourMachinefileName | none
 debugMode="off"          # on | off| strict
-reconstuctCase="true"    # true | false
+reconstructCase="true"   # true | false
 testHarnessPath="$CFDEM_TEST_HARNESS_PATH"
 #--------------------------------------------------------------------------------#
 
 #- call function to run a parallel CFD-DEM case
-parCFDDEMrun $logpath $logfileName $casePath $headerText $solverName $nrProcs $machineFileName $debugMode $reconstuctCase
+parCFDDEMrun $logpath $logfileName $casePath $headerText $solverName $nrProcs $machineFileName $debugMode $reconstructCase
 

tutorials/cfdemSolverPiso/settlingTestMPI/parCFDDEMrun.sh

@@ -22,6 +22,7 @@ solverName="cfdemSolverPiso"
 nrProcs="2"
 machineFileName="none"   # yourMachinefileName | none
 debugMode="off"          # on | off| strict
+reconstructCase="true"   # true | false
 testHarnessPath="$CFDEM_TEST_HARNESS_PATH"
 runOctave="true"
 cleanUp="true"
@@ -29,7 +30,7 @@ postproc="false"
 #--------------------------------------------------------------------------------#
 
 #- call function to run a parallel CFD-DEM case
-parCFDDEMrun $logpath $logfileName $casePath $headerText $solverName $nrProcs $machineFileName $debugMode "true"
+parCFDDEMrun $logpath $logfileName $casePath $headerText $solverName $nrProcs $machineFileName $debugMode $reconstructCase
 
 if [ $runOctave == "true" ]
     then