diff --git a/src/lagrangian/cfdemParticle/subModels/energyModel/YagiWallHT/YagiWallHT.C b/src/lagrangian/cfdemParticle/subModels/energyModel/YagiWallHT/YagiWallHT.C
index 5dbb605a..9f088239 100644
--- a/src/lagrangian/cfdemParticle/subModels/energyModel/YagiWallHT/YagiWallHT.C
+++ b/src/lagrangian/cfdemParticle/subModels/energyModel/YagiWallHT/YagiWallHT.C
@@ -103,17 +103,15 @@ namespace Foam
   velFieldName_(propsDict_.lookupOrDefault<word>("velFieldName","U")),
   U_(sm.mesh().lookupObject<volVectorField> (velFieldName_)),
   UsFieldName_(propsDict_.lookup("granVelFieldName")),
-	Us_(sm.mesh().lookupObject<volVectorField> (UsFieldName_)),
+  Us_(sm.mesh().lookupObject<volVectorField> (UsFieldName_)),
   densityFieldName_(propsDict_.lookupOrDefault<word>("densityFieldName","rho")),
   rho_(sm.mesh().lookupObject<volScalarField> (densityFieldName_)),
   WallTempName_(propsDict_.lookup("WallTempName")),
   WallTemp_(sm.mesh().lookupObject<volScalarField> (WallTempName_)),
   partRe_(NULL),
   multiphase_(propsDict_.lookupOrDefault<bool>("multiphase",false)),
-
   kfFieldName_(propsDict_.lookupOrDefault<word>("kfFieldName",voidfractionFieldName_)), // use voidfractionField as dummy to prevent lookup error when not using multiphase
   kfField_(sm.mesh().lookupObject<volScalarField> (kfFieldName_)),
-
   CpFieldName_(propsDict_.lookupOrDefault<word>("CpFieldName",voidfractionFieldName_)), // use voidfractionField as dummy to prevent lookup error when not using multiphase
   CpField_(sm.mesh().lookupObject<volScalarField> (CpFieldName_))
   {
@@ -169,8 +167,6 @@ namespace Foam
     const volScalarField mufField = particleCloud_.turbulence().nu()*rho_;
     #endif
 
-
-
     // calc La based heat flux
     scalar voidfraction(1);
     vector Ufluid(0,0,0);
@@ -180,7 +176,6 @@ namespace Foam
     scalar muf(0);
     scalar magUr(0);
     scalar Rep(0);
-    scalar H(0);
 
     interpolationCellPoint<scalar> voidfractionInterpolator_(voidfraction_);
     interpolationCellPoint<vector> UInterpolator_(U_);
@@ -239,30 +234,47 @@ namespace Foam
 
             if (WallTemp_.boundaryField().types()[patchi] == "fixedValue")
             {
-                forAll(curPatch, facei)
-                {
-                    label faceCelli = curPatch.faceCells()[facei];
-                    // calculate Urel
-                    scalar magG = mag(U_[faceCelli]-Us_[faceCelli])*voidfraction_[faceCelli]*rho_[faceCelli];
-                    // calculate H
-                    //H = 0.2087 * (pow(ReField_[faceCelli] , -0.20)) * CpField_[faceCelli] * magG / (pow(PrField_[faceCelli] , (2/3)));
-                    H = 1;
-                    // get delta T
-                    //scalar deltaT = WallTemp_.boundaryField()[patchi][facei] - tempField_.boundaryField()[patchi][facei];
-                    scalar deltaT = 100;
-                    // get facei.area
-                    scalar area = curPatch.magSf()[facei];
-                    scalar volume = U_.mesh().V()[faceCelli];
+	        if(tempField_.boundaryField().types()[patchi] == "fixedGradient")
+		{
+		    fixedGradientFvPatchField<scalar>& tempGradField(refCast<fixedGradientFvPatchField<scalar>>(tempField_.boundaryFieldRef()[patchi]));
 
-                    //Info << "faceID: " << facei << " area: " << area << " volume: " << volume << endl;
+		    forAll(curPatch, facei)
+		    {							
+			label faceCelli = curPatch.faceCells()[facei];
+			// calculate Urel
+			scalar magG = mag(U_[faceCelli]-Us_[faceCelli])*voidfraction_[faceCelli]*rho_[faceCelli];
 
-                    QWallFluid_[faceCelli] = H*deltaT*area/volume;
-                }
+			// calculate H
+			scalar H = 0.2087 * (pow(ReField_[faceCelli]+SMALL, -0.20)) * CpField_[faceCelli] * magG / (pow(PrField_[faceCelli] , (2/3))+SMALL);
+
+			// get delta T
+			scalar deltaT = WallTemp_.boundaryField()[patchi][facei] - tempField_.boundaryField()[patchi][facei];
+
+			// calculate heat flux
+			scalar heatFlux = H*deltaT;
+			scalar TGrad = -heatFlux/kfField_[faceCelli];
+			tempGradField.gradient()[facei] = TGrad;
+
+			Info << "####################" << endl;
+			Info << "G : " << magG << endl;
+			Info << "Re: " << ReField_[faceCelli] << endl;
+			Info << "Pr: " << PrField_[faceCelli] << endl;
+			Info << "Cp: " << CpField_[faceCelli] << endl;
+			Info << "kf: " << kfField_[faceCelli] << endl;
+			Info << "H : " << H << endl;
+			Info << "dT: " << deltaT << endl;
+			Info << "q': " << heatFlux << endl;
+			Info << "gradT: " << TGrad << endl;
+		    }
+		}
+		else
+		{
+		    FatalError << "YagiWallHT requires fixedGradient BC for temperature field" << endl;
+		}
+		    
             }
         }
 
-//    h = 0.2087 * (pow(partRe_ , -0.20)) * CpField_ * magUr * rho_ / (pow(PrField_ , (2/3)));
-
     // limit source term
     forAll(QWallFluid_,cellI)
     {
diff --git a/src/lagrangian/cfdemParticle/subModels/energyModel/YagiWallHT/YagiWallHT.H b/src/lagrangian/cfdemParticle/subModels/energyModel/YagiWallHT/YagiWallHT.H
index c535f538..40d25a9b 100644
--- a/src/lagrangian/cfdemParticle/subModels/energyModel/YagiWallHT/YagiWallHT.H
+++ b/src/lagrangian/cfdemParticle/subModels/energyModel/YagiWallHT/YagiWallHT.H
@@ -28,6 +28,7 @@ License
 #include "fvCFD.H"
 #include "cfdemCloudEnergy.H"
 #include "energyModel.H"
+#include "fixedGradientFvPatchFields.H"
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
 namespace Foam
@@ -61,7 +62,7 @@ protected:
 
     word tempFieldName_;
 
-    const volScalarField& tempField_;
+    volScalarField tempField_;
 
     word voidfractionFieldName_;
 
@@ -83,7 +84,7 @@ protected:
 
     word WallTempName_;
 
-    const volScalarField& WallTemp_;                          // Lagrangian array
+    const volScalarField& WallTemp_;
 
     mutable double **partRe_;