diff --git a/applications/solvers/rcfdemSolverRhoSteadyPimple/UEqn.H b/applications/solvers/rcfdemSolverRhoSteadyPimple/UEqn.H
index 3b296ac8..9e75193e 100644
--- a/applications/solvers/rcfdemSolverRhoSteadyPimple/UEqn.H
+++ b/applications/solvers/rcfdemSolverRhoSteadyPimple/UEqn.H
@@ -26,6 +26,7 @@ else
 }
 
 //U.relax();
+#include "limitU.H"
 
 fvOptions.correct(U);
 
diff --git a/applications/solvers/rcfdemSolverRhoSteadyPimple/createFields.H b/applications/solvers/rcfdemSolverRhoSteadyPimple/createFields.H
index 6429ef5f..9df05258 100644
--- a/applications/solvers/rcfdemSolverRhoSteadyPimple/createFields.H
+++ b/applications/solvers/rcfdemSolverRhoSteadyPimple/createFields.H
@@ -189,6 +189,39 @@ Info<< "Reading thermophysical properties\n" << endl;
         )
     );
 
+    dimensionedScalar pMax
+    (
+        dimensionedScalar::lookupOrDefault
+        (
+            "pMax",
+            pimple.dict(),
+            dimPressure,
+            GREAT
+        )
+    );
+
+    dimensionedScalar pMin
+    (
+        dimensionedScalar::lookupOrDefault
+        (
+            "pMin",
+            pimple.dict(),
+            dimPressure,
+            -GREAT
+        )
+    );
+
+    dimensionedScalar UMax
+    (
+        dimensionedScalar::lookupOrDefault
+        (
+            "UMax",
+            pimple.dict(),
+            dimVelocity,
+            -1.0
+        )
+    );
+
     Info<< "Creating turbulence model\n" << endl;
     autoPtr<compressible::turbulenceModel> turbulence
     (
@@ -265,4 +298,6 @@ Info<< "Reading thermophysical properties\n" << endl;
         ),
         Us
     );
+
+
 //===============================
diff --git a/applications/solvers/rcfdemSolverRhoSteadyPimple/pEqn.H b/applications/solvers/rcfdemSolverRhoSteadyPimple/pEqn.H
index b1f352f1..13323354 100644
--- a/applications/solvers/rcfdemSolverRhoSteadyPimple/pEqn.H
+++ b/applications/solvers/rcfdemSolverRhoSteadyPimple/pEqn.H
@@ -60,6 +60,8 @@ else
 // Explicitly relax pressure for momentum corrector
 p.relax();
 
+#include "limitP.H"
+
 // Recalculate density from the relaxed pressure
 rho = thermo.rho();
 rho = max(rho, rhoMin);
@@ -77,6 +79,8 @@ else
     U = HbyA - rAU*(fvc::grad(p)-Ksl*UsRec);
 }
 
+#include "limitU.H"
+
 U.correctBoundaryConditions();
 fvOptions.correct(U);
 K = 0.5*magSqr(U);
diff --git a/applications/solvers/rcfdemSolverRhoSteadyPimple/rcfdemSolverRhoSteadyPimple.C b/applications/solvers/rcfdemSolverRhoSteadyPimple/rcfdemSolverRhoSteadyPimple.C
index 65a5efb8..6d2fe6e1 100644
--- a/applications/solvers/rcfdemSolverRhoSteadyPimple/rcfdemSolverRhoSteadyPimple.C
+++ b/applications/solvers/rcfdemSolverRhoSteadyPimple/rcfdemSolverRhoSteadyPimple.C
@@ -85,6 +85,11 @@ int main(int argc, char *argv[])
     scalar recTimeStep = recurrenceBase.recM().recTimeStep();
     scalar startTime = runTime.startTime().value();
 
