Numerical Study of Powder Flow Nozzle for Laser-Assisted Metal Deposition

: Metal additive manufacturing has received much attention in the past few decades, and it offers a variety of technologies for three-dimensional object production. One of such technologies, allowing large-sized object production, is laser-assisted metal deposition, the limits of which are determined by the capabilities of the positioning system. The already-existing nozzles have either a relatively low build rate or a poor resolution. The goal of this work is to develop a new nozzle with a centered particle beam at high velocity for the laser-assisted metal additive manufacturing technologies. Scientiﬁc challenges are addressed with regards to the ﬂuid dynamics, the particle-substrate contact, and tracking of the thermodynamic state during contact. In this paper, two nozzles based on the de Laval geometry with Witoszynski and Bicubic curves of convergence zone were designed; the results showed that the average ﬂow velocity in a Bicubic outlet curve nozzle is around 615 m/s and in Witoszynski this is 435 m/s. Investigation of particle beam formation for the Bicubic curve geometry revealed that small particles have the highest velocity and the lowest total force at the nozzle outlet. Fine particles have a shorter response time, and therefore, a smaller dispersion area. The elasto-plastic particle-surface contact showed that particles of diameter limited up to 3 µ m are able to reach experimentally obtained critical velocity without additional heating. For particle sizes above 10 µ m, additional heating is needed for deposition. The maximum coefﬁcient of restitution (COR) is achieved with a particle size of 30 µ m; smaller particles are characterized by the values of COR, which are lower due to a relatively high velocity. Particles larger than 30 µ m are scalable, characterized by a small change in velocity and a rise in temperature as their mass increases.


Introduction
In recent decades, the additive manufacturing (AM) technology, used to fabricate physical models, prototypes, or functional gradient components directly from 3D computeraided design (CAD) data using the layered approach [1,2], has been evolving rapidly.In particular, the laser additive manufacturing (LAM) [3] demonstrates the advantage of printing large metal components and having a complex internal structure, as well as opens many new applications.The most popular LAM methods include selective laser sintering/melting (SLS [4]/SLM [5]) and directed energy deposition (DED) [6].Compared to SLS/SLM methods, DED represents the potential of different processes, especially in the repair of worn components [7] and in the preparation of functional gradient materials [8].
One of the AM classes belonging to the DED class is the Laser Powder Fed Additive Manufacturing Process (LPF-AM).LPF-AM functionality has attracted several industries to customize its coating, repair, and component manufacturing features from one or more materials.In LPF-AM, the powder nozzle feeds the powder into a melt pool supported by the laser beam.As the laser beam moves forward along with the powder jet, a weld path is formed, which forms a coating or 3D portion along with several other overlapping weld paths.This process is also known as laser metal deposition (LMD).
LMD technology allows printing 3D objects in large sizes, the limits of which are determined by the capabilities of the positioning system.Using LMD technologies, the printing of 3D objects can be combined with metal processing processes and products with places that are not accessible for processing.Due to their complex shape, surface treatment can be performed at the selected printing stage.Another advantage of LMD technology is the recovery of damaged metal objects and the gradual coating of the surface with other metals [9].
LMD systems use three types of printheads, namely off-axial nozzle, discrete coaxial nozzle, and continuous coaxial nozzle.In the last two, the flow of powder fed from the sides interacts with the laser beam and attenuates the intensity of the beam, resulting in more laser power being used than required by the metal melting process.Laser-heated powder particles in such heads fall into the melt pool at different angles, which ensures freedom of the printing directions of the head but reduces the deposition rate of the metal due to the peculiarities of the metal powder flow.The first type of LMD print head off-axial nozzle has a high metal deposition speed, which leads not only to time savings, but also to a better adaptability in the repair of damaged metal parts.The main disadvantages of such heads are single-direction print and lower resolution than those of off-axial nozzle heads.If the first problem can be solved by simply tilting the head with a print direction, then nozzle optimization is required to increase the resolution [10,11].
In the off-axial nozzle heads used, the flow of metal powder expands as it exits the nozzle and is in the shape of an expanding cone [12].Moreover, in such nozzles, the metal powder is evenly mixed with the carrier gas and as an abrasive destroys the nozzle output, which changes the printing parameters over time [13].
The additional gas flow around the metal powder flow helps to regulate and maintain the powder flow evenly along the axis of fall [14].Such a solution is suitable for low gas flow rates, and de Laval nozzles are used to increase the gas and powder flow rate.Usually, de Laval nozzles are used in Cold spray systems, where the powder flow rate is a critical parameter for effective surface covering [15].
The particles are usually fed into a laser-heated spot, where a melt is formed, which suddenly hardens as the laser beam travels onward [16].One of the main disadvantages is the low resolution of the melted material and the formed track, which results in a high surface roughness.This can be reduced by using a smaller melt pool and a thinner layer thickness, but this can have a negative effect on the speed of the process.Another negative aspect is the relatively low efficiency of the trapped powder, which in some cases can be as low as 8% [17].Recently, several concepts of metal particle positioning heads have been introduced to increase particle deposition efficiency, reduce particle scattering, etc. [18].These are usually coaxial flow heads in which the laser radiation passes through the central axis of the head and the particles are fed from the sides [19,20].Commonly, they consist of two matched cones that form an annular powder outlet.However, such already-existing nozzles have either a relatively low build rate or a poor resolution.Better resolution cannot be achieved because a better focused particle beam cannot be formed due to the limitations of the geometry (the laser beam is positioned in the center of the nozzle).
Heads with a coaxial particle flow and a non-centered laser beam exist as well [21].The main differences are how the particles are fed-whether separate sprays are used [22], whether it is coaxial injection, or whether there are additional gas streams that better focus the particle flow and laser radiation.However, there is no nozzle with a centered particle beam at high velocity.Our approach is different from the already-existing nozzles, namely, the high velocity particle flow is formed in the center of the nozzle with a few laser beams positioned around it.In this way, better resolution can be achieved.
The goal of this work is to develop a new nozzle with a centered particle beam at high velocity for the laser-assisted metal additive manufacturing technologies.For that, we attain the following aims:

