commit cd318a2e1b19e390b11bd19caf6a08b5da597446 Author: ignis Date: Thu Apr 19 18:46:28 2018 +0900 Initial commit (reactingFoam) 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; +} + + +// ************************************************************************* //