ion absorption/neutralization flux at bouundaries
This commit is contained in:
ignis
parent
3eec143b5a
commit
150202d4e9
2 changed files with 180 additions and 43 deletions
|
|
@ -49,6 +49,81 @@ tmp<fv::convectionScheme<scalar> > 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<fv::convectionScheme<scalar> > 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<fvScalarMatrix> electronR(
|
||||
|
|
@ -89,8 +165,8 @@ tmp<fv::convectionScheme<scalar> > 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<fv::convectionScheme<scalar> > 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)
|
||||
|
|
|
|||
|
|
@ -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<scalarList> reflexes (ions.size());
|
||||
PtrList<surfaceScalarField::GeometricBoundaryField> ionFluxBFs (ions.size());
|
||||
PtrList<surfaceScalarField::GeometricBoundaryField> 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_);
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue