From 33b24e6599e40b4b03fb76d09bd91270ad4d6e11 Mon Sep 17 00:00:00 2001 From: andy Date: Mon, 25 Oct 2010 18:28:25 +0100 Subject: [PATCH] ENH: Added new steadyReactingParcelFoam solver --- .../steadyReactingParcelFoam/Make/files | 3 + .../steadyReactingParcelFoam/Make/options | 43 ++++ .../steadyReactingParcelFoam/UEqn.H | 19 ++ .../steadyReactingParcelFoam/YEqn.H | 47 ++++ .../steadyReactingParcelFoam/chemistry.H | 27 +++ .../steadyReactingParcelFoam/createClouds.H | 9 + .../createExplicitSources.H | 27 +++ .../steadyReactingParcelFoam/createFields.H | 149 ++++++++++++ .../createPorousZones.H | 3 + .../steadyReactingParcelFoam/hsEqn.H | 23 ++ .../steadyReactingParcelFoam/pEqn.H | 60 +++++ .../readAdditionalSolutionControls.H | 7 + .../readChemistryProperties.H | 23 ++ .../readTimeControls.H | 71 ++++++ .../steadyReactingParcelFoam/rhoEqn.H | 45 ++++ .../setPressureWork.H | 10 + .../steadyReactingParcelFoam.C | 123 ++++++++++ .../steadyReactingParcelFoam/timeScales.H | 229 ++++++++++++++++++ 18 files changed, 918 insertions(+) create mode 100644 applications/solvers/lagrangian/steadyReactingParcelFoam/Make/files create mode 100644 applications/solvers/lagrangian/steadyReactingParcelFoam/Make/options create mode 100644 applications/solvers/lagrangian/steadyReactingParcelFoam/UEqn.H create mode 100644 applications/solvers/lagrangian/steadyReactingParcelFoam/YEqn.H create mode 100644 applications/solvers/lagrangian/steadyReactingParcelFoam/chemistry.H create mode 100644 applications/solvers/lagrangian/steadyReactingParcelFoam/createClouds.H create mode 100644 applications/solvers/lagrangian/steadyReactingParcelFoam/createExplicitSources.H create mode 100644 applications/solvers/lagrangian/steadyReactingParcelFoam/createFields.H create mode 100644 applications/solvers/lagrangian/steadyReactingParcelFoam/createPorousZones.H create mode 100644 applications/solvers/lagrangian/steadyReactingParcelFoam/hsEqn.H create mode 100644 applications/solvers/lagrangian/steadyReactingParcelFoam/pEqn.H create mode 100644 applications/solvers/lagrangian/steadyReactingParcelFoam/readAdditionalSolutionControls.H create mode 100644 applications/solvers/lagrangian/steadyReactingParcelFoam/readChemistryProperties.H create mode 100644 applications/solvers/lagrangian/steadyReactingParcelFoam/readTimeControls.H create mode 100644 applications/solvers/lagrangian/steadyReactingParcelFoam/rhoEqn.H create mode 100644 applications/solvers/lagrangian/steadyReactingParcelFoam/setPressureWork.H create mode 100644 applications/solvers/lagrangian/steadyReactingParcelFoam/steadyReactingParcelFoam.C create mode 100644 applications/solvers/lagrangian/steadyReactingParcelFoam/timeScales.H diff --git a/applications/solvers/lagrangian/steadyReactingParcelFoam/Make/files b/applications/solvers/lagrangian/steadyReactingParcelFoam/Make/files new file mode 100644 index 0000000000..a64e1ee786 --- /dev/null +++ b/applications/solvers/lagrangian/steadyReactingParcelFoam/Make/files @@ -0,0 +1,3 @@ +steadyReactingParcelFoam.C + +EXE = $(FOAM_APPBIN)/steadyReactingParcelFoam diff --git a/applications/solvers/lagrangian/steadyReactingParcelFoam/Make/options b/applications/solvers/lagrangian/steadyReactingParcelFoam/Make/options new file mode 100644 index 0000000000..6e117fc63b --- /dev/null +++ b/applications/solvers/lagrangian/steadyReactingParcelFoam/Make/options @@ -0,0 +1,43 @@ +EXE_INC = \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I${LIB_SRC}/meshTools/lnInclude \ + -I$(LIB_SRC)/turbulenceModels/compressible/turbulenceModel \ + -I$(LIB_SRC)/lagrangian/basic/lnInclude \ + -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \ + -I$(LIB_SRC)/lagrangian/coalCombustion/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/pdfs/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/specie/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/liquids/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/liquidMixture/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/solids/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/solidMixture/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/thermophysicalFunctions/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/reactionThermo/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/SLGThermo/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/chemistryModel/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/radiation/lnInclude \ + -I$(LIB_SRC)/ODE/lnInclude \ + -I$(LIB_SRC)/surfaceFilmModels/lnInclude + +EXE_LIBS = \ + -lfiniteVolume \ + -lmeshTools \ + -lcompressibleTurbulenceModel \ + -lcompressibleRASModels \ + -lcompressibleLESModels \ + -llagrangian \ + -llagrangianIntermediate \ + -lspecie \ + -lbasicThermophysicalModels \ + -lliquids \ + -lliquidMixture \ + -lsolids \ + -lsolidMixture \ + -lthermophysicalFunctions \ + -lreactionThermophysicalModels \ + -lSLGThermo \ + -lchemistryModel \ + -lradiation \ + -lODE \ + -lsurfaceFilmModels diff --git a/applications/solvers/lagrangian/steadyReactingParcelFoam/UEqn.H b/applications/solvers/lagrangian/steadyReactingParcelFoam/UEqn.H new file mode 100644 index 0000000000..fee0fe1a68 --- /dev/null +++ b/applications/solvers/lagrangian/steadyReactingParcelFoam/UEqn.H @@ -0,0 +1,19 @@ + fvVectorMatrix UEqn + ( +// fvm::ddt(rho, U) + pZones.ddt(rho, U) + + fvm::div(phi, U) + + turbulence->divDevRhoReff(U) + == + rho.dimensionedInternalField()*g + + parcels.SU(U) + + momentumSource.Su() + ); + + pZones.addResistance(UEqn); + + if (momentumPredictor) + { + solve(UEqn == -fvc::grad(p)); + } + diff --git a/applications/solvers/lagrangian/steadyReactingParcelFoam/YEqn.H b/applications/solvers/lagrangian/steadyReactingParcelFoam/YEqn.H new file mode 100644 index 0000000000..4e2f9766a8 --- /dev/null +++ b/applications/solvers/lagrangian/steadyReactingParcelFoam/YEqn.H @@ -0,0 +1,47 @@ + +tmp > mvConvection +( + fv::convectionScheme::New + ( + mesh, + fields, + phi, + mesh.divScheme("div(phi,Yi_h)") + ) +); + + +if (solveSpecies) +{ + label inertIndex = -1; + volScalarField Yt = 0.0*Y[0]; + + forAll(Y, i) + { + if (Y[i].name() != inertSpecie) + { + volScalarField& Yi = Y[i]; + solve + ( + fvm::ddt(rho, Yi) + + mvConvection->fvmDiv(phi, Yi) + - fvm::laplacian(turbulence->muEff(), Yi) + == + parcels.SYi(i, Yi) + + kappa*chemistry.RR(i)().dimensionedInternalField() + + massSource.Su(i), + 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/steadyReactingParcelFoam/chemistry.H b/applications/solvers/lagrangian/steadyReactingParcelFoam/chemistry.H new file mode 100644 index 0000000000..50fc7f575b --- /dev/null +++ b/applications/solvers/lagrangian/steadyReactingParcelFoam/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/steadyReactingParcelFoam/createClouds.H b/applications/solvers/lagrangian/steadyReactingParcelFoam/createClouds.H new file mode 100644 index 0000000000..954b74e069 --- /dev/null +++ b/applications/solvers/lagrangian/steadyReactingParcelFoam/createClouds.H @@ -0,0 +1,9 @@ +Info<< "\nConstructing reacting cloud" << endl; +basicReactingMultiphaseCloud parcels +( + "reactingCloud1", + rho, + U, + g, + slgThermo +); diff --git a/applications/solvers/lagrangian/steadyReactingParcelFoam/createExplicitSources.H b/applications/solvers/lagrangian/steadyReactingParcelFoam/createExplicitSources.H new file mode 100644 index 0000000000..0f2b52e536 --- /dev/null +++ b/applications/solvers/lagrangian/steadyReactingParcelFoam/createExplicitSources.H @@ -0,0 +1,27 @@ +Info<< "Creating mass source\n" << endl; +scalarTimeActivatedExplicitSourceList massSource +( + "mass", + mesh, + dimMass/dimTime/dimVolume, + composition.species() +); + + +Info<< "Creating momentum source\n" << endl; +vectorTimeActivatedExplicitSourceList momentumSource +( + "momentum", + mesh, + dimMass*dimVelocity/dimTime/dimVolume, + "U" +); + +Info<< "Creating energy source\n" << endl; +scalarTimeActivatedExplicitSourceList energySource +( + "energy", + mesh, + dimEnergy/dimTime/dimVolume, + "h" +); diff --git a/applications/solvers/lagrangian/steadyReactingParcelFoam/createFields.H b/applications/solvers/lagrangian/steadyReactingParcelFoam/createFields.H new file mode 100644 index 0000000000..5f4aaa2427 --- /dev/null +++ b/applications/solvers/lagrangian/steadyReactingParcelFoam/createFields.H @@ -0,0 +1,149 @@ + Info<< "Reading thermophysical properties\n" << endl; + + autoPtr pChemistry + ( + rhoChemistryModel::New(mesh) + ); + rhoChemistryModel& chemistry = pChemistry(); + + hsReactionThermo& thermo = chemistry.thermo(); + + SLGThermo slgThermo(mesh, thermo); + + basicMultiComponentMixture& composition = thermo.composition(); + PtrList& Y = composition.Y(); + + const word inertSpecie(thermo.lookup("inertSpecie")); + + if (!composition.contains(inertSpecie)) + { + FatalErrorIn(args.executable()) + << "Specified inert specie '" << inertSpecie << "' not found in " + << "species list. Available species:" << composition.species() + << exit(FatalError); + } + + volScalarField& p = thermo.p(); + volScalarField& hs = thermo.hs(); + const volScalarField& T = thermo.T(); + const volScalarField& psi = thermo.psi(); + + 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" + + DimensionedField kappa + ( + IOobject + ( + "kappa", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + mesh, + dimensionedScalar("zero", dimless, 0.0) + ); + + dimensionedScalar rhoMax + ( + mesh.solutionDict().subDict("PIMPLE").lookup("rhoMax") + ); + + dimensionedScalar rhoMin + ( + mesh.solutionDict().subDict("PIMPLE").lookup("rhoMin") + ); + + Info<< "Creating turbulence model\n" << endl; + autoPtr turbulence + ( + compressible::turbulenceModel::New + ( + rho, + U, + phi, + thermo + ) + ); + + Info<< "Creating multi-variate interpolation scheme\n" << endl; + multivariateSurfaceInterpolationScheme::fieldTable fields; + + forAll(Y, i) + { + fields.add(Y[i]); + } + fields.add(hs); + + DimensionedField chemistrySh + ( + IOobject + ( + "chemistry::Sh", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh, + dimensionedScalar("chemistrySh", dimEnergy/dimTime/dimVolume, 0.0) + ); + + volScalarField invTauFlow + ( + IOobject + ( + "invTauFlow", + runTime.timeName(), + mesh, + IOobject::READ_IF_PRESENT, + IOobject::AUTO_WRITE + ), + mesh, + dimensionedScalar("one", dimless/dimTime, 1), + zeroGradientFvPatchScalarField::typeName + ); + + Info<< "Creating field DpDt\n" << endl; + volScalarField DpDt + ( + IOobject + ( + "DpDt", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE + ), + mesh, + dimensionedScalar("zero", dimPressure/dimTime, 0.0) + ); + + #include "setPressureWork.H" diff --git a/applications/solvers/lagrangian/steadyReactingParcelFoam/createPorousZones.H b/applications/solvers/lagrangian/steadyReactingParcelFoam/createPorousZones.H new file mode 100644 index 0000000000..90506856d2 --- /dev/null +++ b/applications/solvers/lagrangian/steadyReactingParcelFoam/createPorousZones.H @@ -0,0 +1,3 @@ + Info<< "Creating porous zones" << nl << endl; + + porousZones pZones(mesh); diff --git a/applications/solvers/lagrangian/steadyReactingParcelFoam/hsEqn.H b/applications/solvers/lagrangian/steadyReactingParcelFoam/hsEqn.H new file mode 100644 index 0000000000..5954b1217e --- /dev/null +++ b/applications/solvers/lagrangian/steadyReactingParcelFoam/hsEqn.H @@ -0,0 +1,23 @@ +{ + fvScalarMatrix hsEqn + ( + fvm::ddt(rho, hs) + + mvConvection->fvmDiv(phi, hs) + - fvm::laplacian(turbulence->alphaEff(), hs) + == + DpDt + + parcels.Sh(hs) + + radiation->Shs(thermo) + + energySource.Su() + + chemistrySh + ); + + hsEqn.solve(); + + thermo.correct(); + + radiation->correct(); + + Info<< "T gas min/max = " << min(T).value() << ", " + << max(T).value() << endl; +} diff --git a/applications/solvers/lagrangian/steadyReactingParcelFoam/pEqn.H b/applications/solvers/lagrangian/steadyReactingParcelFoam/pEqn.H new file mode 100644 index 0000000000..6ae2712ea4 --- /dev/null +++ b/applications/solvers/lagrangian/steadyReactingParcelFoam/pEqn.H @@ -0,0 +1,60 @@ +{ + rho = thermo.rho(); + + // Thermodynamic density needs to be updated by psi*d(p) after the + // pressure solution - done in 2 parts. Part 1: + thermo.rho() -= psi*p; + + volScalarField rAU = 1.0/UEqn.A(); + U = rAU*UEqn.H(); + + if (pZones.size() > 0) + { + // ddtPhiCorr not well defined for cases with porosity + phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf()); + } + else + { + phi = + fvc::interpolate(rho) + *( + (fvc::interpolate(U) & mesh.Sf()) + + fvc::ddtPhiCorr(rAU, rho, U, phi) + ); + } + + for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) + { + fvScalarMatrix pEqn + ( + fvc::ddt(rho) + psi*correction(fvm::ddt(p)) + + fvc::div(phi) + - fvm::laplacian(rho*rAU, p) + == + parcels.Srho() + + massSource.SuTot() + ); + + pEqn.solve(); + + if (nonOrth == nNonOrthCorr) + { + phi += pEqn.flux(); + } + } + + // Second part of thermodynamic density update + thermo.rho() += psi*p; + + #include "rhoEqn.H" // NOTE: flux and time scales now inconsistent + #include "compressibleContinuityErrs.H" + + U -= rAU*fvc::grad(p); + U.correctBoundaryConditions(); + + rho = thermo.rho(); + rho = max(rho, rhoMin); + rho = min(rho, rhoMax); + + #include "setPressureWork.H" +} diff --git a/applications/solvers/lagrangian/steadyReactingParcelFoam/readAdditionalSolutionControls.H b/applications/solvers/lagrangian/steadyReactingParcelFoam/readAdditionalSolutionControls.H new file mode 100644 index 0000000000..f0763d8421 --- /dev/null +++ b/applications/solvers/lagrangian/steadyReactingParcelFoam/readAdditionalSolutionControls.H @@ -0,0 +1,7 @@ +dictionary additional = mesh.solutionDict().subDict("additional"); + +bool eWork = additional.lookupOrDefault("eWork", true); +bool hWork = additional.lookupOrDefault("hWork", true); + +// flag to activate solve transport for each specie (Y vector) +bool solveSpecies = additional.lookupOrDefault("solveSpecies", true); diff --git a/applications/solvers/lagrangian/steadyReactingParcelFoam/readChemistryProperties.H b/applications/solvers/lagrangian/steadyReactingParcelFoam/readChemistryProperties.H new file mode 100644 index 0000000000..f0bcf7597f --- /dev/null +++ b/applications/solvers/lagrangian/steadyReactingParcelFoam/readChemistryProperties.H @@ -0,0 +1,23 @@ +Info<< "Reading chemistry properties\n" << endl; + +IOdictionary chemistryProperties +( + IOobject + ( + "chemistryProperties", + runTime.constant(), + mesh, + IOobject::MUST_READ_IF_MODIFIED, + 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/steadyReactingParcelFoam/readTimeControls.H b/applications/solvers/lagrangian/steadyReactingParcelFoam/readTimeControls.H new file mode 100644 index 0000000000..b75a191d6a --- /dev/null +++ b/applications/solvers/lagrangian/steadyReactingParcelFoam/readTimeControls.H @@ -0,0 +1,71 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2010-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 + +\*---------------------------------------------------------------------------*/ + +// Maximum flow Courant number +scalar maxCo(readScalar(runTime.controlDict().lookup("maxCo"))); + +// Maximum time scale +scalar maxDeltaT = readScalar(runTime.controlDict().lookup("maxDeltaT")); + +// Smoothing parameter (0-1) when smoothing iterations > 0 +scalar alphaTauSmooth +( + runTime.controlDict().lookupOrDefault("alphaTauSmooth", 0.1) +); + +// Maximum change in cell density per iteration (relative to previous value) +scalar alphaTauRho +( + runTime.controlDict().lookupOrDefault("alphaTauRho", 0.05) +); + +// Maximum change in cell velocity per iteration (relative to previous value) +scalar alphaTauU +( + runTime.controlDict().lookupOrDefault("alphaTauU", 0.05) +); + +// Maximum change in cell temperature per iteration (relative to previous value) +scalar alphaTauTemp +( + runTime.controlDict().lookupOrDefault("alphaTauTemp", 0.05) +); + +// Max specie mass fraction that can be consumed/gained per chemistry +// integration step +scalar alphaTauSpecie +( + runTime.controlDict().lookupOrDefault("alphaTauSpecie", 0.05) +); + +// Maximum unboundedness allowed (fraction of 1) +scalar specieMaxUnbound +( + runTime.controlDict().lookupOrDefault("specieMaxUnbound", 0.01) +); + + +// ************************************************************************* // diff --git a/applications/solvers/lagrangian/steadyReactingParcelFoam/rhoEqn.H b/applications/solvers/lagrangian/steadyReactingParcelFoam/rhoEqn.H new file mode 100644 index 0000000000..33b7943cad --- /dev/null +++ b/applications/solvers/lagrangian/steadyReactingParcelFoam/rhoEqn.H @@ -0,0 +1,45 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 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 . + +Global + rhoEqn + +Description + Solve the continuity for density. + +\*---------------------------------------------------------------------------*/ + +{ + fvScalarMatrix rhoEqn + ( + fvm::ddt(rho) + + fvc::div(phi) + == + parcels.Srho(rho) + + massSource.SuTot() + ); + + rhoEqn.solve(); +} + +// ************************************************************************* // diff --git a/applications/solvers/lagrangian/steadyReactingParcelFoam/setPressureWork.H b/applications/solvers/lagrangian/steadyReactingParcelFoam/setPressureWork.H new file mode 100644 index 0000000000..69d9a21731 --- /dev/null +++ b/applications/solvers/lagrangian/steadyReactingParcelFoam/setPressureWork.H @@ -0,0 +1,10 @@ +DpDt == dimensionedScalar("zero", DpDt.dimensions(), 0.0); + +if (eWork) +{ + DpDt += -p*fvc::div(phi/fvc::interpolate(rho)); +} +if (hWork) +{ + DpDt += fvc::div(phi/fvc::interpolate(rho)*fvc::interpolate(p)); +} diff --git a/applications/solvers/lagrangian/steadyReactingParcelFoam/steadyReactingParcelFoam.C b/applications/solvers/lagrangian/steadyReactingParcelFoam/steadyReactingParcelFoam.C new file mode 100644 index 0000000000..a6aa1290a3 --- /dev/null +++ b/applications/solvers/lagrangian/steadyReactingParcelFoam/steadyReactingParcelFoam.C @@ -0,0 +1,123 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 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 + porousExplicitSourceReactingParcelFoam + +Description + Transient PISO solver for compressible, laminar or turbulent flow with + reacting multiphase Lagrangian parcels for porous media, including explicit + sources for mass, momentum and energy + + The solver includes: + - reacting multiphase parcel cloud + - porous media + - mass, momentum and energy sources + + Note: ddtPhiCorr not used here when porous zones are active + - not well defined for porous calculations + +\*---------------------------------------------------------------------------*/ + +#include "fvCFD.H" +#include "hReactionThermo.H" +#include "turbulenceModel.H" +#include "basicReactingMultiphaseCloud.H" +#include "rhoChemistryModel.H" +#include "chemistrySolver.H" +#include "radiationModel.H" +#include "porousZones.H" +#include "timeActivatedExplicitSource.H" +#include "SLGThermo.H" +#include "fvcSmooth.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +int main(int argc, char *argv[]) +{ + #include "setRootCase.H" + + #include "createTime.H" + #include "createMesh.H" + #include "readGravitationalAcceleration.H" + #include "readTimeControls.H" + #include "readAdditionalSolutionControls.H" + #include "createFields.H" + #include "createRadiationModel.H" + #include "createClouds.H" + #include "createExplicitSources.H" + #include "createPorousZones.H" + #include "initContinuityErrs.H" + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + Info<< "\nStarting time loop\n" << endl; + + while (runTime.run()) + { + #include "readSIMPLEControls.H" + #include "readChemistryProperties.H" + #include "readAdditionalSolutionControls.H" + #include "readTimeControls.H" + + runTime++; + + Info<< "Time = " << runTime.timeName() << nl << endl; + + p.storePrevIter(); + + // --- Pressure-velocity corrector + { + parcels.evolve(); + + #include "chemistry.H" + + #include "timeScales.H" + + #include "rhoEqn.H" + #include "UEqn.H" + #include "YEqn.H" + #include "hsEqn.H" + + #include "pEqn.H" + + turbulence->correct(); + } + + if (runTime.write()) + { + chemistry.dQ()().write(); + } + + Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" + << " ClockTime = " << runTime.elapsedClockTime() << " s" + << nl << endl; + } + + Info<< "End\n" << endl; + + return(0); +} + + +// ************************************************************************* // diff --git a/applications/solvers/lagrangian/steadyReactingParcelFoam/timeScales.H b/applications/solvers/lagrangian/steadyReactingParcelFoam/timeScales.H new file mode 100644 index 0000000000..35feac5a30 --- /dev/null +++ b/applications/solvers/lagrangian/steadyReactingParcelFoam/timeScales.H @@ -0,0 +1,229 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | + \\ / A nd | Copyright (C) 2010-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 + +\*---------------------------------------------------------------------------*/ + +Info<< "Time scales min/max:" << endl; + + +{ + // Cache old time scale field + tmp tinvTauFlow0 + ( + new volScalarField + ( + IOobject + ( + "invTauFlow0", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ), + invTauFlow + ) + ); + const volScalarField& invTauFlow0 = tinvTauFlow0(); + + + // Flow time scale + // ~~~~~~~~~~~~~~~ + { + invTauFlow = + fvc::surfaceSum + ( + mag(phi)*mesh.deltaCoeffs()/(maxCo*mesh.magSf()) + ) + /rho; + + invTauFlow.max(1.0/maxDeltaT); + + Info<< " Flow = " << gMin(1/invTauFlow.internalField()) << ", " + << gMax(1/invTauFlow.internalField()) << endl; + } + + + // Mass source time scale + // ~~~~~~~~~~~~~~~~~~~~~~ + + { + scalarField tau = + runTime.deltaTValue()*mag(parcels.Srho() + massSource.SuTot()); + + tau = alphaTauRho*rho/(tau + ROOTVSMALL); + + Info<< " Density = " << min(maxDeltaT, gMin(tau)) << ", " + << min(maxDeltaT, gMax(tau)) << endl; + + invTauFlow.internalField() = max(invTauFlow.internalField(), 1/tau); + } + + + // Momentum source time scale + // ~~~~~~~~~~~~~~~~~~~~~~~~~~ + + { +/* + // Method 1 - mag(U) limit using 'small' nominal velocity + scalarField tau = + runTime.deltaTValue() + *mag + ( + rho.dimensionedInternalField()*g + + parcels.UTrans()/(mesh.V()*runTime.deltaT()) + + momentumSource.Su() + ) + /rho; + + const scalar nomMagU(dimensionedScalar("1", dimVelocity, 1)); + tau = alphaTauU*(nomMagU + mag(U))/(tau + ROOTVSMALL); +*/ +/* + // Method 2 - based on fluxes and Co-like limit + volVectorField UEqnRhs + ( + IOobject + ( + "UEqnRhs", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::NO_WRITE, + false + ), + mesh, + dimensionedVector("zero", dimDensity*dimAcceleration, vector::zero), + zeroGradientFvPatchVectorField::typeName + ); + + UEqnRhs.internalField() = + rho.dimensionedInternalField()*g + + parcels.UTrans()/(mesh.V()*runTime.deltaT()) + + momentumSource.Su(); + + surfaceScalarField phiSU + ( + "phiSU", + fvc::interpolate(runTime.deltaT()*UEqnRhs) & mesh.Sf() + ); + + scalarField tau = + alphaTauU*rho + /fvc::surfaceSum + ( + mag(phi + phiSU)*mesh.deltaCoeffs()/mesh.magSf() + + dimensionedScalar("SMALL", dimDensity/dimTime, ROOTVSMALL) + ); + +*/ +/* + Info<< " Momentum = " << min(maxDeltaT, gMin(tau)) << ", " + << min(maxDeltaT, gMax(tau)) << endl; + + invTauFlow.internalField() = max(invTauFlow.internalField(), 1/tau); +*/ + } + + + // Temperature source time scale + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + + { + scalarField tau = + runTime.deltaTValue() + *mag + ( + DpDt + + parcels.hsTrans()/(mesh.V()*runTime.deltaT()) + + energySource.Su() + + chemistrySh + ) + /rho; + + tau = alphaTauTemp*thermo.Cp()*T/(tau + ROOTVSMALL); + + Info<< " Temperature = " << min(maxDeltaT, gMin(tau)) << ", " + << min(maxDeltaT, gMax(tau)) << endl; + + invTauFlow.internalField() = max(invTauFlow.internalField(), 1/tau); + } + + + // Specie source time scale + // ~~~~~~~~~~~~~~~~~~~~~~~~ + + { + scalarField tau(mesh.nCells(), ROOTVGREAT); + forAll(Y, fieldI) + { + const volScalarField& Yi = Y[fieldI]; + const scalarField deltaYi = + runTime.deltaTValue() + *mag + ( + kappa*chemistry.RR(fieldI)() + + massSource.Su(fieldI) + + parcels.Srho(fieldI) + ) + /rho; + tau = + min + ( + tau, + alphaTauSpecie + /(deltaYi/(Yi + specieMaxUnbound) + ROOTVSMALL) + ); + } + + Info<< " Specie = " << min(maxDeltaT, gMin(tau)) << ", " + << min(maxDeltaT, gMax(tau)) << endl; + + invTauFlow.internalField() = max(invTauFlow.internalField(), 1/tau); + } + + + // Limit rate of change of time scale + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + // - reduce as much as required for flow, but limit source contributions + const dimensionedScalar deltaTRamp("deltaTRamp", dimless, 1/(1 + 0.2)); + invTauFlow = max(invTauFlow, invTauFlow0*deltaTRamp); + tinvTauFlow0.clear(); + + + // Limit the largest time scale + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + invTauFlow.max(1/maxDeltaT); + + + // Spatially smooth the time scale field + // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + fvc::smooth(invTauFlow, alphaTauSmooth); + + Info<< " Overall = " << min(1/invTauFlow).value() + << ", " << max(1/invTauFlow).value() << nl << endl; +} + + +// ************************************************************************* //