Comparing Relative Bond Characteristics Between Spherical and Elongated Morphologies for Cold Spray Process Using SPH Simulation

Under cold spray conditions, the modified Johnson–Cook model was adopted to perform single and multiple particle simulation for spherical and elongated aluminum alloy- Al-6061 feedstock particles. The splat formations were realistically presented; the temperature evolution throughout the deposition process stayed below the melting point of Al-6061, and the feedstock particles exhibited restitution for impact velocities lower than 200 m/s. Feedstock particles with elongated morphology experienced a lower elastic strain energy level than spherical morphology after impact, which implied the relative bond strength was higher for elongated particles than spherical particles. The displacement curves in single particle simulations for both morphologies suggested a spherical particle experienced a greater shock than the elongated particle upon impact. The relative bond strength achieved by multiple particle impact was lower than the single particle impact, even though the displacement curves showed the feedstock particles were individually embedded in the substrate.


Introduction
Cold spray is a surface coating technology that has become increasingly popular in metal 3D printing industries. It works by heating and accelerating process gas stream through a De Laval (diverging-converging) nozzle, which carries and propels the feedstock powders at velocities varying from Mach 1 to Mach 4 (Ref 1). This high-velocity powder spray plume is then directed at a substrate, so that particles collide onto the surface, causing them to plastically deform and interlock with the substrate surface. Throughout the impact, the temperature of the feedstock powders remains below the melting point, and therefore the entire bonding process is carried out in solid state (Ref 2).
To achieve sufficient bonding, the strain rates during impact must approach * 10 9 s À1 (Ref 3, 4), requiring extremely high particle impact velocities. During powdersubstrate impact, most of the powder kinetic energy is converted into plastic dissipation energy. The plastic dissipation energy is converted into heat energy rapidly within tens of nanoseconds timeframe (e.g., 40 ns for copper powder of 25 lm in diameter) ( Ref 5). The amount of heat generated within this timeframe experiences a lag in dissipation due to a significant difference between the rates at which heat is generated and dissipated. The rate at which heat is generated by deformation is much greater than the heat that is dissipated either to the surrounding environment or throughout the rest of the substrate. The unstable heat transfer consequently localized the heat energy to the specific point of contact, and this further promotes thermal softening. Since the localisation of heat energy happens parallel to the impact itself, the sudden decrease in the particle momentum during impact results in a significant shear mechanism that pushes the deforming particle outward to form a circular jet around the impact crater ( Ref 5).
The bonding process in cold spray could be compared to explosion welding or explosive powder compaction (Ref 4,6,7), for which a critical powder impact velocity value is defined to characterize the likelihood of successful bonding. The feedstock powders will bond with the substrate if the impact velocity is greater than the critical velocity, otherwise, they will bounce back ( Ref 8).
Researchers have been investigating the cold spray bonding mechanism via different models and approaches, exploring phenomena and theories such as adiabatic shear instability ( Ref 2,4,9,10,11,12,13,14), Although empirical data has partly supported different bonding models to a certain degree associated with their own mechanistic framework, none had so far been able to fundamentally explain the cold spray bonding mechanism. Each model is not mutually exclusive to the others, to say, several models could be adopted in explaining the cold spray bonding mechanism depending on the materials involved and setup conditions (see Table 1). However, since adiabatic shear instability criterion steers toward a pragmatic prediction on the critical velocity as a function of processing parameters and material properties for common cold spray metals, it is deemed the prevailing bonding model propounded in cold spray research community (Ref 8).
Although it is possible to study the bonding behavior through the performance of micro-ballistic tests (Ref 22), due to the short lived timescale (* 10 ns range) of deposition activity, measuring the overall stress/strain temporal evolution and the substrate temperature distribution upon impact remains challenging. Therefore, the use of numerical simulations to study cold spray bonding mechanism remains a more effective approach to tackle these challenges.
In addition to this, numerical simulations have already been used to build partial empirical models for predicting both the deposition behavior (Ref 3, 23) as well as mechanical properties such as porosity and residual stress ( Ref 9,24,25). In review, Yin et al. (Ref 26) had systematically examined the capability of cold spray numerical modeling using Lagrangian, Eulerian and Smoothed Particle Hydrodynamics (SPH) methods. It was found that when using the pure SPH method, the feedstock particles did not rebound after impact, which could be advantageous in capturing the dynamic features along the impact interface. A combination of SPH model for feedstock particles and Lagrangian model for substrate was also simulated by Yin et al. (Ref 26), showing a more significant crater on the substrate flat surface after impact, and with more intensive deformation observed. However, when using the combination of Lagrangian and SPH model to simulate cold spray impact process, a contact algorithm had to be implemented. Saleh et al. (Ref 27) had carried out neutron diffraction residual stress measurement and obtained the thickness stress profiles of cold spray coatings. Their experimental results were also successfully reproduced through the use of SPH modeling, which demonstrated the strength of SPH to be used as a predictive tool for the cold spray industry.
As cold spray technology is still a relatively new technique used in additive manufacturing, the existing deposition data is insufficient for engineers to construct a generalized predictive framework. Therefore, whenever a new feedstock material is to be cold sprayed onto substrates with dissimilar material properties, extensive testing would have to be conducted to validate the material combination, which is a costly and time-consuming process. To save time and cost, there is an increasing need for a simulation model that could accurately predict the bonding behavior based soley on the material properties of both feedstock particles and substrate, the environmental setup (initial conditions) and the geometrical information associated with the feedstock morphologies.

