Prediction and characterization of heat-affected zone formation in tin-bismuth alloys due to nickel-aluminum multilayer foil reaction

Reactive multilayer foils have the potential to be used as local high intensity heat sources for a variety of applications. Most of the past research effort concerning these materials have focused on understanding the structure-property relationships of the foils that govern the energy released during a reaction. To improve the ability of researchers to more rapidly develop technologies based on reactive multilayer foils, a deeper and more predictive understanding of the relationship between the heat released from the foil and microstructural evolution in the neighboring materials is needed. This work describes the development of a numerical model for the purpose of predicting heat affected zone size in substrate materials. The model is experimentally validated using a commercially available Ni-Al multilayer foils and alloys from the Sn-Bi binary system. To accomplish this, phenomenological models for predicting the variation of physical properties (i.e., thermal conductivity, density, and heat capacity) ...


I. INTRODUCTION
2][3][4][5][6][7] These materials utilize alternating submicron layers of reacting elements to generate heat.The alternating layers provide short diffusion lengths for mixing of the reactive elements to occur.Once energy is added to the multilayer, the enthalpy of formation creates a selfsustaining exothermic reaction, the result of which is a rapid release of energy that propagates quickly.The amount of energy released and the rate of reaction are dictated by the structure and chemistry of the multilayer. 8This allows for a great degree of control of the RMF reactivity and overall behavior making these materials ideal for a number of applications, such a joining, igniters and power sources. 92][3][4][5][6][7] Little is understood regarding the effect of RMFs on neighboring materials (such as a substrate) and the relationship between heat transport in the bulk material and its consequential effect on the microstructure evolution of the ensemble.(Throughout this work, the word "substrate" refers to the bulk Sn-Bi alloy to which the multilayer is bonded through reaction.This is not the same substrate material used for sputter deposition of reactive multilayers.) This work seeks to develop a better understanding of how RMFs affect the underlying material microstructure by investigating the heat flow and its coupling to phase transformations that may occur.To this end, a numerical model was developed based on the combination of foil and substrate properties.This model was used to predict the scale of the heat-affected zone (HAZ) in four compositions of the Sn-Bi alloy system.The Sn-Bi system was chosen due to its low melting temperature and simple phase diagram, which would allow for larger HAZs, easing microstructure analysis. 21hrough concurrent computational and experimental characterization of the substrate's HAZ, a greater understanding of the underlying physics governing these reactions was obtained.

II. MODEL DEVELOPMENT
The purpose of the numerical model was to create a means by which the HAZ depth can be predicted chiefly based on foil and substrate properties.Additionally, it was desirable for the final model not be computationally intensive to allow for broader utilization.In Secs.II A and II B, the different components of the model are presented and their roles are discussed individually.
A. Heat flow model

Thermal transport in the substrate
Heat flow in the substrate material was handled with a one dimensional (1-D) formulation of the heat equation (Eq.( 1)) 22 @T @t ¼ a @ 2 T @x 2 ; (1) where a) In Eq. ( 1), the variables T, x, a, and t are the temperature, distance, thermal diffusivity, and time, respectively.a is comprised of the thermal conductivity (k) in W/m-K, specific heat capacity (C p ) in J/kg-K, and density (q) in kg/m 3 .The 1-D heat equation was selected in order to simplify the calculation.Additionally, no microstructural constituents were included in the calculation and the substrate material was treated as a continuum (e.g., homogeneous) with the same properties as the bulk alloy.Thus, the only phase transformations captured in this simplified model is the onset of melting and solidification, which occur over the range of temperatures between the solidus (T solidus ) and liquidus (T liquidus ) temperatures.
For implementation into the code, Eq. ( 1) was approximated using a forward in time central difference finite difference approximation (Eq.( 3)).The variable, m, represents the time step that the value of T is taken from and n is the node index Equation ( 3) was further simplified to take the form shown below where The values of Dt and Dx are determined by the time step and the distance between nodes used in the calculation.For Eq. ( 4) to accurately model the heat flow, the R must be less than 0.5, 23 otherwise the solutions will diverge from those based on the differential form of the heat equation.This will result in non-physical temperature predictions within the substrate.For the simulations discussed in this work, the values for Dt and Dx were taken to be 1 ls and 10 lm, respectively.These values were chosen in order to balance the spatial resolution of the output HAZ size, calculation stability, and calculation time.

