From 41accee91d415c5a011a33d98d15fa7be9237047 Mon Sep 17 00:00:00 2001
From: Yeongdo Park <ignis@postech.ac.kr>
Date: Sun, 17 Feb 2019 09:13:51 -0500
Subject: [PATCH] energy diffusion using h gradient, Lorentz force to momentum,
 drift to Y equation

---
 EEqn.H                |  2 +-
 UEqn.H                |  3 ++-
 YEqn.H                | 22 +++++++++-------------
 createElectricField.H | 15 ---------------
 createFields.H        |  8 --------
 5 files changed, 12 insertions(+), 38 deletions(-)

diff --git a/EEqn.H b/EEqn.H
index 70a923c..a0d4151 100644
--- a/EEqn.H
+++ b/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()
diff --git a/UEqn.H b/UEqn.H
index 1058faa..d965fd5 100644
--- a/UEqn.H
+++ b/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();
 
diff --git a/YEqn.H b/YEqn.H
index 98fcefb..cbbcf98 100644
--- a/YEqn.H
+++ b/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)
diff --git a/createElectricField.H b/createElectricField.H
index 02bd427..033202d 100644
--- a/createElectricField.H
+++ b/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
 (
diff --git a/createFields.H b/createFields.H
index b69d6a9..1044431 100644
--- a/createFields.H
+++ b/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