From 44b3dbfabd11de766ccb10ed8eecb3354a99a952 Mon Sep 17 00:00:00 2001
From: ignis <ignis@postech.ac.kr>
Date: Sun, 27 Nov 2016 13:08:38 +0900
Subject: [PATCH] plasmaReactingFoam subcycling

---
 .../combustion/plasmaReactingFoam/PhiEqn.H    |  2 -
 .../combustion/plasmaReactingFoam/UEqn.H      |  2 +
 .../combustion/plasmaReactingFoam/YEqn.H      | 89 ------------------
 .../plasmaReactingFoam/neControls.H           |  4 +
 .../combustion/plasmaReactingFoam/neEqn.H     | 92 +++++++++++++++++++
 .../plasmaReactingFoam/neEqnSubCycle.H        | 22 +++++
 .../plasmaReactingFoam/numberDensity.H        |  3 +
 .../plasmaReactingFoam/plasmaReactingFoam.C   | 11 +++
 8 files changed, 134 insertions(+), 91 deletions(-)
 create mode 100644 applications/solvers/combustion/plasmaReactingFoam/neControls.H
 create mode 100644 applications/solvers/combustion/plasmaReactingFoam/neEqn.H
 create mode 100644 applications/solvers/combustion/plasmaReactingFoam/neEqnSubCycle.H
 create mode 100644 applications/solvers/combustion/plasmaReactingFoam/numberDensity.H

diff --git a/applications/solvers/combustion/plasmaReactingFoam/PhiEqn.H b/applications/solvers/combustion/plasmaReactingFoam/PhiEqn.H
index dbb68e70..6a86a1ea 100644
--- a/applications/solvers/combustion/plasmaReactingFoam/PhiEqn.H
+++ b/applications/solvers/combustion/plasmaReactingFoam/PhiEqn.H
@@ -9,8 +9,6 @@
 
     E = -fvc::grad(Phi);
 
-    ng = p / T * (NA / R);
-
     tmp<volScalarField> tMagE (mag(E));
     const volScalarField &magE = tMagE();
 
diff --git a/applications/solvers/combustion/plasmaReactingFoam/UEqn.H b/applications/solvers/combustion/plasmaReactingFoam/UEqn.H
index aa678de6..900648e7 100644
--- a/applications/solvers/combustion/plasmaReactingFoam/UEqn.H
+++ b/applications/solvers/combustion/plasmaReactingFoam/UEqn.H
@@ -20,3 +20,5 @@
         fvOptions.correct(U);
         K = 0.5*magSqr(U);
     }
