energy diffusion using h gradient, Lorentz force to momentum, drift to Y equation

EEqn.H

@@ -62,7 +62,7 @@
         fvm::ddt(rho, he) + hConvection->fvmDiv(phi, he)
       + fvc::ddt(rho, K) + fvc::div(phi, K)
       - dpdt
-      - fvc::laplacian(diff.k(), T)
+      - fvm::laplacian(diff.k()/thermo.Cp(), he)
       + fvc::div(hDiffusionSrc)
      ==
         reaction->Sh()

UEqn.H

@@ -13,7 +13,8 @@ tmp<fvVectorMatrix> tUEqn
      - fvm::laplacian(diff.mu(), U)
     )
  ==
-    fvOptions(rho, U)
+    rho * thermo.composition().Qc() * E
+  + fvOptions(rho, U)
 );
 fvVectorMatrix& UEqn = tUEqn.ref();
 

YEqn.H

@@ -1,28 +1,19 @@
 Vc = dimensionedVector("zero", dimLength/dimTime, vector(Zero));
 
-phiC = dimensionedScalar("zero", dimDensity*dimVelocity*dimArea, 0.0);
+phiC = phi;
 
 forAll(Y, i)
 {
     Vc += diff.D(i) * fvc::grad(Y[i]);
+    // Vc -= mu_i * E;
 }
 
 forAll(Y, i)
 {
     phiC += linearInterpolate(rho * diff.D(i)) * fvc::snGrad(Y[i]) * mesh.magSf();
+    // phiC -= linearInterpolate(rho * mu_i) * snE * mesh.magSf();
 }
 
-tmp<fv::convectionScheme<scalar>> mvConvection
-(
-    fv::convectionScheme<scalar>::New
-    (
-        mesh,
-        fields,
-        phi + phiC,
-        mesh.divScheme("div(phi,Yi_h)")
-    )
-);
-
 {
     reaction->correct();
     dQ = reaction->dQ();
@@ -38,7 +29,12 @@ tmp<fv::convectionScheme<scalar>> mvConvection
             fvScalarMatrix YiEqn
             (
                 fvm::ddt(rho, Yi)
-              + mvConvection->fvmDiv(phi + phiC, Yi)
+              + (
+                  thermo.composition().z(i) != 0
+                    ? fvm::div(phiC + (linearInterpolate(rho*diff.mu(i,T)*E) & mesh.Sf()),
+			             Yi, "div(phi,Yi_h)")
+                    : fvm::div(phiC, Yi, "div(phi,Yi_h)")
+                )
               - fvm::laplacian(rho * diff.D(i), Yi)
              ==
                 reaction->R(Yi)

createElectricField.H

@@ -40,21 +40,6 @@ volScalarField Phi
 );
 
 
-Info<< "Reading field rhoc\n" << endl;
-volScalarField rhoc
-(
-    IOobject
-    (
-        "rhoc",
-        runTime.timeName(),
-        mesh,
-        IOobject::MUST_READ,
-        IOobject::AUTO_WRITE
-    ),
-    mesh
-);
-
-
 Info<< "Creating field electric field\n" << endl;
 volVectorField E
 (

createFields.H

@@ -116,14 +116,6 @@ Info<< "Creating field kinetic energy K\n" << endl;
 volScalarField K("K", 0.5*magSqr(U));
 
 
-multivariateSurfaceInterpolationScheme<scalar>::fieldTable fields;
-
-forAll(Y, i)
-{
-    fields.add(Y[i]);
-}
-fields.add(thermo.he());
-
 volScalarField dQ
 (
     IOobject