From 747b4a0f4643d44dd8890216eb7d71240a9eb54c Mon Sep 17 00:00:00 2001
From: ignis <ignis@postech.ac.kr>
Date: Fri, 8 Apr 2016 14:37:30 +0900
Subject: [PATCH] plasmaReactingFoam electro static momentum source

---
 .../combustion/plasmaReactingFoam/PhiEqn.H    |  8 ++++---
 .../combustion/plasmaReactingFoam/UEqn.H      |  2 ++
 .../plasmaReactingFoam/createFields.H         | 22 +++++++++++++++++++
 3 files changed, 29 insertions(+), 3 deletions(-)

diff --git a/applications/solvers/combustion/plasmaReactingFoam/PhiEqn.H b/applications/solvers/combustion/plasmaReactingFoam/PhiEqn.H
index 71e4e1df..b73c80e9 100644
--- a/applications/solvers/combustion/plasmaReactingFoam/PhiEqn.H
+++ b/applications/solvers/combustion/plasmaReactingFoam/PhiEqn.H
@@ -1,8 +1,10 @@
 {
-    volScalarField& qc = thermo.qc();
-
     solve
     (
-        fvm::laplacian(Phi) + fvc::Sp(rho, qc)/epsilon0
+        fvm::laplacian(Phi)
+        + fvc::Sp(rho, qc)/epsilon0
+        + fvc::Sp(eCharge, ne)/epsilon0
     );
+
+    E = -fvc::grad(Phi);
 }
diff --git a/applications/solvers/combustion/plasmaReactingFoam/UEqn.H b/applications/solvers/combustion/plasmaReactingFoam/UEqn.H
index 60c2a8a7..36eb9e6a 100644
--- a/applications/solvers/combustion/plasmaReactingFoam/UEqn.H
+++ b/applications/solvers/combustion/plasmaReactingFoam/UEqn.H
@@ -5,6 +5,8 @@
       + turbulence->divDevRhoReff(U)
      ==
         rho*g
+      + rho*qc*E
+      + eCharge*ne*E
       + fvOptions(rho, U)
     );
 
diff --git a/applications/solvers/combustion/plasmaReactingFoam/createFields.H b/applications/solvers/combustion/plasmaReactingFoam/createFields.H
index 899a21bd..ea969cd1 100644
--- a/applications/solvers/combustion/plasmaReactingFoam/createFields.H
+++ b/applications/solvers/combustion/plasmaReactingFoam/createFields.H
@@ -37,6 +37,9 @@ volScalarField Phi
     mesh
 );
 
+Info<< "Creating field electric field\n" << endl;
+volVectorField E("E", -fvc::grad(Phi));
+
 Info<< "Creating reaction model\n" << endl;
 
 autoPtr<combustionModels::psiCombustionModel> reaction
@@ -47,7 +50,12 @@ autoPtr<combustionModels::psiCombustionModel> reaction
 psiReactionThermo& thermo = reaction->thermo();
 thermo.validate(args.executable(), "h", "e");
 
+volScalarField& qc = thermo.qc();
 basicMultiComponentMixture& composition = thermo.composition();
+
+//- Elementary charge (default in [C])
+const dimensionedScalar& eCharge = constant::electromagnetic::e;
+
 PtrList<volScalarField>& Y = composition.Y();
 
 word inertSpecie(thermo.lookup("inertSpecie"));
@@ -140,3 +148,17 @@ volScalarField dQ
     mesh,
     dimensionedScalar("dQ", dimEnergy/dimTime, 0.0)
 );
+
+Info<< "Creating field electron number density\n" << endl;
+volScalarField ne
+(
+    IOobject
+    (
+        "ne",
+        runTime.timeName(),
+        mesh,
+        IOobject::MUST_READ,
+        IOobject::AUTO_WRITE
+    ),
+    mesh
+);