# Self-Propagating Reactive Fronts in Compacts of Multilayered Particles

Authors: Ihab Sraj; Manav Vohra; Leen Alawieh; Timothy P. Weihs; Omar M. Knio
Journal: Journal of Nanomaterials
 (2013)

---

## Abstract

Reactive multilayered foils in the form of thin films have gained interest in various applications such as joining, welding, and ignition. Typically, thin film multilayers support self-propagating reaction fronts with speeds ranging from 1 to 20 m/s. In some applications, however, reaction fronts with much smaller velocities are required. This recently motivated Fritz et al. (2011) to fabricate compacts of regular sized/shaped multilayered particles and demonstrate self-sustained reaction fronts having much smaller velocities than thin films with similar layering. In this work, we develop a simplified numerical model to simulate the self-propagation of reactive fronts in an idealized compact, comprising identical Ni/Al multilayered particles in thermal contact. The evolution of the reaction in the compact is simulated using a two-dimensional transient model, based on a reduced description of mixing, heat release, and thermal transport. Computed results reveal that an advancing reaction front can be substantially delayed as it crosses from one particle to a neighboring particle, which results in a reduced mean propagation velocity. A quantitative analysis is thus conducted on the dependence of these phenomena on the contact area between the particles, the thermal contact resistance, and the arrangement of the multilayered particles.


---

## Body

## 1. Introduction

Reactive multilayered materials have recently gained increasing interest in various applications, including joining, brazing, sealing, and ignition of secondary reactions [1–13]. Typically, these materials are fabricated in the form of multilayered foils [13–20] using vapor deposition techniques [13, 14, 21, 22], ball-milling [23–25], and rolling [26]. The former result in films with precisely engineering microstructures, with layering ranging from nanometers to microns. The layers alternate between materials that mix exothermically; their structure is characterized in terms of the bilayer period, 2(1+γ)δ, defined as the thickness of a pair of chemically distinct layers, where δ is the Al layer half-thickness and γ is the ratio of Ni layer thickness to Al layer thickness (Figure 1). For instance, Ni/Al multilayered foils, which fall within the scope of the present study, are routinely made with δ in the range of 10–250 nm. For such bilayer values, sustained reactions can be initiated via a localized stimulus (e.g., electric or laser discharge), resulting in self-propagating fronts travelling at speeds ranging from 1 to 20 m/s [13, 21, 27]Figure 1
Schematic illustration of an idealized compact of five multilayered particles. Note that each particle contains multiple bilayers, but only one is illustrated. The individual layers are assumed to be separated by a premixed region, labeled NiAl, of thickness,w. The particles are identical of fixed length L=1 mm. The contact area, Ac, between the particles and thickness d are the same across the compact, but different values are considered in the computations. The flame is initiated at the leftmost particle, with a thermal “spark” of temperature Ts and width ws. The grid is schematically shown on the second and third particles. Similar gridding is applied to all particles but is not shown.One of the advantages of multilayered foils concerns their controlled microstructure, which enables a reaction front with a well-defined velocity, which in many cases can be made insensitive to environmental conditions or boundary conditions imposed by the application [28]. One drawback, however, concerns the fact that, with multilayered foils, it is difficult to achieve sustained reaction fronts with slow propagation velocity, at least without the risk of quenching the reaction. Slow moving fronts may be required for certain applications, including chemical time delays, or in situations where the heating provided by the reactions must be sustained over extended durations. To overcome this drawback, Fritz et al. [29] investigated the self-propagation of exothermic formation reactions within loose compacts of Ni/Al multilayered particles. The particles were fabricated by DC magnetron sputtering onto nylon mesh substrates. The sputtered multilayered coating was broken into particles that matched the size of the mesh elements by bending the mesh under water. The loose particles were then collected and loosely packed into a glass tube. This fabrication method resulted in compacts supporting a self-propagation velocity, that is, substantially smaller than in foils with similar multilayering, and that can be controlled by varying the packing density.The present study aims at ultimately developing computational models that can predict the behavior of reaction fronts in compacts of multilayered particles and characterize their dependence on the particle distribution and on the layering within individual particles. To this end, in this paper we consider an idealized two-dimensional model of a compact comprising rectangular multilayered particles in thermal contact. The model is used to examine the evolution of self-propagating reactions, particularly as the front traverses neighboring particles. The computations are then used to analyze the dependence of the front velocity on the contact area, the thermal contact resistance, and the number of particles within the idealized compact.


## 2. Formulation

### 2.1. Idealized Particle Compact

We will focus primarily on a simplified two-dimensional compact comprising rectangularly shaped, multilayered particles that are in thermal contact. A schematic of a typical configuration is shown in Figure1. Note that each particle contains multiple bilayers, though only one is shown in the schematic. The particles are assumed to be identically shaped, namely, with width L and thickness d, and to have the same internal microstructure. The latter consists of uniform layers alternating between chemically distinct reactant layers that are separated by a thin premixed product region; within each bilayer, a 1 : 1 molar ratio of Al and Ni is assumed. The contact area, Ac, between the neighboring particles is also assumed to be uniform across the compact. This contact area allows the conduction of heat, and if possible the transfer of the reactive front from one particle to the next. In the schematic of Figure 1, five particles are shown for illustration purposes, though a smaller or larger number can be considered.


### 2.2. Reaction Model

