diff --git a/YEqn.H b/YEqn.H
index 8cb5664..2552408 100644
--- a/YEqn.H
+++ b/YEqn.H
@@ -15,6 +15,14 @@ tmp<fv::convectionScheme<scalar>> mvConvection
     label inertIndex = -1;
     volScalarField Yt(0.0*Y[0]);
 
+    // volVectorField Vc(U);
+    Vc = dimensionedVector("zero", dimLength/dimTime, U[0]*0.0);
+
+    forAll(Y, i)
+    {
+        Vc += diff.D(i) * fvc::grad(Y[i]);
+    }
+
     forAll(Y, i)
     {
         if (Y[i].name() != inertSpecie)
diff --git a/createFields.H b/createFields.H
index da1cfdb..a2cb408 100644
--- a/createFields.H
+++ b/createFields.H
@@ -40,6 +40,20 @@ volVectorField U
     mesh
 );
 
+// Info<< "Reading field U\n" << endl;
+volVectorField Vc
+(
+    IOobject
+    (
+        "Vc",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    mesh,
+    dimensionedVector("zero", dimLength/dimTime, U[0]*0.0)
+);
 
 volScalarField& p = thermo.p();
 
diff --git a/diffusivityModel/diffusivityModel/diffusivityModel.H b/diffusivityModel/diffusivityModel/diffusivityModel.H
index 1104c60..ca9c128 100644
--- a/diffusivityModel/diffusivityModel/diffusivityModel.H
+++ b/diffusivityModel/diffusivityModel/diffusivityModel.H
@@ -174,6 +174,10 @@ public:
 
         // Access
 
+        inline volScalarField& D(const label i);
+
+        inline const volScalarField& D(const label i) const;
+
         inline volScalarField& mu();
 
         inline const volScalarField& mu() const;
diff --git a/diffusivityModel/diffusivityModel/diffusivityModelI.H b/diffusivityModel/diffusivityModel/diffusivityModelI.H
index 2edf73c..c0036f4 100644
--- a/diffusivityModel/diffusivityModel/diffusivityModelI.H
+++ b/diffusivityModel/diffusivityModel/diffusivityModelI.H
@@ -37,6 +37,18 @@ License
 // * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //
 
 
+Foam::volScalarField& Foam::diffusivityModel::D(const label i)
+{
+    return D_[i];
+}
+
+
+const Foam::volScalarField& Foam::diffusivityModel::D(const label i) const
+{
+    return D_[i];
+}
+
+
 Foam::volScalarField& Foam::diffusivityModel::mu()
 {
     return mu_;