From 98a77f44ca2aabb54ed9f84d74cd415b979e24f2 Mon Sep 17 00:00:00 2001 From: andy Date: Thu, 8 Apr 2010 16:56:50 +0100 Subject: [PATCH] ENH: Added new reactingParcelFilmFoam solver - Takes reactingParcelFoam solver and adds surface film capability --- .../reactingParcelFilmFoam/Make/files | 3 + .../reactingParcelFilmFoam/Make/options | 38 ++++++ .../lagrangian/reactingParcelFilmFoam/UEqn.H | 16 +++ .../lagrangian/reactingParcelFilmFoam/YEqn.H | 44 ++++++ .../reactingParcelFilmFoam/chemistry.H | 27 ++++ .../reactingParcelFilmFoam/createClouds.H | 9 ++ .../reactingParcelFilmFoam/createFields.H | 122 +++++++++++++++++ .../createSurfaceFilmModel.H | 6 + .../lagrangian/reactingParcelFilmFoam/hsEqn.H | 21 +++ .../lagrangian/reactingParcelFilmFoam/pEqn.H | 72 ++++++++++ .../reactingParcelFilmFoam.C | 127 ++++++++++++++++++ .../readChemistryProperties.H | 23 ++++ .../reactingParcelFilmFoam/rhoEqn.H | 43 ++++++ 13 files changed, 551 insertions(+) create mode 100644 applications/solvers/lagrangian/reactingParcelFilmFoam/Make/files create mode 100644 applications/solvers/lagrangian/reactingParcelFilmFoam/Make/options create mode 100644 applications/solvers/lagrangian/reactingParcelFilmFoam/UEqn.H create mode 100644 applications/solvers/lagrangian/reactingParcelFilmFoam/YEqn.H create mode 100644 applications/solvers/lagrangian/reactingParcelFilmFoam/chemistry.H create mode 100644 applications/solvers/lagrangian/reactingParcelFilmFoam/createClouds.H create mode 100644 applications/solvers/lagrangian/reactingParcelFilmFoam/createFields.H create mode 100644 applications/solvers/lagrangian/reactingParcelFilmFoam/createSurfaceFilmModel.H create mode 100644 applications/solvers/lagrangian/reactingParcelFilmFoam/hsEqn.H create mode 100644 applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H create mode 100644 applications/solvers/lagrangian/reactingParcelFilmFoam/reactingParcelFilmFoam.C create mode 100644 applications/solvers/lagrangian/reactingParcelFilmFoam/readChemistryProperties.H create mode 100644 applications/solvers/lagrangian/reactingParcelFilmFoam/rhoEqn.H diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/Make/files b/applications/solvers/lagrangian/reactingParcelFilmFoam/Make/files new file mode 100644 index 0000000000..a76d2acc65 --- /dev/null +++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/Make/files @@ -0,0 +1,3 @@ +reactingParcelFilmFoam.C + +EXE = $(FOAM_APPBIN)/reactingParcelFilmFoam diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/Make/options b/applications/solvers/lagrangian/reactingParcelFilmFoam/Make/options new file mode 100644 index 0000000000..8b15c3e9ae --- /dev/null +++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/Make/options @@ -0,0 +1,38 @@ +DEV_PATH=./../.. + +EXE_INC = \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I${LIB_SRC}/meshTools/lnInclude \ + -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \ + -I$(LIB_SRC)/thermophysicalModels/pdfs/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/solids/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/solidMixture/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/liquids/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/liquidMixture/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/thermophysicalFunctions/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \ + -I$(LIB_SRC)/surfaceFilmModels/lnInclude \ + -I$(LIB_SRC)/lagrangian/basic/lnInclude \ + -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \ + -I$(LIB_SRC)/ODE/lnInclude + +EXE_LIBS = \ + -lfiniteVolume \ + -lmeshTools \ + -lcompressibleRASModels \ + -lcompressibleLESModels \ + -lspecie \ + -lbasicThermophysicalModels \ + -lsolids \ + -lsolidMixture \ + -lthermophysicalFunctions \ + -lreactionThermophysicalModels \ + -lchemistryModel \ + -lradiation \ + -lsurfaceFilmModels \ + -llagrangianIntermediate \ + -lODE diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/UEqn.H b/applications/solvers/lagrangian/reactingParcelFilmFoam/UEqn.H new file mode 100644 index 0000000000..3c4a927091 --- /dev/null +++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/UEqn.H @@ -0,0 +1,16 @@ + fvVectorMatrix UEqn + ( + fvm::ddt(rho, U) + + fvm::div(phi, U) + + turbulence->divDevRhoReff(U) + == + rho.dimensionedInternalField()*g + + parcels.SU() + ); + + UEqn.relax(); + + if (momentumPredictor) + { + solve(UEqn == -fvc::grad(p)); + } diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/YEqn.H b/applications/solvers/lagrangian/reactingParcelFilmFoam/YEqn.H new file mode 100644 index 0000000000..5cf44cb09f --- /dev/null +++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/YEqn.H @@ -0,0 +1,44 @@ +tmp > mvConvection +( + fv::convectionScheme::New + ( + mesh, + fields, + phi, + mesh.divScheme("div(phi,Yi_h)") + ) +); + + +{ + label inertIndex = -1; + volScalarField Yt = 0.0*Y[0]; + + for (label i=0; ifvmDiv(phi, Yi) + - fvm::laplacian(turbulence->muEff(), Yi) + == + parcels.Srho(i) + + kappa*chemistry.RR(i)().dimensionedInternalField(), + mesh.solver("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/lagrangian/reactingParcelFilmFoam/chemistry.H b/applications/solvers/lagrangian/reactingParcelFilmFoam/chemistry.H new file mode 100644 index 0000000000..12ed815d1d --- /dev/null +++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/chemistry.H @@ -0,0 +1,27 @@ +{ + Info << "Solving chemistry" << endl; + + chemistry.solve + ( + runTime.value() - runTime.deltaTValue(), + runTime.deltaTValue() + ); + + // turbulent time scale + if (turbulentReaction) + { + DimensionedField tk = + Cmix*sqrt(turbulence->muEff()/rho/turbulence->epsilon()); + DimensionedField tc = + chemistry.tc()().dimensionedInternalField(); + + // Chalmers PaSR model + kappa = (runTime.deltaT() + tc)/(runTime.deltaT() + tc + tk); + } + else + { + kappa = 1.0; + } + + chemistrySh = kappa*chemistry.Sh()(); +} diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/createClouds.H b/applications/solvers/lagrangian/reactingParcelFilmFoam/createClouds.H new file mode 100644 index 0000000000..f3e043143c --- /dev/null +++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/createClouds.H @@ -0,0 +1,9 @@ +Info<< "\nConstructing reacting cloud" << endl; +thermoReactingCloud parcels +( + "reactingCloud1", + rho, + U, + g, + thermo +); diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/createFields.H b/applications/solvers/lagrangian/reactingParcelFilmFoam/createFields.H new file mode 100644 index 0000000000..f7c2caf8c6 --- /dev/null +++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/createFields.H @@ -0,0 +1,122 @@ + Info<< "Reading thermophysical properties\n" << endl; + + autoPtr pChemistry + ( + psiChemistryModel::New(mesh) + ); + psiChemistryModel& chemistry = pChemistry(); + + hsCombustionThermo& thermo = chemistry.thermo(); + + basicMultiComponentMixture& composition = thermo.composition(); + PtrList& Y = composition.Y(); + + word inertSpecie(thermo.lookup("inertSpecie")); + + volScalarField& p = thermo.p(); + volScalarField& hs = thermo.hs(); + const volScalarField& T = thermo.T(); + const volScalarField& psi = thermo.psi(); + + Info<< "Creating field rho\n" << endl; + volScalarField rho + ( + IOobject + ( + "rho", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + thermo.rho() + ); + + Info<< "\nReading field U\n" << endl; + volVectorField U + ( + IOobject + ( + "U", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh + ); + + #include "compressibleCreatePhi.H" + + Info<< "Creating field kappa\n" << endl; + DimensionedField kappa + ( + IOobject + ( + "kappa", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + mesh, + dimensionedScalar("zero", dimless, 0.0) + ); + + Info<< "Creating turbulence model\n" << endl; + autoPtr turbulence + ( + compressible::turbulenceModel::New + ( + rho, + U, + phi, + thermo + ) + ); + + Info<< "Creating field DpDt\n" << endl; + volScalarField DpDt + ( + "DpDt", + fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p) + ); + + multivariateSurfaceInterpolationScheme::fieldTable fields; + + forAll (Y, i) + { + fields.add(Y[i]); + } + fields.add(hs); + + IOdictionary additionalControlsDict + ( + IOobject + ( + "additionalControls", + runTime.constant(), + mesh, + IOobject::MUST_READ, + IOobject::NO_WRITE + ) + ); + + Switch solvePrimaryRegion + ( + additionalControlsDict.lookup("solvePrimaryRegion") + ); + + DimensionedField chemistrySh + ( + IOobject + ( + "chemistry::Sh", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh, + dimensionedScalar("chemistrySh", dimEnergy/dimTime/dimVolume, 0.0) + ); diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/createSurfaceFilmModel.H b/applications/solvers/lagrangian/reactingParcelFilmFoam/createSurfaceFilmModel.H new file mode 100644 index 0000000000..15046e3d49 --- /dev/null +++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/createSurfaceFilmModel.H @@ -0,0 +1,6 @@ +Info<< "\nConstructing surface film model" << endl; + +autoPtr + tsurfaceFilm(surfaceFilmModels::surfaceFilmModel::New(mesh, g)); +surfaceFilmModels::surfaceFilmModel& surfaceFilm = tsurfaceFilm(); + diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/hsEqn.H b/applications/solvers/lagrangian/reactingParcelFilmFoam/hsEqn.H new file mode 100644 index 0000000000..156d5505ab --- /dev/null +++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/hsEqn.H @@ -0,0 +1,21 @@ +{ + fvScalarMatrix hsEqn + ( + fvm::ddt(rho, hs) + + mvConvection->fvmDiv(phi, hs) + - fvm::laplacian(turbulence->alphaEff(), hs) + == + DpDt + + parcels.Sh() + + radiation->Shs(thermo) + + chemistrySh + ); + + hsEqn.relax(); + + hsEqn.solve(); + + thermo.correct(); + + radiation->correct(); +} diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H b/applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H new file mode 100644 index 0000000000..d27afcce8a --- /dev/null +++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/pEqn.H @@ -0,0 +1,72 @@ +rho = thermo.rho(); + +volScalarField rUA = 1.0/UEqn.A(); +U = rUA*UEqn.H(); + +if (transonic) +{ + surfaceScalarField phid + ( + "phid", + fvc::interpolate(psi) + *( + (fvc::interpolate(U) & mesh.Sf()) + + fvc::ddtPhiCorr(rUA, rho, U, phi) + ) + ); + + for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) + { + fvScalarMatrix pEqn + ( + fvm::ddt(psi, p) + + fvm::div(phid, p) + - fvm::laplacian(rho*rUA, p) + == + parcels.Srho() + ); + + pEqn.solve(); + + if (nonOrth == nNonOrthCorr) + { + phi == pEqn.flux(); + } + } +} +else +{ + phi = + fvc::interpolate(rho) + *( + (fvc::interpolate(U) & mesh.Sf()) + + fvc::ddtPhiCorr(rUA, rho, U, phi) + ); + + for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) + { + fvScalarMatrix pEqn + ( + fvm::ddt(psi, p) + + fvc::div(phi) + - fvm::laplacian(rho*rUA, p) + == + parcels.Srho() + ); + + pEqn.solve(); + + if (nonOrth == nNonOrthCorr) + { + phi += pEqn.flux(); + } + } +} + +#include "rhoEqn.H" +#include "compressibleContinuityErrs.H" + +U -= rUA*fvc::grad(p); +U.correctBoundaryConditions(); + +DpDt = fvc::DDt(surfaceScalarField("phiU", phi/fvc::interpolate(rho)), p); diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/reactingParcelFilmFoam.C b/applications/solvers/lagrangian/reactingParcelFilmFoam/reactingParcelFilmFoam.C new file mode 100644 index 0000000000..f39d8f070a --- /dev/null +++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/reactingParcelFilmFoam.C @@ -0,0 +1,127 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2009-2010 OpenCFD Ltd. + \\/ 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 + reactingParcelFilmFoam + +Description + Transient PISO solver for compressible, laminar or turbulent flow with + reacting Lagrangian parcels, and surface film modelling. + +\*---------------------------------------------------------------------------*/ + +#include "fvCFD.H" +#include "hCombustionThermo.H" +#include "turbulenceModel.H" +#include "BasicReactingCloud.H" +#include "surfaceFilmModel.H" +#include "psiChemistryModel.H" +#include "chemistrySolver.H" +#include "radiationModel.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +int main(int argc, char *argv[]) +{ + #include "setRootCase.H" + + #include "createTime.H" + #include "createMesh.H" + #include "readChemistryProperties.H" + #include "readGravitationalAcceleration.H" + #include "createFields.H" + #include "createClouds.H" + #include "createRadiationModel.H" + #include "createSurfaceFilmModel.H" + #include "initContinuityErrs.H" + #include "readTimeControls.H" + #include "compressibleCourantNo.H" + #include "setInitialDeltaT.H" + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + Info<< "\nStarting time loop\n" << endl; + + while (runTime.run()) + { + #include "readTimeControls.H" + #include "readPISOControls.H" + #include "compressibleCourantNo.H" + #include "setDeltaT.H" + + runTime++; + + Info<< "Time = " << runTime.timeName() << nl << endl; + + parcels.evolve(); + + surfaceFilm.evolve(); + + if (solvePrimaryRegion) + { + #include "chemistry.H" + #include "rhoEqn.H" + + // --- PIMPLE loop + for (int ocorr=1; ocorr<=nOuterCorr; ocorr++) + { + #include "UEqn.H" + #include "YEqn.H" + + // --- PISO loop + for (int corr=1; corr<=nCorr; corr++) + { + #include "hsEqn.H" + #include "pEqn.H" + } + + Info<< "T gas min/max = " << min(T).value() << ", " + << max(T).value() << endl; + } + + turbulence->correct(); + + rho = thermo.rho(); + + if (runTime.write()) + { + chemistry.dQ()().write(); + } + } + else + { + runTime.write(); + } + + Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" + << " ClockTime = " << runTime.elapsedClockTime() << " s" + << nl << endl; + } + + Info<< "End" << endl; + + return(0); +} + + +// ************************************************************************* // diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/readChemistryProperties.H b/applications/solvers/lagrangian/reactingParcelFilmFoam/readChemistryProperties.H new file mode 100644 index 0000000000..ab51afe283 --- /dev/null +++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/readChemistryProperties.H @@ -0,0 +1,23 @@ +Info<< "Reading chemistry properties\n" << endl; + +IOdictionary chemistryProperties +( + IOobject + ( + "chemistryProperties", + runTime.constant(), + mesh, + IOobject::MUST_READ, + IOobject::NO_WRITE, + false + ) +); + +Switch turbulentReaction(chemistryProperties.lookup("turbulentReaction")); + +dimensionedScalar Cmix("Cmix", dimless, 1.0); + +if (turbulentReaction) +{ + chemistryProperties.lookup("Cmix") >> Cmix; +} diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/rhoEqn.H b/applications/solvers/lagrangian/reactingParcelFilmFoam/rhoEqn.H new file mode 100644 index 0000000000..4e711c31ec --- /dev/null +++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/rhoEqn.H @@ -0,0 +1,43 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2008-2010 OpenCFD Ltd. + \\/ 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 2 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, write to the Free Software Foundation, + Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +Global + rhoEqn + +Description + Solve the continuity for density. + +\*---------------------------------------------------------------------------*/ + +{ + solve + ( + fvm::ddt(rho) + + fvc::div(phi) + == + parcels.Srho() + ); +} + +// ************************************************************************* //