+    const IOdictionary& couplingProps = particleCloud.couplingProperties();
+    label nEveryFlow(couplingProps.lookupOrDefault<label>("nEveryFlow",1));
+    Info << "Solving flow equations every " << nEveryFlow << " steps.\n" << endl;
+    label stepcounter = 0;
+
     Info<< "\nStarting time loop\n" << endl;
 
     while (runTime.run())
@@ -127,33 +132,36 @@ int main(int argc, char *argv[])
 
         particleCloud.clockM().start(26,"Flow");
         volScalarField rhoeps("rhoeps",rho*voidfractionRec);
-        while (pimple.loop())
+        if (stepcounter%nEveryFlow==0)
         {
-            // if needed, perform drag update here
-            if (pimple.nCorrPIMPLE() <= 1)
+            while (pimple.loop())
             {
-                #include "rhoEqn.H"
-            }
+                // if needed, perform drag update here
+                if (pimple.nCorrPIMPLE() <= 1)
+                {
+                    #include "rhoEqn.H"
+                }
 
+                // --- Pressure-velocity PIMPLE corrector loop
 
-            // --- Pressure-velocity PIMPLE corrector loop
+                #include "UEqn.H"
+                #include "EEqn.H"
 
-            #include "UEqn.H"
-            #include "EEqn.H"
+                // --- Pressure corrector loop
+                while (pimple.correct())
+                {
+                    // besides this pEqn, OF offers a "pimple consistent"-option
+                    #include "pEqn.H"
+                    rhoeps=rho*voidfractionRec;
+                }
 
-            // --- Pressure corrector loop
-            while (pimple.correct())
-            {
-                // besides this pEqn, OF offers a "pimple consistent"-option
-                #include "pEqn.H"
-                rhoeps=rho*voidfractionRec;
-            }
-
-            if (pimple.turbCorr())
-            {
-                turbulence->correct();
+                if (pimple.turbCorr())
+                {
+                    turbulence->correct();
+                }
             }
         }
+        stepcounter++;
         particleCloud.clockM().stop("Flow");
 
         particleCloud.clockM().start(31,"postFlow");
diff --git a/src/lagrangian/cfdemParticle/Make/files b/src/lagrangian/cfdemParticle/Make/files
index 2f572ada..22cad185 100644
--- a/src/lagrangian/cfdemParticle/Make/files
+++ b/src/lagrangian/cfdemParticle/Make/files
@@ -41,7 +41,7 @@ $(chemistryModels)/reactantPerParticle/reactantPerParticle.C
 $(energyModels)/energyModel/energyModel.C
 $(energyModels)/energyModel/newEnergyModel.C
 $(energyModels)/heatTransferGunn/heatTransferGunn.C
-$(energyModels)/heatTransferGunnPartField/heatTransferGunnPartField.C
+$(energyModels)/heatTransferGranConduction/heatTransferGranConduction.C
 $(energyModels)/reactionHeat/reactionHeat.C
 
 $(thermCondModels)/thermCondModel/thermCondModel.C
