diff --git a/applications/solvers/cfdemSolverPiso/cfdemSolverPisoTemp/TEqn.H b/applications/solvers/cfdemSolverPiso/cfdemSolverPisoTemp/TEqn.H
index a551f90e..f77031d0 100644
--- a/applications/solvers/cfdemSolverPiso/cfdemSolverPisoTemp/TEqn.H
+++ b/applications/solvers/cfdemSolverPiso/cfdemSolverPisoTemp/TEqn.H
@@ -3,15 +3,16 @@
 	
         // solve temperature transport equation assuming const. density and heat capacity
 	
-	thDiff=voidfraction*particleCloud.thermCondM().thermDiff();
-	thDiff.correctBoundaryConditions();
+ 	thDiff=particleCloud.thermCondM().thermDiff();
+        // BC for thDiff are calculated from its mathematical expression
+        // to enforce user-defined values, change thDiff to MUST_READ and do
+        // thDiff.correctBoundaryConditions();
 		
         fvScalarMatrix TEqn
         (
            fvm::ddt(voidfraction,T) + fvm::div(phi, T)
          - fvm::Sp(fvc::ddt(voidfraction),T) - fvm::Sp(fvc::div(phi),T)
-       //  - fvm::laplacian(voidfraction*particleCloud.thermCondM().thermDiff(), T)
-	 - fvm::laplacian(thDiff, T)
+	 - fvm::laplacian(voidfraction*thDiff, T)
          ==
            Qsource/(rho*Cp)
         );
diff --git a/applications/solvers/cfdemSolverPiso/cfdemSolverPisoTemp/createFields.H b/applications/solvers/cfdemSolverPiso/cfdemSolverPisoTemp/createFields.H
index 65cd259a..790914b6 100644
--- a/applications/solvers/cfdemSolverPiso/cfdemSolverPisoTemp/createFields.H
+++ b/applications/solvers/cfdemSolverPiso/cfdemSolverPisoTemp/createFields.H
@@ -127,10 +127,11 @@
             "thDiff",
             runTime.timeName(),
             mesh,
-            IOobject::MUST_READ,
+            IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
-        mesh
+        mesh,
+        dimensionedScalar("zero", dimensionSet(0,2,-1,0,0,0,0), 0.0)
     );
     
     Info<< "\nCreating temperature gradient field\n" << endl;
diff --git a/applications/solvers/cfdemSolverRhoPimple/EEqn.H b/applications/solvers/cfdemSolverRhoPimple/EEqn.H
index b90025f8..f6933aff 100644
--- a/applications/solvers/cfdemSolverRhoPimple/EEqn.H
+++ b/applications/solvers/cfdemSolverRhoPimple/EEqn.H
@@ -5,8 +5,10 @@
     volScalarField& he = thermo.he();
     particleCloud.energyContributions(Qsource);
 
-    thDiff=particleCloud.thermCondM().thermDiff();
-    thDiff.correctBoundaryConditions();
+    thDiff=particleCloud.thermCondM().thermDiff();   
+    // BC for thDiff are calculated from its mathematical expression
+    // to enforce user-defined values, change thDiff to MUST_READ and do
+    // thDiff.correctBoundaryConditions();
 
     fvScalarMatrix EEqn
     (
@@ -25,7 +27,7 @@
      // net heat transfer from particles to fluid
       - Qsource
      // thermal conduction of the fluid with effective conductivity
-      - fvc::laplacian(rhoeps*thDiff,he) 
+      - fvm::laplacian(rhoeps*thDiff,he)
           
       // + particle-fluid energy transfer due to work
       // + fluid energy dissipation due to shearing
diff --git a/applications/solvers/cfdemSolverRhoPimple/cfdemSolverRhoPimpleChem/createFields.H b/applications/solvers/cfdemSolverRhoPimple/cfdemSolverRhoPimpleChem/createFields.H
index d2e5aa3b..e9655a6e 100644
--- a/applications/solvers/cfdemSolverRhoPimple/cfdemSolverRhoPimpleChem/createFields.H
+++ b/applications/solvers/cfdemSolverRhoPimple/cfdemSolverRhoPimpleChem/createFields.H
@@ -106,10 +106,11 @@
             "thDiff",
             runTime.timeName(),
             mesh,
-            IOobject::MUST_READ,
+            IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
-        mesh
+        mesh,
+        dimensionedScalar("zero", dimensionSet(0,2,-1,0,0,0,0), 0.0)
     );
 
     #ifndef compressibleCreatePhi_H
diff --git a/applications/solvers/cfdemSolverRhoPimple/createFields.H b/applications/solvers/cfdemSolverRhoPimple/createFields.H
index b9b61054..8c352c28 100644
--- a/applications/solvers/cfdemSolverRhoPimple/createFields.H
+++ b/applications/solvers/cfdemSolverRhoPimple/createFields.H
@@ -76,10 +76,11 @@ Info<< "Reading thermophysical properties\n" << endl;
             "thDiff",
             runTime.timeName(),
             mesh,
-            IOobject::MUST_READ,
+            IOobject::NO_READ,
             IOobject::AUTO_WRITE
         ),
-        mesh
+        mesh,
+        dimensionedScalar("zero", dimensionSet(0,2,-1,0,0,0,0), 0.0)
     );
 
     #ifndef compressibleCreatePhi_H