OpenFOAM-2.4.x/applications/solvers/combustion/plasmaReactingFoam/createFields.H

Info<< "Reading physicalProperties\n" << endl;

IOdictionary physicalProperties
(
    IOobject
    (
        "physicalProperties",
        runTime.constant(),
        mesh,
        IOobject::MUST_READ_IF_MODIFIED,
        IOobject::NO_WRITE
    )
);

dimensionedScalar epsilon0
(
    physicalProperties.lookup("epsilon0")
);

// Convert E/n in SI unit to table unit.
// Default V m^2 => Td
scalar EnFac (
        physicalProperties.lookupOrDefault("EnFac", 1.0e21)
        );

CSV<scalar> mueN("mueN", physicalProperties);
// Convert mu_e * n_g value from table into SI unit.
scalar mueNFac (
        physicalProperties.lookupOrDefault("mueNFac", 1.0)
        );

CSV<scalar> DeN("DeN", physicalProperties);
// Convert D_e * n_g value from table into SI unit.
scalar DeNFac (
        physicalProperties.lookupOrDefault("DeNFac", 1.0)
        );

CSV<scalar> TeOfEn("TeOfEn", physicalProperties);
// 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("E", -fvc::grad(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;

//- Pi
const scalar pi = constant::mathematical::pi;

PtrList<volScalarField>& Y = composition.Y();

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"

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

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;
surfaceScalarField ve
(
    IOobject
    (
        "ve",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::NO_WRITE
    ),
    - linearInterpolate(mue*E/ng) & mesh.Sf()
);