diff --git a/UEqn.H b/UEqn.H
index e4caa1c..849fe89 100644
--- a/UEqn.H
+++ b/UEqn.H
@@ -2,11 +2,15 @@
 
 MRF.correctBoundaryVelocity(U);
 
+
 tmp<fvVectorMatrix> tUEqn
 (
     fvm::ddt(rho, U) + fvm::div(phi, U)
   + MRF.DDt(rho, U)
-  + turbulence->divDevRhoReff(U)
+  + (
+     - fvc::div(diff.mu()*Foam::dev2(Foam::T(fvc::grad(U))))
+     - fvm::laplacian(diff.mu(), U)
+    )
  ==
     fvOptions(rho, U)
 );
diff --git a/diffusivityModel/diffusivityModel/diffusivityModel.C b/diffusivityModel/diffusivityModel/diffusivityModel.C
index 43f45c3..3c9176d 100644
--- a/diffusivityModel/diffusivityModel/diffusivityModel.C
+++ b/diffusivityModel/diffusivityModel/diffusivityModel.C
@@ -287,6 +287,32 @@ Foam::diffusivityModel::diffusivityModel(const psiReactionThermo& thermo)
 :
     thermo_(thermo),
     D_(thermo_.composition().species().size()),
+    mu_
+    (
+        IOobject
+        (
+            "mu",
+            thermo_.composition().Y(0).mesh().time().timeName(),
+            thermo_.composition().Y(0).mesh(),
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        thermo_.composition().Y(0).mesh(),
+        dimensionedScalar("zero", dimPressure*dimTime, 0.0)
+    ),
+    k_
+    (
+        IOobject
+        (
+            "lambda",
+            thermo_.composition().Y(0).mesh().time().timeName(),
+            thermo_.composition().Y(0).mesh(),
+            IOobject::NO_READ,
+            IOobject::AUTO_WRITE
+        ),
+        thermo_.composition().Y(0).mesh(),
+        dimensionedScalar("zero", dimArea/dimTime, 0.0)
+    ),
     neutrals_(thermo_.composition().species().size()),
     ions_(thermo_.composition().species().size()),
     ecs_(thermo_.composition().species().size()),
@@ -481,40 +507,6 @@ Foam::diffusivityModel::diffusivityModel(const psiReactionThermo& thermo)
         }
     }
 
-    mu_.set
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                "mu",
-                mesh.time().timeName(),
-                mesh,
-                IOobject::NO_READ,
-                IOobject::AUTO_WRITE
-            ),
-            mesh,
-            dimensionedScalar("zero", dimArea/dimTime, 0.0)
-        )
-    );
-
-    k_.set
-    (
-        new volScalarField
-        (
-            IOobject
-            (
-                "lambda",
-                mesh.time().timeName(),
-                mesh,
-                IOobject::NO_READ,
-                IOobject::AUTO_WRITE
-            ),
-            mesh,
-            dimensionedScalar("zero", dimArea/dimTime, 0.0)
-        )
-    );
-
     forAll(species_, i)
     {
         D_.set
@@ -540,7 +532,9 @@ Foam::diffusivityModel::diffusivityModel(const psiReactionThermo& thermo)
 
 Foam::diffusivityModel::diffusivityModel(const diffusivityModel& dm)
 :
-    thermo_(dm.thermo_)
+    thermo_(dm.thermo_),
+    mu_(dm.mu_),
+    k_(dm.k_)
 {}
 
 
@@ -678,8 +672,8 @@ void Foam::diffusivityModel::correct()
             D_[i][celli] = Di[i];
         }
 
-        mu_()[celli] = mixAvgMu(muI, localX, Wpure);
-        k_()[celli] = mixAvgK(kI, localX) / rhoi / Cpi;
+        mu_[celli] = mixAvgMu(muI, localX, Wpure);
+        k_[celli] = mixAvgK(kI, localX) / rhoi / Cpi;
     }
 
 
@@ -755,8 +749,8 @@ void Foam::diffusivityModel::correct()
                 D_[i].boundaryFieldRef()[patchi][facei] = Di[i];
             }
 
-            mu_().boundaryFieldRef()[patchi][facei] = mixAvgMu(muI, localX, Wpure);
-            k_().boundaryFieldRef()[patchi][facei] = mixAvgK(kI, localX) / rhoi / Cpi;
+            mu_.boundaryFieldRef()[patchi][facei] = mixAvgMu(muI, localX, Wpure);
+            k_.boundaryFieldRef()[patchi][facei] = mixAvgK(kI, localX) / rhoi / Cpi;
         }
     }
 
diff --git a/diffusivityModel/diffusivityModel/diffusivityModel.H b/diffusivityModel/diffusivityModel/diffusivityModel.H
index b866fd0..1104c60 100644
--- a/diffusivityModel/diffusivityModel/diffusivityModel.H
+++ b/diffusivityModel/diffusivityModel/diffusivityModel.H
@@ -73,10 +73,10 @@ class diffusivityModel
         PtrList<volScalarField> D_;
 
         //- Mixture averaged viscosity field
-        autoPtr<volScalarField> mu_;
+        volScalarField mu_;
 
         //- Mixture averaged thermal conductivity field
-        autoPtr<volScalarField> k_;
+        volScalarField k_;
 
         //- Electron object
         autoPtr<Electron> electron_;
@@ -174,6 +174,15 @@ public:
 
         // Access
 
+        inline volScalarField& mu();
+
+        inline const volScalarField& mu() const;
+
+        inline volScalarField& k();
+
+        inline const volScalarField& k() const;
+
+
         // Check
 
         // Edit
diff --git a/diffusivityModel/diffusivityModel/diffusivityModelI.H b/diffusivityModel/diffusivityModel/diffusivityModelI.H
index eca8d9e..2edf73c 100644
--- a/diffusivityModel/diffusivityModel/diffusivityModelI.H
+++ b/diffusivityModel/diffusivityModel/diffusivityModelI.H
@@ -37,6 +37,30 @@ License
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 
+Foam::volScalarField& Foam::diffusivityModel::mu()
+{
+    return mu_;
+}
+
+
+const Foam::volScalarField& Foam::diffusivityModel::mu() const
+{
+    return mu_;
+}
+
+
+Foam::volScalarField& Foam::diffusivityModel::k()
+{
+    return k_;
+}
+
+
+const Foam::volScalarField& Foam::diffusivityModel::k() const
+{
+    return k_;
+}
+
+
 // * * * * * * * * * * * * * * * Member Operators  * * * * * * * * * * * * * //