eReactingFoam-2.4.x/neEqn.H
2017-05-17 20:44:45 +09:00

113 lines
3 KiB
C

{
// Electron swarm parameter
if (updateElectronTransport)
{
Te = T;
Te.correctBoundaryConditions();
if (!constETransport)
{
scalarField &W (csSpeciesW);
forAll (Te, cid)
{
scalar eEps = (3./2.) * kB.value() * Te[cid] / eCharge.value();
scalar sumX = 0.0;
sigma_m[cid] = 0;
forAll (csSpecies, isp)
{
scalar Xi = Y[csSpeciesI[isp]][cid]/W[isp];
sumX += Xi;
sigma_m[cid] += Xi * csList[isp].value(eEps);
}
sigma_m[cid] /= sumX;
}
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);
}