# FireFoam

fireFoam

   Transient solver for fires and turbulent diffusion flames with reacting
particle clouds, surface film and pyrolysis modelling.


## 1 Solution Strategy

The source code can be found in fireFoam.C



/*---------------------------------------------------------------------------*\
=========                 |
\\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
\\    /   O peration     |
\\  /    A nd           | Copyright (C) 2011-2017 OpenFOAM Foundation
\\/     M anipulation  |
-------------------------------------------------------------------------------
This file is part of OpenFOAM.

OpenFOAM is free software: you can redistribute it and/or modify it
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 <http://www.gnu.org/licenses/>.

Application
fireFoam

Group
grpCombustionSolvers

Description
Transient solver for fires and turbulent diffusion flames with reacting
particle clouds, surface film and pyrolysis modelling.

\*---------------------------------------------------------------------------*/

#include "fvCFD.H"
#include "turbulentFluidThermoModel.H"
#include "basicReactingCloud.H"
#include "surfaceFilmModel.H"
#include "pyrolysisModelCollection.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[])
{
(
"Transient solver for fires and turbulent diffusion flames"
" with reacting particle clouds, surface film and pyrolysis modelling."
);

#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 "createRegionControls.H"

turbulence->validate();

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

Info<< "\nStarting time loop\n" << endl;

while (runTime.run())
{
#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();

runTime.printExecutionTime(Info);
}

Info<< "End" << endl;

return 0;
}

// ************************************************************************* //


## 2 Equations

### 2.1 Equations Liquid Film

The equation describing the evolution of the liquid film can be found in http://www.ilasseurope.org/ICLASS/iclass2012_Heidelberg/Contributions/Paper-pdfs/Contribution1294_b.pdf and https://www.witpress.com/Secure/elibrary/papers/MPF11/MPF11020FU1.pdf. A description of the method used to derive the equation can be found in https://www.researchgate.net/profile/Petr_Vita/publication/313655371_Thin_Film_Fluid_Flow_Simulation_on_Rotating_Discs/links/58a1a592a6fdccf5e9707665/Thin-Film-Fluid-Flow-Simulation-on-Rotating-Discs.pdf.

 $\frac{ \partial (\rho \delta )}{\partial t} + \frac{\partial}{\partial x_j} \left( \rho {U}_{j} \delta \right) = S_{\delta}$ (x)

$U$ represents the mean film velocity, $\delta$ the film height and $S_\delta$ mass source term resulting from the impingement of the liquid droplets on the film surface. The momentum conservation reads:

 $\frac{ \partial (\rho \delta U)}{\partial t} + \frac{\partial}{\partial x_j} \left( \rho U_i {U}_{j} \delta \right) = -\delta \frac{\partial p}{\partial x_i} + \tau + S_{U\delta}$ (x)

$p$ represents the pressure, $\tau$ the stress like components of the forces acting on the filme and $S_{U\delta}$ represents the contribution to the momentum by the impinging droplets. The pressure contribution is divided in three components: The contribution due to capillary effects $p_\sigma = -\sigma\frac{\partial \delta }{\partial x_i \partial x_i}$, the hydrostatic contribution $p_\delta = - \rho x_i g_i \delta$ and the contribution of the pressure originated by the surrounding gas $p_g$ .

 $\frac{ \partial (\rho \delta h )}{\partial t} + \frac{\partial}{\partial x_j} \left( \rho {U}_{j} \delta h \right) = S_{h\delta}$ (x)

$h$ represents the mean film enthalpy and $S_{h\delta}$ represents the contribution of the impinging droplets to the energy equation.

### 2.2 Equations Fluid

#### 2.2.1 Momentum conservation

The equation of motion are written for a moving frame of reference. They are however formulated for the absolute velocity (the derivation of the equations of motion can be found in https://openfoamwiki.net/index.php/See_the_MRF_development and also in https://diglib.tugraz.at/download.php?id=581303c7c91f9&location=browse. Some additional information can be found in https://pingpong.chalmers.se/public/pp/public_courses/course07056/published/1497955220499/resourceId/3711490/content/UploadedResources/HakanNilssonRotatingMachineryTrainingOFW11-1.pdf):

 $\frac{ \partial (\rho {u}_i)}{\partial t} + \frac{\partial}{\partial x_j} \left( \rho {u}_{rj} u_i \right) + \rho\epsilon_{ijk}\omega_i u_j= - \frac{\partial p_{rgh}} {\partial{x_i}} - \frac{\partial \rho g_j x_j}{\partial x_i} + \frac{\partial}{\partial x_j} \left( \tau_{ij} + \tau_{t_{ij}} \right) + F_{pi}$ (2)

$u$ represent the velocity, $u_r$ the relative veloicty, $g_i$ the gravitational acceleration, $p_{rgh} = p - \rho g_j x_j$ the pressure minus the hydrostatic pressure and $\tau_{ij}$ and $\tau_{t_{ij}}$ are the viscose and turbulent stresses. Note that since the relative velocity $u_r$ appears in the divergence term, the face flux $\phi$ appearing in the finite volume discretization of the momentum equation should be calculated with the relative velocity. $F_{pi}$ represents the force excreted by the particles on the fluid.

The source code can be found in Ueqn.H:



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
(
(