From 67c6ff5c64fb54562af90c43607e97434d783864 Mon Sep 17 00:00:00 2001 From: ignis Date: Tue, 24 Nov 2015 14:26:10 +0900 Subject: [PATCH] plasmaReactingFoam initiated from copy of reactingFoam --- .../combustion/plasmaReactingFoam/Allwmake | 7 ++ .../combustion/plasmaReactingFoam/EEqn.H | 36 +++++++ .../combustion/plasmaReactingFoam/Make/files | 3 + .../plasmaReactingFoam/Make/options | 27 +++++ .../combustion/plasmaReactingFoam/UEqn.H | 21 ++++ .../combustion/plasmaReactingFoam/YEqn.H | 53 +++++++++ .../plasmaReactingFoam/createFields.H | 101 +++++++++++++++++ .../combustion/plasmaReactingFoam/pEqn.H | 90 ++++++++++++++++ .../plasmaReactingFoam/plasmaReactingFoam.C | 102 ++++++++++++++++++ 9 files changed, 440 insertions(+) create mode 100755 applications/solvers/combustion/plasmaReactingFoam/Allwmake create mode 100644 applications/solvers/combustion/plasmaReactingFoam/EEqn.H create mode 100644 applications/solvers/combustion/plasmaReactingFoam/Make/files create mode 100644 applications/solvers/combustion/plasmaReactingFoam/Make/options create mode 100644 applications/solvers/combustion/plasmaReactingFoam/UEqn.H create mode 100644 applications/solvers/combustion/plasmaReactingFoam/YEqn.H create mode 100644 applications/solvers/combustion/plasmaReactingFoam/createFields.H create mode 100644 applications/solvers/combustion/plasmaReactingFoam/pEqn.H create mode 100644 applications/solvers/combustion/plasmaReactingFoam/plasmaReactingFoam.C diff --git a/applications/solvers/combustion/plasmaReactingFoam/Allwmake b/applications/solvers/combustion/plasmaReactingFoam/Allwmake new file mode 100755 index 00000000..1d9987ed --- /dev/null +++ b/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 diff --git a/applications/solvers/combustion/plasmaReactingFoam/EEqn.H b/applications/solvers/combustion/plasmaReactingFoam/EEqn.H new file mode 100644 index 00000000..9267c9a9 --- /dev/null +++ b/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; +} diff --git a/applications/solvers/combustion/plasmaReactingFoam/Make/files b/applications/solvers/combustion/plasmaReactingFoam/Make/files new file mode 100644 index 00000000..dc1e7c43 --- /dev/null +++ b/applications/solvers/combustion/plasmaReactingFoam/Make/files @@ -0,0 +1,3 @@ +plasmaReactingFoam.C + +EXE = $(FOAM_APPBIN)/plasmaReactingFoam diff --git a/applications/solvers/combustion/plasmaReactingFoam/Make/options b/applications/solvers/combustion/plasmaReactingFoam/Make/options new file mode 100644 index 00000000..15cf68cb --- /dev/null +++ b/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 diff --git a/applications/solvers/combustion/plasmaReactingFoam/UEqn.H b/applications/solvers/combustion/plasmaReactingFoam/UEqn.H new file mode 100644 index 00000000..60c2a8a7 --- /dev/null +++ b/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); + } diff --git a/applications/solvers/combustion/plasmaReactingFoam/YEqn.H b/applications/solvers/combustion/plasmaReactingFoam/YEqn.H new file mode 100644 index 00000000..c8d34cf1 --- /dev/null +++ b/applications/solvers/combustion/plasmaReactingFoam/YEqn.H @@ -0,0 +1,53 @@ +tmp > mvConvection +( + fv::convectionScheme::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); +} diff --git a/applications/solvers/combustion/plasmaReactingFoam/createFields.H b/applications/solvers/combustion/plasmaReactingFoam/createFields.H new file mode 100644 index 00000000..32bdd372 --- /dev/null +++ b/applications/solvers/combustion/plasmaReactingFoam/createFields.H @@ -0,0 +1,101 @@ +Info<< "Creating reaction model\n" << endl; + +autoPtr reaction +( + combustionModels::psiCombustionModel::New(mesh) +); + +psiReactionThermo& thermo = reaction->thermo(); +thermo.validate(args.executable(), "h", "e"); + +basicMultiComponentMixture& composition = thermo.composition(); +PtrList& 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 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) +); diff --git a/applications/solvers/combustion/plasmaReactingFoam/pEqn.H b/applications/solvers/combustion/plasmaReactingFoam/pEqn.H new file mode 100644 index 00000000..f7a3004a --- /dev/null +++ b/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); +} diff --git a/applications/solvers/combustion/plasmaReactingFoam/plasmaReactingFoam.C b/applications/solvers/combustion/plasmaReactingFoam/plasmaReactingFoam.C new file mode 100644 index 00000000..831df5c9 --- /dev/null +++ b/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 . + +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; +} + + +// ************************************************************************* //