eReactingFoam-2.4.x/neEqn.H

112 lines
2.7 KiB
C

{
// Electron swarm parameter
if (updateElectronTransport)
{
Te = T;
Te.correctBoundaryConditions();
scalarField &W (csSpeciesW);
Xsum = 0.0;
sigma_m = 0.0;
meanEps = (3./2.) * kB.value() / eCharge.value() * Te.internalField();
forAll (csSpecies, isp)
{
label k = csSpeciesI[isp];
scalarField X (Y[k].internalField() / W[isp]);
Xsum += X;
sigma_m += X * csList[isp].value(meanEps);
}
sigma_m /= Xsum;
Info << min(sigma_m) << " / " << max(sigma_m) << endl;
/*
mue.internalField() = eCharge.value()
/ (sqrt(eMass.value() * kB.value() * Te) * sigma_m);
mue.correctBoundaryConditions();
*/
mue.internalField() = mueNFac * eCharge.value() / 3.0
/ (sqrt(pi * eMass.value() * kB.value() * Te / 8.0) * sigma_m);
mue.correctBoundaryConditions();
De = mue * Te * (kB / eCharge);
De.correctBoundaryConditions();
updateElectronTransport = false;
mue.writeMinMax(Info);
De.writeMinMax(Info);
}
Udrift = -(mue/ng)*E;
if (relaxDrift < 1.0 && relaxDrift > 0.0)
{
Udrift *= relaxDrift;
}
// Uthermal = -((De/ng/Te)*fvc::grad(Te));
Uthermal = -(fvc::grad(De/ng));
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);
}