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

Working plasmaChemFoam

Browse source 
YEqn - electron number density equation
createBaseFields - corrected dimension of vsf ne
createFields
    - Reads description of EnProfile ( units and time step size )
    - vsf En added
    - Some constatns (me, e, NA, R, pi, kB, …)
    - Reads name list of fields to write to post file
output
    - verbose info to stdout log
    - Write many fields to post field
updateEn - update En value according to give time varying profile (also Te)
solveChemistry - chemistryModel.solve to .calculate (electron reaction heat release correction)
setDeltaT - change time step size after the pulse phase
This commit is contained in:
ignis 2016-09-28 17:44:13 +09:00
parent 85596642b8
commit 70a5b812e5
9 changed files with 213 additions and 4 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

30
applications/solvers/combustion/chemFoam/plasmaChemFoam/YEqn.H Normal file
View file

@ -0,0 +1,30 @@
{
forAll(Y, specieI)
{
volScalarField& Yi = Y[specieI];
if (Yi.name() == "E-")
{
tmp<fvScalarMatrix> electronR(
new fvScalarMatrix(ne, dimless/dimTime));
electronR->source() -= mesh.V()*chemistry.RR(specieI).field();
solve
(
// fvm::ddt(ne) - chemistry.RR(specieI),
fvm::ddt(ne) - electronR,
mesh.solver("Yi")
);
}
else
{
solve
(
fvm::ddt(rho, Yi) - chemistry.RR(specieI),
mesh.solver("Yi")
);
}
}
}

2
applications/solvers/combustion/chemFoam/plasmaChemFoam/createBaseFields.H
View file

@ -66,7 +66,7 @@ Info<< "Creating base fields for time " << runTime.timeName() << endl;
false
),
mesh,
dimensionedScalar("ne", dimTemperature, ne0)
dimensionedScalar("ne", dimless/dimVolume, ne0)
);
ne.write();

101
applications/solvers/combustion/chemFoam/plasmaChemFoam/createFields.H
View file

@ -26,6 +26,23 @@
// read CSV
CSV<scalar> EnProfile("EnProfile", initialConditions);
// Unit conversion factor of time in E/n profile table (to seconds)
// 1e-9 if time values in the table are nanoseconds
scalar EnProfileXScale = readScalar(initialConditions.lookup("EnProfileXScale"));
// Unit conversion factor of reduced electric field in E/n profile table (to V m2)
// 1e-21 if time values in the table are Townsends
scalar EnProfileYScale = readScalar(initialConditions.lookup("EnProfileYScale"));
// Time when the pulse ends (last time entry)
scalar EnProfileEndTime = EnProfile.x()->last() * EnProfileXScale;
// Time step size used during the pulse and afterglow
scalar EnDeltaT1 = readScalar(initialConditions.lookup("EnDeltaT1"));
// scalar EnDeltaT2 = readScalar(initialConditions.lookup("EnDeltaT2"));
Info<< "Pulse ends at = " << EnProfileEndTime << endl;
#include "createBaseFields.H"
volScalarField ne
@ -54,6 +71,56 @@
mesh
);
volScalarField En
(
IOobject
(
"En",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
mesh,
dimensionedScalar("En", dimPower/dimCurrent/dimLength*dimVol, EnProfile.value(runTime.timeOutputValue()/EnProfileXScale))
);
Info<< "Reading physicalProperties\n" << endl;
IOdictionary physicalProperties
(
IOobject
(
"physicalProperties",
runTime.constant(),
mesh,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE
)
);
CSV<scalar> TeOfEn("TeOfEn", physicalProperties);
//- Electron mass (default in [kg])
const dimensionedScalar eMass = constant::atomic::me;
//- Elementary charge (default in [C])
const dimensionedScalar eCharge = constant::electromagnetic::e;
//- Avogadro number (default in [1/mol])
const dimensionedScalar NA = constant::physicoChemical::NA;
//- Universal gas constant (default in [J/mol/K])
const dimensionedScalar R = constant::physicoChemical::R;
//- Pi
const scalar pi = constant::mathematical::pi;
//- Boltzmann constant
const scalar kB = constant::physicoChemical::k.value();
Info<< nl << "Reading thermophysicalProperties" << endl;
autoPtr<psiChemistryModel> pChemistry(psiChemistryModel::New(mesh));
@ -129,8 +196,36 @@
)
);
OFstream post(args.path()/"chemFoam.out");
post<< "# Time" << token::TAB << "Temperature [K]" << token::TAB
<< "Pressure [Pa]" << endl;
wordList watchList (initialConditions.lookup("watchList"));
wordList RRWatchList (initialConditions.lookup("RRWatchList"));
OFstream post(args.path()/"plasmaChemFoam.out");
{
post << "t" << token::TAB
<< "T" << token::TAB
<< "p" << token::TAB
<< "ne" << token::TAB;
forAll (watchList, iWatch)
{
word specieName = watchList[iWatch];
// test whether species name is valid.
// const volScalarField &cvf = composition.Y(specieName);
const volScalarField &cvf = mesh.lookupObject<volScalarField>(specieName);
post << cvf.name() << token::TAB;
}
forAll (RRWatchList, iWatch)
{
word specieName = RRWatchList[iWatch];
// test whether species name is valid.
// const volScalarField &cvf = composition.Y(specieName);
const label specieI = composition.species()[specieName];
post << chemistry.RR(specieI).name() << token::TAB;
}
post << endl;
}

