Percolation Phenomena in Filtration Combustion - Industrial

3008

Ind. Eng. Chem. Res. 2004, 43, 3008-3018

Percolation Phenomena in Filtration Combustion Chuan Lu and Yannis C. Yortsos* Department of Chemical Engineering, University of Southern California, Los Angeles, California 90089-1211

Using a pore-network model, we study the effect of randomness in the solid fuel distribution on the propagation of filtration combustion fronts in porous media. Both a complex pore-network model for full filtration combustion and a more simplified gasless combustion model are used. We find that the minimum threshold (akin to a “percolation threshold”) for the fuel distribution is a function of the kinetic and heat-transfer parameters and can vary over a wide range. The balance between heat transfer and reaction kinetics endows the system with a behavior more complex than simple percolation, which was earlier suggested as a model for such systems. The long-range effect of heat conduction acts to move the behavior away from true percolation behavior. Conversely, the effect of heat losses is sufficient to lead to extinction (thus, to “nonpercolation”) even for systems with fully occupied fuel sites. In either case, the availability of the fuel and its distribution in space are important parameters in the process. 1. Introduction Filtration combustion (FC) is an area of significant interest in various applications, from the processing of materials to the recovery of oil from oil reservoirs.1 It involves the combustion of a solid fuel deposited at the pore walls of a porous medium, through the injection of an oxidizing agent. When ignition occurs at the gas inlet, a forward FC process develops, in which reaction and thermal fronts propagate in the same direction as the injected gas. When ignition occurs at the opposite side, the combustion front propagates in the direction opposite to the gas flow, and the process is denoted as reverse FC. A number of studies have addressed various features of FC. Analytical investigations using continuum models of forward FC have described properties of the propagation front2,3 and the effects of heat losses, including multiple steady-state behavior4 and front instabilities.5 Reverse FC has been explored experimentally6 and by some simple numerical models.6,7 In a recently completed work,8 we investigated the properties of FC at the pore-network scale. Although effective continuum models are very valuable, investigations at the pore-network scale provide information that cannot be captured by continuum models. This is particularly important when fronts are relatively thin, in which case an effective medium approach might be questionable. It is also important for understanding the effect of pore-scale heterogeneity, when this heterogeneity is dynamically coupled to the process. In a series of companion papers, we have addressed many of these issues by analyzing forward FC9 or reverse FC10 using random pore-network models. Pore-network representations of porous media have become the most commonly accepted approaches for effectively capturing important geometrical and topological features of the microstructure, and they are currently widely used. In filtration combustion, an important effect of the pore-scale heterogeneity is the fuel distribution. Depending on the method by which it was generated or deposited on the pore walls, the fuel is likely to be heterogeneously distributed in the pore space. For * To whom correspondence should be addressed. E-mail: [email protected].

example, in in situ combustion in porous media, the fuel is generated in situ by oxidation reactions preceding the main combustion fronts.4 This process is controlled by a combination of transport and reaction processes and, thus, by pore-scale heterogeneity. Specifically, fuel formation occurs through a combination of low-temperature oxidation reactions that precede the propagating combustion front. The process is quite complex and involves a number of coupled processes (see Prats1 for further details). Of specific importance to this work is the distribution of the fuel as a binary variable between pores, in which case percolation-like phenomena are possible.11 Although a rather limiting case, such a distribution is useful for addressing more generic issues. We are specifically concerned with ignition and extinction phenomena, which are determined by nontrivial balances involving rates of heat transfer and reaction, as well as by the pore microstructure. The application of percolation-type approaches to combustion processes is not new and has been introduced before in the context of forest fires. The spreading of these fires is affected by the amount and distribution of vegetation; the convection strength (the wind and the field slope); the moisture; and other parameters, such as radiation. Experimental evidence has suggested that a minimum vegetation density must exist for the fire to spread. The understanding of this parameter can be of significance for effective fire-fighting strategies and has been pursued in the literature before. Although a number of simplified models (e.g., Morvan and Dupuy12) have been developed, the full problem presently remains unsolved. However, considerable progress has been made using percolation approaches, in which an analogy is made with phase transitions. In this analogy, the density of fuel is the control parameter, and the condition of whether fire can spread represents the order parameter. A number of studies have been proposed to relate the phase transition in fire spreading with ordinary site percolation, as discussed below. Albinet et al.13 conducted simulations in a 2-D square array, randomly occupied by combustible material, that is subsequently ignited at one end. Ignition of blocks containing fuel followed simple ignition rules (first or second neighbor). To interpret the results, scaling

10.1021/ie0306372 CCC: $27.50 © 2004 American Chemical Society Published on Web 03/12/2004

Ind. Eng. Chem. Res., Vol. 43, No. 12, 2004 3009

relations were proposed according to ordinary site percolation theory. For example, the characteristic length for the spatial extent of the “fire spreading” cluster was set equal to that of the well-studied percolation cluster

ξ ∝ ||-ν

(1)

and the characteristic time for growth was taken equal to

θ ∝ ||-τ

(2)

Here, ) (p - pc)/p is the control parameter, p is the fraction of sites occupied by fuel, pc is the percolation threshold, and ν and τ correspond to critical exponents from ordinary percolation theory.11 Using this model, the percolation threshold for a 2-D system is 0.593 or 0.407 for the case of first-neighbor or first- and-secondneighbor ignition, respectively. Equations 1 and 2 further imply that the mean front position and the total number of burned sites would scale with time as

xj ∝ tv/τ

(3)

M ∝ t(v-β)/τ

(4)

and

