An enthalpy-based pyrolysis model for charring and non-charring materials in case of fire

TL;DR: In this article, a simplified enthalpy-based pyrolysis model is proposed to simulate flame spread in a simulation of a developing fire, and the model is extended to multi-dimensional solid-phase treatments.

About: This article is published in Combustion and Flame.The article was published on 2010-04-01 and is currently open access. It has received 33 citations till now. The article focuses on the topics: Charring & Char.

Summary

2.1 Thermodynamics: Introduction

• The authors approach is largely based on [16] , one of the first theoretical papers on this topic.
• The approach is to consider five constituents: virgin solid material, char, volatiles, liquid water and water vapour.
• In [17, 18] , the same thermodynamic description for the enthalpy, in terms of temperature and local composition, is followed, in combination with Arrhenius type expressions to include pyrolysis process kinetics.
• Wide ranges of values are found in the literature for certain materials.

2.2 Thermodynamics: Relation between enthalpy and temperature

• The local relation between enthalpy and temperature is crucial.
• In principle, equations ( 1) and ( 2) completely define the relation between the specific enthalpy and the temperature, when the local composition of the material is known.
• This equation determines ΔQ pyr , but, as mentioned above, since the constitution of the virgin material, the char material and the pyrolysis gases is typically not exactly known, their formation enthalpies are unknown.
• Indeed, only two constants can be chosen freely, while the third one is determined by eq. ( 4).

a. T > T pyr

• This is only possible in charring materials.
• Under this condition, only char, pyrolysis gases and water vapour can be present.

2.3 Formation enthalpies: discussion

• In the previous section, the authors illustrated a direct link between the formation enthalpies of the different constituents and the 'heat of gasification' or 'heat of pyrolysis'.
• Indeed, starting from the basic expression (2), it becomes clear how ΔQ pyr depends on the temperature where it is defined.
• As mentioned before, one typically does not dispose of the knowledge, required to compute ΔQ pyr .
• This is equivalent to stating that only enthalpy differences need to be considered and the formation enthalpy of virgin material can be chosen equal to zero.
• The discussion is similar for the water vaporisation process, which is considered independent of the pyrolysis process.

2.4 Model description -Enthalpy equation

• In their model, the solid material is divided into a number of control volumes, which are kept fixed during the simulations (fixed computational mesh).
• Provided that a piecewise linear temperature field is used, this mesh need not be extremely fine and time steps in the solution procedure for the enthalpy equation need not be extremely small [15] .
• Discretisation issues and a sensitivity study are discussed below.
• For a fixed (sub-)volume 'V', the energy equation reads: EQUATION with '' q n ⋅ the heat flux out of the volume 'V' through its boundary S, as n is the unit normal vector, pointing outward with respect to the volume.
• In the solid material, though, the authors ignore the effect from static pressure.

2.4.1 Heat fluxes through the volume boundaries

• The right hand side of equation ( 11) consists of conduction and convection heat fluxes, i.e.: '' '' '' cond conv EQUATION.
• Conduction is modelled by Fourier's law: EQUATION.
• The value might depend on temperature and local composition at the cell face.
• In these pores, conduction in the gas phase can take place, possibly in combination with natural convection and, at sufficiently high temperatures, radiation.
• The convective fluxes due to transport of pyrolysis gases and water vapour (and possibly water liquid) are given by: EQUATION '' i m denotes the mass flux (kg/m 2 s) of constituent i leaving the volume, as determined below.

2.4.2 Motion of pyrolysis and evaporation front

• Thus, expression (15) determines the pyrolysis front motion.
• Note that, when there is no pyrolysis, e.g. in a cooling phase, expression (15) merely relates the temperature derivatives to the local thermal conductivities (change of material type over the pyrolysis front).
• At T = 373K, the following expression gives the motion of the evaporation front: EQUATION.
• On the other hand, v,wet means 'wet virgin', i.e. the solid virgin material, containing an amount of water.

2.5 Discussion: Relation with existing pyrolysis models

• The theoretical concepts, elaborated above from basic thermodynamic principles, have already been described in e.g. [16, 17, 24] .
• As such, this theoretical concept is not new.
• Note however, that multiple fronts (at different temperatures) can easily be introduced, so that the pyrolysis process could be modeled by means of several fronts.
• The heat flux onto the solid material is then obtained from the gas phase CFD calculations, while the present model gives mass flow rates and temperatures as boundary conditions to the CFD package.

