plasmaReactingFoam - snGrad of E for boundary drift calc / Disable mvConvection RTS / Separate surfaceFields for ion fluxes / Write electron flux and ions fluxes

This commit is contained in:
ignis 2016-12-14 20:49:28 +09:00
parent ea92642329
commit 6370bb0f5d
4 changed files with 70 additions and 18 deletions

View file

@ -8,6 +8,7 @@
);
E = -fvc::grad(Phi);
snE = -fvc::snGrad(Phi);
tmp<volScalarField> tMagE (mag(E));
const volScalarField &magE = tMagE();

View file

@ -1,3 +1,4 @@
/*
tmp<fv::convectionScheme<scalar> > mvConvection
(
fv::convectionScheme<scalar>::New
@ -8,6 +9,7 @@ tmp<fv::convectionScheme<scalar> > mvConvection
mesh.divScheme("div(phi,Yi_h)")
)
);
*/
{
label inertIndex = -1;
@ -35,9 +37,11 @@ tmp<fv::convectionScheme<scalar> > mvConvection
// Adding drift flux to boundary patches
forAll (bfIonFlux, pidx)
{
Info << "Adding drift flux to boundary patches" << pidx << endl;
bfIonFlux[pidx] +=
(E.boundaryField()[pidx]
& sf.boundaryField()[pidx])
snE.boundaryField()[pidx]
* msf.boundaryField()[pidx]
* rho.boundaryField()[pidx]
* Di.boundaryField()[pidx]
/ T.boundaryField()[pidx]
@ -113,22 +117,22 @@ tmp<fv::convectionScheme<scalar> > mvConvection
if (nCharge != 0)
{
phi_drift = phi;
phi_drift += fvc::interpolate((rho*Di/T*(eCharge*z/kB))*E) & mesh.Sf();
phis[i] = phi;
phis[i] += 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];
// phis[i] updated
phis[i].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];
// update phis[i]
phis[i].internalField() = phi.internalField();
phis[i].boundaryField() = neutralFluxBFs[ibc];
}
@ -137,9 +141,9 @@ tmp<fv::convectionScheme<scalar> > mvConvection
fvm::ddt(rho, Yi)
+
( nCharge != 0
? fvm::div(phi_drift, Yi, "div(phi,Yi_h)")
? fvm::div(phis[i], Yi, "div(phi,Yi_h)")
: ( neutrals.contains(Y[i].name())
? fvm::div(phi_neutral, Yi, "div(phi,Yi_h)")
? fvm::div(phis[i], Yi, "div(phi,Yi_h)")
: fvm::div(phi, Yi, "div(phi,Yi_h)")
)
)

View file

@ -92,6 +92,20 @@ volVectorField E
-fvc::grad(Phi)
);
surfaceScalarField snE
(
IOobject
(
"snE",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
-fvc::snGrad(Phi)
);
Info<< "Creating reaction model\n" << endl;
autoPtr<combustionModels::psiCombustionModel> reaction
@ -363,8 +377,8 @@ surfaceScalarField ve
"ve",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
linearInterpolate(Udrift+Uthermal) & mesh.Sf()
);
@ -377,7 +391,7 @@ surfaceScalarField phi_drift
"phi_drift",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::READ_IF_PRESENT,
IOobject::NO_WRITE
),
phi
@ -391,13 +405,42 @@ surfaceScalarField phi_neutral
"phi_neutral",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::READ_IF_PRESENT,
IOobject::NO_WRITE
),
phi
);
PtrList<surfaceScalarField> phis (composition.species().size());
forAll (composition.species(), isp)
{
const scalar z(composition.z(isp));
const label nCharge(z);
if (nCharge != 0)
{
phis.set(isp,
new surfaceScalarField
(
IOobject
(
"phi." + composition.species()[isp],
runTime.timeName(),
mesh,
IOobject::READ_IF_PRESENT,
IOobject::AUTO_WRITE
),
phi
)
);
}
}
// plasmaWallFluxes
Info<< "Reading plasma wall flux bc control\n" << endl;
// electron wall flux
dictionary wallElectronFlux
@ -418,6 +461,8 @@ forAll (wallPatcheNames, pi)
= mesh.boundaryMesh().findPatchID(patchName);
}
Info<< "plasma walls are \n" << wallPatcheNames << endl;
// ion wall flux
dictionary wallIonFluxes
@ -432,6 +477,8 @@ PtrList<hashedWordList> targetList (ions.size());
PtrList<scalarList> reflexes (ions.size());
PtrList<surfaceScalarField::GeometricBoundaryField> ionFluxBFs (ions.size());
Info<< ions.size() << " ions are \n" << ions << endl;
label nTargets = 0;
forAll (ions, iidx)
@ -485,7 +532,6 @@ for (label nidx = 0; nidx < nTargets; nidx++)
neutrals.append(n);
nNeutrals += 1;
}
}
PtrList<surfaceScalarField::GeometricBoundaryField> neutralFluxBFs (neutrals.size());

View file

@ -72,13 +72,14 @@ int main(int argc, char *argv[])
runTime++;
Info<< "Time = " << runTime.timeName() << nl << endl;
#include "numberDensity.H"
#include "PhiEqn.H"
#include "rhoEqn.H"
while (pimple.loop())
{
#include "numberDensity.H"
#include "PhiEqn.H"
#include "UEqn.H"
reaction->correct();