Boundary conditions
The heat evolved from the multilayer was calculated based on the structure of the foil, assuming discrete transitions in composition between layers after Michaelsen (Eq.( 6)). 8 Energy output (DH) is based on the enthalpy of formation for the compound being formed during the reaction (DH 0 ), the bilayer thickness (k) and the thickness of the intermixed region between two adjoining layers (w).In this work, DH 0 was taken to equal 1383 J/g assuming the formation of an equiatomic NiAl compound 24 The energy from the RMF reaction was converted to an initial foil temperature value by assuming that the energy was reversibly absorbed by the reacted foil, resulting in an increase in temperature. 25Energy transport across the foilsubstrate interface was calculated utilizing a form of Fourier's law, 22 Eq.( 7).This equation was also approximated using a finite difference approximation In the above equation, the variables q f , x f , and k s represent the density of the foil, the foil thickness, and the thermal conductivity of the substrate, respectively.Heat was then conducted through the substrate using Eq. ( 4).The opposite boundary was held constant at 300 K, due to the relative thickness of the substrate material compared to the amount of heat released by the foil.Additionally, it is assumed that the only loss of heat from the foil is due to conduction into the substrate with negligible losses due to convection in the atmosphere and radiation.
Since the reaction of the foils can generate extreme local temperatures melting had to be accounted for in the substrate.This was handled by modifying the C p value used in the calculations depending on temperature.This is shown in Eq. ( 8), 26 where the enthalpy of fusion (DH fus ) and the melting range (DT m ) are used when a node associated with the substrate reached a temperature within the alloy melting range.Here, the subscripts s and l represent the solid and liquid phases, respectively.Due to the rapid heating and cooling conditions created by the RMF reaction, the HAZ was said to extend to the maximum distance, where T liquidus was achieved in the substrate.This was thought to be necessary since the kinetics in the liquid state would allow for sufficient mass diffusion in order to facilitate the redistribution of the primary phase, which creates the observable HAZ structure

B. Modeling physical properties of substrate
Composition and temperature dependent material properties were estimated using a phenomenological solution model that accounted for deviations from ideality using the Redlich-Kister polynomial. 27The Redlich-Kister polynomial (Eq.9][30] This methodology was chosen due to its ability to provide a consistent means to model all of the properties of interest.Additionally, this approach can be easily expanded and applied to more complex systems (i.e., ternary systems and beyond), which allows for greater utility in the future.The model took the basic form of below where In the above formulae, P is the property being modeled, and P 1 and P 2 are the contribution to the property from pure component 1 and the contribution from pure component 2, respectively.The properties of interest, P, in this work are the substrate thermal conductivity and density.The contribution to the overall property is modulated by the atomic fraction of each component in the system, X 1 and X 2 .The deviations from ideality are accounted for with the interaction parameters (L i ) in the P excess term, which can be a positive or negative constant or have temperature dependence depending on the property being fit.Overall, the model implementation is represented schematically in Fig. 1.

Density
A relationship for the variation of volume as a function of temperature is well established 25 (Eq.( 11)).In the equation below, a CTE , V, and P are the volumetric coefficient of thermal expansion, volume, and pressure, respectively When a constant mass assumption is made, the relation can be rewritten to model changes in density as a function of temperature, where q 0 is the density at a reference temperature T 0 Literature data were used to fit the relationship between temperature and density (Eq.( 12)) for the pure elements 24,[31][32][33][34] and subsequently, this information was leveraged to model the variation in density that can be expected as a function of both temperature and composition in the Sn-Bi system based on available literature from different alloys. 35,36Table I contains the coefficients for the models of density for the Sn-Bi system in the solid and liquid states.

Thermal conductivity
It is known that the variation of thermal conductivity in a pure metal with temperature can be modeled using Eq. ( 13). 22Data taken from literature were used to account for this temperature dependence for both Sn 24,37-39 and Bi 24,38,40 kðTÞ ¼ a þ bT: Once this behavior was modeled, additional literature data 38,40 were used to create a model for the calculation of the temperature and compositional dependence of thermal conductivity for liquid Sn-Bi alloys (Table II).The R 2 -value is presented along with the model coefficients to quantify the fit.The residuals from the fit as a function of temperature and composition (not shown) were also analyzed to eliminate artifacts in the model that could originate from over or under fitting the data.The resultant model was implemented into the code to account for changes in thermal conductivity as a result of temperature fluctuations.

Heat capacity
C p as a function of temperature (Eq.( 14)) was modeled using the standard equation for molar heat capacity 41 divided by the molar mass of the alloy (M).The coefficients used for the Sn-Bi calculations are shown below (Table III).It has been established that C p as a function of composition can be well modeled by using an ideal solution model 42,43 (i.e., where the excess term is equal to zero).Therefore, this was used for calculating C p in the simulation