•
To design and to investigate two nozzles (based on the de Laval geometry with Witoszynski and Bicubic curves of convergence zone) having a centered capillary for the particle beam and sheath compressed gas flow around the capillary.

•
To determine optimal Witoszynski and Bicubic curve parameters for the nozzle, which are characterized by high gas velocity (>500 m/s), lower pressure at the end of the capillary, and the smallest gas flow velocity fluctuations in the first 20 mm after exiting the outlet of the nozzle.

•
To investigate particle beam formation for optimal printing resolution.

•
To determine relationship between particle sizes and coefficients of restitution, by performing particle-substrate thermo-mechanical elastic-plastic contact analysis.
In this work, scientific challenges are addressed with regards to the fluid dynamics (chapters 3.1 and 3.3), the particle-substrate contact (chapter 3.4), and tracking of the thermodynamic state during contact (chapter 3.5).Such an approach leads to the optimal design of a nozzle.Our technology works a little bit differently than most directed energy deposition (DED) technologies and has similarities with the cold spray technology.In our case, we are developing a technology where the powder particles attach to a solid surface (substrate or an already-existing layer) while still being in solid state.Only after their sticking to the existing layer, the newly deposited particles are melted by a laser beam.For this reason, the calculations of the described impact effect of a solid particle against a solid surface are significant in our study.We expect our technology to be more precise than the already-existing DED, LMD, or laser cladding technologies.

Methods
De Laval nozzles consist of four parts, namely the straight part, the convergent part, the throat, and the spread part (Figure 1).The flow of compressed gas (Table 1) has a high degree of turbulence, which is not conducive to gas shrinkage in the convergent section of the nozzle.Therefore, a straight section is established in front of the convergent part of the nozzle to ensure an even distribution of the flow inlet velocity.
Mathematics 2021, 9, x FOR PEER REVIEW 3 of 17 the particle flow and laser radiation.However, there is no nozzle with a centered particle beam at high velocity.Our approach is different from the already-existing nozzles, namely, the high velocity particle flow is formed in the center of the nozzle with a few laser beams positioned around it.In this way, better resolution can be achieved.The goal of this work is to develop a new nozzle with a centered particle beam at high velocity for the laser-assisted metal additive manufacturing technologies.For that, we attain the following aims:

•
To design and to investigate two nozzles (based on the de Laval geometry with Witoszynski and Bicubic curves of convergence zone) having a centered capillary for the particle beam and sheath compressed gas flow around the capillary.

•
To determine optimal Witoszynski and Bicubic curve parameters for the nozzle, which are characterized by high gas velocity (>500 m/s), lower pressure at the end of the capillary, and the smallest gas flow velocity fluctuations in the first 20 mm after exiting the outlet of the nozzle.

•
To investigate particle beam formation for optimal printing resolution.

•
To determine relationship between particle sizes and coefficients of restitution, by performing particle-substrate thermo-mechanical elastic-plastic contact analysis.
In this work, scientific challenges are addressed with regards to the fluid dynamics (chapters 3.1 and 3.3), the particle-substrate contact (chapter 3.4), and tracking of the thermodynamic state during contact (chapter 3.5).Such an approach leads to the optimal design of a nozzle.Our technology works a little bit differently than most directed energy deposition (DED) technologies and has similarities with the cold spray technology.In our case, we are developing a technology where the powder particles attach to a solid surface (substrate or an already-existing layer) while still being in solid state.Only after their sticking to the existing layer, the newly deposited particles are melted by a laser beam.For this reason, the calculations of the described impact effect of a solid particle against a solid surface are significant in our study.We expect our technology to be more precise than the already-existing DED, LMD, or laser cladding technologies.

Methods
De Laval nozzles consist of four parts, namely the straight part, the convergent part, the throat, and the spread part (Figure 1).The flow of compressed gas (Table 1) has a high degree of turbulence, which is not conducive to gas shrinkage in the convergent section of the nozzle.Therefore, a straight section is established in front of the convergent part of the nozzle to ensure an even distribution of the flow inlet velocity.

Convergent Section
The shrinkage part is an important part of the supersonic nozzle, which evenly accelerates air flow, reduces turbulence, and improves the stability of gas flow.There are currently many types of shrinkage design curves, the most common types of curves being a double parametric curve and its modifications, a Witoszynski curve and its modifications, and a straight line.The best parameters of speed and stability have been shown by the double parametric curve and the Witoszynski curve [23,24].

Bicubic Parametric Curve
The Bicubic parametric curve is a family of curves with the parameter x m .By selecting the x m position, a number of contraction surfaces can be obtained.The tapering profile R is calculated by the Equation (1): where x m is the coordinate of the point of intersection of the two curves (defining the end of the first curve and the beginning of the second curve).The value of x m in this paper is in the range of 0.2-0.9(Table 2).The choice is explained in Section 3.1.