diff --git a/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGranConduction/heatTransferGranConduction.C b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGranConduction/heatTransferGranConduction.C
new file mode 100644
index 00000000..0669fdb3
--- /dev/null
+++ b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGranConduction/heatTransferGranConduction.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 "heatTransferGranConduction.H"
+#include "addToRunTimeSelectionTable.H"
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+// * * * * * * * * * * * * * * Static Data Members * * * * * * * * * * * * * //
+
+defineTypeNameAndDebug(heatTransferGranConduction, 0);
+
+addToRunTimeSelectionTable(energyModel, heatTransferGranConduction, dictionary);
+
+// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //
+
+// Construct from components
+heatTransferGranConduction::heatTransferGranConduction
+(
+    const dictionary& dict,
+    cfdemCloudEnergy& sm
+)
+:
+    energyModel(dict,sm),
+    propsDict_(dict.subDict(typeName + "Props")),
+    multiTypes_(false),
+    verbose_(propsDict_.lookupOrDefault<bool>("verbose",false)),
+    QPartPartName_(propsDict_.lookupOrDefault<word>("QPartPartName","QPartPart")),
+    QPartPart_
+    (   IOobject
+        (
+            QPartPartName_,
+            sm.mesh().time().timeName(),
+            sm.mesh(),
+            IOobject::READ_IF_PRESENT,
+            IOobject::NO_WRITE
+        ),
+        sm.mesh(),
+        dimensionedScalar("zero", dimensionSet(1,-1,-3,0,0,0,0), 0.0)
+    ),
+    partEffThermCondField_
+    (
+        IOobject
+        (
+            "partEffThermCondField",
+            sm.mesh().time().timeName(),
+            sm.mesh(),
+            IOobject::READ_IF_PRESENT,
+            IOobject::NO_WRITE
+        ),
+        sm.mesh(),
+        dimensionedScalar("one", dimensionSet(1, 1, -3, -1,0,0,0), 1.0),
+	"zeroGradient"
+    ),
+    partThermCondField_
+    (
+        IOobject
+        (
+            "partThermCondField",
+            sm.mesh().time().timeName(),
+            sm.mesh(),
+            IOobject::READ_IF_PRESENT,
+            IOobject::NO_WRITE
+        ),
+        sm.mesh(),
+        dimensionedScalar("one", dimensionSet(1, 1, -3, -1,0,0,0), 1.0),
+	"zeroGradient"
+    ),
+    partTempField_(sm.mesh().lookupObject<volScalarField>("partTemp")),
+    voidfractionFieldName_(propsDict_.lookupOrDefault<word>("voidfractionFieldName","voidfraction")),
+    voidfraction_(sm.mesh().lookupObject<volScalarField> (voidfractionFieldName_)),
+    partHeatFluxName_(propsDict_.lookupOrDefault<word>("partHeatFluxName","conductiveHeatFlux")),
+    partHeatFlux_(NULL),
+    typePartThermCond_(propsDict_.lookupOrDefault<scalarList>("thermalConductivities",scalarList(1,-1.0))),
+    partThermCond_(NULL)
+{
+    allocateMyArrays();
+
+    if (typePartThermCond_[0] < 0.0)
+    {
+        FatalError << "heatTransferGranConduction: provide list of thermal conductivities." << abort(FatalError);
+    }
+
+    if (typePartThermCond_.size() > 1)
+    {
+        multiTypes_ = true;
+    }
+
+    if (multiTypes_ = true && !particleCloud_.getParticleTypes())
+    {
+        FatalError << "heatTransferGranConduction needs data for more than one type, but types are not communicated." << abort(FatalError);
+    }
+
+    if (verbose_)
+    {
+        QPartPart_.writeOpt() = IOobject::AUTO_WRITE;
+        QPartPart_.write();
+    }
+}
+
+
+// * * * * * * * * * * * * * * * * Destructor  * * * * * * * * * * * * * * * //
+
+heatTransferGranConduction::~heatTransferGranConduction()
+{
+    particleCloud_.dataExchangeM().destroy(partHeatFlux_,1);
+    particleCloud_.dataExchangeM().destroy(partThermCond_,1);
+}
+
+// * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
+void heatTransferGranConduction::allocateMyArrays() const
+{
+    // get memory for 2d arrays
+    double initVal=0.0;
+    particleCloud_.dataExchangeM().allocateArray(partHeatFlux_,initVal,1);
+    particleCloud_.dataExchangeM().allocateArray(partThermCond_,initVal,1);
+}
+// * * * * * * * * * * * * * * * * Member Fct  * * * * * * * * * * * * * * * //
+
+void heatTransferGranConduction::calcEnergyContribution()
+{
+   // realloc the arrays
+    allocateMyArrays();
+
+    calcPartEffThermCond();
+
+    QPartPart_ = fvc::laplacian(partEffThermCondField_,partTempField_);
+
+    label cellI=0;
+    scalar partVolume(0);
+    scalar QPartPart(0);
+    scalar voidfraction(1);
+
+    for(int index = 0;index < particleCloud_.numberOfParticles(); ++index)
+    {
+        cellI = particleCloud_.cellIDs()[index][0];
+        if(cellI >= 0)
+        {
+            voidfraction = voidfraction_[cellI];
+            if (voidfraction < 0.01)
+                voidfraction = 0.01;
+
+            partVolume = particleCloud_.particleVolume(index);
+            QPartPart = QPartPart_[cellI];
+
+            heatFlux(index, partVolume, voidfraction, QPartPart);
+        }
+    }
+}
+
+void heatTransferGranConduction::calcPartEffThermCond()
+{
+    calcPartThermCond();
+
+    scalar volFrac = 0.0;
+
+    forAll(partEffThermCondField_, cellI)
+    {
+        volFrac = 1.0 - voidfraction_[cellI];
+        if (volFrac < 0.334)
+        {
+            partEffThermCondField_[cellI] = 0.0;
+        }
+        else
+        {
+            partEffThermCondField_[cellI] = 0.5 * (3*volFrac - 1) * partThermCondField_[cellI];
+        }
+    }
+}
+
+void heatTransferGranConduction::calcPartThermCond()
+{
+    label cellI=0;
+    label partType = 1;
+    for(int index = 0;index < particleCloud_.numberOfParticles(); ++index)
+    {
+            cellI = particleCloud_.cellIDs()[index][0];
+            if (cellI >= 0)
+            {
+                if (multiTypes_)
+                {
+                    partType = particleCloud_.particleType(index);
+                }
+                // LIGGGGHTS counts types 1, 2, ..., C++ array starts at 0
+                partThermCond_[index][0] = typePartThermCond_[partType - 1];
+            }
+    }
+
+    partThermCondField_.primitiveFieldRef() = 0.0;
+    particleCloud_.averagingM().resetWeightFields();
+    particleCloud_.averagingM().setScalarAverage
+    (
+        partThermCondField_,
+        partThermCond_,
+        particleCloud_.particleWeights(),
+        particleCloud_.averagingM().UsWeightField(),
+        NULL
+    );
+}
+
+void heatTransferGranConduction::heatFlux(label index, scalar vol, scalar voidfraction, scalar QPartPart)
+{
+        partHeatFlux_[index][0] = vol / (1.0 - voidfraction) * QPartPart;
+}
+
+void heatTransferGranConduction::giveData()
+{
+    Info << "total conductive particle-particle heat flux [W] (Eulerian) = " << gSum(QPartPart_*1.0*QPartPart_.mesh().V()) << endl;
+
+    particleCloud_.dataExchangeM().giveData(partHeatFluxName_,"scalar-atom", partHeatFlux_);
+}
+
+void heatTransferGranConduction::postFlow()
+{
+    giveData();
+}
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// ************************************************************************* //
+
diff --git a/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGranConduction/heatTransferGranConduction.H b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGranConduction/heatTransferGranConduction.H
new file mode 100644
index 00000000..1dd07050
--- /dev/null
+++ b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGranConduction/heatTransferGranConduction.H
@@ -0,0 +1,129 @@
+/*---------------------------------------------------------------------------*\
+License
+
+    This is free software: you can redistribute it and/or modify it
+    under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    This code is distributed in the hope that it will be useful, but WITHOUT
+    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+    for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this code.  If not, see <http://www.gnu.org/licenses/>.
+
+    Copyright (C) 2015- Thomas Lichtenegger, JKU Linz, Austria
+
+    Description
+    heat transfer from granular conduction according to the
+    Effective Medium Theory discussed e.g. by
+    Carson et al. Int. J. Heat Mass Transfer 48 (2005)
+
+\*---------------------------------------------------------------------------*/
+
+#ifndef heatTransferGranConduction_H
+#define heatTransferGranConduction_H
+
+#include "fvCFD.H"
+#include "cfdemCloudEnergy.H"
+#include "energyModel.H"
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+namespace Foam
+{
+
+/*---------------------------------------------------------------------------*\
+                           Class heatTransferGranConduction Declaration
+\*---------------------------------------------------------------------------*/
+
+class heatTransferGranConduction
+:
+    public energyModel
+{
+protected:
+
+    dictionary propsDict_;
+
+    bool multiTypes_;
+
+    bool verbose_;
+
+    word QPartPartName_;
+
+    volScalarField QPartPart_;
+
+    const volScalarField& partTempField_;
+
+    volScalarField partEffThermCondField_;
+
+    volScalarField partThermCondField_;
+
+    word voidfractionFieldName_;
+
+    const volScalarField& voidfraction_;
+  
+    word partHeatFluxName_;
+
+    mutable double **partHeatFlux_;
+
+    scalarList typePartThermCond_;
+
+    mutable double **partThermCond_;
+
+    void allocateMyArrays() const;
+
+    void calcPartEffThermCond();
+
+    void calcPartThermCond();
+   
+    virtual void giveData();
+
+    virtual void heatFlux(label, scalar, scalar, scalar);
+
+
+
+public:
+
+    //- Runtime type information
+    TypeName("heatTransferGranConduction");
+
+    // Constructors
+
+        //- Construct from components
+        heatTransferGranConduction
+        (
+            const dictionary& dict,
+            cfdemCloudEnergy& sm
+        );
+
+
+    // Destructor
+
+        virtual ~heatTransferGranConduction();
+
+
+    // Member Functions
+
+        void addEnergyContribution(volScalarField&) const {}
+
+        void addEnergyCoefficient(volScalarField&) const {}
+
+	void calcEnergyContribution();
+
+	void postFlow();
+
+};
+
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+} // End namespace Foam
+
+// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
+
+#endif
+
+// ************************************************************************* //
+
diff --git a/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunn/heatTransferGunn.C b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunn/heatTransferGunn.C
index 779a54bf..2760792e 100644
--- a/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunn/heatTransferGunn.C
+++ b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunn/heatTransferGunn.C
@@ -43,6 +43,7 @@ heatTransferGunn::heatTransferGunn
 :
     energyModel(dict,sm),
     propsDict_(dict.subDict(typeName + "Props")),
