diff --git a/ReleaseNotes-dev b/ReleaseNotes-dev index 2418343bb7..7231b2b1d1 100644 --- a/ReleaseNotes-dev +++ b/ReleaseNotes-dev @@ -51,6 +51,18 @@ As 1.4792e-06; Ts 116; } +*** Lagrangian intermediate library + Extensively updated + *Updated* input format + Extended to include steady cloud tracking + *New* collision modelling + *Coupled* to new surface film modelling library + *New* sub-models + + NonSphereDrag: drag model to account for non-spherical particles + + ParticleTracks: post-processing model to generate track data, typically + during steady calculations + *Updated* sub-models + + Devolatilisation models: now act on a per-specie basis *** DSMC *** Dynamic Mesh @@ -106,6 +118,8 @@ *** *New* Solvers + =reactingParcelFilmFoam=: Lagrangian cloud and film transport in a reacting gas phase system + + =steadyReactingParcelFoam=: Steady solution of cloud and reacting systems + using local time stepping methods *** Modifications to multiphase and buoyant solvers + ... *** Modifications to solvers for sensible enthalpy @@ -149,6 +163,8 @@ field data. + =singleCellMesh=: new utility to convert mesh and fields to a single cell mesh. Great for postprocessing. + + =steadyParticleTracks=: Generates VTK tracks from the output of the cloud + =ParticleTracks= post-processing sub-model + Function objects: + =residualControl=: new function object to allow users to terminate steady state calculations when the defined residual levels are achieved diff --git a/applications/solvers/combustion/reactingFoam/chemistry.H b/applications/solvers/combustion/reactingFoam/chemistry.H index 8bfa872d7b..d7faf86b0c 100644 --- a/applications/solvers/combustion/reactingFoam/chemistry.H +++ b/applications/solvers/combustion/reactingFoam/chemistry.H @@ -11,7 +11,7 @@ if (turbulentReaction) { volScalarField tk = - Cmix*sqrt(turbulence->muEff()/rho/turbulence->epsilon()); + Cmix*sqrt(turbulence->muEff()/rho/turbulence->epsilon()); volScalarField tc = chemistry.tc(); // Chalmers PaSR model diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C index 36b3f1c3b0..0f9ff2f7f2 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionFoam.C @@ -31,7 +31,7 @@ Description \*---------------------------------------------------------------------------*/ #include "fvCFD.H" -#include "basicPsiThermo.H" +#include "basicRhoThermo.H" #include "turbulenceModel.H" #include "fixedGradientFvPatchFields.H" #include "regionProperties.H" @@ -121,7 +121,7 @@ int main(int argc, char *argv[]) << nl << endl; } - Info << "End\n" << endl; + Info<< "End\n" << endl; return 0; } diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/chtMultiRegionSimpleFoam.C b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/chtMultiRegionSimpleFoam.C index 17c4b5386c..0c590f781b 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/chtMultiRegionSimpleFoam.C +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/chtMultiRegionSimpleFoam.C @@ -30,7 +30,7 @@ Description \*---------------------------------------------------------------------------*/ #include "fvCFD.H" -#include "basicPsiThermo.H" +#include "basicRhoThermo.H" #include "turbulenceModel.H" #include "fixedGradientFvPatchFields.H" #include "regionProperties.H" diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/createFluidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/createFluidFields.H index 167484a7c7..75191acff9 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/createFluidFields.H +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/createFluidFields.H @@ -1,5 +1,5 @@ // Initialise fluid field pointer lists - PtrList thermoFluid(fluidRegions.size()); + PtrList thermoFluid(fluidRegions.size()); PtrList rhoFluid(fluidRegions.size()); PtrList KFluid(fluidRegions.size()); PtrList UFluid(fluidRegions.size()); @@ -28,7 +28,7 @@ thermoFluid.set ( i, - basicPsiThermo::New(fluidRegions[i]).ptr() + basicRhoThermo::New(fluidRegions[i]).ptr() ); Info<< " Adding to rhoFluid\n" << endl; diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H index b0cbbf80bf..2f339b87d7 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/pEqn.H @@ -1,5 +1,4 @@ { - // From buoyantSimpleFoam rho = thermo.rho(); rho = max(rho, rhoMin[i]); rho = min(rho, rhoMax[i]); @@ -13,6 +12,8 @@ phi = fvc::interpolate(rho)*(fvc::interpolate(U) & mesh.Sf()); bool closedVolume = adjustPhi(phi, U, p_rgh); + dimensionedScalar compressibility = fvc::domainIntegrate(psi); + bool compressible = (compressibility.value() > SMALL); surfaceScalarField buoyancyPhi = rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf(); phi -= buoyancyPhi; @@ -25,7 +26,11 @@ fvm::laplacian(rhorAUf, p_rgh) == fvc::div(phi) ); - p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell)); + p_rghEqn.setReference + ( + pRefCell, + compressible ? getRefCellValue(p_rgh, pRefCell) : pRefValue + ); p_rghEqn.solve(); @@ -50,10 +55,10 @@ // For closed-volume cases adjust the pressure level // to obey overall mass continuity - if (closedVolume) + if (closedVolume && compressible) { - p += (initialMass - fvc::domainIntegrate(psi*p)) - /fvc::domainIntegrate(psi); + p += (initialMass - fvc::domainIntegrate(thermo.rho())) + /compressibility; p_rgh = p - rho*gh; } diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/setRegionFluidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/setRegionFluidFields.H index cac7750e97..754f67c52c 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/setRegionFluidFields.H +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/chtMultiRegionSimpleFoam/fluid/setRegionFluidFields.H @@ -1,6 +1,6 @@ const fvMesh& mesh = fluidRegions[i]; - basicPsiThermo& thermo = thermoFluid[i]; + basicRhoThermo& thermo = thermoFluid[i]; volScalarField& rho = rhoFluid[i]; volScalarField& K = KFluid[i]; volVectorField& U = UFluid[i]; diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleCourantNo.C b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleCourantNo.C index 8895151817..606cd811ea 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleCourantNo.C +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/compressibleCourantNo.C @@ -34,16 +34,13 @@ Foam::scalar Foam::compressibleCourantNo const surfaceScalarField& phi ) { - scalar CoNum = 0.0; - scalar meanCoNum = 0.0; - scalarField sumPhi = fvc::surfaceSum(mag(phi))().internalField() /rho.internalField(); - CoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue(); + scalar CoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue(); - meanCoNum = + scalar meanCoNum = 0.5*(gSum(sumPhi)/gSum(mesh.V().field()))*runTime.deltaTValue(); Info<< "Region: " << mesh.name() << " Courant Number mean: " << meanCoNum diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H index b0a7f95912..012426af6a 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/createFluidFields.H @@ -1,5 +1,5 @@ // Initialise fluid field pointer lists - PtrList thermoFluid(fluidRegions.size()); + PtrList thermoFluid(fluidRegions.size()); PtrList rhoFluid(fluidRegions.size()); PtrList KFluid(fluidRegions.size()); PtrList UFluid(fluidRegions.size()); @@ -23,7 +23,7 @@ thermoFluid.set ( i, - basicPsiThermo::New(fluidRegions[i]).ptr() + basicRhoThermo::New(fluidRegions[i]).ptr() ); Info<< " Adding to rhoFluid\n" << endl; diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H index 686a8edcee..854dc9e2cc 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/pEqn.H @@ -1,5 +1,7 @@ { bool closedVolume = p_rgh.needReference(); + dimensionedScalar compressibility = fvc::domainIntegrate(psi); + bool compressible = (compressibility.value() > SMALL); rho = thermo.rho(); @@ -19,34 +21,48 @@ phi = phiU - rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf(); - for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) { - fvScalarMatrix p_rghEqn + fvScalarMatrix p_rghDDtEqn ( - fvm::ddt(psi, p_rgh) + fvc::ddt(psi, rho)*gh + fvc::ddt(rho) + psi*correction(fvm::ddt(p_rgh)) + fvc::div(phi) - - fvm::laplacian(rhorAUf, p_rgh) ); - p_rghEqn.solve - ( - mesh.solver + // Thermodynamic density needs to be updated by psi*d(p) after the + // pressure solution - done in 2 parts. Part 1: + thermo.rho() -= psi*p_rgh; + + for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++) + { + fvScalarMatrix p_rghEqn ( - p_rgh.select + p_rghDDtEqn + - fvm::laplacian(rhorAUf, p_rgh) + ); + + p_rghEqn.solve + ( + mesh.solver ( + p_rgh.select ( - oCorr == nOuterCorr-1 - && corr == nCorr-1 - && nonOrth == nNonOrthCorr + ( + oCorr == nOuterCorr-1 + && corr == nCorr-1 + && nonOrth == nNonOrthCorr + ) ) ) - ) - ); + ); - if (nonOrth == nNonOrthCorr) - { - phi += p_rghEqn.flux(); + if (nonOrth == nNonOrthCorr) + { + phi += p_rghEqn.flux(); + } } + + // Second part of thermodynamic density update + thermo.rho() += psi*p_rgh; } // Correct velocity field @@ -66,10 +82,10 @@ // For closed-volume cases adjust the pressure and density levels // to obey overall mass continuity - if (closedVolume) + if (closedVolume && compressible) { - p += (initialMass - fvc::domainIntegrate(psi*p)) - /fvc::domainIntegrate(psi); + p += (initialMass - fvc::domainIntegrate(thermo.rho())) + /compressibility; rho = thermo.rho(); p_rgh = p - rho*gh; } diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H index 89aaec4737..9cb2450952 100644 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H +++ b/applications/solvers/heatTransfer/chtMultiRegionFoam/fluid/setRegionFluidFields.H @@ -1,6 +1,6 @@ fvMesh& mesh = fluidRegions[i]; - basicPsiThermo& thermo = thermoFluid[i]; + basicRhoThermo& thermo = thermoFluid[i]; volScalarField& rho = rhoFluid[i]; volScalarField& K = KFluid[i]; volVectorField& U = UFluid[i]; diff --git a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/readSolidMultiRegionPISOControls.H b/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/readSolidMultiRegionPISOControls.H deleted file mode 100644 index 0c965a8322..0000000000 --- a/applications/solvers/heatTransfer/chtMultiRegionFoam/solid/readSolidMultiRegionPISOControls.H +++ /dev/null @@ -1,5 +0,0 @@ - const dictionary& piso = solidRegions[i].solutionDict().subDict("PISO"); - - const int nNonOrthCorr = - piso.lookupOrDefault("nNonOrthogonalCorrectors", 0); - 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/solvers/multiphase/interFoam/LTSInterFoam/setrDeltaT.H b/applications/solvers/multiphase/interFoam/LTSInterFoam/setrDeltaT.H index be8068c288..e84eef1023 100644 --- a/applications/solvers/multiphase/interFoam/LTSInterFoam/setrDeltaT.H +++ b/applications/solvers/multiphase/interFoam/LTSInterFoam/setrDeltaT.H @@ -103,7 +103,7 @@ fvc::sweep(rDeltaT, alpha1, nAlphaSweepIter, alphaSpreadDiff); } - Info<< "Flow time scale min/max = " + Info<< "Smoothed flow time scale min/max = " << gMin(1/rDeltaT.internalField()) << ", " << gMax(1/rDeltaT.internalField()) << endl; @@ -116,13 +116,12 @@ && runTime.timeIndex() > runTime.startTimeIndex() + 1 ) { - Info<< "Damping rDeltaT" << endl; rDeltaT = rDeltaT0*max(rDeltaT/rDeltaT0, 1.0 - rDeltaTDampingCoeff); - } - Info<< "Flow time scale min/max = " - << gMin(1/rDeltaT.internalField()) - << ", " << gMax(1/rDeltaT.internalField()) << endl; + Info<< "Damped flow time scale min/max = " + << gMin(1/rDeltaT.internalField()) + << ", " << gMax(1/rDeltaT.internalField()) << endl; + } label nAlphaSubCycles ( diff --git a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict index dbfcbd12d5..1a27956dbc 100644 --- a/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict +++ b/applications/utilities/mesh/generation/snappyHexMesh/snappyHexMeshDict @@ -1,7 +1,7 @@ /*--------------------------------*- C++ -*----------------------------------*\ | ========= | | | \\ / F ield | OpenFOAM: The Open Source CFD Toolbox | -| \\ / O peration | Version: 1.6 | +| \\ / O peration | Version: 1.7.1 | | \\ / A nd | Web: www.OpenFOAM.com | | \\/ M anipulation | | \*---------------------------------------------------------------------------*/ @@ -106,7 +106,8 @@ castellatedMeshControls // Explicit feature edge refinement // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ - // Specifies a level for any cell intersected by its edges. + // Specifies a level for any cell intersected by explicitly provided + // edges. // This is a featureEdgeMesh, read from constant/triSurface for now. features ( 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