+
+    q = linearInterpolate(U) & mesh.Sf();
diff --git a/applications/solvers/combustion/plasmaReactingFoam/YEqn.H b/applications/solvers/combustion/plasmaReactingFoam/YEqn.H
index 639bb567..5ed001d5 100644
--- a/applications/solvers/combustion/plasmaReactingFoam/YEqn.H
+++ b/applications/solvers/combustion/plasmaReactingFoam/YEqn.H
@@ -10,48 +10,11 @@ tmp<fv::convectionScheme<scalar> > mvConvection
 );
 
 {
-    reaction->correct();
-    dQ = reaction->dQ();
     label inertIndex = -1;
     volScalarField Yt(0.0*Y[0]);
 
     composition.calculateDiffusivities(p, T);
 
-    EnTd = En.internalField();
-    EnTd *= EnToTableUnit;
-
-    Te.internalField() = TeOfEn.value(EnTd) * TeFac;
-    forAll(rho, celli)
-    {
-        Te[celli] = max(Te[celli], T[celli]);
-    }
-    Te.correctBoundaryConditions();
-
-    if (mobility_f_of_Te)
-    {
-        EnTd = Te.internalField();
-        EnTd *= TeToTableUnit;
-    }
-
-    mue.internalField() = mueN.value(EnTd) * mueNFac;
-
-    if (calculateDe)
-    {
-        De = mue * Te * (kB / eCharge);
-    }
-    else
-    {
-        De.internalField() = DeN.value(EnTd) * DeNFac;
-    }
-
-    mue.correctBoundaryConditions();
-    De.correctBoundaryConditions();
-
-    q = linearInterpolate(U) & mesh.Sf();
-
-
-
-
     const surfaceScalarField &msf = mesh.magSf();
     const surfaceVectorField &sf = mesh.Sf();
 
@@ -132,58 +95,6 @@ tmp<fv::convectionScheme<scalar> > mvConvection
 
         if (Y[i].name() == electronSpecie)
         {
-            Udrift = -(mue/ng)*E;
-            Uthermal = -((De/ng/Te)*fvc::grad(Te));
-            ve = (linearInterpolate(Udrift+Uthermal) & mesh.Sf()) + q;
-
-            // Wall electron flux correction
-            forAll (wallPatcheIDs, pidx)
-            {
-                label patchID = wallPatcheIDs[pidx];
-                // Probability of electron reflex
-                scalar pReflex = wallReflexes[pidx];
-                pReflex = max(min(pReflex,1.0),0.0);
-
-                fvsPatchScalarField &wallFlux = ve.boundaryField()[patchID];
-                const fvsPatchScalarField &wallMSf = msf.boundaryField()[patchID];
-                const fvPatchScalarField &wallTe = Te.boundaryField()[patchID];
-
-                scalarField vt(sqrt(8.0*kB.value()/pi/eMass.value()*wallTe) / 4.0);
-
-                // remove negative wallFlux value (flux from wall)
-                wallFlux = max(wallFlux, 0.0);
-
-                // add flux by thermal velocity
-                wallFlux += vt * wallMSf;
-                wallFlux *= (1.0-pReflex);
-            }
-
-            tmp<fvScalarMatrix> electronR(
-                    new fvScalarMatrix(ne,
-                        ne.dimensions()*dimVol/dimTime));
-            electronR->source() = reaction->R(Yi)->source();
-
-            fvScalarMatrix neEqn
-                (
-                 fvm::ddt(ne)
-                 + fvm::div(ve, ne)
-                 - fvm::laplacian(De/ng, ne)
-                 ==
-                 electronR
-                 + fvOptions(ne)
-                );
-
-            neEqn.relax();
-
-            fvOptions.constrain(neEqn);
-
-            neEqn.solve(mesh.solver("ne"));
-
-            fvOptions.correct(ne);
-
-            ne.writeMinMax(Info);
-
-            ne.max(0.0);
         }
         else if (Y[i].name() != inertSpecie)
         {
diff --git a/applications/solvers/combustion/plasmaReactingFoam/neControls.H b/applications/solvers/combustion/plasmaReactingFoam/neControls.H
new file mode 100644
index 00000000..dc1b0b96
--- /dev/null
+++ b/applications/solvers/combustion/plasmaReactingFoam/neControls.H
@@ -0,0 +1,4 @@
+const dictionary& neControls = mesh.solverDict(ne.name());
+
+label nNeSubCycles(readLabel(neControls.lookup("nNeSubCycles")));
+
diff --git a/applications/solvers/combustion/plasmaReactingFoam/neEqn.H b/applications/solvers/combustion/plasmaReactingFoam/neEqn.H
new file mode 100644
index 00000000..5c551986
--- /dev/null
+++ b/applications/solvers/combustion/plasmaReactingFoam/neEqn.H
@@ -0,0 +1,92 @@
+{
+    // Electron swarm parameter
+
+    EnTd = En.internalField();
+    EnTd *= EnToTableUnit;
+
+    Te.internalField() = TeOfEn.value(EnTd) * TeFac;
+    forAll(rho, celli)
+    {
+        Te[celli] = max(Te[celli], T[celli]);
+    }
+    Te.correctBoundaryConditions();
+
+    if (mobility_f_of_Te)
+    {
+        EnTd = Te.internalField();
+        EnTd *= TeToTableUnit;
+    }
+
+    mue.internalField() = mueN.value(EnTd) * mueNFac;
+
+    if (calculateDe)
+    {
+        De = mue * Te * (kB / eCharge);
+    }
+    else
+    {
+        De.internalField() = DeN.value(EnTd) * DeNFac;
+    }
+
+    mue.correctBoundaryConditions();
+    De.correctBoundaryConditions();
+
+
+    Udrift = -(mue/ng)*E;
+    Uthermal = -((De/ng/Te)*fvc::grad(Te));
+    ve = (linearInterpolate(Udrift+Uthermal) & mesh.Sf()) + q;
+
+    const surfaceScalarField &msf = mesh.magSf();
+
+    // Wall electron flux correction
+    forAll (wallPatcheIDs, pidx)
+    {
+        label patchID = wallPatcheIDs[pidx];
+        // Probability of electron reflex
+        scalar pReflex = wallReflexes[pidx];
+        pReflex = max(min(pReflex,1.0),0.0);
+
+        fvsPatchScalarField &wallFlux = ve.boundaryField()[patchID];
+        const fvsPatchScalarField &wallMSf = msf.boundaryField()[patchID];
+        const fvPatchScalarField &wallTe = Te.boundaryField()[patchID];
+
+        scalarField vt(sqrt(8.0*kB.value()/pi/eMass.value()*wallTe) / 4.0);
+
+        // remove negative wallFlux value (flux from wall)
+        wallFlux = max(wallFlux, 0.0);
+
+        // add flux by thermal velocity
+        wallFlux += vt * wallMSf;
+        wallFlux *= (1.0-pReflex);
+    }
+
+
+    volScalarField& Yi = composition.Y(electronSpecie);
+
+    tmp<fvScalarMatrix> electronR(
+            new fvScalarMatrix(ne, ne.dimensions()*dimVol/dimTime)
+            );
+    electronR->source() = reaction->R(Yi)->source();
+
+    fvScalarMatrix neEqn
+        (
+         fvm::ddt(ne)
+         + fvm::div(ve, ne)
+         - fvm::laplacian(De/ng, ne)
+         ==
+         electronR
+         + fvOptions(ne)
+        );
+
+    neEqn.relax();
+
+    fvOptions.constrain(neEqn);
+
+    neEqn.solve(mesh.solver("ne"));
+
+    fvOptions.correct(ne);
+
+    ne.writeMinMax(Info);
+
+    ne.max(0.0);
+}
diff --git a/applications/solvers/combustion/plasmaReactingFoam/neEqnSubCycle.H b/applications/solvers/combustion/plasmaReactingFoam/neEqnSubCycle.H
new file mode 100644
index 00000000..a3e7bba6
--- /dev/null
+++ b/applications/solvers/combustion/plasmaReactingFoam/neEqnSubCycle.H
@@ -0,0 +1,22 @@
+if (nNeSubCycles > 1)
+{
+    dimensionedScalar totalDeltaT = runTime.deltaT();
+
+    for
+    (
+        subCycle<volScalarField> neSubCycle(ne, nNeSubCycles);
+        !(++neSubCycle).end();
+    )
+    {
+        #include "neEqn.H"
+        #include "PhiEqn.H"
+        // rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi;
+    }
+
+    // rhoPhi = rhoPhiSum;
+}
+else
+{
+    #include "neEqn.H"
+    #include "PhiEqn.H"
+}
diff --git a/applications/solvers/combustion/plasmaReactingFoam/numberDensity.H b/applications/solvers/combustion/plasmaReactingFoam/numberDensity.H
new file mode 100644
index 00000000..c38cfb9a
--- /dev/null
+++ b/applications/solvers/combustion/plasmaReactingFoam/numberDensity.H
@@ -0,0 +1,3 @@
+{
+    ng = p / T * (NA / R);
+}
diff --git a/applications/solvers/combustion/plasmaReactingFoam/plasmaReactingFoam.C b/applications/solvers/combustion/plasmaReactingFoam/plasmaReactingFoam.C
index dac1623c..ed1062c6 100644
--- a/applications/solvers/combustion/plasmaReactingFoam/plasmaReactingFoam.C
+++ b/applications/solvers/combustion/plasmaReactingFoam/plasmaReactingFoam.C
@@ -30,6 +30,7 @@ Description
 \*---------------------------------------------------------------------------*/
 
 #include "fvCFD.H"
+#include "subCycle.H"
 #include "turbulenceModel.H"
 #include "psiCombustionModel.H"
 #include "multivariateScheme.H"
@@ -71,12 +72,22 @@ int main(int argc, char *argv[])
         runTime++;
         Info<< "Time = " << runTime.timeName() << nl << endl;
 
+        #include "numberDensity.H"
+
         #include "PhiEqn.H"
         #include "rhoEqn.H"
 
         while (pimple.loop())
         {
             #include "UEqn.H"
+
+            reaction->correct();
+            dQ = reaction->dQ();
+
+            #include "neControls.H"
+            #include "neEqnSubCycle.H"
+            // #include "neEqn.H"
+
             #include "YEqn.H"
             #include "EEqn.H"