research/OpenFOAM-2.4.x
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

plasmaReactingFoam subcycling

Browse source
This commit is contained in:
ignis 2016-11-27 13:08:38 +09:00
parent 4406424391
commit 44b3dbfabd
8 changed files with 134 additions and 91 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

2
applications/solvers/combustion/plasmaReactingFoam/PhiEqn.H
View file

@ -9,8 +9,6 @@
E = -fvc::grad(Phi);
ng = p / T * (NA / R);
tmp<volScalarField> tMagE (mag(E));
const volScalarField &magE = tMagE();

2
applications/solvers/combustion/plasmaReactingFoam/UEqn.H
View file

@ -20,3 +20,5 @@
fvOptions.correct(U);
K = 0.5*magSqr(U);
}
q = linearInterpolate(U) & mesh.Sf();

89
applications/solvers/combustion/plasmaReactingFoam/YEqn.H
View file

@ -10,48 +10,11 @@ tmp<fv::convectionScheme<scalar> > mvConvection
);
{
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();
@ -132,58 +95,6 @@ tmp<fv::convectionScheme<scalar> > mvConvection
if (Y[i].name() == electronSpecie)
{
Udrift = -(mue/ng)*E;
Uthermal = -((De/ng/Te)*fvc::grad(Te));
ve = (linearInterpolate(Udrift+Uthermal) & 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)
{

4
applications/solvers/combustion/plasmaReactingFoam/neControls.H Normal file
View file

@ -0,0 +1,4 @@
const dictionary& neControls = mesh.solverDict(ne.name());
label nNeSubCycles(readLabel(neControls.lookup("nNeSubCycles")));

92
applications/solvers/combustion/plasmaReactingFoam/neEqn.H Normal file
View file

@ -0,0 +1,92 @@
{
// Electron swarm parameter
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();
Udrift = -(mue/ng)*E;
Uthermal = -((De/ng/Te)*fvc::grad(Te));
ve = (linearInterpolate(Udrift+Uthermal) & mesh.Sf()) + q;
const surfaceScalarField &msf = mesh.magSf();
// 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);
}
volScalarField& Yi = composition.Y(electronSpecie);
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);
}

22
applications/solvers/combustion/plasmaReactingFoam/neEqnSubCycle.H Normal file
View file

@ -0,0 +1,22 @@
if (nNeSubCycles > 1)
{
dimensionedScalar totalDeltaT = runTime.deltaT();
for
(
subCycle<volScalarField> neSubCycle(ne, nNeSubCycles);
!(++neSubCycle).end();
)
{
#include "neEqn.H"
#include "PhiEqn.H"
// rhoPhiSum += (runTime.deltaT()/totalDeltaT)*rhoPhi;
}
// rhoPhi = rhoPhiSum;
}
else
{
#include "neEqn.H"
#include "PhiEqn.H"
}

3
applications/solvers/combustion/plasmaReactingFoam/numberDensity.H Normal file
View file

@ -0,0 +1,3 @@
{
ng = p / T * (NA / R);
}

11
applications/solvers/combustion/plasmaReactingFoam/plasmaReactingFoam.C
View file

@ -30,6 +30,7 @@ Description
\*---------------------------------------------------------------------------*/
#include "fvCFD.H"
#include "subCycle.H"
#include "turbulenceModel.H"
#include "psiCombustionModel.H"
#include "multivariateScheme.H"
@ -71,12 +72,22 @@ 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 "UEqn.H"
reaction->correct();
dQ = reaction->dQ();
#include "neControls.H"
#include "neEqnSubCycle.H"
// #include "neEqn.H"
#include "YEqn.H"
#include "EEqn.H"