diff --git a/applications/solvers/combustion/plasmaReactingFoam/YEqn.H b/applications/solvers/combustion/plasmaReactingFoam/YEqn.H index 722f6721..0cf07576 100644 --- a/applications/solvers/combustion/plasmaReactingFoam/YEqn.H +++ b/applications/solvers/combustion/plasmaReactingFoam/YEqn.H @@ -49,6 +49,81 @@ tmp > mvConvection 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]; @@ -60,25 +135,26 @@ tmp > mvConvection Udrift = - linearInterpolate(mue*E/ng); ve = (Udrift & mesh.Sf()) + q; - const surfaceScalarField &msf = mesh.magSf(); - // Wall electron flux correction forAll (wallPatcheIDs, pidx) { label patchID = wallPatcheIDs[pidx]; - scalar pReflex = wallReflexes[pidx]; // Probability of electron reflex + // 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)); + 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 * max((1.0-max(pReflex,1.0)),1.0) / 4.0; + wallFlux += vt * wallMSf; + wallFlux *= (1.0-pReflex); } tmp electronR( @@ -89,8 +165,8 @@ tmp > mvConvection fvScalarMatrix neEqn ( fvm::ddt(ne) - + mvConvection->fvmDiv(ve, ne) - - mvConvection->fvmDiv(fvc::interpolate(De/ng/Te*fvc::grad(Te)) & mesh.Sf(), ne) + + fvm::div(ve, ne) + // - mvConvection->fvmDiv(fvc::interpolate(De/ng/Te*fvc::grad(Te)) & mesh.Sf(), ne) - fvm::laplacian(De/ng, ne) == electronR @@ -120,13 +196,31 @@ tmp > mvConvection 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, Yi) - : mvConvection->fvmDiv(phi_drift, 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) diff --git a/applications/solvers/combustion/plasmaReactingFoam/createFields.H b/applications/solvers/combustion/plasmaReactingFoam/createFields.H index 0ace61f9..ebc80d0a 100644 --- a/applications/solvers/combustion/plasmaReactingFoam/createFields.H +++ b/applications/solvers/combustion/plasmaReactingFoam/createFields.H @@ -64,38 +64,6 @@ scalar TeFac ( ); -dictionary wallElectronFlux -( - physicalProperties.subDict("wallElectronFlux") -); - -word TeName(wallElectronFlux.lookup("TeName")); - -wordList wallPatcheNames (wallElectronFlux.lookup("wallPatches")); -labelList wallPatcheIDs (wallPatcheNames.size(), 0); -scalarList wallReflexes (wallElectronFlux.lookup("wallReflexes")); - -forAll (wallPatcheNames, pi) -{ - word patchName = wallPatcheNames[pi]; - label patchID = mesh.boundaryMesh().findPatchID(patchName); - wallPatcheIDs[pi] = patchID; - -/* - Info<< patchName << patchID << endl; - std::cout << Pstream::myProcNo() << patchName << patchID << std::endl; - - OStringStream temp_ss; - temp_ss << Pstream::myProcNo() << mesh.boundaryMesh(); - std::cout << temp_ss.str() << endl; -*/ -} - -Info<< TeName << endl; -Info<< wallPatcheNames << endl; -Info<< wallPatcheIDs << endl; - - Info<< "Reading field Phi\n" << endl; volScalarField Phi ( @@ -111,7 +79,17 @@ volScalarField Phi ); Info<< "Creating field electric field\n" << endl; -volVectorField E("E", -fvc::grad(Phi)); +volVectorField E( + IOobject + ( + "E", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + -fvc::grad(Phi) +); Info<< "Creating reaction model\n" << endl; @@ -379,3 +357,68 @@ surfaceScalarField phi_drift ), phi ); + + +surfaceScalarField phi_neutral +( + IOobject + ( + "phi_neutral", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + phi +); + +// plasmaWallFluxes + +// electron wall flux +dictionary wallElectronFlux +( + physicalProperties.subDict("wallElectronFlux") +); + +word TeName(wallElectronFlux.lookup("TeName")); + +wordList wallPatcheNames (wallElectronFlux.lookup("wallPatches")); +labelList wallPatcheIDs (wallPatcheNames.size(), 0); +scalarList wallReflexes (wallElectronFlux.lookup("wallReflexes")); + +forAll (wallPatcheNames, pi) +{ + const word patchName = wallPatcheNames[pi]; + wallPatcheIDs[pi] + = mesh.boundaryMesh().findPatchID(patchName); +} + +// ion wall flux + +dictionary wallIonFluxes +( + physicalProperties.subDict("wallIonFluxes") +); + +const hashedWordList ions(wordList(wallIonFluxes.lookup("ions"))); +wordList neutrals_ (ions.size(), ""); +PtrList reflexes (ions.size()); +PtrList ionFluxBFs (ions.size()); +PtrList neutralFluxBFs (ions.size()); + +forAll (ions, iidx) +{ + const dictionary &wallIonFlux = wallIonFluxes.subDict(ions[iidx]); + + neutrals_[iidx] = word(wallIonFlux.lookup("neutral")); + reflexes.set(iidx, + new scalarList(wallIonFlux.lookup("wallReflexes"))); + ionFluxBFs.set(iidx, + new surfaceScalarField::GeometricBoundaryField + (phi.boundaryField())); + neutralFluxBFs.set(iidx, + new surfaceScalarField::GeometricBoundaryField + (phi.boundaryField())); +} + +const hashedWordList neutrals(neutrals_);