eReactingFoam copy of plasmaReactingFoam

Allwmake

@@ -0,0 +1,7 @@
+#!/bin/sh
+cd ${0%/*} || exit 1    # run from this directory
+set -x
+
+wmake
+
+# ----------------------------------------------------------------- end-of-file

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;
+}

Make/files

@@ -0,0 +1,3 @@
+eReactingFoam.C
+
+EXE = $(FOAM_USER_APPBIN)/eReactingFoam

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

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<volScalarField> tMagE (mag(E));
+    const volScalarField &magE = tMagE();
+
+    magE.writeMinMax(Info);
+
+    En = magE / (ng);
+}

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();

YEqn.H

@@ -0,0 +1,181 @@
+/*
+tmp<fv::convectionScheme<scalar> > mvConvection
+(
+    fv::convectionScheme<scalar>::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);
+}

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<scalar>::New("mueN", physicalProperties));
+const DataEntry<scalar> &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<scalar>::New("DeN", physicalProperties));
+const DataEntry<scalar> &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<scalar>::New ("TeOfEn", physicalProperties));
+const DataEntry<scalar> &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<combustionModels::psiCombustionModel> 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<volScalarField>& Y = composition.Y();
+const PtrList<volScalarField>& 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<compressible::turbulenceModel> 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<scalar>::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<surfaceScalarField> 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<hashedWordList> targetList (ions.size());
+PtrList<scalarList> reflexes (ions.size());
+PtrList<surfaceScalarField::GeometricBoundaryField> 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<surfaceScalarField::GeometricBoundaryField> neutralFluxBFs (neutrals.size());
+
+forAll (neutrals, iidx)
+{
+    neutralFluxBFs.set(iidx,
+            new surfaceScalarField::GeometricBoundaryField (phi.boundaryField()));
+}

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 <http://www.gnu.org/licenses/>.
+
+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;
+}
+
+
+// ************************************************************************* //

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 <http://www.gnu.org/licenses/>.
+
+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;
+
+// ************************************************************************* //

neControls.H

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

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<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);
+}

neEqnSubCycle.H

@@ -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"
+}

numberDensity.H

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

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);
+}

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 <http://www.gnu.org/licenses/>.
+
+Global
+    readTimeControls
+
+Description
+    Read the control parameters used by setDeltaT
+
+\*---------------------------------------------------------------------------*/
+
+const bool adjustTimeStep =
+    runTime.controlDict().lookupOrDefault("adjustTimeStep", false);
+
+scalar maxCo =
+    runTime.controlDict().lookupOrDefault<scalar>("maxCo", 1.0);
+
+scalar maxDeltaT =
+    runTime.controlDict().lookupOrDefault<scalar>("maxDeltaT", GREAT);
+
+
+
+
+const dictionary pulseTimeStepDict = runTime.controlDict().subOrEmptyDict("pulseTimeStep");
+
+const bool pulseTimeStep = !pulseTimeStepDict.empty();
+
+const scalar pulseFrequency = // default 20 kHz
+    pulseTimeStepDict.lookupOrDefault<scalar>("pulseFrequency", 20e3);
+
+const scalar pulseDuration =  // default 12 ns
+    pulseTimeStepDict.lookupOrDefault<scalar>("pulseDuration", 12e-9);
+
+const label nStepPulse =  // default 1
+    ceil(pulseTimeStepDict.lookupOrDefault<scalar>("nStepPulse", 1));
+
+const label nStepIdle =  // default 1
+    ceil(pulseTimeStepDict.lookupOrDefault<scalar>("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;
+
+{
+}
+
+// ************************************************************************* //

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 <http://www.gnu.org/licenses/>.
+
+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;
+}
+
+// ************************************************************************* //