- Email: [email protected]
Sensitivity analysis and parameter adjustment in a simplified physical wildland fire model
Advances in Engineering Software 90 (2015) 98–106
Contents lists available at ScienceDirect
Advances in Engineering Software journal homepage: www.elsevier.com/locate/advengsoft
Sensitivity analysis and parameter adjustment in a simplified physical wildland fire model D. Prieto b, M.I. Asensio a,b,∗, L. Ferragut a,b, J.M. Cascón a,c a
Instituto Universitario de Física Fundamental y Matemáticas, Casas del Parque 1, Salamanca, 37008, Spain Departamento de Matemática Aplicada, Universidad de Salamanca, Casas del Parque 2, Salamanca, 37008, Spain c Departamento de Economía e Historia Económica, Universidad de Salamanca, Edificio FES, Campus Miguel de Unamuno, Salamanca, 37007, Spain b
a r t i c l e
i n f o
Article history: Received 10 April 2015 Revised 16 June 2015 Accepted 2 August 2015
Keywords: Simplified physical model Fire model validation Laboratory fire experiments Global sensitivity analysis Parameter adjustment Wildland fire model
a b s t r a c t A global sensitivity analysis and parameter adjustment of a simplified physical fire model applied to a well measured experimental example is developed in order to validate the model. The fire model is a simplified physical 2D wildland fire model with some 3D effects that takes into account the wind, the slope of the orography, the fuel load and type, the moisture content, the energy lost in the vertical direction and the radiation from the flames. The simplicity of the model and the numerical techniques proposed allow very competitive computational times.
1. Introduction The aim of a wildland fire model is to take fuel, topography, weather conditions, ignition and fire suppression as input data and predict the spread of a fire in a time shorter than the real time scale of the fire spread so that it can be a useful tool by firefighters. In recent years, advances in computational power and the increase in the capabilities of spatial information technologies offer great potential for the effective modeling of the behavior of wildland fires. Usually most of reviews in fire behavior modeling classify fire models according to the nature of their construction as being physical, semiphysical or empirical [1–5]. According to the physical system modeled, wildland fire mathematical models may be classified as surface fires, crown fires, spotting and ground fires models [2]. Surface fire spread modeling has been one of the most important task carried out in wildland fire research. In the literature there are two regimes in the propagation of surface fires, wind-driven fires and plume-dominated fires [6]. The fire model we summarize here is a simplified physical model for wind-driven surface fires. It is based on the fundamental physics of combustion and fire spread, together with some appropriate assumptions, that by making use of efficient numerical and computa-
∗
Corresponding author. Tel.: +34923294400; fax: +34627443088. E-mail addresses: [email protected] (D. Prieto), [email protected] (M.I. Asensio), [email protected] (L. Ferragut), [email protected] (J.M. Cascón). http://dx.doi.org/10.1016/j.advengsoft.2015.08.001 0965-9978/© 2015 Published by Elsevier Ltd.
© 2015 Published by Elsevier Ltd.
tional tools, allow that the simulation time be shorter than the real time. This simplified physical model follows the ideas proposed in [7–9], and previous fire models developed by the authors in [10] where a two phase model is presented and in [11–13] where several modifications and ideas were developed. The current simplified physical model is described with detail in [14,15] where is applied to a real case. The model takes into account the radiation from the flames, the wind, the slope of the orography, the fuel load and type, the moisture content, and the energy lost in the vertical direction. This simplified physical model depends on several input variables and three model parameters. In order to validate the model, a global sensitivity analysis [16,17] of the input factors, input variables and model parameters, was carried out, allowing us to determine which input factors are the most influential on model outputs. In the analysis proposed, the model outputs are the rate of spread and the fire thickness. The results of the sensitivity analysis allows to select the most relevant parameters, so that, the parameter adjustment can be carried out in a more efficient way. The sensitivity analysis and the adjustment of parameters have been done using the set of laboratory experiments reported in [18]. The global sensitivity analysis and the parameter adjustment involve the solution of the problem many times and consequently an efficient numerical solution of the direct problem is needed. The numerical solution proposed attempts to reduce the computational time by defining the active nodes and making use of parallel computation.
D. Prieto et al. / Advances in Engineering Software 90 (2015) 98–106
In order to take into account some three-dimensional effects, specifically the radiation from the flames above the surface S, we shall consider the following three-dimensional domain,
Nomenclature Physical quantities T temperature (K) (T∞ reference temperature) E enthalpy (J m−2 ) M fuel load (kg m−2 ) (M0 initial fuel load) Mv moisture content (kg of water/kg of dry fuel) ρ density (kg m−3 ) C heat capacity (J K−1 kg−1 ) H natural convection coefficient (J s−1 m−2 K−1 ) R thermal radiation (J s−1 m−2 ) I radiation intensity (J s−1 m−2 sr−1 ) a mean absorption coefficient (m−1 ) v latent heat of evaporation (J kg−1 ) F flame height (m) half-life time of combustion (s) t1/2 [l], [t] spatial and temporal scales (m, s) Dimensionless quantities u dimensionless temperature (u = e
dimensionless enthalpy (e =
c
mass fraction of solid fuel (c
α r
λv uv up
γ
T −T∞ T∞ )
E MCT∞ ) = MM ) 0
dimensionless natural convection coefficient (α = H[t] MC )
dimensionless thermal radiation (r =
[t] MCT∞ R)
v) dimensionless latent heat of evaporation (λv = MCTv ∞ dimensionless temperature of water evaporation dimensionless temperature of solid fuel pyrolysis dimensionless coefficient of fuel mass variation by pyrolysis ( lnt 2[t] ) 1/2
v
β
99
dimensionless wind velocity correction factor
The outline of the paper is as follows. In Section 2 we recap the fire model, in Section 3 we explain the numerical scheme. In Section 4 we address the sensitivity analysis of the input factors, and in Section 5 the parameter adjustment, both for the same experimental example for which there are available data [18]. Finally, some conclusions are given in Section 6. 2. The fire model In this section we describe briefly the mathematical equations corresponding to the one phase model which follows the works developed by the authors in [11–13]. In particular, a model with local radiation is reported in [11], the idea of a multivalued operator for enthalpy is presented for the first time in [12], and a model with no local radiation is developed in [13]. Additionally, in [14] the model was tested modeling a real fire and in [15] we consider the model with data assimilation. Let d = [0, lx ] × [0, ly ] ⊂ R2 be a rectangle representing the projection of the surface S where the fire occurs, defined by the mapping
S : d −→ R3
(x, y) −→ (x, y, h(x, y)) where h(x, y) is a known function representing the orography of the surface S. We shall assume that vegetation can be represented by a given fuel load M (kg m−2 ) together with a moisture content Mv , (kg of water/kg of dry fuel). M and Mv are scalar functions defined on d. Moreover, we shall assume that the height, F, of the flames in a particular fire is bounded by δ .
D = {(x, y, z): x, y ∈ d, h(x, y) < z < h(x, y) + δ}. The equations governing the fire model are based on the energy and mass conservation equations on the surface, S, and the radiation equation in D. The non-dimensional simplified equations for the fire model are,
∂t e + β v · ∇ e + α u = r in S t ∈ (0, tmax ),
(1)
e ∈ G(u)
(2)
in S t ∈ (0, tmax ),
∂t c = −g(u)c in S t ∈ (0, tmax ).
(3)
We complete the problem with homogeneous Dirichlet boundary conditions and the following initial conditions,
u(x, y, 0) = u0 (x, y)
in S,
c(x, y, 0) = c0 (x, y)
in S.
The unknowns, e =
(4) (5)
E MCT∞ ,
dimensionless enthalpy, u =
T −T∞ T∞ ,
the
dimensionless temperature of the solid fuel and c = MM , the mass 0 fraction of solid fuel, are bidimensional variables defined in S × (0, tmax ), where tmax is the time t of study. The physical quantities E (J m−2 ), T (K) and M (kg m−2 ) are enthalpy, the temperature of solid fuel and the fuel load respectively, and C (J K−1 kg−1 ), the heat capacity of solid fuel, T∞ (K), a reference temperature, and M0 (kg m−2 ), the initial solid fuel load. The influence of moisture content is modeled through a multivalued operator representing the enthalpy. The dimensionless enthalpy e is thus an element of the following multivalued maximal monotone operator G [19],
G(u) =
⎧ ⎪ ⎨u
[uv , uv + λv ] ⎪ ⎩u + λv [u p + λv , ∞]
if if if if
u < uv , u = uv , uv < u < u p , u = up,
where uv and up are the dimensionless evaporation temperature of the water and the dimensionless pyrolysis temperature of the solid fuel, respectively. The quantity λv is the dimensionless evaporation heat related to the latent heat of evaporation v (J kg−1 ) and the fuel moisture content Mv (kg of water/kg of dry fuel)
λv =
Mv v . CT∞
It should be noted that in the burned zone the multivalued operator does not represent the physical phenomena exactly, since the water vapor is no longer in the porous medium. This drawback can be avoided by setting λv = 0 in the burned area. For further details about this multivalued operators see [12]. The convective term, β v · ∇ e, represents the energy convected by the gas pyrolyzed through the elementary control volume, where the surface wind velocity, v, is re-scaled by a correction factor β . In order to deeply understand the meaning of the parameter β we should explain how this one phase model (mainly solid phase model) can be considered a simplification of a gas-solid two phase model where we have retained the enthalpy transported in the gas phase. The transported enthalpy in the gas phase is
φρgCgV.∇ Tg where ρ g , Cg and Tg are the density, the specific heat capacity and the temperature of the gas phase, φ is the volume fraction occupied by the gas phase and V is the velocity of the gas phase. We estimate
100
D. Prieto et al. / Advances in Engineering Software 90 (2015) 98–106 Table 1 Input variables.
Fig. 1. Scheme of the domain, the flame, and the incident ray.
Input variables
Symbol
Units
Dimensionless wind velocity Height of the surface Fuel load Moisture content Reference temperature Flame temperature Pyrolysis temperature Latent heat of evaporation Half-life time of combustion Flame height Heat capacity
v h M Mv T∞ Tf Tp
− m kg/m2 kg of water/kg of dry fuel K K K J kg−1 s m J K−1 kg−1
Table 2 Model parameters.
the fraction of transported heat that must be taking into account in Eq. (1) as
β=
φρgCg Tg (1 − φ)ρsCs Ts
(6)
where ρ s , Cs and Ts are the density, the specific heat capacity and the temperature of the solid phase, and 1 − φ is the volume fraction occupied by the solid phase. Considering the fuel layer as being a porous medium and that the enthalpy of this porous medium is an average of the enthalpy in the gas and in the solid, the fraction φ /(1 − φ) represents the fraction of transported enthalpy that is retained by the solid. To estimate the value of β we take typical values of the heat capacity for the air (ρ g Cg ) and for the wood (ρ s Cs ). These values are ρgCg ≈ 103 J m−3 K−1 and ρsCs ≈ 106 J m−3 K−1 . Then assuming, e.g., φ = 0.8 we can estimate the value of β , taking for instance inside the flame the approximate value, Tg /Ts = 1500/500 and away from the φT
T
g flame Tg /Ts ≈ 1, β ≈ (1−φ) 10−3 = 4 × 10−3 Tgs . Then β ≈ 1.2 × 10−2 Ts inside the flame and β ≈ 4 × 10−3 away from the flame. The surface wind velocity v can be considered as given data or can be computed for example by means of the model developed by the authors in [20,21]. The coupling of this convection model with the wind model was detailed by the authors in [12]. The term α u represents the energy lost by natural convection in the vertical direction. The parameter α is related to physical quan−1 m−2 K−1 ) is the natural convection tities by α = H[t] MC , where H (J s coefficient and [t] is a time scale. The right hand side of Eq. (3) represents the loss of solid fuel due to combustion,
g(u) =
0
γ
if if
u < up, u = up,
i.e., the loss of solid fuel is null if the temperature is below the pyrolysis temperature, and it remains constant when the temperature of pyrolysis is reached. This constant value is inversely proportional to the solid fuel half-life time t1/2 of the combustion of each type of fuel,
γ=
ln 2[t] . t1/2
The right hand side of Eq. (1) describes the thermal radiation reaching the surface S from the flame above the layer,
r=
v
t1/2 F C
[t] R. MCT∞
R represents the incident energy at a point x = (x, y, h(x, y)) of the surface S due to radiation from the flame above the surface per unit time and per unit area, obtained by summing up the contribution of
Model parameters
Symbol
Units
Mean absorption coefficient Natural convection coefficient Correction factor of convective term
a H
m−1 J s−1 m−2 K−1 −
β
all directions ; that is
R(x) =
2π
ω=0
I(x, ) · N dω,
(7)
where ω is the solid angle and N is the unit outer vector normal to the surface S. We only considered the hemisphere above the fuel layer, and each contribution depends on the flame height F (see Fig. 1). I is the total radiation intensity, i.e., the integral over all wavelengths of the radiation energy passing through an area per unit time, per unit of projected area, and per unit of solid angle. The differential equation describing the total radiation intensity I at any position along a given path s in a gray medium may be written neglecting scattering as
dI + a(s)I(s) = a(s)Ib (s), ds
(8)
where Ib is the black body total radiation intensity and is governed by the Stefan–Boltzmann law, corresponding to the integral over all wavelengths of the emissive power of a black body
Ib =
σ 4 T , π
(9)
where σ = 5.6704 × 10−8 J s−1 m−2 K−4 is the Stefan–Boltzmann constant and the temperature T reaches the flame temperature denoted by Tf . a is the radiation absorption coefficient inside the flame. For further details about the radiation computation see [15]. For basic definitions concerning radiation theory see [22]. 2.1. Some considerations about the simplifications of the model In the following we differentiate between input variables and model parameters. Input variables are magnitudes that can be measured with more or less precision, and they are summarised in Table 1. Model parameters are unknowns, although their physical meaning can provide an approximate idea of their range, as explained above for β . The three parameters of our fire model are listed in Table 2. Note that the wind velocity v and the height of the surface h, as well as the initial value of the fuel load and temperature can considered as given data. In our experiment, we took the reference temperature T∞ as the ambient temperature and the latent heat of evaporation v = 2.25 × 106 (J kg−1 ), which is the latent heat of evaporation of water. The other input variables, moisture content, Mv , flame temperature, Tf , pyrolysis temperature, Tp , half-life time of combustion, t1/2 , flame height, F, and heat capacity, C, depend on each kind of fuel.
D. Prieto et al. / Advances in Engineering Software 90 (2015) 98–106
3. Numerical method
if b˜ <
We use a P1 finite element approximation on a regular mesh for spatial discretization and a Crank–Nicolson finite difference scheme time discretization of the Eqs. (1–3). It improves the semi-implicit Euler method considered in [15].
if
α t
1+ 1+
Let t = t n+1 − t n be a time step and let cn , en and un denote approximations at time step tn to the exact solution c, e and u, respectively. We consider a Crank–Nicolson finite difference scheme by discretization of the total derivative [23],
1 n+1 (e − e¯n ), t
if
1+
˜ uv < b< 1+
rn+1 + rn en+1 − e¯n un+1 + un +α = , t 2 2
(10)
en+1 ∈ G(un+1 ),
(11)
The multivalued operator in Eq. (11) is maximal monotone, and hence its resolvent Jμ = (Id + μG)−1 for any μ > 0 is a well defined univalued operator. Moreover, the Yosida approximation of G, Gμ =
c
n+1
=
en+1 = Gμ (un+1 + μen+1 ),
(13)
then
b˜ − λv 1+
u p + μv < b˜ < ∞
2
u˜ =
then
α t 2
u˜ = u p .
un+1 + un rn+1 + rn + t , 2 2
un+1 +un un+12+un cn .
1 − 2t g 1 + 2t g
(17)
(18)
2
The non-local radiation term r strongly depends on the temperature u and on the fuel mass c. Therefore, the problem given by Eqs. (15), (17) and (18) remains implicit. We propose a fixed-point iteration; more specifically, using un+1,0 = un , cn+1,0 = cn as the starting point, we compute for i ≥ 0
2
α t
en+1,i = ten − α t
c
n+1,i
=
e¯n − un +
rn+1,i−1 + rn
α
,
un+1,i−1 + un rn+1,i−1 + rn + t , 2 2
un+1,i−1 +un 2 un+1,i−1
cn , +un
1 − 2t g 1 + 2t g
(19)
(20)
(21)
2
where rn+1,i−1 represents the radiation term computed for the temperature un+1,i−1 and the fuel mass cn+1,i−1 . In practice, Eqs. (19)–(21) are solved only in a neighborhood of the fire flame front.
or
un+1 = Jμ (un+1 + μen+1 ).
(14)
un+1 +
2
α t
en+1 =
2
α t
e¯n − un +
rn+1 + rn
α
,
2 then taking μ = α in Eq. (14) we obtain t
un+1 = J2/α t
2 rn+1 + rn . e¯n − un + α t α
(15)
It remains to explain how to calculate un+1 in Eq. (15). That is, for n+1 n 2 e¯n − un + r α+r , computing u˜ = J2/α t (b) is equivaa given b = α t lent to solving
α t 2
Id + G u˜ b˜ =
Thus, u˜ is given by
α t 2
b.
3.2. Numerical solution of the radiation equation As radiation essentially comes from the flames we shall consider a physical model in which the gases produced by pyrolysis burn above the fuel layer, producing a flame over that layer. The flame emits radiation reaching the points ahead of the flame heating the non burned fuel around it, allowing the propagation of the fire. To evaluate radiation from flames we distinguish the case without wind (vertical flame with rectangular section) and the case with wind (tilted flame).
On the other hand, rearranging the terms in Eq. (10), we have
(16)
3.2.1. Vertical flame. In [15] we obtain for a rectangular flame the energy per unit area at the point x which is given by
R(x) =
aσ T f4
π
,
2
Finally, the problem to solve is given by Eq. (15) and the following rearranged expressions of Eq. (10) and (12),
Id−Jμ
μ is a Lipschitz operator and the inclusion in Eq. (11) is equivalent for all μ > 0 to the equation (see [24])
α t
then u˜ = uv ,
(12)
u v + μv
u p + μv
2
un+1,i = J2/α t
cn+1 − cn un+1 + un cn+1 + cn = −g . t 2 2
2
en+1 = t e¯n − α t
where e¯n = en ◦ X n , and X n (x) = X (x, t n+1 , t n ) ≈ x − β vt is the position at time tn of the particle that is at position x at time t n+1 . Thus, at each time step, we solve,
α t
1+
uv + μv < b˜
α t
α t
b˜
then u˜ =
uv
2
1+
<
2
α t
3.1. Time integration and numerical solution of the multivalued operator
α t
2
if
∂t e + β v · ∇ e ≈
1+
101
f
F 2 + ||x − x˜ ||2 − ||x − x˜ ||
˜ dA, ||x − x˜ || F 2 + ||x − x˜ ||2
(22)
where f is the burning surface, Tf is the temperature inside the flame, F is the flame height and a is the radiation absorption coefficient.
,
102
D. Prieto et al. / Advances in Engineering Software 90 (2015) 98–106
3.2.2. Tilted flame. For a more general tilted flame we consider the flame split in several layers, and we sum up the contribution of each layer. For a flame with L layers the energy per unit area in x due to layer = 0, . . . , L − 1 is [15]
R (x) =
aσ T f4
π
×
f
(l +1)2 F 2 +L2 ||x − x˜ ||2 − l 2 F 2 +L2 ||x − x˜ ||2 ˜ dA. l 2 F 2 + L2 ||x − x˜ ||2 (l +1)2 F 2 +L2 ||x− x˜ ||2
(23)
3.2.3. Discretization of the radiation term. As we have seen in the previous paragraphs for a rectangular flame (resp. for each layer of a tilted flame), the radiation term R (resp. R ) has the form
R(x) =
aσ T f4
π
˜ χ(x˜ ) f (x − x˜ )dA,
(24)
where χ is the characteristic function of the flames, i.e., χ(x) = 1 if T (x) = Tp and χ(x) = 0 otherwise, where Tp is the pyrolysis temperature. Representing the function f with the help of a finite element basis, for each node xi of the mesh, the radiation term can be written as
R(xi ) =
aσ T f4
π
Ri j A j χ(x j )
j
where R = (Ri j = f j (xi )), is the radiation matrix. In practice only the terms in the neighborhood of the flame are calculated. When wind is present, a radiation matrix R corresponding to each layer of the tilted flame is considered. The radiation matrix is computed only once (outside the time loop). 4. A global sensitivity analysis of the model A first step to validate any mathematical model is the global sensitivity analysis [16,17] of its factors (input variables and model parameters) (X1 , . . . , Xn ) in order to determine which are most influential in model outputs Y. The aim of sensitivity analysis is precisely to assess the importance of each factor in the values of the output variable of the model. This is crucial when both the input variables (available data) and model parameters are affected by uncertainties, and fire propagation models are a good example of this situation. Global sensitivity analysis also allows the significance of interactions between the input variables to be quantified, which is frequently of great importance, in particular for parameter adjustment. The sensitivity analysis is not only critical for model validation but also serves to guide future research efforts by focusing on the factors mainly responsible for the uncertainty in the model. Global sensitivity analysis is widely used in environmental models [25] and, particularly, in forest fire propagation models [26]. There are many ways to perform global sensitivity analysis [27]. One of the most effective way to implement global sensitivity analysis is based on the decomposition of the variance of the output Y in a unique additive series of the variance of independent functions of increasing dimensionality
V [Y (X1 , . . . , Xn )] =
i
+
V [Ai (Xi )] +
V [Bi j (Xi , X j )]
i< j
V [Bi jk (Xi , X j , Xk )] + · · ·
i< j
+ V [H (X1 , . . . , Xn )].
(25)
The terms V[Ai (Xi )] are the variances due to the main effects for each of the input variables, the terms V[Bij (Xi , Xj )] are the variances due to the pairwise interactions between all pairs of input variables, and
so on, up to the term considering the interaction effects of all n input variables together V [H (X1 , . . . , Xn )]. This decomposition permits the effect on the model output uncertainty of each input variable, Xi , to be assigned through the variances of the main effects, and of any group of input variables through higher-order variances (interaction effects). The sum of the main effect and all the variances in Eq. (25), where a certain variable Xi is present, is the total effect. A small value for the total effect will be enough to state that the variable analyzed is not important, while the main effects will be useful to quantify the direct effect of variables on the output factor, although small values of the main effects are not sufficient guarantee to reject a variable since it may have a strong influence on the output variable through the interactions. The variances calculated for the main and total effects divided by the total model output variance give the sensitivity indices of first-order and total respectively, which allow the input variables to be ranked quantitatively. Global sensitivity analysis assigns the uncertainty of the output variable to the uncertainty of input factors by the sampling approach of probability density functions (p.d.f.) associated with input factors: each of the input factors is varied simultaneously in a given interval following a given p.d.f., and the sensitivities are calculated over the whole range of variation of the input factors. Two widely used global sensitivity analysis methods are the Extended Fourier Amplitude Sensitivity Test (FAST) [28] and Sobol’s method [31], both included in the Simlab program (version 2.2) [32] used in this work. Both methods take into account the effect of the p.d.f. of each factor, consider the effect of the simultaneous variation of all factors, do not require that the model be additive or linear, and can treat grouped factors. Both Extended FAST and Sobol’s method provide first-order and total sensitivity indices. The aim of the sensitivity analysis of our fire propagation model is on the one side, to verify the adequacy of the model, and on the other side, to design the parameter adjustment as efficiently as possible. For both, sensitivity analysis and parameter adjustment, the most expensive part in terms of computational time is the model evaluation. With the goal of reduce the computational cost, the sensitivity analysis and the parameter adjustment of our fire propagation model was carried out in this paper for a set of laboratory experiments reported in [18]. Experimental results are always necessary to adjust and validate any type of model. In the case of wildland fire models, there are experimental results with very different scales: from real wildland fires to experimental fires in a laboratory. Although the aim of the model described here is to simulate real wildland fires, experimental work at the laboratory scale is significant: in a real wildland fire it is costly and extremely difficult to obtain measurements of many of the variables whereas more accurate measurements can be made in laboratory experiments. In addition, from the point of view of the computational cost of the numerical simulations to be performed, the experimental scale is much more affordable. The experiments reported in [18] were carried out in a low speed wind tunnel, on beds of Pinus pinaster needles. The wind tunnel has a bench of 3.00 m × 1.20 m where the last 2.00 m are movable between 0° and ±15°. The fuel bed consisted of a layer of Pinus pinaster needles of 2.0 m × 0.70 m on the center of the bench, with a depth of approximately 40 mm and an oven dry load of 0.5 kg m−2 . They supply measurements of flame angle, flame height, rate of spread and temperatures. We selected data from the experiments carried out for the following values of wind velocity: 0 m s−1 , 1 m s−1 , 2 m s−1 and 3 m s−1 ; and slope: 0°, 5°, 10° and 15° (up-slope). We did not analyze neither the slope nor the wind velocity as input variables since those are well-known data and obviously the rate of spread depends strongly on both factors. The initial fuel load was the same for all the experiments, M0 = 0.5 kg m−2 , as well as the ambient temperature, T∞ = 300 K. The pyrolysis temperature, Tp = 500 K, the heat capacity, C = −1 1600 J K−1 kg , and the half-life of combustion, t1/2 = 15 s are values
D. Prieto et al. / Advances in Engineering Software 90 (2015) 98–106
current in the literature. Taking into account data in [18] we propose for the flame temperature Tf a range from 900 K to 1200 K, and for the moisture content Mv from 0.08 to 0.2. The ranges selected for the model parameters in Table 2 are [0.5, 15] for the natural convection coefficient H, [0.5, 2] for the mean absorption coefficient a, and [0.001, 0.01] for the correction factor of convective term β . The propagation fire model was implemented in C++, using the API OpenMP [29] in order to take advantage of the multiprocessor platforms to reduce the computational time. It was compiled using the GNU Compiler (GCC) version 4.9. These simulations were computed on a Dell Precision T7610 workstation, equipped with 2 Intel Xeon E5 2670 v2 processors (10 cores, each one working at a frequency of 2.50 GHz; 20 cores altogether) and 64GB RAM. With the aims described above and the characteristics of the experimental example selected in mind, we consider the following input variables and parameters for the sensitivity analysis of our model,
103
Table 3 Output and input dimensionless factors. Input factor
Symbol
1
ROS √ C T F T Ht1/2 √ C T Ht1/2 M0C√ at1/2 C T F √ t1/2 C T
1
2 3 4 5 6 7
Mv t1/2 σ T f4
Lower range
Upper range
–
–
–
–
9.3750 × 10−3
2.8125 × 10−1
4500
18, 000
2.7778 × 10−5
1.1111 × 10−4
0
0.25
M0C T
3.1000
9.7975
β
0.001
0.01
M0 , C, t1/2 , T = Tp − T∞ , H, a, F, Mv , σ T f4 , β summarized in Tables 1 and 2, and two output variables, the rate of spread ROS, computed as the slope of the linear regression line used to fit the position of the fire front at each instant of time of the simulation, and the fire thickness FTH computed as the mean fire thickness measured at each instant of time of the simulation, both measured along the last 2.00 m. Prior to the sensitivity analysis, we organize the variables into sets of dimensionless parameters from the given variables using the Buckingham theorem [30], which provides a method for computing sets of dimensionless parameters, even if the form of the equation is still unknown. The aims of this dimensional analysis are to minimize the number of input factors to be considered and the interaction effects between them, and to get results as independent as possible of the dimension. We assume that the output variable Y (ROS or FTH) can be written as a function of input variables and model parameters (input factors), more precisely,
Y = f (M0 , C, t1/2 , T, H, a, F, Mv , σ T f4 , β) Then there are 11 variables involved, including Y. In our case the dependent variable is the output variable Y and we will chose as “repeating” variables M0 , C, t1/2 , T. The number of fundamental dimensions to be considered are 4, i.e., mass, length, time and temperature. So from the Buckingham theorem, we can group all the variables into 7 dimensionless groups in the general case. Then we can put the output variable Y as a function of 6 dimensionless parameters i or less, depending wether we consider or not the wind, water content, etc.,
ROS = FTH =
CTp f1 (2 , 3 , 4 , 5 , 6 , 7 )
CTp
t1/2
f2 (2 , 3 , 4 , 5 , 6 , 7 )
(26)
or in the case without wind,
ROS = FTH =
CTp f1 (2 , 3 , 4 , 5 , 6 ) CTp
f2 (2 , 3 , 4 , 5 , 6 )
[32]. Since there is no additional information, in a first approximation a uniform probability density function was selected for each of the dimensionless parameters. Both methods, Extended FAST and Sobol, provide very similar results, so for simplicity all figures correspond to the Extended FAST method. In Fig. 2, the first-order (black bars) and total (white bars) sensitivity indices computed with the Extended FAST method for the rate of spread ROS are represented for several slopes and no wind. The total indices for ROS are close to the first-order indices, so the correlation is minimal. Independently of the slope, without wind, the most relevant parameter is√the one associated to the mean absorption coefficient, 3 = at1/2 C T , of the radiation term, followed by t1/2 σ T 4
t1/2
Fig. 2. Sensitivity indices obtained with the Extended FAST method for ROS, no wind and several slopes.
(27)
A straightforward computation gives the expressions of the dimensionless factors i shown in Table 3, as well as their ranges computed from the values and ranges of the variables and parameters described above for the experimental example. The sensitivity analysis was carried out for the 6 (resp. 5) dimensionless groups in (26) (resp. (27)), with a sample of size N = 4998 for the Extended FAST method and a sample of size N = 7168 for Sobol’s method in order to ensure a good estimation of the sensitivity indices [16,17], using the Simlab (version 2.2) software for sensitivity analysis
the corresponding to the flame temperature, 6 = MCTf , and the moisture content, 5 = Mv . This shows that in low wind conditions the dominating mechanism of heat transfer in a wildland fire is radiation, as could be expected. For the second output variable, the fire thickness FTH, the correlation is higher as can be seen in Fig. 3 where the values of first-order (black bars) and total (white bars) sensitivity indices are different. The results again, do not depend on the slope, and the most relevant parameter is also√the one associated to the mean absorption coefficient, 3 = at1/2 C T , but closely followed by the one related to the Ht
natural convection coefficient 2 = M1/2 . This is consistent with the 0C fact the higher the loss by natural convection, the shorter the fire is.
104
D. Prieto et al. / Advances in Engineering Software 90 (2015) 98–106
Fig. 3. Sensitivity indices obtained with the Extended FAST method for FFW, no wind and several slopes.
Fig. 5. Sensitivity indices obtained with the Extended FAST method for FTH, slope 5°, several winds.
5. Parameter adjustment
Fig. 4. Sensitivity indices obtained with the Extended FAST method for ROS, slope 5°, several winds.
Again, another factor closely related with the radiation term, as the t1/2 σ T 4
one corresponding to the flame temperature, 6 = MCTf , is very important, as could be expected in a no wind fire. Under windy conditions, the correction factor of the convective term 7 = β becomes the most relevant input factor for ROS, independently of the slope, but depends strongly on the wind magnitude, since when wind is important wildland fire spread is dominated by convection. The order of importance of the other five input factors remains equal to conditions without wind, see Fig. 4. The correction factor of the convective term 7 = β is also a very important input factor for FTH, with a growing importance with wind magniHt
1/2 remains very tude, but the natural convection coefficient 2 = MC significant, see Fig. 5. The sensitivity analysis allows us to state that our model reflect the most important fire heat transfer in each situation and provides insights about how to devise the parameter adjustment.
The global sensitivity analysis developed allows us to conclude that in conditions without wind the most influential parameter in terms of the rate of spread is the mean absorption coefficient of the radiation term a, and in terms of the fire thickness we can not discard the natural convection coefficient H. Under windy conditions, the correction factor of the convective term β is critical. The experimental examples selected allow us to fix the values of the other input factors with more accuracy and focus on the adjustment of the unknown parameters a, H and β of the model, first a and H without wind, and then β under windy conditions, in order to validate the model proposed and test out the order of magnitude of the model parameters. For all the experiments we assume T∞ = 300 K, Tp = 500 K and −1 C = 1600 J K−1 kg . We distinguish the two different fuel moisture contents reported in [18]: 10%, for which Mv = 0.1 and T f = 1150 K, and 18%, for which Mv = 0.18 and T f = 1050 K. The flame height in our model is the effective flame height in terms of radiation and we approximate it by twice the video recorded flame height data reported in [18]. We adjust the parameters for the same experiments selected for the sensitivity analysis, that is, slope: 0°, 5°, 10° and 15° (up-slope); and wind velocity: 0 m s−1 , 1 m s−1 , 2 m s−1 and 3 m s−1 . An initial proposal for the objective function could be in terms of the two output variables used in the global sensitivity analysis but this has proved very weak due to two basics reasons, the lack of data about the fire thickness in [18] and the fact that the richer the objective function, the better the adjustment is. The objective function we proposed takes into account the time ti that the fire front needs to arrive to a distance xi from the initial fire line, from x1 = 1.10 cm to x18 = 2.80 cm every 10 cm. We do not consider neither first meter of the domain where there is no slope, nor the spurious data of the boundary. For the eight cases without wind, that is, two moisture content for each of the four slopes, the objective function to adjust a and H is,
J1 (a, H ) =
18 i=1
(tic − tie )2
(28)
D. Prieto et al. / Advances in Engineering Software 90 (2015) 98–106 Table 4 Input factors, experimental ROS [18], computed ROS and adjusted a and H for experiments without wind. Mv
Slope
ROSe (m s−1 )
ROSc (m s−1 )
a
H
0.1 0.1 0.1 0.1 0.18 0.18 0.18 0.18
0 5 10 15 0 5 10 15
0.00268 0.00211 0.00252 0.00286 0.00243 0.00183 0.00169 0.00217
0.00268 0.00210 0.00252 0.00285 0.00242 0.00183 0.00168 0.00216
1.41718 1.16641 1.33359 1.52523 1.68769 1.33359 1.30166 1.55244
7.25059 8.55806 8.15403 8.50949 6.94194 7.25059 8.24941 6.94194
Table 5 Input factors, experimental ROS [18], computed ROS and adjusted β for different wind conditions. The values of a and H are the corresponding to each case adjusted without wind. Mv 0.1
Slope 0
5
10
15
0.18
0
5
10
15
v (m s−1 )
ROSe (m s−1 )
ROSc (m s−1 )
β
1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3
0.01044 0.02077 0.04921 0.01198 0.02216 0.04817 0.01424 0.02522 0.06087 0.01424 0.02170 0.07231 0.00536 0.01065 0.01476 0.00719 0.01444 0.02093 0.00788 0.01358 0.02085 0.00899 0.01763 0.02652
0.01040 0.02067 0.04800 0.01199 0.02201 0.04628 0.01291 0.02311 0.046281 0.01367 0.02171 0.05721 0.00534 0.01069 0.01472 0.00715 0.01442 0.02093 0.00785 0.01356 0.02093 0.00893 0.01764 0.02650
0.0066667 0.0084028 0.0147917 0.0088194 0.0094444 0.0144444 0.0093750 0.0098611 0.0143056 0.0097917 0.0089583 0.0179697 0.0024306 0.0036806 0.0037500 0.0045833 0.0057639 0.0060417 0.0053472 0.0054167 0.0061111 0.0059722 0.0072222 0.0079167
Once a and H are adjusted, we adjust β for each of the eight cases and each of the three wind velocities
J2 (β) =
18
(tic − tie )2
(29)
i=1
The superscript c represents the data computed with our model for a given value of a and H (resp. β ) and the superscript e is the data provided in [18]. We used a direct search method to solve the optimization problems (28) and (29). These methods do not require any information about the gradient of the objective function, and search a set of points around the current point, looking for one where the value of the objective function is lower than the value at the current point. This kind of methods are widely used, also in fire modeling, [33]. In particular, for problem (28), we employed a Hooke and Jeeves method [34] consisting of an application of the cyclic coordinate method in addition to a pattern search. In addition, a golden section search method is implemented to minimize along each search direction. For problem (29), we use a Fibonacci method. The parameters should be positive, so it is a constrained optimization problem, for the convergence of direct search algorithms for constrained minimization see [35]. In Table 4 we summarize: (1) the rate of spread (ROSe ) measured in the experiments of [18] for the different slopes and moisture contents mentioned above, in conditions without wind; (2) the adjusted values of a and H computed with the described Hooke and Jeeves method performed in C++, with a mean computational cost of
105
748 s and one iteration of the Hooke and Jeeves algorithm; (3) the corresponding rate of spread (ROSc ) computed for these adjusted values with the model. As might be expected, the adjusted values of the mean absorption coefficient of the radiation term a do not vary significantly neither with the slope nor with the moisture content and flame temperature, always within the same order of magnitude. The Pearson variation coefficient of a, that is, the relative standard deviation (rate of standard deviation and mean expressed as a percentage) is 11%, so the variation of the adjusted values of a is small. Even though the lack of data about the fire thickness in [18] limits our options, we obtain a good adjustment of the natural convection coefficient H within the same order of magnitude, and a small Pearson variation coefficient 9%. In the case the wind is present (Table 5), we adjust the correction factor of the convective term β minimizing (29) with a Fibonacci method implemented in C++, with a mean computational cost less than 30 s. The variability of β is higher but the obtained values remain within an appropriate interval from the physical point of view (see (6) and its explanation). 6. Conclusions Some attempts have been made to validate the simplified physical fire model developed by authors. The model is simple enough to depend on a few parameters but sufficiently precise to reflect important phenomena for the evolution of a fire, taking into account the orography, the fuel load and type, the temperature, the humidity and the wind, and the most important mechanisms of heat transfer: convection and radiation. The model has been applied to a set of experimental examples for which data were available. The experimental scale allows to obtain measurements of many of the involved variables more feasible and accurate than in real fires, and reduces the computational cost of the numerical simulations to be performed. A global sensitivity analysis has been carried out over the simplified physical fire model at experimental scale, focus the analysis on the model parameters and the input variables difficult to be measured. The results obtained from this analysis allow us to conclude that the model properly reflects the importance of radiation in no wind conditions, and convection in windy conditions. The adjustment of the model parameters has been designed according the results of the global sensitivity analysis and available data. Each of the three model parameters a, H and β remains in the same order of magnitude independently of fuel moisture, slope and wind conditions, which gives validity to the model. Acknowledgments This work has been partially supported by Secretaría de Estado de Investigación, Desarrollo e Innovación and Centro para el Desarrollo Tecnológico Industrial of the Ministerio de Economía y Competitividad of the Spanish Government, grant contract: CGL2011-29396-C0302 and by Conserjería de Educación of the Junta de Castilla y León, grant contract: SA266A12-2. References [1] Perry GLW. Current approach to modelling the spread of wildland fire: a review. Prog Phys Geogr 1998;22(2):222–45. [2] Pastor E, Zarate L, Planas E, Arnaldos J. Mathematical models and calculation systems for the study of wildland fire behaviour. Prog Energy Combust Sci 2003;29(2):139–53. [3] Sullivan AL. Wildland surface fire spread modelling, 1990–2007. 1: physical and quasi-physical models. Int J Wildland Fire 2009;18:349–68. [4] Sullivan AL. Wildland surface fire spread modelling, 1990-2007. 2: empirical and quasi-empirical models. Int J Wildland Fire 2009;18:369–86. [5] Sullivan AL. Wildland surface fire spread modelling, 1990–2007. 3: simulation and mathematical analogue models. Int J Wildland Fire 2009c;18:387–403.
106
D. Prieto et al. / Advances in Engineering Software 90 (2015) 98–106
[6] Morvan D. Physical phenomena and length scales governing the behaviour of wildfires: a case for physical modelling. Fire Technol 2011;47:437–60. [7] Cox G. Combustion Fundamentals of Fire. London: Academic Press; 1995. [8] Weber RO. Analytical models of fire spread due to radiation. Combust Flame 1989;78:398–408. [9] Margerit J, Séro-Guillaume O. Modelling forest fires. Part II: reduction to twodimensional models and simulation of propagation. Int J Heat Mass Transf 2002;45:1723–37. [10] Asensio MI, Ferragut L, Simon J. Modelling of convective phenomena in forest fire. Rev R Acad Cien Serie A Mat 2001;95(1):13–27. [11] Asensio MI, Ferragut L. On a wildland fire model with radiation. Int J Numerical Methods Eng 2002;54:137–57. [12] Ferragut L, Asensio MI, Monedero S. A numerical method for solving convectionreaction-diffusion multivalued equations in fire spread modelling. Adv Eng Softw 2006;18:366–71. [13] Ferragut L, Asensio MI, Monedero S. Modelling radiation and moisture content in fire spread. Commun Numerical Methods Eng 2007;23:819–33. [14] Ferragut L, Asensio MI, Cascón JM, Prieto D. A simplified wildland fire model applied to a real case. Adv Differential Equ Appl, SEMA SIMAI Springer Ser 2015;4:155–67. [15] Ferragut L, Asensio MI, Cascón JM, Prieto D. A wildland fire physical model well suited to data assimilation. Pure Appl Geophys 2015;172:121–39. [16] Saltelli A, et al. The primer. England: John Wiley & Sons Ltd; 2008. [17] Saltelli A, et al. Sentitivity analysis in practice: a guide to assessing scientific models. England: John Wiley & Sons Ltd; 2004. [18] Mendes-Lopes JM, Ventura JMP, Amaral JMP. Flame characteristics, temperaturetime curves, and rate of spread in fires propagating in a bed of Pinus pinaster needles. Int J Wildland Fire 2003;12:64–84. [19] Brézis H. Operateurs maximaux monotones et semi-groupes de contractions dans les espaces de Hilbert. N-Holl Math Stud; 1973. [20] Asensio MI, Ferragut L, Simon J. A convection model for fire spread simulation. Appl Math Lett 2005;18:673–7.
[21] Ferragut L, Asensio MI, Simon J. High definition local adjustment model for 3D wind fields performing only 2D computations. Int J Numerical Methods Biomed Eng 2011;27:510–23. [22] Siegel R, Howell JR. Thermal radiation heat transfer. New York: McGraw-Hill,Inc; 1972. [23] Pironneau O. On the transport-diffusion algorithm and its applications to the Navier–Stokes equations. Numerische Math 1982;38:309–32. [24] Bermúdez A, Moreno C. Duality methods for solving variational inequalities. Comput Math Appl 1981;7:43–58. [25] Campolongo F, Saltelli A. Sensitivity analysis of an environmental model: an application of different analysis methods. Reliability Eng Syst Saf 1997;57:49–69. [26] Salvador R, Piñol J, Tarantola S, Pla E. Global sensitivity analysis and scale effects of a fire propagation model used over Mediterranean shrublands. Ecol Model 2001;136:175–89. [27] Hamby DM. A review of techniques for parameter sensitivity analysis of environmental models. Environ Monit Assess 1994;32:135–54. [28] Saltelli A, Tarantola S, Chank K. A quantitative, model independent method for global sensitivity analysis of model output. Technometrics 1999;41(1):39–56. [29] Chapman B, Jost G, Van Der Pas R. Using OpenMP: portable shared memory parallel programming. Cambridge: MIT Press; 2007. [30] Buckingham E. On physically similar systems; Illustrations of the use of dimensional equations. Phys Rev 1914;4(4):345–76. [31] Sobol IM. Sensitivity estimates for non-linear mathematical models. Math Model Comput Exp 1993;1(4):407–14. [32] 2009. SIMLAB Version 2.2 Simulation environment for uncertainty and sensitivity analysis, developed by the Joint Research Centre of the European Commission. [33] Chetehouna K, Séro-Guillaume O, Sochet I, Degiovanni A. On the experimental determination of flame front positions and of propagation parameters for a fire. Int J Therm Sci 2008;47:1148–57. [34] Mokhtar SB, Hanif DS, Shetty CM. Nonlinear programming: theory and algorithms. 3rd ed. Wiley; 2006. [35] Michael LR, Torczon V. Pattern search algorithms for bound constrained minimization. SIAM J Optim 1999;9(4):1082–99.
Contents lists available at ScienceDirect
Advances in Engineering Software journal homepage: www.elsevier.com/locate/advengsoft
Sensitivity analysis and parameter adjustment in a simplified physical wildland fire model D. Prieto b, M.I. Asensio a,b,∗, L. Ferragut a,b, J.M. Cascón a,c a
Instituto Universitario de Física Fundamental y Matemáticas, Casas del Parque 1, Salamanca, 37008, Spain Departamento de Matemática Aplicada, Universidad de Salamanca, Casas del Parque 2, Salamanca, 37008, Spain c Departamento de Economía e Historia Económica, Universidad de Salamanca, Edificio FES, Campus Miguel de Unamuno, Salamanca, 37007, Spain b
a r t i c l e
i n f o
Article history: Received 10 April 2015 Revised 16 June 2015 Accepted 2 August 2015
Keywords: Simplified physical model Fire model validation Laboratory fire experiments Global sensitivity analysis Parameter adjustment Wildland fire model
a b s t r a c t A global sensitivity analysis and parameter adjustment of a simplified physical fire model applied to a well measured experimental example is developed in order to validate the model. The fire model is a simplified physical 2D wildland fire model with some 3D effects that takes into account the wind, the slope of the orography, the fuel load and type, the moisture content, the energy lost in the vertical direction and the radiation from the flames. The simplicity of the model and the numerical techniques proposed allow very competitive computational times.
1. Introduction The aim of a wildland fire model is to take fuel, topography, weather conditions, ignition and fire suppression as input data and predict the spread of a fire in a time shorter than the real time scale of the fire spread so that it can be a useful tool by firefighters. In recent years, advances in computational power and the increase in the capabilities of spatial information technologies offer great potential for the effective modeling of the behavior of wildland fires. Usually most of reviews in fire behavior modeling classify fire models according to the nature of their construction as being physical, semiphysical or empirical [1–5]. According to the physical system modeled, wildland fire mathematical models may be classified as surface fires, crown fires, spotting and ground fires models [2]. Surface fire spread modeling has been one of the most important task carried out in wildland fire research. In the literature there are two regimes in the propagation of surface fires, wind-driven fires and plume-dominated fires [6]. The fire model we summarize here is a simplified physical model for wind-driven surface fires. It is based on the fundamental physics of combustion and fire spread, together with some appropriate assumptions, that by making use of efficient numerical and computa-
∗
Corresponding author. Tel.: +34923294400; fax: +34627443088. E-mail addresses: [email protected] (D. Prieto), [email protected] (M.I. Asensio), [email protected] (L. Ferragut), [email protected] (J.M. Cascón). http://dx.doi.org/10.1016/j.advengsoft.2015.08.001 0965-9978/© 2015 Published by Elsevier Ltd.
© 2015 Published by Elsevier Ltd.
tional tools, allow that the simulation time be shorter than the real time. This simplified physical model follows the ideas proposed in [7–9], and previous fire models developed by the authors in [10] where a two phase model is presented and in [11–13] where several modifications and ideas were developed. The current simplified physical model is described with detail in [14,15] where is applied to a real case. The model takes into account the radiation from the flames, the wind, the slope of the orography, the fuel load and type, the moisture content, and the energy lost in the vertical direction. This simplified physical model depends on several input variables and three model parameters. In order to validate the model, a global sensitivity analysis [16,17] of the input factors, input variables and model parameters, was carried out, allowing us to determine which input factors are the most influential on model outputs. In the analysis proposed, the model outputs are the rate of spread and the fire thickness. The results of the sensitivity analysis allows to select the most relevant parameters, so that, the parameter adjustment can be carried out in a more efficient way. The sensitivity analysis and the adjustment of parameters have been done using the set of laboratory experiments reported in [18]. The global sensitivity analysis and the parameter adjustment involve the solution of the problem many times and consequently an efficient numerical solution of the direct problem is needed. The numerical solution proposed attempts to reduce the computational time by defining the active nodes and making use of parallel computation.
D. Prieto et al. / Advances in Engineering Software 90 (2015) 98–106
In order to take into account some three-dimensional effects, specifically the radiation from the flames above the surface S, we shall consider the following three-dimensional domain,
Nomenclature Physical quantities T temperature (K) (T∞ reference temperature) E enthalpy (J m−2 ) M fuel load (kg m−2 ) (M0 initial fuel load) Mv moisture content (kg of water/kg of dry fuel) ρ density (kg m−3 ) C heat capacity (J K−1 kg−1 ) H natural convection coefficient (J s−1 m−2 K−1 ) R thermal radiation (J s−1 m−2 ) I radiation intensity (J s−1 m−2 sr−1 ) a mean absorption coefficient (m−1 ) v latent heat of evaporation (J kg−1 ) F flame height (m) half-life time of combustion (s) t1/2 [l], [t] spatial and temporal scales (m, s) Dimensionless quantities u dimensionless temperature (u = e
dimensionless enthalpy (e =
c
mass fraction of solid fuel (c
α r
λv uv up
γ
T −T∞ T∞ )
E MCT∞ ) = MM ) 0
dimensionless natural convection coefficient (α = H[t] MC )
dimensionless thermal radiation (r =
[t] MCT∞ R)
v) dimensionless latent heat of evaporation (λv = MCTv ∞ dimensionless temperature of water evaporation dimensionless temperature of solid fuel pyrolysis dimensionless coefficient of fuel mass variation by pyrolysis ( lnt 2[t] ) 1/2
v
β
99
dimensionless wind velocity correction factor
The outline of the paper is as follows. In Section 2 we recap the fire model, in Section 3 we explain the numerical scheme. In Section 4 we address the sensitivity analysis of the input factors, and in Section 5 the parameter adjustment, both for the same experimental example for which there are available data [18]. Finally, some conclusions are given in Section 6. 2. The fire model In this section we describe briefly the mathematical equations corresponding to the one phase model which follows the works developed by the authors in [11–13]. In particular, a model with local radiation is reported in [11], the idea of a multivalued operator for enthalpy is presented for the first time in [12], and a model with no local radiation is developed in [13]. Additionally, in [14] the model was tested modeling a real fire and in [15] we consider the model with data assimilation. Let d = [0, lx ] × [0, ly ] ⊂ R2 be a rectangle representing the projection of the surface S where the fire occurs, defined by the mapping
S : d −→ R3
(x, y) −→ (x, y, h(x, y)) where h(x, y) is a known function representing the orography of the surface S. We shall assume that vegetation can be represented by a given fuel load M (kg m−2 ) together with a moisture content Mv , (kg of water/kg of dry fuel). M and Mv are scalar functions defined on d. Moreover, we shall assume that the height, F, of the flames in a particular fire is bounded by δ .
D = {(x, y, z): x, y ∈ d, h(x, y) < z < h(x, y) + δ}. The equations governing the fire model are based on the energy and mass conservation equations on the surface, S, and the radiation equation in D. The non-dimensional simplified equations for the fire model are,
∂t e + β v · ∇ e + α u = r in S t ∈ (0, tmax ),
(1)
e ∈ G(u)
(2)
in S t ∈ (0, tmax ),
∂t c = −g(u)c in S t ∈ (0, tmax ).
(3)
We complete the problem with homogeneous Dirichlet boundary conditions and the following initial conditions,
u(x, y, 0) = u0 (x, y)
in S,
c(x, y, 0) = c0 (x, y)
in S.
The unknowns, e =
(4) (5)
E MCT∞ ,
dimensionless enthalpy, u =
T −T∞ T∞ ,
the
dimensionless temperature of the solid fuel and c = MM , the mass 0 fraction of solid fuel, are bidimensional variables defined in S × (0, tmax ), where tmax is the time t of study. The physical quantities E (J m−2 ), T (K) and M (kg m−2 ) are enthalpy, the temperature of solid fuel and the fuel load respectively, and C (J K−1 kg−1 ), the heat capacity of solid fuel, T∞ (K), a reference temperature, and M0 (kg m−2 ), the initial solid fuel load. The influence of moisture content is modeled through a multivalued operator representing the enthalpy. The dimensionless enthalpy e is thus an element of the following multivalued maximal monotone operator G [19],
G(u) =
⎧ ⎪ ⎨u
[uv , uv + λv ] ⎪ ⎩u + λv [u p + λv , ∞]
if if if if
u < uv , u = uv , uv < u < u p , u = up,
where uv and up are the dimensionless evaporation temperature of the water and the dimensionless pyrolysis temperature of the solid fuel, respectively. The quantity λv is the dimensionless evaporation heat related to the latent heat of evaporation v (J kg−1 ) and the fuel moisture content Mv (kg of water/kg of dry fuel)
λv =
Mv v . CT∞
It should be noted that in the burned zone the multivalued operator does not represent the physical phenomena exactly, since the water vapor is no longer in the porous medium. This drawback can be avoided by setting λv = 0 in the burned area. For further details about this multivalued operators see [12]. The convective term, β v · ∇ e, represents the energy convected by the gas pyrolyzed through the elementary control volume, where the surface wind velocity, v, is re-scaled by a correction factor β . In order to deeply understand the meaning of the parameter β we should explain how this one phase model (mainly solid phase model) can be considered a simplification of a gas-solid two phase model where we have retained the enthalpy transported in the gas phase. The transported enthalpy in the gas phase is
φρgCgV.∇ Tg where ρ g , Cg and Tg are the density, the specific heat capacity and the temperature of the gas phase, φ is the volume fraction occupied by the gas phase and V is the velocity of the gas phase. We estimate
100
D. Prieto et al. / Advances in Engineering Software 90 (2015) 98–106 Table 1 Input variables.
Fig. 1. Scheme of the domain, the flame, and the incident ray.
Input variables
Symbol
Units
Dimensionless wind velocity Height of the surface Fuel load Moisture content Reference temperature Flame temperature Pyrolysis temperature Latent heat of evaporation Half-life time of combustion Flame height Heat capacity
v h M Mv T∞ Tf Tp
− m kg/m2 kg of water/kg of dry fuel K K K J kg−1 s m J K−1 kg−1
Table 2 Model parameters.
the fraction of transported heat that must be taking into account in Eq. (1) as
β=
φρgCg Tg (1 − φ)ρsCs Ts
(6)
where ρ s , Cs and Ts are the density, the specific heat capacity and the temperature of the solid phase, and 1 − φ is the volume fraction occupied by the solid phase. Considering the fuel layer as being a porous medium and that the enthalpy of this porous medium is an average of the enthalpy in the gas and in the solid, the fraction φ /(1 − φ) represents the fraction of transported enthalpy that is retained by the solid. To estimate the value of β we take typical values of the heat capacity for the air (ρ g Cg ) and for the wood (ρ s Cs ). These values are ρgCg ≈ 103 J m−3 K−1 and ρsCs ≈ 106 J m−3 K−1 . Then assuming, e.g., φ = 0.8 we can estimate the value of β , taking for instance inside the flame the approximate value, Tg /Ts = 1500/500 and away from the φT
T
g flame Tg /Ts ≈ 1, β ≈ (1−φ) 10−3 = 4 × 10−3 Tgs . Then β ≈ 1.2 × 10−2 Ts inside the flame and β ≈ 4 × 10−3 away from the flame. The surface wind velocity v can be considered as given data or can be computed for example by means of the model developed by the authors in [20,21]. The coupling of this convection model with the wind model was detailed by the authors in [12]. The term α u represents the energy lost by natural convection in the vertical direction. The parameter α is related to physical quan−1 m−2 K−1 ) is the natural convection tities by α = H[t] MC , where H (J s coefficient and [t] is a time scale. The right hand side of Eq. (3) represents the loss of solid fuel due to combustion,
g(u) =
0
γ
if if
u < up, u = up,
i.e., the loss of solid fuel is null if the temperature is below the pyrolysis temperature, and it remains constant when the temperature of pyrolysis is reached. This constant value is inversely proportional to the solid fuel half-life time t1/2 of the combustion of each type of fuel,
γ=
ln 2[t] . t1/2
The right hand side of Eq. (1) describes the thermal radiation reaching the surface S from the flame above the layer,
r=
v
t1/2 F C
[t] R. MCT∞
R represents the incident energy at a point x = (x, y, h(x, y)) of the surface S due to radiation from the flame above the surface per unit time and per unit area, obtained by summing up the contribution of
Model parameters
Symbol
Units
Mean absorption coefficient Natural convection coefficient Correction factor of convective term
a H
m−1 J s−1 m−2 K−1 −
β
all directions ; that is
R(x) =
2π
ω=0
I(x, ) · N dω,
(7)
where ω is the solid angle and N is the unit outer vector normal to the surface S. We only considered the hemisphere above the fuel layer, and each contribution depends on the flame height F (see Fig. 1). I is the total radiation intensity, i.e., the integral over all wavelengths of the radiation energy passing through an area per unit time, per unit of projected area, and per unit of solid angle. The differential equation describing the total radiation intensity I at any position along a given path s in a gray medium may be written neglecting scattering as
dI + a(s)I(s) = a(s)Ib (s), ds
(8)
where Ib is the black body total radiation intensity and is governed by the Stefan–Boltzmann law, corresponding to the integral over all wavelengths of the emissive power of a black body
Ib =
σ 4 T , π
(9)
where σ = 5.6704 × 10−8 J s−1 m−2 K−4 is the Stefan–Boltzmann constant and the temperature T reaches the flame temperature denoted by Tf . a is the radiation absorption coefficient inside the flame. For further details about the radiation computation see [15]. For basic definitions concerning radiation theory see [22]. 2.1. Some considerations about the simplifications of the model In the following we differentiate between input variables and model parameters. Input variables are magnitudes that can be measured with more or less precision, and they are summarised in Table 1. Model parameters are unknowns, although their physical meaning can provide an approximate idea of their range, as explained above for β . The three parameters of our fire model are listed in Table 2. Note that the wind velocity v and the height of the surface h, as well as the initial value of the fuel load and temperature can considered as given data. In our experiment, we took the reference temperature T∞ as the ambient temperature and the latent heat of evaporation v = 2.25 × 106 (J kg−1 ), which is the latent heat of evaporation of water. The other input variables, moisture content, Mv , flame temperature, Tf , pyrolysis temperature, Tp , half-life time of combustion, t1/2 , flame height, F, and heat capacity, C, depend on each kind of fuel.
D. Prieto et al. / Advances in Engineering Software 90 (2015) 98–106
3. Numerical method
if b˜ <
We use a P1 finite element approximation on a regular mesh for spatial discretization and a Crank–Nicolson finite difference scheme time discretization of the Eqs. (1–3). It improves the semi-implicit Euler method considered in [15].
if
α t
1+ 1+
Let t = t n+1 − t n be a time step and let cn , en and un denote approximations at time step tn to the exact solution c, e and u, respectively. We consider a Crank–Nicolson finite difference scheme by discretization of the total derivative [23],
1 n+1 (e − e¯n ), t
if
1+
˜ uv < b< 1+
rn+1 + rn en+1 − e¯n un+1 + un +α = , t 2 2
(10)
en+1 ∈ G(un+1 ),
(11)
The multivalued operator in Eq. (11) is maximal monotone, and hence its resolvent Jμ = (Id + μG)−1 for any μ > 0 is a well defined univalued operator. Moreover, the Yosida approximation of G, Gμ =
c
n+1
=
en+1 = Gμ (un+1 + μen+1 ),
(13)
then
b˜ − λv 1+
u p + μv < b˜ < ∞
2
u˜ =
then
α t 2
u˜ = u p .
un+1 + un rn+1 + rn + t , 2 2
un+1 +un un+12+un cn .
1 − 2t g 1 + 2t g
(17)
(18)
2
The non-local radiation term r strongly depends on the temperature u and on the fuel mass c. Therefore, the problem given by Eqs. (15), (17) and (18) remains implicit. We propose a fixed-point iteration; more specifically, using un+1,0 = un , cn+1,0 = cn as the starting point, we compute for i ≥ 0
2
α t
en+1,i = ten − α t
c
n+1,i
=
e¯n − un +
rn+1,i−1 + rn
α
,
un+1,i−1 + un rn+1,i−1 + rn + t , 2 2
un+1,i−1 +un 2 un+1,i−1
cn , +un
1 − 2t g 1 + 2t g
(19)
(20)
(21)
2
where rn+1,i−1 represents the radiation term computed for the temperature un+1,i−1 and the fuel mass cn+1,i−1 . In practice, Eqs. (19)–(21) are solved only in a neighborhood of the fire flame front.
or
un+1 = Jμ (un+1 + μen+1 ).
(14)
un+1 +
2
α t
en+1 =
2
α t
e¯n − un +
rn+1 + rn
α
,
2 then taking μ = α in Eq. (14) we obtain t
un+1 = J2/α t
2 rn+1 + rn . e¯n − un + α t α
(15)
It remains to explain how to calculate un+1 in Eq. (15). That is, for n+1 n 2 e¯n − un + r α+r , computing u˜ = J2/α t (b) is equivaa given b = α t lent to solving
α t 2
Id + G u˜ b˜ =
Thus, u˜ is given by
α t 2
b.
3.2. Numerical solution of the radiation equation As radiation essentially comes from the flames we shall consider a physical model in which the gases produced by pyrolysis burn above the fuel layer, producing a flame over that layer. The flame emits radiation reaching the points ahead of the flame heating the non burned fuel around it, allowing the propagation of the fire. To evaluate radiation from flames we distinguish the case without wind (vertical flame with rectangular section) and the case with wind (tilted flame).
On the other hand, rearranging the terms in Eq. (10), we have
(16)
3.2.1. Vertical flame. In [15] we obtain for a rectangular flame the energy per unit area at the point x which is given by
R(x) =
aσ T f4
π
,
2
Finally, the problem to solve is given by Eq. (15) and the following rearranged expressions of Eq. (10) and (12),
Id−Jμ
μ is a Lipschitz operator and the inclusion in Eq. (11) is equivalent for all μ > 0 to the equation (see [24])
α t
then u˜ = uv ,
(12)
u v + μv
u p + μv
2
un+1,i = J2/α t
cn+1 − cn un+1 + un cn+1 + cn = −g . t 2 2
2
en+1 = t e¯n − α t
where e¯n = en ◦ X n , and X n (x) = X (x, t n+1 , t n ) ≈ x − β vt is the position at time tn of the particle that is at position x at time t n+1 . Thus, at each time step, we solve,
α t
1+
uv + μv < b˜
α t
α t
b˜
then u˜ =
uv
2
1+
<
2
α t
3.1. Time integration and numerical solution of the multivalued operator
α t
2
if
∂t e + β v · ∇ e ≈
1+
101
f
F 2 + ||x − x˜ ||2 − ||x − x˜ ||
˜ dA, ||x − x˜ || F 2 + ||x − x˜ ||2
(22)
where f is the burning surface, Tf is the temperature inside the flame, F is the flame height and a is the radiation absorption coefficient.
,
102
D. Prieto et al. / Advances in Engineering Software 90 (2015) 98–106
3.2.2. Tilted flame. For a more general tilted flame we consider the flame split in several layers, and we sum up the contribution of each layer. For a flame with L layers the energy per unit area in x due to layer = 0, . . . , L − 1 is [15]
R (x) =
aσ T f4
π
×
f
(l +1)2 F 2 +L2 ||x − x˜ ||2 − l 2 F 2 +L2 ||x − x˜ ||2 ˜ dA. l 2 F 2 + L2 ||x − x˜ ||2 (l +1)2 F 2 +L2 ||x− x˜ ||2
(23)
3.2.3. Discretization of the radiation term. As we have seen in the previous paragraphs for a rectangular flame (resp. for each layer of a tilted flame), the radiation term R (resp. R ) has the form
R(x) =
aσ T f4
π
˜ χ(x˜ ) f (x − x˜ )dA,
(24)
where χ is the characteristic function of the flames, i.e., χ(x) = 1 if T (x) = Tp and χ(x) = 0 otherwise, where Tp is the pyrolysis temperature. Representing the function f with the help of a finite element basis, for each node xi of the mesh, the radiation term can be written as
R(xi ) =
aσ T f4
π
Ri j A j χ(x j )
j
where R = (Ri j = f j (xi )), is the radiation matrix. In practice only the terms in the neighborhood of the flame are calculated. When wind is present, a radiation matrix R corresponding to each layer of the tilted flame is considered. The radiation matrix is computed only once (outside the time loop). 4. A global sensitivity analysis of the model A first step to validate any mathematical model is the global sensitivity analysis [16,17] of its factors (input variables and model parameters) (X1 , . . . , Xn ) in order to determine which are most influential in model outputs Y. The aim of sensitivity analysis is precisely to assess the importance of each factor in the values of the output variable of the model. This is crucial when both the input variables (available data) and model parameters are affected by uncertainties, and fire propagation models are a good example of this situation. Global sensitivity analysis also allows the significance of interactions between the input variables to be quantified, which is frequently of great importance, in particular for parameter adjustment. The sensitivity analysis is not only critical for model validation but also serves to guide future research efforts by focusing on the factors mainly responsible for the uncertainty in the model. Global sensitivity analysis is widely used in environmental models [25] and, particularly, in forest fire propagation models [26]. There are many ways to perform global sensitivity analysis [27]. One of the most effective way to implement global sensitivity analysis is based on the decomposition of the variance of the output Y in a unique additive series of the variance of independent functions of increasing dimensionality
V [Y (X1 , . . . , Xn )] =
i
+
V [Ai (Xi )] +
V [Bi j (Xi , X j )]
i< j
V [Bi jk (Xi , X j , Xk )] + · · ·
i< j
+ V [H (X1 , . . . , Xn )].
(25)
The terms V[Ai (Xi )] are the variances due to the main effects for each of the input variables, the terms V[Bij (Xi , Xj )] are the variances due to the pairwise interactions between all pairs of input variables, and
so on, up to the term considering the interaction effects of all n input variables together V [H (X1 , . . . , Xn )]. This decomposition permits the effect on the model output uncertainty of each input variable, Xi , to be assigned through the variances of the main effects, and of any group of input variables through higher-order variances (interaction effects). The sum of the main effect and all the variances in Eq. (25), where a certain variable Xi is present, is the total effect. A small value for the total effect will be enough to state that the variable analyzed is not important, while the main effects will be useful to quantify the direct effect of variables on the output factor, although small values of the main effects are not sufficient guarantee to reject a variable since it may have a strong influence on the output variable through the interactions. The variances calculated for the main and total effects divided by the total model output variance give the sensitivity indices of first-order and total respectively, which allow the input variables to be ranked quantitatively. Global sensitivity analysis assigns the uncertainty of the output variable to the uncertainty of input factors by the sampling approach of probability density functions (p.d.f.) associated with input factors: each of the input factors is varied simultaneously in a given interval following a given p.d.f., and the sensitivities are calculated over the whole range of variation of the input factors. Two widely used global sensitivity analysis methods are the Extended Fourier Amplitude Sensitivity Test (FAST) [28] and Sobol’s method [31], both included in the Simlab program (version 2.2) [32] used in this work. Both methods take into account the effect of the p.d.f. of each factor, consider the effect of the simultaneous variation of all factors, do not require that the model be additive or linear, and can treat grouped factors. Both Extended FAST and Sobol’s method provide first-order and total sensitivity indices. The aim of the sensitivity analysis of our fire propagation model is on the one side, to verify the adequacy of the model, and on the other side, to design the parameter adjustment as efficiently as possible. For both, sensitivity analysis and parameter adjustment, the most expensive part in terms of computational time is the model evaluation. With the goal of reduce the computational cost, the sensitivity analysis and the parameter adjustment of our fire propagation model was carried out in this paper for a set of laboratory experiments reported in [18]. Experimental results are always necessary to adjust and validate any type of model. In the case of wildland fire models, there are experimental results with very different scales: from real wildland fires to experimental fires in a laboratory. Although the aim of the model described here is to simulate real wildland fires, experimental work at the laboratory scale is significant: in a real wildland fire it is costly and extremely difficult to obtain measurements of many of the variables whereas more accurate measurements can be made in laboratory experiments. In addition, from the point of view of the computational cost of the numerical simulations to be performed, the experimental scale is much more affordable. The experiments reported in [18] were carried out in a low speed wind tunnel, on beds of Pinus pinaster needles. The wind tunnel has a bench of 3.00 m × 1.20 m where the last 2.00 m are movable between 0° and ±15°. The fuel bed consisted of a layer of Pinus pinaster needles of 2.0 m × 0.70 m on the center of the bench, with a depth of approximately 40 mm and an oven dry load of 0.5 kg m−2 . They supply measurements of flame angle, flame height, rate of spread and temperatures. We selected data from the experiments carried out for the following values of wind velocity: 0 m s−1 , 1 m s−1 , 2 m s−1 and 3 m s−1 ; and slope: 0°, 5°, 10° and 15° (up-slope). We did not analyze neither the slope nor the wind velocity as input variables since those are well-known data and obviously the rate of spread depends strongly on both factors. The initial fuel load was the same for all the experiments, M0 = 0.5 kg m−2 , as well as the ambient temperature, T∞ = 300 K. The pyrolysis temperature, Tp = 500 K, the heat capacity, C = −1 1600 J K−1 kg , and the half-life of combustion, t1/2 = 15 s are values
D. Prieto et al. / Advances in Engineering Software 90 (2015) 98–106
current in the literature. Taking into account data in [18] we propose for the flame temperature Tf a range from 900 K to 1200 K, and for the moisture content Mv from 0.08 to 0.2. The ranges selected for the model parameters in Table 2 are [0.5, 15] for the natural convection coefficient H, [0.5, 2] for the mean absorption coefficient a, and [0.001, 0.01] for the correction factor of convective term β . The propagation fire model was implemented in C++, using the API OpenMP [29] in order to take advantage of the multiprocessor platforms to reduce the computational time. It was compiled using the GNU Compiler (GCC) version 4.9. These simulations were computed on a Dell Precision T7610 workstation, equipped with 2 Intel Xeon E5 2670 v2 processors (10 cores, each one working at a frequency of 2.50 GHz; 20 cores altogether) and 64GB RAM. With the aims described above and the characteristics of the experimental example selected in mind, we consider the following input variables and parameters for the sensitivity analysis of our model,
103
Table 3 Output and input dimensionless factors. Input factor
Symbol
1
ROS √ C T F T Ht1/2 √ C T Ht1/2 M0C√ at1/2 C T F √ t1/2 C T
1
2 3 4 5 6 7
Mv t1/2 σ T f4
Lower range
Upper range
–
–
–
–
9.3750 × 10−3
2.8125 × 10−1
4500
18, 000
2.7778 × 10−5
1.1111 × 10−4
0
0.25
M0C T
3.1000
9.7975
β
0.001
0.01
M0 , C, t1/2 , T = Tp − T∞ , H, a, F, Mv , σ T f4 , β summarized in Tables 1 and 2, and two output variables, the rate of spread ROS, computed as the slope of the linear regression line used to fit the position of the fire front at each instant of time of the simulation, and the fire thickness FTH computed as the mean fire thickness measured at each instant of time of the simulation, both measured along the last 2.00 m. Prior to the sensitivity analysis, we organize the variables into sets of dimensionless parameters from the given variables using the Buckingham theorem [30], which provides a method for computing sets of dimensionless parameters, even if the form of the equation is still unknown. The aims of this dimensional analysis are to minimize the number of input factors to be considered and the interaction effects between them, and to get results as independent as possible of the dimension. We assume that the output variable Y (ROS or FTH) can be written as a function of input variables and model parameters (input factors), more precisely,
Y = f (M0 , C, t1/2 , T, H, a, F, Mv , σ T f4 , β) Then there are 11 variables involved, including Y. In our case the dependent variable is the output variable Y and we will chose as “repeating” variables M0 , C, t1/2 , T. The number of fundamental dimensions to be considered are 4, i.e., mass, length, time and temperature. So from the Buckingham theorem, we can group all the variables into 7 dimensionless groups in the general case. Then we can put the output variable Y as a function of 6 dimensionless parameters i or less, depending wether we consider or not the wind, water content, etc.,
ROS = FTH =
CTp f1 (2 , 3 , 4 , 5 , 6 , 7 )
CTp
t1/2
f2 (2 , 3 , 4 , 5 , 6 , 7 )
(26)
or in the case without wind,
ROS = FTH =
CTp f1 (2 , 3 , 4 , 5 , 6 ) CTp
f2 (2 , 3 , 4 , 5 , 6 )
[32]. Since there is no additional information, in a first approximation a uniform probability density function was selected for each of the dimensionless parameters. Both methods, Extended FAST and Sobol, provide very similar results, so for simplicity all figures correspond to the Extended FAST method. In Fig. 2, the first-order (black bars) and total (white bars) sensitivity indices computed with the Extended FAST method for the rate of spread ROS are represented for several slopes and no wind. The total indices for ROS are close to the first-order indices, so the correlation is minimal. Independently of the slope, without wind, the most relevant parameter is√the one associated to the mean absorption coefficient, 3 = at1/2 C T , of the radiation term, followed by t1/2 σ T 4
t1/2
Fig. 2. Sensitivity indices obtained with the Extended FAST method for ROS, no wind and several slopes.
(27)
A straightforward computation gives the expressions of the dimensionless factors i shown in Table 3, as well as their ranges computed from the values and ranges of the variables and parameters described above for the experimental example. The sensitivity analysis was carried out for the 6 (resp. 5) dimensionless groups in (26) (resp. (27)), with a sample of size N = 4998 for the Extended FAST method and a sample of size N = 7168 for Sobol’s method in order to ensure a good estimation of the sensitivity indices [16,17], using the Simlab (version 2.2) software for sensitivity analysis
the corresponding to the flame temperature, 6 = MCTf , and the moisture content, 5 = Mv . This shows that in low wind conditions the dominating mechanism of heat transfer in a wildland fire is radiation, as could be expected. For the second output variable, the fire thickness FTH, the correlation is higher as can be seen in Fig. 3 where the values of first-order (black bars) and total (white bars) sensitivity indices are different. The results again, do not depend on the slope, and the most relevant parameter is also√the one associated to the mean absorption coefficient, 3 = at1/2 C T , but closely followed by the one related to the Ht
natural convection coefficient 2 = M1/2 . This is consistent with the 0C fact the higher the loss by natural convection, the shorter the fire is.
104
D. Prieto et al. / Advances in Engineering Software 90 (2015) 98–106
Fig. 3. Sensitivity indices obtained with the Extended FAST method for FFW, no wind and several slopes.
Fig. 5. Sensitivity indices obtained with the Extended FAST method for FTH, slope 5°, several winds.
5. Parameter adjustment
Fig. 4. Sensitivity indices obtained with the Extended FAST method for ROS, slope 5°, several winds.
Again, another factor closely related with the radiation term, as the t1/2 σ T 4
one corresponding to the flame temperature, 6 = MCTf , is very important, as could be expected in a no wind fire. Under windy conditions, the correction factor of the convective term 7 = β becomes the most relevant input factor for ROS, independently of the slope, but depends strongly on the wind magnitude, since when wind is important wildland fire spread is dominated by convection. The order of importance of the other five input factors remains equal to conditions without wind, see Fig. 4. The correction factor of the convective term 7 = β is also a very important input factor for FTH, with a growing importance with wind magniHt
1/2 remains very tude, but the natural convection coefficient 2 = MC significant, see Fig. 5. The sensitivity analysis allows us to state that our model reflect the most important fire heat transfer in each situation and provides insights about how to devise the parameter adjustment.
The global sensitivity analysis developed allows us to conclude that in conditions without wind the most influential parameter in terms of the rate of spread is the mean absorption coefficient of the radiation term a, and in terms of the fire thickness we can not discard the natural convection coefficient H. Under windy conditions, the correction factor of the convective term β is critical. The experimental examples selected allow us to fix the values of the other input factors with more accuracy and focus on the adjustment of the unknown parameters a, H and β of the model, first a and H without wind, and then β under windy conditions, in order to validate the model proposed and test out the order of magnitude of the model parameters. For all the experiments we assume T∞ = 300 K, Tp = 500 K and −1 C = 1600 J K−1 kg . We distinguish the two different fuel moisture contents reported in [18]: 10%, for which Mv = 0.1 and T f = 1150 K, and 18%, for which Mv = 0.18 and T f = 1050 K. The flame height in our model is the effective flame height in terms of radiation and we approximate it by twice the video recorded flame height data reported in [18]. We adjust the parameters for the same experiments selected for the sensitivity analysis, that is, slope: 0°, 5°, 10° and 15° (up-slope); and wind velocity: 0 m s−1 , 1 m s−1 , 2 m s−1 and 3 m s−1 . An initial proposal for the objective function could be in terms of the two output variables used in the global sensitivity analysis but this has proved very weak due to two basics reasons, the lack of data about the fire thickness in [18] and the fact that the richer the objective function, the better the adjustment is. The objective function we proposed takes into account the time ti that the fire front needs to arrive to a distance xi from the initial fire line, from x1 = 1.10 cm to x18 = 2.80 cm every 10 cm. We do not consider neither first meter of the domain where there is no slope, nor the spurious data of the boundary. For the eight cases without wind, that is, two moisture content for each of the four slopes, the objective function to adjust a and H is,
J1 (a, H ) =
18 i=1
(tic − tie )2
(28)
D. Prieto et al. / Advances in Engineering Software 90 (2015) 98–106 Table 4 Input factors, experimental ROS [18], computed ROS and adjusted a and H for experiments without wind. Mv
Slope
ROSe (m s−1 )
ROSc (m s−1 )
a
H
0.1 0.1 0.1 0.1 0.18 0.18 0.18 0.18
0 5 10 15 0 5 10 15
0.00268 0.00211 0.00252 0.00286 0.00243 0.00183 0.00169 0.00217
0.00268 0.00210 0.00252 0.00285 0.00242 0.00183 0.00168 0.00216
1.41718 1.16641 1.33359 1.52523 1.68769 1.33359 1.30166 1.55244
7.25059 8.55806 8.15403 8.50949 6.94194 7.25059 8.24941 6.94194
Table 5 Input factors, experimental ROS [18], computed ROS and adjusted β for different wind conditions. The values of a and H are the corresponding to each case adjusted without wind. Mv 0.1
Slope 0
5
10
15
0.18
0
5
10
15
v (m s−1 )
ROSe (m s−1 )
ROSc (m s−1 )
β
1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3 1 2 3
0.01044 0.02077 0.04921 0.01198 0.02216 0.04817 0.01424 0.02522 0.06087 0.01424 0.02170 0.07231 0.00536 0.01065 0.01476 0.00719 0.01444 0.02093 0.00788 0.01358 0.02085 0.00899 0.01763 0.02652
0.01040 0.02067 0.04800 0.01199 0.02201 0.04628 0.01291 0.02311 0.046281 0.01367 0.02171 0.05721 0.00534 0.01069 0.01472 0.00715 0.01442 0.02093 0.00785 0.01356 0.02093 0.00893 0.01764 0.02650
0.0066667 0.0084028 0.0147917 0.0088194 0.0094444 0.0144444 0.0093750 0.0098611 0.0143056 0.0097917 0.0089583 0.0179697 0.0024306 0.0036806 0.0037500 0.0045833 0.0057639 0.0060417 0.0053472 0.0054167 0.0061111 0.0059722 0.0072222 0.0079167
Once a and H are adjusted, we adjust β for each of the eight cases and each of the three wind velocities
J2 (β) =
18
(tic − tie )2
(29)
i=1
The superscript c represents the data computed with our model for a given value of a and H (resp. β ) and the superscript e is the data provided in [18]. We used a direct search method to solve the optimization problems (28) and (29). These methods do not require any information about the gradient of the objective function, and search a set of points around the current point, looking for one where the value of the objective function is lower than the value at the current point. This kind of methods are widely used, also in fire modeling, [33]. In particular, for problem (28), we employed a Hooke and Jeeves method [34] consisting of an application of the cyclic coordinate method in addition to a pattern search. In addition, a golden section search method is implemented to minimize along each search direction. For problem (29), we use a Fibonacci method. The parameters should be positive, so it is a constrained optimization problem, for the convergence of direct search algorithms for constrained minimization see [35]. In Table 4 we summarize: (1) the rate of spread (ROSe ) measured in the experiments of [18] for the different slopes and moisture contents mentioned above, in conditions without wind; (2) the adjusted values of a and H computed with the described Hooke and Jeeves method performed in C++, with a mean computational cost of
105
748 s and one iteration of the Hooke and Jeeves algorithm; (3) the corresponding rate of spread (ROSc ) computed for these adjusted values with the model. As might be expected, the adjusted values of the mean absorption coefficient of the radiation term a do not vary significantly neither with the slope nor with the moisture content and flame temperature, always within the same order of magnitude. The Pearson variation coefficient of a, that is, the relative standard deviation (rate of standard deviation and mean expressed as a percentage) is 11%, so the variation of the adjusted values of a is small. Even though the lack of data about the fire thickness in [18] limits our options, we obtain a good adjustment of the natural convection coefficient H within the same order of magnitude, and a small Pearson variation coefficient 9%. In the case the wind is present (Table 5), we adjust the correction factor of the convective term β minimizing (29) with a Fibonacci method implemented in C++, with a mean computational cost less than 30 s. The variability of β is higher but the obtained values remain within an appropriate interval from the physical point of view (see (6) and its explanation). 6. Conclusions Some attempts have been made to validate the simplified physical fire model developed by authors. The model is simple enough to depend on a few parameters but sufficiently precise to reflect important phenomena for the evolution of a fire, taking into account the orography, the fuel load and type, the temperature, the humidity and the wind, and the most important mechanisms of heat transfer: convection and radiation. The model has been applied to a set of experimental examples for which data were available. The experimental scale allows to obtain measurements of many of the involved variables more feasible and accurate than in real fires, and reduces the computational cost of the numerical simulations to be performed. A global sensitivity analysis has been carried out over the simplified physical fire model at experimental scale, focus the analysis on the model parameters and the input variables difficult to be measured. The results obtained from this analysis allow us to conclude that the model properly reflects the importance of radiation in no wind conditions, and convection in windy conditions. The adjustment of the model parameters has been designed according the results of the global sensitivity analysis and available data. Each of the three model parameters a, H and β remains in the same order of magnitude independently of fuel moisture, slope and wind conditions, which gives validity to the model. Acknowledgments This work has been partially supported by Secretaría de Estado de Investigación, Desarrollo e Innovación and Centro para el Desarrollo Tecnológico Industrial of the Ministerio de Economía y Competitividad of the Spanish Government, grant contract: CGL2011-29396-C0302 and by Conserjería de Educación of the Junta de Castilla y León, grant contract: SA266A12-2. References [1] Perry GLW. Current approach to modelling the spread of wildland fire: a review. Prog Phys Geogr 1998;22(2):222–45. [2] Pastor E, Zarate L, Planas E, Arnaldos J. Mathematical models and calculation systems for the study of wildland fire behaviour. Prog Energy Combust Sci 2003;29(2):139–53. [3] Sullivan AL. Wildland surface fire spread modelling, 1990–2007. 1: physical and quasi-physical models. Int J Wildland Fire 2009;18:349–68. [4] Sullivan AL. Wildland surface fire spread modelling, 1990-2007. 2: empirical and quasi-empirical models. Int J Wildland Fire 2009;18:369–86. [5] Sullivan AL. Wildland surface fire spread modelling, 1990–2007. 3: simulation and mathematical analogue models. Int J Wildland Fire 2009c;18:387–403.
106
D. Prieto et al. / Advances in Engineering Software 90 (2015) 98–106
[6] Morvan D. Physical phenomena and length scales governing the behaviour of wildfires: a case for physical modelling. Fire Technol 2011;47:437–60. [7] Cox G. Combustion Fundamentals of Fire. London: Academic Press; 1995. [8] Weber RO. Analytical models of fire spread due to radiation. Combust Flame 1989;78:398–408. [9] Margerit J, Séro-Guillaume O. Modelling forest fires. Part II: reduction to twodimensional models and simulation of propagation. Int J Heat Mass Transf 2002;45:1723–37. [10] Asensio MI, Ferragut L, Simon J. Modelling of convective phenomena in forest fire. Rev R Acad Cien Serie A Mat 2001;95(1):13–27. [11] Asensio MI, Ferragut L. On a wildland fire model with radiation. Int J Numerical Methods Eng 2002;54:137–57. [12] Ferragut L, Asensio MI, Monedero S. A numerical method for solving convectionreaction-diffusion multivalued equations in fire spread modelling. Adv Eng Softw 2006;18:366–71. [13] Ferragut L, Asensio MI, Monedero S. Modelling radiation and moisture content in fire spread. Commun Numerical Methods Eng 2007;23:819–33. [14] Ferragut L, Asensio MI, Cascón JM, Prieto D. A simplified wildland fire model applied to a real case. Adv Differential Equ Appl, SEMA SIMAI Springer Ser 2015;4:155–67. [15] Ferragut L, Asensio MI, Cascón JM, Prieto D. A wildland fire physical model well suited to data assimilation. Pure Appl Geophys 2015;172:121–39. [16] Saltelli A, et al. The primer. England: John Wiley & Sons Ltd; 2008. [17] Saltelli A, et al. Sentitivity analysis in practice: a guide to assessing scientific models. England: John Wiley & Sons Ltd; 2004. [18] Mendes-Lopes JM, Ventura JMP, Amaral JMP. Flame characteristics, temperaturetime curves, and rate of spread in fires propagating in a bed of Pinus pinaster needles. Int J Wildland Fire 2003;12:64–84. [19] Brézis H. Operateurs maximaux monotones et semi-groupes de contractions dans les espaces de Hilbert. N-Holl Math Stud; 1973. [20] Asensio MI, Ferragut L, Simon J. A convection model for fire spread simulation. Appl Math Lett 2005;18:673–7.
[21] Ferragut L, Asensio MI, Simon J. High definition local adjustment model for 3D wind fields performing only 2D computations. Int J Numerical Methods Biomed Eng 2011;27:510–23. [22] Siegel R, Howell JR. Thermal radiation heat transfer. New York: McGraw-Hill,Inc; 1972. [23] Pironneau O. On the transport-diffusion algorithm and its applications to the Navier–Stokes equations. Numerische Math 1982;38:309–32. [24] Bermúdez A, Moreno C. Duality methods for solving variational inequalities. Comput Math Appl 1981;7:43–58. [25] Campolongo F, Saltelli A. Sensitivity analysis of an environmental model: an application of different analysis methods. Reliability Eng Syst Saf 1997;57:49–69. [26] Salvador R, Piñol J, Tarantola S, Pla E. Global sensitivity analysis and scale effects of a fire propagation model used over Mediterranean shrublands. Ecol Model 2001;136:175–89. [27] Hamby DM. A review of techniques for parameter sensitivity analysis of environmental models. Environ Monit Assess 1994;32:135–54. [28] Saltelli A, Tarantola S, Chank K. A quantitative, model independent method for global sensitivity analysis of model output. Technometrics 1999;41(1):39–56. [29] Chapman B, Jost G, Van Der Pas R. Using OpenMP: portable shared memory parallel programming. Cambridge: MIT Press; 2007. [30] Buckingham E. On physically similar systems; Illustrations of the use of dimensional equations. Phys Rev 1914;4(4):345–76. [31] Sobol IM. Sensitivity estimates for non-linear mathematical models. Math Model Comput Exp 1993;1(4):407–14. [32] 2009. SIMLAB Version 2.2 Simulation environment for uncertainty and sensitivity analysis, developed by the Joint Research Centre of the European Commission. [33] Chetehouna K, Séro-Guillaume O, Sochet I, Degiovanni A. On the experimental determination of flame front positions and of propagation parameters for a fire. Int J Therm Sci 2008;47:1148–57. [34] Mokhtar SB, Hanif DS, Shetty CM. Nonlinear programming: theory and algorithms. 3rd ed. Wiley; 2006. [35] Michael LR, Torczon V. Pattern search algorithms for bound constrained minimization. SIAM J Optim 1999;9(4):1082–99.