+    multiTypes_(false),
     expNusselt_(propsDict_.lookupOrDefault<bool>("expNusselt",false)),
     interpolation_(propsDict_.lookupOrDefault<bool>("interpolation",false)),
     verbose_(propsDict_.lookupOrDefault<bool>("verbose",false)),
@@ -58,8 +59,7 @@ heatTransferGunn::heatTransferGunn
             IOobject::AUTO_WRITE
         ),
         sm.mesh(),
-        dimensionedScalar("zero", dimensionSet(1,-1,-3,0,0,0,0), 0.0),
-	"zeroGradient"
+        dimensionedScalar("zero", dimensionSet(1,-1,-3,0,0,0,0), 0.0)
     ),
     QPartFluidCoeffName_(propsDict_.lookupOrDefault<word>("QPartFluidCoeffName","QPartFluidCoeff")),
     QPartFluidCoeff_
@@ -85,7 +85,7 @@ heatTransferGunn::heatTransferGunn
         ),
         sm.mesh(),
         dimensionedScalar("zero", dimensionSet(0,0,0,1,0,0,0), 0.0),
-	"zeroGradient"
+        "zeroGradient"
     ),
     partRelTempField_
     (   IOobject
@@ -97,8 +97,7 @@ heatTransferGunn::heatTransferGunn
             IOobject::NO_WRITE
         ),
         sm.mesh(),