Material Deformation Model
The material model plays a crucial role in the simulation of the cold spray deposition process, therefore, the correct choice of an appropriate plasticity model is fundamental for ensuring the accuracy of such simulations. As the strain rates during impact could go up to * 10 9 s À1 ( Ref 3,4) in cold spray, the plasticity model of choice must be able to capture the deposition behavior throughout the extensive deformation process.

Standard Johnson-Cook
The Johnson-Cook plasticity model is one of the most widely used material models in the simulation studies of cold spray; it specified the strain rate sensitivity and the hardening behavior of the yield stress ( Ref 28). The original form of Johnson-Cook model had the following form: And, where A is the initial yield strength of the substrate material at reference strain rate. B, C, n, m are material constants which would require tuning for each material, r JC is the plastic flow stress, e p is the equivalent plastic strain (PEEQ), _ e p is the normalized equivalent plastic strain rate, _ e 0 is the reference strain rate, T Ã is the homologous temperature, T ref is the reference temperature, and T m is the melting temperature of the substrate.
This model is easy to implement, and material constants exist for many common metals. However, the original form of Johnson-Cook model does not capture the strain rate sensitivity at strain rates higher than 10 4 s À1 (Ref 29,30,31). This is because in the original Johnson-Cook model the flow stress is linearly dependent on the strain rate, and therefore at high strain rates it would produce unrealistic distortion in the results ( Ref 32,33).

