Citation |

- Permanent Link:
- http://digital.auraria.edu/AA00003432/00001
## Material Information- Title:
- Reaction diffusion equations and numerical wildland fire models
- Creator:
- Kim, Minjeong
- Publication Date:
- 2011
- Language:
- English
- Physical Description:
- xvii, 139 leaves : ; 28 cm
## Subjects- Subjects / Keywords:
- Reaction-diffusion equations ( lcsh )
Wildfires -- Mathematical models ( lcsh ) Differential equations, Partial ( lcsh ) Differential equations, Partial ( fast ) Reaction-diffusion equations ( fast ) Wildfires -- Mathematical models ( fast ) - Genre:
- bibliography ( marcgt )
theses ( marcgt ) non-fiction ( marcgt )
## Notes- Bibliography:
- Includes bibliographical references (leaves 131-139).
- General Note:
- Department of Mathematical and Statistical Sciences
- Statement of Responsibility:
- by Minjeong Kim.
## Record Information- Source Institution:
- |University of Colorado Denver
- Holding Location:
- |Auraria Library
- Rights Management:
- All applicable rights reserved by the source institution and holding location.
- Resource Identifier:
- 761029679 ( OCLC )
ocn761029679 - Classification:
- LD1193.L622 2011d K55 ( lcc )
## Auraria Membership |

Full Text |

