commit 9949acea5edb33c440867133c9ef08debff37b42 Author: ignis Date: Mon May 15 12:26:05 2017 +0900 eReactingFoam copy of plasmaReactingFoam 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; +} + +// ************************************************************************* //