From cd318a2e1b19e390b11bd19caf6a08b5da597446 Mon Sep 17 00:00:00 2001 From: ignis Date: Thu, 19 Apr 2018 18:46:28 +0900 Subject: [PATCH] Initial commit (reactingFoam) --- EEqn.H | 36 +++++++++++++ Make/files | 3 ++ Make/options | 28 ++++++++++ UEqn.H | 25 +++++++++ YEqn.H | 53 +++++++++++++++++++ createFieldRefs.H | 2 + createFields.H | 125 ++++++++++++++++++++++++++++++++++++++++++++ eReactingFoam.C | 128 ++++++++++++++++++++++++++++++++++++++++++++++ pEqn.H | 106 ++++++++++++++++++++++++++++++++++++++ pcEqn.H | 127 +++++++++++++++++++++++++++++++++++++++++++++ setRDeltaT.H | 128 ++++++++++++++++++++++++++++++++++++++++++++++ 11 files changed, 761 insertions(+) create mode 100644 EEqn.H create mode 100644 Make/files create mode 100644 Make/options create mode 100644 UEqn.H create mode 100644 YEqn.H create mode 100644 createFieldRefs.H create mode 100644 createFields.H create mode 100644 eReactingFoam.C create mode 100644 pEqn.H create mode 100644 pcEqn.H create mode 100644 setRDeltaT.H diff --git a/EEqn.H b/EEqn.H new file mode 100644 index 0000000..9267c9a --- /dev/null +++ b/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/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..6363f6d --- /dev/null +++ b/Make/options @@ -0,0 +1,28 @@ +EXE_INC = \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude \ + -I$(LIB_SRC)/sampling/lnInclude \ + -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \ + -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \ + -I$(LIB_SRC)/transportModels/compressible/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 \ + -lturbulenceModels \ + -lcompressibleTurbulenceModels \ + -lreactionThermophysicalModels \ + -lspecie \ + -lcompressibleTransportModels \ + -lfluidThermophysicalModels \ + -lchemistryModel \ + -lODE \ + -lcombustionModels diff --git a/UEqn.H b/UEqn.H new file mode 100644 index 0000000..e4caa1c --- /dev/null +++ b/UEqn.H @@ -0,0 +1,25 @@ +// Solve the Momentum equation + +MRF.correctBoundaryVelocity(U); + +tmp tUEqn +( + fvm::ddt(rho, U) + fvm::div(phi, U) + + MRF.DDt(rho, U) + + turbulence->divDevRhoReff(U) + == + fvOptions(rho, U) +); +fvVectorMatrix& UEqn = tUEqn.ref(); + +UEqn.relax(); + +fvOptions.constrain(UEqn); + +if (pimple.momentumPredictor()) +{ + solve(UEqn == -fvc::grad(p)); + + fvOptions.correct(U); + K = 0.5*magSqr(U); +} diff --git a/YEqn.H b/YEqn.H new file mode 100644 index 0000000..8cb5664 --- /dev/null +++ b/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/createFieldRefs.H b/createFieldRefs.H new file mode 100644 index 0000000..6e1184a --- /dev/null +++ b/createFieldRefs.H @@ -0,0 +1,2 @@ +const volScalarField& psi = thermo.psi(); +const volScalarField& T = thermo.T(); diff --git a/createFields.H b/createFields.H new file mode 100644 index 0000000..da22c1a --- /dev/null +++ b/createFields.H @@ -0,0 +1,125 @@ +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(); + +#include "compressibleCreatePhi.H" + +dimensionedScalar rhoMax +( + dimensionedScalar::lookupOrDefault + ( + "rhoMax", + pimple.dict(), + dimDensity, + GREAT + ) +); + +dimensionedScalar rhoMin +( + dimensionedScalar::lookupOrDefault + ( + "rhoMin", + pimple.dict(), + dimDensity, + 0 + ) +); + +mesh.setFluxRequired(p.name()); + +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) +); + +#include "createMRF.H" diff --git a/eReactingFoam.C b/eReactingFoam.C new file mode 100644 index 0000000..efcc5eb --- /dev/null +++ b/eReactingFoam.C @@ -0,0 +1,128 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2011-2016 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 "turbulentFluidThermoModel.H" +#include "psiCombustionModel.H" +#include "multivariateScheme.H" +#include "pimpleControl.H" +#include "fvOptions.H" +#include "localEulerDdtScheme.H" +#include "fvcSmooth.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +int main(int argc, char *argv[]) +{ + #include "postProcess.H" + + #include "setRootCase.H" + #include "createTime.H" + #include "createMesh.H" + #include "createControl.H" + #include "createTimeControls.H" + #include "createRDeltaT.H" + #include "initContinuityErrs.H" + #include "createFields.H" + #include "createFieldRefs.H" + #include "createFvOptions.H" + + turbulence->validate(); + + if (!LTS) + { + #include "compressibleCourantNo.H" + #include "setInitialDeltaT.H" + } + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + Info<< "\nStarting time loop\n" << endl; + + while (runTime.run()) + { + #include "readTimeControls.H" + + if (LTS) + { + #include "setRDeltaT.H" + } + else + { + #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()) + { + if (pimple.consistent()) + { + #include "pcEqn.H" + } + else + { + #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/pEqn.H b/pEqn.H new file mode 100644 index 0000000..ac7107a --- /dev/null +++ b/pEqn.H @@ -0,0 +1,106 @@ +rho = thermo.rho(); +rho = max(rho, rhoMin); +rho = min(rho, rhoMax); +rho.relax(); + +volScalarField rAU(1.0/UEqn.A()); +surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU)); +volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p)); + +if (pimple.nCorrPISO() <= 1) +{ + tUEqn.clear(); +} + +if (pimple.transonic()) +{ + surfaceScalarField phid + ( + "phid", + fvc::interpolate(psi) + *( + fvc::flux(HbyA) + + rhorAUf*fvc::ddtCorr(rho, U, phi)/fvc::interpolate(rho) + ) + ); + + MRF.makeRelative(fvc::interpolate(psi), phid); + + while (pimple.correctNonOrthogonal()) + { + fvScalarMatrix pEqn + ( + fvm::ddt(psi, p) + + fvm::div(phid, p) + - fvm::laplacian(rhorAUf, p) + == + fvOptions(psi, p, rho.name()) + ); + + pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter()))); + + if (pimple.finalNonOrthogonalIter()) + { + phi == pEqn.flux(); + } + } +} +else +{ + surfaceScalarField phiHbyA + ( + "phiHbyA", + ( + fvc::flux(rho*HbyA) + + rhorAUf*fvc::ddtCorr(rho, U, phi) + ) + ); + + MRF.makeRelative(fvc::interpolate(rho), phiHbyA); + + // Update the pressure BCs to ensure flux consistency + constrainPressure(p, rho, U, phiHbyA, rhorAUf, MRF); + + while (pimple.correctNonOrthogonal()) + { + fvScalarMatrix pEqn + ( + fvm::ddt(psi, p) + + fvc::div(phiHbyA) + - fvm::laplacian(rhorAUf, p) + == + fvOptions(psi, p, rho.name()) + ); + + pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter()))); + + if (pimple.finalNonOrthogonalIter()) + { + phi = phiHbyA + pEqn.flux(); + } + } +} + +#include "rhoEqn.H" +#include "compressibleContinuityErrs.H" + +// Explicitly relax pressure for momentum corrector +p.relax(); + +// Recalculate density from the relaxed pressure +rho = thermo.rho(); +rho = max(rho, rhoMin); +rho = min(rho, rhoMax); +rho.relax(); +Info<< "rho max/min : " << max(rho).value() + << " " << min(rho).value() << endl; + +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/pcEqn.H b/pcEqn.H new file mode 100644 index 0000000..713f443 --- /dev/null +++ b/pcEqn.H @@ -0,0 +1,127 @@ +rho = thermo.rho(); +rho = max(rho, rhoMin); +rho = min(rho, rhoMax); +rho.relax(); + +volScalarField rAU(1.0/UEqn.A()); +volScalarField rAtU(1.0/(1.0/rAU - UEqn.H1())); +volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p)); + +if (pimple.nCorrPISO() <= 1) +{ + tUEqn.clear(); +} + +if (pimple.transonic()) +{ + surfaceScalarField phid + ( + "phid", + fvc::interpolate(psi) + *( + fvc::flux(HbyA) + + fvc::interpolate(rho*rAU)*fvc::ddtCorr(rho, U, phi) + /fvc::interpolate(rho) + ) + ); + + MRF.makeRelative(fvc::interpolate(psi), phid); + + surfaceScalarField phic + ( + "phic", + fvc::interpolate(rho*(rAtU - rAU))*fvc::snGrad(p)*mesh.magSf() + ); + + HbyA -= (rAU - rAtU)*fvc::grad(p); + + volScalarField rhorAtU("rhorAtU", rho*rAtU); + + while (pimple.correctNonOrthogonal()) + { + fvScalarMatrix pEqn + ( + fvm::ddt(psi, p) + + fvm::div(phid, p) + + fvc::div(phic) + - fvm::laplacian(rhorAtU, p) + == + fvOptions(psi, p, rho.name()) + ); + + pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter()))); + + if (pimple.finalNonOrthogonalIter()) + { + phi == phic + pEqn.flux(); + } + } +} +else +{ + surfaceScalarField phiHbyA + ( + "phiHbyA", + ( + fvc::flux(rho*HbyA) + + fvc::interpolate(rho*rAU)*fvc::ddtCorr(rho, U, phi) + ) + ); + + MRF.makeRelative(fvc::interpolate(rho), phiHbyA); + + phiHbyA += fvc::interpolate(rho*(rAtU - rAU))*fvc::snGrad(p)*mesh.magSf(); + HbyA -= (rAU - rAtU)*fvc::grad(p); + + volScalarField rhorAtU("rhorAtU", rho*rAtU); + + // Update the pressure BCs to ensure flux consistency + constrainPressure(p, rho, U, phiHbyA, rhorAtU, MRF); + + while (pimple.correctNonOrthogonal()) + { + fvScalarMatrix pEqn + ( + fvm::ddt(psi, p) + + fvc::div(phiHbyA) + - fvm::laplacian(rhorAtU, p) + == + fvOptions(psi, p, rho.name()) + ); + + pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter()))); + + if (pimple.finalNonOrthogonalIter()) + { + phi = phiHbyA + pEqn.flux(); + } + } +} + +#include "rhoEqn.H" +#include "compressibleContinuityErrs.H" + +// Explicitly relax pressure for momentum corrector +p.relax(); + +U = HbyA - rAtU*fvc::grad(p); +U.correctBoundaryConditions(); +fvOptions.correct(U); +K = 0.5*magSqr(U); + +if (thermo.dpdt()) +{ + dpdt = fvc::ddt(p); +} + +// Recalculate density from the relaxed pressure +rho = thermo.rho(); +rho = max(rho, rhoMin); +rho = min(rho, rhoMax); + +if (!pimple.transonic()) +{ + rho.relax(); +} + +Info<< "rho max/min : " << max(rho).value() << " " << min(rho).value() << endl; diff --git a/setRDeltaT.H b/setRDeltaT.H new file mode 100644 index 0000000..5c24294 --- /dev/null +++ b/setRDeltaT.H @@ -0,0 +1,128 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2013-2016 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 . + +\*---------------------------------------------------------------------------*/ + +{ + volScalarField& rDeltaT = trDeltaT.ref(); + + const dictionary& pimpleDict = pimple.dict(); + + // Maximum flow Courant number + scalar maxCo(readScalar(pimpleDict.lookup("maxCo"))); + + // Maximum time scale + scalar maxDeltaT(pimpleDict.lookupOrDefault("maxDeltaT", GREAT)); + + // Smoothing parameter (0-1) when smoothing iterations > 0 + scalar rDeltaTSmoothingCoeff + ( + pimpleDict.lookupOrDefault("rDeltaTSmoothingCoeff", 0.1) + ); + + // Damping coefficient (1-0) + scalar rDeltaTDampingCoeff + ( + pimpleDict.lookupOrDefault("rDeltaTDampingCoeff", 1) + ); + + // Maximum change in cell temperature per iteration + // (relative to previous value) + scalar alphaTemp(pimpleDict.lookupOrDefault("alphaTemp", 0.05)); + + + Info<< "Time scales min/max:" << endl; + + // Cache old reciprocal time scale field + volScalarField rDeltaT0("rDeltaT0", rDeltaT); + + // Flow time scale + { + rDeltaT.ref() = + ( + fvc::surfaceSum(mag(phi))()() + /((2*maxCo)*mesh.V()*rho()) + ); + + // Limit the largest time scale + rDeltaT.max(1/maxDeltaT); + + Info<< " Flow = " + << gMin(1/rDeltaT.primitiveField()) << ", " + << gMax(1/rDeltaT.primitiveField()) << endl; + } + + // Reaction source time scale + if (alphaTemp < 1.0) + { + volScalarField::Internal rDeltaTT + ( + mag(reaction->Sh())/(alphaTemp*rho*thermo.Cp()*T) + ); + + Info<< " Temperature = " + << gMin(1/(rDeltaTT.field() + VSMALL)) << ", " + << gMax(1/(rDeltaTT.field() + VSMALL)) << endl; + + rDeltaT.ref() = max + ( + rDeltaT(), + rDeltaTT + ); + } + + // Update tho boundary values of the reciprocal time-step + rDeltaT.correctBoundaryConditions(); + + // Spatially smooth the time scale field + if (rDeltaTSmoothingCoeff < 1.0) + { + fvc::smooth(rDeltaT, rDeltaTSmoothingCoeff); + } + + // Limit rate of change of time scale + // - reduce as much as required + // - only increase at a fraction of old time scale + if + ( + rDeltaTDampingCoeff < 1.0 + && runTime.timeIndex() > runTime.startTimeIndex() + 1 + ) + { + rDeltaT = max + ( + rDeltaT, + (scalar(1.0) - rDeltaTDampingCoeff)*rDeltaT0 + ); + } + + // Update tho boundary values of the reciprocal time-step + rDeltaT.correctBoundaryConditions(); + + Info<< " Overall = " + << gMin(1/rDeltaT.primitiveField()) + << ", " << gMax(1/rDeltaT.primitiveField()) << endl; +} + + +// ************************************************************************* //