REACTION DIFFUSION EQUATIONS AND NUMERICAL WILDLAND FIRE
MODELS by Minjeong Kim M.S., Kyungpook National University, South Korea, 2001 A thesis submitted to the University of Colorado Denver in partial fulfillment of the requirements for the degree of Doctor of Philosophy Applied Mathematics 2011 This thesis for the Doctor of Philosophy degree by Minjeong Kim has been approved by Jan Mandel Loren Cobb Jiangguo (James) Liu Weldon Lodwick Date Kim, Minjeong (Ph.D., Applied Mathematics) Reaction Diffusion Equations and Numerical Wildland Fire Models Thesis directed by Professor Jan Mandel ABSTRACT In this thesis, we formulate a mathematical model based upon partial differential equations that describes wildland fires. We first investigate the reaction rate and then the PDE-based wildland fire model is discussed. The traveling wave solutions of the wildland fire model and measured temperature profile are presented. We also discuss the identification of coefficients, and a dimensionless form of wildland fire model. Standard numerical schemes are used for the numerical simulations, and various numerical tests are performed in order to examine the properties of the numerical model and numerical discretization. An efficient implementation of a finite differ- ence scheme in 2D is described and discussed. Furthermore, we present a fireline propagation model by using the level set method, and observed problems in the level set method are discussed. Lastly, we analyze a reaction-diffusion model in ID in phase space. We present some concepts of dynamic systems and prove the existence of a traveling wave in phase space under some assumptions. We obtain rigorous inequalities which show the existence of a traveling wave solution in phase space, and we finally present numerical solutions of the traveling wave in phase space. This abstract accurately represents the content of the candidates thesis. I recommend its publication. Signed Jan Mandel ACKNOWLEDGMENT I would like to give my special thank to Professor Jan Mandel as my advisor, and as a mentor. Without his support, guidance and endless patience, my dissertation would have been impossible. I am also indebted to Professor Lynn Bennethum and Professor Loren Cobb for their advices. I thank to Professor Sangdong Kim for motivating me to study abroad. I would also like thank the faculty and my fellow students at the Department of Mathematical and Statistical Sciences for providing great working conditions. In particular, thank to Dr. Bedrich Sousedik and Dr. Christopher Harder for studying and discussions. Lastly, I thank to my family, and in particular, my husband and my parents. Contents 1 Introduction 2 2 PDE-based model 7 2.1 Introduction.............................................. 7 2.2 Formulation of the model................................. 10 2.3 Derivation of the model ................................. 13 2.3.1 Heat and fuel supply balance equations............. 13 2.3.2 Reaction rate...................................... 17 2.3.3 Heat balance with full fuel supply..................20 2.4 Identification of coefficients .....................23 2.4.1 Reaction rate coefficients..........................23 2.4.2 Cooling coefficient ................................26 2.4.3 Scales and non-dimensional coefficients.............29 2.4.4 Related work....................................... 36 2.5 Calibration of coefficients in one dimension..............37 vi 3 Numerical methods for the PDE model 41 3.1 Efficient numerical implementation............................41 3.2 Computational tests ..........................................48 3.2.1 Line ignition with oscillatory wind ...................50 3.2.2 Spot fire without wind.................................50 3.2.3 Effect of the wind speed...............................52 3.2.4 Some observations .................................... 56 3.3 Numerical convergence tests.................................. 58 3.4 Numerical difficulties........................................61 4 The fireline propagation model 68 4.1 Level set based wildland fire model.......................... 70 4.1.1 Some basics about the level set method................ 71 4.2 Numerical schemes for fireline propagation model............. 74 4.2.1 Numerical schemes for the Level Set equation...........74 4.2.2 Observations and future work...........................77 5 Phase-space analysis of the traveling wave 83 5.1 Traveling wave in phase space.................................84 5.1.1 Parameterized system of ODEs in phase space............85 5.2 Properties of dynamical systems...............................93 5.3 Traveling wave in the fire model as a dynamical system .... 100 5.3.1 Invariant manifolds...................................100 5.3.2 Nullclines............................................103 vii 5.4 Inequalities for the speed of the traveling wave..........Ill 5.5 Existence of a traveling wave speed.......................118 5.6 Numerical solution in phase space.........................123 6 Conclusion 129 viii List of Figures 2.1 Sample reaction heat balance function f(T) from equation (2.21)............................................................21 2.2 Reaction heat balance potential U.................................22 2.3 Solution with the Arrhenius reaction rate. Due to nonzero reaction rate at ambient temperature, fuel starts disappearing and thus a propagating combustion wave does not develop. The coefficients are k = 2.1360 x 10_1m2s_1 A*3, A = 1.8793 x 102Ks~\ B = 5.5849 x 102A, C = 4.8372 x lO^A-1, and Cs = 1-6250 x lO^s"1. Tc = 1200A and T{ = 670A...............24 2.4 Solution with Arrhenius reaction rate modified by temperature offset Ta to force zero reaction rate at ambient temperature. A propagating combustion wave develops. The coefficients are k = 2.1360 x 10-1m2s-1 A-3, A = 1.8793 x 102As-\ B = 5.5849 xlO2 A, C = 4.8372x10- A"1, and Cs = 1.6250 x lO-V1. Tc = 1200A and Tj = 670A.............................. 27 IX 2.5 Temperature profile of a traveling wave. The wave moved about 398m from the initial position in 2300s...............28 2.6 The fuel fraction remaining after the combustion wave as a function of dimensionless parameters A and (3...............32 2.7 Ratio of the leading and the trailing edge of combustion wave with respect to A and (3 from (2.31). The width is measured as the length of the part of the wave above temperature 500A7 33 2.8 Width of the trailing slope as a function of A and /?. The rela- tion between the wave speed and (diffusion coefficient)0,5, i.e. \fk. We see approximately linear relations.........................35 2.9 Periodic patterns of the solutions by varied (3. (a) (3 = 0.4172, (b) (3 = 0.4636, (c) (3 = 0.4945, (d) /3 = 0.5718. Each graph shows the temperature profile by timeframe...........38 2.10 Time-temperature profile (dotted line) measured in a grass wildland fire at a fixed sensor location, digitized from Kremens et al. [2003], and a computed profile (solid line) from simulation. 39 3.1 Five-point stencil........................................... 44 3.2 Pentadiagonal matrix Jtt by using standard five-point finite difference discretization.....................................45 3.3 Diagonal matrix Jts, Jst, and Jss............................46 x 3.4 Column storage of the sparse Jacobian matrix. Five columns, labelled as i m, i 1, i, i + 1, and i + m, are from the pentadiagonal matrix Jtt The last three columns store the diagonal matrices Jts> Jst, and Jss........................47 3.5 North to south straight line ignition and a sinusoidal wind with speed 3m/s after 400s. The dark areas are burnt regions. 49 3.6 Spot fire without wind simulations by Section 2.2 equations (2.2-2.3). The mesh size is 4m and time step is Is after 600s. (a) is the numerical fire evolution and (b) is a contour of (a). (c) is fuel consumption in (a).............................51 3.7 Effect of wind speed 0.005m/s fire simulation of Section 2.2 equations (2.2-2.3) after 800s. Time step is Is and mesh size is 4m. With a small wind speed the fire also moves more slowly and it is closer to circular shape, (b) is a contour plots of (a). (c) is a fuel consumption in (a)...........................53 3.8 Effect of wind speed 0.05m/s fire simulation of Section 2.2 equations (2.2-2.3) after 400s. The numerical solution moves from southwest to northeast. Time step is Is and mesh size is 8m. (b) is a contour plots of (a), (c) is a fuel consumption in (a)........................................................ 54 xi 3.9 Effect of wind speed 3m/s fire simulation of Section 2.2 equa- tions (2.2-2.3). Fire moves from southwest to northeast. Time step is Is and mesh size is 16m. (b) is contour plots of (a), and (c) is a fuel consumption of (a). Note that the sharp front edges in (a).................................................. 55 3.10 Wind effect of fuel consumption, (a) Fuel consumption with domain size 256m x 256m and wind speed 3m/s after 20s. (b) Fuel consumption with wind speed 5m/s. The dark area means fuel consumption and light yellow is unburnt region. We see that the fire spreads faster with a higher wind speed. 57 3.11 Errors from Table 3.1 in the logarithmic scale(left), and the convergence norm ratios(right). Errors from Table 3.2 in the logarithmic scale(left), and the convergence norm ratio(right). 60 3.12 Errors from Table 3.3 in the logarithmic scale(left), and the convergence norm ratios (right). Errors from Table 3.4 in the logarithmic scale(left), and the convergence norm ratio(right). 62 3.13 Fire simulation with large mesh size of 80m and wind speed of 0.005m/s. The time step is 2s, and ignition area is 400m by 320m. (a) is the initial plot of (b).......................65 3.14 Fire simulation with mesh size of 2.5m after 600s. All other conditions are the same in Figure 3.13........................66 xii 3.15 Fire simulation with large mesh size of 80m and wind speed of 1 m/s. While the case of small wind speed (3.13), a traveling fire develops............................................................66 4.1 the normal to the curve(4.6)...................................... 72 4.2 A circular fire propagation with constant wind. The contour is a fireline...................................................... 77 4.3 (a) A trace of circle movement along wind direction. Initial value circle, (b) The circle moved to the south-east. As time is passing, the circle keep shrinking. The re-initialization avoids the shrinking problem in (c)..................................... 78 4.4 A trace of motion in a rotation velocity field. Re-initialization method is added in (a) and (b)................................... 79 4.5 (a) Points on the zero level set on the mesh grid lines are stored as a set A and points on the zero level set function, which is not on the mesh grid lines, are stored as a set B. (b) Closest points to the grid points are stored as a set C..................80 4.6 Motion of a circle in a rotational velocity field. The Narrow Band scheme with re-initialization is added in (a) and (b). There still exists some mass loss while the circle is rotating and eventually the circle does not match at the end.............. 82 5.1 Intersection of the solution of equations (5.35) (5.36) with the plane V = 0 at (zq, Cq) xiii 98 5.2 The unstable, stable, and center subspaces, Eu, Es, and Ec, and manifolds, WU(Xo), VT5(Xo), and WC(Xo) respectively. 102 5.3 Nullclines Vz = 0 of (5.32) with different wave speed for 5 = 0, 5 = 0.5, and 5=1 and the traveling wave speed c = 0.45 in (a), c = 0.4 in (b), c = 0.35 in (c). The parameters are A = 0.25 and (3 = 0.2...........................................104 5.4 Nullclines V = Â£(Se~7 XT) and V = 0, for a fixed 5 with direction fields. Scaled to unit length in the (V,T) projection. 105 5.5 Direction field of equation (5.35) in the plane 5=1 with the nullcline Vz = 0 and tangent line to the nullcline at the origin in (a), in the plane 5 = 0.8 in (b), and in the plane 5 = 0.5 in (c) also shown. In (a) and (b), Tig is the ignition temperature. The parameters are c = 0.4, A = 0.25, and (3 = 0.2..............106 5.6 The solutions of the systems (5.35)-(5.36), projected on the (T, V) plane. The solution is not a traveling solution because the solution does not crosses the nullcline Vz 0..............107 5.7 The solution of equations (5.35)-(5.36). A traveling wave so- lution is shown. The view is rotated around the V-axis. The parameters are c = 0.7992, A = 0.06, and (3 0.15671... 109 xiv 5.8 Nullcline Vz = 0 and the solution of equation (5.35). The view is rotated around the V-axis. The maximum value of the nullcline is negative and the solution keeps going upward from the critical point. Thus, no traveling wave can exist. The 5.9 Solution of equations (5.35)-(5.36) in the plane 5=1 and S = 0.9. We see that A) in the region where V > ^(e_1 A)T. The solution of equations (5.35)-(5.36) keeps going upward, and this solution is not the traveling wave solution. The parameters are c = 0.4, A = 0.25, and (3 0.4. 113 5.10 The solution of equation (5.35)-(5.36) in 3D. In (a), c is too large, and the solution with the plane T = 0. In (b), c is too small and the solution with the plane V = 0. In (c), the end point of the solution goes to the limit point (5.37), and this solution is the traveling wave solution. The parameters parameters are c = 0.075, A = e *, and (3 = 0.15671 110 A = 0.06 and (3 = 0.15671 114 5.11 Function e t AT in Theorem 8, which is denoted by 0(T) here. The parameters are c = 0.6 and (3 = 0.4................ 115 5.12 (a) p with varied A in equation (5.50). (b) with varied A. The speed of the traveling wave c = 0.07992 and c < (c) The traveling wave solution of equations (5.35)-(5.36) with the plane (T, V) . 116 xv 5.13 (a) hfi- with varied A and the parameter c = 0.6. We see y 1 that the parameter c > (b) The solution of equations (5.35) -(5.36) with the plane (T, V) and the parameter c = 0.6. The solution of equations (5.35)-(5.36) keeps increasing, and clearly it is not the traveling wave solution of the equation (5.35) -(5.36). The parameters are 0 = 0.4 and A = 0.254. 5.14 Function <&(c) is continuous............................. 117 121 xvi List of Tables 3.1 Convergence ratio for decreasing time step, wind speed 2m/s, the mesh size h = 2m....................................59 3.2 Convergence rate for decreasing time step, wind speed 4m/s, the mesh size h = 4m....................................61 3.3 Numerical convergence behavior for both the time step dt and the mesh size h, wind speed 0.5m/s......................61 3.4 Numerical convergence behavior for both the time step dt and the mesh size h. Wind speed is 8m/s.....................63 xvn Chapter 1 Introduction Reaction-diffusion equations arise in a large area of mathematical models involving chemical reactions. As a practical application, there is a wide interest in such equations, with numerous recent studies being performed on, in particular, biological models such as the epidemic model and bacterial colonies. We present a highly simplified wildland fire model as an application of reaction diffusion equations. This model has two physical variables, the fuel temperature and the fuel mass. The speed of propagation of the fire is a property of the solution of a system of differential equations, and it is thus determined by the coefficients of the equations. One reason for considering a PDE-based model is that even simple reaction- diffusion equations are capable of complex nonlinear, unsteady behavior such as pulsation and bifurcation. We note that the physical behavior can be 2 achieved by very simple models that can reproduce qualitative behavior of some features of fire. On the other hand, it is not easy to determine the coefficients of the PDEs. In this thesis, properties of the combustion wave are used to calibrate the parameters of the model. Exact solutions to the PDE wildland fire model are not known. Thus, computational simulation is important. In addition these PDE models can- not exactly replicate the physical phenomena. Furthermore, numerical so- lutions can exhibit errors related to the implementation itself. Therefore, numerical tests and observations on the physical correctness of numerical so- lutions are necessary. Efficient implementations can dramatically reduce the costs of the simulation. Although the reaction-diffusion-advection equations with enriched finite element methods are designed and error estimates are given in Asensio et al. [2004], Codina [1998], Franca et al. [2005], and Franca et al. [2006], these works do not account for nonlinear effects at the leading edge. Error esti- mates are provided in Asensio and Ferragut [2000] and Ferragut and Asensio [2002] in the context of mixed finite elements, which is applied to a single species combustion equation and more general reaction-diffusion problems, but these works do not consider a fuel-balance equation. In this thesis, we present a nonlinear reaction-diffusion equation with fuel consumption, and we obtain traveling wave solutions. Although it is known that coupling with the atmosphere is essential for 3 wildland fire behavior, such interaction is not considered in this thesis. We instead limit our attention to investigating the wildland fire models by them- selves and finding numerical solutions. Because a common feature of the solutions of reaction-diffusion equations is the development of a sharp wave, with the solution being almost constant elsewhere, interface tracking techniques such as level set methods (Sethian [1999] and Sussman et al. [1994]) are relevant here as well. This method is a basis of the implementation in WRF-FIRE (Mandel et al. [2009], Mandel et al. [2011]), and this empirical model developed by NCAR (Clark et al. [2004], Coen [2005a]). In this empirical model, the fire propagation speed is normal to the fireline and is a function of wind and terrain slope, and an exponential decay of fuel from the time of ignition is assumed. Thus, the fire propagation speed can be input directly from known fire behavior for various fuels. We observe that the numerical solution of PDEs is a traveling wave so- lution. This observation gives us a motivation to study the traveling wave solution of the equations. Thus, we are interested in the traveling wave solu- tions of reaction-diffusion equations. The traveling wave solution has a form in connection with a time and a speed of wave, and this form of solution allows us to have ODEs instead of PDEs. In general, ODEs are easy to analyze, and it often allows us to find an explicit solution of the equation. The parameters of the equations determine a behavior of the solution in phase space. The exact solution of reaction-diffusion equations should 4 represent the behavior of the model or physical properties. The analysis of the direction field in phase space clearly indicates the properties of the solution in phase space which depend on parameters. These dependencies provide rigorous properties of traveling wave solutions in phase space. This thesis is organized as follows. In Chapter 2, we formulate a simple PDE-based wildland fire model and identify the coefficients instead of taking these from physical and chemical properties of fuel. In Chapter 3, we first introduce a familiar discretization of unsteady state reaction-diffusion-convection equations, and then provide an efficient imple- mentation. Various computational tests and numerical convergence tests will be performed. Finally, we will discuss observed numerical difficulties. In Chapter 4, we introduce the level set method and apply it to the wild- land fire problem. An observed numerical drawback, the mass loss problem, is introduced from some numerical experiments. We then apply simple res- olutions to cure this difficulty, and then it turns out that preventing this difficulty is a interesting area in studying the level set method. In Chapter 5, we turn our attention to the phase-space method for the reaction-diffusion equations. The phase-space method is often used to reduce second order equations to first order equations. Then the first order equations become autonomous dynamic systems. We introduce the basic concepts of dynamics and investigate a phase-space direction field, which is 3D here. We prove the existence of a traveling wave under suitable assumptions. The authors personal contributions in this work are as follows. Chapter 5 2: Derivation of the simple PDE-based model (2.2)-(2.3) in (2.9)-(2.12) and identification of coefficients. Analysis of dimensionless parameters. Chapter 3: Matlab and Fortran implementation of equations (2.2)-(2.3). Compu- tational tests of numerical convergence and observations of numerical diffi- culties. Furthermore, the efficient implementation of equations (2.2)-(2.3) allowed for faster data assimilation in the DDDAS(Dynamic Data Driven Application Systems) project. Chapter 4: Survey of level set method and the investigation of problems of level set method in experiments. Contri- butions to the implementation of the method and fuel burning in the fire code in WRF-Fire. Chapter 5: Rigorous results in equations (5.35)-(5.37). Phase-space analysis of the traveling wave using invariant manifolds. Joint work with advisor: Proofs of Lemmas and Theorems showing that a traveling wave does not exist under certain conditions, and thus proving inequalities for the parameters of the traveling wave. Proof of the existence of traveling wave solutions in phase space under some assumptions. 6 Chapter 2 PDE-based model 2.1 Introduction Modeling wildland fire requires us to understand the properties of fire be- haviors. This modeling contains complexities of fuel characteristics and heat transfer mechanism with the effect of wind. We formulate a highly simple mathematical wildland fire model. It is a plausible mathematical model of wildland fire which helps to understand the physical meanings of a compli- cated real wildland fire. Various aspects of special cases of the reaction-diffusion equations (2.2)- (2.3) in Section 2.2 below have been studied in a number of papers. A formal expansion in an Arrhenius reaction model to get the wave speed and a pre- diction of whether a small fire will or will not spread was used by Weber [1991]. The ignition wave speed and the extinction wave speed were com- 7 puted numerically by Mercer and Weber [1995]. The speed of combustion waves were analyzed by asymptotic expansion and the stability of the travel- ing wave front was obtained by linearized system of ODEs and its eigenvalues by Gubernov et al. [2003]. The speed of the traveling wave in a simplified model with the reaction started by ignition at a given temperature was derived by Campos et al. [2004] under the assumption that a temperature-independent reaction rate is proportional to the fuel remaining, similar to the model by Balbi et al. [1999]. This temperature-independent reaction rate is often used in a PDE- based model, and the reaction rate of the model has a proportionality to the fuel remaining. The speed of a combustion traveling wave fronts was investigated by using relaxation and shooting methods by Gubernov et al. [2003] and Gubernov et al. [2004] from the solid to gaseous fuel (i.e., also with fuel diffusion instead of just temperature diffusion) for different values of parameter, j3 (the ratio of activation to heat release). They found that the large /? and gaseous fuel give a good correspondence between analytical and numerical results. A nonlinear eigenvalue problem for a traveling wave in a different combus- tion problem, with fuel reaction, was derived by Norbury and Stuart [1988a] and Norbury and Stuart [1988b], The equation was solved numerically by the shooting method, and the existence and stability of the traveling wave solution is studied. For an analytical study of the evolution of waves to a traveling waveform, 8 see Sherratt [1998], and for a numerical study, see Gazdag and Canosa [1974] and Zhao and Wei [2003]. Proofs of the existence of a solution and attractors for reaction-diffusion equations are given by [Robinson, 2001, Ch. 8 and 11], although traveling waves are not mentioned. The existence of a solution for a finite time is proved by Asensio and Ferragut [2002] for the reaction-diffusion equation of (2.2) without convection, but the nonlinear diffusion term was replaced by V (kT3VT). Again, traveling waves and fuel depletion were not considered. Equation (2.2) without fuel depletion (i.e., with constant S) and without wind (i.e., v = 0) is a special case of the nonlinear reaction-diffusion equation, V (VT) + / (T). (2.1) Reaction-diffusion equations of the form of equation (2.1) are known to pos- sess traveling wave solutions, which switch values close to stationary states given by / (T) = 0 (Gilding and Kersner [2004] and Infeld and Rowlands [2000]). The simplest model problem is Fishers equation with f (T) = T(l T), for which the existence of a traveling wave solution and a for- mula for its speed were proven by Kolmogorov et al. [1937]. Two reaction models (solid and pyrolysis gas) are considered by Seron et al. [2005], who argue that the modeling of pyrolysis by a separate reaction is essential for capturing realistic fire behavior. For more complicated theo- retical models of this type that include other components, e.g., water vapor, 9 see Grishin [1996]. For a similar system, but with fuel diffusion, the existence and the speed of traveling waves were obtained by asymptotic methods already in classi- cal works, summarized by the monograph Zeldovich et al. [1985]. For the system (2.2)-(2.3), an approximate combustion wave speed by asymptotic methods was obtained by Weber et al. [1997]. However, the heat loss term, C (T Ta), was not included in the approximation. No heat loss case is equivalent to our setting C = 0 in equations (2.2)-(2.4). The proposed model in this chapter is based on a reaction-diffusion- convection equation (2.2)-(2.4) in Section 2.2. A more detailed derivation of the model from physical principles such as the conservation of mass and energy is contained by Section 2.3. In Section 2.4, we identify coefficients, using a dimensionless form of reaction-diffusion-convection equation (2.29)-(2.30). Non-dimensional pa- rameters presented in this section are used to identify coefficients. This chapter is partly based on Mandel et al. [2008]. 2.2 Formulation of the model We consider a model of a fire in a layer just above the ground. First, we define the following terms: T (K) is the temperature of the fire layer (K is degrees Kelvin), 10 S [0,1] is the fuel supply mass fraction (the relative amount of fuel re- maining: the mass of remaining fuel divided by the total mass of fuel), k (m2s 1) is the thermal diffusivity, A (Ks 1) is the temperature rise per second at the maximum burning rate with full initial fuel load and no cooling present, B (K) is the proportionality coefficient in the modified Arrhenius law below, C (K-1) is the scaled coefficient of the heat transfer to the environment, Cs (s-1) is the fuel relative disappearance rate, Ta (K) is the ambient temperature, and v (ms-1) is the wind speed given by atmospheric data or model. The model equations consist of the conservation of energy and balance of fuel: ^ = V (kVT) v VT + A (5e-B/(T-T) C (T Ta)), (2.2) (2.3) with the initial values S (t0) = 1 and T (t0) = Tinit. (2.4) 11 The diffusion term V (&VT) models short-range heat transfer by radia- tion in a semi-permeable medium, v VT models heat advected by the wind, Se-B/(T~Ta) jg ra^e fuej jg consumed due to burning, and the fuel deple- tion gains heat when the fire is ignited, i.e., T > Ta. AC (T Ta) models the convective heat lost to the atmosphere. The reaction rate e~B^T~Ta'> is obtained by modifying the reaction rate e~B!T from the Arrhenius law by an offset to force zero reaction at ambient temperature, with the resulting reaction rate smoothly depending on temperature. The reaction rate from the Arrhenius law will be discussed in Section 2.3.2. It is known that systems of a form similar to (2.2 2.3) admit traveling wave solutions. An approximation to the temperature reaction equation gives the size of the reaction zone and the slope of the temperature curve. The temperature in the traveling wave has a sharp leading edge, followed by an exponentially decaying cool-down trailing edge. This was observed numerically but we were not able to find a rigorous proof in the literature for exactly this (2.2 2.3) case. Models that do not guarantee zero combustion at ambient temperature suffer from the cold boundary difficulty: by the time a combustion wave gets to a given location, the fuel at that location is depleted by the ongoing reaction at ambient temperature. So, no perpetual traveling combustion waves can exist, and there are only pseudo-waves that travel only for a finite time and then vanish (Berestycki et al. [1991], Mercer et al. [1996], and Zeldovich et al. [1985]). 12 2.3 Derivation of the model We consider fire in a ground layer of some unspecified finite small thickness h. The fire layer consists of the fuel and air just above the fuel. All mod- eled quantities are treated as two dimensional, homogenized in the vertical direction over the ground layer. We will not attempt to substitute the coef- ficients from material properties because of the degree of simplification and the uncertainty present in the homogenization. Instead, physical laws will be used to derive the form of the equations and the coefficients will be identified later from the dynamic behavior of the solution. We first derive the system of PDEs based on conservation of energy and fuel reaction rate in Section 2.3.1 and then discuss the choice of the reaction term in Section 2.3.2 2.3.1 Heat and fuel supply balance equations In this section, we write the equations in a form suitable for identification of the coefficients, which will be formally done in Section 2.4 The chemical reactions are a heat source and the chemical reaction of burning is generating the heat. We model the burning as a reaction where the reaction rate only depends on the temperature. Then the reaction rate with a coefficient of proportionality, Cs (1 /s), is Csr(T), and r is a dimensionless reaction rate. Let Fc > 0 (kgm~2) be the concentration of fuel remaining. Then the fuel consumption rate is proportional to the reaction rate where the amount 13 of fuel available, -^ = -F'Csr(T) We introduce the mass fraction of fuel by (2.5) S Fc Fo where F0 is the initial fuel quantity. The heat generated per unit surface area is then proportional to the fuel lost, qg = AiFcCsr (T) (Wm~2), (2.6) where Ay (J kg~x) is the heat released per unit mass of fuel. Heat transfer is due to radiation and convection to the atmosphere. The short-range heat transfer due to radiation and turbulence is modeled by dif- fusion. The two dimensional heat flux vector then is qr = AqVT (Wm"1). (2.7) Equation (2.7) models the heat transfer from the area of high temperature to the area of low temperature through the unit length, and equation (2.7) guarantees this heat transfer. The constant Aq (WK~l) will be identified later. Heat per unit area lost due to natural convection to the atmosphere and 14 to buoyancy is given by Newtons law of cooling, qc = Ca(T Ta) (Wm~2) , (2.8) where Ta is the ambient temperature (K) and Ca (Wm~2K~x) is the heat transfer coefficient. In this model, it is assumed that the convective heat transfer is dominant, and so the effect of radiation into the atmosphere is included in (2.8). The material time derivative, which is a derivative taken along a path following the flow with velocity field, v, of the temperature is given by DT ~Dt dT = Ti+'"vt- (2.9) Here dT/dt is the Eulerian (or spatial) time derivative of temperature, v is the (homogenized) velocity of the air, p is the homogenized surface density of the fire layer (kg m-3), and cp is the homogenized specific heat of the fire layer (J kg~1K~l). Again, none of the coefficients h, p, or Cp can be assumed to be actually known. The units of the product hpcp are JK~lm~2. The velocity vector v is obtained from the state of the atmosphere as data, or in future work by coupling with an atmospheric model. From the divergence theorem, we now obtain the conservation of energy in the fire layer as DT hpCp~Dt = V Qr + Q9 ~ Qc ^21^ 15 The goal is to obtain a system of equations with a minimal number of coefficients and in as simple form as possible. In addition, we wish to relate the behavior of the solution to the non-dimensional parameters in Section 2.4.3 below, rather than to material and physical properties of the medium, which are in general unknown. Substituting in the appropriate values for the heat sources and fluxes, (2.7), (2.6), and (2.8), into (2.9) and (2.10), and some simple algebra, we obtain the energy balance and fuel reaction rate equations = V-(fcVT)-v VT + A(Sr(T)-C0(T-Ta)), (2.11) ^ = -CsSr(T), (2.12) with k = ki/(hpcp), A = AiCs/{hcpp), and C0-Ca/Ai. Alternatively, we could have started with the disappearance of fuel on the left hand side of the heat balance equation (2.10) (fuel that has burned does not need to be heated), which would lead to an equation of the form (1 + CjS) ^ = V (kVT) v VT + A (Sr (T) -C0(T- Ta)) at instead of (2.11). We have chosen not to do so since our goal is to work with the simplest possible model, whose coefficients can be identified. The fuel dis- 16 appearance in the heat equation affects the temperature profile significantly only in the reaction zone, which is the highest part of the temperature curve. Before the ignition, there is a full fuel load, 5=1, and after a fairly short reaction time, well before most of the cooling takes place, the remaining fuel settles to some residual value which then remains constant. The effect of the decreased heat capacity of the remaining fuel is then absorbed into the cooling term ACq. 2.3.2 Reaction rate In this section we discuss the reaction rate with an Arrhenius type of equa- tion, and we derive a modified equation by introducing dimensionless tem- perature and parameter. The Arrhenius equation from physical chemistry is given by r{T) = Ae(~-tr\ (2.13) where r is a rate constant of reaction and T is a temperature, and the reaction rate falls exponentially with decreasing in the temperature. This equation gives the reaction rate as a function of temperature, where A is the pre- exponential factor, Ra is the gas constant, Ea is an activation energy, and Ea/Ra has units K. Equation (2.13) is valid only for a single reaction and lower temperature regime, but it is used here as an approximation form. The activation energy Ea is mostly obtained empirically. One consequence of (2.13) is that the reaction has a nonzero rate at some temperature above 17 absolute zero. As already explained in Section 2.2, the cold boundary problem is caused by the fact that the time scale of burning is much smaller than the oxidation rate at ambient temperature (Williams [1985]). A dimensionless parameter e* of equation (2.17) below is used in Mercer et al. [1996], Asensio and Ferragut [2002], Mercer and Weber [1995], and the numerically investigated results in the dependency of dimensionless parameters with respect to a wave speed can be found in Mercer et al. [1996] and Mercer and Weber [1995]. However, none of these papers explains why they used equation (2.16) as the reaction rate, and they did not mention why the dimensionless parameter e* needs to be introduced. We summarize the assumptions and a process of deriving the dimension- less parameter e* (Zeldovich et al. [1985]). The idea is that the parameter e* is approximating the reaction rate in a narrow range of temperature. The narrow range of temperature, AT = T T, where the temperature T* is taken from a neighborhood where the reaction occurs. This temperature T* can be Ta in a study of spontaneous ignition, or it can be the combus- tion temperature in a study of propagating combustion. In the problem, T, is determined during problem solving since T* is the unknown. Here, the temperature T is the ambient temperature Ta. By the method of expanding the exponent, which is originally suggested 18 in Williams [1985], we introduce the dimensionless temperature difference u, u = Eg RaT? (T T*). (2.14) From the identity we obtain 1 +at 1 T. = 1 - AT T. Ea 1 Ea RgT ~ T(T + AT) TaT, 1 + ^ ~ RaTt ( ) (2.15) Equation (2.15) allows us to rewrite the reaction rate as ___ Eg __ Eg U c =6 KaT* gHTeTtTj (2.16) where u is given by (2.14) and the parameter e. RaTt Eg ' (2.17) We simply rewrite the Arrhenius reaction rate (2.16) as r(T) = e-B'T, (2.18) where the coefficient B has units K. This equation is valid only for gas fuel premixed with a sufficient supply of oxygen. This approximation ignores fuel 19 surface effects but it is widely used nonetheless. We establish the reaction rate by neglecting the initial temperature, i.e., by assuming that the reaction rate is equal to zero at the ambient temperature and modify (2.18) so that no oxidation occurs below some fixed temperature, To (Section 2.4), and take instead Note that the fuel consumption rate is a smooth function of T, which is favorable for a numerical solution, unlike the rate in Asensio and Ferragut [2002], where a cutoff function was used. 2.3.3 Heat balance with full fuel supply Consider first the hypothetical case in which T is constant in space, so that only heat due to the reaction (burning) and natural convection contribute non-zero terms in the heat equation, (2.2), and essentially the full initial fuel supply Fo is present at all times (the rate of fuel consumption is negligible, Cs ~ 0, so S ~ 1). Then, (2.11) reduces to Constant values of temperature which are solutions to (2.20) are called equi- librium points, and at these points the heat produced by the reaction equals (2.19) = A (e-B/(T-To) -C(T- Ta)). (2.20) 20 T =300 T =1200 T=670 a c i Figure 2.1: Sample reaction heat balance function f(T) from equation (2.21). the heat lost to the environment, f(T) = e-BI(T-To) -C{T-Ta) = 0. (2.21) Equation (2.21) has at most three roots (Frank-Kamenetskii [1955]), see for example, Figure 2.1. The first zero, denoted as Tp, called the lower temperature regime by Frank-Kamenetskii [1955], is a stable equilibrium temperature. If the tem- perature goes below this temperature then the heat generated from the re- action dominates and the temperature rises, since some reaction is present even at ambient. If the temperature goes above this temperature then cool- ing dominates and the temperature decreases. The middle zero, 7), is an 21 10 T =300 T =1200 T.=670 a c i Figure 2.2: Reaction heat balance potential U unstable equilibrium point. If the temperature goes above Ta and below then cooling dominates, and the temperature decreases. Above 7), the heat due to chemical reactions dominate and the temperature increases. We refer to Ti as the auto-ignition temperature, the temperature above which the re- action is self-sustaining (Quintiere [1998]). The stable equilibrium Tc is the maximum stable combustion temperature, assuming replenishing of the sup- ply of fuel and oxygen. The temperature Tc is called the high temperature regime by Frank-Kamenetskii [1955]. The stability properties of the equi- librium points are also clear from the graph of the potential U(T), defined by U'(T) = f(T). The stable equilibrium points are local minima of the potential, while the auto-ignition temperature is a local maximum and thus an unstable equilibrium (Figure 2.2). 22 2.4 Identification of coefficients We wish to calibrate the coefficients from physically observable quantities of the fire rather than physical material properties. It is not simple to obtain a reasonable behavior of the solution from substituting physical coefficients values into the equations. Further, as explained in Section 2.3, because of a number of simplifying assumptions employed and because of the homoge- nization of coefficients over a fire layer of unspecified thickness, it is not quite clear what the material properties should be anyway. We first use basic reaction dynamics and a reduced model to find rough approximate values of the coefficients that produce a reasonable solution. After that, we transform the equations to a non-dimensional form, which allows us to separate the coefficients into those that determine the qualitative behavior of the solution and those that determine the scales. We use the approximate coefficients obtained from the reduced models as initial values for the identification of the coefficients by the nondimensionalization method to match observed temperature profiles. 2.4.1 Reaction rate coefficients While the coefficients B and C/F0 in (2.21) are generally unknown, they can be found from the equilibrium temperature points. Suppose that two roots Ti and Tc of f(T) from (2.21) are given such that T0 23 Fuel supply mass fraction Figure 2.3: Solution with the Arrhenius reaction rate. Due to nonzero reaction rate at ambient temperature, fuel starts disappearing and thus a propagating combustion wave does not develop. The coefficients are k = 2.1360 x 10-1m2s-1 A-3, A = 1.8793 x 102Ks~\ B = 5.5849 x 102A, C = 4.8372 x 10-5K~\ and Cs 1.6250 x 10_1s_1. Tc = 1200K and Ti = 670K. 24 / (Ti) = / (Tc) = 0, thus e~B/Ti (7 (T< Ta) = 0, e-B/Tc-C(Tc-Ta) = 0. Eliminating C gives 0 = e~BtTc (T{ Ta) e-B/T (Tc-Ta), which is equivalent to Tj-Ta Tc-Ta _ e-B(l/Ti-l/Tc) = Q5 hence Tc-T0 Ti-To -B/(Ti-T0) and C = ------ Ti-Ta (2.22) When T{ is too close to T, the resulting energy balance equation f(T) = 0 has only two roots. This, however, does not occur for the values of Ta, Tt, and Tc of interest. First consider the solution of (2.11-2.12) and the reaction rate (2.19) with T0 = 0. Then the reaction is the Arrhenius rate known from chemistry and there is a nonzero reaction rate at T = Ta. This results in fuel loss everywhere, and, in our computational experiments, no traveling combustion 25 wave developed in Figure 2.3 since, after a relatively short time, there was not enough fuel to sustain combustion due to the cold boundary difficulty (see Section 2.2). Since the values of B and C are obtained from realistic Ti and Tc, the fuel disappearance happens quickly and the traveling wave is sustained (Figure 2.3), we force the reaction rate, r(T), to be zero at ambient temperature by choosing the offset T0 = Ta. Using the offset by Ta is essentially the same as assuming that the ambient temperature is absolute zero, as commonly done in combustion literature Weber et al. [1997]. In this case, a propagating combustion wave develops on Figures 2.4 and 2.5. Therefore, we use T0 = Ta. ft should be noted that traveling combustion waves, such as in Figure 2.5, are caused by the combined effect of reaction and diffusion; convection does not play a role in this section. The reaction heat diffuses forward on the leading edge, heating the fuel ahead of the wave, until the reaction ahead of the wave can sustain itself, thus causing the combustion to spread. On the trailing edge, the reaction drops off due to fuel depletion, after which temperature decays due to cooling. 2.4.2 Cooling coefficient The coefficients B and C in the modified reaction form (2.21) have been determined from reasonable values of Tc and Ti by (2.22). Now we want to determine the remaining coefficient, A, which can be done if we know the characteristic cooling time. Consider the trailing edge of a traveling combus- 26 reaction Figure 2.4: Solution with Arrhenius reaction rate modified by temperature offset Ta to force zero reaction rate at ambient temperature. A propagating combustion wave develops. The coefficients are k 2.1360 x 10_1m2s_1 K~3, A = 1.8793 x 102Ks~\ B = 5.5849 x 102A, C = 4.8372 x lO"5^-1, and Cs = 1.6250 x lO^s-1. Tc = 1200K and T* = 670A. 27 Traveling combustion waves Figure 2.5: Temperature profile of a traveling wave. The wave moved about 398m from the initial position in 2300s. tion wave, after all or most of the fuel has been depleted, temperature drops, and heat generated by the reaction and diffusion drop to an insignificant level. From that point on, the temperature approximately satisfies dT - = -AC(T-Ta). Thus, at the trailing edge, given T at some time to, we have T(t) = Ta + (T(t0) Ta) e~AC^ and we can define the characteristic cooling time, tc, (see similar discussions in Balbi et al. [1999]), to be the time which the fire layer takes to cool by a 28 factor of 1/e, i.e., T(t0 + tc) Ta = (T(t0) Ta), which gives that ACtc 1, or A (2.23) 2.4.3 Scales and non-dimensional coefficients We now write the model in terms of nondimensional variables, which control the qualitative behavior of the system (2.2-2.3). Again, we do not consider the wind here yet, and so = V (kVT) + A (^Se~T^ C (T Ta)j , (2.24) fjQ B = -SCse~T=K, T > Ta. (2.25) at Given time temperature scale Tj, space scale xx, and time scale t,x, con- sider the substitution T=1 x = t = -. (2.26) Tj xi 11 29 The equation (2.24)-(2.25) becomes ^2 = (kVf) + A (Se 5% CTiT) , *i dt x\ \ ) \ )' ^ = -SCse~7n, T > Ta. at Multiplying through by tx/Tx and tx respectively, we have the equations in the dimensionless form, dT t\k ~ /~ ~\ txA - = -Vv VT J + ASe dt xf \ ) Ti dS -JL. -^ = -hSCse TT^, dt T > 0, ACtiT, which depends on the six dimensionless coefficients, h Cs, (2.27) 2 = d'V (Vf) + c3Se~c*/f cbT, 2 = -c&Se~Cilf, f > 0. dt _ tik _Ta txA B ^1 2 ^2 rp i Qj rji ) ^4 rj~) 1 ^5 Cg ^1 Tl Tj 7l so (2.24)-(2.25) becomes We choose the scales aq, fi, and 7\ from what we expect the dominant mechanisms to be and to simplify the equations. It is natural to expect the temperature scale to be governed by the reaction rate, thus 7\ = B, and so 30 C4 = 1. Next we scale the time so that c3 = 1, hence t\ T\/A = B/A. Finally, we expect that the spatial scale is governed by the diffusion term and scale to make c\ = 1, which gives aq = (kt\)^2 = kxl2Bxl2 A~1!2. Then c5 = ACt\ ACB/A and c$ tiCs = BCs/A. In conclusion, the scaling (2.26) with Tj = B, Xl = kl/2B^2A~^2, h = BA~\ (2.28) transforms (2.24-2.25) into the non-dimensional form = V (vf) + Se~1/f AT, dt V / S = -(3Se-Vf, T> 0, dt with two dimensionless coefficients A = CB, (2.29) (2.30) (2.31) Therefore, the qualitative behavior of the solution is determined only by the nondimensional coefficients A and (3, which can be varied independently. The nondimensional form (2.29-2.30) suggests a strategy for identification of the coefficients k, A, B, C, Cs first match nondimensional properties of the traveling combustion wave, such as the remaining fuel fraction, or the ratio of the width of the leading edge to the trailing edge. As we see from Figures 2.6 and 2.7, the relationship between nondimensional parameters, A 31 x 10' 0.06 x 10' Figure 2.6: The fuel fraction remaining after the combustion wave as a func- tion of dimensionless parameters A and 0. 32 0.05 1 0.01 0.02 0.06 Figure 2.7: Ratio of the leading and the trailing edge of combustion wave with respect to A and /3 from (2.31). The width is measured as the length of the part of the wave above temperature 500/C. 33 and ft, and nondimensional properties (e.g. remaining fuel fraction, the ratio of the width of the leading edge to the trailing edge) of the combustion wave is approximately linear. This relation helps to get an initial temperature curve. Once we have found the initial temperature curve, we can optimize the two curves by con- structing a distance function between two curves. We used FMINSEARCH in MATLAB, which minimizes the distance function starting at an initial value with parameters. The optimized measured temperature curve can then be found by this process. In Figure 2.6, the narrow range, which was perturbed by very small amounts from the optimized nondimensional parameters, of A and /?, was used in order to have an approximately linear relation. Thus, we see that the fuel remaining is very small. The width of the wave can be measured e.g. as the distance of the points where the temperature equals 50% of the maximum with respect to A and /3, and this measured width is used in Figure 2.7. The nondimensional traveling wave solution T(t,x), S(t,x) has some (nondimensional) maximal temperature Tmax, width w, and speed v, while the data (a measured temperature profile) has the maximal temperature Tmax, the width ui, and the speed v of the traveling wave. This determines the scales T, = w vw X\ = , and 11 = ^3. w vw 34 3.8-, Figure 2.8: Width of the trailing slope as a function of A and j3. The rela- tion between the wave speed and (diffusion coefficient)0'5, i.e. \fk. We see approximately linear relations. 35 The substitution ~ x x = and into the system (2.24-2.25) with the coefficients A = Ti/ti, B = Tlt C = X/TU Cs P/h, and k = x\/ (Tfh) which has the desired nondimensional properties as well as the correct max- imal temperature, width, and speed of a traveling combustion wave. 2.4.4 Related work A dimensionless system similar to (2.29-2.30) was studied in Weber et al. [1997] when A = 0, i.e., combustion insulated against heat loss. The speed of the traveling wave as a function of /? was determined numerically and by an asymptotic expansion in Weber et al. [1997]. The observed existence of a traveling combustion wave only for small values of (1 is also shown in Weber et al. [1997]. When increasing /?, the solution becomes periodic, then the period doubles in Figure 2.9. We have observed that increasing /? has a (2.32) admits the scaled solution and 36 similar effect for equations (2.29-2.30) when A > 0. Also, we have observed that a sustained combustion wave is possible only when A is small enough. The dimensionless parameters e, A, and (3 were introduced in Mercer and Weber [1995], Mercer et al. [1996], Weber et al. [1997], and Asensio and Ferragut [2002], as we discussed in Section (2.3.2). The derived dimension- less equation can still have the ubiquitous cold boundary problem due to a natural cooling term, depending on how the dimensionless equation was de- rived (Weber et al. [1997]). By setting the ambient temperature to absolute zero in the natural cooling term, the cold boundary difficulty can be avoided. However, none of the paper analyzed with the natural cooling term at their work. In a different way, the derived dimensionless equation of (Weber et al. [1997] and Mercer and Weber [1997]) avoid the cold boundary problem in the heat loss term (Mercer and Weber [1995], Mercer et al. [1996], and Asensio and Ferragut [2002]), but it still requires extra work in dealing with new dimensionless parameter e, besides A and (3. Our model only has two di- mensionless parameters, A and (3. Having fewer parameters in the equations is a huge benefit for understanding the model. 2.5 Calibration of coefficients in one dimen- sion We have found the initial values B = 5.5849 x 104K, and C = 5.9739 x 10-4A-1 by using the values Tj = 670A and Tc = 1200A' in equation (2.22), 37 (a) (b) (c) (d) Figure 2.9: Periodic patterns of the solutions by varied (3. (a) (3 = 0.4172, (b) (3 = 0.4636, (c) (3 = 0.4945, (d) (3 = 0.5718. Each graph shows the temperature profile by timeframe. 38 1000 Figure 2.10: Time-temperature profile (dotted line) measured in a grass wild- land fire at a fixed sensor location, digitized from Kremens et al. [2003], and a computed profile (solid line) from simulation. and then the value A = 1.5217 x lO1/^-1 from (2.23), using the value tc 110s from Kremens et al. [2003]. Not every initial condition gives rise to a traveling combustion wave (Mercer et al. [1996]). Inspired by Mercer and Weber [1997], we use an initial condition of the form T(x,t0) = Tce~^x~x^+ Ta, where xq is in the center of the interval and a = 10\/2m.. This initial condition is smooth. Thus, it does not excite possible numerical artifacts. It has numerically local support and for a mod- est a provides ignition sufficient to develop into two sustained combustion waves traveling from the center. We then found empirically suitable values of k and Cg that result in traveling combustion waves, computed the nondimensional coefficients A and /3, adjusted them using Figure 2.10 (dotted line), and scaled using (2.32) to 39 match the maximal temperature and the width of the wave in Figure 2.10 (dotted line), and the speed of the traveling combustion wave, 0.17m/s, from Kremens et al. [2003]. There was a small amount of wind in Kremens et al. [2003], however we have not considered the wind here. The resulting coefficients are k = 2.1360 x 10~1m2s~1 K~3, A = 1.8793 x 102/fs-1, B 5.5849 x 102K, C = 4.8372 x lO"5^"1, and Cs = 1.6250 x lO-1*-1. The corresponding nondimensional coefficients were A = 2.7000 x 10-2 and 0 = 4.829 x 10_1. The computed traveling combustion wave (Figures 2.4, 2.5 and 2.10, solid line) is a reasonable match with the observation (Figure 2.10, dotted line). The trailing edge of the computed temperature profile (Figure 2.10, solid line) was not so well matched but this model is quite limited and other matches reported in the literature are similar (Balbi et al. [1999]). The real data looks like the superposition of two exponential decay modes, possibly the fast one from cooling and the long one from the heat stored in water in the ground. We have also noted that when the ratio A/Cs increases, the temperature in the traveling combustion wave increases, increasing the thermal diffusiv- ity coefficient k increases the width and the speed of the combustion wave and that the maximum temperature in the traveling wave decreases if Cs increases. Sufficiently small value of Cs is needed for sustained combustion. We have noted that the numerical solution by finite differences becomes un- stable when the ratio k/h, where h is the mesh size, is too small. 40 Chapter 3 Numerical methods for the PDE model 3.1 Efficient numerical implementation There are several approaches for solving systems of nonlinear PDEs. We used a standard finite difference method, with an upwind scheme in the space variables for stability purposes. The trapezoidal method in time with Newtons method was implemented. An explicit method, such as Runge- Kutta, would be much more expensive than the trapezoidal method because a very small time step is required due to the diffusion and reaction terms. In Newtons method, the GMRES method with FFT as a preconditioner was also used. The discretized systems of equations (2.2-2.3) in space is ^ 41 where H~1----S+Ul + ^) (max(t;x,0)uM + min(t;x, Oju^) B (max(vy, 0)ui j + min^, 0)rt^) + Au2je U1~T AC{u\j Ta),
-Csu2e
speed R(x,y,t), the spread rate, in the normal direction. The speed may
= IWr
distance function, but this property is not preserved in general by advancing
method, which is a central difference approximation in space, [Osher and
)-ty(
Ax
and VT
Ax Here G\ and G2 are partial derivatives of G with respect to ipx and
spectively. Then we have the advection equation in two dimensions, G\ = VX
n+1 0 . |