Feasibility of numerical simulation methods on the Cold Gas Dynamic Spray (CGDS) Deposition process for ductile materials

The techniques of cold gas dynamic spray (CGDS) coating involve the deposition of solid, high speed micron to nano particles onto a substrate. In contrast to a thermal spray, CGDS does not melt particles to retain their physico-chemical properties. There have been many advantages in developing microscopic analysis of deformation mechanisms with numerical simulation methods. Therefore, this study focuses on four cardinal numerical methods of analysis which are: Lagrangian, Smoothed Particles Hydrodynamics (SPH), Arbitrary Lagrangian-Eulerian (ALE), and Coupled Eulerian-Lagrangian (CEL) to examine the Cold Gas Dynamic Spray (CGDS) deposition system by simulating and analyzing the contact/impact problem at deformation zone using ductile materials. The details of these four numerical approaches are explained with some aspects of analysis procedure, model description, material model, boundary conditions, contact algorithm and mesh refinement. It can be observed that the material of the particle greatly influences the deposition and the deformation than the material of the substrate. Concerning the particle, a higher-density material such as Cu has a higher initial kinetic energy, which leads to a larger contact area, a longer contact time and, therefore, better bonding between the particle and the substrate. All the numerical methods studied, however, can be used to analyze the contact/impact problem at deformation zone during cold gas dynamic spray process.


