diff --git a/applications/solvers/combustion/plasmaReactingFoam/YEqn.H b/applications/solvers/combustion/plasmaReactingFoam/YEqn.H index 50850456..5caaedf2 100644 --- a/applications/solvers/combustion/plasmaReactingFoam/YEqn.H +++ b/applications/solvers/combustion/plasmaReactingFoam/YEqn.H @@ -92,11 +92,23 @@ tmp > 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) == diff --git a/applications/solvers/combustion/plasmaReactingFoam/createFields.H b/applications/solvers/combustion/plasmaReactingFoam/createFields.H index c0192ac3..8eb34ed6 100644 --- a/applications/solvers/combustion/plasmaReactingFoam/createFields.H +++ b/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 +); diff --git a/src/thermophysicalModels/basic/heThermo/heThermo.C b/src/thermophysicalModels/basic/heThermo/heThermo.C index d6c99c81..2cca5709 100644 --- a/src/thermophysicalModels/basic/heThermo/heThermo.C +++ b/src/thermophysicalModels/basic/heThermo/heThermo.C @@ -90,6 +90,15 @@ void Foam::heThermo::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(); + } + } } diff --git a/src/thermophysicalModels/reactionThermo/mixtures/SpecieMixture/SpecieMixture.C b/src/thermophysicalModels/reactionThermo/mixtures/SpecieMixture/SpecieMixture.C index 938244b5..552f9fe5 100644 --- a/src/thermophysicalModels/reactionThermo/mixtures/SpecieMixture/SpecieMixture.C +++ b/src/thermophysicalModels/reactionThermo/mixtures/SpecieMixture/SpecieMixture.C @@ -65,6 +65,16 @@ Foam::scalar Foam::SpecieMixture::W } +template +Foam::scalar Foam::SpecieMixture::z +( + const label speciei +) const +{ + return this->getLocalThermo(speciei).z(); +} + + template Foam::scalar Foam::SpecieMixture::Cp ( diff --git a/src/thermophysicalModels/reactionThermo/mixtures/SpecieMixture/SpecieMixture.H b/src/thermophysicalModels/reactionThermo/mixtures/SpecieMixture/SpecieMixture.H index a954119d..a7a9329b 100644 --- a/src/thermophysicalModels/reactionThermo/mixtures/SpecieMixture/SpecieMixture.H +++ b/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 diff --git a/src/thermophysicalModels/reactionThermo/mixtures/basicMultiComponentMixture/basicMultiComponentMixture.H b/src/thermophysicalModels/reactionThermo/mixtures/basicMultiComponentMixture/basicMultiComponentMixture.H index fdfae416..970dad5d 100644 --- a/src/thermophysicalModels/reactionThermo/mixtures/basicMultiComponentMixture/basicMultiComponentMixture.H +++ b/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 diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/BolsigTableReactionRate/BolsigTableReactionRate.H b/src/thermophysicalModels/specie/reaction/reactionRate/BolsigTableReactionRate/BolsigTableReactionRate.H index 05b6297d..9805541f 100644 --- a/src/thermophysicalModels/specie/reaction/reactionRate/BolsigTableReactionRate/BolsigTableReactionRate.H +++ b/src/thermophysicalModels/specie/reaction/reactionRate/BolsigTableReactionRate/BolsigTableReactionRate.H @@ -63,7 +63,16 @@ class BolsigTableReactionRate // scalar beta_; // scalar Ta_; - autoPtr > 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 > ke_; public: diff --git a/src/thermophysicalModels/specie/reaction/reactionRate/BolsigTableReactionRate/BolsigTableReactionRateI.H b/src/thermophysicalModels/specie/reaction/reactionRate/BolsigTableReactionRate/BolsigTableReactionRateI.H index a1d4e6b8..1de353f0 100644 --- a/src/thermophysicalModels/specie/reaction/reactionRate/BolsigTableReactionRate/BolsigTableReactionRateI.H +++ b/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 ("k", dict)) + cvFac_(dict.lookupOrDefault("convertToSI", 1.0)), + unitInNumber_(dict.lookupOrDefault("unitInNumber", true)), + numberToMole_(unitInNumber_ ? 6.0221417930e+23 : 1.0), + ke_(DataEntry::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; diff --git a/src/thermophysicalModels/specie/transport/cantera/canteraTransport.C b/src/thermophysicalModels/specie/transport/cantera/canteraTransport.C index a90d6584..8747bc0a 100644 --- a/src/thermophysicalModels/specie/transport/cantera/canteraTransport.C +++ b/src/thermophysicalModels/specie/transport/cantera/canteraTransport.C @@ -164,7 +164,11 @@ void Foam::canteraTransport::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)