Modified Johnson-Cook
To mitigate the limitation posed by the original form of Johnson-Cook model, several modified versions of Johnson-Cook plasticity models have previously been proposed, each tailored to address different applications. Lemiale et al. (Ref 32) investigated the combined effect of having the correct initial temperature and strain rate hardening by modeling copper/copper impact in SPH, and found that ignoring both the temperature and strain rate dependance resulted in poor predictions in terms of the splats formation profiles. When temperature was accounted for, but the sensitivity of strain rates remained unchanged, this had led to even poorer predictions. Therefore, to better account for the strain rate dependance, the standard Johnson-Cook model was adjusted to include the following expression for strain rates, _ e p greater than the critical value, _ e c (Which is $ 10 4 s À1 for copper).
where C is the material constant associated with the sensitivity of strain rate below the critical value, and D is that above the critical value. Tuazon et al. (Ref 34) studied the effectiveness of this modification by fitting the stress-strain curves of three highstrength steels (HSA800, Hi-Mn and AISI 4340) with bodycentered cubic crystal configuration into both the standard Johnson-Cook as well as a modified version which had the strain hardening term expanded under logarithmic power-law. The adjusted form of expression for the strain hardening term in the modified version of Johnson-Cook model was as follows.
where the material constant C governs the sensitivity of the strain rate and a is that in the power form. The yield strength of the materials were determined at different strain rates at 0.2% offset. As a result, the dependency between the yield strength of the materials and their respective logarithmic strain rates varied at different level of strain rates. As the relationship between plastic flow stress and the strain rate was non-linear, the original Johnson-Cook model did not fit well with the experimental data. In contrast to the standard Johnson-Cook model, the logarithmic power form of the modified Johnson-Cook model fits relatively well for all steels investigated in the study, with the downside that nonphysical stress was exhibited at low strain rates. In addition to Tuazon et al. (Ref 34)'s study, a quadratic form of modification made to the strain rate term was proposed by Huh and Kang (Ref 35), and when compared fits nearly as well. However, while the model fits well (down to * 5.3% error) with steel materials, when fit against copper the errors in fitting increased up to 41.4%. Although this modification more accurately maps the substantial increase in yield strain at high strain rates, it also exhibits an increase in yield stress at quasi static range, which is believed to be non-physical.
In this study, the modified Johnson-cook model from Chakrabarty and Song (Ref 36) was selected. Chakrabarty and Song (Ref 36) developed a modified power-law form of the Johnson-Cook model. The modification took viscous regimes into account, which was critical for accurate modeling of high strain-rate deformation. This modified model had the following form: where, And _ e c ¼ ys À1 . D became nonzero when the strain rate was equal to or greater than the critical strain rate determined by function _ e c y ð Þ. The parameters D and _ e c were determined by curve fitting from experimental data, which was performed using the MATLAB curve fitting tool. The other original material constants from the standard Johnson-Cook model (A, B, n, C, m and _ e 0 ) remained unchanged.

Numerical Simulation
In this study, both spherical aluminum particles and elongated aluminum particles were modeled as representations of powder with regular and irregular morphologies, respectively. The spherical particle was internally modeled in ABAQUS part module (Ref 37) with an 80 lm diameter, while the geometry of elongated particles were first parametrically modeled in Open-SCAD, and then imported into ABAQUS. To avoid mesh distortions during simulations, the smoothed particle hydrodynamics method was used in the numerical modeling. To ensure reasonable comparisons were able to be made, the maximum particle diameter for each morphology was maintained at 80 lm.
The impact simulation of single spherical aluminum powder was first simulated under the standard form of Johnson-Cook model, and then this was repeated using the modified Johnson-Cook model, to enable direct comparison of results. The modified Johnson-Cook model was then used exclusively throughout the rest of the simulations. Then, a set of simulations were conducted of single-particle powder impact on substrate were conducted, with a comparison carried out between single powder particles of both spherical and of elongated morphologies, impacting at velocities of 200, 400, 600, 700 and 800 m/s, in order to study the difference in deposition behaviors. Finally, multiple-particle simulations were conducted, with simulation of five elongated particles impacting simultaneously, and this was compared against the previous single-particle results.
The results would serve as the basis to allow comparison between the simulated deposition behaviors of single and multiple particle system. For each case simulation, the overall initial temperature was set to 25°C, and the size of the substrate was also constant.

