calculation of correction flux

YEqn.H

@@ -18,11 +18,18 @@ tmp<fv::convectionScheme<scalar>> mvConvection
     // volVectorField Vc(U);
     Vc = dimensionedVector("zero", dimLength/dimTime, U[0]*0.0);
 
+    phiC = dimensionedScalar("zero", dimDensity*dimVelocity*dimArea, 0.0);
+
     forAll(Y, i)
     {
         Vc += diff.D(i) * fvc::grad(Y[i]);
     }
 
+    forAll(Y, i)
+    {
+        phiC += linearInterpolate(rho * diff.D(i)) * fvc::snGrad(Y[i]) * mesh.magSf();
+    }
+
     forAll(Y, i)
     {
         if (Y[i].name() != inertSpecie)

createFields.H

@@ -141,3 +141,17 @@ volScalarField dQ
 #include "createMRF.H"
 
 diff.correct();
+
+
+surfaceScalarField phiC
+(
+    IOobject
+    (
+        "phiC",
+        runTime.timeName(),
+        mesh,
+        IOobject::READ_IF_PRESENT,
+        IOobject::AUTO_WRITE
+    ),
+    linearInterpolate(rho*U) & mesh.Sf()
+);