-        dimensionedScalar("zero", dimensionSet(0,0,0,0,0,0,0), 0.0),
-	"zeroGradient"
+        dimensionedScalar("zero", dimensionSet(0,0,0,0,0,0,0), 0.0)
     ),
     ReField_
     (   IOobject
@@ -110,22 +109,19 @@ heatTransferGunn::heatTransferGunn
             IOobject::NO_WRITE
         ),
         sm.mesh(),
-        dimensionedScalar("zero", dimensionSet(0,0,0,0,0,0,0), 0.0),
-	"zeroGradient"
+        dimensionedScalar("zero", dimensionSet(0,0,0,0,0,0,0), 0.0)
     ),
-    NusseltFieldName_(propsDict_.lookupOrDefault<word>("NusseltFieldName","NuField")),
     NuField_
     (   IOobject
         (
-            NusseltFieldName_,
+            "NuField",
             sm.mesh().time().timeName(),
             sm.mesh(),
-            IOobject::READ_IF_PRESENT,
+            IOobject::NO_READ,
             IOobject::NO_WRITE
         ),
         sm.mesh(),
-        dimensionedScalar("zero", dimensionSet(0,0,0,0,0,0,0), 0.0),
-	"zeroGradient"
+        dimensionedScalar("zero", dimensionSet(0,0,0,0,0,0,0), 0.0)
     ),
     partRefTemp_("partRefTemp", dimensionSet(0,0,0,1,0,0,0), 0.0),
     calcPartTempField_(propsDict_.lookupOrDefault<bool>("calcPartTempField",false)),
