From 9949acea5edb33c440867133c9ef08debff37b42 Mon Sep 17 00:00:00 2001 From: ignis Date: Mon, 15 May 2017 12:26:05 +0900 Subject: [PATCH] eReactingFoam copy of plasmaReactingFoam --- Allwmake | 7 + EEqn.H | 36 +++ Make/files | 3 + Make/options | 27 +++ PhiEqn.H | 19 ++ UEqn.H | 24 ++ YEqn.H | 181 +++++++++++++++ createFields.H | 548 ++++++++++++++++++++++++++++++++++++++++++++ eReactingFoam.C | 120 ++++++++++ electronCourantNo.H | 51 +++++ neControls.H | 4 + neEqn.H | 96 ++++++++ neEqnSubCycle.H | 22 ++ numberDensity.H | 3 + pEqn.H | 90 ++++++++ readTimeControls.H | 70 ++++++ setDeltaT.H | 85 +++++++ 17 files changed, 1386 insertions(+) create mode 100755 Allwmake create mode 100644 EEqn.H create mode 100644 Make/files create mode 100644 Make/options create mode 100644 PhiEqn.H create mode 100644 UEqn.H create mode 100644 YEqn.H create mode 100644 createFields.H create mode 100644 eReactingFoam.C create mode 100644 electronCourantNo.H create mode 100644 neControls.H create mode 100644 neEqn.H create mode 100644 neEqnSubCycle.H create mode 100644 numberDensity.H create mode 100644 pEqn.H create mode 100644 readTimeControls.H create mode 100644 setDeltaT.H diff --git a/Allwmake b/Allwmake new file mode 100755 index 0000000..1d9987e --- /dev/null +++ b/Allwmake @@ -0,0 +1,7 @@ +#!/bin/sh +cd ${0%/*} || exit 1 # run from this directory +set -x + +wmake + +# ----------------------------------------------------------------- end-of-file diff --git a/EEqn.H b/EEqn.H new file mode 100644 index 0000000..9262b91 --- /dev/null +++ b/EEqn.H @@ -0,0 +1,36 @@ +{ + volScalarField& he = thermo.he(); + + fvScalarMatrix EEqn + ( + fvm::ddt(rho, he) + fvm::div(phi, he, "div(phi,Yi_h)") + + fvc::ddt(rho, K) + fvc::div(phi, K) + + ( + he.name() == "e" + ? fvc::div + ( + fvc::absolute(phi/fvc::interpolate(rho), U), + p, + "div(phiv,p)" + ) + : -dpdt + ) + - fvm::laplacian(turbulence->alphaEff(), he) + == + reaction->Sh() + + fvOptions(rho, he) + ); + + EEqn.relax(); + + fvOptions.constrain(EEqn); + + EEqn.solve(); + + fvOptions.correct(he); + + thermo.correct(); + + Info<< "min/max(T) = " + << min(T).value() << ", " << max(T).value() << endl; +} diff --git a/Make/files b/Make/files new file mode 100644 index 0000000..de06199 --- /dev/null +++ b/Make/files @@ -0,0 +1,3 @@ +eReactingFoam.C + +EXE = $(FOAM_USER_APPBIN)/eReactingFoam diff --git a/Make/options b/Make/options new file mode 100644 index 0000000..15cf68c --- /dev/null +++ b/Make/options @@ -0,0 +1,27 @@ +EXE_INC = \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/fvOptions/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude \ + -I$(LIB_SRC)/sampling/lnInclude \ + -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \ + -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \ + -I$(LIB_SRC)/ODE/lnInclude \ + -I$(LIB_SRC)/combustionModels/lnInclude + +EXE_LIBS = \ + -lfiniteVolume \ + -lfvOptions \ + -lmeshTools \ + -lsampling \ + -lcompressibleTurbulenceModel \ + -lcompressibleRASModels \ + -lcompressibleLESModels \ + -lreactionThermophysicalModels \ + -lspecie \ + -lfluidThermophysicalModels \ + -lchemistryModel \ + -lODE \ + -lcombustionModels diff --git a/PhiEqn.H b/PhiEqn.H new file mode 100644 index 0000000..8a912c7 --- /dev/null +++ b/PhiEqn.H @@ -0,0 +1,19 @@ +{ + rhoq = (rho * qc) - (eCharge * ne); + + solve + ( + fvm::laplacian(Phi) + + rhoq/epsilon0 + ); + + E = -fvc::grad(Phi); + snE = -fvc::snGrad(Phi); + + tmp tMagE (mag(E)); + const volScalarField &magE = tMagE(); + + magE.writeMinMax(Info); + + En = magE / (ng); +} diff --git a/UEqn.H b/UEqn.H new file mode 100644 index 0000000..6dccc4d --- /dev/null +++ b/UEqn.H @@ -0,0 +1,24 @@ + fvVectorMatrix UEqn + ( + fvm::ddt(rho, U) + + fvm::div(phi, U) + + turbulence->divDevRhoReff(U) + == + rho*g + + rhoq*E + + fvOptions(rho, U) + ); + + UEqn.relax(); + + fvOptions.constrain(UEqn); + + if (pimple.momentumPredictor()) + { + solve(UEqn == -fvc::grad(p)); + + fvOptions.correct(U); + K = 0.5*magSqr(U); + } + + q = linearInterpolate(U) & mesh.Sf(); diff --git a/YEqn.H b/YEqn.H new file mode 100644 index 0000000..ac1e073 --- /dev/null +++ b/YEqn.H @@ -0,0 +1,181 @@ +/* +tmp > mvConvection +( + fv::convectionScheme::New + ( + mesh, + fields, + phi, + mesh.divScheme("div(phi,Yi_h)") + ) +); +*/ + +{ + label inertIndex = -1; + volScalarField Yt(0.0*Y[0]); + + composition.calculateDiffusivities(p, T); + + const surfaceScalarField &msf = mesh.magSf(); + const surfaceVectorField &sf = mesh.Sf(); + + forAll(ions, k) // ion-neutral pair + { + const word nIon(ions[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]; + bfIonFlux = phi.boundaryField(); + + // Adding drift flux to boundary patches + forAll (bfIonFlux, pidx) + { + Info << "Adding drift flux to boundary patches" << pidx << endl; + + bfIonFlux[pidx] += + snE.boundaryField()[pidx] + * msf.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 MIon(WIon / NA.value() / 1000.0); + + const volScalarField& Yion = composition.Y(nIon); + + const hashedWordList &targets = targetList[k]; + + forAll (targets, tidx) + { + const word &nNeu(targets[tidx]); + + surfaceScalarField::GeometricBoundaryField & + bfNeuFlux = neutralFluxBFs[neutrals[nNeu]]; + bfNeuFlux = phi.boundaryField(); + + const scalar WNeu(composition.W(composition.species()[nNeu])); + const scalar MNeu(WNeu / NA.value() / 1000.0); + + const volScalarField& Yneu = composition.Y(nNeu); + + forAll(wallPatcheIDs, pidx) // loop over wall patches + { + const label patchID = wallPatcheIDs[pidx]; + + // Probability of ion reflex + const scalar pReflex = max(min(rK[pidx],1.0),0.0); + + scalarField &wallFluxIon = bfIonFlux[patchID]; + + const scalarField &wallMSf = msf.boundaryField()[patchID]; + const scalarField &wallT = T.boundaryField()[patchID]; + const scalarField &wallYion = Yion.boundaryField()[patchID]; + + const 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 + scalarField &wallFluxNeu = bfNeuFlux[patchID]; + + const scalarField &wallYneu = Yneu.boundaryField()[patchID]; + + wallFluxNeu -= wallFluxIon * wallYion / wallYneu + / (WIon / WNeu); + } + } + } + + forAll(Y, i) + { + volScalarField& Yi = Y[i]; + + const volScalarField& Di = D[i]; + + if (Y[i].name() == electronSpecie) + { + } + else if (Y[i].name() != inertSpecie) + { + const scalar z(composition.z(i)); + const label nCharge(z); + + if (nCharge != 0) + { + phis[i] = phi; + phis[i] += fvc::interpolate((rho*Di/T*(eCharge*z/kB))*E) & mesh.Sf(); + + if (relaxDrift < 1.0 && relaxDrift > 0.0) + { + phis[i] *= relaxDrift; + } + } + + if (ions.contains(Y[i].name())) + { + const label ibc = ions[Y[i].name()]; + // phis[i] updated + phis[i].boundaryField() = ionFluxBFs[ibc]; + } + else if (neutrals.contains(Y[i].name())) + { + const label ibc = neutrals[Y[i].name()]; + // update phis[i] + phis[i].internalField() = phi.internalField(); + phis[i].boundaryField() = neutralFluxBFs[ibc]; + } + + + fvScalarMatrix YiEqn + ( + fvm::ddt(rho, Yi) + + + ( nCharge != 0 + ? fvm::div(phis[i], Yi, "div(phi,Yi_h)") + : ( neutrals.contains(Y[i].name()) + ? fvm::div(phis[i], Yi, "div(phi,Yi_h)") + : fvm::div(phi, Yi, "div(phi,Yi_h)") + ) + ) + // - fvm::laplacian(turbulence->muEff(), Yi) + - fvm::laplacian(rho*Di, Yi) + == + reaction->R(Yi) + + fvOptions(rho, Yi) + ); + + YiEqn.relax(); + + fvOptions.constrain(YiEqn); + + YiEqn.solve(mesh.solver("Yi")); + + fvOptions.correct(Yi); + + Yi.max(0.0); + Yt += Yi; + } + else + { + inertIndex = i; + } + } + + Y[inertIndex] = scalar(1) - Yt; + Y[inertIndex].max(0.0); +} diff --git a/createFields.H b/createFields.H new file mode 100644 index 0000000..1bfe560 --- /dev/null +++ b/createFields.H @@ -0,0 +1,548 @@ +Info<< "Reading physicalProperties\n" << endl; + +IOdictionary physicalProperties +( + IOobject + ( + "physicalProperties", + runTime.constant(), + mesh, + IOobject::MUST_READ_IF_MODIFIED, + IOobject::NO_WRITE + ) +); + +scalar relaxDrift +( + physicalProperties.lookupOrDefault("relaxDrift", 1.0) +); + +dimensionedScalar epsilon0 +( + physicalProperties.lookup("epsilon0") +); + +Switch mobility_f_of_Te = physicalProperties.lookupOrDefault("mobility_f_of_Te", false); + +Switch calculateDe = physicalProperties.lookupOrDefault("calculateDe", false); + +// Convert E/n in SI unit to table unit. +// Default V m^2 => Td +scalar TeToTableUnit ( + physicalProperties.lookupOrDefault("TeToTableUnit", 1.0) + ); + +// Convert E/n in SI unit to table unit. +// Default V m^2 => Td +scalar EnToTableUnit ( + physicalProperties.lookupOrDefault("EnToTableUnit", 1.0e21) + ); + + +autoPtr< DataEntry< scalar > > pmueN ( + DataEntry::New("mueN", physicalProperties)); +const DataEntry &mueN(pmueN()); + +// Convert mu_e * n_g value from table into SI unit. +scalar mueNFac ( + physicalProperties.lookupOrDefault("mueNFac", 1.0) + ); + + +autoPtr< DataEntry< scalar > > pDeN ( + DataEntry::New("DeN", physicalProperties)); +const DataEntry &DeN(pDeN()); + +// Convert D_e * n_g value from table into SI unit. +scalar DeNFac ( + physicalProperties.lookupOrDefault("DeNFac", 1.0) + ); + + +autoPtr< DataEntry< scalar > > pTeOfEn( + DataEntry::New ("TeOfEn", physicalProperties)); +const DataEntry &TeOfEn(pTeOfEn()); + +// Convert T_e value from table into SI unit. +scalar TeFac ( + physicalProperties.lookupOrDefault("TeFac", 1.0) + ); + + +Info<< "Reading field Phi\n" << endl; +volScalarField Phi +( + IOobject + ( + "Phi", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh +); + +Info<< "Creating field electric field\n" << endl; +volVectorField E +( + IOobject + ( + "E", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + -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 reaction +( + combustionModels::psiCombustionModel::New(mesh) +); + +psiReactionThermo& thermo = reaction->thermo(); +thermo.validate(args.executable(), "h", "e"); + +volScalarField& qc = thermo.qc(); +basicMultiComponentMixture& composition = thermo.composition(); + +//- 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; + +//- Boltzmann constant +const dimensionedScalar kB = constant::physicoChemical::k; + +//- Pi +const scalar pi = constant::mathematical::pi; + +PtrList& Y = composition.Y(); +const PtrList& D = composition.D(); + +word inertSpecie(thermo.lookup("inertSpecie")); +word electronSpecie("E-"); + +volScalarField rho +( + IOobject + ( + "rho", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + thermo.rho() +); + +Info<< "Reading field U\n" << endl; +volVectorField U +( + IOobject + ( + "U", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh +); + + +volScalarField& p = thermo.p(); +const volScalarField& psi = thermo.psi(); +const volScalarField& T = thermo.T(); + +#include "compressibleCreatePhi.H" + +surfaceScalarField q +( + IOobject + ( + "q", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + linearInterpolate(U) & mesh.Sf() +); + +Info << "Creating turbulence model.\n" << nl; +autoPtr turbulence +( + compressible::turbulenceModel::New + ( + rho, + U, + phi, + thermo + ) +); + +// Set the turbulence into the reaction model +reaction->setTurbulence(turbulence()); + + +Info<< "Creating field dpdt\n" << endl; +volScalarField dpdt +( + IOobject + ( + "dpdt", + runTime.timeName(), + mesh + ), + mesh, + dimensionedScalar("dpdt", p.dimensions()/dimTime, 0) +); + +Info<< "Creating field kinetic energy K\n" << endl; +volScalarField K("K", 0.5*magSqr(U)); + + +multivariateSurfaceInterpolationScheme::fieldTable fields; + +forAll(Y, i) +{ + fields.add(Y[i]); +} +fields.add(thermo.he()); + +volScalarField dQ +( + IOobject + ( + "dQ", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + mesh, + dimensionedScalar("dQ", dimEnergy/dimTime, 0.0) +); + +Info<< "Creating field electron number density\n" << endl; +volScalarField ne +( + IOobject + ( + "ne", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh +); + + +Info<< "Creating field charge density\n" << endl; +volScalarField rhoq +( + IOobject + ( + "rhoq", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh + // (rho * qc) - (eCharge * ne) +); + +Info<< "Creating field gas number density\n" << endl; +volScalarField ng("ng", p / R / T * NA); + +Info<< "Creating field reduced electric field\n" << endl; +volScalarField En ("En", mag(E) / (ng)); +scalarField EnTd(En.internalField() * EnToTableUnit); + +Info<< "Creating field electron mobility\n" << endl; +volScalarField mue +( + IOobject + ( + "mue", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh +); + +Info<< "Creating field electron diffusivity\n" << endl; +volScalarField De +( + IOobject + ( + "De", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh +); + +Info<< "Creating field electron temperature\n" << endl; +volScalarField Te +( + IOobject + ( + "Te", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh +); + + +/* +BE_IX::Bolos bolos; + +bolos.readCollisions("./LXCat-June2013.txt"); + +bolos.presolve(); +{ + bolos.maxwellian(2); + forAll(rho, celli) + { + De[celli] = bolos.diffusivity(); + mue[celli] = bolos.mobility(); + } + De.correctBoundaryConditions(); + mue.correctBoundaryConditions(); +} +*/ + +Info<< "Calculating face flux field ve\n" << endl; +volVectorField Udrift +( + IOobject + ( + "Udrift", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + - (mue*E/ng) +); + +volVectorField Uthermal +( + IOobject + ( + "Uthermal", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + - ((De/ng/Te)*fvc::grad(Te)) +); + +surfaceScalarField ve +( + IOobject + ( + "ve", + runTime.timeName(), + mesh, + IOobject::READ_IF_PRESENT, + IOobject::AUTO_WRITE + ), + linearInterpolate(Udrift+Uthermal) & mesh.Sf() +); + + +surfaceScalarField phi_drift +( + IOobject + ( + "phi_drift", + runTime.timeName(), + mesh, + IOobject::READ_IF_PRESENT, + IOobject::NO_WRITE + ), + phi +); + + +surfaceScalarField phi_neutral +( + IOobject + ( + "phi_neutral", + runTime.timeName(), + mesh, + IOobject::READ_IF_PRESENT, + IOobject::NO_WRITE + ), + phi +); + + +PtrList 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 +( + 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); +} + +Info<< "plasma walls are \n" << wallPatcheNames << endl; + +// ion wall flux + +dictionary wallIonFluxes +( + physicalProperties.subDict("wallIonFluxes") +); + +const label nMaxTargets = 5; +const hashedWordList ions(wordList(wallIonFluxes.lookup("ions"))); +wordList neutrals_ (ions.size()*nMaxTargets); +PtrList targetList (ions.size()); +PtrList reflexes (ions.size()); +PtrList ionFluxBFs (ions.size()); + +Info<< ions.size() << " ions are \n" << ions << endl; + +label nTargets = 0; + +forAll (ions, iidx) +{ + const dictionary &wallIonFlux = wallIonFluxes.subDict(ions[iidx]); + + const Switch noTarget (wallIonFlux.lookupOrDefault("noTarget", false)); + + const wordList targets (wallIonFlux.lookupOrDefault("neutrals", wordList())); + + if (noTarget) + { + targetList.set(iidx, new hashedWordList()); + } + else if (targets.empty()) + { + word neutralName(wallIonFlux.lookup("neutral")); + neutrals_[nTargets] = neutralName; + nTargets += 1; + + targetList.set(iidx, new hashedWordList(wordList(1, neutralName))); + } + else + { + wordList::subList sub (neutrals_, targets.size(), nTargets); + forAll (targets, tidx) + { + sub[tidx] = targets[tidx]; + } + nTargets += targets.size();; + + targetList.set(iidx, new hashedWordList(targets)); + } + + reflexes.set(iidx, new scalarList(wallIonFlux.lookup("wallReflexes"))); + + ionFluxBFs.set(iidx, + new surfaceScalarField::GeometricBoundaryField(phi.boundaryField())); +} + +hashedWordList neutrals; + +label nNeutrals = 0; + +for (label nidx = 0; nidx < nTargets; nidx++) +{ + const word n(neutrals_[nidx]); + + if (!neutrals.contains(n)) + { + neutrals.append(n); + nNeutrals += 1; + } +} + +PtrList neutralFluxBFs (neutrals.size()); + +forAll (neutrals, iidx) +{ + neutralFluxBFs.set(iidx, + new surfaceScalarField::GeometricBoundaryField (phi.boundaryField())); +} diff --git a/eReactingFoam.C b/eReactingFoam.C new file mode 100644 index 0000000..e5a6cdd --- /dev/null +++ b/eReactingFoam.C @@ -0,0 +1,120 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011-2013 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +Application + reactingFoam + +Description + Solver for combustion with chemical reactions. + +\*---------------------------------------------------------------------------*/ + +#include "fvCFD.H" +#include "subCycle.H" +#include "turbulenceModel.H" +#include "psiCombustionModel.H" +#include "multivariateScheme.H" +#include "pimpleControl.H" +#include "fvIOoptionList.H" + +#include "CSV.H" + +// #include "bolos.h" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +int main(int argc, char *argv[]) +{ + #include "setRootCase.H" + #include "createTime.H" + #include "createMesh.H" + #include "readGravitationalAcceleration.H" + #include "createFields.H" + #include "createFvOptions.H" + #include "initContinuityErrs.H" + #include "readTimeControls.H" + #include "compressibleCourantNo.H" + #include "setInitialDeltaT.H" + + pimpleControl pimple(mesh); + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + Info<< "\nStarting time loop\n" << endl; + + while (runTime.run()) + { + #include "readTimeControls.H" + #include "compressibleCourantNo.H" + #include "electronCourantNo.H" + #include "setDeltaT.H" + + runTime++; + Info<< "Time = " << runTime.timeName() << nl << endl; + + + #include "rhoEqn.H" + + while (pimple.loop()) + { + #include "numberDensity.H" + #include "PhiEqn.H" + + #include "UEqn.H" + + reaction->correct(); + dQ = reaction->dQ(); + + #include "neControls.H" + #include "neEqnSubCycle.H" + // #include "neEqn.H" + + #include "YEqn.H" + #include "EEqn.H" + + // --- Pressure corrector loop + while (pimple.correct()) + { + #include "pEqn.H" + } + + if (pimple.turbCorr()) + { + turbulence->correct(); + } + } + + runTime.write(); + + Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" + << " ClockTime = " << runTime.elapsedClockTime() << " s" + << nl << endl; + } + + Info<< "End\n" << endl; + + return 0; +} + + +// ************************************************************************* // diff --git a/electronCourantNo.H b/electronCourantNo.H new file mode 100644 index 0000000..44728c0 --- /dev/null +++ b/electronCourantNo.H @@ -0,0 +1,51 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +Global + CourantNo + +Description + Calculates and outputs the mean and maximum Courant Numbers from electron flux. + +\*---------------------------------------------------------------------------*/ + +scalar CoNumVe = 0.0; +scalar meanCoNumVe = 0.0; + +if (mesh.nInternalFaces()) +{ + scalarField sumPhi + ( + fvc::surfaceSum(mag(ve))().internalField() + ); + + CoNumVe = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue(); + + meanCoNumVe = + 0.5*(gSum(sumPhi)/gSum(mesh.V().field()))*runTime.deltaTValue(); +} + +Info<< "Electron Courant Number mean: " << meanCoNumVe + << " max: " << CoNumVe << endl; + +// ************************************************************************* // diff --git a/neControls.H b/neControls.H new file mode 100644 index 0000000..dc1b0b9 --- /dev/null +++ b/neControls.H @@ -0,0 +1,4 @@ +const dictionary& neControls = mesh.solverDict(ne.name()); + +label nNeSubCycles(readLabel(neControls.lookup("nNeSubCycles"))); + diff --git a/neEqn.H b/neEqn.H new file mode 100644 index 0000000..dd1978f --- /dev/null +++ b/neEqn.H @@ -0,0 +1,96 @@ +{ + // 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; + if (relaxDrift < 1.0 && relaxDrift > 0.0) + { + Udrift *= relaxDrift; + } + 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 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); +} diff --git a/neEqnSubCycle.H b/neEqnSubCycle.H new file mode 100644 index 0000000..a3e7bba --- /dev/null +++ b/neEqnSubCycle.H @@ -0,0 +1,22 @@ +if (nNeSubCycles > 1) +{ + dimensionedScalar totalDeltaT = runTime.deltaT(); + + for + ( + subCycle 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" +} diff --git a/numberDensity.H b/numberDensity.H new file mode 100644 index 0000000..c38cfb9 --- /dev/null +++ b/numberDensity.H @@ -0,0 +1,3 @@ +{ + ng = p / T * (NA / R); +} diff --git a/pEqn.H b/pEqn.H new file mode 100644 index 0000000..f7a3004 --- /dev/null +++ b/pEqn.H @@ -0,0 +1,90 @@ +rho = thermo.rho(); + +volScalarField rAU(1.0/UEqn.A()); +surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU)); + +volVectorField HbyA("HbyA", U); +HbyA = rAU*UEqn.H(); + +if (pimple.transonic()) +{ + surfaceScalarField phid + ( + "phid", + fvc::interpolate(psi) + *( + (fvc::interpolate(rho*HbyA) & mesh.Sf()) + + rhorAUf*fvc::ddtCorr(rho, U, phi) + )/fvc::interpolate(rho) + ); + + fvOptions.makeRelative(fvc::interpolate(psi), phid); + + while (pimple.correctNonOrthogonal()) + { + fvScalarMatrix pEqn + ( + fvm::ddt(psi, p) + + fvm::div(phid, p) + - fvm::laplacian(rho*rAU, p) + == + fvOptions(psi, p, rho.name()) + ); + + fvOptions.constrain(pEqn); + + pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter()))); + + if (pimple.finalNonOrthogonalIter()) + { + phi == pEqn.flux(); + } + } +} +else +{ + surfaceScalarField phiHbyA + ( + "phiHbyA", + ( + (fvc::interpolate(rho*HbyA) & mesh.Sf()) + + rhorAUf*fvc::ddtCorr(rho, U, phi) + ) + ); + + fvOptions.makeRelative(fvc::interpolate(rho), phiHbyA); + + while (pimple.correctNonOrthogonal()) + { + fvScalarMatrix pEqn + ( + fvm::ddt(psi, p) + + fvc::div(phiHbyA) + - fvm::laplacian(rho*rAU, p) + == + fvOptions(psi, p, rho.name()) + ); + + fvOptions.constrain(pEqn); + + pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter()))); + + if (pimple.finalNonOrthogonalIter()) + { + phi = phiHbyA + pEqn.flux(); + } + } +} + +#include "rhoEqn.H" +#include "compressibleContinuityErrs.H" + +U = HbyA - rAU*fvc::grad(p); +U.correctBoundaryConditions(); +fvOptions.correct(U); +K = 0.5*magSqr(U); + +if (thermo.dpdt()) +{ + dpdt = fvc::ddt(p); +} diff --git a/readTimeControls.H b/readTimeControls.H new file mode 100644 index 0000000..78522b1 --- /dev/null +++ b/readTimeControls.H @@ -0,0 +1,70 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +Global + readTimeControls + +Description + Read the control parameters used by setDeltaT + +\*---------------------------------------------------------------------------*/ + +const bool adjustTimeStep = + runTime.controlDict().lookupOrDefault("adjustTimeStep", false); + +scalar maxCo = + runTime.controlDict().lookupOrDefault("maxCo", 1.0); + +scalar maxDeltaT = + runTime.controlDict().lookupOrDefault("maxDeltaT", GREAT); + + + + +const dictionary pulseTimeStepDict = runTime.controlDict().subOrEmptyDict("pulseTimeStep"); + +const bool pulseTimeStep = !pulseTimeStepDict.empty(); + +const scalar pulseFrequency = // default 20 kHz + pulseTimeStepDict.lookupOrDefault("pulseFrequency", 20e3); + +const scalar pulseDuration = // default 12 ns + pulseTimeStepDict.lookupOrDefault("pulseDuration", 12e-9); + +const label nStepPulse = // default 1 + ceil(pulseTimeStepDict.lookupOrDefault("nStepPulse", 1)); + +const label nStepIdle = // default 1 + ceil(pulseTimeStepDict.lookupOrDefault("nStepIdle", 1)); + +const scalar pulsePeriod = 1 / pulseFrequency; +const scalar pulseTime = 2*pulseDuration; +const scalar idleTime = pulsePeriod - pulseTime; + +const scalar deltaTpulse = pulseTime / nStepPulse; +const scalar deltaTidle = idleTime / nStepIdle; + +{ +} + +// ************************************************************************* // diff --git a/setDeltaT.H b/setDeltaT.H new file mode 100644 index 0000000..88a215d --- /dev/null +++ b/setDeltaT.H @@ -0,0 +1,85 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +Global + setDeltaT + +Description + Reset the timestep to maintain a constant maximum courant Number. + Reduction of time-step is immediate, but increase is damped to avoid + unstable oscillations. + +\*---------------------------------------------------------------------------*/ + +if (pulseTimeStep) +{ + const scalar t = runTime.timeOutputValue(); + const scalar T = pulsePeriod; + const scalar neg_saw_t = (t - floor(t/T+1.)*T); + + const bool inPulse = neg_saw_t >= -pulseTime; + + if (inPulse) + { + runTime.setDeltaT (deltaTpulse); + Info<< "deltaT to proceed the pulse time " << endl; + } + else + { + const scalar eta = - pulseTime - neg_saw_t; + + if (eta >= deltaTidle) + { + runTime.setDeltaT (deltaTidle); + Info<< "deltaT to proceed the idle time " << endl; + } + else + { + runTime.setDeltaT (deltaTpulse); + Info<< "deltaT to the start of the pulse " << endl; + } + + } + + + Info<< "deltaT = " << runTime.deltaTValue() << endl; +} +else if (adjustTimeStep) +{ + CoNum = max(CoNum, CoNumVe); + scalar maxDeltaTFact = maxCo/(CoNum + SMALL); + scalar deltaTFact = min(min(maxDeltaTFact, 1.0 + 0.1*maxDeltaTFact), 1.2); + + runTime.setDeltaT + ( + min + ( + deltaTFact*runTime.deltaTValue(), + maxDeltaT + ) + ); + + Info<< "deltaT = " << runTime.deltaTValue() << endl; +} + +// ************************************************************************* //