3. Implementation and solution procedure

• As described in the previous section, the model considers enthalpy as the basic variable, for which a transport equation is solved.
• The latter is related to the motion of the pyrolysis front, which is assumed infinitely thin in the present model formulation.
• This is done, using the relationships between these variables, as described above.
• This is the simplest possible transport model.
• The vapour and volatiles are further assumed to take the local temperature of the solid material (local thermal equilibrium throughout the solid).

3.1 Solution procedure

• The flow chart in Figure 3 illustrates the solution procedure for dry charring materials.
• Starting from the initial conditions, physical time steps Stepping from time t n to t n+1 occurs in an iterative manner.
• In fact, only the central node is treated point-implicitly in the conductive fluxes in the subiterations described below.
• From the new enthalpy field, the temperature field, the position of the pyrolysis front and the pyrolysis front temperature must be reconstructed.
• From that point onwards, the authors proceed as just described.

3.3 Convective fluxes

• Equation (18) shows that the fluxes through the cell faces determine the update in enthalpy of the cell.
• Besides heat transfer by means of conduction, energy is also transported with the movement of water vapour and volatiles out of the solid material.
• The convective flux through cell face i-1/2 reads: EQUATION ).
• The mass flow rates are determined as: EQUATION ( ) EQUATION EQUATION EQUATION.
• Also note that the authors only consider uni-directional flow of the pyrolysis gases towards the side where the external heat flux is imposed.

3.4 Reconstruction of temperature field and fronts' position

• From equation (18) , the enthalpy update can be calculated.
• For the construction of the fluxes, the knowledge of the temperature and fronts' position is needed, though.
• Therefore, from the updated enthalpy values, the temperature field and fronts' position must be reconstructed.
• Two constraints determine the relation between these variables: the enthalpy is a function of temperature and local composition and the motion of the front correlates the temperature gradients on both sides of the fronts.

3.4.1 Constraint 1: enthalpy as a function of temperature and composition

• The authors introduce a function F enth,i , for each computational cell, to express the relationship between the enthalpy value and the temperature.
• As the authors assume here that water vapour and volatiles leave the solid material as soon as they are formed, the mass fractions α w,v and α g are zero.
• The authors consider the value for enthalpy as the averaged value over the computational cell.
• As such, a few possible configurations can be distinguished, for which different expressions can be formulated, all derived from the general form: EQUATION EQUATION ).
• The authors distinguish between the following possible cell configurations.

3.4.2 Constraint 2: motion of the front

• For the cells containing a front, the authors require an extra constraint.
• For non-charring materials, ρ c equals zero and the conduction term in the char is replaced by the external heat flux.
• The authors use the notation k v,dry here for clarity.

3.4.3 Inversion of the constraints to determine temperature and fronts' position.

• Using expressions ( 23) -( 25), for a given enthalpy field, the corresponding temperature field and front positions can be found.
• Since an iterative procedure is adopted anyway, this coupling is not taken into account for the inversion of the system and only the temperature of cell i is considered as an unknown in the functions.
• The authors focus now on a cell that contains the pyrolysis front.
• The notation means that the quantities behind the semi-colon are 'known', whereas the variables ahead of the semi-colon are to be computed.
• The convergence check has been reported above, in section 3.1.

3.5 Treatment of the boundaries

• It is clear that the authors cannot use expressions (24) , since these functions return zero by construction of the temperature profile.
• Hence, the alternative constraint at the boundaries.

3.6 Discussion: zero-th order temperature field representation

• In [25] , a zero-th order representation is adopted for the temperature field, i.e. the temperature is uniform in each of the computational cells.
• The essential difference to the model formulation as presented above is that, when the mushy cell is pyrolysing, its temperature is kept fixed, equal to the pyrolysis temperature.
• This has serious consequences on the evolution of the pyrolysis gases mass flow rates and pyrolysis front motion: when the mushy cell has just become pure char, the next cell to pyrolyse must first still heat up to T pyr and during this period, '' pyr m drops to zero, which is unacceptable.
• This was already recognised in [8] , but the problem was not really solved there.
• A dual mesh technique was introduced, effectively reducing the mentioned undesired phenomenon, but not solving the problem.

