251 lines
7 KiB
C
251 lines
7 KiB
C
tmp<fv::convectionScheme<scalar> > mvConvection
|
|
(
|
|
fv::convectionScheme<scalar>::New
|
|
(
|
|
mesh,
|
|
fields,
|
|
phi,
|
|
mesh.divScheme("div(phi,Yi_h)")
|
|
)
|
|
);
|
|
|
|
{
|
|
reaction->correct();
|
|
dQ = reaction->dQ();
|
|
label inertIndex = -1;
|
|
volScalarField Yt(0.0*Y[0]);
|
|
|
|
composition.calculateDiffusivities(p, T);
|
|
|
|
EnTd = En.internalField();
|
|
EnTd *= EnToTableUnit;
|
|
|
|
Te.internalField() = TeOfEn.value(EnTd) * TeFac;
|
|
forAll(rho, celli)
|
|
{
|
|
Te[celli] = max(Te[celli], T[celli]);
|
|
}
|
|
Te.correctBoundaryConditions();
|
|
|
|
if (mobility_f_of_Te)
|
|
{
|
|
EnTd = Te.internalField();
|
|
EnTd *= TeToTableUnit;
|
|
}
|
|
|
|
mue.internalField() = mueN.value(EnTd) * mueNFac;
|
|
|
|
if (calculateDe)
|
|
{
|
|
De = mue * Te * (kB / eCharge);
|
|
}
|
|
else
|
|
{
|
|
De.internalField() = DeN.value(EnTd) * DeNFac;
|
|
}
|
|
|
|
mue.correctBoundaryConditions();
|
|
De.correctBoundaryConditions();
|
|
|
|
q = linearInterpolate(U) & mesh.Sf();
|
|
|
|
|
|
|
|
|
|
const surfaceScalarField &msf = mesh.magSf();
|
|
const surfaceVectorField &sf = mesh.Sf();
|
|
|
|
forAll(ions, k) // ion-neutral pair
|
|
{
|
|
const word nIon(ions[k]);
|
|
const word nNeu(neutrals[k]);
|
|
|
|
const volScalarField& Di = composition.D(nIon);
|
|
const scalar z(composition.z(composition.species()[nIon]));
|
|
|
|
// P_Reflex list for the ion
|
|
const scalarList &rK = reflexes[k];
|
|
|
|
surfaceScalarField::GeometricBoundaryField &bfIonFlux
|
|
= ionFluxBFs[k];
|
|
surfaceScalarField::GeometricBoundaryField &bfNeuFlux
|
|
= neutralFluxBFs[k];
|
|
|
|
bfIonFlux = phi.boundaryField();
|
|
bfNeuFlux = phi.boundaryField();
|
|
|
|
// Adding drift flux to boundary patches
|
|
forAll (bfIonFlux, pidx)
|
|
{
|
|
bfIonFlux[pidx] +=
|
|
(E.boundaryField()[pidx]
|
|
& sf.boundaryField()[pidx])
|
|
* rho.boundaryField()[pidx]
|
|
* Di.boundaryField()[pidx]
|
|
/ T.boundaryField()[pidx]
|
|
* (eCharge*z/kB).value();
|
|
}
|
|
|
|
const scalar WIon(composition.W(composition.species()[nIon]));
|
|
const scalar WNeu(composition.W(composition.species()[nNeu]));
|
|
|
|
const scalar MIon(WIon / NA.value() / 1000.0);
|
|
const scalar MNeu(WNeu / NA.value() / 1000.0);
|
|
|
|
const volScalarField& Yion = composition.Y(nIon);
|
|
const volScalarField& Yneu = composition.Y(nNeu);
|
|
|
|
forAll(wallPatcheIDs, pidx) // loop over wall patches
|
|
{
|
|
label patchID = wallPatcheIDs[pidx];
|
|
|
|
// Probability of ion reflex
|
|
const scalar pReflex = max(min(rK[pidx],1.0),0.0);
|
|
|
|
scalarField &wallFluxIon = bfIonFlux[patchID];
|
|
scalarField &wallFluxNeu = bfNeuFlux[patchID];
|
|
|
|
const scalarField &wallMSf = msf.boundaryField()[patchID];
|
|
const scalarField &wallT = T.boundaryField()[patchID];
|
|
const scalarField &wallYion = Yion.boundaryField()[patchID];
|
|
const scalarField &wallYneu = Yneu.boundaryField()[patchID];
|
|
|
|
scalarField vt(sqrt(8.0*kB.value()/pi/MIon*wallT) / 4.0);
|
|
|
|
// remove negative wallFlux value (flux from wall)
|
|
wallFluxIon = max(wallFluxIon, 0.0);
|
|
|
|
// add flux by thermal velocity
|
|
wallFluxIon += vt * wallMSf;
|
|
wallFluxIon *= (1.0 - pReflex);
|
|
|
|
// add flux by ion neutralization
|
|
wallFluxNeu -= wallFluxIon * wallYion / wallYneu / (WIon / WNeu);
|
|
}
|
|
}
|
|
|
|
forAll(Y, i)
|
|
{
|
|
volScalarField& Yi = Y[i];
|
|
|
|
const volScalarField& Di = D[i];
|
|
|
|
if (Y[i].name() == electronSpecie)
|
|
{
|
|
Udrift = - linearInterpolate
|
|
((mue/ng)*E + ((De/ng/Te)*fvc::grad(Te)));
|
|
ve = (Udrift & mesh.Sf()) + q;
|
|
|
|
// 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);
|
|
}
|
|
|
|
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);
|
|
}
|
|
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();
|
|
}
|
|
|
|
if (ions.contains(Y[i].name()))
|
|
{
|
|
const label ibc = ions[Y[i].name()];
|
|
// phi_drift updated
|
|
phi_drift.boundaryField() = ionFluxBFs[ibc];
|
|
}
|
|
else if (neutrals.contains(Y[i].name()))
|
|
{
|
|
const label ibc = neutrals[Y[i].name()];
|
|
// update phi_neutral
|
|
phi_neutral.internalField() = phi.internalField();
|
|
phi_neutral.boundaryField() = neutralFluxBFs[ibc];
|
|
}
|
|
|
|
|
|
fvScalarMatrix YiEqn
|
|
(
|
|
fvm::ddt(rho, Yi)
|
|
+
|
|
( nCharge != 0
|
|
? mvConvection->fvmDiv(phi_drift, Yi)
|
|
: ( neutrals.contains(Y[i].name())
|
|
? mvConvection->fvmDiv(phi_neutral, Yi)
|
|
: mvConvection->fvmDiv(phi, Yi)
|
|
)
|
|
)
|
|
// - fvm::laplacian(turbulence->muEff(), Yi)
|
|
- fvm::laplacian(rho*Di, Yi)
|
|
==
|
|
reaction->R(Yi)
|
|
+ fvOptions(rho, Yi)
|
|
);
|
|
|
|
YiEqn.relax();
|
|
|
|
fvOptions.constrain(YiEqn);
|
|
|
|
YiEqn.solve(mesh.solver("Yi"));
|
|
|
|
fvOptions.correct(Yi);
|
|
|
|
Yi.max(0.0);
|
|
Yt += Yi;
|
|
}
|
|
else
|
|
{
|
|
inertIndex = i;
|
|
}
|
|
}
|
|
|
|
Y[inertIndex] = scalar(1) - Yt;
|
|
Y[inertIndex].max(0.0);
|
|
}
|