respectively, where β corresponds to the exponent of the percolation probability in site percolation.11 Tephany et al.14 conducted experiments to mimic the above simulations. The occupied sites were filled with wood shavings; the unoccupied sites left empty or filled with fiber cement, wood shavings, treated with fire retardant. The percolation threshold and the dependence of the total number of burned sites on time agreed with the above. However, previous experimental work did not show a similar agreement. Beer15 conducted experiments on arrays of matches vertically inserted into a wire mesh. In the resulting fire, strong heattransfer effects were observed. More than two adjacent burning sites could ignite a second or even a third neighbor, where a single, isolated burning site was incapable of doing so. The conclusion was reached that the combination of percolation theory and simple ignition rules was perhaps too simplistic to model fire propagation in a heterogeneous system when heattransfer effects are important. The fuel percolation threshold in FC can be determined by detailed simulations using a full pore-network model, such as one of the type developed by Lu and Yortsos,9 as shown below. Because of the computational demands in investigating the effects of various parameters, however, it is more instructive to simplify the description by using a simpler network model. In this study, we consider such a simpler model and investigate the comparative effects of heat conduction, but not convection, and of reaction kinetics on percolation-type phenomena. We show that, although a percolation-like threshold in the fuel density exists, the process is not a typical percolation process. Rather, depending on the ignition characteristics of the system, the process can be dictated by longer-range interactions between sites. In such a process, the balance between kinetics and heat transfer is an important consideration. Percolation is then the sustained propagation of a combustion front.

Figure 1. Thermally coupled pore network and inert solid matrix: (A) schematic of the pore network, (B) the coupling between the pore network and the solid lattice. White denotes the pore network, black the inert solid. The dashed lines indicate pore-solid interactions. The one-to-one assignment between solid sites and pore sites is made for simplicity of the calculations.9

The paper is organized as follows: First, we use a full pore-network model for the forward FC process under typical conditions, but with fuel distributed in a binary fashion, to illustrate the existence of a percolation threshold and to obtain its dependence on thermal conductivity. Because of the high computational cost, a simpler model in which convection is not a factor is subsequently introduced. The comparative effects of heat conduction, heat losses, and chemical kinetics on the percolation threshold are analyzed using models of increasing complexity. We close with some remarks on the percolation nature of the process. 2. Percolation Threshold of Forward FC in a Pore Network To demonstrate the existence of a percolation threshold in the fuel fraction, we used the pore-network model of Lu and Yortsos.9 The model consists of a network of pores and throats, representing the pore space, that is embedded in another network of solid sites, representing the solid matrix.16 The pores are the places where reaction occurs and solid fuel exists. They are interconnected via throats (bonds) that control the transport of mass, momentum, and heat. Heat transfer between pore and solid sites leads to the coupling of the two networks of solid and pore space, as shown in Figure 1. In the 2-D simulations reported below, both networks are square lattices. Important features of mass and heat transport include the following: 1. In two dimensions, each solid site communicates with four solid sites (Figure 1B) and one pore site (the dashed line between solid and pore sites in Figure 1B) through heat conduction. The one-to-one assignment between solid sites and pore sites (for example, contrast with Figure 1A) is made for simplicity in the calculations. 2. Injection of a mixture of oxidant and inert gas occurs at one end, where the pressure, temperature, and composition are specified. The outlet end is at a constant pressure. Flow in the pore throats is governed by Poiseuille’s law. The oxidant in the gas phase is transported by gas-phase diffusion and convection. Ignition occurs at the gas inlet.

3010 Ind. Eng. Chem. Res., Vol. 43, No. 12, 2004

3. The following one-step heterogeneous reaction is assumed

solid fuel + µ(gas oxidant) f solid product + µ(gas product) (5) where µ is a stoichiometric coefficient and we assume no net gas production. Implicit in the above expression is the assumption that the oxidant is O2, the fuel C, and the gas product CO2. Of course, we are using reaction 5 to more generally refer to a lumped system. The reaction rate is expressed by the kinetic model

( )

R ) krArPXoH(Vf) exp -

Ea RT

(6)

where H is the step function, Vf is the volume of fuel, P is the total pressure, kr is a kinetic constant, Ar is the interfacial area in an individual pore site, Xo is the mass fraction of the oxidant, and Ea is the activation energy. 4. Within a pore site, the concentrations, pressure, and temperature are uniform. This implies that fuel and gas are at thermal equilibrium. However, heat transfer does take place between adjacent pore sites, between pore and solid sites, and between adjacent solid sites. The governing equations in the pore-network model are discretizations of the relevant conservation equations, and they are detailed in Lu and Yortsos.9 Along with dimensionless similarity geometric parameters, a number of dimensionless parameters control the process. These include the thermal Thiele modulus, Th ) kl2/R; the mass Peclet number, Pe ) ul/D; the Nusselt number, Nu ) hl/λ; the Lewis number, Le ) R/D; the Arrhenius number, γ ) Ea/RT0; and the dimensionless heat of reaction, q ) φfFf∆H/[(1 - φ)FCpT0] (Akkutlu and Yortsos4). Here, k is the reaction constant, l is a local pore-length scale, R is the solid thermal diffusivity, u is the injection velocity, D is the mass diffusivity, h is the heat-transfer coefficient to the surroundings, λ is the thermal conductivity of the solid, R is the initial temperature, ∆H is the heat of reaction, FCp is the mass heat capacity of the rock, Ff expresses the mass density of the fuel, and φf is the volumetric fraction of the fuel (generally a small number). The Thiele modulus here expresses the ratio between thermal diffusion and kinetic time. The Nusselt number expresses heat losses to the surroundings. Because of the 2-D geometry, the heat-loss model corresponds to a thermally thin problem, in which the heat transfer across the thickness of the medium is sufficiently fast compared to the rate of heat loss . For the general problem, a 3-D description is more appropriate. To study percolation phenomena, the fuel was randomly distributed in the sites of the pore network based on a binary scheme (present or absent, as shown in Figure 2, where the pore throats connecting adjacent sites have been omitted for simplicity). Each site o Ff. Thus, the measure contains the same fuel mass, Vi,f of fuel availability is expressed through the occupancy fraction p. In the simulations below, a 2-D pore-network of size 64 × 64 was used. The system was ignited by assigning a sufficiently high temperature (the adiabatic front temperature, which in rescaled notation using T0 is equal to (1 + q) to the first column at the upstream end. The percolation threshold pc was determined by monitoring whether the reaction front reaches the other end (Figure 3A) or becomes extinct before doing so

Figure 2. Schematic of the fuel distribution in the pore network. Grey indicates the presence of fuel, white the absence of fuel. The pore throats connecting adjacent sites have been omitted for simplicity.

Figure 3. Conversion patterns in forward FC processes in the pore network with distributed fuel. Examples of (A) a percolating front (p> pc), (B) a front that extinguishes before reaching the end (p < pc). Gray indicates burned fuel, white indicates inert sites, and black indicates unburned fuel.

(Figure 3B). Thus, strictly speaking, it does not correspond to a simple percolation process, but rather, it also involves issues of extinction as a result of heat transfer. The threshold was determined by a bisection method, where, to maintain accuracy, 10 bisections were run to determine the threshold corresponding to the same seed for the random realization. The numerical


Figure 4. Fuel percolation threshold (minimal fuel occupy fraction) plotted vs the thermal conductivity for a nonadiabatic forward FC process in a 64 × 64 2-D square lattice using the full pore-network model. The parameters from Table 1 were used with Nu ) 0.04. Table 1. List of the Parameter Values Used parameter

