plasmaReactingFoam initiated from copy of reactingFoam

applications/solvers/combustion/plasmaReactingFoam/Allwmake

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

applications/solvers/combustion/plasmaReactingFoam/EEqn.H

@@ -0,0 +1,36 @@
+{
+    volScalarField& he = thermo.he();
+
+    fvScalarMatrix EEqn
+    (
+        fvm::ddt(rho, he) + mvConvection->fvmDiv(phi, he)
+      + 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;
+}

applications/solvers/combustion/plasmaReactingFoam/Make/files

@@ -0,0 +1,3 @@
+plasmaReactingFoam.C
+
+EXE = $(FOAM_APPBIN)/plasmaReactingFoam

applications/solvers/combustion/plasmaReactingFoam/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

applications/solvers/combustion/plasmaReactingFoam/UEqn.H

@@ -0,0 +1,21 @@
+    fvVectorMatrix UEqn
+    (
+        fvm::ddt(rho, U)
+      + fvm::div(phi, U)
+      + turbulence->divDevRhoReff(U)
+     ==
+        rho*g
+      + fvOptions(rho, U)
+    );
+
+    UEqn.relax();
+
+    fvOptions.constrain(UEqn);
+
+    if (pimple.momentumPredictor())
+    {
+        solve(UEqn == -fvc::grad(p));
+
+        fvOptions.correct(U);
+        K = 0.5*magSqr(U);
+    }

applications/solvers/combustion/plasmaReactingFoam/YEqn.H

@@ -0,0 +1,53 @@
+tmp<fv::convectionScheme<scalar> > mvConvection
+(
+    fv::convectionScheme<scalar>::New
+    (
+        mesh,
+        fields,
+        phi,
+        mesh.divScheme("div(phi,Yi_h)")
+    )
+);
+
+{
+    reaction->correct();
+    dQ = reaction->dQ();
+    label inertIndex = -1;
+    volScalarField Yt(0.0*Y[0]);
+
+    forAll(Y, i)
+    {
+        if (Y[i].name() != inertSpecie)
+        {
+            volScalarField& Yi = Y[i];
+
+            fvScalarMatrix YiEqn
+            (
+                fvm::ddt(rho, Yi)
+              + mvConvection->fvmDiv(phi, Yi)
+              - fvm::laplacian(turbulence->muEff(), 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);
+}

applications/solvers/combustion/plasmaReactingFoam/createFields.H

@@ -0,0 +1,101 @@
+Info<< "Creating reaction model\n" << endl;
+
+autoPtr<combustionModels::psiCombustionModel> reaction
+(
+    combustionModels::psiCombustionModel::New(mesh)
+);
+
+psiReactionThermo& thermo = reaction->thermo();
+thermo.validate(args.executable(), "h", "e");
+
+basicMultiComponentMixture& composition = thermo.composition();
+PtrList<volScalarField>& Y = composition.Y();
+
+word inertSpecie(thermo.lookup("inertSpecie"));
+
+volScalarField rho
+(
+    IOobject
+    (
+        "rho",
+        runTime.timeName(),
+        mesh
+    ),
+    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)
+);

applications/solvers/combustion/plasmaReactingFoam/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);
+}

applications/solvers/combustion/plasmaReactingFoam/plasmaReactingFoam.C

@@ -0,0 +1,102 @@
+/*---------------------------------------------------------------------------*\
+  =========                 |
+  \\      /  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 "turbulenceModel.H"
+#include "psiCombustionModel.H"
+#include "multivariateScheme.H"
+#include "pimpleControl.H"
+#include "fvIOoptionList.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 "setDeltaT.H"
+
+        runTime++;
+        Info<< "Time = " << runTime.timeName() << nl << endl;
+
+        #include "rhoEqn.H"
+
+        while (pimple.loop())
+        {
+            #include "UEqn.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;
+}
+
+
+// ************************************************************************* //