Introduction
The cold gas dynamic spray (CGDS) mechanism is based on the solid-state deposition technique. CGDS is suitable for various engineering applications, which include composites, ceramics, metals and polymers. The gradual shear and plastic deformation generated by impact velocity of the accelerated particle is accomplished by the expansion of pressurized gases through a nozzle thus, the metallurgical coalescence is produced between the particle and the substrate [1,2]. Plastic deformation of particles occurs as it impacts on the target surface to form a uniform layer. Only when the velocity of the sprayed material reaches a predefined velocity called critical velocity under given operating conditions can the particles/substrate bonding occur [1,3]. Nanoscale cold spraying is a potential technology for depositing or coating nanostructured materials on the surface of the substrate without affecting its properties or structure significantly [3][4][5]38]. The technology finds enormous applications including metal matrix composite (MMC), metallic, ceramic, and plastic coatings in vital engineering fields [3,[5][6][7]. The suitability of the material for CGDS depends on its physical and mechanical properties, including density, melting temperature, melting hardness and material hardness [4,5,8]. Relatively low yield material such as zinc, aluminium and copper are considered desirable because they display comparatively greater softening at high temperature [5,9,10], whereas high resistance materials are not perfect for CGDS, because of the lack of adequate energy available for deformation.
The technique of finite element analysis has been the focus of many researchers among several models of material impact phenomena developed to investigate the deformation mechanism due to its modelling capacity of material models, complex geometries and contact algorithms. In cold gas dynamic spraying, the particle/substrate impact can be termed a high-velocity impact process which can be handled by the Lagrangian method [11]. The Lagrangian numerical algorithm has been the focus of many researchers. 3D and axisymmetric models for particle and substrate impact was first used by Assadi et al. [12], to establish the impact phenomenon by using ABAQUS/Explicit version 6.2-1. Yin et al. [13], Li et al. [14], and Grujicic et al. [15] also investigated the behaviour of particles and substrate during impact with the application of Lagrangian analysis model.
Li et al. [14] were the first among the researchers that incorporated material damage mechanisms in the model as well as Lagrangian adaptive mesh domains (Lagrangian-Eulerian method) to control the excessive element distortion and mesh size control respectively. Xie et al. [16], provides the basic understanding of particle/substrate impact during cold spray by explicitly examined the different numerical model by modeling high-velocity impacts of spherical particles onto a flat substrate under various conditions. For the first time, they proposed the coupled Eulerian-Lagrangian (CEL) numerical approach as a means of solving the high-strain rate deformation problem.
Research to study the particles and substrate impact behaviour was carried out with Cu as the particle material using Eulerian Formulation, another ABAQUS/Explicit model [17]. At about 290 m/s minimum velocity, a jet was discovered to have formed and had maximum reached plastic strain (PEEQ). Therefore, at a velocity below 290 m/s, no jet could occur. This velocity was assumed to be a critical velocity. Jet formation discontinued and the material splashed at a speed greater than 290-400 m/s. A critical velocity could, therefore, be predicted in accordance with the theoretical analysis of jet morphology as a prediction tool by the Eulerian model [18][19][20].
Yildirim et al. [21] use several model reference domains in ABAQUS/Explicit package to systematically study the impact of a single particle on the semi-infinite substrate at an initial velocity of 100-700 m/s. The 2-dimensional axisymmetric structure was used in the case of Lagrangian and Arbitrary Lagrangian-Eulerian (ALE) model for the impact process, and one-quarter symmetry 3D model was used to study the material failure mechanism. The temperature was initiated at 293K (room temperature). The fictional formulation with surface-to-surface interaction was also incorporated in the 2D model and general contact algorithm for the 3D model with the friction coefficient of 0.3 at the particle and substrate interface. It was discovered that interpolation error occurs in the ALE adaptive remeshing technique with a significant decrease in equivalent plastic strain.
Another numerical approach in investigating the impact behaviour of cold spray particles is Smoothed Particle Hydrodynamics (SPH) [22]. Manap et al. [23] and Yildirim et al. [21] used the SPH approach for critical velocity prediction during the CS deposition process. Furthermore, the SPH method is appropriate for the multiple-particles impact process due to the appropriate solution techniques regarding the interface contact problem and its unique meshless feature [24]. However, the problem of tensile instability in SPH approach and the lack of interaction between the particles can lead to large tensile deformation which is a significant numerical problem.
With the introduction of a modern numerical system and complex representations of finite elements, the damage was induced to aeronautical structures by Smojver and Ivančevic [25] to predict the damages induced of a bird strike. Coupled Eulerian-Lagrangian (CEL), a modern finite element technique was however used to model and solved the soft body impacts. For modelled bird replacement material hydrodynamic reactions, material volumetric force and pressure-density ratio by the material equation of state (EOS) were used. The lagrangian bird model and experimental results are used to verify observations from the bird model with the CEL approach.
Gang et al. [26] explored the potential of the CEL numerical method to address geotechnical problems. Their studies have shown that CEL is capable of resolving difficult problems which FEM considers difficult to solve. In order to further explore the capacity of the CEL, a pile installation was simulated by a CEL approach and it was discovered that the CEL is better suitable for studying the pile influence on the soil-preceding structure's relationship, namely that the friction values are high when the simulation results and measuring data are considered. Therefore, they believe that due to the quality of parallelization, the CEL method achieves successful results.
This study presents four cardinal numerical methods of analysis which are: Lagrangian, Smoothed Particles Hydrodynamics (SPH), Arbitrary Lagrangian-Eulerian (ALE), and Coupled Eulerian-Lagrangian (CEL) to examine the Cold Gas Dynamic Spray (CGDS) deposition system by simulating and analyzing the contact/impact problem at deformation zone using ductile materials. This was done with the aim of accomplishing a qualitative understanding of when the particle deforms plastically during the cold gas dynamic spray process and find a way of addressing high-strain-rate dynamic problems of cold sprayed particles and substrate.

Problem description
This study investigates the feasibility of four numerical analysis models to simulate both the single and multiple particle impact on a deformable substrate during the CGDS process. Based on Abaqus Analysis User's Manual [27], an explicit-finite element analysis program was adopted for the analysis. A 3D model was established for deformable Copper (Cu) particle and deformable Aluminum (Al) substrate. The 500 m/s initial impact velocity used for this simulation of Cu/Al impact is below the critical velocity of copper/aluminum system. Note that 507 m/s is the approximated critical velocity for Cu/Al by using shear localization analysis [15]. Because of SEM observation, the morphology of copper particle was taken to be spherical for the numerical model ( Fig. 1) and its corresponding mean particle size was taken to be 10 mm in the CGDS process. The penalty formulation and general contact explicit available in ABAQUS/Explicit FEA program were used to describe the relative motion between the particle and substrate surface during cold gas dynamic spray process. The coefficient of friction was taken to be 0.3 for all the analyses [10,21]. The initial temperature is set to be 25°C in all the cases of Cu/Al impacts [28] with a fixed simulation time of 60 ns which is enough to study the contact process [10,16,29,30]. The cylindrical model was used for the substrate, the height and the radius were 8R and 16R respectively, where R is the diameter of the particle.
The Johnson-Cook plasticity model offers a definition of material movement for both the particles and substrate [29]. The flow stress (s) functions as illustrated in equation (1) are the strain hardening, strain rate hardening and temperature softening where the working hardening exponent (n) and plastic strain is denoted by _ e, the dimensionless plastic strain rate is the ratio of _ e=_ e 0 , _ e 0 = 1.0 s À1 and the substance constants A, B, C and m are shown in Table 1. T, T m and T 0 are the measured, melting and reference temperature respectively The thermal reaction study is carried out using the thermal conductivity properties and specific heat. Table 1 shows the properties of the materials used for the analysis [21,31,32] 3 Computational procedure 3

.1 Lagrangian model description
Abaqus/explicit program [27], was used to model the impact behaviour of Cu particle upon Al substrate in the Lagrangian simulation. The penalty formulation and general contact explicit available in ABAQUS/Explicit FEA program were used to describe the relative motion between the particle and substrate surface during cold gas dynamic spray process. The following analysis was also used for the Lagrangian domain; an 8-node thermally coupled brick (C3D8RT), reduced integration, hourglass control, trilinear displacement and temperature [27]. The mesh size of 0.0003 mm was used, i.e. the resolution of 1/100 particle diameter. Hexahedral meshing elements have been used. The application of boundary conditions with respect to the x-z plane involves symmetry, the bottom of the substrate with zero displacements, and the boundary conditions Lagrangian parts are shown in Figure 1c. Where particle initial velocity (v = 500 m/s), nodal rotation (r) and nodal displacement (u) are defined with the coordinates x, y z. The particle and substrate initial temperature of 25°C was used for all the calculations.
The conservative equations outlined in equations (4)-(7) represent the equation of mass, energy and momentum derived by the spatial time derivative approach [33]. Where r, s, u, E and e are the density, Cauchy stress, material velocity, total energy per unit volume and the internal energy respectively. The addition of the internal energy (e) and kinetic energy gives the total energy E (Eq. 7).

Smoothed particle hydrodynamics (SPH) model description
Another numerical model used in this study is Smoothed particle hydrodynamics (SPH). These methods belong to the meshless (or mesh-free)-family in which elements and nodes are not defined in the domain against the normal practice in the analysis of the finite element method; instead, the given body is represented by a collection of points. These nodes are generally called pseudo-particles or particles in smoothed particle hydrodynamics. Smoothed particle hydrodynamics is a modelling scheme in a fully Lagrangian domain in which continuum equations of a prescribed set are discretized by directly interpolating the properties over the solution region at a discrete set of points without necessarily define a spatial mesh. All the theory of energy conservation such as conservation of momentum, conservation of mass and conservation of energy are all defined with the Lagrangian domain. The interaction that occurs between the neighbouring pseudo-particles and current pseudo-particle over time causes the material to change and the field approximation at each step is done on the basis of locally distributed neighbouring particles. As the word 'particle' might suggest, the discrete particles (spheres) collide in compression with each other or in tension exhibiting cohesive-like behaviour. The SPH method at its core does not base on such phenomenon, rather, a method of discretization of continuum partial differential equations. In this analysis, the SPH reference frame was used only to model the particle because its deformation is much more than that of the substrate, whereas the substrate was modelled using the Lagrangian reference frame. The Explicit dynamic stress-displacement analysis was used in this method because the coupled mechanical Àthermal procedure in Lagrangian domain does not support the PC3D element type SPH approach. The element type for the substrate is an 8-node thermally coupled brick (C3D8RT), reduced integration, hourglass control, trilinear displacement and temperature, while the PC3D unique ABAQUS/Explicit element type was used for the particle in the SPH analysis. The interpolation theory is the foundation of the smoothed particle hydrodynamics (SPH) numerical approach. The interpolation theory transforms the continuum fluid dynamics conservation laws into an integral equation. Smoothing kernels (W) can be used at a certain position to obtain the kernel approximation of a function f(R 0 ) by integration over the computational domain [34].
where weighting function is W with respect to support scale (h). Equation (10) is the Delta-function property when the limit of smoothing length is tending to zero. Kernel functions have many possible choices. The third-order B-spline function was selected for this analysis to reduce the code frequency and the number of interactions of the particles. If the relative displacement (R) is defined between points Rnd R 0 , r = |R À R 0 |/h then equation (11) gives the B-spine function. The function of normalization N (d) is given to be {3/2, 7/10p, p, 31/5p 2 } ; d = 1, ..., 5. Herefore, the fluid flow momentum equation is given by equation (12): were the particle velocity, viscosity, density and pressure are respectively represented by v a , m a , r a , and p a for particle, a. For particle b, the velocity, viscosity, density, mass and pressure are represented by v b , m b , r b , m b and p b respectively. The position vector and interpolation kernel are denoted by R ab = R a À R b form particle a to b and W ab = W (R ab , h) respectively. Gravity vector and viscous term factor are denoted by g and j. At R ab = 0, smoothing of singularity is h. 3D SPH model was employed in this analysis. There is no special treatment for contact boundary condition in the SPH approach. The bottom surface of substrate is constrained by the PINNED boundary condition (the substrate is constraint in both x, y and z-direction). The condition of SPH particles geometric proximity is detected automatically to calculate the contact process of particle/substrate interface. The SPH interact with each other by meeting this condition following the requirement of boundary compatibility.

Arbitrary Lagrangian-Eulerian (ALE) model description
The Arbitrary Lagrangian-Eulerian (ALE) analysis is used to study large deformations of the transient problem by using Lagrangian adaptive mesh domains. In Figure 1a, the domain of Lagrangian adaptive mesh is created by the blue line region so that the orientation of the material present in this domain will follow the material flow path, which validates the most structural analyses of physical interpretation. The computational cost will be reduced by defining the domains of adaptive mesh as a fraction of the entire domain. The material direction is always proportional to the mesh and normal to the boundary on the Lagrangian domain boundary so that at all times the material domain is covered by the mesh. A smoother, new mesh is made in an adaptive meshing increase by iteratively sweeping over the domain of the adaptive mesh. Element distortion is reduced as each mesh sweep in the domain by relocating the nodes based on the current location of the neighbouring nodes as well as the elements. In each increment of adaptive meshing, the adaptive meshing intensity increases when the number of sweeps increases. The sole objective of the mesh smoothing method in the adaptive mesh domain is to improve element aspect ratios by minimizing mesh distortion expense of diffusing initial mesh gradation. Adaptive meshing robustness in ABAQUS/Explicit is achieved by adopting an enhanced algorithm based on the geometry of the evolving element. 2D axisymmetric ALE model was employed in this analysis. CAX4RT: a 4-node bilinear displacement, viscoelastic hourglass control, axisymmetric quadrilateral, reduced integration, thermally coupled and the temperature was used for the ALE domain. Here, 0.0001 mm was used as the mesh size. Quadrilateral elements were used for the meshing. The number of mesh sweeps and frequency for this analysis is 5 and 15 respectively.

Coupled Eulerian-Lagrangian model description
Another numerical model used in this study is the Coupled Eulerian-Lagrangian approach. This method is also used to study large deformations of transient problems. Coupled temperature-displacement elements for the Eulerian domain (EC3D8RT) and fully coupled thermal-stress analysis for the Lagrangian domain (C3D8RT) are used in this analysis. The volume-of-fluid method is the foundation of Coupled Eulerian-Lagrangian model for the implementation of the Eulerian part in ABAQUS/ Explicit. Within each element, this method computes the Eulerian Volume-Fraction (EVF) as the material flows and tracked in the mesh. Generally, one is the volume fraction of an element if it is completely filled with material and if the element contains no material, its volume fraction is zero. More than one material can be present in a Eulerian element simultaneously. If all the volume fractions of the material in an element are sum together and is less than one, automatically, the remaining elements will be occupied by 'void' material. There is no mass and weight for void material. The 3D model is implemented for this numerical approach. The geometry, boundary conditions, analysis procedure, mesh and interaction are the same as that of the pure Lagrangian method. In the case of Eulerian frame, prescribed velocity along x, y and z-direction are constraint.

Lagrangian numerical approach for single-particle impact model
At 50 m/s initial impact velocity, the equivalent plastic strain evolution (PEEQ) for Cu/Al impact after the simulation time of 30ns is shown in Figure 2. The peak value of PEEQ during the calculation is found to be at the interfacial zone between the particle and substrate.
As predicted by the Lagrangian method, the single spherical particle impacts the flat substrate, becomes flattens, and generate a crater at the edge of the contact region. At the outer interfacial region, there is an observation of an intensive plastic deformation between the Cu particle and Al substrate, where the value of PEEQ exceeds 2.0.
The temperature (TEMP) distribution at 500 m/s initial impact velocity for Cu/Al impact after the simulation time of 30 ns is shown in Figure 3. The process of deformation in cold gas dynamic spray is purely adiabatic [35]. Temperature distribution of material is dependent on the plastic deformation, and since the deformation of the particle is largely observed with a maximum plastic strain close to the interfacial region, the maximum TEMP is also found to be located near the interface. The particle and substrate materials formed jet at the interfacial region as the deformation of the particle and substrate proceeds. The appearance of the deformed particle is now lens-like shape. As the material of the particle and the substrate deform further, then, more jet builds up. Figures 2a-d and 3a-d shows the evolution of the jet impact time histories of 30 ns, and Figure 4a and b displays the TEMP and PEEQ distribution for the corresponding single-particle The material true time history impact at the impact time of 60 ns. The shear strain and temperature evolution is not enough to clearly predict adiabatic shear instability compared with [16,36] from these figures. As shown in Figure 4, the substrate equivalent plastic strain evolution is higher than that of the particle, thus higher temperature is obtained.

Smoothed particle hydrodynamics (SPH)
numerical approach for single-particle impact model SPH deformation pattern for Cu/Al impact is different from Lagrangian and ALE deformation pattern of particle/ substrate impact. At the edge of the interfacial zone, there is no material jet of a particle, and the particle penetration in the substrate is deeper.
At 500 m/s initial impact velocity, the evolution of equivalent plastic strain (PEEQ) for Cu/Al impact after the simulation time of 30 ns is shown in Figure 5. In these figures, the particle aspect ratio is decreased, with an increase in the collision time of the contact process, and the width and depth of the substrate increase. The plastic strain observed in the substrate increases greatly after 5 ns, whereas, the maximum equivalent plastic strain continues to approach a horizontal asymptote. A viscous-like resistance created by the excessive deformation can hinder the further deformation process under high strain rate and high impact velocity. This effect can be attributed to the second term Johnson-cook of the stress-strain law which describes the strain rate hardening. Figure 6b shows the stress evolution at 25 ns impact time which is however different from the PEEQ curve in Figure 6a. When the substrate temperature increases, there is low resistance of the material to shear flow when thermal softening is considered. This means that if any amount of shear stress is applied as the material approaching melting temperature, the shear strength of the material will be lost, and it will experience excessive deformation. However, S.T. Oyinbo and T.-C. Jen: Manufacturing Rev. 7, 24 (2020) this analysis is isotherm, and the temperature is kept constant of the whole system at 25°C. Moreover, material damage model is not employed in the material. The Dynamic-Explicit procedure used in this analysis produces much higher PEEQ as clearly observed in Figure 6. The conduction of heat is obvious in the case of Dynamic-Temperature Displacement-Explicit procedure from the impact zone and conducted into the inner part of the particle and substrate. Although for large models' analysis with extremely discontinuous events and relatively short times response, Dynamic-Explicit procedure is recommended because of its computational efficiency.

Arbitrary Lagrangian-Eulerian (ALE) numerical approach for single-particle impact model
The ALE numerical approach is also used to solve the single-impact problem. The equivalent plastic strain (PEEQ) evolution and temperature distribution at 500 m/s initial impact velocity for Cu/Al impact during the simulation time of 30 ns are shown in Figures 7 and 8 respectively. Excessive distortion of mesh does not occur by using this numerical method. The material jet formed in this analysis is smoother at the interfacial region instead of acute and thin. The value of frequency applied in the analysis is 15 and remeshing sweeps per increment are 5. The absence of adhesion model makes the particle rebound to occur after 30 ns of the impact time. It is observed that the temperature at the outer region of the particle and substrate impact is higher than what the inner part experiences because of the locally formed plastic deformation around the region surrounding the interface.
The evolution of the particle and substrate equivalent plastic strain (PEEQ) are shown in Figure 9b at 60 ns impact time. The plastic strain observed in the substrate increases unreasonably after 10 ns, whereas, the maximum equivalent plastic strain continues to approach a horizontal asymptote. However, the monotonical increase in the history of PEEQ until it reaches the plateau is another feature in a pure Lagrangian approach. The sudden decrease of the PEEQ after reaching the peak value is an unrealistic behaviour and can be attributed to errors of interpolation caused in the remapping point and surrounding interface of high strain gradients by the adaptive meshing algorithm [16,37]. The deviation of material points and integration points is due to the usage of adaptive meshing with some features similar to CEL approach. Therefore, the adaptive meshing and material motion composite effect is represented by the motion of the interior mesh of an adaptive mesh domain. Figure 9a shows the temperature (TEMP) distribution at 60 ns after impact. Before 10 ns, the TEMP slowly increases, but because of the high strain rate, it increases sharply between 10 and 20 ns, thereafter, there is a gradual decrease as the calculation continues. Because the substrate experiences larger plastic deformation, its temperature is higher than that of the particle. Adaptive meshing frequency and intensity are the key factors that affect the computational cost and simulation results in this analysis.

Coupled Eulerian-Lagrangian (CEL) numerical approach for single-particle impact model
The problem of mesh excessive distortion and abnormal deformation can be avoided using the Coupled Eulerian-Lagrangian numerical approach. CEL particle model enhances the modelling of the fluid-like particle as shown in Figure 10. Although there is still an occurrence of the material jet in this analysis and it has no effect on the completion of the analysis. The interpretation of the results from CEL should be different from that of a Lagrangian method (Figs. 10 and 11). Since the Eulerian part in CEL analysis is rigid and fixed, any nodal displacements result is      meaningless. Therefore, the calculation of PEEQ of all materials within the element is based on volume fraction weighted average (PEEQVAVG). The value of this volume average is significantly lower whether the adaptive remeshing is used or not than that of the Lagrangian analysis. The same explanation can as well be applied to the particle temperature (TEMPMAVG) distribution as shown in Figure 11. The equivalent plastic strain (PEEQ) evolution and temperature distribution at 500 m/s initial impact velocity for Cu/Al impact during the simulation time of 20 ns for particle and substrate are shown in Figure 12. The PEEQ evolution in the substrate increases rapidly until it reaches the peak value of 5. The PEEQAVG for the substrate is significantly higher than that of the particle. One of the CEL analysis shortcomings is the inability to trace the materials history behaviour.

Comparison of four numerical methods
At 500 m/s initial impact velocity, the evolution of equivalent plastic strain (PEEQ) for Cu/Al impact after the simulation time of 60 ns is shown in Figure 13 calculated by Lagrangian, SPH, ALE, and CEL numerical model. As predicted by the Lagrangian method (Fig. 13a), the single spherical particle impacts the flat substrate and becomes flattens and generate a crater at the edge of the contact region. The particle and substrate materials formed jet at the interfacial region as the deformation of the particle and substrate proceeds. The appearance of the deformed particle is now lens-like shape. As the material of the particle and the substrate deform further, then, more jet builds up. For the SPH numerical method, as indicated in Figure 13b, there is no material jet of the particle at the interfacial region edge. The PEEQ distribution in this analysis, when compared to the Lagrangian result, was similar. Excessive distortion of mesh does not occur by using ALE numerical method. The material jet formed in this analysis is smoother at the interfacial region instead of acute and thin (Fig. 13c). CEL particle model enhances the modelling of the fluid-like particle (Fig. 13d). Although the occurrence of the material jet in this analysis has no effect on the completion of the analysis, and the particle penetration in the substrate is deeper than other numerical methods. Figure 14a shows the normalised kinetic energy over the period of 60 ns impact time using Cu/Al impact for the four numerical approaches under consideration. In the deposition of Cold Spray process, energy from initial kinetic energy (ALLKE) is converted into energy saved in    the particle/substrate (ALLSE), the energy that deforms the material plastically (ALLPD) and energy that propagates the stress wave. After plastic deformation of the material, the process is irreversible, and the kinetic energy stored, is the energy which can be recovered during the restitution. This energy is known as rebound kinetic energy. The initial kinetic energy is much higher than the rebound kinetic energy. That means over 98% of the kinetic energy is converted into internal energy and 2% is converted into rebound kinetic energy. The kinetic energy which is then produced by rebounded particles has not been monotonous, because, after 20 ns, the kinetic energy damped periodically. The recoverable elastic energy matches this energy. The pattern of the normalized kinetic energy in all cases of numerical approaches agrees well except for the CEL approach. The deviation is due to the different approach used in the analysis.
The removal of hourglass control (singular modes) is usually done by the associated 'artificial' strain energy. The strain energy that will control hourglass deformation will be too much if there is excessive artificial strain energy during the process. Figure 14b shows how to determine whether the artificial strain energy is excessive or not by comparing the internal energy (ALLIE) with that of artificial strain energy (ALLAE). Generally, the proportion of when ALLAE divides ALLIE should not be up to 5%.
The initial ratios close to 0.0 s can be neglected because it is basically noise produced when a very small number is divided by another, all four numerical approaches have a ratio of less than 4%. The problem of the hourglass can be sufficiently prevented by intrinsic hourglass control. The fine mesh size used in this analysis is enough to stop zeroenergy modes propagation which potentially yields inaccurate results. Table 2 presents a comparison between the four numerical approaches in term of computational costs. The most efficient approach among the four is the pure Lagrangian approach, and the SPH consumes more time.

Conclusion
The Finite Element analysis has been carried out using four different numerical approaches -Lagrangian, Smoothed Particles Hydrodynamics (SPH), Arbitrary Lagrangian-Eulerian (ALE), and Coupled Eulerian-Lagrangian (CEL), to examine the Cold Gas Dynamic Spray (CGDS) deposition system, through simulating and analyzing the contacts/impacts at the deformation zones. It can be observed that the particle material has greater influences on the deposition process and the deformations than the substrate material does. Regarding the particle, a material  with higher densitysuch as Cu has a higher initial kinetic energy, leading to a larger deformation area and a longer contacttime, and hence, suggesting a better bonding between the particle and the substrate. The study suggested that all the numerical methods tested could be used to analyze the contact/impact problems at the deformation zones in cold gas dynamic spray processes. The higher computational efficiency of the Lagrangian approach and its ability to incorporate a complex material model into the simulation, nevertheless, makes it one of the most suitable numerical methods. However, the severe distortion of the mesh structure in the deformed area could result in non-convergence of simulation and inaccuracy of the calculated result, due to the numerical backslide effects.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.