value

Q, heat of reaction Ea, activation energy k0, kinetic constant T0, initial temperature u, injection velocity (SS) Pi, injection pressure Yi, injected mole fraction D, gas-phase diffusivity λ, maximum conductivity φ, porosity Fg,i, gas density cg, gas mass heat capacity Fs ) Ff, solid and fuel densities, respectively cs ) cf, solid and fuel capacities, respectively φf, volumetric fuel fraction µg, stoichiometric coefficient

39 452 kJ/kg 73501.6 kJ/kmol 498 s-1 atm-1 373 K 10 mm/s 10 atm 0.23 2 × 10-6 m2/s 0.8655 J/m‚s‚K 0.3 1.2254 kg/m3 306.8386 kJ/kg‚K 2.5 × 103 kg/m3 0.35 kJ/kg‚K 0.0256 1.00

values obtained were then averaged over a number of different seeds. This resulted in a high computational cost. Because of the high computational requirements, as well as the complexity of the process, we studied only one effect, namely, the dependence of the percolation threshold on the thermal conductivity. The results are shown in Figure 4. These results correspond to one realization only. However, they are indicative of a general trend, as also discussed below. The parameters of Table 1 were used under conditions of local heat conduction control (Pe ) 0.01, Nu ) 0.04). It is clear that, for this set of parameter values, a definite percolation threshold exists that ranges between 0.4 and 0.9 in the parameter space investigated. The threshold reflects the balance required among the rate of heat transfer, the rate of heat losses, and the rate of heat release for combustion to be sustained throughout the

domain. We note that the threshold has a nonmonotonic dependence on the conductivity, increasing at low and high values of the thermal conductivity and displaying a generally wide plateau at intermediate values. At low conductivities, increasing the thermal conductivity value increases the probability that neighboring sites containing fuel will be ignited by an adjacent burning site. Thus, the threshold initially decreases with the thermal conductivity. As the conductivity becomes too large, however, the thermal energy of a burning site diffuses quickly, such that adjacent sites might not reach ignition, and the process rapidly becomes extinct. It is rather surprising that percolation thresholds as high as 0.85 are found. As discussed below, such high values are due to the nonadiabatic conditions assumed. It appears that, at high thermal conductivities, sites not containing fuel act as a thermal sink of such strength that this effect is sufficient to delay and ultimately extinct the propagation of a combustion front. As a result, the spacing between fuel sites must be reduced for sustained propagation to occur. The general behavior depicted here also occurs in the computationally simpler gasless system described below. It should be noted that, in the full model, extinction can occur only if heat losses are included. Heat losses have been shown to lead to extinction under certain conditions in 1-D continuum models,4 where fuel exists continuously in the system (i.e., even for pc ) 1). One concludes that FC in a system where the fuel is distributed in the pores in a binary, random fashion is not a simple percolation problem, in which only first or second nearest neighbors are ignited. Rather, heat transfer, heat losses, and reaction kinetics combine to endow the system with a more complex behavior. Significantly, the thermal conductivity plays an important role in determining the percolation threshold. To analyze the process in more detail, we examine in the following section a simpler model in which the reaction is controlled by the fuel availability and is independent of the oxidant supply. 3. Percolation Threshold of Forward FC under Solid Fuel Control (Gasless Reaction) We consider a simplified model of FC controlled by the solid fuel (gasless reaction). In this model, we assume that the gas-phase oxidant is present when needed and at the required supply rate and that the reaction kinetics are zeroth order, independent of the oxidant concentration. The objective is to decouple the oxidant and gas-phase balances from the energy balance, which here is controlled only by fuel availability and heat conduction. Fluid flow and transport in the gas phase are not issues, and consequently, the process simulates “gasless combustion”.17 Strictly speaking, therefore, this approach represents only a limiting case of FC. We use a lattice model, as before. However, the pore lattice contains only an immobile fuel, and the assumption is made that the heat capacity of a site changes little after combustion. Clearly a simplification, this assumption allows us nonetheless to reach some useful insights into the process. This situation is similar to the least significant bit (LSB) model of Tephany et al.,14 where heterogeneity is obtained by mixing combustible and noncombustible blocks. The combustible blocks were made of sawdust with paraffin, and the noncombustible consisted of fiber cement slabs. The blocks were placed randomly on a square lattice.


The governing equations for this system are simpler than those required in the complete pore-network FC model of Lu and Yortsos,9 described above. The energy balance in a pore site i is given by

∆Ei ∆t

)-

k l

∑j Aij (Ti - Tj) - hLAi(Ti - T0)

(7)