Witoszynski Parametric Curve
The Witoszynski curve is used for an ideal, uncompressed, and axially symmetric flow that can be used to obtain a well-distributed and stable gas flow.The Witoszynski curve is obtained according to the described Equation ( 2): where R in is the radius of the inlet, R c is the radius of the nozzle throat, and R x denotes the radius of the nozzle cross-section in any x along the convergence section.L 1 is the length of the convergence section and coefficients k 1 and k 2 are used to describe the sharpness of the narrowing of the Witoszynski curve.In this paper, the values of k 1 and k 2 are in the range of 10-900 (Table 2).The choice is explained in Section 3.1.

Divergent Section
To select the geometry of the divergence section, a Quadratic Bezier curve was used to describe it [25].The Quadratic Bezier curve can uniformly describe the most commonly used geometry of a divergence section and their variations in a convex parabola and cone.
Quadratic Bezier curve equation is defined by a set of control points P 0 , P 1 , and P 2 through P n , where n refers to the order of the curve (n = 1 for linear, 2 for quadratic, etc.).The first and last control points are always the endpoints of the curve; however, the intermediate control points generally do not lie on the curve.The equation of the quadratic Bezier curve is Equation (3): where: P 0 = (x 0 , y 0 ), P 1 = (x 1 , y 1 ), P 2 = (x 2 , y 2 ), t = (∈ 0, 1) Given the coordinates of control points P n , the first control point has coordinates P 1 = (x 1 , y 1 ), the second P 2 = (x 2 , y 2 ), etc., the curve coordinates are described by the equation that depends on the parameter t from the segment P 0 = (0, R c ), P 1 = (4, 2), and P 3 = (L 2 , R out ).

Parametrization of the Nozzle Geometry
The de Laval nozzle is used to accelerate the powder coming out of the capillary, so the central part of the nozzle has a holder fixating the capillary tube in the central axis.The geometry calculations were carried out using COMSOL Multiphysics 5.6 software applying the formulas listed above.
The High Mach Number flow, Spalart-Allmaras (RANS model) interface is used for the proposed geometries calculation.Fully compressible Navier-Stokes equations are used for a simple compressible fluid.The physics interface assumes that the fluid is an ideal gas.This is necessary for the formulation of the Consistent Inlet and Outlet Conditions.The total pressure of the internal inlet is 1 bar, while the total pressure of the external inlet is 13 bar.Nitrogen dynamic viscosity, specific gas constant, and heat capacity at constant pressure are given as a piecewise function dependent on the temperature.Density as an analytic function is temperature-dependent.The viscosity and thermal conductivity of an ideal gas can be accurately approximated using Sutherland's Law.Sutherland's formula is an approximation for how the viscosity of gases depends on the temperature.This law is based on an idealized intermolecular-force potential.Nitrogen's Sutherland constants are given in Table 1.The nozzle parameters used for the calculations and their limits are presented in Table 2.The parameter ranges listed in Table 2 were determined by the geometry limits.Parameter values out of the chosen ranges would cause geometry overlaps.The values for L 1 , L 2 , and R in were chosen as they ensure the production of the most compact nozzle.
Variation steps for k 1 and k 2 were 10, for P 1x -1 mm, for R c , P 1y , and R out -0.1 mm for each test of the Witoszynski curve geometry.The total number of tests was 124.For each test of the Bicubic curve geometry, the variation steps were as follows: X m -0.1, P 1x -1 mm, R c , P 1y , and R out -0.1 mm.The total number of tests was 124.
By changing the sets of parameters listed in Table 2, different gas flow parameters (velocity and pressure along the central axis of the nozzle) were calculated, as illustrated graphically in Figure 2. Figure 2 qualitatively illustrates the influence of different parameters of Witoszynski and Bicubic nozzle geometry on gas flow pressure and velocity.It can be observed from the graphs that Witoszynski and Bicubic geometries react differently to the change of parameters.Each variable is presented in Figure 2 by a family of curves.
Comparing the gas flow pressure shift of Witoszynski (Figure 2A) and Bicubic (Figure 2C) curve nozzle geometry, we see that the Witoszynski geometry responds to changes in parameters more than Bicubic.The Bicubic curve geometry proves to be very sensitive to the radius of the nozzle throat parameter, while the Witoszynski curve geometry shows a larger gas flow pressure drop in the capillary tube, a smoother transition between parameters, and lower fluctuations in the nozzle outlet.By comparing the velocity parameter shift of Witoszynski (Figure 2B) and Bicubic (Figure 2D) curve nozzle geometry, we see that Witoszynski geometry shows higher velocities in the capillary tube, but lower velocities in the divergent section and a sharper decay than the Bicubic curve geometry.By changing the sets of parameters listed in Table 2, different gas flow parameters (velocity and pressure along the central axis of the nozzle) were calculated, as illustrated graphically in Figure 2.  Comparing the gas flow pressure shift of Witoszynski (Figure 2A) and Bicubic (Figure 2C) curve nozzle geometry, we see that the Witoszynski geometry responds to changes in parameters more than Bicubic.The Bicubic curve geometry proves to be very sensitive to the radius of the nozzle throat parameter, while the Witoszynski curve geometry shows a larger gas flow pressure drop in the capillary tube, a smoother transition between parameters, and lower fluctuations in the nozzle outlet.By comparing the velocity parameter shift of Witoszynski (Figure 2B) and Bicubic (Figure 2D) curve nozzle geometry, we see that Witoszynski geometry shows higher velocities in the capillary tube, While trying to determine the optimal parameters for the convergent part of the Witoszynski curve nozzle, it was established that the coefficient k 1 has a bigger impact on the gas flow velocity compared with k 2 .The optimal value k 1 = 100 was chosen in consideration of sufficient gas flow velocity in the powder delivery capillary and little flow velocity fluctuations after exiting the nozzle outlet.When the said value is below 100, the flow velocity starts increasing inside the nozzle and decreases drastically upon exiting the outlet.When k 1 > 100, the gas flow in the capillary decreases rapidly resulting in less efficient powder delivery.Similar gas flow dependence on the value of the coefficient k 2 can be observed, and the optimal value for k 2 is the same as k 1 (k 2 = 100).
While trying to determine the optimal parameters for the convergent part of the Bicubic curve nozzle, it was established that the coefficient X m and the coefficients k 1 and k 2 in the Witoszynski case have a similar effect on the gas flow velocity.The optimal value for the coefficient X m is 0.7.When X m < 0.7, gas flow increases inside the capillary but decreases in the divergent zone and at the outlet.Similar results are observed when X m > 0.7.
The radius of the thinnest part of the de Laval geometry (R c ) is 1 mm.The same value was chosen for both Witoszynski and Bicubic curve nozzles.When R c < 1, gas flow increases inside the capillary but decreases at the outlet.When R c > 1, gas flow velocity maintains stability for a longer time upon exiting the outlet.However, flow velocity in the capillary decreases rapidly and negatively affects the efficiency of powder delivery.
The parameters for the divergence part of the nozzle were chosen the same for Witoszynski and Bicubic curve nozzles, except for the parameter P 1y .It was determined from CFD calculations that the gas flow is less turbulent when the P 1y value is 1.7 mm for Witoszynski and 2.0 mm for Bicubic.
It is suspected that the difference is due to an uneven narrowing of the convergent part in the Witoszynski and Bicubic curves; in the case of Witoszynski it has a constant tapered form, in the case of the Bicubic it has a sudden narrowing.The radius of the nozzle outlet (R out ) determines gas flow velocity of the divergent part of the nozzle.The highest flow velocity is observed when R out = 2 mm.When R out < 2 mm, velocity increases inside the nozzle but sharply decreases at the outlet.When R out > 2 mm, the drop in gas flow velocity inside the capillary is significantly larger compared to the optimal R out value.