48
applications/solvers/combustion/chemFoam/plasmaChemFoam/output.H Normal file
View file

@ -0,0 +1,48 @@
runTime.write();
Info<< "Sh = " << Sh
<< ", T = " << thermo.T()[0]
<< ", p = " << thermo.p()[0]
<< ", E/n = " << En[0]
<< ", Te = " << Te[0]
<< ", ne = " << ne[0]
<< ", " << composition.Y("CH4").name()
<< " = " << composition.Y("CH4")[0]
<< ", " << composition.Y("O2").name()
<< " = " << composition.Y("O2")[0]
<< ", " << composition.Y("CH3").name()
<< " = " << composition.Y("CH3")[0]
<< ", " << composition.Y("O").name()
<< " = " << composition.Y("O")[0]
<< ", " << composition.Y("H").name()
<< " = " << composition.Y("H")[0]
<< endl;
{
post << runTime.value() << token::TAB
<< thermo.T()[0] << token::TAB
<< thermo.p()[0] << token::TAB
<< ne[0] << token::TAB;
forAll (watchList, iWatch)
{
word specieName = watchList[iWatch];
// test whether species name is valid.
// const volScalarField &cvf = composition.Y(specieName);
const scalarField &csf = mesh.lookupObject<volScalarField>(specieName);
post << csf[0] << token::TAB;
}
forAll (RRWatchList, iWatch)
{
word specieName = RRWatchList[iWatch];
// test whether species name is valid.
const label specieI = composition.species()[specieName];
const scalarField &csf = chemistry.RR(specieI);
post << csf[0] << token::TAB;
}
post << endl;
}

2
applications/solvers/combustion/chemFoam/plasmaChemFoam/plasmaChemFoam.C
View file

@ -69,7 +69,9 @@ int main(int argc, char *argv[])
runTime++;
Info<< "Time = " << runTime.timeName() << nl << endl;
Info<< "E/n = " << EnProfile.value(runTime.timeOutputValue()/EnProfileXScale) << nl << endl;
#include "updateEn.H"
#include "solveChemistry.H"
#include "YEqn.H"
#include "hEqn.H"

10
applications/solvers/combustion/chemFoam/plasmaChemFoam/readControls.H Normal file
View file

@ -0,0 +1,10 @@
if (runTime.controlDict().lookupOrDefault("suppressSolverInfo", false))
{
lduMatrix::debug = 0;
}
Switch adjustTimeStep(runTime.controlDict().lookup("adjustTimeStep"));
scalar maxDeltaT(readScalar(runTime.controlDict().lookup("maxDeltaT")));
scalar deltaT(readScalar(runTime.controlDict().lookup("deltaT")));

11
applications/solvers/combustion/chemFoam/plasmaChemFoam/setDeltaT.H Normal file
View file

@ -0,0 +1,11 @@
if (runTime.timeOutputValue() < 2 * EnProfileEndTime)
{
runTime.setDeltaT(EnDeltaT1);
}
else
{
runTime.setDeltaT(deltaT);
}
Info<< "deltaT = " << runTime.deltaT().value() << endl;

4
applications/solvers/combustion/chemFoam/plasmaChemFoam/solveChemistry.H Normal file
View file

@ -0,0 +1,4 @@
// dtChem = chemistry.solve(runTime.deltaT().value());
chemistry.calculate();
scalar Sh = chemistry.Sh()()[0]/rho[0];
integratedHeat += Sh*runTime.deltaT().value();

9
applications/solvers/combustion/chemFoam/plasmaChemFoam/updateEn.H Normal file
View file

@ -0,0 +1,9 @@
const scalar newEn = EnProfile.value(runTime.timeOutputValue()/EnProfileXScale);
forAll(En, celli)
{
En[celli] = newEn * EnProfileYScale;
Te[celli] = eCharge.value() / kB * TeOfEn.value(newEn);
Te[celli] = max(Te[celli], thermo.T()[celli]);
}
En.correctBoundaryConditions();
Te.correctBoundaryConditions();