III. EXPERIMENTAL PROCEDURES
In order to validate the model, four experimental alloys from the Sn-Bi system were used: Sn with 15, 39, 57, and 78 wt.% Bi (Fig. 2).The alloys were cast from Sn shot (99.8% metal basis) and Bi needles (99.99% metal basis) sourced from Alfa Aesar.Each alloy was melted under ambient atmosphere in a graphite crucible located inside a resistance box furnace at 400 C. The melt was held at temperature for 20 min; they were then removed from the furnace and stirred with a graphite rod.After the alloys were stirred, they were placed back in the furnace for approximately 5 min before casting.Each alloy was cast into a room temperature copper mold with dimensions of approximately 19 mm Â 19 mm Â 100 mm.
To verify composition, a sample from each casting was removed and analyzed using X-ray fluorescence spectroscopy (XRF).The analysis was performed on an EDAX Eagle III with an accelerating voltage of 40 kV and a spot size of 90 lm.Three spots were collected from each alloy and used to determine the average composition of the casting.Prior to analysis, the surface of each of the samples was dry-polished with 1200 grit paper to remove any surface oxide that may have formed.
The solid state thermal conductivity values of the experimental alloys were determined following ASTM standard E1461-07 44 on an Anter Flashline 4010 system, which employed the laser flash method to determine thermal diffusivity.Specimens from each alloy were lathed into cylinders with diameters between 10.4 and 10.8 mm and ground to thicknesses of 1.0-2.0mm per the testing fixture geometry requirements.The samples were tested at three temperatures (40 C, 60 C, and 80 C) with five measurements taken at each condition.Five shots of each sample were performed with a 1300 W laser, and the thermal diffusivity was calculated utilizing the Clark and Taylor correction method. 45Then along with experimentally measured C p values relative to a Mo standard (also measured by the Anter Flashline 4010 system), thermal conductivity was calculated through Eq. ( 2).
Differential scanning calorimetry (DSC) was utilized to determine the solidus and liquidus phase transformation temperatures, and DH fus values.The experiments were conducted on a Perkin Elmer DSC 8000.Prior to testing, the machine was calibrated using In and Zn standards to verify that DH fus and transformation temperatures were accurate.Samples weighing between 150 and 300 mg were placed in sealed Al crucibles and the temperature was ramped from 30 C to 400 C for each of the alloys.All of the experiments were conducted at a ramp rate of 20 C/min.Each of the alloys was analyzed three times to ensure consistency in the measurement.Values for DH fus and transformation temperatures were determined using the area under the heat flow curve and tangent line intersection, respectively.
The structure of the as-received RMFs was characterized with transmission electron microscopy (TEM).Samples were made using a FEI Dual-Beam Strata DB235 focusedion beam microscope.A protective layer of Pt was deposited on the surface and a series of milling operations were carried out with the ion beam until a foil <300 nm in thickness was created.The sample was then lifted out ex situ using micromanipulators and placed on a Cu sample holder.A JEOL 200CX TEM was used to characterize the foil structure.
To experimentally characterize the HAZ size, free standing multilayers were reacted in contact with each of the alloys.The sputter deposited Ni(V)-Al multilayer foils were sourced from Indium Corporation.Each multilayer was approximately 12.5 mm Â 12.5 mm Â 60 lm thick and had a 1 lm thick layer of InCuSil V R -ABA TM (Ag 59 Cu 27.25 In 12.5 Ti 1.25 ), a braze alloy, on the surface.The experimental set up of the multilayer-substrate reaction is shown schematically in Fig. 3.This reaction geometry was chosen to improve the contact between the substrate and the free standing foils.Preliminary experimentation with the free standing foils indicated that the unbound foils tended to distort their shape significantly as a result of the reaction.Therefore, it was determined that placing the foil beneath the substrate would allow for the most consistent results.Prior to the reaction, the substrate materials were prepared by sectioning the casting in to square specimens with thickness of 6-8 mm.This allowed for sufficient thickness for virtually all of the heat from the multilayer to be absorbed into the sample without the need to consider transport effects on the opposite side of the sample.Additionally, it allowed for many samples to be studied from the individual castings, which maximized the data obtainable from the limited material.The samples were then planed using a lapping fixture on a polishing wheel with 800 grit grinding paper.After the samples were planed, the face to be in contact with the foil was polished down to a 0.04 lm colloidal silica final polish using standard metallographic techniques.The surface was prepared in this manner to reduce the influence of contact resistance due to surface roughness. 46Once the substrate materials were polished, the reaction was immediately carried out to mitigate any effects of surface oxide formation.
Once the specimen was in position, the reaction was initiated using a 20 V DC power supply.An arc was created by connecting the power supply to two solid core copper wires that were brought into contact with the foil.For consistency between reactions, the power supply was set to its maximum of 20 V for every reaction.
HAZ thickness measurements were performed using optical microscopy due to a combination of the scale of the reaction layer and the contrast between the phases present in the microstructure.Images of the HAZ were collected along the entire length of all samples with the exception of the extreme ends of the foils.These regions were removed from the calculation of the thickness due to edge effects from local material flow into these regions during the reaction, exacerbating the measurements taken from the sample.Analysis of the thickness was performed using ImageJ. 47,48For each image, five measurements were taken of the HAZ.The measurements were taken by drawing a line perpendicular to the foil-substrate interface to the edge of the HAZ.Higher resolution microstructural analysis was performed on a Phenom Pro scanning electron microscope (SEM).

