#### 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.

#### 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, , 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]

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 Figure 1. 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 and thickness , 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, , 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: where is the region-averaged enthalpy, is the area of the region, is the region-averaged heat release term, is the mean specific heat, is the nominal flame temperature, and is the (negative) heat of mixing. Note that we have assumed that exhibits a quadratic dependence on the concentration, , according to [33, 34] where is a dimensionless concentration defined so that in the reactants, and for the product. Also note that the temperature, , can be recovered from the enthalpy , following the methodology outlined in [34]; namely, where K, K, and K denote the melting temperatures of Al, Ni, and NiAl, respectively, , , and are the corresponding heats of fusion (per unit mole), represents the fraction of pure (unmixed) Al, , , , , , , , , , and are the densities of Al and Ni, and and are the corresponding heat capacities, whereas , , and denote the atomic weights of Al, Ni, and NiAl, respectively.

The evolution of the mean concentration and its second moment, and , 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
where denotes atomic diffusivity and is the half-thickness of an Al layer. The diffusivity is assumed to follow an Arrhenius dependence on , according to
where m^{2}/s is the preexponent, kJ/mol is the activation energy, and is the universal gas constant. The values of and 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 flux which is given by Fourier’s law: where is the thermal conductivity. Generally, 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: where and 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 domain that is discretized with a uniform Cartesian grid of mesh size and along the and directions, respectively; see Figure 1. The two field variables, 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 where is the augmented value of the thermal conductivity at the solid-solid boundaries, and 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, mm. The microstructure of the particles is also held fixed, with the Al layer half-thickness nm and the premix width nm. Attention is thus focused on the effects of contact area, , and thermal contact resistance, .

The self-propagating reactions are initiated using a thermal spark of temperature K and width *μ*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:
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 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

Figure 2 shows the flame position (-direction) versus time for a compact comprising five layered particles of thickness and contact area . Once established, the front initially propagates with uniform velocity in the first particle, ms. 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 ms (at 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 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.

To further investigate the mechanism of front crossover, we plot in Figure 3 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 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, K, and the width of the reaction zone is around m. At , 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, K. This is accompanied by a dramatic decrease in the peak instantaneous heat release rate, from about 7 MW/m^{3} at ms to about MW/m^{3} at ms, and a substantial widening of the reaction zone width. Figure 3(c) shows that, at ms, the reactive front is still located within the contact area. The peak temperature has dropped to 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 ms and ms, though the peak heat release rate has increased to MW/m^{3}. At 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/m^{3}) once again reach values associated with a steadily propagating front.

**(a)**

**(b)**

**(c)**

**(d)**

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 Figure 5(a), which clearly shows that the heat release rate drops to very small values during the heating phase (). Another perspective is provided from section-averaged temperature profiles plotted in Figure 5(b). It is interesting to note that at 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.

**(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, , on bilayer and premix parameters [30]. These studies have shown that for steadily propagating fronts, 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 of by taking the ratio of the total length of the compact, , by the total reaction time, , defined as the time needed for reaction front to reach the end of the last particle. Here, 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 of on the number of particles present in the compact, namely, by systematically increasing in the range . The results (not shown) indicated that as becomes large, tends to an asymptotic value. Specifically, when , falls within 5% of each other, in the entire range of that is considered. Consequently, we have used the estimates of for in all the computations presented below.

###### 3.2.1. Effect of Contact Area

Simulations were conducted by systematically varying in the range 10–200 *μ*m; in all cases, the particle thickness was held fixed at *μ*m. Note that for large values of , 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 *μ*m, 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 .

Figure 6 shows the dependence of on for the case of perfect thermal contact. It is seen that monotonically decreases as increases. For *μ*m, the computed average velocity, m/s, coincides with the flame speed of a continuous multilayered film of the same thickness *μ*m. Thus, for this small value of , 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 *μ*m. For *μ*m and *μ*m, the impact of the ignition delay is significant, with predicted values m/s and 0.15 m/s, respectively. These correspond to about a twofold or threefold decrease in compared to continuous multilayered film with uniform thickness.

To further illustrate the impact of the heating time, , we plot in Figure 7 the flame position versus time; curves are generated for different values of . The results indicate that during the front propagation phase, the slopes of the individual curves are essentially identical. Thus, for small values of , 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 that occurs as 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 , the ignition delay essentially stabilizes after the first crossover. For *μ*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 *μ*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.

The analysis above was repeated for a particle thickness *μ*m. The results (not shown) revealed similar trends for the dependence of on . 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, , 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 9 shows curves of flame position versus time for different values of . The results are obtained using the single-row particle model. Also shown is a curve obtained using the periodic configuration model with *μ*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 increases. However, the variations of with 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.

Also shown on Figure 9 are the results obtained with the periodic stack model using *μ*m. Note that the curve, obtained with the periodic model with *μ*m, lies between *μ*m and *μ*m curves obtained using the single-row model. Thus, the periodic arrangement with *μ*m yields similar results to the single-row model with *μ*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 fall in a wide range, between and m^{2} 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 m^{2} K/W and as before vary in the range 10–200 *μ*m.

The predicted values of are plotted against in Figure 10; also shown for comparison is the curve obtained for perfect thermal contact . The results indicate that, unlike the case of perfect thermal contact, for m^{2} K/W, the as a function of is no longer monotonic. Specifically, peaks for *μ*m. For lower values of , the velocity rapidly decreases as is reduced, whereas for higher values, the rate of decay is substantially smaller.

To investigate the origin of this trend, we plot in Figure 11 the front position versus time for the different values of . The results indicate that the observed trends in 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 *μ*m, the predictions for and m^{2} K/W become close to each other. Thus, in this range, the impact of the thermal contact resistance becomes substantially weaker.

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 ) 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 *μ*m and thermal contact resistance m^{2} 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 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 as 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 by .(4)The thermal contact resistance can have a substantial impact on the average front velocity, especially for low values of . In this regime, incorporation of a thermal contact resistance leads to further reduction in . At larger values of , however, the impact of becomes much weaker.(5)The comparison of computed predictions with experimental measurements of Fritz et al. reveals a favorable agreement with predictions obtained using high ; 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.

#### Acknowledgments

This work was supported by the Office of Naval Research through Award N00014-07-1-0740 and by the Defense Threat Reduction Agency, Basic Research Award no. HDTRA1-11-1-0063, to Johns Hopkins University.