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) with temperature and composition in the Sn-Bi system were utilized using literature data.
I. INTRODUCTION
Reactive multilayer foils (RMF) have garnered significant attention because of their potential use as localized high intensity heat sources.1–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 self-sustaining 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.8 This 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.9
The majority of the literature available for this class of materials is largely focused on elucidating the physics of the foil reaction and how structure-property relationships of the foil affect the energy output.10–20 By comparison, fewer investigations have been undertaken to develop applications for the RMFs.1–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.21 Through 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
1. 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
where
In Eq. (1), the variables T, x, α, and t are the temperature, distance, thermal diffusivity, and time, respectively. α is comprised of the thermal conductivity (k) in W/m-K, specific heat capacity (Cp) in J/kg-K, and density (ρ) in kg/m3. 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 (Tsolidus) and liquidus (Tliquidus) 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 Δt and Δx 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 Δt and Δx were taken to be 1 μs and 10 μm, respectively. These values were chosen in order to balance the spatial resolution of the output HAZ size, calculation stability, and calculation time.
2. 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 (ΔH) is based on the enthalpy of formation for the compound being formed during the reaction (ΔH0), the bilayer thickness (λ) and the thickness of the intermixed region between two adjoining layers (w). In this work, ΔH0 was taken to equal 1383 J/g assuming the formation of an equiatomic NiAl compound24
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.25 Energy transport across the foil-substrate 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 ρf, xf, and ks 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 Cp value used in the calculations depending on temperature. This is shown in Eq. (8),26 where the enthalpy of fusion (ΔHfus) and the melting range (ΔTm) 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 Tliquidus 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.27 The Redlich-Kister polynomial (Eq. (10)) has been used extensively to model temperature and composition dependent properties for solution phases.28–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 P1 and P2 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, X1 and X2. The deviations from ideality are accounted for with the interaction parameters (Li) in the Pexcess 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.
1. Density
A relationship for the variation of volume as a function of temperature is well established25 (Eq. (11)). In the equation below, α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 ρ0 is the density at a reference temperature T0
Literature data were used to fit the relationship between temperature and density (Eq. (12)) for the pure elements24,31–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,36 Table I contains the coefficients for the models of density for the Sn-Bi system in the solid and liquid states.
. | ρ0 Bi . | T0 Bi . | αCTE Bi . | ρ0 Sn . | T0 Sn . | αCTE Sn . | L0 . | L1 . | R2-value . |
---|---|---|---|---|---|---|---|---|---|
Solid Sn-Bi | 9790 | 300 | 3.04 × 10−5 | 7260 | 300 | 5.43 × 10−5 | 1054.6 + 1.2 T | 2108.7 − 0.3 T | 0.958 |
Liquid Sn-Bi | 10050 | 544 | 1.25 × 10−4 | 6980 | 505 | 1.01 × 10−4 | 2409.1 − 3.0 T | 2011.7 + 0.3 T | 0.988 |
. | ρ0 Bi . | T0 Bi . | αCTE Bi . | ρ0 Sn . | T0 Sn . | αCTE Sn . | L0 . | L1 . | R2-value . |
---|---|---|---|---|---|---|---|---|---|
Solid Sn-Bi | 9790 | 300 | 3.04 × 10−5 | 7260 | 300 | 5.43 × 10−5 | 1054.6 + 1.2 T | 2108.7 − 0.3 T | 0.958 |
Liquid Sn-Bi | 10050 | 544 | 1.25 × 10−4 | 6980 | 505 | 1.01 × 10−4 | 2409.1 − 3.0 T | 2011.7 + 0.3 T | 0.988 |
2. Thermal conductivity
It is known that the variation of thermal conductivity in a pure metal with temperature can be modeled using Eq. (13).22 Data taken from literature were used to account for this temperature dependence for both Sn24,37–39 and Bi24,38,40
Once this behavior was modeled, additional literature data38,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 R2-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.
3. Heat capacity
Cp as a function of temperature (Eq. (14)) was modeled using the standard equation for molar heat capacity41 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 Cp as a function of composition can be well modeled by using an ideal solution model42,43 (i.e., where the excess term is equal to zero). Therefore, this was used for calculating Cp 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 μm. 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-0744 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.0 mm 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.45 Then along with experimentally measured Cp 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 ΔHfus 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 ΔHfus 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 ΔHfus 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 focused-ion 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 μm thick and had a 1 μm thick layer of InCuSil®-ABA™ (Ag59Cu27.25In12.5Ti1.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 μm 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.46 Once 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,48 For 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).
IV. RESULTS AND DISCUSSION
A. Experimentally determined model parameters
The compositions of the four alloys determined using XRF are shown in Table IV.
. | . | Composition . | |
---|---|---|---|
Alloy designation . | Sn . | wt. % Bi . | at. % Bi . |
15Bi | Bal. | 15.0 (14.5 ± 0.3) | 9.1 (8.8 ± 0.2) |
39Bi | Bal. | 39.0 (34.5 ± 0.4) | 26.6 (23.1 ± 0.3) |
57Bi | Bal. | 57.0 (52.0 ± 1.2) | 43.0 (38.1 ± 1.2) |
78Bi | Bal. | 78.0 (71.0 ± 0.8) | 66.8 (58.2 ± 0.9) |
. | . | Composition . | |
---|---|---|---|
Alloy designation . | Sn . | wt. % Bi . | at. % Bi . |
15Bi | Bal. | 15.0 (14.5 ± 0.3) | 9.1 (8.8 ± 0.2) |
39Bi | Bal. | 39.0 (34.5 ± 0.4) | 26.6 (23.1 ± 0.3) |
57Bi | Bal. | 57.0 (52.0 ± 1.2) | 43.0 (38.1 ± 1.2) |
78Bi | Bal. | 78.0 (71.0 ± 0.8) | 66.8 (58.2 ± 0.9) |
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.38 The 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.
Alloy . | Temperature ( °C) . | Thermal conductivity (W/m-K) . |
---|---|---|
15Bi | 36 | 32.6 ± 1.1 |
56 | 32.0 ± 0.5 | |
76 | 32.1 ± 0.3 | |
39Bi | 37 | 30.9 ± 0.4 |
56 | 29.7 ± 0.2 | |
77 | 28.4 ± 0.6 | |
57Bi | 37 | 14.1 ± 0.1 |
56 | 14.2 ± 0.1 | |
77 | 13.5 ± 0.1 | |
78Bi | 37 | 9.1 ± 0.2 |
56 | 9.0 ± 0.04 | |
78 | 8.4 ± 0.1 |
Alloy . | Temperature ( °C) . | Thermal conductivity (W/m-K) . |
---|---|---|
15Bi | 36 | 32.6 ± 1.1 |
56 | 32.0 ± 0.5 | |
76 | 32.1 ± 0.3 | |
39Bi | 37 | 30.9 ± 0.4 |
56 | 29.7 ± 0.2 | |
77 | 28.4 ± 0.6 | |
57Bi | 37 | 14.1 ± 0.1 |
56 | 14.2 ± 0.1 | |
77 | 13.5 ± 0.1 | |
78Bi | 37 | 9.1 ± 0.2 |
56 | 9.0 ± 0.04 | |
78 | 8.4 ± 0.1 |
The values for the transition temperatures (Tsolidus and Tliquidus) and ΔHfus 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.21 The 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 λ for the Ni-Al RMFs was determined to be 47.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.24 As 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.49 As 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.50 This 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®-ABA™ 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 Tliquidus 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 Tm 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 Tm, 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.
ACKNOWLEDGMENTS
This work was support by the Sandia Campus Executive Fellowship. Sandia National Laboratories is a multi-program laboratory managed and operated by Sandia Corporation, a wholly owned subsidiary of Lockheed Martin Company, for the United States Department of Energy's National Nuclear Security Administration under Contract No. DE-AC04-94AL85000. P. M. Johns was supported by the U.S. DoE Nuclear Energy University Program Fellowship and C. G. Davis was supported by the U.S. NSF Graduate Research Fellowship Program (No. DGE-082270).