diff --git a/applications/solvers/cfdemSolverRhoPimple/EEqn.H b/applications/solvers/cfdemSolverRhoPimple/EEqn.H
index fc770162..0f1e6aff 100644
--- a/applications/solvers/cfdemSolverRhoPimple/EEqn.H
+++ b/applications/solvers/cfdemSolverRhoPimple/EEqn.H
@@ -39,7 +39,7 @@
      ==
         fvOptions(rho, he)
     );
-    
+
 
     EEqn.relax();
 
@@ -52,4 +52,9 @@
     thermo.correct();
 
     Info<< "T max/min : " << max(T).value() << " " << min(T).value() << endl;
+
+    particleCloud.clockM().start(31,"energySolve");
+    particleCloud.solve();
+    particleCloud.clockM().stop("energySolve");
+
 }
diff --git a/applications/solvers/cfdemSolverRhoSimple/UEqn.H b/applications/solvers/cfdemSolverRhoSimple/UEqn.H
index c4a58ad2..45b78abe 100644
--- a/applications/solvers/cfdemSolverRhoSimple/UEqn.H
+++ b/applications/solvers/cfdemSolverRhoSimple/UEqn.H
@@ -19,14 +19,21 @@ fvOptions.constrain(UEqn);
 if (modelType=="B" || modelType=="Bfull")
 {
     solve(UEqn == -fvc::grad(p)+ Ksl*Us);
-
-    fvOptions.correct(U);
-    K = 0.5*magSqr(U);
 }
 else
 {
     solve(UEqn == -voidfraction*fvc::grad(p)+ Ksl*Us);
-
-    fvOptions.correct(U);
-    K = 0.5*magSqr(U);
 }
+
+fvOptions.correct(U);
+
+// limit vel
+
+forAll(U,cellI)
+{
+    scalar mU = mag(U[cellI]);
+    if (mU > 200.0) U[cellI] *= 200.0/mU;
+}
+
+
+K = 0.5*magSqr(U);
diff --git a/applications/solvers/cfdemSolverRhoSimple/cfdemSolverRhoSimple.C b/applications/solvers/cfdemSolverRhoSimple/cfdemSolverRhoSimple.C
index e960b40b..b505252e 100644
--- a/applications/solvers/cfdemSolverRhoSimple/cfdemSolverRhoSimple.C
+++ b/applications/solvers/cfdemSolverRhoSimple/cfdemSolverRhoSimple.C
@@ -108,12 +108,14 @@ int main(int argc, char *argv[])
         // Pressure-velocity SIMPLE corrector
 
         #include "UEqn.H"
-        #include "EEqn.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");
diff --git a/applications/solvers/cfdemSolverRhoSimple/pEqn.H b/applications/solvers/cfdemSolverRhoSimple/pEqn.H
index a85e7f44..c032cbeb 100644
--- a/applications/solvers/cfdemSolverRhoSimple/pEqn.H
+++ b/applications/solvers/cfdemSolverRhoSimple/pEqn.H
@@ -75,6 +75,16 @@ else
 {
     U = HbyA - rAU*(fvc::grad(p)-Ksl*Us);
 }
+
+// limit vel
+
+forAll(U,cellI)
+{
+    scalar mU = mag(U[cellI]);
+    if (mU > 200.0) U[cellI] *= 200.0/mU;
+}
+
+
 U.correctBoundaryConditions();
 fvOptions.correct(U);
 K = 0.5*magSqr(U);
diff --git a/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunnPartField/heatTransferGunnPartField.C b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunnPartField/heatTransferGunnPartField.C
index 17be0c17..2dac2222 100644
--- a/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunnPartField/heatTransferGunnPartField.C
+++ b/src/lagrangian/cfdemParticle/subModels/energyModel/heatTransferGunnPartField/heatTransferGunnPartField.C
@@ -81,6 +81,8 @@ heatTransferGunnPartField::heatTransferGunnPartField
         FatalError << "heatTransferGunnPartField: provide list of specific heat capacities." << abort(FatalError);
     }
 
+    partTempField_.writeOpt() = IOobject::AUTO_WRITE;
+
     allocateMyArrays();
 }
 
@@ -173,18 +175,29 @@ void heatTransferGunnPartField::postFlow()
 
 void heatTransferGunnPartField::solve()
 {
-    Info << "patch types of partTemp boundary: " << partTempField_.boundaryField().types() << endl;
-    volScalarField Qsource = QPartFluidCoeff_*tempField_;
+    // 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(QPartFluidCoeff_, partTempField_)
+      - fvm::Sp(correctedQPartFluidCoeff, partTempField_)
       - fvm::laplacian(thCondEff,partTempField_)
      ==
         fvOptions(partCpEff, partTempField_)
@@ -199,8 +212,15 @@ void heatTransferGunnPartField::solve()
 
     partTEqn.solve();
 
+    partTempField_.relax();
+
     fvOptions.correct(partTempField_);
 
+    dimensionedScalar pTMin("pTMin",dimensionSet(0,0,0,1,0,0,0),293.0);
+    dimensionedScalar pTMax("pTMin",dimensionSet(0,0,0,1,0,0,0),3000.0);
+    partTempField_ = max(partTempField_, pTMin);
+    partTempField_ = min(partTempField_, pTMax);
+
     Info<< "partT max/min : " << max(partTempField_).value() << " " << min(partTempField_).value() << endl;
 
 }