Development of palladium nanoparticles deposition on a copper substrate using a molecular dynamic (MD) simulation: a cold gas dynamic spray process

. The objective of this study is to create an ultra-thin palladium foil with a molecular dynamic (MD) simulation technique on a copper substrate surface. The layer formed onto the surface consists of a singular 3D palladium (Pd) nanoparticle structure which, by the cold gas dynamic spray (CGDS) technique, is especially incorporated into the low-cost copper substrate. Pd and Cu have been chosen for their possible hydrogen separation technology applications. The nanoparticles were deposited to the substrate surface with an initial velocity ranging from 500 to 1500 m/s. The particle radius was 1 to 4nm and an angle of impact of 90 ° at room temperature of 300 K, in order to evaluate changes in the conduct of deformation caused by effects of size. The deformation mechanisms study revealed that the particle and substrate interface is subject to the interfacial jet formation and adiabatic softening resulting in a uniform layering. However, shear instabilities at high impact speeds were con ﬁ rmed by the evolution of von Mises shear strain, temperature evolution and plastic strain. The results of this study can be used to further our existing knowledge in the complex spraying processes of cold gas dynamic spray technology.


Introduction
Hydrogen is an ideal energy carrier due to its low environmental impact combustion enthalpy. Hydrogen fuel can be used in aircraft, vehicles, ships, or other machinery directly. Although many issues are to be overcome in the use of hydrogen technology, hydrogen is generally believed to become the main energy carrier in this century as part of the Fourth Industrial Revolution (4IR) [1]. To improve the pureness of hydrogen, consideration is given to the separation membrane as the ultimate and effective technologies. This includes dense metal membranes, ion conductive membranes, and microfiltration membranes [2]. Of these options, dense metal membranes produce the purest stream of hydrogen [3]. Among pure metal membranes, palladium is the best candidate since it has excellent hydrogen dissociation capabilities and sufficiently highly hydrogen diffusivity [4][5][6]. The membranes from palladium and palladium alloys were the subjects of past research [3,[7][8][9]. Except for chemical vapour deposition and sputtering deposition [10][11][12], which are normally used as separation technique, a new surface treatment technique, the cold gas dynamic spraying (CGDS), can also be used to produce the Pd membrane. CGDS is preferable because this advanced technology ensures the generation of layers without unnecessary particle heating which preserves the original physical and chemical features. The velocity of the particles is increased to supersonic velocity by a high-speed gas stream through CGDS. Nevertheless, the characteristic interaction of the deposited particle with factors that influences the deposition process has difficulties, including the elimination of defects in the deposited substrate, as well as the preservation of chemical and mechanical stability in working conditions. Chemical stability is required to reduce surface fouling and poisoning to avoid inter distributions of elements between the deposited layer and substrate. In order to reduce stresses caused by the different heat expansion factor of the deposited layer and the substrate, mechanical stability is also required [13][14][15][16][17]. The performance of the membrane produced by cold gas dynamic can be enhanced by studying the process of Pd nanoparticles deposition on a metal substrate.
The computational modelling can be an alternative means of analysing the deposition process compared to the larger experimental studies. Typically, the molecular dynamics (MD) approach was used to examine the transient phenomena occurring the process of impact and to investigate the microscopic development of substrate film growth [18][19][20][21][22]. Some studies also documented single particle impact simulation using MD. The impact behaviour of ductile materials such as copper and aluminium after deposition was investigated by Oyinbo and Jen [23] using the effects of various processing variables such as preheating temperature, initial impact velocity, coefficients of friction, and different material combinations. They discovered that the increase in initial kinetic energy increases with increasing substrate temperature and equivalent plastic strain. Au and Cu particles were modelled by Gao et al. [24,25] with a couple of hundred atoms to study the influence of particle deposition variables in cold spray, they found that higher temperature, initial impact velocity, and larger particle size increased the particle/substrate overall bonding strength. Malama et al. [26] modelled the Ti and Ni impact on a Ti substrate, they discovered that higher bond strength is achieved by increasing particle size. Their study examined about 2000 atoms for the analysis. The effects of spray parameters on small clusters of up to 400 atoms as model systems were studied by Joshi and James and it was demonstrated that particles seemed to rebound out of the substrate at impact angles exceeding 90 degrees [24,27,28]. Daneshian and Assadi [29] modelled on the interatomic impact of constitutively ductile nanomaterials and identified the particle size effect and speed of impact in the fragmentation, deformation, and rebounding of ductile particles.
Experimental and numerical investigations were conducted by Del Popolo et al. [30,31] with Pd clusters and Au substrate. The orientation of the substrate crystal was taken to be in [111] direction. They discovered that layer by layer dissolution of the clusters does not occur but rather at the Pd cluster edges. Another investigation was carried out by Guo et al. [32] on the Pd islands growth on the surface of Nickel substrate with orientation [111] direction. They found out that with 0.4 monolayers or less of the Pd coverage, 2-dimensional mechanism of the growth appeared while 3-dimensional islands occur at the higher monolayer. The stability and structure of Pd particles were also studied by Rojas et al. [33] on the Pt [h k l ] and Au [h k l ] surfaces at different degrees and discovered that along the crystallographic orientation of the substrate, the Pd adlayers grew pseudo morphically and epitaxially.
Therefore, this study tailors the effect of many variables to describe the atomistic interaction between Pd nanoparticles and Cu substrate using molecular dynamics. Our research focuses on the fundamental systems of plasticity and metal hardening for the characterization of Pd nanoparticle deformation behaviour at a fundamental level. The particle/substrate interface is also taken into account, which leads to a further deformation mechanism due to the high localized temperature.

Molecular dynamic simulation
In order to carry out a series of MD simulations in this study, the CGDS method is simulated using a "Large-scale atomic/molecular massively parallel simulator" (LAMMPS) [34]. For visualizing and analysing atomistic simulation data, open visualization tool (OVITO) [35] together with the algorithm of dislocation extraction (DXA) [36][37][38] were adopted. The research considers the effect of nanoparticles on a three-dimensional (3D) metal substrate surface. The simulation system includes both palladium nanoparticles (Pd) and copper (Cu) atoms as the substrate. The Pd nanoparticles consist of 2,310 atoms while the Cu substrate with fcc lattice structure consists of 77,326 atoms. The substrate is considered rectangular, single-crystal block with dimensions of 15 nm Â 15 nm Â 4.0 nm and lattice spacing of 0.361 nm. The Pd nanoparticles are spherical with a lattice spacing of 0.389 nm arranged in fcc lattice structure. Figure 1 shows the simulation model used for the CGDS deposition process and the simulation conditions for the analysis are presented in Table 1. The model atomic interaction is described by the Embedded Atom Method (EAM) potential function [39,40] as described below Here, the atoms in the model are labelled i and j where (j ≠ i). r ij and r b are the distance and electron density respectively between the atom i and j in the solid. The density, r i is the summation of individual atomic densities, f (r ij )with respect to atom i. The potential function is denoted by Fab, and the embedding function, F is the energy that is needed to move atom i into the electron cloud of type a.
Periodic boundary conditions were used in the x and y directions and non-periodic boundary conditions in the z-direction. The Newton motion equation was integrated with the Verlet algorithm for each atom [41,42]. 1fs integration time-step was adopted. To prevent interaction between the nanoparticle and substrate, the particle was positioned 1.5 nm over the substrate surface. Nosé-Hoover thermostat [43] were used before the actual deposition of the particle, to both equilibrate the substrate and the particle for 10 ps with NVT ensemble at 300 K. The initial particle velocity of 900 m/s was set after equilibration in the z-direction. For adhesion to occur, this initial velocity should be greater than the critical velocity of copper and below particle erosion velocity [44]. The study carried out by Goel et al. [45] with varying initial temperature revealed that the velocity used in this analysis is less than the velocity of which particle melting is anticipated to occur. This CGDS deposition process takes place in the NVE ensemble. The thermal coupling was removed, and the entire system was assumed to be thermally insulated for 20 ps. The MD simulation of particle impact was assumed to be an adiabatic process in this study.

Palladium nanoparticle's deposition morphology on the substrate at different time periods
First, the morphology of the deposited powder was carefully studied and the Pallaidum nanoparticle morphological evolution deposited at 900 m/s initial impact velocity on the Cu substrate were shown in Figure 2 small cross-sections of the system were shown in the figures and the Pd nanoparticle atoms are showed in red balls and the blue atoms are the Cu substrate atoms. The entire process of deposition is of three phases: absorption, impact and relaxation according to the evolution of the morphology. At absorption, the Pd nanoparticle with an initial velocity accelerate towards the surface of the substrate from its original location slowly and retained almost the same spherical structure. At impact, the nanoparticle/substrate collision occurs at the interfacial zone. The Pd atoms structure changes as some of the atoms spread on the surface of the substrate and some penetrated into the substrate. Meanwhile, some atoms of the substrate were forced to form jet on the interface and buildup around the particle. The interpenetration of the nanoparticle/ substrate atoms on the surface occurs at this stage. Both the nanoparticle and substrate at the interfacial zone deform plastically. Some atoms at the impacted region of the substrate rebound and the nanoparticle were rearranged. Full expansion of the system was then achieved, and the stage of relaxation came. The energy minimization was performed to aid the atoms structural rearrangement. As shown in Figure 3, the substrate lattice arrangement was followed by the reconstruction of the particle structure after 10 ps. The kinetic energy of the system assumed to approach a minimum stable value. The kinetic energy of the system first increased to a maximum value at around 2 ps then decreased and later experienced a stable minimum value. At the initial stage, the nanoparticle retained the initial crystal structure, then the atoms at the lower part of the particle restructure after impact follow the arrangement of the substrate atoms, and eventually, the whole particle reconstructed from the lower part to the upper part. The degree of reconstruction depends on the incident energy of the nanoparticle. The rearrangement of the deposited Pd coating lattice structure and the atoms of the substrate is advantageous for the hydrogen separation membrane to increase the membrane efficiency by  substantially reducing hydrogen atoms resistance and the diffusion path inside the membrane because of the lattice mismatch.
After the deposition, the entire system tended towards a balanced state in the course of the deposition. The nanoparticle morphological evolution was accompanied by the structural development observed at certain points by the nanoparticle RDF (radial distribution functions) as shown in Figure 4. During the first phase of deposition, there were small projections in the RDFs, showing the short-ranged and long-ranged order of the crystal-structural properties. The cluster crystallinity subsequently weakened significantly, and characteristics liquid-like nature emerged. The fluctuation of RDFs which tended to 1 showed the characteristics of the long-ranged disorder and short-ranged order. On the relaxation point,  in the RDFs, some small peaks took place and the crystal structure again emerged. Therefore, the nanoparticle transforms to liquid-like state from the crystal structure state then evolved back into its crystal form after the impact. Figure 5 shows cross-sections of the splat after it has solidified at the end of the simulation. The splats topography is specifically illustrated by the figures. The dark line in the cross-sections reflects the substrate's top surface, and below that line reflect splat propagation in the substrate. To this extent, the cross-section shows, therefore, the depth of penetration, flattening, and subsurface penetration dimensions simulated at different temperatures and velocities, and the maximum splat diameter and height are presented in Table 2. For clarification of measurements, slicing the x-x section of the coated surface and disengage substrate atoms was used to obtain cross-sections of the splats. The x-x section was named for all particles.
The influence of impact velocity is a determining factor for the splat topography as shown in Figure 5. With a 500 m/s velocity and at 500 K initial temperature, the topography is significantly hemispheric, pyramidal at 1000 m/s, and flat over the surface of the substrate like a mushroom at 1500 m/s. In addition, the splat diameter decreases with an increase in the velocity of impact, and so does the degree of splat dissemination in the sub-surface. The flattening diameter of the nanoparticles was lower at preheating temperature of 700 K compared with 500 K at 500 m/s initial velocity. The liquidity behaviour in relation to temperature is therefore nonlinear. The liquidity has improved considerably at higher velocities: the diameter of the splat is gradually flattened at higher velocities (1000 m/s À1500 m/s). The propagation of subsurface also appears to be decreasing with increasing particle temperature. The height ratio of splat below and above the impact surface is 83.2%, 73.6%, and 63.4% for 300, 500, and 700 K preheating temperatures respectively at 1500 m/s impact velocity.
According to the investigation carried out by Escure et al. [46], at temperature of 300 K, the diameter of aluminium splats was less than the splats obtained at 1000 K. This phenomenon has only occurred at higher impact velocities in this current study, i.e. splat diameter increased at 1000 m/s to1500 m/s higher impact velocities but was reduced at relatively lower velocity of 500 m/s when the preheating temperature increases. There is a need for an appropriate and higher initial impact velocity to overcome the inertial forces that obstruct the spreading of the nanoparticle on the substrate. The particle kinetic energy is converted into three components of work-surface energy, viscous deformation and substrate in-depth propagation. The strength of the intermediate contact is also dependent on the pressure of the impact, which depends greatly on the surface of contact, and therefore the quality of the splat.

Flattening ratios
In this work, flattening aspect ratio and flattening diameter ratio [47] were utilised as a quantitative indicator and their variation with impact velocities and preheating  temperatures are shown in Figure 6. Flattening aspect ratio is defined as the ratio of post-impact particle maximum diameter to the splat total height after impact. Whereas, flattening diameter ratio is the ratio of the post-impact maximum particle diameter to the original diameter of the particle. In the earlier stages of thermal coating work, a theoretical analysis of the flattening mechanism was conducted in order to quantitatively understand the process of flattening. Nevertheless, the degree of subsurface substrate penetration due to the impacting nanoparticle is not considered. Thus, in this study, flattening aspect ratio is implemented as a better parameter for measuring both horizontal and vertical deformations. From Figure 6, it is clear from that flattening aspect ratio is consistent at temperatures up to 500 K in comparison to the flattening diameter ratio and at velocities below 1000 m/s. Well, that is because the overall deformation of the system is considered in flattening aspect ratio, and the effect of the vertical deformation is not discarded unlike the flattening diameter ratio. Since the degree of the splat subsurface propagation has a greater influence on the overall flatness, the measurements of splat topography are in line with the explanation above. It is recommended therefore that a more appropriate measure of the flattening aspect ratio can be used to assess the overall deformation mechanisms than the flattening diameter ratio. This parameter changed quickly over and above 1000 m/s to 1500 m/s impact velocity.

Evolution of von Mises shear strain
To examine the recrystallization behaviour and to describe the extent of deviatoric deformation, it is important to evaluate the value of a single characteristic strain, von Mises shear strain was used to calculate this value. The atoms plastic deformation was quantify using Shimizu et al. [48] proposed algorithms. A comparison of two atomic configurations (start configuration and during impact) using OVITO was carried out for each atom i to calculate the atomic local shear strain (von Mises strain). The first approach was to calculate the local lagrangian strain matrix, h i ¼ 1 2 JJ T i À I À Á using a local deformation matrix J i , von Mises shear strain was then computed for each atoms i. The von Mises shear strain h mises i is a good measure of inelastic local deformation and is shown in equation (3) h mises Using this algorithm at different initial impact velocities, the von Mises shear strains is shown in Figure 7 at 500 K initial temperature. Figure 9 indicates that the highest shear strain is centred directly below the subsurface with higher impact velocities while the average strain is apparent at lower velocities at the interface or at peripheral atoms. This illustrates how the material's recrystallization patterns are affected by different impacting velocities. This also explains why the residual strains are located on the subsurface of the substrate resulting from cold spray systems at higher spray rates. Residual strain caused by such phenomena was shown to cause an orthopaedic implant to fail. The non-destructive neutron diffraction method was also proposed [49] to Fig. 6. Variation in (a) flattening diameter ratio (b) flattening aspect ratio with respect to particle impact velocity at preheating temperatures of 500 K, 700 K and 1000 K. characterize the substrate residual strains. MD is therefore useful for studying these sub-surface phenomena that require complex experimental setups otherwise in order to examine these details.

Effect of the velocity of impact on palladium nanoparticle deposition
In accordance with the morphological evolution, the Pd nanoparticle temperature evolution after impacted on Cu substrate at a different initial velocity of the incident was shown in Figure 8a. As nanoparticles were adsorbed on the surface of the substrate, the rise in the binding energy led to increasing the temperature of the nanoparticles. The particle kinetic energy was converted to thermodynamic energy of the substrate during the collision process and thus the temperature of the particle dropped rapidly and continuing approximately for about 2-6 ps. The particle temperature with high impacted kinetic energy decreases rapidly than the lower cases. During relaxation, the temperature of the particle declined in the form of a damping oscillation and approached the thermal temperature of the substrate-temperature-control layer. The nanoparticle with high initial impact velocity at each point of the deposition cycle took less time and held higher temperature until equilibrium because of the higher initial energy than the low initial impact velocity.

Particle size effect on palladium nanoparticle deposition
Another important parameter for optimal cold gas dynamic spray (CGDS) coatings for separation membrane is the size of the deposited particle. Figure 9 demonstrates the result obtained by varying particle sizes with an initial velocity of 500 m/s at an impact angle of 90°. The result obtained from this analysis shows that the particle size should be no less than 1 nm in radius, as the substrate cannot be coated uniformly. At 2 nm, thick and uniform coatings were observed. The deposition height increases as the particle size increase to 2 nm, as shown in Figure 9b. Above 2 nm particle size, no significant differences are observed in the deposition indicating that optimal coating for conditions used in this analysis is obtained with particle sizes varying in the range of 1 nm to 2 nm. Figure 10 indicates the changes in the temperature and plastic strain S.T. Oyinbo et al.: Manufacturing Rev. 7, 29 (2020) of the nanoparticle with respect to time after impact at the interfacial zone. The nanoparticle with higher particle size at each point of the deposition cycle held a higher temperature until equilibrium because of the higher initial energy and larger contact particle/substrate surface area than the particle with a small radius. Also, the plastic strain increases by increasing the particle size as the deformation progress overtime in the interfacial zone, leading to the softer and plastic deformation of the Pd nanoparticle in the contact region because of wide contact areas between the particle and the substrate material, varying from plastic nature to viscous flow at deformation Fig. 9. Cross-section view of particle/substrate impact at particle size of (a) 1 nm, (b) 2 nm, (c) 3 nm, (d) 4 nm diameter at initial temperature of 300 K. zone. This phenomena reflects the thermal relaxation in this analysis and can be attributed to the parameters in the Johnson-Cook model [50].

Comparison with previous experimental and numerical results
The results of this study of molecular dynamic (MD) simulation comply with previously published experimental and numerical results [14,[51][52][53][54]. In the MD report, the flattening increases with the increase in impact speed. The flattening ratio increases and then starts to decrease, indicating an optimal initial velocity for high-quality deposition. There is a similarity with Assadi et al. [14] experimental observations, where particles of copper are deposited with the cold spray process on the copper substrate. In addition, the outcome of the MD simulation is consistent with previously reported experimental results on von Mises shear strain variations [14]. Jodoin et al. [52] experimental study show that the effect of particle morphology on particle velocities is more pronounced for the larger range of particle size. According to Goldbaum et al. [53], the adhesive force of splats varies in the deposition pass with their particle size and particle location. The particles deposited at the centre of the deposition pass displayed a greater adhesion relative to the particles on the boundary of the deposition pass. These results are in agreement with the present MD simulation analysis. Also, there is a similarity between the results of the MD simulation and the findings of Ghelichi et al. on the numerical analysis for the cold spray Cu/Cu coating process [55]. The estimates of MD simulation of the time variance of the plastic strain during the cold spray process are consistent with the results of the Grujicic et al. [56] and Joshi et al. [27] FE simulation research. In both analyses, the strains for initial velocities greater than critical velocity are unexpectedly subjected to high shear strain-rate deformation and material softening. The variance in the strains values in relation to the approaching distance in this study showed similar trends in the MD simulation and findings of thermo-mechanical finite element simulation of the cold gas dynamic spray process [57].

Conclusions
A detailed analysis of the deposition, morphology, energy, and surface characteristics of the Pd nanoparticle at the specific initial velocity of 900 m/s was performed at room temperature by MD simulations. The flattening ratio, the local maximum substrate temperature, the amount of the particle embedded, and the reconstructed structure were examined in detail. There are three steps in the entire deposition process: absorption, collision and relaxation. The deposition method revealed at each point strikingly distinct morphological properties and processes of energy transformation.
The topography and splats measurement were studied explicitly using an automatic dislocation extraction (DXA) algorithm including penetration depth calculation, overall shape following sub-surface penetration and flattening measures at different temperatures and impact velocities. The topography of the splat depends heavily on impact velocity. With a 500 m/s velocity and at 500 K initial temperature, the topography is significantly hemispheric, pyramidal at 1000 m/s, and flat over the surface of the substrate like a mushroom at 1500 m/s. In addition, the splat diameter decreases with an increase in the velocity of S.T. Oyinbo et al.: Manufacturing Rev. 7, 29 (2020) impact, and so does the degree of splat dissemination in the sub-surface. Impact velocity has been discovered for optimal dense and uniform coating within a medium-range in this analysis as 500 to 1100 m/s. In order to ascertain a deposition of good quality during the cold gas dynamic spray process, it is important to maintain impact velocity within an optimum range. As the velocity of impact increases, the kinetic energy decreases, the temperature of the particle decreases, and the plastic strains increase considerably.
Analysis of the effect of particle size revealed that the height of the particle deposited increases as the particle size increased to 2 nm. Above 2 nm particle size, no significant differences are observed in the deposition indicating that optimal coating conditions used in this analysis are obtained with particle sizes varying in the range of 1 nm to 2 nm.
The analysis of von Mises shear strain, temperature distribution and plastic strain at higher initial velocities established shear instabilities of the system. The results of this study of molecular dynamic (MD) simulation comply with previously published experimental and numerical results. Our current knowledge of cold gas dynamic spray processes can also be enhanced with the results of this study.