diff --git a/Allwmake b/Allwmake index 6d6eb7932..574729e3c 100755 --- a/Allwmake +++ b/Allwmake @@ -24,5 +24,6 @@ wmake applications/solvers/df0DFoam wmake applications/solvers/dfLowMachFoam wmake applications/solvers/dfHighSpeedFoam wmake applications/solvers/dfSprayFoam +wmake applications/solvers/dfBuoyancyFoam wmake applications/utilities/flameSpeed diff --git a/applications/solvers/dfBuoyancyFoam/Make/files b/applications/solvers/dfBuoyancyFoam/Make/files new file mode 100644 index 000000000..179d90a27 --- /dev/null +++ b/applications/solvers/dfBuoyancyFoam/Make/files @@ -0,0 +1,3 @@ +dfBuoyancyFoam.C + +EXE = $(DF_APPBIN)/dfBuoyancyFoam diff --git a/applications/solvers/dfBuoyancyFoam/Make/options b/applications/solvers/dfBuoyancyFoam/Make/options new file mode 100644 index 000000000..24a198dc9 --- /dev/null +++ b/applications/solvers/dfBuoyancyFoam/Make/options @@ -0,0 +1,76 @@ +-include $(GENERAL_RULES)/mplibType + +EXE_INC = -std=c++14 \ + -g \ + -fopenmp \ + -Wno-unused-variable \ + -Wno-unused-but-set-variable \ + -Wno-old-style-cast \ + -I. \ + $(PFLAGS) $(PINC) \ + -I$(LIB_SRC)/finiteVolume/lnInclude \ + -I$(LIB_SRC)/meshTools/lnInclude \ + -I$(LIB_SRC)/sampling/lnInclude \ + -I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \ + -I$(LIB_SRC)/TurbulenceModels/compressible/lnInclude \ + -I$(LIB_SRC)/transportModels/compressible/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/basic/lnInclude \ + -I$(LIB_SRC)/regionModels/regionModel/lnInclude \ + -I$(LIB_SRC)/regionModels/surfaceFilmModels/lnInclude \ + -I$(LIB_SRC)/regionModels/pyrolysisModels/lnInclude \ + -I$(LIB_SRC)/lagrangian/basic/lnInclude \ + -I$(DF_SRC)/lagrangian/intermediate/lnInclude \ + -I$(LIB_SRC)/lagrangian/intermediate/lnInclude \ + -I$(DF_SRC)/lagrangian/spray/lnInclude \ + -I$(LIB_SRC)/lagrangian/spray/lnInclude \ + -I$(LIB_SRC)/lagrangian/distributionModels/lnInclude \ + -I$(LIB_SRC)/ODE/lnInclude \ + -I$(DF_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \ + -I$(LIB_SRC)/thermophysicalModels/thermophysicalProperties/lnInclude \ + -I$(DF_SRC)/thermophysicalModels/SLGThermo/lnInclude \ + -I$(LIB_SRC)/Pstream/mpi \ + -I$(DF_SRC)/dfCanteraMixture/lnInclude \ + -I$(DF_SRC)/dfChemistryModel/lnInclude \ + -I$(DF_SRC)/dfCombustionModels/lnInclude \ + -I$(CANTERA_ROOT)/include \ + $(if $(LIBTORCH_ROOT),-I$(LIBTORCH_ROOT)/include,) \ + $(if $(LIBTORCH_ROOT),-I$(LIBTORCH_ROOT)/include/torch/csrc/api/include,) \ + $(PYTHON_INC_DIR) \ + $(if $(AMGX_DIR), -I$(DF_ROOT)/src_gpu,) \ + $(if $(AMGX_DIR), -I/usr/local/cuda/include,) \ + $(if $(AMGX_DIR), -I$(AMGX_DIR)/include,) \ + $(if $(ODE_GPU_SOLVER), -I$(OPENCC_PATH)/include,) \ + $(if $(ODE_GPU_SOLVER), -DODE_GPU_SOLVER,) + +EXE_LIBS = \ + -lfiniteVolume \ + -lmeshTools \ + -lsampling \ + -lturbulenceModels \ + -lcompressibleTransportModels \ + -llagrangian \ + -lregionModels \ + -L$(DF_LIBBIN) \ + -ldfSurfaceFilmModels \ + -ldfSLGThermo \ + -ldfLagrangianIntermediate \ + -ldfLagrangianTurbulence \ + -ldfLagrangianSpray \ + -ldfFluidThermophysicalModels \ + -ldfCompressibleTurbulenceModels \ + -ldfThermophysicalProperties \ + -ldfCanteraMixture \ + -ldfChemistryModel \ + -ldfCombustionModels \ + $(CANTERA_ROOT)/lib/libcantera.so \ + $(if $(LIBTORCH_ROOT),$(LIBTORCH_ROOT)/lib/libtorch.so,) \ + $(if $(LIBTORCH_ROOT),$(LIBTORCH_ROOT)/lib/libc10.so,) \ + $(if $(LIBTORCH_ROOT),-rdynamic,) \ + $(if $(LIBTORCH_ROOT),-lpthread,) \ + $(if $(LIBTORCH_ROOT),$(DF_SRC)/dfChemistryModel/DNNInferencer/build/libDNNInferencer.so,) \ + $(if $(PYTHON_LIB_DIR),$(PYTHON_LIB_DIR),) \ + $(if $(AMGX_DIR), /usr/local/cuda/lib64/libcudart.so,) \ + $(if $(AMGX_DIR), /usr/local/cuda/lib64/libnccl.so,) \ + $(if $(AMGX_DIR), $(DF_ROOT)/src_gpu/build/libdfMatrix.so,) \ + $(if $(AMGX_DIR), $(AMGX_DIR)/build/libamgxsh.so,) \ + $(if $(ODE_GPU_SOLVER), $(ODE_GPU_SOLVER)/lib/libopencc.so,) diff --git a/applications/solvers/dfBuoyancyFoam/UEqn.H b/applications/solvers/dfBuoyancyFoam/UEqn.H new file mode 100644 index 000000000..5314a0ff9 --- /dev/null +++ b/applications/solvers/dfBuoyancyFoam/UEqn.H @@ -0,0 +1,34 @@ + MRF.correctBoundaryVelocity(U); + + fvVectorMatrix UEqn + ( + fvm::ddt(rho, U) + fvm::div(phi, U) + + MRF.DDt(rho, U) + + turbulence->divDevRhoReff(U) + == + parcels.SU(U) + // + fvOptions(rho, U) + ); + + UEqn.relax(); + + // fvOptions.constrain(UEqn); + + if (pimple.momentumPredictor()) + { + solve + ( + UEqn + == + fvc::reconstruct + ( + ( + - ghf*fvc::snGrad(rho) + - fvc::snGrad(p_rgh) + )*mesh.magSf() + ) + ); + + // fvOptions.correct(U); + K = 0.5*magSqr(U); + } diff --git a/applications/solvers/dfBuoyancyFoam/YEEqn.H b/applications/solvers/dfBuoyancyFoam/YEEqn.H new file mode 100644 index 000000000..245aeda0d --- /dev/null +++ b/applications/solvers/dfBuoyancyFoam/YEEqn.H @@ -0,0 +1,90 @@ +tmp> mvConvection +( + fv::convectionScheme::New + ( + mesh, + fields, + phi, + mesh.divScheme("div(phi,Yi_h)") + ) +); +{ + // radiation->correct(); + combustion->correct(); + + HRR = combustion->Qdot(); + + volScalarField Yt(0.0*Y[0]); + + forAll(Y, i) + { + if (i != inertIndex) + { + volScalarField& Yi = Y[i]; + + fvScalarMatrix YiEqn + ( + fvm::ddt(rho, Yi) + + mvConvection->fvmDiv(phi, Yi) + - fvm::laplacian(turbulence->alphaEff(), Yi) + == + parcels.SYi(i, Yi) + // + surfaceFilm.Srho(i) + + combustion->R(Yi) + // + fvOptions(rho, Yi) + ); + YiEqn.relax(); + + // fvOptions.constrain(YiEqn); + + YiEqn.solve("Yi"); + + // fvOptions.correct(Yi); + + Yi.max(0.0); + Yt += Yi; + } + } + + Y[inertIndex] = scalar(1) - Yt; + Y[inertIndex].max(0.0); + + volScalarField& he = thermo.he(); + + fvScalarMatrix EEqn + ( + fvm::ddt(rho, he) + mvConvection->fvmDiv(phi, he) + + fvc::ddt(rho, K) + fvc::div(phi, K) + + ( + he.name() == "e" + ? fvc::div + ( + fvc::absolute(phi/fvc::interpolate(rho), U), + p, + "div(phiv,p)" + ) + : -dpdt + ) + - fvm::laplacian(turbulence->alphaEff(), he) + == + // combustion->Qdot() // ha Eqn has no reaction source + // + radiation->Sh(thermo, he) + parcels.Sh(he) + // + surfaceFilm.Sh() + // + fvOptions(rho, he) + ); + + EEqn.relax(); + + // fvOptions.constrain(EEqn); + + EEqn.solve(); + + // fvOptions.correct(he); + + // thermo.correct(); + chemistry->correctThermo(); + + Info<< "min/max(T) = " + << min(T).value() << ", " << max(T).value() << endl; +} diff --git a/applications/solvers/dfBuoyancyFoam/createClouds.H b/applications/solvers/dfBuoyancyFoam/createClouds.H new file mode 100644 index 000000000..c568be12a --- /dev/null +++ b/applications/solvers/dfBuoyancyFoam/createClouds.H @@ -0,0 +1,9 @@ +Info<< "\nConstructing reacting cloud" << endl; +basicReactingCloud parcels +( + "reactingCloud1", + rho, + U, + g, + slgThermo +); diff --git a/applications/solvers/dfBuoyancyFoam/createFieldRefs.H b/applications/solvers/dfBuoyancyFoam/createFieldRefs.H new file mode 100644 index 000000000..6e9886e92 --- /dev/null +++ b/applications/solvers/dfBuoyancyFoam/createFieldRefs.H @@ -0,0 +1,5 @@ +const volScalarField& psi = thermo.psi(); +const volScalarField& T = thermo.T(); +// regionModels::surfaceFilmModel& surfaceFilm = tsurfaceFilm(); +const word inertSpecie(chemistry->lookup("inertSpecie")); +const label inertIndex(chemistry->species()[inertSpecie]); diff --git a/applications/solvers/dfBuoyancyFoam/createFields.H b/applications/solvers/dfBuoyancyFoam/createFields.H new file mode 100644 index 000000000..6b9464c8f --- /dev/null +++ b/applications/solvers/dfBuoyancyFoam/createFields.H @@ -0,0 +1,168 @@ +Info<< "Reading thermophysical properties\n" << endl; +//autoPtr pThermo(psiReactionThermo::New(mesh)); +//psiReactionThermo& thermo = pThermo(); +//thermo.validate(args.executable(), "h", "e"); + +CanteraMixture::setEnergyName("ha"); +// fluidThermo* pThermo = new heRhoThermo(mesh, word::null); +fluidThermo* pThermo = new heRhoThermo(mesh, word::null); +fluidThermo& thermo = *pThermo; +SLGThermo slgThermo(mesh, thermo); + +// basicSpecieMixture& composition = thermo.composition(); +// PtrList& Y = composition.Y(); + +// const word inertSpecie(thermo.lookup("inertSpecie")); +// if (!composition.species().found(inertSpecie)) +// { +// FatalIOErrorIn(args.executable().c_str(), thermo) +// << "Inert specie " << inertSpecie << " not found in available species " +// << composition.species() +// << exit(FatalIOError); +// } + +Info<< "Creating field rho\n" << endl; +volScalarField rho +( + IOobject + ( + "rho", + runTime.timeName(), + mesh, + IOobject::NO_READ, + IOobject::AUTO_WRITE + ), + thermo.rho() +); + +volScalarField& p = thermo.p(); + +Info<< "\nReading field U\n" << endl; +volVectorField U +( + IOobject + ( + "U", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh +); + +#include "compressibleCreatePhi.H" + +#include "createMRF.H" + + +Info<< "Creating turbulence model\n" << endl; +autoPtr turbulence +( + compressible::turbulenceModel::New + ( + rho, + U, + phi, + thermo + ) +); + +Info<< "Creating combustion model\n" << endl; +autoPtr> combustion +( + CombustionModel::New(thermo, turbulence()) +); + +dfChemistryModel* chemistry = combustion->chemistry(); +PtrList& Y = chemistry->Y(); + +#include "readGravitationalAcceleration.H" +#include "readhRef.H" +#include "gh.H" +#include "readpRef.H" + +volScalarField p_rgh +( + IOobject + ( + "p_rgh", + runTime.timeName(), + mesh, + IOobject::MUST_READ, + IOobject::AUTO_WRITE + ), + mesh +); + +mesh.setFluxRequired(p_rgh.name()); + +#include "phrghEqn.H" + + +multivariateSurfaceInterpolationScheme::fieldTable fields; + +forAll(Y, i) +{ + fields.add(Y[i]); +} +fields.add(thermo.he()); + +volScalarField HRR +( + IOobject + ( + "HRR", + runTime.timeName(), + mesh, + IOobject::READ_IF_PRESENT, + IOobject::AUTO_WRITE + ), + mesh, + dimensionedScalar(dimEnergy/dimVolume/dimTime, 0) +); + +IOdictionary additionalControlsDict +( + IOobject + ( + "additionalControls", + runTime.constant(), + mesh, + IOobject::MUST_READ_IF_MODIFIED, + IOobject::NO_WRITE + ) +); + +Switch solvePrimaryRegion +( + additionalControlsDict.lookup("solvePrimaryRegion") +); + +Switch solvePyrolysisRegion +( + additionalControlsDict.lookupOrDefault("solvePyrolysisRegion", true) +); + + +Info<< "Creating field dpdt\n" << endl; +volScalarField dpdt +( + IOobject + ( + "dpdt", + runTime.timeName(), + mesh + ), + mesh, + dimensionedScalar(p.dimensions()/dimTime, 0) +); + +Info<< "Creating field kinetic energy K\n" << endl; +volScalarField K("K", 0.5*magSqr(U)); + +#include "createClouds.H" +#include "createSurfaceFilmModel.H" +// #include "createPyrolysisModel.H" +// #include "createRadiationModel.H" +// #include "createFvOptions.H" diff --git a/applications/solvers/dfBuoyancyFoam/createPyrolysisModel.H b/applications/solvers/dfBuoyancyFoam/createPyrolysisModel.H new file mode 100644 index 000000000..4c93ad2c6 --- /dev/null +++ b/applications/solvers/dfBuoyancyFoam/createPyrolysisModel.H @@ -0,0 +1,3 @@ +Info<< "Creating pyrolysis model" << endl; + +regionModels::pyrolysisModels::pyrolysisModelCollection pyrolysis(mesh); diff --git a/applications/solvers/dfBuoyancyFoam/createSurfaceFilmModel.H b/applications/solvers/dfBuoyancyFoam/createSurfaceFilmModel.H new file mode 100644 index 000000000..81995c09a --- /dev/null +++ b/applications/solvers/dfBuoyancyFoam/createSurfaceFilmModel.H @@ -0,0 +1,6 @@ +Info<< "\nConstructing surface film model" << endl; + +autoPtr tsurfaceFilm +( + regionModels::surfaceFilmModel::New(mesh, g) +); diff --git a/applications/solvers/dfBuoyancyFoam/dfBuoyancyFoam.C b/applications/solvers/dfBuoyancyFoam/dfBuoyancyFoam.C new file mode 100644 index 000000000..b32a006b6 --- /dev/null +++ b/applications/solvers/dfBuoyancyFoam/dfBuoyancyFoam.C @@ -0,0 +1,136 @@ +/*---------------------------------------------------------------------------*\ + ========= | + \\ / F ield | OpenFOAM: The Open Source CFD Toolbox + \\ / O peration | Website: https://openfoam.org + \\ / A nd | Copyright (C) 2011-2018 OpenFOAM Foundation + \\/ M anipulation | +------------------------------------------------------------------------------- +License + This file is part of OpenFOAM. + + OpenFOAM is free software: you can redistribute it and/or modify it + under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OpenFOAM is distributed in the hope that it will be useful, but WITHOUT + ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or + FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License + for more details. + + You should have received a copy of the GNU General Public License + along with OpenFOAM. If not, see . + +Application + fireFoam + +Description + Transient solver for fires and turbulent diffusion flames with reacting + particle clouds, surface film and pyrolysis modelling. + +\*---------------------------------------------------------------------------*/ + +#include "stdlib.h" +#include "dfChemistryModel.H" +#include "CanteraMixture.H" +#include "dfSingleStepReactingMixture.H" +// #include "hePsiThermo.H" +#include "heRhoThermo.H" + +#include "fvCFD.H" +#include "turbulentFluidThermoModel.H" +#include "basicReactingCloud.H" +#include "surfaceFilmModel.H" +#include "pyrolysisModelCollection.H" +// #include "radiationModel.H" +#include "SLGThermo.H" +// #include "solidChemistryModel.H" +// #include "psiReactionThermo.H" +#include "CombustionModel.H" +#include "pimpleControl.H" +#include "fvOptions.H" + +// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + +int main(int argc, char *argv[]) +{ + #include "postProcess.H" + + #include "setRootCaseLists.H" + #include "createTime.H" + #include "createMesh.H" + #include "createControl.H" + #include "createFields.H" + #include "createFieldRefs.H" + #include "initContinuityErrs.H" + #include "createTimeControls.H" + #include "compressibleCourantNo.H" + #include "setInitialDeltaT.H" + // #include "readPyrolysisTimeControls.H" + + turbulence->validate(); + + // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // + + Info<< "\nStarting time loop\n" << endl; + + while (runTime.run()) + { + #include "readTimeControls.H" + #include "compressibleCourantNo.H" + // #include "solidRegionDiffusionNo.H" + #include "setMultiRegionDeltaT.H" + #include "setDeltaT.H" + + runTime++; + + Info<< "Time = " << runTime.timeName() << nl << endl; + + parcels.evolve(); + + // surfaceFilm.evolve(); + + // if(solvePyrolysisRegion) + // { + // pyrolysis.evolve(); + // } + + if (solvePrimaryRegion) + { + #include "rhoEqn.H" + + // --- PIMPLE loop + while (pimple.loop()) + { + #include "UEqn.H" + #include "YEEqn.H" + + // --- Pressure corrector loop + while (pimple.correct()) + { + #include "pEqn.H" + } + + if (pimple.turbCorr()) + { + turbulence->correct(); + } + } + + rho = thermo.rho(); + } + + runTime.write(); + + Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s" + << " ClockTime = " << runTime.elapsedClockTime() << " s" + << nl << endl; + } + + Info<< "End" << endl; + + return 0; +} + + +// ************************************************************************* // diff --git a/applications/solvers/dfBuoyancyFoam/pEqn.H b/applications/solvers/dfBuoyancyFoam/pEqn.H new file mode 100644 index 000000000..38db9b17e --- /dev/null +++ b/applications/solvers/dfBuoyancyFoam/pEqn.H @@ -0,0 +1,60 @@ +rho = thermo.rho(); + +volScalarField rAU(1.0/UEqn.A()); +surfaceScalarField rhorAUf("rhorAUf", fvc::interpolate(rho*rAU)); +volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p)); + +surfaceScalarField phig("phig", -rhorAUf*ghf*fvc::snGrad(rho)*mesh.magSf()); + +surfaceScalarField phiHbyA +( + "phiHbyA", + ( + fvc::flux(rho*HbyA) + + MRF.zeroFilter(rhorAUf*fvc::ddtCorr(rho, U, phi)) + ) + + phig +); + +MRF.makeRelative(fvc::interpolate(rho), phiHbyA); + +// Update the pressure BCs to ensure flux consistency +constrainPressure(p_rgh, rho, U, phiHbyA, rhorAUf, MRF); + +while (pimple.correctNonOrthogonal()) +{ + fvScalarMatrix p_rghEqn + ( + fvm::ddt(psi, p_rgh) + + fvc::ddt(psi, rho)*gh + + fvc::ddt(psi)*pRef + + fvc::div(phiHbyA) + - fvm::laplacian(rhorAUf, p_rgh) + == + parcels.Srho() + // + surfaceFilm.Srho() + // + fvOptions(psi, p_rgh, rho.name()) + ); + + p_rghEqn.solve(); + + if (pimple.finalNonOrthogonalIter()) + { + phi = phiHbyA + p_rghEqn.flux(); + U = HbyA + rAU*fvc::reconstruct((p_rghEqn.flux() + phig)/rhorAUf); + U.correctBoundaryConditions(); + // fvOptions.correct(U); + } +} + +p = p_rgh + rho*gh + pRef; + +#include "rhoEqn.H" +#include "compressibleContinuityErrs.H" + +K = 0.5*magSqr(U); + +if (thermo.dpdt()) +{ + dpdt = fvc::ddt(p); +} diff --git a/applications/solvers/dfBuoyancyFoam/phrghEqn.H b/applications/solvers/dfBuoyancyFoam/phrghEqn.H new file mode 100644 index 000000000..3c5175f41 --- /dev/null +++ b/applications/solvers/dfBuoyancyFoam/phrghEqn.H @@ -0,0 +1,62 @@ +if (pimple.dict().lookupOrDefault("hydrostaticInitialization", false)) +{ + volScalarField& ph_rgh = regIOobject::store + ( + new volScalarField + ( + IOobject + ( + "ph_rgh", + "0", + mesh, + IOobject::MUST_READ, + IOobject::NO_WRITE + ), + mesh + ) + ); + + if (equal(runTime.value(), 0)) + { + p = ph_rgh + rho*gh + pRef; + thermo.correct(); + rho = thermo.rho(); + + label nCorr + ( + pimple.dict().lookupOrDefault