@@ -142,11 +138,13 @@ heatTransferGunn::heatTransferGunn
     rho_(sm.mesh().lookupObject<volScalarField> (densityFieldName_)),
     partTempName_(propsDict_.lookup("partTempName")),
     partTemp_(NULL),
-    partHeatFluxName_(propsDict_.lookup("partHeatFluxName")),
+    partHeatFluxName_(propsDict_.lookupOrDefault<word>("partHeatFluxName","convectiveHeatFlux")),
     partHeatFlux_(NULL),
     partHeatFluxCoeff_(NULL),
     partRe_(NULL),
-    partNu_(NULL)
+    partNu_(NULL),
+    scaleDia_(1.),
+    typeCG_(propsDict_.lookupOrDefault<scalarList>("coarseGrainingFactors",scalarList(1,1.0)))
 {
     allocateMyArrays();
 
@@ -186,13 +184,14 @@ heatTransferGunn::heatTransferGunn
           FatalError <<"Cannot read and create NuField at the same time!\n" << abort(FatalError);
         }
     }
-    
-    if (expNusselt_)
+
+    if (propsDict_.found("scale") && typeCG_.size()==1)
     {
-        NuField_.writeOpt() = IOobject::AUTO_WRITE;
-        NuField_.write();
-        Info <<  "Using predefined Nusselt number field." << endl;
+        // if "scale" is specified and there's only one single type, use "scale"
+        scaleDia_=scalar(readScalar(propsDict_.lookup("scale")));
+        typeCG_[0] = scaleDia_;
     }
+    else if (typeCG_.size()>1) multiTypes_ = true;
 }
 
 
@@ -229,68 +228,7 @@ 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()
 {
@@ -301,9 +239,26 @@ void heatTransferGunn::calcEnergyContribution()
     QPartFluid_.primitiveFieldRef() = 0.0;
 
     // get DEM data
-    particleCloud_.clockM().start(29,"getDEM_Tdata");
     particleCloud_.dataExchangeM().getData(partTempName_,"scalar-atom",partTemp_);
-    particleCloud_.clockM().stop("getDEM_Tdata");
+
+    if(calcPartTempField_)
+    {
+        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());
+        partRelTempField_ = (partTempField_ - aveTemp) / (aveTemp - partRefTemp_);
+        Info << "heatTransferGunn: average part. temp = " << aveTemp.value() << endl;
+    }
 
     #ifdef compre
        const volScalarField mufField = particleCloud_.turbulence().mu();
@@ -311,6 +266,15 @@ void heatTransferGunn::calcEnergyContribution()
        const volScalarField mufField = particleCloud_.turbulence().nu()*rho_;
     #endif
 