Determination of the Nozzle Geometry
The criteria for the optimum of a nozzle are the maximum gas flow velocity at the critical distance (from the nozzle outlet to the substrate) and the minimum pressure at the end of the capillary.One Witoszynski and one Bicubic curve corresponding to these criteria were used to choose an optimal geometry for the nozzle development.Curve parameter values obtained from numerical calculations are listed in Table 3.To choose which nozzle geometry to use, a comparison of velocity and pressure has been carried out.Calculation results of the central nozzle axis are provided in Figure 3.The flow is accelerated in the convergent section of the nozzle and a velocity of 180-200 m/s is reached in Witoszynski and Bicubic curve nozzles (Figure 3A).Upon entering the divergent section, the flow is accelerated to 640 m/s in the Witoszynski and 700 m/s in the Bicubic curve nozzle.Velocities in the same range (600-1000 m/s) are used in extremely high speed laser metal deposition technologies [9].When the flow reaches the nozzle outlet, a decrease of its velocity can be observed from the graph.The critical distance from the nozzle outlet to the substrate is 20 mm and it is crucial to maintain an accelerated and stable particle flow until the particles reach the surface.To choose which nozzle geometry to use, a comparison of velocity and pressure has been carried out.Calculation results of the central nozzle axis are provided in Figure 3.The flow is accelerated in the convergent section of the nozzle and a velocity of 180-200 m/s is reached in Witoszynski and Bicubic curve nozzles (Figure 3A).Upon entering the divergent section, the flow is accelerated to 640 m/s in the Witoszynski and 700 m/s in the Bicubic curve nozzle.Velocities in the same range (600-1000 m/s) are used in extremely high speed laser metal deposition technologies [9].When the flow reaches the nozzle outlet, a decrease of its velocity can be observed from the graph.The critical distance from the nozzle outlet to the substrate is 20 mm and it is crucial to maintain an accelerated and stable particle flow until the particles reach the surface.The average flow velocity in the critical distance after exiting the outlet of the Bicubic curve nozzle is around 615 m/s.For the Witoszynski curve nozzle, the velocity is around 435 m/s.Fluctuations in velocity can be observed in both nozzles upon exiting the outlet.Despite being more apparent in the Bicubic curve nozzle, the average flow velocity is maintained higher in the distance from the outlet to the critical 20 mm point compared to the Witoszynski curve nozzle.
To ensure high flow velocity and avoid particle clogging, it is crucial to create subatmospheric pressure.Said subatmospheric pressure is achieved at the end of the capillary, as seen in Figure 3B.After this point, a rapid increase of pressure helps accelerate the particle flow in the divergent section of the nozzle.It can be seen from the graphs in Figure 3B that the pressure drop for the Witoszynski curve geometry is larger compared with the Bicubic.However, the pressure drop difference is not as significant as the difference in gas flow velocity.Thus, the Bicubic curve geometry is more suitable for the laser-assisted metal deposition application.