Smoothed Particle Hydrodynamics (SPH)
The SPH method can handle extensive deformation problems such as the feedstock particles/substrate impact process in cold spray. The SPH elements are responsible for carrying the property data and the field variable such as pressure and temperature, and not to be confused with feedstock powders particles. Unlike standard grid-based elements that only serve as the interpolation data points, SPH elements propagate with the characteristics of an assigned material, they act similarly as real particles in space ( Ref 39).
In SPH simulation the resolution of the model constructed is governed by the amount of the SPH elements and their size. In comparison to the traditional tetrahedron mesh, the points in SPH are not connected by edges. At the beginning of the simulation mesh is still used as a means of distributing SPH elements, the mesh is not required after time evolution started. The meshless nature of SPH also helps prevent mesh distortion throughout the simulation.
The SPH field function is integrally represented by Kernel approximation as the following equation: where f is an arbitrary function of a positional vector x in 3D space, V is the volume, x À x 0 ð Þ is the length between the particle of interest x and a reference particle x 0 in V, h is the smoothing length and W is the Kernel function. Typically, the conventional continuum finite elements are defined first and then called to convert into SPH elements at the start of, or during the analysis.
The SPH approach is a total Lagrangian modeling method, which allows for the discretisation of an ordered group of continuum equations by interpolating the property data at a set of finite points discretely over the support domain without the need for structural mesh definition. The summation form of Kernel approximation function, which is mathematically more applicable in finite analysis could be expressed as the following equation.
where p j , m j , N and h are the density of particle j, mass of particle j, the total number of particles in the domain and the smoothing length, respectively (Ref 8).
The SPH method could be applied to standard materials in ABAQUS/Explicit, as well as for user defined materials, so long as the initial and boundary conditions are pre-specified as per standard Lagrangian model conventions. The SPH method is also applicable to a large range of applications, because it also allows for the contact interactions with multiple Lagrangian models. However, it is generally not as accurate as the Lagrangian finite element analysis under low deformation scenarios, or the coupled Eulerian-Lagrangian (CEL) analysis under higher deformation regimes ( Ref 40). Cubic polynomials were used as the SPH interpolator in this research as they are the default in ABAQUS/Explicit.

Setting Up Simulation Environment
Due to the nature of the deposition process, the actual number of elements are not always equal to the expected number of elements within the volume of each SPH particle. Therefore, a convergence study was carried out solely based on varying mesh density to optimize the selection of seed size for each feature. The smoothing length computed during the start of each analysis was kept at a level whereby the average number of particles associated with each element was between 50 and 60. As a result of the convergence study, the following subsections described the optimized global seed size as well as the partition requirement for each feedstock particle and substrate.

Spherical Particle
The 80 lm spherical particle was modeled in the ABAQUS CAE part module, cell-type partitioned by the 3 principal datum planes XY, YZ and XZ as shown in Fig. 1, and meshed at 3 Â10 À6 global size. 8-node linear brick elements were used with reduced integration and hourglass control.

Elongated Particle
The elongated particle shown in Fig. 2. was parametrically modeled in OpenSCAD. A degree 5 Bezier function was Fig. 1 The spherical particle is partitioned by XY, YZ and XZ principal datum planes for mesh control. The maximum particle diameter was kept at 80 lm Fig. 2 The elongated particle was meshed with ABAQUS default free technique. Due to its irregular surface, tri meshing on bounding faces is used where appropriate. The maximum particle diameter was kept at 80 lm used as the main curve to give a 2D profile for the rotatory extrusion, with the following form: where p 0 ; p 1 ; p 2 ; p 3 ; p 4 and p 5 are the coefficients. A closed 2D profile was obtained by joining the two end points of the Bezier curve by a polyline. Rotatory extrusion function was applied to compile a symmetrically elongated 3D model, and then a hull operation was used to combine two separate additional spheres with the main 3D model, in order to break up the symmetry. The elongated particles were then freely meshed with 4-node linear tetrahedron elements at 5 Â 10 À6 global size.

Substrate
The substrate was modeled in the ABAQUS CAE part module. It was modeled as a cylindrical plate with thickness t ¼ 2:5 Â d 0 , where d 0 was the maximum diameter of the particle in the single particle system, or that of the particle cluster in multiple particle system. To achieve an evenly distributed mesh a cylindrical profile, a smaller circular cross section was sketched on the surface of the substrate and facetype partitioned into two concentric cylinders by applying the extrude/sweep edges method on the sketched circular cross section. The 2 principal planes XY and YZ were used to partition the entire substrate. A global mesh size of 3.5 Â10 À6 was applied to the overall seeds and an element size of 3 Fig. 3 The cylindrical substrate is partitioned into two concentric cylinders, and principal datum planes XY, YZ and XZ for mesh control Â10 À6 was applied as the local seeds of the middle cylindrical section of the substrate. Encastre boundary condition was applied to the side and bottom face of the substrate so that only the top surface was subject to deformation (see Fig. 3).