where Ei ) [(1 - φ)FCps]Tj + Ffφf∆rH is the energy of site i, l is the (pore) length between two neighboring sites, A stands for the heat-exchange area, hl is a heatexchange coefficient between a pore site and the environment, the subscript j stands for nearest-neighbor pore sites, and the subscript 0 stands for the environment. In the above formulation, heat conduction occurs directly between pore sites (although physically through the solid); hence, we used only a one-lattice network (rather than the two-lattice network used in the full model above). Assuming zeroth-order kinetics, the fuel mass balance equation for site i is expressed as

( )

Ea ∆Vi,f kr Ar ) exp H(Vi,f) ∆t Ff RTi

(8)

with various sizes showed that this lattice size was sufficiently large. A set of 10 different realizations were taken for each parameter set, and the percolation threshold was determined as the arithmetic average of these realizations. In the simulations, the process is ignited by assigning to the first column a temperature slightly higher than the ignition threshold, e.g., (1 + a)∆Ta, where a > 0 is a small parameter (typically, a ) 0.05). At any given value of p, we determine the minimum value of the ignition temperature, or more precisely of the ratio ∆Ta/∆Tr ≡ ∆Ta/qT0, needed to propagate the front. Thus, we can either fix the ratio ∆Ta/∆Tr between the ignition and reaction temperatures and vary p until no further propagation is possible or fix p and reduce ∆Ta/∆Tr until the process becomes extinct. An estimate of the upper bound of the ignition temperature can be obtained by considering the limit in which fuel occupies the entire lattice (i.e., pc ) 1). An analytic result can be obtained if we consider the effective continuum limit. Assuming that ignition occurs in the entire first row, the problem is essentially that of 1-D heat conduction in an infinite medium, described by

and, in terms of the reaction fractional conversion η ) o 1 - Vi,f /Vi,f , as

∂Θ ∂2Θ ) 2 ∂t ∂x

Ea ∆η ) k exp H(1 - η) ∆t RTi

with vanishing temperature in the far field and with the initial conditions

( )

(9)

where we lumped the reaction area, the fuel density, and the kinetics into one overall kinetic constant, k ) o kr Ar/V i,f Ff, with units of inverse time. Here, Vi,f is the fuel volume in a site, and the superscript o denotes initial conditions. All variables except for the fuel, which was randomly distributed, are uniform. Ignition is imposed at the inlet row only. To investigate the percolation properties of this system, we considered the following specific cases. 3.1. Infinite Reaction Rate. The simplest model case is one in which, following ignition (that is, after the temperature exceeds a threshold), the kinetics is infinitely fast and does not control the process at all. To model this situation, we applied the following ignition rule: If a site contains fuel, it will be ignited when its temperature exceeds a specified ignition threshold ∆Ta + T0. Upon ignition, the temperature of the site is instantly increased by the adiabatic reaction temperature ∆Tr ) φfFf∆H/(1 - φ)FCp, as noted above (see also Akkutlu and Yortsos4), and simultaneously the fuel amount in the site instantly becomes zero, because of the assumed infinitely fast kinetics. The ignited site, then, dissipates its temperature through heat conduction to the other sites and to the surroundings. This process is described by two dimensionless parameters, the dimensionless ratio between the ignition and reaction temperatures, ∆Ta/∆Tr, and the heat losses to the surroundings, which are expressed through the Nusselt number. (i) Effect of the Ignition Temperature (∆Ta/∆Tr) in an Adiabatic System (Nu ) 0). We conducted simulations of this problem in an adiabatic system in a 128 × 128 square lattice, using no-flux boundary conditions on the sides. The lattice size was selected on the basis of two considerations: limiting the finite size effect but also making computations feasible. Simulations

Θ)

{

∆Ta 1 1 + 1 for - < x < ∆Tr 2 2 elsewhere

(1 + a)

0

(10)

(11)

In the above equation, the temperature was made dimensionless with the dimensionless reaction temperature ∆Tr, distance with the pore length l, and time with the conductive time. In this notation, the adiabatic temperature is equal to unity. Equation 11 expresses the fact that, once a site is ignited (in the interval -1/2 < x < 1/2), its temperature instantly increases by the adiabatic temperature. The solution of eqs 10 and 11 is readily found to be

[

][

Θ(x,t) )

( ) ( )]

∆Ta 1 1 +1 -x x+ ∆Tr 2 2 erf + erf 2 x2t x2t (12)

(1 + a)

By evaluating at the nearest neighboring site (x ) 1), we can then determine whether the temperature at the nearest neighbor will surpass the ignition threshold, which is a necessary and sufficient condition for the propagation of a combustion front in the porous medium. The solution Θ(1,t) in eq 12 has a nonmonotonic behavior, reaching a maximum where

∂Θ(1,t) ≡ ∂t

[

][

∆Ta +1 ∆Tr

(1 + a)

4x2πt

3/2

( 8t1 ) - 3 exp(- 8t9 )] ) 0 (13)

exp -


Figure 5. Comparison between the analytical and discrete results for determining the upper bound on the ignition temperature for the case p ) 1. The discrete version shows a slightly higher temperature response at the nearest-neighbor site at x ) 1, when ignition occurs at site x ) 0. Case of an infinite reaction rate in an adiabatic gasless system.

This occurs at the time t* ) 1/(ln 3). For successful ignition of the nearest neighbor, the corresponding maximum must exceed the ignition temperature, that is, the following inequality must be satisfied

(

1+a+

)[ (

) ( )]

∆Tr xln 3 3xln 3 erf - erf ∆Ta 2x2 2x2

g 2 (14)

For small a, this requires the condition ∆Ta/∆Tr < 0.3247. Therefore, we anticipate that there would be an upper bound on the ratio ∆Ta/∆Tr above which sustained ignition is not possible, even for a fully occupied lattice. The existence of an upper bound was verified in the numerical results reported below. However, its actual value was somewhat different, due to the fact that the simulation was on a discrete, rather than continuum medium. In this case, the relevant equations are different, the equation having the matrix form