The evolution of the reaction within the particle compact is analyzed using the reduced model developed in [30–32]. This model, originally developed for a single multilayered solid, is adapted to the present setup by treating individual particles as separate regions, while accounting for thermal conduction through the contact areas.Within individual particles, the reaction is described in terms of a coarse-grained continuum model, coupling the conservation of energy equation to a mixture evolution equation. Dividing the particle into a finite number of regions or cells, conservation of energy within each cell is expressed in terms of the region-averaged enthalpy equation:(1)∂H∂t=-1|Ω|∫Ω‍(∇·q)dΩ+∂Q¯∂t,
where H is the region-averaged enthalpy, |Ω| is the area of the region,
(2)∂Q¯∂t=-ρcp¯ΔTf∂C2¯(t)∂t
is the region-averaged heat release term, ρcp¯ is the mean specific heat, ΔTf≡ΔHrx/ρcp¯ is the nominal flame temperature, and ΔHrx is the (negative) heat of mixing. Note that we have assumed that Q exhibits a quadratic dependence on the concentration, C, according to [33, 34]
(3)Q¯=-ρcp¯ΔTfC2¯,
where C is a dimensionless concentration defined so that |C|=1 in the reactants, and C=0 for the product. Also note that the temperature, T, can be recovered from the enthalpy H, following the methodology outlined in [34]; namely,
(4)T={T0+Hρcp¯ifH<H1TmAlifH1<H<H2TmAl+(H-H2)ρcp¯ifH2<H<H3TmNiifH3<H<H4TmNi+(H-H4)ρcp¯ifH4<H<H5TmNiAlifH5<H<H6TmNiAl+(H-H6)ρcp¯ifH6<H,
where TmAl=933 K, TmNi=1728 K, and TmNiAl=1912 K denote the melting temperatures of Al, Ni, and NiAl, respectively, ℋfAl, ℋfNi, and ℋfNiAl are the corresponding heats of fusion (per unit mole), β≡C¯/(1+γ) represents the fraction of pure (unmixed) Al, ΔHfAl≡ρAlℋfAl/MAl, ΔHfNi≡ρNiℋfNi/MNi, ΔHfNiAl≡ρ¯ℋfNiAl/MNiAl, H1≡ρcp¯(TmAl-T0), H2≡H1+βΔHfAl, H3≡H2+ρcp¯(TmNi-TmAl), H4≡H3+βγΔHfNi, H5≡H4+ρcp¯(TmNiAl-TmNi), H6≡H5+(1-C¯)ΔHfNiAl,
(5)γ≡ρAlρNiMNiMAl,ρ¯≡(ρAl+γρNi)(1+γ),ρcp¯≡(ρAlcpAl+γρNicpNi)(1+γ)ρAl and ρNi are the densities of Al and Ni, and cpAl and cpNi are the corresponding heat capacities, whereas MAl, MNi, and MNiAl denote the atomic weights of Al, Ni, and NiAl, respectively.The evolution of the mean concentration and its second moment,C¯ and C2¯, are obtained through the reduced formalism introduced in [30]. Specifically, canonical solutions are used to relate the evolution of these moments to that of a stretched time variable τ that is governed by
(6)∂τ∂t=D(T)δ2,
where D denotes atomic diffusivity and δ is the half-thickness of an Al layer. The diffusivity is assumed to follow an Arrhenius dependence on T, according to
(7)D(T)=D0exp(-ERT),
where D0=2.18×10-6 m2/s is the preexponent, E=137 kJ/mol is the activation energy, and R is the universal gas constant. The values of D0 and E are obtained as best fits based on experimental measurements of the velocity of axially propagating fronts [35].In the computations below, we assume the local heat fluxq which is given by Fourier’s law:
(8)q=-k∇T,
where k is the thermal conductivity. Generally, k may depend on the composition and temperature, and possibly on microstructure as well [36]. However, since our primary focus in the present study has been on the compact structure, we have restricted our attention to a uniform conductivity model:
(9)k=k¯≡kAl+γkNi1+γ,
where kAl and kNi are the thermal conductivities of Al and Ni, respectively. We finally note that, apart from interparticle thermal transport, heat losses to the environment have also been ignored.


### 2.3. Simulation

The coupled system, (1) and (6), is solved numerically using the finite difference scheme adapted from [31]. In the present implementation, each layered particle is considered as a 2D domain that is discretized with a uniform Cartesian grid of mesh size Δx and Δy along the x and y directions, respectively; see Figure 1. The two field variables, H and τ, are discretized at cell centers, and conservative second-order differences are used to compute fluxes. All other local physical quantities can be readily obtained based on these two field variables. Interparticle mass diffusion is ignored, and adiabatic conditions are imposed on all domain boundaries, expect for contact areas where continuity of temperature and heat flux is imposed. These relationships may be readily generalized by incorporating a simplified model for thermal contact resistance. Specifically, with the present second-order centered difference discretization, the impact of thermal resistance is accounted for by modifying the value of the thermal conductivity at the contact surfaces, namely; according to
(10)kmod=k1+k×Rc/Δy,
where kmod is the augmented value of the thermal conductivity at the solid-solid boundaries, and Rc is the magnitude of the thermal contact resistance.


## 2.1. Idealized Particle Compact

We will focus primarily on a simplified two-dimensional compact comprising rectangularly shaped, multilayered particles that are in thermal contact. A schematic of a typical configuration is shown in Figure1. Note that each particle contains multiple bilayers, though only one is shown in the schematic. The particles are assumed to be identically shaped, namely, with width L and thickness d, and to have the same internal microstructure. The latter consists of uniform layers alternating between chemically distinct reactant layers that are separated by a thin premixed product region; within each bilayer, a 1 : 1 molar ratio of Al and Ni is assumed. The contact area, Ac, between the neighboring particles is also assumed to be uniform across the compact. This contact area allows the conduction of heat, and if possible the transfer of the reactive front from one particle to the next. In the schematic of Figure 1, five particles are shown for illustration purposes, though a smaller or larger number can be considered.


## 2.2. Reaction Model

The evolution of the reaction within the particle compact is analyzed using the reduced model developed in [30–32]. This model, originally developed for a single multilayered solid, is adapted to the present setup by treating individual particles as separate regions, while accounting for thermal conduction through the contact areas.Within individual particles, the reaction is described in terms of a coarse-grained continuum model, coupling the conservation of energy equation to a mixture evolution equation. Dividing the particle into a finite number of regions or cells, conservation of energy within each cell is expressed in terms of the region-averaged enthalpy equation:(1)∂H∂t=-1|Ω|∫Ω‍(∇·q)dΩ+∂Q¯∂t,
where H is the region-averaged enthalpy, |Ω| is the area of the region,
(2)∂Q¯∂t=-ρcp¯ΔTf∂C2¯(t)∂t
is the region-averaged heat release term, ρcp¯ is the mean specific heat, ΔTf≡ΔHrx/ρcp¯ is the nominal flame temperature, and ΔHrx is the (negative) heat of mixing. Note that we have assumed that Q exhibits a quadratic dependence on the concentration, C, according to [33, 34]
(3)Q¯=-ρcp¯ΔTfC2¯,
where C is a dimensionless concentration defined so that |C|=1 in the reactants, and C=0 for the product. Also note that the temperature, T, can be recovered from the enthalpy H, following the methodology outlined in [34]; namely,
(4)T={T0+Hρcp¯ifH<H1TmAlifH1<H<H2TmAl+(H-H2)ρcp¯ifH2<H<H3TmNiifH3<H<H4TmNi+(H-H4)ρcp¯ifH4<H<H5TmNiAlifH5<H<H6TmNiAl+(H-H6)ρcp¯ifH6<H,
where TmAl=933 K, TmNi=1728 K, and TmNiAl=1912 K denote the melting temperatures of Al, Ni, and NiAl, respectively, ℋfAl, ℋfNi, and ℋfNiAl are the corresponding heats of fusion (per unit mole), β≡C¯/(1+γ) represents the fraction of pure (unmixed) Al, ΔHfAl≡ρAlℋfAl/MAl, ΔHfNi≡ρNiℋfNi/MNi, ΔHfNiAl≡ρ¯ℋfNiAl/MNiAl, H1≡ρcp¯(TmAl-T0), H2≡H1+βΔHfAl, H3≡H2+ρcp¯(TmNi-TmAl), H4≡H3+βγΔHfNi, H5≡H4+ρcp¯(TmNiAl-TmNi), H6≡H5+(1-C¯)ΔHfNiAl,
(5)γ≡ρAlρNiMNiMAl,ρ¯≡(ρAl+γρNi)(1+γ),ρcp¯≡(ρAlcpAl+γρNicpNi)(1+γ)ρAl and ρNi are the densities of Al and Ni, and cpAl and cpNi are the corresponding heat capacities, whereas MAl, MNi, and MNiAl denote the atomic weights of Al, Ni, and NiAl, respectively.The evolution of the mean concentration and its second moment,C¯ and C2¯, are obtained through the reduced formalism introduced in [30]. Specifically, canonical solutions are used to relate the evolution of these moments to that of a stretched time variable τ that is governed by
(6)∂τ∂t=D(T)δ2,
where D denotes atomic diffusivity and δ is the half-thickness of an Al layer. The diffusivity is assumed to follow an Arrhenius dependence on T, according to
(7)D(T)=D0exp(-ERT),
where D0=2.18×10-6 m2/s is the preexponent, E=137 kJ/mol is the activation energy, and R is the universal gas constant. The values of D0 and E are obtained as best fits based on experimental measurements of the velocity of axially propagating fronts [35].In the computations below, we assume the local heat fluxq which is given by Fourier’s law:
(8)q=-k∇T,
where k is the thermal conductivity. Generally, k may depend on the composition and temperature, and possibly on microstructure as well [36]. However, since our primary focus in the present study has been on the compact structure, we have restricted our attention to a uniform conductivity model:
(9)k=k¯≡kAl+γkNi1+γ,
where kAl and kNi are the thermal conductivities of Al and Ni, respectively. We finally note that, apart from interparticle thermal transport, heat losses to the environment have also been ignored.


## 2.3. Simulation

The coupled system, (1) and (6), is solved numerically using the finite difference scheme adapted from [31]. In the present implementation, each layered particle is considered as a 2D domain that is discretized with a uniform Cartesian grid of mesh size Δx and Δy along the x and y directions, respectively; see Figure 1. The two field variables, H and τ, are discretized at cell centers, and conservative second-order differences are used to compute fluxes. All other local physical quantities can be readily obtained based on these two field variables. Interparticle mass diffusion is ignored, and adiabatic conditions are imposed on all domain boundaries, expect for contact areas where continuity of temperature and heat flux is imposed. These relationships may be readily generalized by incorporating a simplified model for thermal contact resistance. Specifically, with the present second-order centered difference discretization, the impact of thermal resistance is accounted for by modifying the value of the thermal conductivity at the contact surfaces, namely; according to
(10)kmod=k1+k×Rc/Δy,
where kmod is the augmented value of the thermal conductivity at the solid-solid boundaries, and Rc is the magnitude of the thermal contact resistance.


## 3. Results

Simulations were conducted to characterize the evolution of reactive fronts within the idealized compacts and to examine their dependence on the properties of the compact. To this end, the particle length is held fixed,L=1 mm. The microstructure of the particles is also held fixed, with the Al layer half-thickness δ=250 nm and the premix width w=0.8 nm. Attention is thus focused on the effects of contact area, Ac, and thermal contact resistance, Rc.The self-propagating reactions are initiated using a thermal spark of temperatureTs=1400 K and width ws=200 μm. We start by characterizing how the self-propagation front crosses from one particle to another and then characterize the dependence of the reaction front speed on the properties of the compact. As in [32], the 2D front position is tracked through the first moment of the heat release rate:
(11)xf=∫Γ‍x(∂Q¯/∂t)dΩ∫Γ‍(∂Q¯/∂t)dΩ,
where Γ denotes the entire computational domain. For the compacts considered in the present study, the reaction front remains predominantly vertical, and thus we focus exclusively on its x coordinate. Differentiating the latter with time then yields its instantaneous velocity. Unless otherwise noted, perfect thermal contact conditions are assumed.


### 3.1. Front Propagation in Particle Compact

Figure2 shows the flame position (x-direction) versus time for a compact comprising five layered particles of thickness 10μm and contact area Ac=200μm. Once established, the front initially propagates with uniform velocity in the first particle, t≤2ms. As it enters the contact region between the first and second particle, its rate of advancement is significantly delayed. The flame transfers to the second particle at about t=12 ms (at xf=0.9 mm); this corresponds to about half the total reaction time. At the second contact position, a similar delay occurs, though its duration is considerably smaller, approximately 1.6 ms. For the subsequent crossovers, the time delay remains nearly the same, approximately 4 ms. It is thus seen that a substantial reduction in the front propagation speed occurs due to the delay in the crossover of the front from one particle to the next. The flame retardation is a manifestation of the nonpremixed structure of the system, and of the fact that the thermal conductivity is several orders of magnitude larger than atomic diffusivity, even at high temperature. This results in a local quenching of the front, as further described below.Figure 2
Flame front position versus time for contact areaAc=200 μm with perfect thermal contact. The dashed lines indicate the region at the contact area between successive particles.To further investigate the mechanism of front crossover, we plot in Figure3 distributions of temperature and heat release rate at selected time instants. Attention is focused on the front crossover between the first and second layered particles. At t=1 ms (Figure 3(a)), the reactive front is clearly located in the middle of the first particle. At this instant, the maximum temperature in the compact corresponds to the flame adiabatic temperature, Tmax=1912 K, and the width of the reaction zone is around 100μm. At t=5ms, the front enters the first contact region, as shown in Figure 3(b). A noticeable rise in the temperature of the second particle can also be observed. As heat is lost from the first particle to the second, a slight reduction in the peak temperature occurs, Tmax=1846 K. This is accompanied by a dramatic decrease in the peak instantaneous heat release rate, from about 7 MW/m3 att=1 ms to about 0.05 MW/m3 at t=5 ms, and a substantial widening of the reaction zone width. Figure 3(c) shows that, at t=12 ms, the reactive front is still located within the contact area. The peak temperature has dropped to Tmax=1296 K. This phenomenon is dependent on the size of the contact area, as illustrated in Figure 4. There appears to be little change in the width of the reaction zone between t=5 ms and t=12 ms, though the peak heat release rate has increased to 0.14 MW/m3. At t=13.6 ms, the flame front has completely transferred to the second particle, as shown in Figure 3(d). At this time, the maximum temperature (1912 K) and heat release rate (15 MW/m3) once again reach values associated with a steadily propagating front.Instantaneous distributions of temperature in K (left) and local heat release rate in MW/m3 (right) at selected times ((a)–(d)). The contact area Ac=200 μm and perfect thermal contact is assumed. Distributions are shown on the front crossover between the first and second particles.
(a)
(b)
(c)
(d)Figure 4
Evolution of the maximum temperature within the compact for different values ofAc, as indicated. The results are obtained assuming perfect thermal contact.It appears that, for the present configuration, the flame crossover from one particle to the next can be characterized by a front propagation phase that essentially leads to the consumption of the reactants in the first particle, followed by an extended heating phase, and followed by a rapid ignition and the formation of a self-propagating front in the neighboring particle. This succession of regimes is further illustrated in Figure5(a), which clearly shows that the heat release rate drops to very small values during the heating phase (5ms≤t≤12ms). Another perspective is provided from section-averaged temperature profiles plotted in Figure 5(b). It is interesting to note that at t=13.6 ms, that is, following the formation of a front in the second particle, heat is being conducted away from the reaction zone in two directions, towards the first particle and towards the third. This phenomenon can be readily appreciated from the nonmonotonic behavior of the corresponding temperature profile; it is due to the sudden rise in heat release rate and the substantial heat loss from the first particle to the second during the heating phase.Instantaneous (a) local heat release rate and (b) temperature profiles at selected time instants forAc=200 μm with perfect thermal contact. The x-axis spans the first two particles in the compact. The overlap region between the two particles is shown as a dashed curve.
(a)
(b)The results obtained for the present configuration indicate that the motion of self-propagating fronts in particle compacts can involve a complex succession of phenomena, and large variations in the flame temperature and in heat release rate. They also point to the possibility of substantial reduction in the mean front propagation velocity due to heating delay. Below, we focus on quantifying this effect and characterizing its dependence on properties of the compact.


### 3.2. Flame Propagation Speed

The average speed of reaction fronts propagating in Ni/Al laminates of uniform thickness has been analyzed in numerous published works [27, 28, 33, 34, 37–39]. The reduced model used in this work has been successful in capturing the dependence of the average front velocity, Vavg, on bilayer and premix parameters [30]. These studies have shown that for steadily propagating fronts, Vavg can be simply determined by differentiating the flame position versus time curve. On the other hand, when the front propagates in an unsteady manner, a sufficiently long time period must be considered, so as to appropriately filter out the impact of velocity oscillations.In the present work, we extend the concept used for uniform thickness foils to the case of particle compacts as follows. Regardless of the details of the configuration, we form an estimate ofVavg by taking the ratio of the total length of the compact, Lt=nL-(n-1)Ac, by the total reaction time, tr, defined as the time needed for reaction front to reach the end of the last particle. Here, n is the total number of particles in the compact.As shown in the previous section, the delay in front crossover may vary appreciably from one contact region to another. To ensure an appropriate estimate of the average velocity, we analyzed the dependence ofVavg on the number of particles present in the compact, namely, by systematically increasing n in the range 1≤n≤21. The results (not shown) indicated that as n becomes large, Vavg tends to an asymptotic value. Specifically, when n≥5, Vavg falls within 5% of each other, in the entire range of Ac that is considered. Consequently, we have used the estimates of Vavg for n=5 in all the computations presented below.


#### 3.2.1. Effect of Contact Area

Simulations were conducted by systematically varyingAc in the range 10–200 μm; in all cases, the particle thickness was held fixed at d=10 μm. Note that for large values of Ac, with the current setup and ignition properties, the flame front quenches as it enters the contact region between the first and second particles. This is illustrated in Figure 4, which shows that for Ac=250 μm, Tmax drops to a low value, and the reaction is not initiated in the second particle. This observation is at the origin of the restricted range of Ac.Figure6 shows the dependence of Vavg on Ac for the case of perfect thermal contact. It is seen that Vavg monotonically decreases as Ac increases. For Ac=10 μm, the computed average velocity, Vavg=0.42 m/s, coincides with the flame speed of a continuous multilayered film of the same thickness d=10 μm. Thus, for this small value of Ac, the ignition delay associated with flame crossover is rather insignificant. As shown in Figure 4, the drop in the maximum temperature of the first particle is small for Ac≤50 μm. For Ac=100 μm and 200 μm, the impact of the ignition delay is significant, with predicted values Vavg=0.22 m/s and 0.15 m/s, respectively. These correspond to about a twofold or threefold decrease in Vavg compared to continuous multilayered film with uniform thickness.Figure 6
Profiles ofVavg versus Ac. The results are obtained assuming perfect thermal contact. The average speed for a continuous film is indicated using a dashed line for reference.To further illustrate the impact of the heating time,th, we plot in Figure 7 the flame position versus time; curves are generated for different values of Ac. The results indicate that during the front propagation phase, the slopes of the individual curves are essentially identical. Thus, for small values of Ac, the reaction time, that is, the time needed for the front to traverse an individual particle, is independent of the contact area. Accordingly, our calculations suggest that the average flame velocity in a compact of particles is close to that of a continuous film when the contact area is of the same order of magnitude as the particle thickness. The results also reveal that the substantial reduction in Vavg that occurs as Ac increases is mainly due to an increase in the particle heating time. It is also interesting to note that for the same configuration, the heating time may vary appreciably from one flame crossover event to another. For small values of Ac, the ignition delay essentially stabilizes after the first crossover. For Ac=200 μm, however, it takes more than a single passage for the ignition delay to stabilize. In particular, a large variation is observed between the first and the second flame crossings. The pronounced delay associated with the first crossing for Ac=200 μm is affected by the substantial heat loss that occurs immediately following ignition, due to the proximity between the ignition spot and the contact between the first and second particles.Figure 7
Flame front position versus time for different values ofAc, as indicated. The results are obtained assuming perfect thermal contact.The analysis above was repeated for a particle thicknessd=50 μm. The results (not shown) revealed similar trends for the dependence of Vavg on Ac. Small quantitative differences were however observed, as further discussed below.


#### 3.2.2. Effect of Particle Thickness

In this section, we briefly examine the impact of the particle thickness,d, on the average flame velocity. We also take advantage of the present analysis to make a connection between the single-row arrangement shown in Figure 1 and a periodic configuration depicted in Figure 8. In the latter, the stack used in Figure 1 is repeated multiple times in the vertical direction. Assuming that the leftmost particles are ignited identically and simultaneously, this is modelled as a periodic arrangement as indicated in Figure 8. Specifically, periodic boundary conditions are used at the horizontal boundaries of the domain, indicated in Figure 8 using dashed lines. Based on symmetry arguments, predictions obtained using the periodic configuration should correspond to those obtained using the single-row model, provided that the same physical parameters are used except for the particle thickness, whose value in the periodic arrangement is twice that of the single-row model. As briefly outlined below, this is in fact observed in the results of the computations.Figure 8
Schematic illustration of particles compact with periodic arrangement. Though the individual particles comprise multiple bilayers, only one is shown in the schematic. Periodicity is simulated by restricting the domain to the region lying between the dotted lines, and imposing periodic boundary conditions on these lines. Other boundary and initial conditions are similar to those depicted in Figure1.Figure9 shows curves of flame position versus time for different values of d. The results are obtained using the single-row particle model. Also shown is a curve obtained using the periodic configuration model with d=50 μm. The results obtained using the single-row model indicate that in the range of particle thicknesses considered, the average flame velocity generally decreases as d increases. However, the variations of Vavg with d are generally small, about 5% or less. Thus, in the present setup (where heat losses are ignored), the average flame speed appears to be less affected by particle thickness than by particle contact area.Figure 9
Flame front position versus time for the single-row configuration, with particle thickness in the ranged= 5–50 μm, as indicated. The symbols depict results obtained for a periodic configuration with d=50 μm.Also shown on Figure9 are the results obtained with the periodic stack model using d=50 μm. Note that the curve, obtained with the periodic model with d=50 μm, lies between d=20 μm and d=30 μm curves obtained using the single-row model. Thus, the periodic arrangement with d=50 μm yields similar results to the single-row model with d~25 μm, which illustrates our earlier claim regarding the correspondence between the two arrangements.


#### 3.2.3. Effect of Thermal Contact Resistance

The thermal contact resistance in a particle compact may depend on a variety of factors, including geometry, applied pressure, as well as surface cleanliness, roughness, and flatness [40–42]. Nominal reported values of Rc fall in a wide range, between 10-7 and 10-5 m2 K/W. In this section, we briefly investigate the effect of thermal contact resistance on the average velocity of the front, namely, by repeating the analysis for a compact of 10 μm thick particles. For brevity, we consider a single value Rc=10-5 m2 K/W and as before vary Ac in the range 10–200 μm.The predicted values ofVavg are plotted against Ac in Figure 10; also shown for comparison is the curve obtained for perfect thermal contact Rc=0. The results indicate that, unlike the case of perfect thermal contact, for Rc=10-5 m2 K/W, the Vavg as a function of Ac is no longer monotonic. Specifically, Vavg peaks for Ac=100 μm. For lower values of Ac, the velocity rapidly decreases as Ac is reduced, whereas for higher values, the rate of decay is substantially smaller.Figure 10
Average flame speed versus contact area. Curves are generated for the case of vanishing thermal contact resistance, and forRc=10-5 m2 K/W. Also shown for comparison are the experimental average flame speed measurements of Fritz et al. [29], for particle compacts with δ in the range 26–375 nm.To investigate the origin of this trend, we plot in Figure11 the front position versus time for the different values of Ac. The results indicate that the observed trends in Vavg are essentially governed by variation of the ignition delay. Specifically, for small contact areas, the presence of a large contact resistance leads to a significant increase in the heating time, which delays the initiation of the front in the adjacent particle. This phenomenon has also been observed through detailed examination of temperature and heating rate distributions (not shown). On the other hand, in the range Ac≥100 μm, the predictions for Rc=0 and Rc=10-5 m2 K/W become close to each other. Thus, in this range, the impact of the thermal contact resistance becomes substantially weaker.Figure 11
Flame front position versus time forRc=10-5 m2 K/W. Curves are generated for different values of Ac, as indicated.We finally compare our present predictions with experimental measurements of Fritz et al. [29], who reported average flame speeds of 0.0058–0.0260 m/s for low-density particle compacts (packing densities of 20%) with different half-layer thickness. The data from Fritz et al. are also plotted in Figure 10. One observes that the range of the measured flame speeds is comparable with the computed average flame speed for a contact area Ac=10 μm and thermal contact resistance Rc=10-5 m2 K/W. Although differences between the computational and experimental results can arise due to a variety of factors that are not accounted for in the model, such as heat losses, particle geometry, and compact packing density, the agreement between reported measurements with predictions obtained with high Rc and small contact area (low packing density) strongly suggests that thermal contact resistance may play a strong role in determining the self-propagation velocity in practical compacts.


## 3.1. Front Propagation in Particle Compact

Figure2 shows the flame position (x-direction) versus time for a compact comprising five layered particles of thickness 10μm and contact area Ac=200μm. Once established, the front initially propagates with uniform velocity in the first particle, t≤2ms. As it enters the contact region between the first and second particle, its rate of advancement is significantly delayed. The flame transfers to the second particle at about t=12 ms (at xf=0.9 mm); this corresponds to about half the total reaction time. At the second contact position, a similar delay occurs, though its duration is considerably smaller, approximately 1.6 ms. For the subsequent crossovers, the time delay remains nearly the same, approximately 4 ms. It is thus seen that a substantial reduction in the front propagation speed occurs due to the delay in the crossover of the front from one particle to the next. The flame retardation is a manifestation of the nonpremixed structure of the system, and of the fact that the thermal conductivity is several orders of magnitude larger than atomic diffusivity, even at high temperature. This results in a local quenching of the front, as further described below.Figure 2
Flame front position versus time for contact areaAc=200 μm with perfect thermal contact. The dashed lines indicate the region at the contact area between successive particles.To further investigate the mechanism of front crossover, we plot in Figure3 distributions of temperature and heat release rate at selected time instants. Attention is focused on the front crossover between the first and second layered particles. At t=1 ms (Figure 3(a)), the reactive front is clearly located in the middle of the first particle. At this instant, the maximum temperature in the compact corresponds to the flame adiabatic temperature, Tmax=1912 K, and the width of the reaction zone is around 100μm. At t=5ms, the front enters the first contact region, as shown in Figure 3(b). A noticeable rise in the temperature of the second particle can also be observed. As heat is lost from the first particle to the second, a slight reduction in the peak temperature occurs, Tmax=1846 K. This is accompanied by a dramatic decrease in the peak instantaneous heat release rate, from about 7 MW/m3 att=1 ms to about 0.05 MW/m3 at t=5 ms, and a substantial widening of the reaction zone width. Figure 3(c) shows that, at t=12 ms, the reactive front is still located within the contact area. The peak temperature has dropped to Tmax=1296 K. This phenomenon is dependent on the size of the contact area, as illustrated in Figure 4. There appears to be little change in the width of the reaction zone between t=5 ms and t=12 ms, though the peak heat release rate has increased to 0.14 MW/m3. At t=13.6 ms, the flame front has completely transferred to the second particle, as shown in Figure 3(d). At this time, the maximum temperature (1912 K) and heat release rate (15 MW/m3) once again reach values associated with a steadily propagating front.Instantaneous distributions of temperature in K (left) and local heat release rate in MW/m3 (right) at selected times ((a)–(d)). The contact area Ac=200 μm and perfect thermal contact is assumed. Distributions are shown on the front crossover between the first and second particles.
(a)
(b)
(c)
(d)Figure 4
Evolution of the maximum temperature within the compact for different values ofAc, as indicated. The results are obtained assuming perfect thermal contact.It appears that, for the present configuration, the flame crossover from one particle to the next can be characterized by a front propagation phase that essentially leads to the consumption of the reactants in the first particle, followed by an extended heating phase, and followed by a rapid ignition and the formation of a self-propagating front in the neighboring particle. This succession of regimes is further illustrated in Figure5(a), which clearly shows that the heat release rate drops to very small values during the heating phase (5ms≤t≤12ms). Another perspective is provided from section-averaged temperature profiles plotted in Figure 5(b). It is interesting to note that at t=13.6 ms, that is, following the formation of a front in the second particle, heat is being conducted away from the reaction zone in two directions, towards the first particle and towards the third. This phenomenon can be readily appreciated from the nonmonotonic behavior of the corresponding temperature profile; it is due to the sudden rise in heat release rate and the substantial heat loss from the first particle to the second during the heating phase.Instantaneous (a) local heat release rate and (b) temperature profiles at selected time instants forAc=200 μm with perfect thermal contact. The x-axis spans the first two particles in the compact. The overlap region between the two particles is shown as a dashed curve.
(a)
(b)The results obtained for the present configuration indicate that the motion of self-propagating fronts in particle compacts can involve a complex succession of phenomena, and large variations in the flame temperature and in heat release rate. They also point to the possibility of substantial reduction in the mean front propagation velocity due to heating delay. Below, we focus on quantifying this effect and characterizing its dependence on properties of the compact.


## 3.2. Flame Propagation Speed

The average speed of reaction fronts propagating in Ni/Al laminates of uniform thickness has been analyzed in numerous published works [27, 28, 33, 34, 37–39]. The reduced model used in this work has been successful in capturing the dependence of the average front velocity, Vavg, on bilayer and premix parameters [30]. These studies have shown that for steadily propagating fronts, Vavg can be simply determined by differentiating the flame position versus time curve. On the other hand, when the front propagates in an unsteady manner, a sufficiently long time period must be considered, so as to appropriately filter out the impact of velocity oscillations.In the present work, we extend the concept used for uniform thickness foils to the case of particle compacts as follows. Regardless of the details of the configuration, we form an estimate ofVavg by taking the ratio of the total length of the compact, Lt=nL-(n-1)Ac, by the total reaction time, tr, defined as the time needed for reaction front to reach the end of the last particle. Here, n is the total number of particles in the compact.As shown in the previous section, the delay in front crossover may vary appreciably from one contact region to another. To ensure an appropriate estimate of the average velocity, we analyzed the dependence ofVavg on the number of particles present in the compact, namely, by systematically increasing n in the range 1≤n≤21. The results (not shown) indicated that as n becomes large, Vavg tends to an asymptotic value. Specifically, when n≥5, Vavg falls within 5% of each other, in the entire range of Ac that is considered. Consequently, we have used the estimates of Vavg for n=5 in all the computations presented below.


### 3.2.1. Effect of Contact Area

Simulations were conducted by systematically varyingAc in the range 10–200 μm; in all cases, the particle thickness was held fixed at d=10 μm. Note that for large values of Ac, with the current setup and ignition properties, the flame front quenches as it enters the contact region between the first and second particles. This is illustrated in Figure 4, which shows that for Ac=250 μm, Tmax drops to a low value, and the reaction is not initiated in the second particle. This observation is at the origin of the restricted range of Ac.Figure6 shows the dependence of Vavg on Ac for the case of perfect thermal contact. It is seen that Vavg monotonically decreases as Ac increases. For Ac=10 μm, the computed average velocity, Vavg=0.42 m/s, coincides with the flame speed of a continuous multilayered film of the same thickness d=10 μm. Thus, for this small value of Ac, the ignition delay associated with flame crossover is rather insignificant. As shown in Figure 4, the drop in the maximum temperature of the first particle is small for Ac≤50 μm. For Ac=100 μm and 200 μm, the impact of the ignition delay is significant, with predicted values Vavg=0.22 m/s and 0.15 m/s, respectively. These correspond to about a twofold or threefold decrease in Vavg compared to continuous multilayered film with uniform thickness.Figure 6
Profiles ofVavg versus Ac. The results are obtained assuming perfect thermal contact. The average speed for a continuous film is indicated using a dashed line for reference.To further illustrate the impact of the heating time,th, we plot in Figure 7 the flame position versus time; curves are generated for different values of Ac. The results indicate that during the front propagation phase, the slopes of the individual curves are essentially identical. Thus, for small values of Ac, the reaction time, that is, the time needed for the front to traverse an individual particle, is independent of the contact area. Accordingly, our calculations suggest that the average flame velocity in a compact of particles is close to that of a continuous film when the contact area is of the same order of magnitude as the particle thickness. The results also reveal that the substantial reduction in Vavg that occurs as Ac increases is mainly due to an increase in the particle heating time. It is also interesting to note that for the same configuration, the heating time may vary appreciably from one flame crossover event to another. For small values of Ac, the ignition delay essentially stabilizes after the first crossover. For Ac=200 μm, however, it takes more than a single passage for the ignition delay to stabilize. In particular, a large variation is observed between the first and the second flame crossings. The pronounced delay associated with the first crossing for Ac=200 μm is affected by the substantial heat loss that occurs immediately following ignition, due to the proximity between the ignition spot and the contact between the first and second particles.Figure 7
Flame front position versus time for different values ofAc, as indicated. The results are obtained assuming perfect thermal contact.The analysis above was repeated for a particle thicknessd=50 μm. The results (not shown) revealed similar trends for the dependence of Vavg on Ac. Small quantitative differences were however observed, as further discussed below.


### 3.2.2. Effect of Particle Thickness

In this section, we briefly examine the impact of the particle thickness,d, on the average flame velocity. We also take advantage of the present analysis to make a connection between the single-row arrangement shown in Figure 1 and a periodic configuration depicted in Figure 8. In the latter, the stack used in Figure 1 is repeated multiple times in the vertical direction. Assuming that the leftmost particles are ignited identically and simultaneously, this is modelled as a periodic arrangement as indicated in Figure 8. Specifically, periodic boundary conditions are used at the horizontal boundaries of the domain, indicated in Figure 8 using dashed lines. Based on symmetry arguments, predictions obtained using the periodic configuration should correspond to those obtained using the single-row model, provided that the same physical parameters are used except for the particle thickness, whose value in the periodic arrangement is twice that of the single-row model. As briefly outlined below, this is in fact observed in the results of the computations.Figure 8
Schematic illustration of particles compact with periodic arrangement. Though the individual particles comprise multiple bilayers, only one is shown in the schematic. Periodicity is simulated by restricting the domain to the region lying between the dotted lines, and imposing periodic boundary conditions on these lines. Other boundary and initial conditions are similar to those depicted in Figure1.Figure9 shows curves of flame position versus time for different values of d. The results are obtained using the single-row particle model. Also shown is a curve obtained using the periodic configuration model with d=50 μm. The results obtained using the single-row model indicate that in the range of particle thicknesses considered, the average flame velocity generally decreases as d increases. However, the variations of Vavg with d are generally small, about 5% or less. Thus, in the present setup (where heat losses are ignored), the average flame speed appears to be less affected by particle thickness than by particle contact area.Figure 9
Flame front position versus time for the single-row configuration, with particle thickness in the ranged= 5–50 μm, as indicated. The symbols depict results obtained for a periodic configuration with d=50 μm.Also shown on Figure9 are the results obtained with the periodic stack model using d=50 μm. Note that the curve, obtained with the periodic model with d=50 μm, lies between d=20 μm and d=30 μm curves obtained using the single-row model. Thus, the periodic arrangement with d=50 μm yields similar results to the single-row model with d~25 μm, which illustrates our earlier claim regarding the correspondence between the two arrangements.


### 3.2.3. Effect of Thermal Contact Resistance

The thermal contact resistance in a particle compact may depend on a variety of factors, including geometry, applied pressure, as well as surface cleanliness, roughness, and flatness [40–42]. Nominal reported values of Rc fall in a wide range, between 10-7 and 10-5 m2 K/W. In this section, we briefly investigate the effect of thermal contact resistance on the average velocity of the front, namely, by repeating the analysis for a compact of 10 μm thick particles. For brevity, we consider a single value Rc=10-5 m2 K/W and as before vary Ac in the range 10–200 μm.The predicted values ofVavg are plotted against Ac in Figure 10; also shown for comparison is the curve obtained for perfect thermal contact Rc=0. The results indicate that, unlike the case of perfect thermal contact, for Rc=10-5 m2 K/W, the Vavg as a function of Ac is no longer monotonic. Specifically, Vavg peaks for Ac=100 μm. For lower values of Ac, the velocity rapidly decreases as Ac is reduced, whereas for higher values, the rate of decay is substantially smaller.Figure 10
Average flame speed versus contact area. Curves are generated for the case of vanishing thermal contact resistance, and forRc=10-5 m2 K/W. Also shown for comparison are the experimental average flame speed measurements of Fritz et al. [29], for particle compacts with δ in the range 26–375 nm.To investigate the origin of this trend, we plot in Figure11 the front position versus time for the different values of Ac. The results indicate that the observed trends in Vavg are essentially governed by variation of the ignition delay. Specifically, for small contact areas, the presence of a large contact resistance leads to a significant increase in the heating time, which delays the initiation of the front in the adjacent particle. This phenomenon has also been observed through detailed examination of temperature and heating rate distributions (not shown). On the other hand, in the range Ac≥100 μm, the predictions for Rc=0 and Rc=10-5 m2 K/W become close to each other. Thus, in this range, the impact of the thermal contact resistance becomes substantially weaker.Figure 11
Flame front position versus time forRc=10-5 m2 K/W. Curves are generated for different values of Ac, as indicated.We finally compare our present predictions with experimental measurements of Fritz et al. [29], who reported average flame speeds of 0.0058–0.0260 m/s for low-density particle compacts (packing densities of 20%) with different half-layer thickness. The data from Fritz et al. are also plotted in Figure 10. One observes that the range of the measured flame speeds is comparable with the computed average flame speed for a contact area Ac=10 μm and thermal contact resistance Rc=10-5 m2 K/W. Although differences between the computational and experimental results can arise due to a variety of factors that are not accounted for in the model, such as heat losses, particle geometry, and compact packing density, the agreement between reported measurements with predictions obtained with high Rc and small contact area (low packing density) strongly suggests that thermal contact resistance may play a strong role in determining the self-propagation velocity in practical compacts.


## 3.2.1. Effect of Contact Area

Simulations were conducted by systematically varyingAc in the range 10–200 μm; in all cases, the particle thickness was held fixed at d=10 μm. Note that for large values of Ac, with the current setup and ignition properties, the flame front quenches as it enters the contact region between the first and second particles. This is illustrated in Figure 4, which shows that for Ac=250 μm, Tmax drops to a low value, and the reaction is not initiated in the second particle. This observation is at the origin of the restricted range of Ac.Figure6 shows the dependence of Vavg on Ac for the case of perfect thermal contact. It is seen that Vavg monotonically decreases as Ac increases. For Ac=10 μm, the computed average velocity, Vavg=0.42 m/s, coincides with the flame speed of a continuous multilayered film of the same thickness d=10 μm. Thus, for this small value of Ac, the ignition delay associated with flame crossover is rather insignificant. As shown in Figure 4, the drop in the maximum temperature of the first particle is small for Ac≤50 μm. For Ac=100 μm and 200 μm, the impact of the ignition delay is significant, with predicted values Vavg=0.22 m/s and 0.15 m/s, respectively. These correspond to about a twofold or threefold decrease in Vavg compared to continuous multilayered film with uniform thickness.Figure 6
Profiles ofVavg versus Ac. The results are obtained assuming perfect thermal contact. The average speed for a continuous film is indicated using a dashed line for reference.To further illustrate the impact of the heating time,th, we plot in Figure 7 the flame position versus time; curves are generated for different values of Ac. The results indicate that during the front propagation phase, the slopes of the individual curves are essentially identical. Thus, for small values of Ac, the reaction time, that is, the time needed for the front to traverse an individual particle, is independent of the contact area. Accordingly, our calculations suggest that the average flame velocity in a compact of particles is close to that of a continuous film when the contact area is of the same order of magnitude as the particle thickness. The results also reveal that the substantial reduction in Vavg that occurs as Ac increases is mainly due to an increase in the particle heating time. It is also interesting to note that for the same configuration, the heating time may vary appreciably from one flame crossover event to another. For small values of Ac, the ignition delay essentially stabilizes after the first crossover. For Ac=200 μm, however, it takes more than a single passage for the ignition delay to stabilize. In particular, a large variation is observed between the first and the second flame crossings. The pronounced delay associated with the first crossing for Ac=200 μm is affected by the substantial heat loss that occurs immediately following ignition, due to the proximity between the ignition spot and the contact between the first and second particles.Figure 7
Flame front position versus time for different values ofAc, as indicated. The results are obtained assuming perfect thermal contact.The analysis above was repeated for a particle thicknessd=50 μm. The results (not shown) revealed similar trends for the dependence of Vavg on Ac. Small quantitative differences were however observed, as further discussed below.


## 3.2.2. Effect of Particle Thickness

In this section, we briefly examine the impact of the particle thickness,d, on the average flame velocity. We also take advantage of the present analysis to make a connection between the single-row arrangement shown in Figure 1 and a periodic configuration depicted in Figure 8. In the latter, the stack used in Figure 1 is repeated multiple times in the vertical direction. Assuming that the leftmost particles are ignited identically and simultaneously, this is modelled as a periodic arrangement as indicated in Figure 8. Specifically, periodic boundary conditions are used at the horizontal boundaries of the domain, indicated in Figure 8 using dashed lines. Based on symmetry arguments, predictions obtained using the periodic configuration should correspond to those obtained using the single-row model, provided that the same physical parameters are used except for the particle thickness, whose value in the periodic arrangement is twice that of the single-row model. As briefly outlined below, this is in fact observed in the results of the computations.Figure 8
Schematic illustration of particles compact with periodic arrangement. Though the individual particles comprise multiple bilayers, only one is shown in the schematic. Periodicity is simulated by restricting the domain to the region lying between the dotted lines, and imposing periodic boundary conditions on these lines. Other boundary and initial conditions are similar to those depicted in Figure1.Figure9 shows curves of flame position versus time for different values of d. The results are obtained using the single-row particle model. Also shown is a curve obtained using the periodic configuration model with d=50 μm. The results obtained using the single-row model indicate that in the range of particle thicknesses considered, the average flame velocity generally decreases as d increases. However, the variations of Vavg with d are generally small, about 5% or less. Thus, in the present setup (where heat losses are ignored), the average flame speed appears to be less affected by particle thickness than by particle contact area.Figure 9
Flame front position versus time for the single-row configuration, with particle thickness in the ranged= 5–50 μm, as indicated. The symbols depict results obtained for a periodic configuration with d=50 μm.Also shown on Figure9 are the results obtained with the periodic stack model using d=50 μm. Note that the curve, obtained with the periodic model with d=50 μm, lies between d=20 μm and d=30 μm curves obtained using the single-row model. Thus, the periodic arrangement with d=50 μm yields similar results to the single-row model with d~25 μm, which illustrates our earlier claim regarding the correspondence between the two arrangements.


## 3.2.3. Effect of Thermal Contact Resistance

The thermal contact resistance in a particle compact may depend on a variety of factors, including geometry, applied pressure, as well as surface cleanliness, roughness, and flatness [40–42]. Nominal reported values of Rc fall in a wide range, between 10-7 and 10-5 m2 K/W. In this section, we briefly investigate the effect of thermal contact resistance on the average velocity of the front, namely, by repeating the analysis for a compact of 10 μm thick particles. For brevity, we consider a single value Rc=10-5 m2 K/W and as before vary Ac in the range 10–200 μm.The predicted values ofVavg are plotted against Ac in Figure 10; also shown for comparison is the curve obtained for perfect thermal contact Rc=0. The results indicate that, unlike the case of perfect thermal contact, for Rc=10-5 m2 K/W, the Vavg as a function of Ac is no longer monotonic. Specifically, Vavg peaks for Ac=100 μm. For lower values of Ac, the velocity rapidly decreases as Ac is reduced, whereas for higher values, the rate of decay is substantially smaller.Figure 10
Average flame speed versus contact area. Curves are generated for the case of vanishing thermal contact resistance, and forRc=10-5 m2 K/W. Also shown for comparison are the experimental average flame speed measurements of Fritz et al. [29], for particle compacts with δ in the range 26–375 nm.To investigate the origin of this trend, we plot in Figure11 the front position versus time for the different values of Ac. The results indicate that the observed trends in Vavg are essentially governed by variation of the ignition delay. Specifically, for small contact areas, the presence of a large contact resistance leads to a significant increase in the heating time, which delays the initiation of the front in the adjacent particle. This phenomenon has also been observed through detailed examination of temperature and heating rate distributions (not shown). On the other hand, in the range Ac≥100 μm, the predictions for Rc=0 and Rc=10-5 m2 K/W become close to each other. Thus, in this range, the impact of the thermal contact resistance becomes substantially weaker.Figure 11
Flame front position versus time forRc=10-5 m2 K/W. Curves are generated for different values of Ac, as indicated.We finally compare our present predictions with experimental measurements of Fritz et al. [29], who reported average flame speeds of 0.0058–0.0260 m/s for low-density particle compacts (packing densities of 20%) with different half-layer thickness. The data from Fritz et al. are also plotted in Figure 10. One observes that the range of the measured flame speeds is comparable with the computed average flame speed for a contact area Ac=10 μm and thermal contact resistance Rc=10-5 m2 K/W. Although differences between the computational and experimental results can arise due to a variety of factors that are not accounted for in the model, such as heat losses, particle geometry, and compact packing density, the agreement between reported measurements with predictions obtained with high Rc and small contact area (low packing density) strongly suggests that thermal contact resistance may play a strong role in determining the self-propagation velocity in practical compacts.


## 4. Conclusions

In this work, we developed a simple numerical model for the simulation of self-propagation of reactive fronts in an idealized compact. The compacts comprise identical Ni/Al multilayered particles in thermal contact. The evolution of the reaction front in the compact was simulated using a two-dimensional transient reduced reaction model. The model was used to investigate the propagation of the reaction front within the compact and to analyze the dependence of its average velocity on the contact area between neighboring particles, the particle thickness, and the thermal contact resistance. Analysis of the computations revealed that(1)
The average velocity of the front can be substantially smaller than that of a continuous multilayered foil of the same microstructure. Consistent with the analysis of [29], the reduction in average flame speed is primarily influenced by a drop in flame temperature and reaction rates as the flame enters the contact region between neighboring particles, which results in a delay in the crossover of the reaction front.(2)
For the case of perfect thermal contact, the average flame velocity decreases asAc increases, in the considered range of contact areas. The same trend is observed for all considered values of the particle thickness.(3)
The average flame velocity decreases slightly as the particle thickness increases. With the present adiabatic computations, the dependence on particle diameter appears insignificant compared to variations induced byAc.(4)
The thermal contact resistance can have a substantial impact on the average front velocity, especially for low values ofAc. In this regime, incorporation of a thermal contact resistance leads to further reduction in Vavg. At larger values of Ac, however, the impact of Rc becomes much weaker.(5)
The comparison of computed predictions with experimental measurements of Fritz et al. reveals a favorable agreement with predictions obtained using highRc; this suggests that the thermal contact resistance plays an important role in determining the front velocity in practical particle compacts.Work is currently underway to extend the present model to account for realistic, three-dimensional microstructures, as well as the effect of radiative and convective heat losses.



---

*Source: 198096-2013-03-24.xml*