Adaptation of Modified JC Model
The implementation of the modified JC model in this study was validated by simulating copper feedstock particle onto copper substrate under the same setup published in Chakrabarty and Song (Ref 36), to ensure that under the same conditions, this model provided consistent results. The feedstock particle and substrate were then changed to Al 6061 for subsequent simulations. Figure 4 showed the von Misses stress distribution of a simulation of an 80 lm Al 6061 feedstock particle impacting (at 700 m/s) on Al 6061 substrate with the standard JC model and modified JC model, after 200 ns. As expected (Ref 36), the jetting was less significant and the maximum von Misses stress was much greater (Increased by 44.2%) in the modified JC model.
As shown in Fig. 5, in both cases, the increase in temperature due to kinetic deformation was observed surrounding the particle impact zone, localized along the region of circular jet formation. The maximum temperature of the overall system after impact under modified JC model was greater than the increase observed by using the standard JC model. The rise in maximum temperature was due to plastic work being done under higher stress conditions. However, the maximum temperature in both cases remained below the melting point of Al6061.The maximum equivalent plastic strain under modified JC model was 14.5% lower than that under standard JC model. The increase in stress experienced and lower strain observed via the modified JC model suggested more energy was  converted into heat energy during deposition than that under the standard JC model. This was supported by the increase of maximum temperature observed in Fig. 6. As described, the modified JC model should achieve a lower equivalent plastic strain, which implied lower deformation at this level of strain rate. Therefore, the circular jet was more realistically represented in the case using the modified JC model than that in the standard JC model, as shown in Fig. 4. In addition to the difference in the extent of deformation between the two cases, the splat flattening ratio in each case was also calculated for comparison. As seen in Fig. 7, the length and height of splats were measured in the ABAQUS viewer by distance queries, and these were recorded in Table 2. Evidently, taking the viscous regime into account in the modified JC model had resulted a 33.1% decrease in the flattening ratio.
The Chakrabarty and Song (Ref 36)'s modified Johnson-Cook model adopted in this paper had a simple mathematical representation and only required two additional parameters to accurately capture the strain rate sensitivity. As expected, it was found the temperature stayed below the melting point of Al6061 for all cases regardless of morphology type and the number of feedstock particles involved.

Spherical Morphology Versus Elongated Morphology
Simulations were conducted under the modified Johnson-Cook model for both the spherical and elongated particles at 200, 400, 600, 700 (Critical velocity) and 800 m/s. Results suggested that as the impact velocity decreased, the likelihood of particles bouncing back increased. Displacement along the direction of impact was used to study the adherence of feedstock particle onto substrate surface. Figure 8 shows that at 200 m/s impact velocity both spherical and elongated particles bounced back and beyond their original location. At 700 m/s, both spherical and  elongated particles eventually converged within 8 to 10 micron range in displacement and stayed embedded in the substrate after 150 ns. However, it should be noted here that this simulation model only considered the plasticity model. Therefore, it cannot definitively predict whether bonding occurred. Nevertheless, future work could be extended to include a cohesive zone model to comprehensively explain the bonding properties. The elastic strain energy of the particle after impact signified the tendency of separation between particle and substrate, therefore studies had suggested the use of recoverable strain energy to assess the relative bond strength between feedstock particles and substrate (Ref 10,41,42).
While the impact kinetic energy plays a critical role in cold spray adhesion in terms of the total energy landscape, the work of Bae et al. (Ref 10) was used as the foundation in this paper. The terminology was drilled down to the level where recoverable strain energy became apparent for both cases of simulations, so that relative bond strength in each case could be indicatively compared.
In contrast, the spherical particle exhibited higher strain energy than that in elongated particle, which implied higher tendency to separate from substrate after impact. However, Fig. 9 shows the impact velocity in each case also affected the level of strain energy.
In each of the simulated morphologies, since the strain energy level at 200 m/s was much lower than that at 800 m/s, when the impact velocity was below the critical velocity the particles still bounced back. For the case of impact velocity lower than the critical velocity, the restitution of feedstock particle resulted in a low change in internal energy, E U upon impact, and this was translated into lower strain energy level, as the energy contributed toward plastic deformation, E p became insignificant in comparison to the recoverable strain energy, E r .
As the energy dissipated due to viscous effects, E v was insignificant in comparison to the recoverable strain energy, E r (Ref 10), it was ignored in this case.
Whereas, for the cases where impact velocity was higher than the critical velocity, the high recoverable strain energy resulted in the breakup of materials, and produced visible splashing upon impact, as shown in Fig. 10(a) and 11(a).
At this point, it is important to note that the boundary reflected waves as a result of insufficiently large substrate size could have an effect on the bonding predictions. In this work, a compromise was required to ensure acceptable computing time, therefore the depth of the substrate used in this study was not large enough to completely prevent the bottom boundary reflection. However, since this study was only intended to indicatively study the effects of deposition behavior in different morphologies, the findings remained unaffected so long as the same substrate size was used for all morphology cases.