Particle Beam Formation
Particle tracing for fluid flow interface is used to simulate the motion of particles in a background fluid.We are specifying the particle density and diameter in our numerical calculation.It is known that the velocity decreases with the increase of particle size.Three different sizes of copper particles are used in our simulation: 1 μm, 10 μm, and 40 μm (Figure 4).The 1 μm particle was chosen as a particle with maximum velocity, 10 μm with The average flow velocity in the critical distance after exiting the outlet of the Bicubic curve nozzle is around 615 m/s.For the Witoszynski curve nozzle, the velocity is around 435 m/s.Fluctuations in velocity can be observed in both nozzles upon exiting the outlet.Despite being more apparent in the Bicubic curve nozzle, the average flow velocity is maintained higher in the distance from the outlet to the critical 20 mm point compared to the Witoszynski curve nozzle.
To ensure high flow velocity and avoid particle clogging, it is crucial to create subatmospheric pressure.Said subatmospheric pressure is achieved at the end of the capillary, as seen in Figure 3B.After this point, a rapid increase of pressure helps accelerate the particle flow in the divergent section of the nozzle.It can be seen from the graphs in Figure 3B that the pressure drop for the Witoszynski curve geometry is larger compared with the Bicubic.However, the pressure drop difference is not as significant as the difference in gas flow velocity.Thus, the Bicubic curve geometry is more suitable for the laser-assisted metal deposition application.

Particle Beam Formation
Particle tracing for fluid flow interface is used to simulate the motion of particles in a background fluid.We are specifying the particle density and diameter in our numerical calculation.It is known that the velocity decreases with the increase of particle size.Three different sizes of copper particles are used in our simulation: 1 µm, 10 µm, and 40 µm (Figure 4).The 1 µm particle was chosen as a particle with maximum velocity, 10 µm with half velocity, and 40 µm with low velocity.Particles larger than 40 µm have similar low velocity as well; therefore, 40 µm was chosen as the maximum particle size.
beam have a dispersion from 3.5 to 7.0 mm for 20 μm particles [14].Here, the researchers developed a nozzle with a centered particle beam; however, the particles enter the nozzle in the convergent zone.As a result, particles reach a higher velocity (471 m/s for 20 μm particles) and can cause erosion at the nozzle throat.Although the nozzle proposed in our paper has a lower particle beam velocity (360 m/s for 10 μm particles), it has a smaller particle trajectory dispersion.Additionally, due to the particles entering the nozzle in the divergent zone, the throat erosion is eliminated.The developed fluid model was applied for the evaluation of particle impact velocities.Particle tracing trajectory model was coupled with a high Mach number flow interface.The results of CFD simulation allowed to generate required velocities for differently sized powder particles.Particle release and propagation formulation is Newtonian.When particles contact with a wall, we can determine the bounce wall condition, where kinetic energy of the particles is conserved and the particles specularly reflect from the wall.Two forces are used on particles in a fluid: drag force and lift force.The drag force is defined as the Stokes drag law: where τ p is the particle velocity response time, v is the velocity of the particle, and u is the fluid velocity as the particle's position.The validity of a given drag law depends on the relative Reynolds number and the particle velocity response time: Re r = ρ|u − v|d p µ (5) The standard drag correlation is used in our simulation with defined drag coefficient as a piecewise function of the relative Reynolds number, which can reach up to 1.6•10 4 in our numerical calculation.The drag coefficient for Reynolds number in the range of 1.2•10 4 < Re r ≤ 4.4•10 4 is defined as: where d p is the particle diameter, ρ p is the density of the fluid, and µ is the dynamic viscosity of the fluid.Fluid flow in our system is nonuniform because we assume to add lift force on particles in a fluid.The wall-induced lift force in Newtonian formulations is applied to account for the effects of the nearby walls as particles move through the channel.The wall-induced lift force F L is given by: where: In Equation (8) r p is the particle radius, u is the fluid velocity, I is the identity matrix, n is the wall normal at the nearest point on the reference wall, D is the distance between the channel walls, s is the distance from the particle to the reference wall, and G 1 and G 2 are functions of the nondimensionalized wall distance s.The number of particles used in the simulation is 3 × 10 3 (three particles are released into the modeling domain every 10 ns).The difference between the highest and lowest values of the particle's size in our investigation is greater than 10 µm (from 1 to 40 µm).When statistical analysis is applied to particle size distributions, it is routinely based on lognormal distributions, where the data range of the particle's size is greater than about 10 µm.There is no real theoretical reason as to why lognormal distribution is used in this case, it is merely empirically the best fit.The number of particles released at the same time is based on the geometric standard deviation of the lognormal distribution.The diameter of the circle channel in our investigation is 410 µm.In the ideal case, nine particles with the diameter of 40 µm can be released at the same time.According to the lognormal distribution, a more statistical case would be six particles at the same time.CFD and particle tracing module used the plane symmetry geometry, as shown in Figure 4.It means that half of the maximum possible particle numbers can be released at the same time.We are assuming three particles for the largest particle size.
The particle velocity magnitude and trajectory dispersion for three different particle sizes, i.e., 1 µm, 10 µm and 40 µm, are shown in Figure 4. Large datasets of particle trajectories are visualized as a phase portrait to plot the particle position on the X-Y axes and colors represent the response time of the particle velocity.Fine particles have a shorter response time, and therefore, a smaller dispersion area.By default, the position of the particle is taken as the distance from the origin (0, 0).
As can be seen in Figure 4, the dispersion of 10 µm particles is about 3.5 mm.For comparison, the already-existing results of cold spray nozzles with a centered particle beam have a dispersion from 3.5 to 7.0 mm for 20 µm particles [14].Here, the researchers developed a nozzle with a centered particle beam; however, the particles enter the nozzle in the convergent zone.As a result, particles reach a higher velocity (471 m/s for 20 µm particles) and can cause erosion at the nozzle throat.Although the nozzle proposed in our paper has a lower particle beam velocity (360 m/s for 10 µm particles), it has a smaller particle trajectory dispersion.Additionally, due to the particles entering the nozzle in the divergent zone, the throat erosion is eliminated.
The developed fluid model was applied for the evaluation of particle impact velocities.Particle tracing trajectory model was coupled with a high Mach number flow interface.The results of CFD simulation allowed to generate required velocities for differently sized powder particles.Figure 5 presents particle velocity before contact with substrate obtained for particles with different diameters (1,5,10,20,30,40,70, and 100 µm).