dΘ h ) AΘ h dt

(15)

where Θ h is the vector of temperature values at the site (of size 128) and A is the (128 × 128) square matrix

[

-1 1 1 -2 1 A) ‚‚‚ ‚‚‚ 1 -1

]

(16)

The solution of this problem gives a slightly higher estimate, as shown in Figure 5. We find that the more appropriate condition for this system is ∆Ta/∆Tr < 0.4627. This result was indeed confirmed in the 2-D results shown below. This analysis shows that, as also suggested in Figure 4, there exist bounds on the various controlling parameters, above (or below) which sustained propagation is not possible for reasons other than fuel scarcity. Strictly speaking, therefore, the behaviors exhibited by these systems are not percolation phenomena. When these bounds are reached, it is possible that the fraction in

Figure 6. log-log plot of the effect on the percolation threshold of the dimensionless ignition parameter ∆Ta/∆Tr, expressing ignition temperature or heat of reaction effects. Case of an infinite reaction rate in an adiabatic gasless system.

the fuel distribution is, in fact, smaller than unity, particularly for an adiabatic system. This situation is further analyzed below. In addition, there is a singular behavior when the bound is reached, in that, to the right of the limiting point, a solution is not possible. Figure 6 shows the effect on the percolation threshold of the parameter ∆Ta/∆Tr in the adiabatic case for different values of ∆Tr. Because of the adiabatic conditions, the percolation threshold should depend only on the dimensionless group ∆Ta/∆Tr. There is a slight effect of ∆Tr, attributed to the finite size of the lattice and the small number of realizations. As expected, the percolation threshold is a monotonically increasing function of the ratio ∆Ta/∆Tr. Indeed, if the solid fuel is more energetic, releasing more heat (higher ∆Tr), the ignition temperature is reached faster at larger distances from the ignited site; hence, a combustion front can be propagated with more sparsely distributed fuel. Clearly, an increase in the heat of reaction will reduce the percolation threshold. Conversely, for a fixed heat of reaction, a decrease in the ignition temperature will reduce the threshold. The two limiting regimes in the figure can be further analyzed. At the limit where no ignition threshold is required, the percolation threshold must be zero, as any temperature increase would be sufficient to ignite the system. Therefore, we expect pc f 0 as ∆Ta/∆Tr f 0, as shown in Figure 6. The opposite limit corresponds to the case where the amount of fuel is sufficient for nearest-neighbor percolation of fuel-occupied sites, but the ignition threshold is too high to allow for successful ignition. An estimate of this limit was obtained above using a 1-D analysis and assuming that all sites are occupied with fuel (pc ) 1). Figure 6 shows that the upper bound estimate of 0.4627 is not violated, even for the more general 2-D disordered system (where p < 1). At the point where this bound is reached, we note that the percolation threshold can be smaller than 1. Even though the 1-D analysis is not expected to apply rigorously to the more sparsely distributed system, it appears that its predictions are quite adequate.

3014 Ind. Eng. Chem. Res., Vol. 43, No. 12, 2004 Table 2. Comparison of the Critical Time and the Upper Bound on the Ratio ∆Ta/∆Tr as a Function of the Nusselt Number for the Case of Infinite Reaction Rate in an Nonadiabatic Gasless Systema 1-D effective discretization continuum with in a lattice with full fuel occupancy full fuel occupancy

Figure 7. Effect of the Nusselt number on the dependence of the percolation threshold on the ratio ∆Ta/∆Tr. Case of an infinite reaction rate in a nonadiabatic gasless system. The various symbols correspond to Nu ) 0.00, 0.05, 0.10, 0.15, 0.20, and 0.25 (see legend). The numerical prediction is indicated by the solid line.

Figure 6 shows the regression of the data in a loglog plot. The resulting curve is close to a straight line with exponent 1. We conclude that, in the case of infinitely fast kinetics, the percolation threshold is a linear function of the ratio ∆Ta/∆Tr, which expresses ignition vs heat release capacity. The threshold vanishes when ignition is not required and approaches an upper limit when the ratio reaches an upper bound. Beyond this bound, sustained propagation is not possible, but for reasons unrelated to the fuel distribution in the lattice (always assuming, of course, that the pore length remains fixed and that the fuel resides only in the sites of the lattice). We should note that the results do not depend on the thermal conductivity, as expected. Indeed, various different values tested showed very little dependence.8 (ii) Effect of the Ignition Temperature (∆Ta/∆Tr) in a Nonadiabatic System. In the presence of heat losses, the overall level of the percolation threshold is expected to increase, as higher rates of heat losses will lead to faster heat dissipation and, thus, more difficult conditions for the ignition of distant sites. Figure 7 shows the effect of ∆Ta/∆Tr on pc for different values of the Nusselt number. As in the previous section, we expect that the percolation threshold will vanish in the absence of an ignition threshold and that it will reach a limit at a specific upper bound value ∆T/a/∆Tr. Above this upper bound, the ignition threshold is too large to be overcome, even for densely packed fuel. Previous studies on FC have shown that heat losses can lead to multiplicity and extinction.4,9 In those works, however, the dominant variable considered was the injection velocity. The limiting value of the percolation threshold, when the upper bound ∆T/a/∆Tr is reached, decreases with the Nusselt number. Such behavior is anticipated in general terms, as larger rates of heat loss impose stronger requirements on the ignition thresholds. Between the two limits, the behavior appears to be close to a linear dependence, with the slope increasing with increasing Nusselt number. We can obtain an approximate estimate on the upper bound by extending the previous 1-D analysis to include heat losses. In the

Nu

t*

upper bound

0.00 0.05 0.10 0.15 0.20 0.25

0.9102 0.8261 0.764 0.7151 0.6751 0.6415

0.3247 0.3065 0.291 0.2774 0.2653 0.2544

simulation of the full problem

t*

upper bound

percolation upper threshold at bound upper bound

1.2240 1.0860 0.9880 0.9120 0.8520 0.8000