Single Particle Simulation Versus Multiple Particle Simulation
To study the difference in terms of strain energy level in multiple particle system in comparison to that in single particle system, the impact of 5 elongated particles and 5 spherical particles were simulated to represent multiple particle systems. Figure 12 shows the position of individual particles in each case. In addition, all elongated particles were oriented identically as in the single particle simulation, with the direction of maximum particle diameter parallel to the direction of impact. All particles were initially 2 lm above the impact interface.
For each morphology case, 5 particles were simulated to hit the Al6061 substrate at 700 m=s simultaneously under the modified Johnson-Cook model. The corresponding strain energy evolution curves in each morphology Fig. 9 The evolution of strain energy level of the spherical particle and elongated particle at different impact velocities presented in Fig. 13. showed the results simulated under multiple particle condition could be 5 times higher than that under single particle condition. As shown in Fig. 14, the displacement curve of individual particle for each morphology were clustered along the same trend. This implied that for multiple particle simulation, the displacement was not significantly affected by the type of morphology associated. Displacement curves showed the evolution of the location of SPH elements as a function of time before and after the deposition process. For feedstock particles that did not deposit upon impact, the displacement curves tended toward zero along the y (displacement) axis. The corresponding stress profile for each morphology during multiple particle impact can be seen in Fig. 15.

Conclusion
When simulating cold spray deposition using the Chakrabarty and Song (Ref 36)'s modified JC model, phenomena were reproduced supporting the notion of critical velocity. Instead of depositing as splats, particles impacting below the critical velocity are expected to bounce back. In simulation, feedstock particles were observed to bounce off the substrate for impact velocities lower than 200 m/s. Feedstock particles with spherical morphology experienced higher elastic strain energy than the elongated morphology, which implied the likelihood for a rebound was greater in the case of spherical morphology. However, it was also shown that the impact velocity could affect the level of elastic strain energy experienced by the embedded particle as well. Therefore, both the impact velocity and the elastic strain energy would have to be considered in deducing the relative bond strength for each case.
The number of feedstock particles also affected the accuracy of a simulation result. Therefore, in determining the relative bond strength among different morphology types, the results generated from multiple particle simulation would give a more realistic representation than from single particle simulations. Due to the nature of cold spray impact process, the SPH elements distribution in the feedstock and substrate are not fully ordered throughout the impact process, i.e., the actual number of elements do not always equal to the expected number of elements within the volume of each SPH particle, throughout the impact. It can be challenging to decide the level of convergence one is after solely based on varying smoothing lengths. Experimental data are required to fine tune simulation model to a point one knows what smoothing length is appropriate. Given this paper is purely simulation, the smoothing error was controlled by fixing a relatively small smoothing length (By default, large enough to scan 50-60 elements), and performing convergence study based on varying mesh density. The idea was to keep the smoothing length fixed for all simulations, and use the result of mesh density convergence study to guide the choice of seed size for each feature (feedstock/substrate), i.e., fixed distance but more neighboring elements. The following was the result of our mesh density convergence study (see Fig. 16): Note: The seed size governs the distance between each nodes, the smaller the seed size the closer the distance between neighboring nodes, i.e., greater number of elements.  Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.
Funding Open Access funding enabled and organized by CAUL and its Member Institutions.