4. Discussion of results

• The authors restrict ourselves to configurations where the externally imposed heat flux is not computed from flame radiation, in order to avoid related uncertainty.
• Unless stated otherwise, all results are obtained with the piecewise linear temperature field representation.

4.1 One-dimensional configuration -charring material

• Consider a one-dimensional configuration, with an external radiative heat flux imposed at one side.
• As mentioned in [5] , the peak mass flow rate corresponds very well to the experimentally reported value.
• The authors now discuss the dependence of the results on the number of computational cells and the physical time step .
• The remainder of the pyrolysis process takes place in dry virgin material, so that very similar mass flow rate profiles are observed.
• With the increase in moisture content in the solid, the inner solid temperature rises more slowly.

4.2 One-dimensional configuration -Non-charring material

• The front surface of the virgin material, which is now also the pyrolysis front, is exposed to the external heat flux.
• Agreement with the data of [21] is very good when the heat losses are accounted for.
• After some time, the mass loss rate attains a quasi-steady state: the incoming heat flux provides energy for the pyrolysis process and for heat conduction into the virgin solid.
• Obviously, this is more pronounced when there is no heat loss at the exposed surface.
• All these issues are well captured with the present model.

4.3 Multi-dimensional configuration -Upward flame spread

• Extensions to more dimensions can be done with different degrees of complexity.
• As a test, consider a vertically oriented sample with the same material properties as described in section 4.1.
• Figure 13 , showing front surface temperature evolutions in time, supports these findings.
• Note that the equilibrium end temperature at y = 5cm and y = 9cm almost exactly matches the equilibrium end temperature: the configuration is practically onedimensional, with negligible net effect of conductive fluxes in the upward direction.

5. Summary and conclusions

• Starting from a basic thermodynamic description of pyrolysis phenomena, a simplified The major assumption is that pyrolysis and evaporation are isothermal processes, taking place at infinitely thin fronts.
• The model has then been applied to dry and wet charring materials, to non-charring materials and to a multi-dimensional configuration, resembling upward flame spread.
• Good agreement was illustrated for mass flow rates of pyrolysis gases in dry charring and non-charring materials (to which the model is directly applicable).
• The model formulation is robust with respect to several numerical aspects: the dependence of the results on the computational mesh cells' size and the physical time step size is small.
• Finally, the authors illustrated that the model can deal with multi-dimensional configurations by means of a test case, resembling upward flame spread.

Figures

Figure 1 summarises the expressions above for charring and non-charring materials. For ease of drawing, constant thermal capacities are assumed. With expressions (1), (6) or (7) and (9), we can now define the local relations between the specific enthalpy and the temperature for the different temperature ranges:

Figure 13. Surface temperature evolution in time at y = 0mm, y = 1cm, y = 5cm and y = 9cm. Top: thickness = 3mm; middle: thickness = 5mm; bottom: thickness = 1cm. 1D: onedimensional results for comparison purposes.

Figure 7: Evolution in time of pyrolysis gases mass flow rate (left) and surface temperature (right). Dashed line: piecewise linear temperature representation; solid line: zero-th order temperature representation.

Figure 10. Mass loss rate evolution in time (non-charring material).

Figure 9. Temperature evolution in time (Depths: surface, 5mm and 10mm). Solid line: dry material; dashed line: 10% moisture content. Squares: experimental data [3].

Figure 11. Temperature (top) and char fraction (bottom) fields after 10s, 50s, 100s and 200s. Material thickness equal to 5mm.

Figure 3: Flow chart: Solution procedure for dry charring materials (mushy cell).

Figure 2. Energy balance at the pyrolysis front.

Figure 4: Piecewise linear temperature distribution.

Figure 12. Left: Mass flow rate evolution in time at y = 0mm, y = 1cm, y = 5cm and y = 9cm. Middle: evolution of yp and yf; Right: early stages of the evolution of yp and yf. Top: thickness = 3mm; middle: thickness = 5mm; bottom: thickness = 10mm. 1D: one-dimensional results for comparison purposes.

Figure 5. Left: Pyrolysis gases mass flow rate evolution in time. Right: Temperature evolution in time (Depths: surface, 5mm and 10mm). Solid line: numerical simulation. Symbols: squares: experimental data [3]; triangles: numerical results [2].

