From 026a8eb2be7530cfbf09215b29d4781074d591df Mon Sep 17 00:00:00 2001
From: ignis <ignis@postech.ac.kr>
Date: Sun, 9 Oct 2016 20:36:33 +0900
Subject: [PATCH] plasmaReactingFoam ne equation flux mean flow component and
 electron aborbing wall BC

---
 .../combustion/plasmaReactingFoam/YEqn.H      | 28 ++++++++++-
 .../plasmaReactingFoam/createFields.H         | 47 ++++++++++++++++++-
 2 files changed, 73 insertions(+), 2 deletions(-)

diff --git a/applications/solvers/combustion/plasmaReactingFoam/YEqn.H b/applications/solvers/combustion/plasmaReactingFoam/YEqn.H
index d870f6bd..69fb948a 100644
--- a/applications/solvers/combustion/plasmaReactingFoam/YEqn.H
+++ b/applications/solvers/combustion/plasmaReactingFoam/YEqn.H
@@ -27,13 +27,39 @@ tmp<fv::convectionScheme<scalar> > mvConvection
     mue.correctBoundaryConditions();
     Te.correctBoundaryConditions();
 
+
+    q = linearInterpolate(U) & mesh.Sf();
+
+
     forAll(Y, i)
     {
         if (Y[i].name() == electronSpecie)
         {
             volScalarField& Yi = Y[i];
 
-            ve = - linearInterpolate(mue*E/ng) & mesh.Sf();
+            Udrift = - linearInterpolate(mue*E/ng);
+            ve = (Udrift & mesh.Sf()) + q;
+
+            const surfaceScalarField &msf = mesh.magSf();
+
+            // Wall electron flux correction
+            forAll (wallPatcheIDs, pidx)
+            {
+                label patchID = wallPatcheIDs[pidx];
+                fvsPatchScalarField &wallFlux = ve.boundaryField()[patchID];
+
+                const fvsPatchScalarField &wallMSf = msf.boundaryField()[patchID];
+
+                const fvPatchScalarField &wallTe = Te.boundaryField()[patchID];
+
+                scalarField vt(sqrt((8.0/pi)*(2.0/3.0)/16.0*eCharge.value()/eMass.value()*wallTe));
+
+                // remove negative wallFlux value (flux from wall)
+                wallFlux = max(wallFlux, 0.0);
+
+                // add flux by thermal velocity
+                wallFlux += vt * wallMSf;
+            }
 
             tmp<fvScalarMatrix> electronR(
                     new fvScalarMatrix(ne, dimless/dimTime));
diff --git a/applications/solvers/combustion/plasmaReactingFoam/createFields.H b/applications/solvers/combustion/plasmaReactingFoam/createFields.H
index 994759cc..b3f3d315 100644
--- a/applications/solvers/combustion/plasmaReactingFoam/createFields.H
+++ b/applications/solvers/combustion/plasmaReactingFoam/createFields.H
@@ -41,6 +41,36 @@ scalar TeFac (
         physicalProperties.lookupOrDefault("TeFac", 1.0)
         );
 
+dictionary wallElectronFlux
+(
+    physicalProperties.subDict("wallElectronFlux")
+);
+
+word TeName(wallElectronFlux.lookup("TeName"));
+
+wordList wallPatcheNames (wallElectronFlux.lookup("wallPatches"));
+labelList wallPatcheIDs (wallPatcheNames.size(), 0);
+
+forAll (wallPatcheNames, pi)
+{
+    word patchName = wallPatcheNames[pi];
+    label patchID = mesh.boundaryMesh().findPatchID(patchName);
+    wallPatcheIDs[pi] = patchID;
+
+/*
+    Info<< patchName << patchID << endl;
+    std::cout << Pstream::myProcNo() << patchName << patchID << std::endl;
+
+    OStringStream temp_ss;
+    temp_ss << Pstream::myProcNo() << mesh.boundaryMesh();
+    std::cout << temp_ss.str() << endl;
+*/
+}
+
+Info<< TeName << endl;
+Info<< wallPatcheNames << endl;
+Info<< wallPatcheIDs << endl;
+
 
 Info<< "Reading field Phi\n" << endl;
 volScalarField Phi
@@ -126,6 +156,19 @@ const volScalarField& T = thermo.T();
 
 #include "compressibleCreatePhi.H"
 
+surfaceScalarField q
+(
+    IOobject
+    (
+        "q",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    linearInterpolate(U) & mesh.Sf()
+);
+
 Info << "Creating turbulence model.\n" << nl;
 autoPtr<compressible::turbulenceModel> turbulence
 (
@@ -263,6 +306,8 @@ bolos.presolve();
 */
 
 Info<< "Calculating face flux field ve\n" << endl;
+surfaceVectorField Udrift ("Udrift", - linearInterpolate(mue*E/ng));
+
 surfaceScalarField ve
 (
     IOobject
@@ -273,6 +318,6 @@ surfaceScalarField ve
         IOobject::NO_READ,
         IOobject::NO_WRITE
     ),
-    - linearInterpolate(mue*E/ng) & mesh.Sf()
+    Udrift & mesh.Sf()
 );