0.4627 0.4253 0.3949 0.3695 0.3478 0.3288

0.4623 0.4252 0.3947 0.3693 0.3478 0.3287

0.82 0.88 0.92 0.93 0.95 0.958

a 1-D effective continuum with p ) 1, 1-D numerical discretization with p ) 1, and the full simulation in a 2-D disordered system (pc < 1) are shown. The numerical results are in good agreement with the full simulation.

case of a 1-D geometry, the heat-transfer equation in dimensionless notation is

∂Θ ∂2Θ ) 2 - NuΘ ∂t ∂x

(17)

which, through the transformation

Θ ) exp(-Nut)φ

(18)

can be mapped to the heat equation

∂φ ∂2φ ) ∂t ∂x2

(19)

The solution of the latter is the same as in the adiabatic case. Now, however, the time for the nearest neighbor to reach ignition also depends on the Nusselt number. After some algebraic manipulations, we find that this time, t*, solves the algebraic equation

( )

exp -

( ) [ (

1 9 - 3 exp ) 8t* 8t* Nu2x2π(t*)3/2 erf

) (

)]

3 1 - erf 2x2t* 2x2t*

(20)

Then, the condition on the ignition temperature becomes

∆T/a B e ∆Tr 2 - (1 + a)B

where

[ (

B ≡ exp(-Nut*) erf

) (

)]

3 1 - erf 2x2t* 2x2t*

(21)

The solution of eq 21 shows that the upper bound decreases as a function of the Nusselt number. Results from the analysis are reported in Table 2. Included in the rightmost columns are also the results from the full simulation of the problem. We note that the effective continuum underpredicts the upper bound by about the same amount as in the adiabatic case, although the trend of a decreasing upper bound with increasing Nusselt number is reproduced. The reason for the difference has to do with the finite discretization of the problem. To correct this discrepancy, we also calculated the upper bound by solving eq 17 numerically, as in the


Figure 8. Effect of reaction kinetics (through the Thiele modulus) on the dependence of the percolation threshold on the ratio ∆Ta/ ∆Tr for finite kinetics (nonadiabatic system with Nu ) 0.05).

previous section. Now, the equivalent matrix is

[

-1 - Nu 1 1 -2 - Nu 1 A) ‚‚‚ ‚‚‚ 1 -1 - Nu

]

(22)

The results of this computation are listed in the middle columns of Table 2. We note the very good agreement, this time, with the results of the full simulation, that is, of the problem where the 2-D disordered system was solved. The trend of decreasing upper bound with increasing heat losses is clearly reproduced. As in the adiabatic case, it is interesting that the 1-D results for the fully occupied lattice are reproduced for the 2-D problem and more sparsely occupied lattice (although, here, the corresponding occupancy fraction is quite close to unity in the extinction limit). 3.2. Finite Reaction Rate. Consider, next, the case of finite reaction kinetics. Now, there is no need to specify an ignition temperature, and the parameter ∆Ta is replaced by ∆Ta ) E/R ) γT0. In addition to ∆Ta/∆Tr and Nu, the problem is also controlled by the thermaldiffusion-based Thiele modulus, Th ) kl2/R, and the ratio q ) ∆Tr/T0. In the results shown below, the ambient temperature was fixed to T0 ) 300 K. Because of the finite reaction kinetics, the behavior of the system becomes more complex. Nonetheless, for the case with a high activation energy, the fronts percolate mostly by near-neighbor sites, as otherwise, the temperature is dissipated before any substantial reaction occurs. The effect of the kinetic parameter Th on the dependence of the percolation threshold on the parameter ∆Ta/∆Tr is shown in Figure 8 for a nonadiabatic system (Nu ) 0.05). In the adiabatic case, the percolation threshold should theoretically be zero for a finite-size lattice. Indeed, even if the reaction is slow, given a sufficiently long time, the system will ultimately reach combustion behavior. Therefore, meaningful percolation thresholds can be found only in the nonadiabatic case, where heat losses might prevent the system from reaching a sustained propagation. In the simulations corresponding to Figure 8, the process was initiated by igniting the first row at one-half the adiabatic temperature level. Different ignition temperatures might affect

Figure 9. Effect of heat losses (through the Nusselt number) on the dependence of the percolation threshold on the ratio ∆Ta/∆Tr for finite kinetics (Th ) 1).

a finite-size system. The parameter ∆Ta/∆Tr was varied by keeping q ) ∆Tr/T0 constant and equal to 1.02 and varying ∆Ta ) γT0. In addition, because of computational requirements, only three realizations were run, in a lattice of 64 × 64. Figure 8 shows that the dependence of pc on ∆Ta/∆Tr is similar to that observed for the previous case of fast kinetics. In the limit of zero activation energy (or zero ∆Ta/∆Tr), the percolation threshold is zero: the reaction is independent of temperature and occurs at a fixed finite rate, ultimately burning all of the fuel. As the activation energy increases, the percolation threshold also increases. The effect of the Thiele modulus is to decrease the percolation threshold. Such behavior is expected, as faster kinetics leads to a faster release of heat and, hence, higher temperatures in the system and the reduced possibility of extinction due to the constant heat losses. In the limit of very fast kinetics (high Thiele modulus), an asymptotic curve is approached. Figure 8 shows that this asymptote is reached when the Thiele modulus is close to an order-1 value. The limiting curve is almost the same as that for the case of infinitely fast kinetics (compare with Figure 7 and also note the different range in the ordinate axis). The effect of the heat losses on the dependence of the percolation threshold on the parameter ∆Ta/∆Tr is shown in Figure 9 for Th ) 1. As expected, the results are very similar to those presented in Figure 7: higher heat losses leading to higher percolation thresholds, with an almost linear dependence in the curves shown. Comparison between the curves corresponding to the same Nu value shows a consistent behavior (again note the different range in the ordinate compared to Figure 6). It should be pointed out that, in the simulations corresponding to both Figures 8 and 9, the range of the parameter ∆Ta/∆Tr was not constrained by an upper bound, which did not appear to exist in this case. We can demonstrate this by paralleling the 1-D analysis of Akkutlu and Yortsos,4 discussed below. Consider the 1-D propagation of a reaction front with full fuel occupancy. The corresponding equations now read