Figure 6. Left: Influence of the number of cells on the evolution of pyrolysis gases mass flow rate in time (physical time step equal to 0.5s). Right: Influence of the physical time step size on the evolution of pyrolysis gases mass flow rate in time (number of cells equal to 40).

Figure 1 summarises the expressions above for charring and non-charring materials. For ease of drawing, constant thermal capacities are assumed. With expressions (1), (6) or (7) and (9), we can now define the local relations between the specific enthalpy and the temperature for the different temperature ranges:

Frequently asked questions

Q1. What are the contributions in "An enthalpy-based pyrolysis model for charring and non-charring materials in case of fire wasan," ?

In this paper, first, the basic thermodynamic description of pyrolysis phenomena is revisited for charring and non-charring materials, possibly containing moisture. Next, numerical issues and implementation are discussed. The authors consider basic test cases with imposed external heat flux to a solid material that can be dry or contain moisture. The authors illustrate that continuous pyrolysis gases mass flow rates are obtained when a piecewise linear representation of the temperature field is adopted on the fixed computational mesh. The authors show that the present model formulation is robust with respect to numerical aspects ( cell size and time step ) and that the model performs well for variable external heat fluxes. By means of a numerical test case, the authors show that the model, when coupled to CFD calculations, is able to simulate flame spread. 

Q2. What is the effect of adding a constant value to the external heat flux?

A constant value of 25kW/m2 is added tothe external heat flux in the region y < yf(t), resembling radiative heat feedback from flames. 

Q3. What is the effect of the positive feedback loop on the pyrolysis front?

As yp approaches the top surface, an acceleration is observed again, as there are no conductive heat losses at the top surface and thus, due to net incoming conductive heat fluxes from below, the pyrolysis front moves more rapidly than in a one-dimensional configuration. 

Q4. What is the mass loss rate after a while?

After some time, the mass loss rate attains a quasi-steady state: the incoming heat flux provides energy for the pyrolysis process and for heat conduction into the virgin solid. 

Q5. What is the calculation of the pyrolysis?

The calculation then consists of solving a conduction problem, with an incoming heat flux and a moving boundary as the pyrolysis takes place. 

Q6. What is the effect of the evaporation front on the mass flow rate?

A sudden drop in the total mass flow rate is observed when the evaporation front reaches the back surface, because the evaporation front reaches this surface with a non-zero velocity. 

Q7. how do the authors calculate the conductive fluxes of cell faces?

The conductive fluxes of cell faces are calculated using Fourier’s law:'' , 1 / 2 1 / 2 / , dT cond i i dx r l i q k± ±= − (19)As mentioned, the local material properties are used: 1 / 2i ck k± = if face i±1/2 is in charmaterial and 1 / 2 , ,i v v w l w lk k kα α± = + if it is in virgin material, with αi the local massfraction of constituent i. 

Q8. What is the role of simulations in pyrolysis?

In numerical simulations, this implies coupling of gas phase CFD (‘Computational Fluid Dynamics’) simulations, including turbulent combustion and radiation, to pyrolysis simulations in the solid material. 

Q9. What is the appealing feature of the present model?

In fact, it is a particularly appealing feature of the present model that the equations are solved on a fixed computational mesh. 

Q10. What is the mass flow rate in the mushy cell?

As explained above, the mass flow rate (21) inevitably always drops to zero when the char fraction in the mushy cell becomes equal to 1. 

Q11. How can the model deal with multi-dimensional configurations?

the authors illustrated that the model can deal with multi-dimensional configurations by means of a test case, resembling upward flame spread. 

Q12. How do the authors calculate the flame height of a material?

As soon as pyrolysis starts, the authors use a correlation for upward flame spread [27] to calculate the flame height yf(t) from the pyrolysis height yp(t) (i.e. the height over which the material has pyrolysed at its front surface or, alternatively, where the front surface temperature exceeds Tpyr): ( ) ( ) ( ) 2 3' '0.0433f p b my t y t Q Q= + + . 

Q13. What is the boundary condition at the front surface?

The back surface is perfectly insulated and the boundary condition at the front surface reads [21]:( )'' 4 2; 46 ( / )f net s amb sx x q h T T T kW mεσ= = − − − .