A. Experimentally determined model parameters
The compositions of the four alloys determined using XRF are shown in Table IV.
The measurements of thermal conductivity are shown in Table V.In general, the experimental values follow the trend shown in literature even though the data are from different temperatures. 38The minimal difference between the magnitudes of thermal conductivity at different temperatures is due to its limited temperature dependence in these alloys.For this reason and due to the limited temperature range between room temperature and the solidus temperature in these alloys (max difference 142 K), the thermal conductivity was averaged over the three temperatures tested and held constant for the simulations involving each alloy.
The values for the transition temperatures (T solidus and T liquidus ) and DH fus were obtained from DSC analysis.Representative heating curves for each of the alloys are presented in Fig. 4. All of the transition temperatures corresponded to those expected based on the phase diagram. 21he initial peak in the 15Bi curve indicated that a small amount of eutectic was retained in the microstructure under the cooling conditions experienced during casting.This behavior is not predicted by the phase diagram (see Fig. 2).
The final parameters needed for the model were related to the foil structure.Using TEM, the value of k for the Ni-Al RMFs was determined to be 47.6 6 1.3 nm.A representative micrograph of the foil structure is presented in Fig. 5. Due to resolution limitations of the microscope, the intermixed region could not be measured.However, in a study conducted by   Wang et al., 3 which studied similar foils, the thickness of the intermixed layer was determined to be 2.3 nm.