+    if (typeCG_.size()>1 || typeCG_[0] > 1)
+    {
+        Info << "heatTransferGunn using scale = " << typeCG_ << endl;
+    }
+    else if (particleCloud_.cg() > 1)
+    {
+        scaleDia_=particleCloud_.cg();
+        Info << "heatTransferGunn using scale from liggghts cg = " << scaleDia_ << endl;
+    }
 
     // calc La based heat flux
     scalar voidfraction(1);
@@ -319,14 +283,17 @@ void heatTransferGunn::calcEnergyContribution()
     label cellI=0;
     vector Us(0,0,0);
     scalar ds(0);
+    scalar ds_scaled(0);
+    scalar scaleDia3 = typeCG_[0]*typeCG_[0]*typeCG_[0];
     scalar muf(0);
     scalar magUr(0);
     scalar Rep(0);
     scalar Pr(0);
     scalar Nup(0);
     scalar Tsum(0.0);
-    scalar Nsum(0.0);
 
+    scalar cg = typeCG_[0];
+    label partType = 1;
 
     interpolationCellPoint<scalar> voidfractionInterpolator_(voidfraction_);
     interpolationCellPoint<vector> UInterpolator_(U_);
@@ -351,38 +318,33 @@ void heatTransferGunn::calcEnergyContribution()
                     Tfluid = tempField_[cellI];
                 }
 
-                
                 if (voidfraction < 0.01)
                     voidfraction = 0.01;
 
+                if (multiTypes_)
+                {
+                    partType = particleCloud_.particleType(index);
+                    cg = typeCG_[partType - 1];
+                    scaleDia3 = cg*cg*cg;
+                }
+
+                // calc relative velocity
+                Us = particleCloud_.velocity(index);
+                magUr = mag(Ufluid - Us);
                 ds = 2.*particleCloud_.radius(index);
+                ds_scaled = ds/cg;
+                muf = mufField[cellI];
+                Rep = ds_scaled * magUr * voidfraction * rho_[cellI]/ muf;
+                Pr = max(SMALL, Cp_ * muf / kf0_);
 
-                if (expNusselt_)
-                {
-                    Nup = NuField_[cellI];
-                    if (Nup < 2.0)
-                        Nup = 2.0;
-                }
-                else
-                {
-                    Us = particleCloud_.velocity(index);
-                    magUr = mag(Ufluid - Us);
-                
-                    muf = mufField[cellI];
-                    Rep = ds * magUr * voidfraction * rho_[cellI]/ muf;
-                    Pr = max(SMALL, Cp_ * muf / kf0_);
+                Nup = Nusselt(voidfraction, Rep, Pr);
 
-                    Nup = Nusselt(voidfraction, Rep, Pr);
-                }
-                
                 Tsum += partTemp_[index][0];
-                Nsum += 1.0;
-
-                scalar h = kf0_ * Nup / ds;
-                scalar As = ds * ds * M_PI; // surface area of sphere
+                scalar h = kf0_ * Nup / ds_scaled;
+                scalar As = ds_scaled * ds_scaled * M_PI; // surface area of sphere
 
                 // calc convective heat flux [W]
-                heatFlux(index, h, As, Tfluid);
+                heatFlux(index, h, As, Tfluid, scaleDia3);
 
                 if(verbose_)
                 {
@@ -405,16 +367,15 @@ void heatTransferGunn::calcEnergyContribution()
                 }
             }
     }
-    
+
     // gather particle temperature sums and obtain average
     if(calcPartTempAve_)
     {
         reduce(Tsum, sumOp<scalar>());
-        reduce(Nsum, sumOp<scalar>());
-        partTempAve_ = Tsum / Nsum;
+        partTempAve_ = Tsum / particleCloud_.numberOfParticles();
         Info << "mean particle temperature = " << partTempAve_ << endl;
     }
-    
+
     if(calcPartTempField_) partTempField();
 
     particleCloud_.averagingM().setScalarSum
@@ -484,6 +445,55 @@ void heatTransferGunn::calcEnergyContribution()
     QPartFluid_.correctBoundaryConditions();
 }
 
