diff --git a/applications/solvers/lagrangian/coalChemistryFoam/UEqn.H b/applications/solvers/lagrangian/coalChemistryFoam/UEqn.H index 5851cdcf95..a6d7127800 100644 --- a/applications/solvers/lagrangian/coalChemistryFoam/UEqn.H +++ b/applications/solvers/lagrangian/coalChemistryFoam/UEqn.H @@ -5,8 +5,8 @@ + turbulence->divDevRhoReff(U) == rho.dimensionedInternalField()*g - + coalParcels.SU() - + limestoneParcels.SU() + + coalParcels.SU(U) + + limestoneParcels.SU(U) ); UEqn.relax(); diff --git a/applications/solvers/lagrangian/coalChemistryFoam/YEqn.H b/applications/solvers/lagrangian/coalChemistryFoam/YEqn.H index abe03f60c0..0b04bdabba 100644 --- a/applications/solvers/lagrangian/coalChemistryFoam/YEqn.H +++ b/applications/solvers/lagrangian/coalChemistryFoam/YEqn.H @@ -25,7 +25,7 @@ tmp > mvConvection + mvConvection->fvmDiv(phi, Yi) - fvm::laplacian(turbulence->muEff(), Yi) == - coalParcels.Srho(i) + coalParcels.SYi(i, Yi) + kappa*chemistry.RR(i)().dimensionedInternalField() ); diff --git a/applications/solvers/lagrangian/coalChemistryFoam/hsEqn.H b/applications/solvers/lagrangian/coalChemistryFoam/hsEqn.H index 40eb26709e..2dcb13c79f 100644 --- a/applications/solvers/lagrangian/coalChemistryFoam/hsEqn.H +++ b/applications/solvers/lagrangian/coalChemistryFoam/hsEqn.H @@ -6,8 +6,8 @@ - fvm::laplacian(turbulence->alphaEff(), hs) == DpDt - + coalParcels.Sh() - + limestoneParcels.Sh() + + coalParcels.Sh(hs) + + limestoneParcels.Sh(hs) + enthalpySource.Su() + radiation->Shs(thermo) + chemistrySh diff --git a/applications/solvers/lagrangian/coalChemistryFoam/rhoEqn.H b/applications/solvers/lagrangian/coalChemistryFoam/rhoEqn.H index 08fbb120f5..5e4bb09d84 100644 --- a/applications/solvers/lagrangian/coalChemistryFoam/rhoEqn.H +++ b/applications/solvers/lagrangian/coalChemistryFoam/rhoEqn.H @@ -35,7 +35,7 @@ Description fvm::ddt(rho) + fvc::div(phi) == - coalParcels.Srho() + coalParcels.Srho(rho) ); } diff --git a/applications/solvers/lagrangian/incompressibleUncoupledKinematicParcelDyMFoam/incompressibleUncoupledKinematicParcelDyMFoam.C b/applications/solvers/lagrangian/incompressibleUncoupledKinematicParcelDyMFoam/incompressibleUncoupledKinematicParcelDyMFoam.C index f5fba006e0..40b642a284 100644 --- a/applications/solvers/lagrangian/incompressibleUncoupledKinematicParcelDyMFoam/incompressibleUncoupledKinematicParcelDyMFoam.C +++ b/applications/solvers/lagrangian/incompressibleUncoupledKinematicParcelDyMFoam/incompressibleUncoupledKinematicParcelDyMFoam.C @@ -51,7 +51,7 @@ int main(int argc, char *argv[]) #include "setRootCase.H" #include "createTime.H" - # include "createDynamicFvMesh.H" + #include "createDynamicFvMesh.H" #include "readGravitationalAcceleration.H" #include "createFields.H" diff --git a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/UEqn.H b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/UEqn.H index e77fe75dab..fee0fe1a68 100644 --- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/UEqn.H +++ b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/UEqn.H @@ -6,7 +6,7 @@ + turbulence->divDevRhoReff(U) == rho.dimensionedInternalField()*g - + parcels.SU() + + parcels.SU(U) + momentumSource.Su() ); diff --git a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/YEqn.H b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/YEqn.H index f54be04bb8..ac369f3df4 100644 --- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/YEqn.H +++ b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/YEqn.H @@ -26,7 +26,7 @@ tmp > mvConvection + mvConvection->fvmDiv(phi, Yi) - fvm::laplacian(turbulence->muEff(), Yi) == - parcels.Srho(i) + parcels.SYi(i, Yi) + kappa*chemistry.RR(i)().dimensionedInternalField() + massSource.Su(i), mesh.solver("Yi") diff --git a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/hsEqn.H b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/hsEqn.H index 73c276c294..f448144f16 100644 --- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/hsEqn.H +++ b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/hsEqn.H @@ -37,7 +37,7 @@ - fvm::laplacian(turbulence->alphaEff(), hs) == pWork() - + parcels.Sh() + + parcels.Sh(hs) + radiation->Shs(thermo) + energySource.Su() + chemistrySh diff --git a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/rhoEqn.H b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/rhoEqn.H index 0dea187eda..b1770a5a9d 100644 --- a/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/rhoEqn.H +++ b/applications/solvers/lagrangian/porousExplicitSourceReactingParcelFoam/rhoEqn.H @@ -35,7 +35,7 @@ Description fvm::ddt(rho) + fvc::div(phi) == - parcels.Srho() + parcels.Srho(rho) + massSource.SuTot() ); diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/UEqn.H b/applications/solvers/lagrangian/reactingParcelFilmFoam/UEqn.H index 3c4a927091..2e8f979be4 100644 --- a/applications/solvers/lagrangian/reactingParcelFilmFoam/UEqn.H +++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/UEqn.H @@ -5,7 +5,7 @@ + turbulence->divDevRhoReff(U) == rho.dimensionedInternalField()*g - + parcels.SU() + + parcels.SU(U) ); UEqn.relax(); diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/YEqn.H b/applications/solvers/lagrangian/reactingParcelFilmFoam/YEqn.H index 0a2c36b7b5..da9a289d69 100644 --- a/applications/solvers/lagrangian/reactingParcelFilmFoam/YEqn.H +++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/YEqn.H @@ -25,7 +25,7 @@ tmp > mvConvection + mvConvection->fvmDiv(phi, Yi) - fvm::laplacian(turbulence->muEff(), Yi) == - parcels.Srho(i) + parcels.SYi(i, Yi) + surfaceFilm.Srho(i) + kappa*chemistry.RR(i)().dimensionedInternalField(), mesh.solver("Yi") diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/hsEqn.H b/applications/solvers/lagrangian/reactingParcelFilmFoam/hsEqn.H index 3c76f1384c..0cb2318252 100644 --- a/applications/solvers/lagrangian/reactingParcelFilmFoam/hsEqn.H +++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/hsEqn.H @@ -6,7 +6,7 @@ - fvm::laplacian(turbulence->alphaEff(), hs) == DpDt - + parcels.Sh() + + parcels.Sh(hs) + surfaceFilm.Sh() + radiation->Shs(thermo) + chemistrySh diff --git a/applications/solvers/lagrangian/reactingParcelFilmFoam/rhoEqn.H b/applications/solvers/lagrangian/reactingParcelFilmFoam/rhoEqn.H index 2c5d41f4fe..617f1df8a6 100644 --- a/applications/solvers/lagrangian/reactingParcelFilmFoam/rhoEqn.H +++ b/applications/solvers/lagrangian/reactingParcelFilmFoam/rhoEqn.H @@ -35,7 +35,7 @@ Description fvm::ddt(rho) + fvc::div(phi) == - parcels.Srho() + parcels.Srho(rho) + surfaceFilm.Srho() ); } diff --git a/applications/solvers/lagrangian/reactingParcelFoam/UEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/UEqn.H index 3c4a927091..2e8f979be4 100644 --- a/applications/solvers/lagrangian/reactingParcelFoam/UEqn.H +++ b/applications/solvers/lagrangian/reactingParcelFoam/UEqn.H @@ -5,7 +5,7 @@ + turbulence->divDevRhoReff(U) == rho.dimensionedInternalField()*g - + parcels.SU() + + parcels.SU(U) ); UEqn.relax(); diff --git a/applications/solvers/lagrangian/reactingParcelFoam/YEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/YEqn.H index c687f2035b..c4a929c449 100644 --- a/applications/solvers/lagrangian/reactingParcelFoam/YEqn.H +++ b/applications/solvers/lagrangian/reactingParcelFoam/YEqn.H @@ -25,7 +25,7 @@ tmp > mvConvection + mvConvection->fvmDiv(phi, Yi) - fvm::laplacian(turbulence->muEff(), Yi) == - parcels.Srho(i) + parcels.SYi(i, Yi) + kappa*chemistry.RR(i)().dimensionedInternalField(), mesh.solver("Yi") ); diff --git a/applications/solvers/lagrangian/reactingParcelFoam/hsEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/hsEqn.H index ce12ec3ad4..7821d340d4 100644 --- a/applications/solvers/lagrangian/reactingParcelFoam/hsEqn.H +++ b/applications/solvers/lagrangian/reactingParcelFoam/hsEqn.H @@ -6,7 +6,7 @@ - fvm::laplacian(turbulence->alphaEff(), hs) == DpDt - + parcels.Sh() + + parcels.Sh(hs) + radiation->Shs(thermo) + chemistrySh ); diff --git a/applications/solvers/lagrangian/reactingParcelFoam/rhoEqn.H b/applications/solvers/lagrangian/reactingParcelFoam/rhoEqn.H index b9b2d6a13f..d4f69d8f60 100644 --- a/applications/solvers/lagrangian/reactingParcelFoam/rhoEqn.H +++ b/applications/solvers/lagrangian/reactingParcelFoam/rhoEqn.H @@ -35,7 +35,7 @@ Description fvm::ddt(rho) + fvc::div(phi) == - parcels.Srho() + parcels.Srho(rho) ); } 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..ec92091de3 --- /dev/null +++ b/applications/solvers/lagrangian/steadyReactingParcelFoam/chemistry.H @@ -0,0 +1,28 @@ +if (chemistry.chemistry()) +{ + Info<< "Solving chemistry" << endl; + + // update reaction rates + chemistry.calculate(); + + // turbulent time scale + if (turbulentReaction) + { + typedef DimensionedField dsfType; + + const dimensionedScalar e0("e0", sqr(dimLength)/pow3(dimTime), SMALL); + + const dsfType tk = + Cmix*sqrt(turbulence->muEff()/rho/(turbulence->epsilon() + e0)); + + const dsfType tc = chemistry.tc()().dimensionedInternalField(); + + kappa = tc/(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..68734ab9a7 --- /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("SIMPLE").lookup("rhoMax") + ); + + dimensionedScalar rhoMin + ( + mesh.solutionDict().subDict("SIMPLE").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..e7d7c55dbb --- /dev/null +++ b/applications/solvers/lagrangian/steadyReactingParcelFoam/pEqn.H @@ -0,0 +1,65 @@ +{ + 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(); + } + } + + // Explicitly relax pressure for momentum corrector + p.relax(); + + Info<< "p min/max = " << min(p).value() << ", " << max(p).value() << endl; + + // 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..e742e9fea7 --- /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..0ab2be096e --- /dev/null +++ b/applications/solvers/lagrangian/steadyReactingParcelFoam/readTimeControls.H @@ -0,0 +1,70 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 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 . + +\*---------------------------------------------------------------------------*/ + +// 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/src/lagrangian/intermediate/submodels/Kinematic/PostProcessingModel/PostProcessingModel/PostProcessingModelI.H b/applications/solvers/lagrangian/steadyReactingParcelFoam/rhoEqn.H similarity index 72% rename from src/lagrangian/intermediate/submodels/Kinematic/PostProcessingModel/PostProcessingModel/PostProcessingModelI.H rename to applications/solvers/lagrangian/steadyReactingParcelFoam/rhoEqn.H index 1004d532b8..e42db72399 100644 --- a/src/lagrangian/intermediate/submodels/Kinematic/PostProcessingModel/PostProcessingModel/PostProcessingModelI.H +++ b/applications/solvers/lagrangian/steadyReactingParcelFoam/rhoEqn.H @@ -21,34 +21,28 @@ License 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. + \*---------------------------------------------------------------------------*/ -template -const Foam::dictionary& Foam::PostProcessingModel::dict() const { - return dict_; + fvScalarMatrix rhoEqn + ( + fvm::ddt(rho) + + fvc::div(phi) + == + parcels.Srho(rho) + + massSource.SuTot() + ); + + rhoEqn.solve(); + + Info<< "rho min/max = " << min(rho).value() << ", " << max(rho).value() + << endl; } - -template -const CloudType& Foam::PostProcessingModel::owner() const -{ - return owner_; -} - - -template -CloudType& Foam::PostProcessingModel::owner() -{ - return owner_; -} - - -template -const Foam::dictionary& Foam::PostProcessingModel::coeffDict() const -{ - return coeffDict_; -} - - // ************************************************************************* // 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..26ebe6fc67 --- /dev/null +++ b/applications/solvers/lagrangian/steadyReactingParcelFoam/timeScales.H @@ -0,0 +1,228 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 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 . + +\*---------------------------------------------------------------------------*/ + +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; +} + + +// ************************************************************************* // diff --git a/applications/utilities/postProcessing/lagrangian/steadyParticleTracks/Make/files b/applications/utilities/postProcessing/lagrangian/steadyParticleTracks/Make/files new file mode 100644 index 0000000000..b72bc1ae83 --- /dev/null +++ b/applications/utilities/postProcessing/lagrangian/steadyParticleTracks/Make/files @@ -0,0 +1,3 @@ +steadyParticleTracks.C + +EXE = $(FOAM_APPBIN)/steadyParticleTracks diff --git a/applications/utilities/postProcessing/lagrangian/steadyParticleTracks/Make/options b/applications/utilities/postProcessing/lagrangian/steadyParticleTracks/Make/options new file mode 100644 index 0000000000..004f0474b8 --- /dev/null +++ b/applications/utilities/postProcessing/lagrangian/steadyParticleTracks/Make/options @@ -0,0 +1,9 @@ +EXE_INC = \ + -I$(LIB_SRC)/lagrangian/basic/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude \ + -I$(LIB_SRC)/finiteVolume/lnInclude + +EXE_LIBS = \ + -llagrangian \ + -lmeshTools \ + -lfiniteVolume diff --git a/applications/utilities/postProcessing/lagrangian/steadyParticleTracks/createFields.H b/applications/utilities/postProcessing/lagrangian/steadyParticleTracks/createFields.H new file mode 100644 index 0000000000..24854ab6b6 --- /dev/null +++ b/applications/utilities/postProcessing/lagrangian/steadyParticleTracks/createFields.H @@ -0,0 +1,16 @@ +word dictName(args.optionLookupOrDefault("dict", "particleTrackDict")); + +IOdictionary propsDict +( + IOobject + ( + dictName, + runTime.constant(), + mesh, + IOobject::MUST_READ_IF_MODIFIED + ) +); + +word cloudName(propsDict.lookup("cloudName")); + +List userFields(propsDict.lookup("fields")); \ No newline at end of file diff --git a/applications/utilities/postProcessing/lagrangian/steadyParticleTracks/particleTrackDict b/applications/utilities/postProcessing/lagrangian/steadyParticleTracks/particleTrackDict new file mode 100644 index 0000000000..32c50127d5 --- /dev/null +++ b/applications/utilities/postProcessing/lagrangian/steadyParticleTracks/particleTrackDict @@ -0,0 +1,23 @@ +/*--------------------------------*- C++ -*----------------------------------*\ +| ========= | | +| \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | +| \\ / O peration | Version: dev | +| \\ / A nd | Web: www.OpenFOAM.org | +| \\/ M anipulation | | +\*---------------------------------------------------------------------------*/ +FoamFile +{ + version 2.0; + format ascii; + class dictionary; + location "constant"; + object particleTrackDict; +} +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +cloudName reactingCloud1Tracks; + +fields ( d U T ); + +// ************************************************************************* // + diff --git a/applications/utilities/postProcessing/lagrangian/steadyParticleTracks/steadyParticleTracks.C b/applications/utilities/postProcessing/lagrangian/steadyParticleTracks/steadyParticleTracks.C new file mode 100644 index 0000000000..d6a0352ae0 --- /dev/null +++ b/applications/utilities/postProcessing/lagrangian/steadyParticleTracks/steadyParticleTracks.C @@ -0,0 +1,326 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / 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 + steadyParticleTracks + +Description + Generates a VTK file of particle tracks for cases that were computed using + a steady-state cloud + NOTE: case must be re-constructed (if running in parallel) before use + +\*---------------------------------------------------------------------------*/ + +#include "argList.H" +#include "Cloud.H" +#include "IOdictionary.H" +#include "fvMesh.H" +#include "Time.H" +#include "timeSelector.H" +#include "OFstream.H" +#include "passiveParticleCloud.H" + +#include "SortableList.H" +#include "IOobjectList.H" +#include "PtrList.H" +#include "Field.H" +#include "steadyParticleTracksTemplates.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +using namespace Foam; + +namespace Foam +{ + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +label validateFields +( + const List& userFields, + const IOobjectList& cloudObjs +) +{ + List ok(userFields.size(), false); + + forAll(userFields, i) + { + ok[i] = ok[i] || fieldOk