plasmaReactingFoam electro static momentum source

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);
 }

applications/solvers/combustion/plasmaReactingFoam/UEqn.H

@@ -5,6 +5,8 @@
       + turbulence->divDevRhoReff(U)
      ==
         rho*g
+      + rho*qc*E
+      + eCharge*ne*E
       + fvOptions(rho, U)
     );
 

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
+);