Self-Propagating Reactive Fronts in Compacts of Multilayered Particles

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 20m/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 selfpropagation 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.


Introduction
Reactive multilayered materials have recently gained increasing interest in various applications, including joining, brazing, sealing, and ignition of secondary reactions [1][2][3][4][5][6][7][8][9][10][11][12][13].Typically, these materials are fabricated in the form of multilayered foils [13][14][15][16][17][18][19][20] using vapor deposition techniques [13,14,21,22], ball-milling [23][24][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] 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 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, .The particles are identical of fixed length  = 1 mm.The contact area,   , between the particles and thickness  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   and width   .The grid is schematically shown on the second and third particles.Similar gridding is applied to all particles but is not shown.
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 selfpropagating 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.

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.

Reaction Model.
The evolution of the reaction within the particle compact is analyzed using the reduced model developed in [30][31][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 || = 1 in the reactants, and  = 0 for the product.Also note that the temperature, , can be recovered from the enthalpy , following the methodology outlined in [34]; namely, where The evolution of the mean concentration and its second moment,  and  2 , 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  0 = 2.18 × 10 −6 m 2 /s is the preexponent,  = 137 kJ/mol is the activation energy, and  is the universal gas constant.The values of  0 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 q 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  Al and  Ni 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 Δ 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  mod is the augmented value of the thermal conductivity at the solid-solid boundaries, and   is the magnitude of the thermal contact resistance.

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,  = 1mm.The microstructure of the particles is also held fixed, with the Al layer half-thickness  = 250 nm and the premix width  = 0.8 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   = 1400 K and width   = 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: 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.

Front Propagation in Particle
Compact. Figure 2 shows the flame position (-direction) versus time for a compact comprising five layered particles of thickness 10 m and contact area   = 200 m.Once established, the front initially propagates with uniform velocity in the first particle,  ≤ 2 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  = 12 ms (at   = 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.
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  = 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,  max = 1912 K, and the width of the reaction zone is around 100 m.At  = 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,  max = 1846 K.This is accompanied by a dramatic decrease in the peak instantaneous heat release rate, from about 7 MW/m 3 at  = 1 ms to about 0.05 MW/m 3 at  = 5 ms, and a substantial widening of the reaction zone width.Figure 3(c) shows that, at  = 12 ms, the reactive front is still located within the contact area.The peak temperature has dropped to  max = 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  = 5 ms and  = 12 ms, though the peak heat release rate has increased to 0.14 MW/m 3 .At  = 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/m 3 ) once again reach values associated with a steadily propagating front.
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 (5 ms ≤  ≤ 12ms).Another perspective is provided from section-averaged temperature profiles plotted in Figure 5(b).It is interesting to note that at  = 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.
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.

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][38][39]. The reduced model used in this work has been successful in capturing the dependence of the average front velocity,  avg , on bilayer and premix parameters [30].These studies have shown that for steadily propagating fronts,  avg 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  avg by taking the ratio of the total length of the compact,   =  − (−1)  , 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  avg on the number of particles present in the compact, namely, by systematically increasing  in the range 1 ≤  ≤ 21.The results (not shown) indicated that as  becomes large,  avg tends to an asymptotic value.Specifically, when  ≥ 5,  avg falls within 5% of each other, in the entire range of   that is considered.Consequently, we have used the estimates of  avg for  = 5 in all the computations presented below.cases, the particle thickness was held fixed at  = 10 m.

Effect of Contact
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   = 250 m,  max 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  avg on   for the case of perfect thermal contact.It is seen that  avg monotonically decreases as   increases.For   = 10 m, the computed average velocity,  avg = 0.42 m/s, coincides with the flame speed of a continuous multilayered film of the same thickness  = 10 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   ≤ 50 m.For   = 100 m and 200 m, the impact of the ignition delay is significant, with predicted values  avg = 0.22 m/s and 0.15 m/s, respectively.These correspond to about a twofold or threefold decrease in  avg 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  avg 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   = 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   = 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.
The analysis above was repeated for a particle thickness  = 50 m.The results (not shown) revealed similar trends for the dependence of  avg on   .Small quantitative differences were however observed, as further discussed below.

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 singlerow particle model.Also shown is a curve obtained using the periodic configuration model with  = 50m.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  avg 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  = 50 m.Note that the curve, obtained with the periodic model with  = 50 m, lies between  = 20 m and  = 30 m curves obtained using the single-row model.Thus, the periodic arrangement with  = 50 m yields similar results to the single-row model with  ∼ 25 m, which illustrates our earlier claim regarding the correspondence between the two arrangements.variety of factors, including geometry, applied pressure, as well as surface cleanliness, roughness, and flatness [40][41][42].

Effect of
Nominal reported values of   fall in a wide range, between 10 −7 and 10 −5 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  The predicted values of  avg are plotted against   in Figure 10; also shown for comparison is the curve   obtained for perfect thermal contact   = 0.The results indicate that, unlike the case of perfect thermal contact, for   = 10 −5 m 2 K/W, the  avg as a function of   is no longer monotonic.Specifically,  avg peaks for   = 100 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  avg 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   ≥ 100 m, the predictions for   = 0 and   = 10 −5 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.0260m/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   = 10 m and thermal contact resistance   = 10 −5 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 selfpropagation velocity in practical compacts.

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 twodimensional 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  avg .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.

Figure 2 :
Figure 2: Flame front position versus time for contact area   = 200 m with perfect thermal contact.The dashed lines indicate the region at the contact area between successive particles.

Figure 3 :
Figure 3: Instantaneous distributions of temperature in K (left) and local heat release rate in MW/m 3 (right) at selected times ((a)-(d)).The contact area   = 200 m and perfect thermal contact is assumed.Distributions are shown on the front crossover between the first and second particles.
Area.Simulations were conducted by systematically varying   in the range 10-200 m; in all

Figure 4 :
Figure 4: Evolution of the maximum temperature within the compact for different values of   , as indicated.The results are obtained assuming perfect thermal contact.
Thermal Contact Resistance.The thermal contact resistance in a particle compact may depend on a ms  = 5 ms  = 12 ms  = 13.6 ms  = 15 ms (b)

Figure 5 :
Figure 5: Instantaneous (a) local heat release rate and (b) temperature profiles at selected time instants for   = 200 m with perfect thermal contact.The -axis spans the first two particles in the compact.The overlap region between the two particles is shown as a dashed curve.

Figure 6 :
Figure 6: Profiles of  avg versus   .The results are obtained assuming perfect thermal contact.The average speed for a continuous film is indicated using a dashed line for reference.

Figure 7 :
Figure 7: Flame front position versus time for different values of   , as indicated.The results are obtained assuming perfect thermal contact.

Figure 8 :
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 Figure 1.

Figure 9 :
Figure 9: Flame front position versus time for the single-row configuration, with particle thickness in the range  = 5-50 m, as indicated.The symbols depict results obtained for a periodic configuration with  = 50 m.

Figure 10 :
Figure 10: Average flame speed versus contact area.Curves are generated for the case of vanishing thermal contact resistance, and for   = 10 −5 m 2 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.

Figure 11 :
Figure 11: Flame front position versus time for   = 10 −5 m 2 K/W.Curves are generated for different values of   , as indicated.