B. Experimental HAZ vs predictions
To validate the assumption that the effect of heat loss due to convective and radiative cooling would be negligible, the temperature of the foil as a function of time was calculated assuming heat loss due solely to convection, radiation, and conduction, respectively.The calculations assumed that heat loss could occur through both sides of the foil, in the case of convective and radiative heat loss, and only through the foil-substrate interface in the case of conduction.Additionally, the conduction calculation was performed using the 78Bi alloy material properties, which had the lowest thermal conductivity of the alloys studied and a low thermal conductivity when compared to most metals and alloys, in general. 24As can be seen (Fig. 6), cooling in the foil occurs much faster in the case of conduction than in the other two cases.This supports the decision to not include the other two heat loss mechanisms in the overall simulation.Once the assumption was validated and the requisite data were collected, simulations for each of the four alloys were performed.The simulation predictions for HAZ size were then compared to the corresponding HAZs experimentally measured in each alloy.
In general, the HAZ microstructure morphology (Fig. 7) can be explained on the basis of rapid solidification front propagation caused by the high cooling rates. 49As the cooling rate is increased, so too is the solidification front velocity.At sufficiently high solidification rates, partitioning of solutes is reduced or even prevented entirely, which can result in a single phase structure. 50This phenomenon is manifested in these alloys in the form of the more homogenous distribution of phases in the HAZ microstructures.
With the exception of the 78Bi alloy, the predicted value fell within one standard deviation of the average measured HAZ size for the alloy (Fig. 8(a)).The large deviations in the HAZ size are hypothesized to be due to the distortion of the foil during the reaction.Even though the foil was held in place by the mass of the substrate in the experimental setup distortion of the foil was still possible.This led to local variations in heat conduction and a highly variable HAZ thickness across the length of the sample (Fig. 8(b)).It is worth noting that the interfacial region of the RMF and the substrate had no evidence of any reaction products from interactions between the foil, InCuSil V R -ABA TM capping layer, and the underlying alloys.There was evidence of melt infiltration between the foil grains local to the interface (not shown), which indicates good adherence between the RMF and substrate.
In order to account for the inconsistent heat flow across the foil-substrate interface into the substrate alloy and capture the resulting variability in HAZ size, the initial foil temperature in the simulation was varied.The simulation was then run 50 times for each alloy (equivalent to the number of individual experimental measurements of the HAZ).Each of the initial temperatures was determined by randomly selecting values from a normal distribution centered on the calculated initial temperature, which was based on the foil structure and chemistry.The standard deviation of the distribution was determined by averaging the standard deviation (as a percentage of the mean) for each of the four experimental alloys.The results of these simulations are presented in Fig. 9.In each of the alloys except 78Bi, the average predicted HAZ size fell within one standard deviation of the average measured value.In addition, the stochastic nature of the heat flow across the foil-substrate interface appears to be captured over the multiple simulations.
With the variability of HAZ size considered, the compositional effect on HAZ size appears to be minimal in the Sn-Bi system.This could be due to the fact that melting is initiated in each of the alloys at the same temperature (recall that eutectic was present in each of the initial alloy microstructures) and the values for T liquidus were relatively similar.Although the physical properties of the alloys did vary with composition, their differences are small on an absolute scale.This could also contribute to the similarities in the observed HAZ sizes across the four experimental alloys.In order for a clear difference in HAZ to be observed, the differences in melting temperature and physical properties would need to be more pronounced between alloys.
The discrepancies between the predicted and measured HAZ sizes could be due to a lack of microstructure considerations in the simulation.This could explain why the HAZ size is over predicted in the alloys with Sn-rich pro-eutectic phase (15Bi and 39Bi) and under predicted in alloy 78Bi, where the pro-eutectic phase is Bi-rich.The differences in thermal transport through the phases present in the alloys could influence the highly localized melting events that are thought to dictate the HAZ microstructure development.In  the case of 15Bi and 39 Bi, the Sn-rich pro-eutectic phase has a higher T m and greater thermal conductivity than the eutectic.Based on Fourier's Law, 22 there would be a higher flux of heat through the pro-eutectic phase than through the eutectic.Combined with the higher T m , the pro-eutectic phase would persist in the microstructure increasing the local thermal conductivity beyond what the continuum treatment of the substrate can account for.The result would be a more efficient conduction of heat away from the foil with less melting, which would produce a smaller HAZ.The opposite would be true for the 78Bi alloy, where more efficient transport of heat would occur through the eutectic and local melting at the pro-eutectic phase-eutectic interface would result in the under prediction of the HAZ size (Fig. 10).
Overall the model was effective at predicting the HAZ size based on the foil and substrate properties within the error of experimental measurement.It also allowed for the prediction of the temperature evolution as a function of time within the substrate to be predicted, which can be used to inform analysis of the microstructure evolution resulting from the RMF reaction.

V. CONCLUSIONS
To conclude, a predictive model was developed to evaluate the effects of RMFs on heating and modifying the microstructure of different substrate materials.The model leveraged a 1-D finite difference approximation of the heat equation in conjunction with phenomenological models for the physical properties of the substrate to predict heat flow and resultant HAZ depth.Experimental validation of the model was conducted using four Sn-Bi alloys with different thermal properties and initial microstructures, specifically different phase morphologies and distributions.The best agreement between the predicted and measured HAZ dimensions occurred when the variations due to foil distortion during the reaction were included in the simulation.Although close in value, a persistent discrepancy existed between the simulated and measured results in the 78Bi alloy.This could be due to a microstructural contribution to the heat flow that was not accounted for in the model.

FIG. 1 .
FIG. 1.A schematic representation of the data flow within the model.

FIG. 2 .
FIG. 2. Calculated Sn-Bi phase diagram using the assessed NIST_Solder database developed for the Pandat TM software package.

4 6 0. 1 FIG. 4 .
FIG. 4. DSC heating curves for each of the Sn-Bi alloys taken from the second heating cycle at a heating rate of 20 C/min.In this figure, the peaks represent endothermic reactions, specifically melting.Each of the curves has been normalized by the samples size to enable direct comparison.

FIG. 8 .
FIG. 8. (a) Comparison between predicted and measured values for the HAZ size.(b) Optical micrograph of 15Bi after the reaction of a foil showing the wavy nature of the foil post-reaction (dashed line represents the boundary of the HAZ).

FIG. 9 .
FIG.9.A graphical comparison of the predicted vs. measured HAZ sizes when the foil distortion is included in the calculation (6 1 standard deviation).

TABLE I .
Coefficients used for modeling density (kg/m 3 ) in both solid and liquid Sn-Bi alloys.

TABLE III .
24efficients used in Eq. (14) for the calculation of specific heat capacity (J/kg-K) in Sn and Bi in the solid and liquid states.24

TABLE V .
Experimental thermal conductivity values for the solid Sn-Bi alloys at different temperatures.