Particle-Substrate Impact Analysis
Deformation behavior of particle impacting substrate is a very important point of the laser-assisted metal deposition technology process.For this reason, we provide a series of numerical experiments with copper particles impacting copper substrate with prescribed velocities, as shown in Figure 5.
A characteristic feature of the contact analysis is that contact between the micrometerscale particles and the flat metal surfaces is investigated at high velocities ranging from 229 to 668 m/s, resulting in characteristic strain rates from 2.3 × 10 6 s −1 to 6.68 × 10 8 s −1 , respectively.
A numerical particle-substrate impact experiment is conducted in the following manner.The micron-sized copper particle is represented by the sphere of radius r striking the plane surface of the solid copper substrate with impact velocity v im .Simulation of a particle impacting the substrate is formulated as an elastic-plastic interaction of two solid bodies under high velocities.

Particle-Substrate Impact Analysis
Deformation behavior of particle impacting substrate is a very important point of the laser-assisted metal deposition technology process.For this reason, we provide a series of numerical experiments with copper particles impacting copper substrate with prescribed velocities, as shown in Figure 5.
A characteristic feature of the contact analysis is that contact between the micrometer-scale particles and the flat metal surfaces is investigated at high velocities ranging from 229 to 668 m/s, resulting in characteristic strain rates from 2.3 × 10 6 s −1 to 6.68 × 10 8 s −1 , respectively.
A numerical particle-substrate impact experiment is conducted in the following manner.The micron-sized copper particle is represented by the sphere of radius r striking the plane surface of the solid copper substrate with impact velocity vim.Simulation of a particle impacting the substrate is formulated as an elastic-plastic interaction of two solid bodies under high velocities.
The yielding of metals is defined by the von Mises yield criterion limiting stress values, where the stress level is defined by a material constant-the yield stress  .In case of high velocities, yield stress is a quantity dependent not only on strain, but also on strain rate and temperature.A review of common plasticity models can be found in [26].The emphasis on strain rate-dependent models is given in review paper [27].
The plasticity model is the product of three functions expressing strain hardening, strain rate hardening, and thermal softening.The first function  () reflects isotropic strain hardening, referred to as the Ludwik plasticity model [28].It represents a classical yield criterion and describes an increase of initial yield strength  due to plastic strain: where B is the strain hardening coefficient,  is the equivalent plastic strain, and n is the strain hardening exponent.
For the evaluation of the strain rate, Cowper-Symonds [29] and Johnson-Cook [30] models are commonly used.The increasing contribution of the strain rate is better described by exponential model; therefore, the Cowper-Symonds model was used: where D is the strain rate hardening coefficient,  is the strain rate,  is the reference strain rate, and k refers to the strain rate hardening exponent.The yielding of metals is defined by the von Mises yield criterion limiting stress values, where the stress level is defined by a material constant-the yield stress σ y .In case of high velocities, yield stress is a quantity dependent not only on strain, but also on strain rate and temperature.A review of common plasticity models can be found in [26].The emphasis on strain rate-dependent models is given in review paper [27].
The plasticity model is the product of three functions expressing strain hardening, strain rate hardening, and thermal softening.The first function σ 0s (ε) reflects isotropic strain hardening, referred to as the Ludwik plasticity model [28].It represents a classical yield criterion and describes an increase of initial yield strength σ 0 due to plastic strain: where B is the strain hardening coefficient, ε is the equivalent plastic strain, and n is the strain hardening exponent.For the evaluation of the strain rate, Cowper-Symonds [29] and Johnson-Cook [30] models are commonly used.The increasing contribution of the strain rate is better described by exponential model; therefore, the Cowper-Symonds model was used: where D is the strain rate hardening coefficient, .ε is the strain rate, .ε 0 is the reference strain rate, and k refers to the strain rate hardening exponent.
An empirical model of Johnson-Cook [30] was used in order to predict the dependence of the material yielding on the material temperature: where m refers to the temperature softening exponent, and T * is a relative temperature change calculated as: where T is material temperature, T melt is the melting temperature, T ref is the reference/room temperature.
Combining the particular models defined by Equations ( 9)-( 12) leads to the final expression of the yield criterion: The relationship in Equation ( 13) expresses a rather general variation of the yield strength that is further applied in the contact analysis.
The mathematical model of the impact problem combines a dynamic equation of motion and heat transfer.The coupled non-linear thermo-mechanical problem is formulated on the basis of the Lagrangian approach.Coupling is governed by the bilateral contribution; additional heat due to plastic dissipation contributes to the heat balance, while the change of mechanical properties due to temperature shifts contributes to the dynamic behavior.Previous research has shown that the fraction η of dissipate energy converted to heat reaches 0.9 in the case of a plastically deformed system [31].
The above-described thermomechanical model was implemented in the COMSOL Multiphysics finite element software.During the simulation, the analysis of the normal contact begins at zero time instant of the local time scale with the impact velocity v imp .It continues to follow numerical motion and deformation of the particle and substrate during the loading until maximum indentation.As a result, the state of the deformed metallic particle and substrate undergoes changes through elastic-plastic deformation until full yielding.Simulation ends after unloading until full detachment.
The copper material parameter values for Ludwik and Johnson-Cook models were selected from particle impact simulations of [32,33].Full list of material parameters used in this paper are presented in Table 4.
The impact problem is described in Cartesian coordinates, the origin of which is located in the center of a spherical particle.The particle moves along the vertical axis Oy and impacts the surface.Since normal impact is considered, particle rotation is ignored and 2D axisymmetric model is used to simplify the simulation of the particle impacting the wall in the normal direction.A mesh consisting of 12k triangular quadratic serendipity elements was used for the simulation.

