access to D fields and correction velocity calculation

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)

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();
 

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;

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_;