- Email: [email protected]
Critical evaluation of an integral model for the pyrolysis of charring materials
ARTICLE IN PRESS
Fire Safety Journal 40 (2005) 121–140 www.elsevier.com/locate/firesaf
Critical evaluation of an integral model for the pyrolysis of charring materials E. Theuns, B. Merci, J. Vierendeels, P. Vandevelde Department of Flow, Heat and Combustion Mechanics, Ghent University, Gent, Belgium Received 2 February 2004; received in revised form 2 August 2004; accepted 27 September 2004 Available online 2 December 2004
Abstract When field models are used to predict fires and flame spread, a solid model is required to calculate the amount of released pyrolysis gases. Here, the integral model for the pyrolysis of charring materials of Moghtaderi et al., extended with a cooling stage is considered. When an incident heat flux measured during a flame spread experiment, is imposed to the solid, the cooling stage is shown to be indispensable to solve the pyrolysis process. Next, the integral model is compared with a moving grid model. The latter uses the same physical model but does not make any assumption on the temperature profile and thus gives the true solution of the physical model. Comparison of both models reveals that the integral model has problems with suddenly varying boundary conditions and produces erroneous results for parameter studies on thickness and rear boundary condition variations. However, when compared to inert pyrolysis experiments, the integral and the moving grid model are comparable in quality. An automatic optimisation technique is described to obtain material fire properties. r 2004 Elsevier Ltd. All rights reserved. Keywords: Integral model; Pyrolysis; Charring materials
1. Introduction Flame spread is an important aspect of fire. It determines the fire growth and thus the intensity and severity of a fire. Most of the current commercial field (CFD) and Corresponding author. Tel. +32 92 64 33 14; fax. +32 92 64 35 86.
E-mail address: [email protected] (B. Merci). 0379-7112/$ - see front matter r 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.firesaf.2004.09.003
ARTICLE IN PRESS E. Theuns et al. / Fire Safety Journal 40 (2005) 121–140
122
Nomenclature a0 ðtÞ; a1 ðtÞ; a2 ðtÞ coefficients of the quadratic char temperature profile b0 ðtÞ; b1 ðtÞ; b2 ðtÞ coefficients of the quadratic virgin temperature profile c specific heat capacity (J/kg K) h convection coefficient (W/m2 K) k thermal conductivity (W/m K) L total thickness of solid (m) _ 00 m mass flux (kg/m2 s) 00 q_ bs back surface heat flux (W/m2) 00 q_ c heat flux from char layer to pyrolysis front (W/m2) 00 q_ ext external heat flux (W/m2) 00 q_ fl flame heat flux (W/m2) 00 q_ net net incident heat flux (W/m2) 00 q_ v heat flux from pyrolysis front to virgin material (W/m2) DQpyr latent pyrolysis heat (J/kg) T temperature (K) t time (s) x distance (m) Greek symbols a d e y r s
thermal diffusivity (m2/s) layer thickness (m) emissivity of solid surface (dimensionless) integrated temperature (m K) density (kg/m3) Stefan–Boltzmann constant (W/m2 K4)
Subscripts bs c fr g 0 pyr s v N
back surface char pyrolysis front pyrolysis gases initial pyrolysis front surface virgin ambient
ARTICLE IN PRESS E. Theuns et al. / Fire Safety Journal 40 (2005) 121–140
123
zone models are not capable of modelling fire spread [1]. Their use is limited to local fires or fires where the growth is defined before the simulation. The fire growth is thus independent of the conditions during the simulation, like temperature or incident heat flux. This is an important drawback of these models. When the fire growth could be predicted from an initiating fire (e.g. a burning cigarette), the utility of fire models would substantially improve. Therefore, there is a need for fire spread modelling. When solving a fire spread problem we have to deal with the gas phase where combustion takes place, and with the solid phase which determines the reaction of the solid material to the conditions in the gas phase. The reaction of the solid material to the incident heat flux is important for modelling fire spread and is described by a pyrolysis model. Different types of models exist to describe the pyrolysis of materials, going from simple analytical equations to complex numerical models. In the simple analytical models the gasification rate is assumed to be proportional to the net absorbed heat flux, while most of the physical and chemical processes are ignored [2]. The more advanced numerical models can be classified in two categories: models where the pyrolysis front is solved in detail, and models where the pyrolysis front is reduced to a surface at constant pyrolysis temperature. One of the initial numerical models in the first category was developed by Kung [3]. The decomposition of wood is described by single-step Arrhenius equation and the solid is regarded as being partly residual char and partly as yet unpyrolysed active material. Similar models, described by partial differential equations, have been developed [3–9]. An overview of more complex and sophisticated models is given by Di Blasi [10]. In the second category of the numerical models, the pyrolysis front divides the solid into a separate virgin and char layer, for charring material. For noncharring materials the solid exists entirely of virgin material and the pyrolysis front coincides with the front surface. Often integral methods are applied to reduce the calculation time. In the models of Chen et al. [11], Delichatsios et al. [12] and Delichatsios and Saito [13] the first two moments of the heat equations are solved by assuming an exponential temperature profile. In the models Moghtaderi et al. [14,15], Quintiere and Iqbal [16], Spearpoint and Quintiere [17] and Staggs [18] the second moment of the heat equation is not considered, and a quadratic temperature profile is assumed. Other methods in this category do not assume a temperature profile, but use a moving grid [19–21] or a Landau transformation [22]. The model of Galgano and Di Blasi [23] has characteristics of the two categories: degradation reactions are described by an Arrhenius rate equation but they take place at an infinitely thin surface. The temperature of the pyrolysis front is variable. The integral model that is evaluated here, is based on the work of Moghtaderi et al. on the pyrolysis of charring materials [14,15]. The integral model has already been used for predicting the ignition behaviour [17], and the pyrolysis of materials [14]. Here, the model is extended with a cooling stage where virgin and char material are both present but the pyrolysis reactions are interrupted. As will be shown, this stage is necessary for fire simulations. The integral model is also compared to a moving
ARTICLE IN PRESS 124
E. Theuns et al. / Fire Safety Journal 40 (2005) 121–140
grid method for several cases. This allows a true evaluation of the simplifications made by the integral model. Both models are also compared to data from inert pyrolysis experiments.
2. Extended integral model When a virgin solid material is heated by an external heat flux, the temperature in the solid will rise. Once the temperature is sufficiently high, endothermic pyrolysis reactions will start. The virgin material is transformed into char and combustible pyrolysis gases. The latter are very important in fire calculations as these gases mix with ambient air and burn. Part of the energy released in the flame is fed back to the solid and assists maintaining the pyrolysis process. At low temperature the rate pyrolysis reactions are negligible, only at high temperatures they become intense. Most often it is assumed that degradation takes place at a constant pyrolysis temperature, which is considered a material property. The pyrolysis zone in the solid is then reduced to a pyrolysis surface at constant temperature, separating the pure virgin and pure char zone. 2.1. Mathematical description of the pyrolysis stage The integral model used here is based on the work of Moghtaderi [14]. For completeness some formulas will be repeated here, more details can be found in the cited work. The governing equations for the integral model are based on the conservation of energy for the complete char and virgin layer. For the pyrolysis stage the conservation of energy for the char layer is given by Z dc d ddc _ 00g cg ðT s T pyr Þ: T pyr ¼ q_ 00net q_ 00c m rc cc T c dx rc cc (1) dt 0 dt The last term in Eq. (1) represents the heating of the pyrolysis gases when flowing from the point where they are produced, i.e. the pyrolysis front, to the solid surface where they are released in the atmosphere. In some models this term is omitted. Depending on the experiment and the stage of pyrolysis, the net incident heat flux consists of an external incident heat flux, a flame heat flux, reradiation, and convective cooling: qT c kc ¼ q_ 00net ¼ q_ 00ext þ q_ 00fl sðT 4s T 41 Þ hs ðT s T 1 Þ: (2) qx x¼0 For the heating stages and inert pyrolysis experiments, no flame is present and thus q_ 00fl ¼ 0: The description of the virgin zone depends on the location of the thermal penetration front ( ¼ influenced virgin thickness). If the front has not reached the rear surface yet, the virgin layer is treated as semi-infinite. The conservation of
ARTICLE IN PRESS E. Theuns et al. / Fire Safety Journal 40 (2005) 121–140
energy for the semi-infinite pyrolysis and cooling stage is then given by Z dc þdv d ddc ddc ddv T pyr rv cv þ rv cv T v dx þ rv cv T 0 ¼ q_ 00v : dt dc dt dt dt
125
(3)
The influenced virgin thickness (dv ) is always variable, while the char layer thickness (dc ) is only a variable for the pyrolysis stage. Once the thermal penetration front has reached the rear surface the virgin layer is treated as finite. There is no longer uninfluenced virgin material and the conservation equation for the virgin zone becomes Z L d ddc T pyr ¼ q_ 00v þ q_ 00bs : r cv T v dx þ rv cv (4) dt dc v dt In the integral model, the shape of the temperature profile in the char and virgin layer is prescribed. Typical profiles are quadratic [14,15,17] or exponential [13]. Here, the temperature distributions in the char and virgin layer are assumed quadratic and are given by T c ðx; tÞ ¼ a0 ðtÞ þ a1 ðtÞðdc xÞ þ a2 ðtÞðdc xÞ2 ; T v ðx; tÞ ¼ b0 ðtÞ þ b1 ðtÞðx dc Þ þ b2 ðtÞðx dc Þ2 :
ð5Þ
The unknown coefficients a0 ; a1 ; a2 ; b0 ; b1 and b2 are determined with conservation equation (1) and (3) or (1) and (4), and the boundary conditions, e.g. for the pyrolysis stages the temperature at the pyrolysis front is equal to the pyrolysis temperature. The finite and semi-infinite pyrolysis model can be reduced to a system of three differential equations in three unknowns. When yc and yv the integrated temperature for the char and virgin layer, respectively, are defined as Z dc Z dv þdc yc ¼ ðT c T pyr Þ dx; yv ¼ ðT v T pyr Þ dx (6) 0
dv
the set of equations for the semi-infinite pyrolysis model reads ! ddc 1 ¼ ðrv rc ÞDQpyr þ rv cv ðT pyr T 0 Þ dt dyc rv cv dyv 00 þ q_ net rc cc ; dt 2 dt ! dyc yc q_ 00 ¼ 3ac net 2 ; dt 2kc dc dyv 8av ðT 0 T pyr Þ2 ddc ¼ : þ 2ðT pyr T 0 Þ dt 3yv dt
ð7Þ
ARTICLE IN PRESS 126
E. Theuns et al. / Fire Safety Journal 40 (2005) 121–140
After transition to the finite stage, the thermal front no longer exists and the pyrolysis process is described by ! ddc 1
¼ dt rv rc DQpyr dyc dyv q_ 00net rc cc rv c v þ q_ 00bs ; dt dt ! 00 dyc yc q_ ¼ 3ac net 2 ; dt 2kc dc 00 dyv yv q_ bs ¼ 3av : ð8Þ dt 2kv ðL dc Þ2 The last set of system equations, Eqs. (8), is slightly different form the formulation in Moghtaderi et al. [14]. Indeed, instead of assuming a perfectly insulated or isothermal rear boundary, a convective boundary condition is applied. The system equations are more general because, by changing the convection coefficient, an insulated and isothermal rear boundary can be simulated as well. The rate of released pyrolysis gases is determined by the advancement of the pyrolysis front. By coupling the virgin and char solution, the mass loss can be calculated from _ 00g ¼ ðrv rc Þ m
ddc ¼ ðkv ðqT v =qxÞx¼dþc dt kc ðqT c =qxÞx¼dc Þ=DQpyr :
ð9Þ
2.2. Mathematical description of the interrupted pyrolysis stage When the external heat flux is constant (and high enough), all the virgin material will transform into char in a continuous way. But when the external heat flux is variable, it is possible that the pyrolysis process is interrupted. When the external heat flux decreases, less heat is provided at the pyrolysis front and thus less heat is available for the endothermic pyrolysis reactions. When an insufficient amount of heat is supplied, reactions will come to a standstill and a cooling stage is introduced. This is possible for both the semi-infinite and finite pyrolysis stage. When the external heat flux rises again the pyrolysis reactions can resume. Depending on the way of heating (external heat flux), the semi-infinite cooling stage can also result in the finite cooling stage. Note that the name cooling stage is chosen because this stage is most probable to occur when the pyrolysing solid is suddenly cooled down, although, this stage can occur as well for a reduced front surface heat flux. The two cooling stages are required for the simulation of real fires, as will be shown later. They are, however, not included in the formulations of the integral model for the pyrolysis of charring materials found in literature [12,14,17]. More details about the cooling stages can be found in Theuns et al. [24].
ARTICLE IN PRESS E. Theuns et al. / Fire Safety Journal 40 (2005) 121–140
127
3. Evaluation of the integral model For certain initial and boundary conditions an analytical solution exists of the integral model, see for example the heating of a semi-infinite solid with a constant net heat flux [25]. For most variable boundary conditions no analytical solutions exist, and therefore, the results of the integral model are often compared with experiments [12,14,17]. A drawback of the latter, is the lack of knowledge with respect to the modelling of the experiment, and particularly, of the incident heat flux. The differences between experimental data and simulation results can be due to the assumptions typical for the integral model (prescribed temperature profile), due to the representation of the experiment (modelling of incident heat flux, material properties), due to the physical model, due to numerical errors or due to a combination of these types of errors. 3.1. Moving grid model In this work the integral model is compared to a moving grid model [24]. Both models are based on the same physical model: the pyrolysis zone is represented by a single surface at pyrolysis temperature where virgin material is entirely transformed into char. The moving grid model, however, does not rely on any assumption for the temperature distribution and hence the moving grid acts as a numerical reference solution for the physical model. When the same material properties are used for the integral and moving grid model the evaluation of the integral model is rigorous, which is a major advantage. Indeed, when the quality of the integral model is judged on the basis of experiments in a reacting atmosphere, the heat feedback of the flame must be calculated. The simulation then becomes more complex and the modelling uncertainty is larger. Also, the determination of the material properties interferes with the evaluation of the model. In a real fire or fire test, the solid is heated mainly by radiation and convection. The incident heat flux can originate from an external heater, e.g. radiating cone in the Cone Calorimeter test, flames at the solid surface, a hot gas layer, etc. In this paper the model is run as a standalone model and therefore, the net incident heat flux is supplied as a boundary condition. The external heat flux is defined before the simulation and is not dependent on the mass release rate. When the solid model is coupled to a CFD model or zone model this net incident heat flux is computed by the code. All simulations for the rigorous numerical comparisons, are performed on particle board, for which the material properties are given in Table 1 (initial). The following cases are examined:
experimentally measured heat flux; sudden variation of external heat flux; variation of solid thickness; variation of rear boundary condition.
Finally, the integral and moving grid model are compared to inert pyrolysis experiments.
ARTICLE IN PRESS E. Theuns et al. / Fire Safety Journal 40 (2005) 121–140
128
Table 1 Initial guess and optimised material properties rc kv kc cv cc T pyr DQpyr rv (kg/m3) (kg/m3) (W/m K) (W/m K) (J/kg K) (J/kg K) (1C) (J/kg) Initial [29] 600 Optimised 670
60 117
0.36 0.15
0.23 0.38
2500 4655
375 280
8.7 105 0 5.1 105 171
1.0 1.0
60
0.015
measured heat flux (kW/m2)
moving grid integral arrhenius
m" (kg/m2.s)
2500 4463
v (—) cg (J/kg K)
0.010
0.005
40
20
0
0.000 0
300
600
time (s)
900
0
300
600
900
time (s)
Fig. 1. (Left) Mass flux of pyrolysis gases and (right) experimentally measured heat flux.
3.2. Experimentally measured heat flux In the first case the importance of the cooling stage is illustrated. An experimentally measured heat flux is imposed to the solid. The incident heat flux was measured in a vertical flame spread experiment on particle board [26], see Fig. 1 (right). The heat flux peaks after about 180 s, and afterwards drops drastically. This kind of heating has of course a strong influence on the pyrolysis behaviour. As a consequence of the drop in the heat flux, the cooling stage occurs. This is indeed observed in Fig. 1 (left): the release of pyrolysis gases drops to zero between 322 and 460 s. If the cooling stage were not included in the integral model, the experiment could only be adequately simulated during the first 322 s. The prediction of the mass release rate of pyrolysis gases is very good when compared to the moving grid model. Only immediately after the first peak (at 175 s), the mass release rate is underestimated with errors up to 24%. The time interval where the pyrolysis reactions are interrupted and the time of resuming of the reactions is predicted very well, with deviations less than 15 s. This indicates that the description of the cooling stage by a quadratic temperature profile in char and virgin layer is valid for a typical flame spread heat flux. In reality it is not probable that the mass release suddenly drops to zero. This can be illustrated by using for example the time-consuming Arrhenius model to predict the mass release. In this model the pyrolysis reactions are described by an Arrhenius
ARTICLE IN PRESS E. Theuns et al. / Fire Safety Journal 40 (2005) 121–140
129
law [4–6,27,28]: qr ¼ ðr rc ÞA expðE a =RTÞ; qt
(10)
where r is the local instantaneous solid density, and rc the final char density, A the pre-exponential factor and E a the activation energy. In the pyrolysis zone the solid consists of a local mixture of virgin and char material. Hence, the density (and material) changes in a continuous way, and not abrupt as in the integral and moving grid model. In this model the pyrolysis zone is finite and pyrolysis reactions proceed at lower temperatures as well, though at a lower rate. Therefore, the mass release rate in the Arrhenius model is still significant during the cooling stage of the integral model, as can be seen in Fig. 1. The temperature in the pyrolysis zone is still high enough to produce a significant mass release of pyrolysis gases. The weakness of the present integral model to predict the mass loss rate when the external heat flux becomes too low is associated with the assumption of a constant temperature on the pyrolysis front. When insufficient heat is provided, the front can no longer hold the pyrolysis temperature and comes to a standstill. An alternative can be found in the work of Galgano and Di Blasi [23]. They coupled the integral method with an Arrhenius reaction rate on the pyrolysis front for wood cylinders. The temperature of the front is no longer constant but is calculated during the simulation. The advancement of the front is given by the Arrhenius reaction rate. Such an approach is thought to improve the results for low external heat fluxes.
3.3. Sudden variation of the external heat flux In enclosed fires sometimes sudden variations of the external heat flux of the solid material can be noticed. For example, when flashover occurs, all exposed combustible items in the enclosure suddenly get involved in the combustion process. This phenomenon is here modelled as a sudden rise in the external heat flux some fixed time after the start the simulation (at 300 s). The heat flux suddenly drops after 400 s, which represents the effect of extinguishment of flames or the activation of sprinklers. In reality, the time of flashover, the extinguishment, and the activation of sprinklers will depend on the enclosure conditions. The imposed external heat flux is shown in Fig. 2 (right). The integral model has remarkable results, as shown in Fig. 2 (left). Immediately after the rise of the external heat flux, at 300 s, the mass flux of pyrolysis gases drops sharply instead of rising. The reason is the unrealistic, direct influence of the net incident heat flux on the complete temperature profile of the solid. As can be seen in Fig. 2, already after 1 s the integral model has recovered and the mass release is normal again. When the change in the external heat flux occurs at the start of pyrolysis instead at 300 s, the mass release behaves normal (not shown here). When the external heat flux suddenly drops, the integral model shows again a bad prediction just after the fall of the external heat flux, see Fig. 2 at 400 s. Again this is due to the unrealistic, direct influence of the net incident heat flux on the complete
ARTICLE IN PRESS E. Theuns et al. / Fire Safety Journal 40 (2005) 121–140
130 0.020
60 moving grid integral
q"ext (kW/m2)
m" (kg/m2.s)
0.015
0.010
40
20
0.005
0.000
0 100
200
300
time (s)
400
100
200
300
400
time (s)
Fig. 2. (Left) Mass flux of pyrolysis gases and (right) external heat flux.
temperature profile. There is a sudden short rise in the mass flux of pyrolysis gases, but less than a few seconds later, the mass flux is normal again. A concluding remark with respect to the occurrence of the erroneous peaks with the integral model is that as the shape of the temperature profile is prescribed (e.g. quadratic), the integral model can only be valid for types of heating which will result in that prescribed temperature profile. In flame spread simulations, the incident heat flux can suddenly rise due to for example ignition of the pyrolysis gases in the gas stage, but it can also suddenly fall due to for example lack of oxygen in the gas stage. This kind of heating can result in temperature profiles that are different from the quadratic ones. A simple example is given by solving the heating of a semi-infinite slab (no pyrolysis) with a quadratic temperature profile. At the thermal penetration depth, the temperature is equal to the initial slab temperature and the slope of the profile is zero because no conductive heat flows to the unheated zone. When the net incident heat flux suddenly changes from positive to negative, no realistic solution can be found for this problem. The temperature profile is indeed entirely determined by the net heat flux (kv ðqT=qxÞx¼0 ¼ q_ 00net ), the zero heat flux at the end of the thermal layer (kv ðqT=qxÞx¼dv ¼ 0), and the temperature at the thermal penetration depth. (T v jx¼dn ¼ T 0 ). These boundary conditions do not allow conservation of energy for a quadratic temperature profile and thus no acceptable solution exists for such event. 3.4. Variation of the solid thickness Next, the integral model is here used for a parameter study of the solid thickness. The thickness is varied from 1 mm, in fires typically thermally thin, to 50 mm, typically thermally thick. The results of the integral model are shown in Fig. 3, the results of the moving grid model in Fig. 4. For all these cases, the rear surface is perfectly insulated and the external heat flux is constant (q_ 00ext ¼ 50 kW=m2 ). For a small thickness the mass release is concentrated in one short peak and the solid is consumed very quickly. For a thick solid the mass release rate contains two
ARTICLE IN PRESS E. Theuns et al. / Fire Safety Journal 40 (2005) 121–140
131
0.04
L = 1 mm
m" (kg /m2.s)
0.03
Integral model
2 mm
0.02
5 mm 10 mm
0.01
20 mm 50 mm 0.00 10
100
1000
10000
time (s)
Fig. 3. Mass flux of pyrolysis gases for several thicknesses (integral model).
0.04
Moving grid
L = 1 mm
m" (kg/m2.s)
0.03
2 mm
0.02
5 mm 10 mm
0.01
20 mm 50 mm 0.00 10
100
1000
10000
time (s)
Fig. 4. Mass flux of pyrolysis gases for several thicknesses (moving grid model).
peaks. This first peak occurs at the start of pyrolysis, the second peak occurs near the end, when the pyrolysis front is reaching the rear surface of the solid. It is due to the back-effect, caused by the insulated rear surface [17,29]. When the pyrolysis front reaches the rear surface, continuously less energy flows from the pyrolysis front to the virgin material. The energy flowing to the virgin material, cannot leave the solid because of the insulated rear surface. Subsequently, more heat is available for the heat absorbing pyrolysis reactions at the pyrolysis front and hence the mass loss rate increases. When, for example, there is a strong convective boundary condition at the rear surface, the second peak does not occur. Also, for real thick materials this second peak is less pronounced. The growing char layer acts as an thermal insulator and causes the decrease of the mass release rate.
ARTICLE IN PRESS 132
E. Theuns et al. / Fire Safety Journal 40 (2005) 121–140
The first peaks in the mass release rates, see Figs. 3 and 4, coincide for thicknesses of 10, 20 and 50 mm. This is the region where the solid is calculated with the semiinfinite equations and no influence of the rear boundary or thickness is felt. Once the pyrolysis process is in the finite region, a thicker material produces less pyrolysis gases for a certain location of the pyrolysis front (for a insulated rear surface and constant incident heat flux). Indeed, for the thick material a significant amount of heat is ‘‘lost’’ at the pyrolysis front to the remaining cold virgin material which leaves less energy for the endothermic pyrolysis reactions. For a thinner material, on the other hand, the remaining virgin material has a higher temperature because there is less virgin material to be heated (no heat loss at rear surface). Consequently, the heat flow from the pyrolysis front to the virgin layer is lower, which leaves more energy for the pyrolysis reactions. As in this case the external heat flux is constant, the mass release rate is indeed lower, as can be seen in Fig. 4. The results for the integral model, in Fig. 3, again reveal a shortcoming. As long as the solid can be treated as semi-finite the mass release curves coincide (e.g. 10 and 20 mm overlap from 26 to 50 s). At 51 s the thinner solid (10 mm) is solved by finite theory, while the thicker (20 mm) is still solved by semi-infinite theory. The mass release rate for the thinner material should be higher, as explained above, but the results show the opposite. Due to the transition to the finite stage for the thinner material, the correct behaviour is no longer produced. Before transition, the thermal front is still advancing through the virgin layer and the pyrolysis process is described by Eqs. (7). After transition the thermal front no longer exists and the pyrolysis process is described by Eqs. (8). At 115 s the correct behaviour, i.e. higher mass release rate for the thinner material, is given again. The problem at transitions in the integral model manifests itself in a more severe error. Looking at the mass release rate of pyrolysis gases for a thickness of 5 and 10 mm in Fig. 3, we see that the integral model predicts the first mass release (often used as the condition for ignition) for the thicker material (10 mm) earlier than for the thin material (5 mm). It is clear that for a solid with insulated rear surface this cannot be the case. The correct behaviour is shown by the moving grid model in Fig. 4. The reason for the flaw in the integral model is the underprediction of the surface temperature after the transition from semi-infinite to finite stage. Immediately after the transition from the semi-infinite to the finite stage, the surface temperature rises slower than for a thicker solid still in the semi-infinite stage. As a consequence, it takes longer for the surface to reach the pyrolysis temperature. For the results with the integral model in Fig. 3, the thin solid (5 mm) is already treated as finite when the first mass release occurs, while the thicker solid (10 mm) is still treated as semi-infinite at the start of mass release. So, again it is the transition from semi-infinite to finite stage that is responsible for the error. 3.5. Variation of rear boundary condition In order to examine the back-effect, the rear boundary condition is varied. In all cases a convective heat loss is applied: q_ 00bs ¼ hbs ðT bs T 1 Þ:
(11)
ARTICLE IN PRESS E. Theuns et al. / Fire Safety Journal 40 (2005) 121–140
133
The convection coefficient hbs is varied from 0 to 20 W/m2 K. In these simulations the thickness of the solid was 20 mm and the external radiative heat flux 50 kW/m2. Results are shown for the integral model in Fig. 5, and for the moving grid model in Fig. 6. When the convection coefficient is raised, the second peak in the mass release rate becomes less pronounced and eventually disappears. In the integral model the pyrolysis process is described by the semi-infinite model for times shorter than 183 s. All mass release curves coincide because the rear boundary condition does not yet appear in the system equations for the integral model. Of course, it is assumed here that the ambient temperature at the rear surface is equal to the initial temperature. If the ambient temperature is variable or different from the initial temperature, the integral model in its present form, is no longer applicable for the semi-infinite stages. At the transition from semi-infinite to finite pyrolysis, it is expected that the higher the convection coefficient, the higher the heat loss, and thus the lower the mass release rate, as is shown by the moving grid in Fig. 6. For the integral model the opposite is true. Higher heat losses give temporarily (between 183 and 400 s) higher mass release rates, see Fig. 5, which is incorrect. The transition from semi-infinite to finite pyrolysis stage is responsible for this error. 3.6. Comparison to experiments The integral and moving grid model are compared with inert pyrolysis experiments on particle board [30]. Due to the inert atmosphere, there is no combustion of pyrolysis gases in the atmosphere. Modelling of flames and flame heat fluxes is avoided, and thus the incident heat flux is well known. Experimental results for particle board are available for two variable external radiant heat flux rates, where the rise was given by 16.6 kW/m2 s (exp16) and 33.3 kW/m2 s (exp33), and one constant radiant heat flux of 31 kW/m2 (exp31). The exposed surface of the samples loses energy by convection and reradiation. During the initial heating stage, thus at low temperatures, the convective cooling is 0.012
m" (kg/m2.s)
Integral model
Detail
0.008
h = 0 W/m2K 0.004
5 W/m2K 20 W/m2K 0.000 10
100
1000
10 W/m2K 10000
time (s)
Fig. 5. Mass flux of pyrolysis gases for several rear boundary conditions (integral model).
ARTICLE IN PRESS E. Theuns et al. / Fire Safety Journal 40 (2005) 121–140
134
0.012
m" (kg/m2.s)
Moving grid
Detail
0.008
h = 0 W/m2K 0.004
5 W/m2K
20 W/m2K
10 W/m2K
0.000 10
100
1000
10000
time (s)
Fig. 6. Mass flux of pyrolysis gases for several rear boundary conditions (moving grid model).
important. Therefore, in the simulations the convection coefficient hs for the natural convective cooling is determined by the correlation [31]: hs ¼ 0:54
l ðGr PrÞ0:25 L
(12)
instead of using a constant value. Hence, the convection coefficient is dependent on the surface temperature. A value of about 6 W/m2 K is typical for a surface temperature of 30 1C, while 16 W/m2 K for 500 1C. Evidently, at high surface temperatures the reradiation becomes more important. A major problem of pyrolysis models, is the determination of the (equivalent) material properties. Not all the material properties are known or can be found in literature. Sometimes some properties are estimated [28], but often the material properties are varied until good agreement with experiments is obtained [29]. The latter is done here as well. 3.6.1. Determination of material properties by error minimisation Optimisation of material properties is done only with the moving grid model. By adjusting the material properties, the simulated mass release rates of exp16 and exp31 are optimised simultaneously, leaving exp33 as a blind calculation. At each time interval (Dt) the simulated mass flux is compared with the experimental value. The errors for exp16 and exp33 are multiplied and the result defines the function that has been minimised: errorðrv ; rc ; kv ; kc ; cv ; cc ; T pyr ; Qpyr ; cg ; v Þ X M M X 00 00 _ 00g; exp 16 ði DtÞ _ 00g; exp 31 ði DtÞ; _ g; sim16 ði DtÞ m _ g; sim31 ði DtÞ m ¼ m m i¼1
i¼1
ð13Þ
ARTICLE IN PRESS E. Theuns et al. / Fire Safety Journal 40 (2005) 121–140
135
where M is large enough to capture the total mass release curve. The first term in the RHS represents the error of the simulation of exp16, the second the error of the simulation of exp31. The result of the two terms in the RHS are a measure for the difference between the calculated and the experimental result and have the dimensions of a mass flux (kg/m2 s). If these errors are multiplied by the time interval (Dt), a measure of the total mass release error (in kg/m2) is obtained. These errors can then be seen as the difference of the area spanned by the experimental and simulated mass flux curves. The minimisation is done by a downhill simplex method in multidimensions [32]. The method is started by N+1 linear independent points which define a simplex in an N-dimensional space. Here N is the number of material properties that are optimised. By reflecting, expanding or contracting points of the simplex, a minimum of the function is found. The iterative optimisation procedure is stopped when the difference between the lowest error and the highest error of the simplex falls below a small number (104). The vertices of the initial simplex were derived from the guessed material properties in Table 1. One vertex was constructed by contracting the initial guess to the origin by 25% in all directions of the N-space, while the other N vertices were created by expanding the initial guess by 25% for only one direction (or material property). This gives the N+1 linear independent points of the starting simplex. During the optimisation procedure the variables are bounded to allow only realistic values. In the calculations presented in this paper the bounds were chosen very wide. As a consequence, during the optimisation there were no restrictions for the material properties. The high number of variable parameters can be reduced by determining some parameters by independent tests. For example, the virgin density can easily be determined by weight and volume measurements. In our method the virgin density was deliberately left free, which allows to check whether the method predicts realistic values. Instead of using independent tests it is also possible to include extra information in the optimisation procedure, for example, surface or interior temperature measurements. The error in Eq. (13) should then include in addition of the mass loss error also a temperature error. 3.6.2. Discussion of the optimised properties Minimisation has been done for 10 variable material properties. The 10 variables are: the virgin and char density, the virgin and char conductivity, the virgin and char specific heat capacity, the pyrolysis temperature, the pyrolysis heat, the volatiles specific heat capacity and the virgin surface emissivity. The optimised results for the moving grid model are given in Table 1. The virgin density is realistic and resembles the value of de Ris and Yan [29]. The other properties strongly differ with the initial ones, e.g. the char properties and the pyrolysis temperature. Differences can be due to different types of particle board used in the work of de Ris and Yan [29] and in the experiments [30]. The high-optimised virgin heat capacity can be explained by the vaporisation of the water in the particle board. This phenomenon is not directly included in the pyrolysis model and will act as an energy sink in the virgin layer. By increasing the virgin heat capacity the effect of the vaporisation is partially included
ARTICLE IN PRESS E. Theuns et al. / Fire Safety Journal 40 (2005) 121–140
136
in the pyrolysis model. The optimised char heat capacity seems rather high as well, but the influence of the char heat capacity in the simulations is rather small. This has been tested by setting the char heat capacity at 2500 J/kg K. The effect of the char heat capacity is in such simulations small because it appears in the heat equation in combination with the char density (rc) which is small. Similar observations have been made by Chen et al. [33]. They noticed that the thermal capacity ratio R ¼ ðrc cc Þ=ðrv cv Þ does not significantly affect the mass loss and surface temperature histories. By bounding some variable parameters during the optimisation unrealistic results can be avoided. Here, the char heat capacity was kept constant (2500 J/kg K), but tight bounds could be used as well. The optimised material properties are only determined by information of the mass release experiments exp16 and exp31. They do well in predicting the three experiments, see Fig. 7. The experimental curves in Fig. 7 are the mean of the two experimental results in [30]. The good results for exp16 and exp31 could be expected as the material properties were optimised for those two experiments. Exp33 on the other hand was not included, and thus represents a blind prediction. In Fig. 7 it is noticed that the mass loss predictions for exp16 and exp33 are less satisfying for the start and end of pyrolysis. In the beginning of the experiment the incident heat flux is still low, for exp16 the heat flux is rising linearly by 1.66 103 W/cm2 s and for exp33 by 3.33 103 W/cm2 s. As a consequence, degradation reactions will proceed at relative low temperature. Models with a constant surface pyrolysis temperature cannot predict this early mass loss because the solid must first attain the pyrolysis temperature before mass can be released. In the work of Chen [33] it is stated that a surface pyrolysis front at constant temperature is valid for the range 10–100 kW/m2. In a real fire the incident heat flux to most materials will rise slowly and sometimes will stay low. For those low heat fluxes the integral and moving grid model may be less suitable and an Arrhenius rate equation for the pyrolysis reactions may be more appropriate. Still fire models that use a surface pyrolysis front at constant temperature in the solid phase, give good 0.009
Q=31kW/m2
m" (kg/m2.s)
Q=33kW/m2
exp inte mov
0.006
Q=16kW/m2 0.003
0.000 0
1000
2000
3000
time (s)
Fig. 7. Mass flux of pyrolysis gases for particle board exposed to several radiation intensities.
ARTICLE IN PRESS E. Theuns et al. / Fire Safety Journal 40 (2005) 121–140
137
results for fire and flame spread calculations [12,28,34]. By including experiments with low heat flux in the optimisation methodology, the material fire properties for the integral and moving grid model are thought to perform better for a wider range of incident heat fluxes (low and high). 3.6.3. Difference between integral and moving grid model The optimised material properties are used for predicting exp16, exp31 and exp33 with the integral and moving grid model. The differences between the integral and moving grid predictions are small for the three experiments, as can be seen in Fig. 7. The difference is only large at the beginning and end of pyrolysis. The reason is that the mass release rises or falls quickly at these points and the ignition and extinguishment times are not predicted to be equal. During pyrolysis it is noticed that at first, the integral model overpredicts the mass release rate, then gives the same value as the moving grid model, and for the remaining part underpredicts the mass release rate. When errors at the beginning and end of pyrolysis are not considered, the maximum difference between integral and moving grid model is about 1.5 103 kg/m2 s or about 15% of the peak mass release rate (exp31). The average difference is for the three experiments less than 1 104 kg/m2 s or about 1% of the peak mass release rate. The error ( ¼ difference between experiment and model) of the integral model is a bit higher than of the moving grid model, see Table 2. When we compare that error with the total experimental mass release, the difference between the two models can be neglected. The error made by approximating the temperature profile by a quadratic function is much smaller than the error made by the physical model. Therefore, we can say that the integral model is a suitable method for calculating the release of pyrolysis gases as described by the physical model.
4. Conclusions An integral model extended with a cooling stage has been evaluated. The cooling stages occur when an insufficient amount of energy is supplied to the solid and hence the pyrolysis reactions are interrupted. When the heat supply rises again, the pyrolysis reactions can resume. The occurrence of the cooling stage was illustrated for an experimental heat flux. The integral model predicts the interrupted pyrolysis well, although models with an Arrhenius rate equation, still predict a low mass release rate during the cooling stage of the integral model. Table 2 Error of mass release rate prediction
Integral (kg/m2) Moving (kg/m2) Total experimental mass release (kg/m2)
Exp16
Exp31
Exp33
1.631 1.593 5.280
0.528 0.496 4.530
1.490 1.391 5.536
ARTICLE IN PRESS 138
E. Theuns et al. / Fire Safety Journal 40 (2005) 121–140
A comparison of the integral model with the moving grid model is performed to test the validity of the prescribed quadratic temperature profile shape. It was illustrated that, due to this assumption the integral model is not always capable of giving correct or realistic results. Sudden variations in the boundary conditions are transmitted immediately and unrealistically high through the entire solid and lead to temporarily unrealistic results. On the other hand, it must be noted that the integral model recovers quickly from this error. A practically measured heat flux or slowly changing boundary conditions do not cause any problems. Next, it was illustrated that the integral model is not suited for parameter studies. At the transition from the semi-infinite to the finite stage, the correct behaviour is not always predicted. The most important example is the wrong prediction of the first release of pyrolysis gases when varying the solid thickness. Often the first mass release is taken as the condition for ignition time or ignition temperature and is as such very important. Also, when the rear boundary condition is varied, unrealistic results are noted. If the integral model is used as a pyrolysis model, the user should be aware of its limits and shortcomings. On the other hand, in the study of inert pyrolysis experiments the differences between the mass release rate of the integral and moving grid model are very small. Only temporarily, differences of 15% of the peak mass release rate are possible. The error of the physical model is much larger than the error introduced by solving the physical model with the integral model. This implies that the extended integral model is useful as pyrolysis model. Finally, it is noted that an automatic optimisation procedure was presented to determine the material properties for the solid pyrolysis experiments. It is based on fitting the simulated mass release curves to combined experimental results. This has been done by an automatic minimisation of the mass release error.
Acknowledgement The second author is a Postdoctoral Fellow of the Fund of Scientific Research— Flanders (Belgium) (FWO-Vlaanderen).
References [1] Olenick SM, Carpenter DJ. An updated international survey of computer models for fire and smoke. J Fire Protect Eng 2003;13:87–110. [2] Tewarson A. SFPE handbook of fire engineering, 2nd ed. Boston, MA: Society of Fire Protection Engineers; 1995. p. 3.53–3.124. [3] Kung HC. A mathematical model of wood pyrolysis. Combust Flame 1972;18:185–95. [4] Di Blasi C. Processes of flames spreading over the surface of charring solid fuels: effects of fuel thickness. Combust Flame 1994;97:225–39. [5] Staggs JEJ. A simple model of polymer pyrolysis including transport of volatiles. Fire Saf J 2000;34:69–80.
ARTICLE IN PRESS E. Theuns et al. / Fire Safety Journal 40 (2005) 121–140
139
[6] Yuen R, Casey R, De Vahl, Davis G, Leonardi E, Yeoh GH, Chandrasekaran V, Grubits SJ. A threedimensional mathematical model for the pyrolysis of wet wood. In: Hasemi Y, editor. Fire safety science—proceedings of the fifth international symposium, 1997. p.189–200. [7] Wichman IS, Atreya A. A simplified model for the pyrolysis of charring materials. Combust Flame 1987;68:231–47. [8] Staggs JEJ, Whiteley RH. Modelling the combustion of solid-phase fuels in cone calorimeter experiments. Fire Mater 1999;23:63–9. [9] Novozhilov V, Moghtaderi B, Fletcher DF, Kent JH. Computational fluid dynamics modelling of wood combustion. Fire Saf J 1996;27:69–84. [10] Di Blasi C. Modeling and simulation of combustion processes of charring and non-charring solid fuels. Prog. Energy Combust 1993;19:71–104. [11] Chen Y, Delichatsios MA, Motevalli V. Material pyrolysis properties, Part I: An integral model for one-dimensional transient pyrolysis of charring and non-charring materials. Combust Sci Technol 1993;88:309–28. [12] Delichatsios MM, Mathews MK, Delichatsios MA. An upward fire spread and growth simulation. In: Cox G, Langford B, editors. Fire safety science—proceedings of the third international symposium. New York: Elsevier Applied Science; 1991. p. 207–16. [13] Delichatsios MA, Saito K. Upward fire spread: key flammability properties, similarity solutions and flammability indices. In: Cox G and Langford B, editors. Fire safety science—proceedings of the third international symposium. New York: Elsevier Applied Science; 1991. p. 217–26. [14] Moghtaderi B, Novozhilov V, Fletcher D, Kent JH. An integral model for the transient pyrolysis of solid materials. Fire Mater 1997;21:7–16. [15] Moghtaderi B, Novozhilov V, Fletcher D, Kent JH. Mathematical modelling of the piloted igniton of wet wood using the heat-balance integral method. J Appl Fire Sci 1997;6(2):91–107. [16] Quintiere J, Iqbal N. An approximate integral model for the burning rate of a thermoplastic-like material. Fire Mater 1994;18:89–98. [17] Spearpoint MJ, Quintiere JG. Predicting the burning of wood using an integral model. Combust Flame 2000;123:308–24. [18] Staggs JEJ. Short communication: approximate solutions of the pyrolysis of char forming and filled polymers under thermally thick conditions. Fire Mater 2000;24:305–8. [19] Jia F, Galea ER, Patel MK. Numerical simulation of the mass loss process in pyrolizing char materials. Fire Mater 1999;23:71–8. [20] Staggs JEJ, Whiteley RH. Modelling the combustion of solid-phase fuels in cone calorimeter experiments. Fire Mater 1999;23:63–9. [21] Theuns E, Vierendeels J, Vandevelde P. A moving grid model for the pyrolysis of charring. Int J Numer Methods Heat Fluid Flow 2002;12(5):541–59. [22] Staggs JEJ. A discussion of modelling idealised ablative materials with particular reference to fire testing. Fire Saf J 1997;28:47–66. [23] Galgano A, Di Blasi C. Modeling wood degradation by the unreacted-core-shrinking approximation. Ind Eng Chem Res 2003;42:2101–11. [24] Theuns E, Merci B, Vierendeels J, Vandevelde P. Extension and evaluation of the integral model for transient pyrolysis of charring materials. Fire Mater, in press. [25] Carslaw HS, Jaeger JC. Conduction of heat in solids. London: Oxford University Press; 1959. p. 510. [26] Kokkala M, Mikkola E, Immonen M, Juutilainen H, Manner P, Parker WJ. Large-scale upward flame spread tests on wood products year. VTT Research Notes 1834, Technical Research Centre of Finland, Espoo, Finland, 1997. [27] Ritchie SJ, Steckler KD, Hamins A, Cleary TG, Yang JC, Kashiwagi T. The effect of sample size on the heat release of charring materials. In: Hasemi Y, editors. Fire safety science—proceedings of the fifth international symposium, 1997. p.177–88. [28] Yan Z, Holmstedt G. CFD and experimental studies of room fire growth on wall lining materials. Fire Saf J 1996;27:201–38. [29] de Ris JL, Yan Z. Modeling ignition and pyrolysis of charring fuels. Proceedings of the fifth international conference on fire and materials ‘98, San Antonio, TX, USA, 1998. p.111–21.
ARTICLE IN PRESS 140
E. Theuns et al. / Fire Safety Journal 40 (2005) 121–140
[30] Vovelle C, Akrich R, Delfau JL. Mass loss rate measurements on solid materials under radiative heating. Comb Sci Technol 1984;36:1–18. [31] Incropera FP, DeWitt DP. Fundamentals of heat and mass transfer. New York: Wiley; 1996. p. 886. [32] Press WH, Teukolsky SA, Vetterling WT, Flannery BP. Numerical recipes in Fortran 77: the art of scientific computation, vol. 1. Cambridge: Cambridge University Press; 1992. p. 992. [33] Chen Y, Motevalli V, Delichatsios MA. Material pyrolysis properties, Part II: Methodology for the derivation of pyrolysis properties for charring materials. Combust Sci Technol 1995;104(4–6):401–25. [34] Jia F, Galea ER, Patel MK. The numerical simulation of fire spread within a compartment using an integrated gas and solid phase combustion model. J Appl Fire Sci 1999;8(4):327–52.
Fire Safety Journal 40 (2005) 121–140 www.elsevier.com/locate/firesaf
Critical evaluation of an integral model for the pyrolysis of charring materials E. Theuns, B. Merci, J. Vierendeels, P. Vandevelde Department of Flow, Heat and Combustion Mechanics, Ghent University, Gent, Belgium Received 2 February 2004; received in revised form 2 August 2004; accepted 27 September 2004 Available online 2 December 2004
Abstract When field models are used to predict fires and flame spread, a solid model is required to calculate the amount of released pyrolysis gases. Here, the integral model for the pyrolysis of charring materials of Moghtaderi et al., extended with a cooling stage is considered. When an incident heat flux measured during a flame spread experiment, is imposed to the solid, the cooling stage is shown to be indispensable to solve the pyrolysis process. Next, the integral model is compared with a moving grid model. The latter uses the same physical model but does not make any assumption on the temperature profile and thus gives the true solution of the physical model. Comparison of both models reveals that the integral model has problems with suddenly varying boundary conditions and produces erroneous results for parameter studies on thickness and rear boundary condition variations. However, when compared to inert pyrolysis experiments, the integral and the moving grid model are comparable in quality. An automatic optimisation technique is described to obtain material fire properties. r 2004 Elsevier Ltd. All rights reserved. Keywords: Integral model; Pyrolysis; Charring materials
1. Introduction Flame spread is an important aspect of fire. It determines the fire growth and thus the intensity and severity of a fire. Most of the current commercial field (CFD) and Corresponding author. Tel. +32 92 64 33 14; fax. +32 92 64 35 86.
E-mail address: [email protected] (B. Merci). 0379-7112/$ - see front matter r 2004 Elsevier Ltd. All rights reserved. doi:10.1016/j.firesaf.2004.09.003
ARTICLE IN PRESS E. Theuns et al. / Fire Safety Journal 40 (2005) 121–140
122
Nomenclature a0 ðtÞ; a1 ðtÞ; a2 ðtÞ coefficients of the quadratic char temperature profile b0 ðtÞ; b1 ðtÞ; b2 ðtÞ coefficients of the quadratic virgin temperature profile c specific heat capacity (J/kg K) h convection coefficient (W/m2 K) k thermal conductivity (W/m K) L total thickness of solid (m) _ 00 m mass flux (kg/m2 s) 00 q_ bs back surface heat flux (W/m2) 00 q_ c heat flux from char layer to pyrolysis front (W/m2) 00 q_ ext external heat flux (W/m2) 00 q_ fl flame heat flux (W/m2) 00 q_ net net incident heat flux (W/m2) 00 q_ v heat flux from pyrolysis front to virgin material (W/m2) DQpyr latent pyrolysis heat (J/kg) T temperature (K) t time (s) x distance (m) Greek symbols a d e y r s
thermal diffusivity (m2/s) layer thickness (m) emissivity of solid surface (dimensionless) integrated temperature (m K) density (kg/m3) Stefan–Boltzmann constant (W/m2 K4)
Subscripts bs c fr g 0 pyr s v N
back surface char pyrolysis front pyrolysis gases initial pyrolysis front surface virgin ambient
ARTICLE IN PRESS E. Theuns et al. / Fire Safety Journal 40 (2005) 121–140
123
zone models are not capable of modelling fire spread [1]. Their use is limited to local fires or fires where the growth is defined before the simulation. The fire growth is thus independent of the conditions during the simulation, like temperature or incident heat flux. This is an important drawback of these models. When the fire growth could be predicted from an initiating fire (e.g. a burning cigarette), the utility of fire models would substantially improve. Therefore, there is a need for fire spread modelling. When solving a fire spread problem we have to deal with the gas phase where combustion takes place, and with the solid phase which determines the reaction of the solid material to the conditions in the gas phase. The reaction of the solid material to the incident heat flux is important for modelling fire spread and is described by a pyrolysis model. Different types of models exist to describe the pyrolysis of materials, going from simple analytical equations to complex numerical models. In the simple analytical models the gasification rate is assumed to be proportional to the net absorbed heat flux, while most of the physical and chemical processes are ignored [2]. The more advanced numerical models can be classified in two categories: models where the pyrolysis front is solved in detail, and models where the pyrolysis front is reduced to a surface at constant pyrolysis temperature. One of the initial numerical models in the first category was developed by Kung [3]. The decomposition of wood is described by single-step Arrhenius equation and the solid is regarded as being partly residual char and partly as yet unpyrolysed active material. Similar models, described by partial differential equations, have been developed [3–9]. An overview of more complex and sophisticated models is given by Di Blasi [10]. In the second category of the numerical models, the pyrolysis front divides the solid into a separate virgin and char layer, for charring material. For noncharring materials the solid exists entirely of virgin material and the pyrolysis front coincides with the front surface. Often integral methods are applied to reduce the calculation time. In the models of Chen et al. [11], Delichatsios et al. [12] and Delichatsios and Saito [13] the first two moments of the heat equations are solved by assuming an exponential temperature profile. In the models Moghtaderi et al. [14,15], Quintiere and Iqbal [16], Spearpoint and Quintiere [17] and Staggs [18] the second moment of the heat equation is not considered, and a quadratic temperature profile is assumed. Other methods in this category do not assume a temperature profile, but use a moving grid [19–21] or a Landau transformation [22]. The model of Galgano and Di Blasi [23] has characteristics of the two categories: degradation reactions are described by an Arrhenius rate equation but they take place at an infinitely thin surface. The temperature of the pyrolysis front is variable. The integral model that is evaluated here, is based on the work of Moghtaderi et al. on the pyrolysis of charring materials [14,15]. The integral model has already been used for predicting the ignition behaviour [17], and the pyrolysis of materials [14]. Here, the model is extended with a cooling stage where virgin and char material are both present but the pyrolysis reactions are interrupted. As will be shown, this stage is necessary for fire simulations. The integral model is also compared to a moving
ARTICLE IN PRESS 124
E. Theuns et al. / Fire Safety Journal 40 (2005) 121–140
grid method for several cases. This allows a true evaluation of the simplifications made by the integral model. Both models are also compared to data from inert pyrolysis experiments.
2. Extended integral model When a virgin solid material is heated by an external heat flux, the temperature in the solid will rise. Once the temperature is sufficiently high, endothermic pyrolysis reactions will start. The virgin material is transformed into char and combustible pyrolysis gases. The latter are very important in fire calculations as these gases mix with ambient air and burn. Part of the energy released in the flame is fed back to the solid and assists maintaining the pyrolysis process. At low temperature the rate pyrolysis reactions are negligible, only at high temperatures they become intense. Most often it is assumed that degradation takes place at a constant pyrolysis temperature, which is considered a material property. The pyrolysis zone in the solid is then reduced to a pyrolysis surface at constant temperature, separating the pure virgin and pure char zone. 2.1. Mathematical description of the pyrolysis stage The integral model used here is based on the work of Moghtaderi [14]. For completeness some formulas will be repeated here, more details can be found in the cited work. The governing equations for the integral model are based on the conservation of energy for the complete char and virgin layer. For the pyrolysis stage the conservation of energy for the char layer is given by Z dc d ddc _ 00g cg ðT s T pyr Þ: T pyr ¼ q_ 00net q_ 00c m rc cc T c dx rc cc (1) dt 0 dt The last term in Eq. (1) represents the heating of the pyrolysis gases when flowing from the point where they are produced, i.e. the pyrolysis front, to the solid surface where they are released in the atmosphere. In some models this term is omitted. Depending on the experiment and the stage of pyrolysis, the net incident heat flux consists of an external incident heat flux, a flame heat flux, reradiation, and convective cooling: qT c kc ¼ q_ 00net ¼ q_ 00ext þ q_ 00fl sðT 4s T 41 Þ hs ðT s T 1 Þ: (2) qx x¼0 For the heating stages and inert pyrolysis experiments, no flame is present and thus q_ 00fl ¼ 0: The description of the virgin zone depends on the location of the thermal penetration front ( ¼ influenced virgin thickness). If the front has not reached the rear surface yet, the virgin layer is treated as semi-infinite. The conservation of
ARTICLE IN PRESS E. Theuns et al. / Fire Safety Journal 40 (2005) 121–140
energy for the semi-infinite pyrolysis and cooling stage is then given by Z dc þdv d ddc ddc ddv T pyr rv cv þ rv cv T v dx þ rv cv T 0 ¼ q_ 00v : dt dc dt dt dt
125
(3)
The influenced virgin thickness (dv ) is always variable, while the char layer thickness (dc ) is only a variable for the pyrolysis stage. Once the thermal penetration front has reached the rear surface the virgin layer is treated as finite. There is no longer uninfluenced virgin material and the conservation equation for the virgin zone becomes Z L d ddc T pyr ¼ q_ 00v þ q_ 00bs : r cv T v dx þ rv cv (4) dt dc v dt In the integral model, the shape of the temperature profile in the char and virgin layer is prescribed. Typical profiles are quadratic [14,15,17] or exponential [13]. Here, the temperature distributions in the char and virgin layer are assumed quadratic and are given by T c ðx; tÞ ¼ a0 ðtÞ þ a1 ðtÞðdc xÞ þ a2 ðtÞðdc xÞ2 ; T v ðx; tÞ ¼ b0 ðtÞ þ b1 ðtÞðx dc Þ þ b2 ðtÞðx dc Þ2 :
ð5Þ
The unknown coefficients a0 ; a1 ; a2 ; b0 ; b1 and b2 are determined with conservation equation (1) and (3) or (1) and (4), and the boundary conditions, e.g. for the pyrolysis stages the temperature at the pyrolysis front is equal to the pyrolysis temperature. The finite and semi-infinite pyrolysis model can be reduced to a system of three differential equations in three unknowns. When yc and yv the integrated temperature for the char and virgin layer, respectively, are defined as Z dc Z dv þdc yc ¼ ðT c T pyr Þ dx; yv ¼ ðT v T pyr Þ dx (6) 0
dv
the set of equations for the semi-infinite pyrolysis model reads ! ddc 1 ¼ ðrv rc ÞDQpyr þ rv cv ðT pyr T 0 Þ dt dyc rv cv dyv 00 þ q_ net rc cc ; dt 2 dt ! dyc yc q_ 00 ¼ 3ac net 2 ; dt 2kc dc dyv 8av ðT 0 T pyr Þ2 ddc ¼ : þ 2ðT pyr T 0 Þ dt 3yv dt
ð7Þ
ARTICLE IN PRESS 126
E. Theuns et al. / Fire Safety Journal 40 (2005) 121–140
After transition to the finite stage, the thermal front no longer exists and the pyrolysis process is described by ! ddc 1
¼ dt rv rc DQpyr dyc dyv q_ 00net rc cc rv c v þ q_ 00bs ; dt dt ! 00 dyc yc q_ ¼ 3ac net 2 ; dt 2kc dc 00 dyv yv q_ bs ¼ 3av : ð8Þ dt 2kv ðL dc Þ2 The last set of system equations, Eqs. (8), is slightly different form the formulation in Moghtaderi et al. [14]. Indeed, instead of assuming a perfectly insulated or isothermal rear boundary, a convective boundary condition is applied. The system equations are more general because, by changing the convection coefficient, an insulated and isothermal rear boundary can be simulated as well. The rate of released pyrolysis gases is determined by the advancement of the pyrolysis front. By coupling the virgin and char solution, the mass loss can be calculated from _ 00g ¼ ðrv rc Þ m
ddc ¼ ðkv ðqT v =qxÞx¼dþc dt kc ðqT c =qxÞx¼dc Þ=DQpyr :
ð9Þ
2.2. Mathematical description of the interrupted pyrolysis stage When the external heat flux is constant (and high enough), all the virgin material will transform into char in a continuous way. But when the external heat flux is variable, it is possible that the pyrolysis process is interrupted. When the external heat flux decreases, less heat is provided at the pyrolysis front and thus less heat is available for the endothermic pyrolysis reactions. When an insufficient amount of heat is supplied, reactions will come to a standstill and a cooling stage is introduced. This is possible for both the semi-infinite and finite pyrolysis stage. When the external heat flux rises again the pyrolysis reactions can resume. Depending on the way of heating (external heat flux), the semi-infinite cooling stage can also result in the finite cooling stage. Note that the name cooling stage is chosen because this stage is most probable to occur when the pyrolysing solid is suddenly cooled down, although, this stage can occur as well for a reduced front surface heat flux. The two cooling stages are required for the simulation of real fires, as will be shown later. They are, however, not included in the formulations of the integral model for the pyrolysis of charring materials found in literature [12,14,17]. More details about the cooling stages can be found in Theuns et al. [24].
ARTICLE IN PRESS E. Theuns et al. / Fire Safety Journal 40 (2005) 121–140
127
3. Evaluation of the integral model For certain initial and boundary conditions an analytical solution exists of the integral model, see for example the heating of a semi-infinite solid with a constant net heat flux [25]. For most variable boundary conditions no analytical solutions exist, and therefore, the results of the integral model are often compared with experiments [12,14,17]. A drawback of the latter, is the lack of knowledge with respect to the modelling of the experiment, and particularly, of the incident heat flux. The differences between experimental data and simulation results can be due to the assumptions typical for the integral model (prescribed temperature profile), due to the representation of the experiment (modelling of incident heat flux, material properties), due to the physical model, due to numerical errors or due to a combination of these types of errors. 3.1. Moving grid model In this work the integral model is compared to a moving grid model [24]. Both models are based on the same physical model: the pyrolysis zone is represented by a single surface at pyrolysis temperature where virgin material is entirely transformed into char. The moving grid model, however, does not rely on any assumption for the temperature distribution and hence the moving grid acts as a numerical reference solution for the physical model. When the same material properties are used for the integral and moving grid model the evaluation of the integral model is rigorous, which is a major advantage. Indeed, when the quality of the integral model is judged on the basis of experiments in a reacting atmosphere, the heat feedback of the flame must be calculated. The simulation then becomes more complex and the modelling uncertainty is larger. Also, the determination of the material properties interferes with the evaluation of the model. In a real fire or fire test, the solid is heated mainly by radiation and convection. The incident heat flux can originate from an external heater, e.g. radiating cone in the Cone Calorimeter test, flames at the solid surface, a hot gas layer, etc. In this paper the model is run as a standalone model and therefore, the net incident heat flux is supplied as a boundary condition. The external heat flux is defined before the simulation and is not dependent on the mass release rate. When the solid model is coupled to a CFD model or zone model this net incident heat flux is computed by the code. All simulations for the rigorous numerical comparisons, are performed on particle board, for which the material properties are given in Table 1 (initial). The following cases are examined:
experimentally measured heat flux; sudden variation of external heat flux; variation of solid thickness; variation of rear boundary condition.
Finally, the integral and moving grid model are compared to inert pyrolysis experiments.
ARTICLE IN PRESS E. Theuns et al. / Fire Safety Journal 40 (2005) 121–140
128
Table 1 Initial guess and optimised material properties rc kv kc cv cc T pyr DQpyr rv (kg/m3) (kg/m3) (W/m K) (W/m K) (J/kg K) (J/kg K) (1C) (J/kg) Initial [29] 600 Optimised 670
60 117
0.36 0.15
0.23 0.38
2500 4655
375 280
8.7 105 0 5.1 105 171
1.0 1.0
60
0.015
measured heat flux (kW/m2)
moving grid integral arrhenius
m" (kg/m2.s)
2500 4463
v (—) cg (J/kg K)
0.010
0.005
40
20
0
0.000 0
300
600
time (s)
900
0
300
600
900
time (s)
Fig. 1. (Left) Mass flux of pyrolysis gases and (right) experimentally measured heat flux.
3.2. Experimentally measured heat flux In the first case the importance of the cooling stage is illustrated. An experimentally measured heat flux is imposed to the solid. The incident heat flux was measured in a vertical flame spread experiment on particle board [26], see Fig. 1 (right). The heat flux peaks after about 180 s, and afterwards drops drastically. This kind of heating has of course a strong influence on the pyrolysis behaviour. As a consequence of the drop in the heat flux, the cooling stage occurs. This is indeed observed in Fig. 1 (left): the release of pyrolysis gases drops to zero between 322 and 460 s. If the cooling stage were not included in the integral model, the experiment could only be adequately simulated during the first 322 s. The prediction of the mass release rate of pyrolysis gases is very good when compared to the moving grid model. Only immediately after the first peak (at 175 s), the mass release rate is underestimated with errors up to 24%. The time interval where the pyrolysis reactions are interrupted and the time of resuming of the reactions is predicted very well, with deviations less than 15 s. This indicates that the description of the cooling stage by a quadratic temperature profile in char and virgin layer is valid for a typical flame spread heat flux. In reality it is not probable that the mass release suddenly drops to zero. This can be illustrated by using for example the time-consuming Arrhenius model to predict the mass release. In this model the pyrolysis reactions are described by an Arrhenius
ARTICLE IN PRESS E. Theuns et al. / Fire Safety Journal 40 (2005) 121–140
129
law [4–6,27,28]: qr ¼ ðr rc ÞA expðE a =RTÞ; qt
(10)
where r is the local instantaneous solid density, and rc the final char density, A the pre-exponential factor and E a the activation energy. In the pyrolysis zone the solid consists of a local mixture of virgin and char material. Hence, the density (and material) changes in a continuous way, and not abrupt as in the integral and moving grid model. In this model the pyrolysis zone is finite and pyrolysis reactions proceed at lower temperatures as well, though at a lower rate. Therefore, the mass release rate in the Arrhenius model is still significant during the cooling stage of the integral model, as can be seen in Fig. 1. The temperature in the pyrolysis zone is still high enough to produce a significant mass release of pyrolysis gases. The weakness of the present integral model to predict the mass loss rate when the external heat flux becomes too low is associated with the assumption of a constant temperature on the pyrolysis front. When insufficient heat is provided, the front can no longer hold the pyrolysis temperature and comes to a standstill. An alternative can be found in the work of Galgano and Di Blasi [23]. They coupled the integral method with an Arrhenius reaction rate on the pyrolysis front for wood cylinders. The temperature of the front is no longer constant but is calculated during the simulation. The advancement of the front is given by the Arrhenius reaction rate. Such an approach is thought to improve the results for low external heat fluxes.
3.3. Sudden variation of the external heat flux In enclosed fires sometimes sudden variations of the external heat flux of the solid material can be noticed. For example, when flashover occurs, all exposed combustible items in the enclosure suddenly get involved in the combustion process. This phenomenon is here modelled as a sudden rise in the external heat flux some fixed time after the start the simulation (at 300 s). The heat flux suddenly drops after 400 s, which represents the effect of extinguishment of flames or the activation of sprinklers. In reality, the time of flashover, the extinguishment, and the activation of sprinklers will depend on the enclosure conditions. The imposed external heat flux is shown in Fig. 2 (right). The integral model has remarkable results, as shown in Fig. 2 (left). Immediately after the rise of the external heat flux, at 300 s, the mass flux of pyrolysis gases drops sharply instead of rising. The reason is the unrealistic, direct influence of the net incident heat flux on the complete temperature profile of the solid. As can be seen in Fig. 2, already after 1 s the integral model has recovered and the mass release is normal again. When the change in the external heat flux occurs at the start of pyrolysis instead at 300 s, the mass release behaves normal (not shown here). When the external heat flux suddenly drops, the integral model shows again a bad prediction just after the fall of the external heat flux, see Fig. 2 at 400 s. Again this is due to the unrealistic, direct influence of the net incident heat flux on the complete
ARTICLE IN PRESS E. Theuns et al. / Fire Safety Journal 40 (2005) 121–140
130 0.020
60 moving grid integral
q"ext (kW/m2)
m" (kg/m2.s)
0.015
0.010
40
20
0.005
0.000
0 100
200
300
time (s)
400
100
200
300
400
time (s)
Fig. 2. (Left) Mass flux of pyrolysis gases and (right) external heat flux.
temperature profile. There is a sudden short rise in the mass flux of pyrolysis gases, but less than a few seconds later, the mass flux is normal again. A concluding remark with respect to the occurrence of the erroneous peaks with the integral model is that as the shape of the temperature profile is prescribed (e.g. quadratic), the integral model can only be valid for types of heating which will result in that prescribed temperature profile. In flame spread simulations, the incident heat flux can suddenly rise due to for example ignition of the pyrolysis gases in the gas stage, but it can also suddenly fall due to for example lack of oxygen in the gas stage. This kind of heating can result in temperature profiles that are different from the quadratic ones. A simple example is given by solving the heating of a semi-infinite slab (no pyrolysis) with a quadratic temperature profile. At the thermal penetration depth, the temperature is equal to the initial slab temperature and the slope of the profile is zero because no conductive heat flows to the unheated zone. When the net incident heat flux suddenly changes from positive to negative, no realistic solution can be found for this problem. The temperature profile is indeed entirely determined by the net heat flux (kv ðqT=qxÞx¼0 ¼ q_ 00net ), the zero heat flux at the end of the thermal layer (kv ðqT=qxÞx¼dv ¼ 0), and the temperature at the thermal penetration depth. (T v jx¼dn ¼ T 0 ). These boundary conditions do not allow conservation of energy for a quadratic temperature profile and thus no acceptable solution exists for such event. 3.4. Variation of the solid thickness Next, the integral model is here used for a parameter study of the solid thickness. The thickness is varied from 1 mm, in fires typically thermally thin, to 50 mm, typically thermally thick. The results of the integral model are shown in Fig. 3, the results of the moving grid model in Fig. 4. For all these cases, the rear surface is perfectly insulated and the external heat flux is constant (q_ 00ext ¼ 50 kW=m2 ). For a small thickness the mass release is concentrated in one short peak and the solid is consumed very quickly. For a thick solid the mass release rate contains two
ARTICLE IN PRESS E. Theuns et al. / Fire Safety Journal 40 (2005) 121–140
131
0.04
L = 1 mm
m" (kg /m2.s)
0.03
Integral model
2 mm
0.02
5 mm 10 mm
0.01
20 mm 50 mm 0.00 10
100
1000
10000
time (s)
Fig. 3. Mass flux of pyrolysis gases for several thicknesses (integral model).
0.04
Moving grid
L = 1 mm
m" (kg/m2.s)
0.03
2 mm
0.02
5 mm 10 mm
0.01
20 mm 50 mm 0.00 10
100
1000
10000
time (s)
Fig. 4. Mass flux of pyrolysis gases for several thicknesses (moving grid model).
peaks. This first peak occurs at the start of pyrolysis, the second peak occurs near the end, when the pyrolysis front is reaching the rear surface of the solid. It is due to the back-effect, caused by the insulated rear surface [17,29]. When the pyrolysis front reaches the rear surface, continuously less energy flows from the pyrolysis front to the virgin material. The energy flowing to the virgin material, cannot leave the solid because of the insulated rear surface. Subsequently, more heat is available for the heat absorbing pyrolysis reactions at the pyrolysis front and hence the mass loss rate increases. When, for example, there is a strong convective boundary condition at the rear surface, the second peak does not occur. Also, for real thick materials this second peak is less pronounced. The growing char layer acts as an thermal insulator and causes the decrease of the mass release rate.
ARTICLE IN PRESS 132
E. Theuns et al. / Fire Safety Journal 40 (2005) 121–140
The first peaks in the mass release rates, see Figs. 3 and 4, coincide for thicknesses of 10, 20 and 50 mm. This is the region where the solid is calculated with the semiinfinite equations and no influence of the rear boundary or thickness is felt. Once the pyrolysis process is in the finite region, a thicker material produces less pyrolysis gases for a certain location of the pyrolysis front (for a insulated rear surface and constant incident heat flux). Indeed, for the thick material a significant amount of heat is ‘‘lost’’ at the pyrolysis front to the remaining cold virgin material which leaves less energy for the endothermic pyrolysis reactions. For a thinner material, on the other hand, the remaining virgin material has a higher temperature because there is less virgin material to be heated (no heat loss at rear surface). Consequently, the heat flow from the pyrolysis front to the virgin layer is lower, which leaves more energy for the pyrolysis reactions. As in this case the external heat flux is constant, the mass release rate is indeed lower, as can be seen in Fig. 4. The results for the integral model, in Fig. 3, again reveal a shortcoming. As long as the solid can be treated as semi-finite the mass release curves coincide (e.g. 10 and 20 mm overlap from 26 to 50 s). At 51 s the thinner solid (10 mm) is solved by finite theory, while the thicker (20 mm) is still solved by semi-infinite theory. The mass release rate for the thinner material should be higher, as explained above, but the results show the opposite. Due to the transition to the finite stage for the thinner material, the correct behaviour is no longer produced. Before transition, the thermal front is still advancing through the virgin layer and the pyrolysis process is described by Eqs. (7). After transition the thermal front no longer exists and the pyrolysis process is described by Eqs. (8). At 115 s the correct behaviour, i.e. higher mass release rate for the thinner material, is given again. The problem at transitions in the integral model manifests itself in a more severe error. Looking at the mass release rate of pyrolysis gases for a thickness of 5 and 10 mm in Fig. 3, we see that the integral model predicts the first mass release (often used as the condition for ignition) for the thicker material (10 mm) earlier than for the thin material (5 mm). It is clear that for a solid with insulated rear surface this cannot be the case. The correct behaviour is shown by the moving grid model in Fig. 4. The reason for the flaw in the integral model is the underprediction of the surface temperature after the transition from semi-infinite to finite stage. Immediately after the transition from the semi-infinite to the finite stage, the surface temperature rises slower than for a thicker solid still in the semi-infinite stage. As a consequence, it takes longer for the surface to reach the pyrolysis temperature. For the results with the integral model in Fig. 3, the thin solid (5 mm) is already treated as finite when the first mass release occurs, while the thicker solid (10 mm) is still treated as semi-infinite at the start of mass release. So, again it is the transition from semi-infinite to finite stage that is responsible for the error. 3.5. Variation of rear boundary condition In order to examine the back-effect, the rear boundary condition is varied. In all cases a convective heat loss is applied: q_ 00bs ¼ hbs ðT bs T 1 Þ:
(11)
ARTICLE IN PRESS E. Theuns et al. / Fire Safety Journal 40 (2005) 121–140
133
The convection coefficient hbs is varied from 0 to 20 W/m2 K. In these simulations the thickness of the solid was 20 mm and the external radiative heat flux 50 kW/m2. Results are shown for the integral model in Fig. 5, and for the moving grid model in Fig. 6. When the convection coefficient is raised, the second peak in the mass release rate becomes less pronounced and eventually disappears. In the integral model the pyrolysis process is described by the semi-infinite model for times shorter than 183 s. All mass release curves coincide because the rear boundary condition does not yet appear in the system equations for the integral model. Of course, it is assumed here that the ambient temperature at the rear surface is equal to the initial temperature. If the ambient temperature is variable or different from the initial temperature, the integral model in its present form, is no longer applicable for the semi-infinite stages. At the transition from semi-infinite to finite pyrolysis, it is expected that the higher the convection coefficient, the higher the heat loss, and thus the lower the mass release rate, as is shown by the moving grid in Fig. 6. For the integral model the opposite is true. Higher heat losses give temporarily (between 183 and 400 s) higher mass release rates, see Fig. 5, which is incorrect. The transition from semi-infinite to finite pyrolysis stage is responsible for this error. 3.6. Comparison to experiments The integral and moving grid model are compared with inert pyrolysis experiments on particle board [30]. Due to the inert atmosphere, there is no combustion of pyrolysis gases in the atmosphere. Modelling of flames and flame heat fluxes is avoided, and thus the incident heat flux is well known. Experimental results for particle board are available for two variable external radiant heat flux rates, where the rise was given by 16.6 kW/m2 s (exp16) and 33.3 kW/m2 s (exp33), and one constant radiant heat flux of 31 kW/m2 (exp31). The exposed surface of the samples loses energy by convection and reradiation. During the initial heating stage, thus at low temperatures, the convective cooling is 0.012
m" (kg/m2.s)
Integral model
Detail
0.008
h = 0 W/m2K 0.004
5 W/m2K 20 W/m2K 0.000 10
100
1000
10 W/m2K 10000
time (s)
Fig. 5. Mass flux of pyrolysis gases for several rear boundary conditions (integral model).
ARTICLE IN PRESS E. Theuns et al. / Fire Safety Journal 40 (2005) 121–140
134
0.012
m" (kg/m2.s)
Moving grid
Detail
0.008
h = 0 W/m2K 0.004
5 W/m2K
20 W/m2K
10 W/m2K
0.000 10
100
1000
10000
time (s)
Fig. 6. Mass flux of pyrolysis gases for several rear boundary conditions (moving grid model).
important. Therefore, in the simulations the convection coefficient hs for the natural convective cooling is determined by the correlation [31]: hs ¼ 0:54
l ðGr PrÞ0:25 L
(12)
instead of using a constant value. Hence, the convection coefficient is dependent on the surface temperature. A value of about 6 W/m2 K is typical for a surface temperature of 30 1C, while 16 W/m2 K for 500 1C. Evidently, at high surface temperatures the reradiation becomes more important. A major problem of pyrolysis models, is the determination of the (equivalent) material properties. Not all the material properties are known or can be found in literature. Sometimes some properties are estimated [28], but often the material properties are varied until good agreement with experiments is obtained [29]. The latter is done here as well. 3.6.1. Determination of material properties by error minimisation Optimisation of material properties is done only with the moving grid model. By adjusting the material properties, the simulated mass release rates of exp16 and exp31 are optimised simultaneously, leaving exp33 as a blind calculation. At each time interval (Dt) the simulated mass flux is compared with the experimental value. The errors for exp16 and exp33 are multiplied and the result defines the function that has been minimised: errorðrv ; rc ; kv ; kc ; cv ; cc ; T pyr ; Qpyr ; cg ; v Þ X M M X 00 00 _ 00g; exp 16 ði DtÞ _ 00g; exp 31 ði DtÞ; _ g; sim16 ði DtÞ m _ g; sim31 ði DtÞ m ¼ m m i¼1
i¼1
ð13Þ
ARTICLE IN PRESS E. Theuns et al. / Fire Safety Journal 40 (2005) 121–140
135
where M is large enough to capture the total mass release curve. The first term in the RHS represents the error of the simulation of exp16, the second the error of the simulation of exp31. The result of the two terms in the RHS are a measure for the difference between the calculated and the experimental result and have the dimensions of a mass flux (kg/m2 s). If these errors are multiplied by the time interval (Dt), a measure of the total mass release error (in kg/m2) is obtained. These errors can then be seen as the difference of the area spanned by the experimental and simulated mass flux curves. The minimisation is done by a downhill simplex method in multidimensions [32]. The method is started by N+1 linear independent points which define a simplex in an N-dimensional space. Here N is the number of material properties that are optimised. By reflecting, expanding or contracting points of the simplex, a minimum of the function is found. The iterative optimisation procedure is stopped when the difference between the lowest error and the highest error of the simplex falls below a small number (104). The vertices of the initial simplex were derived from the guessed material properties in Table 1. One vertex was constructed by contracting the initial guess to the origin by 25% in all directions of the N-space, while the other N vertices were created by expanding the initial guess by 25% for only one direction (or material property). This gives the N+1 linear independent points of the starting simplex. During the optimisation procedure the variables are bounded to allow only realistic values. In the calculations presented in this paper the bounds were chosen very wide. As a consequence, during the optimisation there were no restrictions for the material properties. The high number of variable parameters can be reduced by determining some parameters by independent tests. For example, the virgin density can easily be determined by weight and volume measurements. In our method the virgin density was deliberately left free, which allows to check whether the method predicts realistic values. Instead of using independent tests it is also possible to include extra information in the optimisation procedure, for example, surface or interior temperature measurements. The error in Eq. (13) should then include in addition of the mass loss error also a temperature error. 3.6.2. Discussion of the optimised properties Minimisation has been done for 10 variable material properties. The 10 variables are: the virgin and char density, the virgin and char conductivity, the virgin and char specific heat capacity, the pyrolysis temperature, the pyrolysis heat, the volatiles specific heat capacity and the virgin surface emissivity. The optimised results for the moving grid model are given in Table 1. The virgin density is realistic and resembles the value of de Ris and Yan [29]. The other properties strongly differ with the initial ones, e.g. the char properties and the pyrolysis temperature. Differences can be due to different types of particle board used in the work of de Ris and Yan [29] and in the experiments [30]. The high-optimised virgin heat capacity can be explained by the vaporisation of the water in the particle board. This phenomenon is not directly included in the pyrolysis model and will act as an energy sink in the virgin layer. By increasing the virgin heat capacity the effect of the vaporisation is partially included
ARTICLE IN PRESS E. Theuns et al. / Fire Safety Journal 40 (2005) 121–140
136
in the pyrolysis model. The optimised char heat capacity seems rather high as well, but the influence of the char heat capacity in the simulations is rather small. This has been tested by setting the char heat capacity at 2500 J/kg K. The effect of the char heat capacity is in such simulations small because it appears in the heat equation in combination with the char density (rc) which is small. Similar observations have been made by Chen et al. [33]. They noticed that the thermal capacity ratio R ¼ ðrc cc Þ=ðrv cv Þ does not significantly affect the mass loss and surface temperature histories. By bounding some variable parameters during the optimisation unrealistic results can be avoided. Here, the char heat capacity was kept constant (2500 J/kg K), but tight bounds could be used as well. The optimised material properties are only determined by information of the mass release experiments exp16 and exp31. They do well in predicting the three experiments, see Fig. 7. The experimental curves in Fig. 7 are the mean of the two experimental results in [30]. The good results for exp16 and exp31 could be expected as the material properties were optimised for those two experiments. Exp33 on the other hand was not included, and thus represents a blind prediction. In Fig. 7 it is noticed that the mass loss predictions for exp16 and exp33 are less satisfying for the start and end of pyrolysis. In the beginning of the experiment the incident heat flux is still low, for exp16 the heat flux is rising linearly by 1.66 103 W/cm2 s and for exp33 by 3.33 103 W/cm2 s. As a consequence, degradation reactions will proceed at relative low temperature. Models with a constant surface pyrolysis temperature cannot predict this early mass loss because the solid must first attain the pyrolysis temperature before mass can be released. In the work of Chen [33] it is stated that a surface pyrolysis front at constant temperature is valid for the range 10–100 kW/m2. In a real fire the incident heat flux to most materials will rise slowly and sometimes will stay low. For those low heat fluxes the integral and moving grid model may be less suitable and an Arrhenius rate equation for the pyrolysis reactions may be more appropriate. Still fire models that use a surface pyrolysis front at constant temperature in the solid phase, give good 0.009
Q=31kW/m2
m" (kg/m2.s)
Q=33kW/m2
exp inte mov
0.006
Q=16kW/m2 0.003
0.000 0
1000
2000
3000
time (s)
Fig. 7. Mass flux of pyrolysis gases for particle board exposed to several radiation intensities.
ARTICLE IN PRESS E. Theuns et al. / Fire Safety Journal 40 (2005) 121–140
137
results for fire and flame spread calculations [12,28,34]. By including experiments with low heat flux in the optimisation methodology, the material fire properties for the integral and moving grid model are thought to perform better for a wider range of incident heat fluxes (low and high). 3.6.3. Difference between integral and moving grid model The optimised material properties are used for predicting exp16, exp31 and exp33 with the integral and moving grid model. The differences between the integral and moving grid predictions are small for the three experiments, as can be seen in Fig. 7. The difference is only large at the beginning and end of pyrolysis. The reason is that the mass release rises or falls quickly at these points and the ignition and extinguishment times are not predicted to be equal. During pyrolysis it is noticed that at first, the integral model overpredicts the mass release rate, then gives the same value as the moving grid model, and for the remaining part underpredicts the mass release rate. When errors at the beginning and end of pyrolysis are not considered, the maximum difference between integral and moving grid model is about 1.5 103 kg/m2 s or about 15% of the peak mass release rate (exp31). The average difference is for the three experiments less than 1 104 kg/m2 s or about 1% of the peak mass release rate. The error ( ¼ difference between experiment and model) of the integral model is a bit higher than of the moving grid model, see Table 2. When we compare that error with the total experimental mass release, the difference between the two models can be neglected. The error made by approximating the temperature profile by a quadratic function is much smaller than the error made by the physical model. Therefore, we can say that the integral model is a suitable method for calculating the release of pyrolysis gases as described by the physical model.
4. Conclusions An integral model extended with a cooling stage has been evaluated. The cooling stages occur when an insufficient amount of energy is supplied to the solid and hence the pyrolysis reactions are interrupted. When the heat supply rises again, the pyrolysis reactions can resume. The occurrence of the cooling stage was illustrated for an experimental heat flux. The integral model predicts the interrupted pyrolysis well, although models with an Arrhenius rate equation, still predict a low mass release rate during the cooling stage of the integral model. Table 2 Error of mass release rate prediction
Integral (kg/m2) Moving (kg/m2) Total experimental mass release (kg/m2)
Exp16
Exp31
Exp33
1.631 1.593 5.280
0.528 0.496 4.530
1.490 1.391 5.536
ARTICLE IN PRESS 138
E. Theuns et al. / Fire Safety Journal 40 (2005) 121–140
A comparison of the integral model with the moving grid model is performed to test the validity of the prescribed quadratic temperature profile shape. It was illustrated that, due to this assumption the integral model is not always capable of giving correct or realistic results. Sudden variations in the boundary conditions are transmitted immediately and unrealistically high through the entire solid and lead to temporarily unrealistic results. On the other hand, it must be noted that the integral model recovers quickly from this error. A practically measured heat flux or slowly changing boundary conditions do not cause any problems. Next, it was illustrated that the integral model is not suited for parameter studies. At the transition from the semi-infinite to the finite stage, the correct behaviour is not always predicted. The most important example is the wrong prediction of the first release of pyrolysis gases when varying the solid thickness. Often the first mass release is taken as the condition for ignition time or ignition temperature and is as such very important. Also, when the rear boundary condition is varied, unrealistic results are noted. If the integral model is used as a pyrolysis model, the user should be aware of its limits and shortcomings. On the other hand, in the study of inert pyrolysis experiments the differences between the mass release rate of the integral and moving grid model are very small. Only temporarily, differences of 15% of the peak mass release rate are possible. The error of the physical model is much larger than the error introduced by solving the physical model with the integral model. This implies that the extended integral model is useful as pyrolysis model. Finally, it is noted that an automatic optimisation procedure was presented to determine the material properties for the solid pyrolysis experiments. It is based on fitting the simulated mass release curves to combined experimental results. This has been done by an automatic minimisation of the mass release error.
Acknowledgement The second author is a Postdoctoral Fellow of the Fund of Scientific Research— Flanders (Belgium) (FWO-Vlaanderen).
References [1] Olenick SM, Carpenter DJ. An updated international survey of computer models for fire and smoke. J Fire Protect Eng 2003;13:87–110. [2] Tewarson A. SFPE handbook of fire engineering, 2nd ed. Boston, MA: Society of Fire Protection Engineers; 1995. p. 3.53–3.124. [3] Kung HC. A mathematical model of wood pyrolysis. Combust Flame 1972;18:185–95. [4] Di Blasi C. Processes of flames spreading over the surface of charring solid fuels: effects of fuel thickness. Combust Flame 1994;97:225–39. [5] Staggs JEJ. A simple model of polymer pyrolysis including transport of volatiles. Fire Saf J 2000;34:69–80.
ARTICLE IN PRESS E. Theuns et al. / Fire Safety Journal 40 (2005) 121–140
139
[6] Yuen R, Casey R, De Vahl, Davis G, Leonardi E, Yeoh GH, Chandrasekaran V, Grubits SJ. A threedimensional mathematical model for the pyrolysis of wet wood. In: Hasemi Y, editor. Fire safety science—proceedings of the fifth international symposium, 1997. p.189–200. [7] Wichman IS, Atreya A. A simplified model for the pyrolysis of charring materials. Combust Flame 1987;68:231–47. [8] Staggs JEJ, Whiteley RH. Modelling the combustion of solid-phase fuels in cone calorimeter experiments. Fire Mater 1999;23:63–9. [9] Novozhilov V, Moghtaderi B, Fletcher DF, Kent JH. Computational fluid dynamics modelling of wood combustion. Fire Saf J 1996;27:69–84. [10] Di Blasi C. Modeling and simulation of combustion processes of charring and non-charring solid fuels. Prog. Energy Combust 1993;19:71–104. [11] Chen Y, Delichatsios MA, Motevalli V. Material pyrolysis properties, Part I: An integral model for one-dimensional transient pyrolysis of charring and non-charring materials. Combust Sci Technol 1993;88:309–28. [12] Delichatsios MM, Mathews MK, Delichatsios MA. An upward fire spread and growth simulation. In: Cox G, Langford B, editors. Fire safety science—proceedings of the third international symposium. New York: Elsevier Applied Science; 1991. p. 207–16. [13] Delichatsios MA, Saito K. Upward fire spread: key flammability properties, similarity solutions and flammability indices. In: Cox G and Langford B, editors. Fire safety science—proceedings of the third international symposium. New York: Elsevier Applied Science; 1991. p. 217–26. [14] Moghtaderi B, Novozhilov V, Fletcher D, Kent JH. An integral model for the transient pyrolysis of solid materials. Fire Mater 1997;21:7–16. [15] Moghtaderi B, Novozhilov V, Fletcher D, Kent JH. Mathematical modelling of the piloted igniton of wet wood using the heat-balance integral method. J Appl Fire Sci 1997;6(2):91–107. [16] Quintiere J, Iqbal N. An approximate integral model for the burning rate of a thermoplastic-like material. Fire Mater 1994;18:89–98. [17] Spearpoint MJ, Quintiere JG. Predicting the burning of wood using an integral model. Combust Flame 2000;123:308–24. [18] Staggs JEJ. Short communication: approximate solutions of the pyrolysis of char forming and filled polymers under thermally thick conditions. Fire Mater 2000;24:305–8. [19] Jia F, Galea ER, Patel MK. Numerical simulation of the mass loss process in pyrolizing char materials. Fire Mater 1999;23:71–8. [20] Staggs JEJ, Whiteley RH. Modelling the combustion of solid-phase fuels in cone calorimeter experiments. Fire Mater 1999;23:63–9. [21] Theuns E, Vierendeels J, Vandevelde P. A moving grid model for the pyrolysis of charring. Int J Numer Methods Heat Fluid Flow 2002;12(5):541–59. [22] Staggs JEJ. A discussion of modelling idealised ablative materials with particular reference to fire testing. Fire Saf J 1997;28:47–66. [23] Galgano A, Di Blasi C. Modeling wood degradation by the unreacted-core-shrinking approximation. Ind Eng Chem Res 2003;42:2101–11. [24] Theuns E, Merci B, Vierendeels J, Vandevelde P. Extension and evaluation of the integral model for transient pyrolysis of charring materials. Fire Mater, in press. [25] Carslaw HS, Jaeger JC. Conduction of heat in solids. London: Oxford University Press; 1959. p. 510. [26] Kokkala M, Mikkola E, Immonen M, Juutilainen H, Manner P, Parker WJ. Large-scale upward flame spread tests on wood products year. VTT Research Notes 1834, Technical Research Centre of Finland, Espoo, Finland, 1997. [27] Ritchie SJ, Steckler KD, Hamins A, Cleary TG, Yang JC, Kashiwagi T. The effect of sample size on the heat release of charring materials. In: Hasemi Y, editors. Fire safety science—proceedings of the fifth international symposium, 1997. p.177–88. [28] Yan Z, Holmstedt G. CFD and experimental studies of room fire growth on wall lining materials. Fire Saf J 1996;27:201–38. [29] de Ris JL, Yan Z. Modeling ignition and pyrolysis of charring fuels. Proceedings of the fifth international conference on fire and materials ‘98, San Antonio, TX, USA, 1998. p.111–21.
ARTICLE IN PRESS 140
E. Theuns et al. / Fire Safety Journal 40 (2005) 121–140
[30] Vovelle C, Akrich R, Delfau JL. Mass loss rate measurements on solid materials under radiative heating. Comb Sci Technol 1984;36:1–18. [31] Incropera FP, DeWitt DP. Fundamentals of heat and mass transfer. New York: Wiley; 1996. p. 886. [32] Press WH, Teukolsky SA, Vetterling WT, Flannery BP. Numerical recipes in Fortran 77: the art of scientific computation, vol. 1. Cambridge: Cambridge University Press; 1992. p. 992. [33] Chen Y, Motevalli V, Delichatsios MA. Material pyrolysis properties, Part II: Methodology for the derivation of pyrolysis properties for charring materials. Combust Sci Technol 1995;104(4–6):401–25. [34] Jia F, Galea ER, Patel MK. The numerical simulation of fire spread within a compartment using an integrated gas and solid phase combustion model. J Appl Fire Sci 1999;8(4):327–52.