Coefficient of Restitution
The main results of numerical simulations are characterized by the coefficient of restitution (COR) and substrate indentation.The COR is defined as a ratio between the rebound and impact velocities (here, the average velocity of particle mass after rebound obtained by CFD simulation is considered as the rebound velocity).The energetic interpretation of COR is used to define the loss of energy during the impact.The experimental results presented in [34] show that the decreasing of COR below 0.1 leads to the particle sticking to the substrate.Our obtained variation of COR of different diameter particles is shown in Figure 6.COR has the highest value for the particle diameter of 30 µm.Smaller particles have a higher velocity, which reduces COR due to the increase in kinetic energy of the impact.Particles larger than 30 µm reduce COR due to their increased mass.
Nondimensional indentation is defined as a ratio of the indentation of substrate and the impacting particle diameter.It shows the scalability of the model for differently sized particles.The main advantage of using this parameter is that it is an experimentally observable parameter independent from particle sticking or rebound [35].Nondimensional indentation shown in Figure 7 closely resembles the particle velocity graph (Figure 5).However, the difference between these graphs is that while particle velocity slightly decreases with particle size, while indentation slightly increases for particle sizes above 40 µm.
Mathematics 2021, 9, x FOR PEER REVIEW 13 of 17 Reference temperature, K Tr 293 The impact problem is described in Cartesian coordinates, the origin of which is located in the center of a spherical particle.The particle moves along the vertical axis Oy and impacts the surface.Since normal impact is considered, particle rotation is ignored and 2D axisymmetric model is used to simplify the simulation of the particle impacting the wall in the normal direction.A mesh consisting of 12k triangular quadratic serendipity elements was used for the simulation.

Coefficient of Restitution
The main results of numerical simulations are characterized by the coefficient of restitution (COR) and substrate indentation.The COR is defined as a ratio between the rebound and impact velocities (here, the average velocity of particle mass after rebound obtained by CFD simulation is considered as the rebound velocity).The energetic interpretation of COR is used to define the loss of energy during the impact.The experimental results presented in [34] show that the decreasing of COR below 0.1 leads to the particle sticking to the substrate.Our obtained variation of COR of different diameter particles is shown in Figure 6.COR has the highest value for the particle diameter of 30 μm.Smaller particles have a higher velocity, which reduces COR due to the increase in kinetic energy of the impact.Particles larger than 30 μm reduce COR due to their increased mass.Nondimensional indentation is defined as a ratio of the indentation of substrate and the impacting particle diameter.It shows the scalability of the model for differently sized particles.The main advantage of using this parameter is that it is an experimentally observable parameter independent from particle sticking or rebound [35].Nondimensional indentation shown in Figure 7 closely resembles the particle velocity graph (Figure 5).However, the difference between these graphs is that while particle velocity slightly decreases with particle size, while indentation slightly increases for particle sizes above 40 μm.The graphs plotted in Figures 6 and 7 allow to distinguish two regions of particle size.The first region varying between 1 to 30 µm illustrates strong variation of indentation dependent on the particle size.The second region varying between 30 to 100 µm is practically insensitive to particle diameter.Consequently, the three impacts of particle diameter, 1, 30, and 100 µm, representing the region boundaries, are further discussed.
The color plot of the induced plastic strain is shown in Figure 8, whereas the resulting temperatures are presented in Figure 9.For the first region, with an increase of particle diameter, a considerable reduction of strain is observed (Figure 8a), whereas the variation of plastic strain in the second region is indistinguishable (Figure 8b,c).The graphs plotted in Figures 6 and 7 allow to distinguish two regions of particle size.The first region varying between 1 to 30 μm illustrates strong variation of indentation dependent on the particle size.The second region varying between 30 to 100 μm is practically insensitive to particle diameter.Consequently, the three impacts of particle diameter, 1, 30, and 100 μm, representing the region boundaries, are further discussed.
The color plot of the induced plastic strain is shown in Figure 8, whereas the resulting temperatures are presented in Figure 9.For the first region, with an increase of particle diameter, a considerable reduction of strain is observed (Figure 8a), whereas the variation of plastic strain in the second region is indistinguishable (Figure 8b,c).The graphs plotted in Figures 6 and 7 allow to distinguish two regions of particle size.The first region varying between 1 to 30 μm illustrates strong variation of indentation dependent on the particle size.The second region varying between 30 to 100 μm is practically insensitive to particle diameter.Consequently, the three impacts of particle diameter, 1, 30, and 100 μm, representing the region boundaries, are further discussed.
The color plot of the induced plastic strain is shown in Figure 8, whereas the resulting temperatures are presented in Figure 9.For the first region, with an increase of particle diameter, a considerable reduction of strain is observed (Figure 8a), whereas the variation of plastic strain in the second region is indistinguishable (Figure 8b,c).The plot shown in Figure 9 illustrates the particle size-dependent temperature distribution.The temperature in small particles has enough time to permeate the material during contact, while larger particles display more localized temperatures that closely resembles strain distribution.
The obtained results were compared to experimental and simulation results presented by [34].Our results indicate stiffer response, characterized by higher strain ratedependent yield strength values.This difference may be explained by using different experimental data for the plasticity model calibration.Our model data are based on the experimental data of Hassani [36], while in Razavipour [34], lower values of yield strength are reported as compared to Tong [37] and Meyers [38].