( )

∆Ta ∆Tr ∂Θ ∂ Θ H(1 - η) ) 2 - NuΘ + Th exp T0 ∂t ∂x Θ+ ∆Tr (23) 2


( )

and

∆Ta ∆Tr ∂η H(1 - η) ) Th exp T0 ∂t Θ+ ∆Tr

(24)

We look for traveling-wave solutions and follow the large-activation-energy approach of Akkutlu and Yortsos,4 in which expansions in terms of a large Zheldovich number are taken. Akkutlu and Yortsos4 considered the full FC problem, which includes convection and transport of the oxidant gas. However, their analysis is also applicable to the present case of gasless combustion. By closely following their approach, it is not difficult to show that the rescaled temperature of the moving front is the solution of the following equation

Θf2 )

1 4qγNu

1+

(

Th(1 + qΘf) exp -

(25)

)

γ 1 + qΘf

The dimensionless front velocity (made dimensionless with the conductive speed R/l), is given by the expression

(

Th γ (1 + qΘf) exp v ) qγ 1 + qΘf 2

)

1 4qγNu

(

Th(1 + b) exp -

(26)

(27) γ 1+b

)

After performing some algebraic manipulations and taking the limit of small b, the upper bound on γ is found to be the solution of the equation

γ exp(γ) )

qTh 4Nub2

on the percolation threshold. The overall dependence of the percolation threshold on λ can be expressed as

∂pc ∂pc ∂pc Th Nu + ) ∂λ ∂Th λ ∂Nu λ

(

)

(

)

(29)

Because ∂pc/∂Th < 0 and ∂pc/∂Nu > 0, the combination of parameters is such that the percolation threshold has the nonmonotonic dependence found. 4. Discussion

To define an extinction condition, we arbitrarily impose the condition that the front temperature is higher than the initial temperature by a small fraction, namely, that Tf ) Θf∆Tr + T0 g(1 + b)T0, where b > 0. This translates into the condition Θf e b/q. Inserting this condition into eq 25 then gives the upper bound condition

b2 ) q2 1 +

Figure 10. Effect of thermal conductivity on the percolation threshold for different values of the ratio ∆Ta/∆Tr for the case of finite kinetics (Th ) 1).

(28)

Given that b is small, the above equation predicts an upper bound well beyond the range shown in Figures 8 and 9. Effect of Thermal Conductivity. It would be of interest to determine whether the present gasless model also leads to the same sensitivity to the thermal conductivity as demonstrated in the full FC model in Figure 4. The results shown in Figure 10 display a behavior qualitatively similar to that of the full FC case. A nonmonotonic trend is exhibited, and the effect is more pronounced at higher values of the ratio ∆Ta/∆Tr ) γ/q. We can interpret this behavior using the previous findings. It was found that, in the nonadiabatic case, the two dimensionless groups Th and Nu, both of which contain the thermal conductivity, have inverse effects

It was noted in the introduction of this paper, that, in some previous works, the problem of propagating combustion fronts was treated as a phase transition with percolation properties. Scaling laws were proposed to examine the behavior of the system, and in some cases, agreement with ordinary percolation theory was noted. Our analysis, however, shows that the behavior of such systems is much more complex. Specifically, we found that the competition among heat transfer, heat losses, and reaction kinetics plays an important role and can lead to extinction, i.e., to “nonpercolating” fronts, even in systems that are fully occupied with fuel. Such lack of percolation is due not to the lack of connectivity of the underlying fuel sites, as typical percolation processes assume, but rather to the intrinsic nonlinearity of the Arrhenius kinetics. In addition, in a disordered lattice, heat transfer through the conducting solid has a long-ranged effect that can lead to the ignition of distant sites, even though they might be disconnected, in the context of fuel occupancy. The result is a wide variability of the percolation threshold (or perhaps more appropriately, the extinction threshold), as a function of the various kinetic, thermodynamic, and heat-transfer variables. To demonstrate such a long-range effect, we consider consecutive snapshots at the same occupation probability for a gasless nonadiabatic process (Figure 11). In the first snapshot, all nearest neighbors to the front do not contain active fuel. The second snapshot shows how the front then advances by igniting a more distant third neighbor, after which the front spreads rapidly (third snapshot in Figure 11). Thus, there exist sites that become activated and burned, even though they are not nearest neighbors of previously burned sites and could not have been reached by a percolation process. Rather, there is a jump across sites not containing fuel. In addition, the corresponding cluster of reacted sites has no apparent similarity to the fractal percolation cluster. Figure 12 shows three patterns corresponding


Figure 11. Consecutive snapshots at the same occupation probability for a gasless nonadiabatic process. In the first snapshot (A), all nearest neighbors to the front do not contain active fuel. The second snapshot (B) shows how the front then advances by igniting a more distant third neighbor (indicated by the arrows), after which the front spreads rapidly [shaded area in third snapshot (C)].

to three different values of the occupation probability with Nu ) 0.05 and Th ) 1, under conditions such that pc ) 0.242. The lack of fractal properties (which roughly translates to the absence of “holes” of all sizes) in the patterns in the second and third panels suggests that a true percolation process did not develop, at least in this example. The underlying reason is the long-range reach of heat conduction, which allows for the possibility of ignition of even distant sites and “homogenizes” the problem. The time evolution of the completely reacted sites for the parameters of Figure 12 is shown in Figure 13 for the three occupation fractions. The simulations demonstrate an almost linear dependence of the total number of burned sites on rescaled time (Figure 13). This is in contrast to the results of Albinet et al.13 or Beer,15 where the dependence on time is a power law with exponents 0.78 and 1.2 for a small array or 1.9 for a large array. Our results show instead a classical behavior with an exponent in the time dependence close to 1 (Figure 13). As we noted above, the reason for this difference is the fact that true percolation behavior does

