From 70a5b812e5c75a6fc60fece19492b334f612e715 Mon Sep 17 00:00:00 2001 From: ignis Date: Wed, 28 Sep 2016 17:44:13 +0900 Subject: [PATCH] =?UTF-8?q?Working=20plasmaChemFoam=20YEqn=20-=20electron?= =?UTF-8?q?=20number=20density=20equation=20createBaseFields=20-=20correct?= =?UTF-8?q?ed=20dimension=20of=20vsf=20ne=20createFields=20=20=20=20=20-?= =?UTF-8?q?=20Reads=20description=20of=20EnProfile=20(=20units=20and=20tim?= =?UTF-8?q?e=20step=20size=20)=20=20=20=20=20-=20vsf=20En=20added=20=20=20?= =?UTF-8?q?=20=20-=20Some=20constatns=20(me,=20e,=20NA,=20R,=20pi,=20kB,?= =?UTF-8?q?=20=E2=80=A6)=20=20=20=20=20-=20Reads=20name=20list=20of=20fiel?= =?UTF-8?q?ds=20to=20write=20to=20post=20file=20output=20=20=20=20=20-=20v?= =?UTF-8?q?erbose=20info=20to=20stdout=20log=20=20=20=20=20-=20Write=20man?= =?UTF-8?q?y=20fields=20to=20post=20field=20updateEn=20-=20update=20En=20v?= =?UTF-8?q?alue=20according=20to=20give=20time=20varying=20profile=20(also?= =?UTF-8?q?=20Te)=20solveChemistry=20-=20chemistryModel.solve=20to=20.calc?= =?UTF-8?q?ulate=20(electron=20reaction=20heat=20release=20correction)=20s?= =?UTF-8?q?etDeltaT=20-=20change=20time=20step=20size=20after=20the=20puls?= =?UTF-8?q?e=20phase?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- .../combustion/chemFoam/plasmaChemFoam/YEqn.H | 30 ++++++ .../plasmaChemFoam/createBaseFields.H | 2 +- .../chemFoam/plasmaChemFoam/createFields.H | 101 +++++++++++++++++- .../chemFoam/plasmaChemFoam/output.H | 48 +++++++++ .../chemFoam/plasmaChemFoam/plasmaChemFoam.C | 2 + .../chemFoam/plasmaChemFoam/readControls.H | 10 ++ .../chemFoam/plasmaChemFoam/setDeltaT.H | 11 ++ .../chemFoam/plasmaChemFoam/solveChemistry.H | 4 + .../chemFoam/plasmaChemFoam/updateEn.H | 9 ++ 9 files changed, 213 insertions(+), 4 deletions(-) create mode 100644 applications/solvers/combustion/chemFoam/plasmaChemFoam/YEqn.H create mode 100644 applications/solvers/combustion/chemFoam/plasmaChemFoam/output.H create mode 100644 applications/solvers/combustion/chemFoam/plasmaChemFoam/readControls.H create mode 100644 applications/solvers/combustion/chemFoam/plasmaChemFoam/setDeltaT.H create mode 100644 applications/solvers/combustion/chemFoam/plasmaChemFoam/solveChemistry.H create mode 100644 applications/solvers/combustion/chemFoam/plasmaChemFoam/updateEn.H diff --git a/applications/solvers/combustion/chemFoam/plasmaChemFoam/YEqn.H b/applications/solvers/combustion/chemFoam/plasmaChemFoam/YEqn.H new file mode 100644 index 00000000..717165a8 --- /dev/null +++ b/applications/solvers/combustion/chemFoam/plasmaChemFoam/YEqn.H @@ -0,0 +1,30 @@ +{ + forAll(Y, specieI) + { + volScalarField& Yi = Y[specieI]; + + if (Yi.name() == "E-") + { + tmp 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") + ); + } + } + +} diff --git a/applications/solvers/combustion/chemFoam/plasmaChemFoam/createBaseFields.H b/applications/solvers/combustion/chemFoam/plasmaChemFoam/createBaseFields.H index 471377b9..d966a6bf 100644 --- a/applications/solvers/combustion/chemFoam/plasmaChemFoam/createBaseFields.H +++ b/applications/solvers/combustion/chemFoam/plasmaChemFoam/createBaseFields.H @@ -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(); diff --git a/applications/solvers/combustion/chemFoam/plasmaChemFoam/createFields.H b/applications/solvers/combustion/chemFoam/plasmaChemFoam/createFields.H index b9d99ef5..c4d63282 100644 --- a/applications/solvers/combustion/chemFoam/plasmaChemFoam/createFields.H +++ b/applications/solvers/combustion/chemFoam/plasmaChemFoam/createFields.H @@ -26,6 +26,23 @@ // read CSV CSV 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 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 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(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; + } diff --git a/applications/solvers/combustion/chemFoam/plasmaChemFoam/output.H b/applications/solvers/combustion/chemFoam/plasmaChemFoam/output.H new file mode 100644 index 00000000..55876b54 --- /dev/null +++ b/applications/solvers/combustion/chemFoam/plasmaChemFoam/output.H @@ -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(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; + } + diff --git a/applications/solvers/combustion/chemFoam/plasmaChemFoam/plasmaChemFoam.C b/applications/solvers/combustion/chemFoam/plasmaChemFoam/plasmaChemFoam.C index e41ba4ea..98f686d1 100644 --- a/applications/solvers/combustion/chemFoam/plasmaChemFoam/plasmaChemFoam.C +++ b/applications/solvers/combustion/chemFoam/plasmaChemFoam/plasmaChemFoam.C @@ -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" diff --git a/applications/solvers/combustion/chemFoam/plasmaChemFoam/readControls.H b/applications/solvers/combustion/chemFoam/plasmaChemFoam/readControls.H new file mode 100644 index 00000000..8570175e --- /dev/null +++ b/applications/solvers/combustion/chemFoam/plasmaChemFoam/readControls.H @@ -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"))); diff --git a/applications/solvers/combustion/chemFoam/plasmaChemFoam/setDeltaT.H b/applications/solvers/combustion/chemFoam/plasmaChemFoam/setDeltaT.H new file mode 100644 index 00000000..dfc3fb43 --- /dev/null +++ b/applications/solvers/combustion/chemFoam/plasmaChemFoam/setDeltaT.H @@ -0,0 +1,11 @@ +if (runTime.timeOutputValue() < 2 * EnProfileEndTime) +{ + runTime.setDeltaT(EnDeltaT1); +} +else +{ + runTime.setDeltaT(deltaT); +} + +Info<< "deltaT = " << runTime.deltaT().value() << endl; + diff --git a/applications/solvers/combustion/chemFoam/plasmaChemFoam/solveChemistry.H b/applications/solvers/combustion/chemFoam/plasmaChemFoam/solveChemistry.H new file mode 100644 index 00000000..6bf4848a --- /dev/null +++ b/applications/solvers/combustion/chemFoam/plasmaChemFoam/solveChemistry.H @@ -0,0 +1,4 @@ + // dtChem = chemistry.solve(runTime.deltaT().value()); + chemistry.calculate(); + scalar Sh = chemistry.Sh()()[0]/rho[0]; + integratedHeat += Sh*runTime.deltaT().value(); diff --git a/applications/solvers/combustion/chemFoam/plasmaChemFoam/updateEn.H b/applications/solvers/combustion/chemFoam/plasmaChemFoam/updateEn.H new file mode 100644 index 00000000..5b1dd850 --- /dev/null +++ b/applications/solvers/combustion/chemFoam/plasmaChemFoam/updateEn.H @@ -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();