Conclusions
Based on the numerical simulation results, the new suggested laser-assisted metal deposition system potential abilities of particle deposition, the following conclusions are formulated: 1.A comparison of the geometry of Bicubic and Witoszynski nozzles showed that the nozzle with the Bicubic curve geometry reaches a higher maximum velocity at the nozzle outlet.Moreover, Bicubic curve geometry has less apparent velocity fluctuations than the nozzle with the Witoszynski curve geometry.For these reasons, the nozzle with the Bicubic curve was most suitable for our further study.2. It was found that the geometry of the nozzle described by the Bicubic curve allows to reach the higher maximum velocity at the nozzle outlet compared to that of the Witoszynski curve.Moreover, the bicubic curve geometry has less apparent velocity fluctuations than the nozzle with the Witoszynski geometry.3. The results of CFD simulation allowed us to evaluate limitations of the designed nozzle to generate required deposition velocities for differently sized powder particles.4. The results obtained revealed that small particles have the highest velocity and the lowest total force at the outlet of the nozzle.Fine particles have a shorter response time, and therefore, a smaller dispersion area. 5.It was determined that particles with a diameter limited up to 3 μm are able to reach experimentally obtained critical velocity without additional heating.The remainder of particle sizes requires additional energy for sticking.These impact velocity results were further applied by simulating the particle-substrate interaction.It was found that for particle sizes above 10 μm, additional heating is needed; therefore, future investigation of laser heating is necessary.6.The COR, indentation, and temperature obtained by the elasto-plastic particle-substrate impact analysis was used to evaluate the potential of adhesion of the different The plot shown in Figure 9 illustrates the particle size-dependent temperature distribution.The temperature in small particles has enough time to permeate the material during contact, while larger particles display more localized temperatures that closely resembles strain distribution.
The obtained results were compared to experimental and simulation results presented by [34].Our results indicate stiffer response, characterized by higher strain rate-dependent yield strength values.This difference may be explained by using different experimental data for the plasticity model calibration.Our model data are based on the experimental data of Hassani [36], while in Razavipour [34], lower values of yield strength are reported as compared to Tong [37] and Meyers [38].

Conclusions
Based on the numerical simulation results, the new suggested laser-assisted metal deposition system potential abilities of particle deposition, the following conclusions are formulated: 1.
A comparison of the geometry of Bicubic and Witoszynski nozzles showed that the nozzle with the Bicubic curve geometry reaches a higher maximum velocity at the nozzle outlet.Moreover, Bicubic curve geometry has less apparent velocity fluctuations than the nozzle with the Witoszynski curve geometry.For these reasons, the nozzle with the Bicubic curve was most suitable for our further study.

2.
It was found that the geometry of the nozzle described by the Bicubic curve allows to reach the higher maximum velocity at the nozzle outlet compared to that of the Witoszynski curve.Moreover, the bicubic curve geometry has less apparent velocity fluctuations than the nozzle with the Witoszynski geometry.

3.
The results of CFD simulation allowed us to evaluate limitations of the designed nozzle to generate required deposition velocities for differently sized powder particles.4.
The results obtained revealed that small particles have the highest velocity and the lowest total force at the outlet of the nozzle.Fine particles have a shorter response time, and therefore, a smaller dispersion area.

5.
It was determined that particles with a diameter limited up to 3 µm are able to reach experimentally obtained critical velocity without additional heating.The remainder of particle sizes requires additional energy for sticking.These impact velocity results were further applied by simulating the particle-substrate interaction.It was found that for particle sizes above 10 µm, additional heating is needed; therefore, future investigation of laser heating is necessary.6.
The COR, indentation, and temperature obtained by the elasto-plastic particle-substrate impact analysis was used to evaluate the potential of adhesion of the different powder
Figure 2 qualitatively illustrates the influence of different parameters of Witoszynski and Bicubic nozzle geometry on gas flow pressure and velocity.It can be observed from the graphs that Witoszynski and Bicubic geometries react differently to the change of parameters.Each variable is presented in Figure 2 by a family of curves.

Figure 2 .
Figure 2. Numerical calculation results using different parameters listed in Table 2. (A) gas flow pressure and (B) gas flow velocity for Witoszynski curve nozzle geometry; (C) gas flow pressure and (D) gas flow velocity for Bicubic curve nozzle geometry.

Figure 2 .
Figure 2. Numerical calculation results using different parameters listed in Table 2. (A) gas flow pressure and (B) gas flow velocity for Witoszynski curve nozzle geometry; (C) gas flow pressure and (D) gas flow velocity for Bicubic curve nozzle geometry.

Figure 3 .
Figure 3.Comparison of the velocity (A) and pressure (B) of Bicubic and Witoszynski curve nozzles.

Figure 3 .
Figure 3.Comparison of the velocity (A) and pressure (B) of Bicubic and Witoszynski curve nozzles.

Figure 5 .
Figure 5. Particle impact velocity for different particle sizes.

Figure 5 .
Figure 5. Particle impact velocity for different particle sizes.

Figure 6 .
Figure 6.Variation of coefficient of restitution against particle size.

Figure 6 .
Figure 6.Variation of coefficient of restitution against particle size.

Figure 7 .
Figure 7. Variation of nondimensional indentations h/d against particle size.

Figure 7 .
Figure 7. Variation of nondimensional indentations h/d against particle size.

Figure 7 .
Figure 7. Variation of nondimensional indentations h/d against particle size.

Table 2 .
Parameters used for numerical calculations.