Figure 12. Patterns of front propagation for a nonadiabatic, gasless FC with finite kinetics and pc ) 0.242, corresponding to p ) 0.22 (case A), p ) 0.25 (case B), and p ) 0.35 (case C). Case A leads to extinction. It is clear that, despite the fact that the system is far below the percolation threshold of ordinary site percolation (equal to 0.59 for a square lattice), front propagation occurs, as a result of the long-ranged effects of heat conduction. The cluster of burned sites indicated by * does not have the fractal properties of a percolation cluster.

not actually become established in the problem, as the heat conduction acts to regularize the problem. Before closing, we reiterate that, in this paper, convection was not included, except for the full simulation case shown in Figure 4. The rest of the simulations did not account for convection. The reason for this simplification was to accelerate the calculations. Nonetheless, it is noted that, for the typical cases considered, the qualitative features of the results obtained are very similar in the corresponding cases (e.g., compare Figures 4 and 10). 5. Conclusions In this work, we used a pore-network model to study the effect of the randomness in the solid fuel distribution


with fully occupied fuel sites. In either case, the availability of the fuel and its distribution in space were found to be important parameters in the process. Acknowledgment This paper was written to honor Professor George R. Gavalas. The senior author of the work (Y.C.Y.) owes a great deal of gratitude to Professor Gavalas for his mentorship, guidance, encouragement, and friendship over the many years of close association with him. Professor Gavalas has left an indelible mark on the author’s career, for which he is deeply grateful. Literature Cited

Figure 13. Evolution of the number of burned sites as a function of time in nonadiabatic gasless FC with finite kinetics (NuL ) 0.05, Th ) 1) for three occupation fractions, corresponding to those in Figure 12. The exponent refers to the time evolution of the number of burned sites, assumed as a power law of time.

on the propagation of filtration combustion fronts in porous media. Both a complex pore-network model for full filtration combustion and a more simplified gasless combustion model were used. When the solid fuel is not continuously distributed, a critical threshold exists in the occupancy fraction of the solid fuel, below which sustained combustion cannot occur in a nonadiabatic system with finite kinetics, or in a general system with a specified ignition temperature. We find that the minimum threshold (akin to a “percolation threshold”) for the fuel distribution is a function of the kinetic and heat-transfer parameters and can vary over a large range. The detailed behavior is controlled by the balance between heat transfer and reaction kinetics, which endows the system with a behavior more complex than simple percolation. Percolation theory and the associated scaling laws were earlier suggested as a model for such systems. Arguments on whether site percolation theory can explain this critical behavior abound in the literature. However, experimental work does not support this contention. Our simulations also show that it is not the case in systems in which heat transfer is important. The long-range effect of heat conduction acts to move the behavior of the system away from true percolation behavior. Observations of the propagation of the fronts confirmed the existence of long-ranged ignition, which is the main cause for the departure from site percolation. Conversely, the effect of heat losses is sufficient to lead to extinction (thus, to “nonpercolation”) even for systems

(1) Prats, M. Thermal Recovery; SPE Monograph Series 7; SPE (Society of Petroleum Engineers): Richardson, TX, 1982. (2) Aldushin, A. P.; Seplyarsky, B. S. Propagation of exothermic reaction in a porous-medium during gas blowing. Dokl. Akad. Nauk SSSR 1978, 241 (1), 72. (3) Aldushin, A. P. New results in the theory of filtration combustion. Combust. Flame 1993, 94, 308. (4) Akkutlu, Y. I.; Yortsos, Y. C. The dynamics of In-situ Combustion Fronts in Porous Media. Combust. Flame 2003, 134, 229. (5) Aldushin, A. P.; Matkowsky, B. J. Instabilities, fingering and Saffman-Taylor problem in filtration combustion. Combust. Sci. Technol. 1998, 133, 293. (6) Zik, O.; Moses, E. Fingering instability in combustion: An extended view. Phys. Rev. E 1999, 60, 518. (7) Conti, M.; Marconi, U. M. B. Fingering in slow combustion. Physica A 2002, 312, 381. (8) Lu, C. The dynamics of filtration combustion at the porenetwork scale. Ph.D. Dissertation, University of Southern California, Los Angeles, CA, 2003. (9) Lu, C.; Yortsos, Y. C. The dynamics of forward filtration combustion at the pore-network scale. AIChE J., manuscript submitted. (10) Lu, C.; Yortsos, Y. C. The dynamics of reverse filtration combustion at the pore-network scale. Phys. Rev. E, manuscript submitted. (11) Stauffer, D. Introduction to Percolation Theory; Taylor & Francis: London, 1992. (12) Morvan, D.; Dupuy, D. L. Modeling of fire spread through a forest fuel bed using a multiphase formulation. Combust. Flame 2001, 127, 1981. (13) Albinet, G.; Searby, G.; Stauffer, D. Fire propagation in a 2-D random medium. J. Phys. (Paris) 1986, 47, 1. (14) Tephany, H.; Nahmias J.; Duarte J. A. M. S. Combustion in heterogeneous media: A critical phenomenon. Physica A 1997, 242, 57. (15) Beer, T. Percolation theory and fire spread. Combust. Sci. Technol. 1990, 72, 297. (16) Satik, C.; Yortsos, Y. C. A pore-network study of bubble growth in porous media driven by heat transfer. J. Heat Transfer 1996, 118, 455. (17) Hwang, S.; Mukasyan, A. S.; Varma, A. Mechanisms of combustion wave propagation in heterogeneous reaction systems. Combust. Flame 1998, 115, 354.

Received for review August 4, 2003 Revised manuscript received December 15, 2003 Accepted December 17, 2003 IE0306372

Percolation Phenomena in Filtration Combustion - Industrial

Recommend Documents