eReactingFoam-2.4.x/YEqn.H

/*
tmp<fv::convectionScheme<scalar> > mvConvection
(
    fv::convectionScheme<scalar>::New
    (
        mesh,
        fields,
        phi,
        mesh.divScheme("div(phi,Yi_h)")
    )
);
*/

{
    label inertIndex = -1;
    volScalarField Yt(0.0*Y[0]);

    composition.calculateDiffusivities(p, T);

    const surfaceScalarField &msf = mesh.magSf();
    const surfaceVectorField &sf = mesh.Sf();

    forAll(ions, k) // ion-neutral pair
    {
        const word nIon(ions[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];
        bfIonFlux = phi.boundaryField();

        // Adding drift flux to boundary patches
        forAll (bfIonFlux, pidx)
        {
            Info << "Adding drift flux to boundary patches" << pidx << endl;

            bfIonFlux[pidx] +=
                snE.boundaryField()[pidx]
                * msf.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 MIon(WIon / NA.value() / 1000.0);

        const volScalarField& Yion = composition.Y(nIon);

        const hashedWordList &targets = targetList[k];

        forAll (targets, tidx)
        {
            const word &nNeu(targets[tidx]);

            surfaceScalarField::GeometricBoundaryField &
                bfNeuFlux = neutralFluxBFs[neutrals[nNeu]];
            bfNeuFlux = phi.boundaryField();

            const scalar WNeu(composition.W(composition.species()[nNeu]));
            const scalar MNeu(WNeu / NA.value() / 1000.0);

            const volScalarField& Yneu = composition.Y(nNeu);

            forAll(wallPatcheIDs, pidx) // loop over wall patches
            {
                const label patchID = wallPatcheIDs[pidx];

                // Probability of ion reflex
                const scalar pReflex = max(min(rK[pidx],1.0),0.0);

                scalarField &wallFluxIon = bfIonFlux[patchID];

                const scalarField &wallMSf  =  msf.boundaryField()[patchID];
                const scalarField &wallT    =    T.boundaryField()[patchID];
                const scalarField &wallYion = Yion.boundaryField()[patchID];

                const 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
                scalarField &wallFluxNeu = bfNeuFlux[patchID];

                const scalarField &wallYneu = Yneu.boundaryField()[patchID];

                wallFluxNeu -= wallFluxIon * wallYion / wallYneu
                                              / (WIon / WNeu);
            }
        }
    }

    forAll(Y, i)
    {
        volScalarField& Yi = Y[i];

        const volScalarField& Di = D[i];

        if (Y[i].name() == electronSpecie)
        {
        }
        else if (Y[i].name() != inertSpecie)
        {
            const scalar z(composition.z(i));
            const label nCharge(z);

            if (nCharge != 0)
            {
                phis[i] = phi;
                phis[i] += fvc::interpolate((rho*Di/T*(eCharge*z/kB))*E) & mesh.Sf();

                if (relaxDrift < 1.0 && relaxDrift > 0.0)
                {
                    phis[i] *= relaxDrift;
                }
            }

            if (ions.contains(Y[i].name()))
            {
                const label ibc = ions[Y[i].name()];
                // phis[i] updated
                phis[i].boundaryField() = ionFluxBFs[ibc];
            }
            else if (neutrals.contains(Y[i].name()))
            {
                const label ibc = neutrals[Y[i].name()];
                // update phis[i]
                phis[i].internalField() = phi.internalField();
                phis[i].boundaryField() = neutralFluxBFs[ibc];
            }


            fvScalarMatrix YiEqn
            (
                fvm::ddt(rho, Yi)
              +
              ( nCharge != 0
              ? fvm::div(phis[i], Yi, "div(phi,Yi_h)")
              : ( neutrals.contains(Y[i].name())
                  ? fvm::div(phis[i], Yi, "div(phi,Yi_h)")
                  : fvm::div(phi, Yi, "div(phi,Yi_h)")
                )
              )
              // - 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);

    updateElectronTransport = true;
}