diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/turbulentDispersion/turbulentDispersion.C b/src/lagrangian/cfdemParticle/subModels/forceModel/turbulentDispersion/turbulentDispersion.C
index 62afa5f8..cb7cb4a8 100644
--- a/src/lagrangian/cfdemParticle/subModels/forceModel/turbulentDispersion/turbulentDispersion.C
+++ b/src/lagrangian/cfdemParticle/subModels/forceModel/turbulentDispersion/turbulentDispersion.C
@@ -74,13 +74,16 @@ turbulentDispersion::turbulentDispersion
         sm.mesh(),
         dimensionedScalar("zero", dimensionSet(0,0,0,0,0,0,0), 0.0)
     ),
-    minTurbKinetcEnergy_(propsDict_.lookupOrDefault<scalar>("minTurbKinetcEnergy", 0.0)),
+    minTurbKineticEnergy_(propsDict_.lookupOrDefault<scalar>("minTurbKineticEnergy", 0.0)),
+    turbKineticEnergyFieldName_(propsDict_.lookupOrDefault<word>("turbKineticEnergyFieldName","")),
+    turbKineticEnergy_(NULL),
+    existTurbKineticEnergyInObjReg_(false),
     voidfractionFieldName_(propsDict_.lookupOrDefault<word>("voidfractionFieldName","voidfraction")),
     voidfraction_(sm.mesh().lookupObject<volScalarField> (voidfractionFieldName_)),
     critVoidfraction_(propsDict_.lookupOrDefault<scalar>("critVoidfraction", 0.9)),
     ranGen_(clock::getTime()+pid())
 {
-    if(ignoreCellsName_ != "none")
+    if (ignoreCellsName_ != "none")
     {
        ignoreCells_.set(new cellSet(particleCloud_.mesh(),ignoreCellsName_));
        Info << type << ": ignoring fluctuations in cellSet " << ignoreCells_().name() <<
@@ -90,7 +93,7 @@ turbulentDispersion::turbulentDispersion
 
     // define a field to indicate if a cell is next to boundary
      label cellI = -1;
-     forAll(mesh_.boundary(),patchI)
+     forAll (mesh_.boundary(),patchI)
      {
          word patchName = mesh_.boundary()[patchI].name();
          if (patchName.rfind("procB",0) == 0) continue;
@@ -102,6 +105,29 @@ turbulentDispersion::turbulentDispersion
          }
      }
 
+    if (turbKineticEnergyFieldName_ != "")
+    {
+        existTurbKineticEnergyInObjReg_ = true;
+        volScalarField& k(const_cast<volScalarField&>(sm.mesh().lookupObject<volScalarField> (turbKineticEnergyFieldName_)));
+        turbKineticEnergy_ = &k;
+    }
+    else
+    {
+        turbKineticEnergy_ = new volScalarField
+        (
+            IOobject
+            (
+                "turbKinEnergy",
+                mesh_.time().timeName(),
+                mesh_,
+                IOobject::NO_READ,
+                IOobject::NO_WRITE
+            ),
+            mesh_,
+            dimensionedScalar("zero", dimensionSet(0,2,-2,0,0), 0)
+        );
+    }
+
     // make sure this is the last force model in list so that fluid velocity does not get overwritten
     label numLastForceModel = sm.nrForceModels();
     word lastForceModel = sm.forceModels()[numLastForceModel-1];
@@ -116,6 +142,7 @@ turbulentDispersion::turbulentDispersion
 
 turbulentDispersion::~turbulentDispersion()
 {
+    if (!existTurbKineticEnergyInObjReg_) delete turbKineticEnergy_;
 }
 
 // * * * * * * * * * * * * * * * private Member Functions  * * * * * * * * * * * * * //
@@ -130,7 +157,10 @@ bool turbulentDispersion::ignoreCell(label cell) const
 
 void turbulentDispersion::setForce() const
 {
-    const volScalarField turbKinetcEnergy(particleCloud_.turbulence().k());
+    if (!existTurbKineticEnergyInObjReg_)
+    {
+        *turbKineticEnergy_ = particleCloud_.turbulence().k()();
+    }
 
     label cellI = -1;
     label patchID = -1;
@@ -142,7 +172,7 @@ void turbulentDispersion::setForce() const
     vector position = vector::zero;
     word patchName("");
 
-    interpolationCellPoint<scalar> turbKinetcEnergyInterpolator_(turbKinetcEnergy);
+    interpolationCellPoint<scalar> turbKineticEnergyInterpolator_(*turbKineticEnergy_);
 
     for(int index = 0;index <  particleCloud_.numberOfParticles(); ++index)
     {
@@ -156,14 +186,14 @@ void turbulentDispersion::setForce() const
                 if (interpolate_)
                 {
                     position = particleCloud_.position(index);
-                    k = turbKinetcEnergyInterpolator_.interpolate(position,cellI);
+                    k = turbKineticEnergyInterpolator_.interpolate(position,cellI);
                 }
                 else
                 {
-                    k = turbKinetcEnergy[cellI];
+                    k = (*turbKineticEnergy_)[cellI];
                 }
 
-                if (k < minTurbKinetcEnergy_) k = minTurbKinetcEnergy_;
+                if (k < minTurbKineticEnergy_) k = minTurbKineticEnergy_;
 
                 flucU=unitFlucDir()*Foam::sqrt(2.0*k);
 
diff --git a/src/lagrangian/cfdemParticle/subModels/forceModel/turbulentDispersion/turbulentDispersion.H b/src/lagrangian/cfdemParticle/subModels/forceModel/turbulentDispersion/turbulentDispersion.H
index 4ca70eb8..8721e73b 100644
--- a/src/lagrangian/cfdemParticle/subModels/forceModel/turbulentDispersion/turbulentDispersion.H
+++ b/src/lagrangian/cfdemParticle/subModels/forceModel/turbulentDispersion/turbulentDispersion.H
@@ -60,7 +60,13 @@ protected:
 
     mutable volScalarField wallIndicatorField_;
 
-    scalar minTurbKinetcEnergy_;
+    scalar minTurbKineticEnergy_;
+
+    word turbKineticEnergyFieldName_;
+
+    volScalarField* turbKineticEnergy_;
+
+    bool existTurbKineticEnergyInObjReg_;
 
     word voidfractionFieldName_;