plasmaReactingFoam ion drift

heThermo init_qc boundary correction
 canteraTransport break after Xi=1 found
 BolsigTableReactionRate DataEntry flexiblity and unit conversion

applications/solvers/combustion/plasmaReactingFoam/YEqn.H

@@ -92,11 +92,23 @@ tmp<fv::convectionScheme<scalar> > mvConvection
         }
         else if (Y[i].name() != inertSpecie)
         {
+            const scalar z(composition.z(i));
+            const label nCharge(z);
+
+            if (nCharge != 0)
+            {
+                phi_drift = phi;
+                phi_drift += fvc::interpolate((rho*Di/T*(eCharge*z/kB))*E) & mesh.Sf();
+            }
 
             fvScalarMatrix YiEqn
             (
                 fvm::ddt(rho, Yi)
-              + mvConvection->fvmDiv(phi, Yi)
+              +
+              ( nCharge == 0
+              ? mvConvection->fvmDiv(phi, Yi)
+              : mvConvection->fvmDiv(phi_drift, Yi)
+              )
               // - fvm::laplacian(turbulence->muEff(), Yi)
               - fvm::laplacian(rho*Di, Yi)
              ==

applications/solvers/combustion/plasmaReactingFoam/createFields.H

@@ -354,3 +354,16 @@ surfaceScalarField ve
     Udrift & mesh.Sf()
 );
 
+
+surfaceScalarField phi_drift
+(
+    IOobject
+    (
+        "phi_drift",
+        runTime.timeName(),
+        mesh,
+        IOobject::NO_READ,
+        IOobject::NO_WRITE
+    ),
+    phi
+);

src/thermophysicalModels/basic/heThermo/heThermo.C

@@ -90,6 +90,15 @@ void Foam::heThermo<BasicThermo, MixtureType>::init_qc()
         qcCells[celli] =
             this->cellMixture(celli).QC();
     }
+
+    forAll(qc_.boundaryField(), patchi)
+    {
+        forAll(qc_.boundaryField()[patchi], facei)
+        {
+            qc_.boundaryField()[patchi][facei]
+                = this->patchFaceMixture(patchi, facei).QC();
+        }
+    }
 }
 
 

src/thermophysicalModels/reactionThermo/mixtures/SpecieMixture/SpecieMixture.C

@@ -65,6 +65,16 @@ Foam::scalar Foam::SpecieMixture<MixtureType>::W
 }
 
 
+template<class MixtureType>
+Foam::scalar Foam::SpecieMixture<MixtureType>::z
+(
+    const label speciei
+) const
+{
+    return this->getLocalThermo(speciei).z();
+}
+
+
 template<class MixtureType>
 Foam::scalar Foam::SpecieMixture<MixtureType>::Cp
 (

src/thermophysicalModels/reactionThermo/mixtures/SpecieMixture/SpecieMixture.H

@@ -80,6 +80,9 @@ public:
             //- Molecular weight [kg/kmol]
             virtual scalar W(const label speciei) const;
 
+            //- Number of elementary charges []
+            virtual scalar z(const label specieI) const;
+
 
         // Per specie thermo properties
 

src/thermophysicalModels/reactionThermo/mixtures/basicMultiComponentMixture/basicMultiComponentMixture.H

@@ -157,6 +157,9 @@ public:
             //- Molecular weight [kg/kmol]
             virtual scalar W(const label specieI) const = 0;
 
+            //- Number of elementary charges []
+            virtual scalar z(const label specieI) const = 0;
+
 
         // Per specie thermo properties
 

src/thermophysicalModels/specie/reaction/reactionRate/BolsigTableReactionRate/BolsigTableReactionRate.H

@@ -63,7 +63,16 @@ class BolsigTableReactionRate
         // scalar beta_;
         // scalar Ta_;
 
-        autoPtr<TableBase<scalar> > ke_;
+        //- Unit conversion factor to m3/s
+        const scalar cvFac_;
+
+        //- True if rate coef is in cm3/s or m3/s units
+        const Switch unitInNumber_;
+
+        //- Avogadro's number if unitInNumber_ or one
+        const scalar numberToMole_;
+
+        autoPtr<DataEntry<scalar> > ke_;
 
 
 public:

src/thermophysicalModels/specie/reaction/reactionRate/BolsigTableReactionRate/BolsigTableReactionRateI.H

@@ -29,6 +29,9 @@ inline Foam::BolsigTableReactionRate::BolsigTableReactionRate
 (
 )
 :
+    cvFac_(1.0),
+    numberToMole_(1.0),
+    unitInNumber_(true),
     ke_(NULL)
 {}
 
@@ -39,6 +42,9 @@ inline Foam::BolsigTableReactionRate::BolsigTableReactionRate
     Istream& is
 )
 :
+    cvFac_(1.0),
+    numberToMole_(1.0),
+    unitInNumber_(true),
     ke_(NULL)
 {
 }
@@ -50,7 +56,10 @@ inline Foam::BolsigTableReactionRate::BolsigTableReactionRate
     const dictionary& dict
 )
 :
-    ke_(new CSV<scalar> ("k", dict))
+    cvFac_(dict.lookupOrDefault("convertToSI", 1.0)),
+    unitInNumber_(dict.lookupOrDefault("unitInNumber", true)),
+    numberToMole_(unitInNumber_ ? 6.0221417930e+23 : 1.0),
+    ke_(DataEntry<scalar>::New("k", dict))
 {}
 
 
@@ -63,12 +72,12 @@ inline Foam::scalar Foam::BolsigTableReactionRate::operator()
     const scalarField&
 ) const
 {
-    const scalar NA = 6.0221417930e+23; // mol^-1
+    const scalar NA = numberToMole_; // mol^-1
     const scalar mol2kmol = 1000.0; //mol/kmol^-1
     const scalar oneTd = 1.0e-21; // 1e-21 Vm^2 = 1 Td
 
     // p is actually E/n in Vm^2 unit
-    scalar coef =  ke_->value(p / oneTd) * NA * mol2kmol;
+    scalar coef =  ke_->value(p / oneTd) * NA * mol2kmol * cvFac_;
     // Info << ke_->value(p * 1e21) << endl;
     return coef;
     // return ke_->value(p / 1e-21) * NA * mol2kmol;

src/thermophysicalModels/specie/transport/cantera/canteraTransport.C

@@ -164,7 +164,11 @@ void Foam::canteraTransport<Thermo>::calculateDiffusivities
     scalar sumX = sum(X_);
     forAll (X_, i)
     {
-        if (X_[i]/sumX == 1.0) idxFracOne = i;
+        if (X_[i]/sumX == 1.0)
+        {
+            idxFracOne = i;
+            break;
+        }
     }
 
     if (idxFracOne < 0)