+void heatTransferGunn::addEnergyContribution(volScalarField& Qsource) const
+{
+    Qsource += QPartFluid_;
+}
+
+void heatTransferGunn::addEnergyCoefficient(volScalarField& Qsource) const
+{
+    if(implicit_)
+    {
+        Qsource += QPartFluidCoeff_;
+    }
+}
+
+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, scalar cg3)
+{
+    scalar hAs = h * As * cg3;
+
+    if (particleCloud_.getParticleEffVolFactors()) 
+    {
+        scalar effVolFac = particleCloud_.particleEffVolFactor(index);
+        hAs *= effVolFac;
+    }
+
+    partHeatFlux_[index][0] = - hAs * partTemp_[index][0];
+    if(!implicit_)
+    {
+        partHeatFlux_[index][0] += hAs * Tfluid;
+    }
+    else
+    {
+        partHeatFluxCoeff_[index][0] = hAs;
+    }
+}
+
+void heatTransferGunn::giveData()
+{
+    Info << "total convective particle-fluid heat flux [W] (Eulerian) = " << gSum(QPartFluid_*1.0*QPartFluid_.mesh().V()) << endl;
+
+    particleCloud_.dataExchangeM().giveData(partHeatFluxName_,"scalar-atom", partHeatFlux_);
+}
+
 void heatTransferGunn::postFlow()
 {
     if(implicit_)
@@ -520,6 +530,23 @@ 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
+    );
+
+    dimensionedScalar aveTemp("aveTemp",dimensionSet(0,0,0,1,0,0,0), partTempAve_);
+    partRelTempField_ = (partTempField_ - aveTemp) / (aveTemp - partRefTemp_);
+}
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 } // End namespace Foam
diff --git a/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunn/heatTransferGunn.H b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunn/heatTransferGunn.H
index fb077759..6139d3c6 100644
--- a/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunn/heatTransferGunn.H
+++ b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunn/heatTransferGunn.H
@@ -45,6 +45,8 @@ protected:
 
     dictionary propsDict_;
 
+    bool multiTypes_;
+
     bool expNusselt_;
 
     bool interpolation_;
@@ -67,16 +69,14 @@ protected:
 
     volScalarField ReField_;
 
-    word NusseltFieldName_;
-
     volScalarField NuField_;
 
     dimensionedScalar partRefTemp_;
 
     bool calcPartTempField_;
-    
+
     bool calcPartTempAve_;
-    
+
     scalar partTempAve_;
 
     word tempFieldName_;
@@ -111,6 +111,10 @@ protected:
 
     mutable double **partNu_;
 
+    mutable scalar scaleDia_;
+
+    scalarList typeCG_;
+
     void allocateMyArrays() const;
 
     void partTempField();
@@ -119,7 +123,7 @@ protected:
 
     virtual void giveData();
 
-    virtual void heatFlux(label, scalar, scalar, scalar);
+    virtual void heatFlux(label, scalar, scalar, scalar, scalar cg3 = 1.0);
 
 public:
 
diff --git a/src/lagrangian/cfdemParticleComp/Make/files b/src/lagrangian/cfdemParticleComp/Make/files
index 077b22cf..008b2935 100644
--- a/src/lagrangian/cfdemParticleComp/Make/files
+++ b/src/lagrangian/cfdemParticleComp/Make/files
@@ -39,7 +39,7 @@ $(chemistryModels)/reactantPerParticle/reactantPerParticle.C
 $(energyModels)/energyModel/energyModel.C
 $(energyModels)/energyModel/newEnergyModel.C
 $(energyModels)/heatTransferGunn/heatTransferGunn.C
-$(energyModels)/heatTransferGunnPartField/heatTransferGunnPartField.C
+$(energyModels)/heatTransferGranConduction/heatTransferGranConduction.C
 $(energyModels)/reactionHeat/reactionHeat.C
 
 $(thermCondModels)/thermCondModel/thermCondModel.C