A parametric simulation model for HVOF coating thickness control

High-velocity oxygen-fuel (HVOF) thermal spraying is a coating process involving multidisciplinary aspects, e.g., fuel–oxidant combustion, flame–particle jet, particle deposition, mass and heat transfer, and even robotic kinematics. Like most coating processes, in HVOF processes, coating thickness is a significant property determining the coating performance; hence, this property should be accurately controlled during the process. In view of green, smart, and digital manufacturing, the coating thickness prediction model is demanded to produce high-quality coatings efficiently. This paper presents an approach to parametrically simulate the coating thickness in HVOF processes via an integrated numerical model. Firstly, an axisymmetric computational fluid dynamics (CFD) model is constructed to compute the behaviors of the fuel–oxidant combustion, flame–particle jet, and particle deposition distribution. Secondly, based on the particle distribution in a 2D axisymmetric model, a 3D single coating thickness profile model is developed by constructing a circular pattern using the axis of the nozzle. Further, this profile is smoothened by a Gaussian model, and its mathematical expression is obtained. Finally, a numerical model couples spray paths with the mathematical expression to model the coating thickness distribution on a substrate surface under industrial scenarios. At the end of this paper, to verify the proposed model’s effectiveness, four sets of operating parameters with a single straight path were experimentally implemented. The width and height of the bead-like shape coating were in good agreement with the simulated results. The normalized root-mean-square errors of the cross-sectional profile heights were around 10%.


Introduction
High-velocity oxygen-fuel (HVOF) deposition is a thermal spraying technology used to coat a protecting layer to improve the workpiece surface performance [1]. It is a multidisciplinary process involving chemical reaction, turbulence, compressible flow, multiphase interactions, subsonic/supersonic transitions, and droplet deformation [2]. During the process, the thoroughly mixed fuel-oxygen gases (typically hydrocarbon-oxygen) and microparticles are fed into the gun chamber where a combustion reaction takes place which in turn generates a hot-sonic/supersonic multiphase gas stream. Within the stream, microparticles of metals, alloys, and/or ceramics are accelerated, heated, and ultimately deposited onto a substrate at high speed [3,4]. When the particles hit a substrate, the sudden deceleration causes a pressure buildup at the particle-substrate interface. The high energy inside the particles forces the molten material to flow laterally or the ductile solid material to deform. The liquid bonding material spreads outward from the point of impact and forms a splat [5]. According to the studies of the coating formation mechanism [3,6,7], the coating formation is predominantly determined by the in-flight particle behavior near the substrate surface, which in turn is affected by the operating parameters, such as the gas flow rate, the fuel/oxygen ratio, and other various operational conditions, such as the spray systems used, substrate shape, and the fuel/oxidant types [7,8]. In contemporary manufacturing, robots and handling systems have been developed to control the movement of the spray torch to deposit the coating on substrate surfaces. The macrostructure of the coating is affected by the parameters relevant to the movement, e.g., scanning speed, scanning step, and spray distance [9,10].
Like most coating processes, the coating thickness is a critical property that should be precisely controlled during the process [11,12]. For the purpose of obtaining a homogeneous coating, the coating thickness modeling is always a hot spot to control this property under various industrial scenarios. In the aspect of HVOF thermal spraying, the modeling is still highly challenging due to its complex system and multidisciplinary process parameters. Its coating thickness models, especially for studying the specific HVOF parameters, are still in relative infancy. Most research works were devoted to the kinematic parameters' effects, e.g., scanning speed, scanning step, and spray distance, on the coating thickness [11,13]. The operating parameters and working conditions specifically used in the HVOF coating, e.g., the gas flow rates and the nozzle geometry, are also critical factors determining the coating quality [14]; however, they have not been considered in the previous research works. Moreover, most of the previous models were empirical models derived from the experimental data. Although these models are reliable for specific applications, they may not be applicable to other systems and working conditions, and these modeling methods based on experiments are time-consuming and expensive. Therefore, this paper proposes a parametric model for HVOF coating thickness control with due regard to the complex spray system and the multidisciplinary parameters by using numerical modeling.
The remainder of this paper is organized as follows: Section 2 reviews the research works related to the coating thickness modeling in thermal spray processes. Section 3 introduces and demonstrates the procedure of our proposed model by simulating a specific scenario. In Section 4, the simulation results and experimental validations under four sets of operating parameters are presented to verify the effectiveness of the model.

Related work
As this paper aims to propose a parametric model for constructing the coating thickness that is specific to HVOF processes, the review of the related works is divided into two parts: numerical modeling of HVOF thermal spraying and coating thickness modeling. In general, the numerical modeling of the HVOF process mostly focused on the simulation of the flame behaviors and in-flight particle characteristics. The works that integrate the numerical models into the coating thickness simulation are very limited. Regarding previous research efforts on the coating thickness modeling, most works were dedicated to the kinematic parameters in robotized processes for thermal spray technologies, e.g., scanning velocity, scanning step, and spray distance, instead of the parameters in specific thermal spray technology, e.g., gas flow rates in HVOF processes.

Numerical modeling of HVOF process
The numerical study of the HVOF process mainly focuses on understanding the physicochemical phenomenon in the process, such as the combustion reaction, gas-flow dynamics, and particle in-flight behavior [8]. In recent years, computational fluid dynamics (CFD) simulation, as an important composition of CAE tools, has reached a high level of sophistication, which makes the parametric integrated study of the HVOF process feasible.
The numerical modeling of HVOF thermal spraying was rapidly developed in the early years of the 21st century with the increasing maturity of the CFD technique. Given a large number of parameters or factors in HVOF processes, the related works investigated four significant aspects affecting the coating quality, nozzle design, combustion reaction, powder particle, and substrate component. Baik et al. [15] created a 2D comprehensive CFD model to explore the effect of nozzle geometries on the performance of the flame by using parametric analysis. By using a similar computational model, Khan and Shamim [16] further explored the effect of nozzle geometries and investigated the physicochemical performance in a dual-stage HVOF system. Emami et al. [17] carried out the research works on the dual-stage HVOF system and discussed the effects of two combustion models on the temperature and velocity fields in this system. Li and Christofides [2] reported a hybrid multiscale model for the process and revealed the effects of chemical reactions on the in-flight behaviors of the particle-gas jet. Their models were both propylene-oxygen combustion and 2D simulation. Tabbara and Gu [18] constructed a 2D kerosene-fueled HVOF model to examine the effects of the liquid fuel droplets on the thermodynamics of the combustion gas flow. Pan et al. [19] also built a 2D liquid-fueled HVOF model but with a kerosene-fueled chemical reaction and analyzed the effects of particle properties on their in-flight behavior. More recently, some researchers concentrated on the numerical modeling of the more advanced suspension HVOF (SHVOF) process. Jadidi et al. [20] built a 3D two-way coupled Eulerian-Lagrangian model to analyze the in-flight behavior of suspension particles in their early publications. Then they further explored the effect of the substrate shape on the particle behavior and analyzed the catch rate for different substrate shapes and standoff distances [7].
In general, most previous research works focused on the numerical analysis of the physicochemical phenomenon in the HVOF from different aspects. The application of these models in the coating thickness simulation is still restricted.

Coating thickness model in robotized thermal spray processes
At the beginning of robotized thermal spray process modeling, some researchers worked on the spray path generation for robotized thermal spray coating application and the modeling method linking the path, thermal history of the workpiece, and coating thickness prediction. Candel and Gadow [21] developed a computer-aided trajectory generation method by using reverse engineering technology. Then they constructed a linking method between the impinging flame and the substrate physics performance by integrating the generated path into the finite element method (FEM). The influence of the generated path on the substrate temperature was further analyzed. Although their method might be applied to the coating thickness simulation, the validation was not involved in their work. Bolot et al. [22] proposed a similar linking model between the nozzle path and the dynamic behavior of the coating (thermal history and coating thickness) by using the user-defined function (UDF) within ANSYS FLUENT [23]. However, their spray patterns for thermal flux and coating thickness were estimated by mathematical derivations based on a specific vacuum plasma spraying system instead of a general approach, which limits the model's applicability to different systems. Lately, Ma et al. [24,25] proposed an iterative coupling method to enhance the accuracy of simulating the dynamic temperature of the substrate. The mass transfer might be applied to their method, but it was not studied in their work.
As the robotized thermal spray process modeling getting sophisticated, the coating thickness prediction and optimization have been brought to the fore. Based on a simulation model, Hegels et al. [26] concentrated on optimizing a given nozzle path to improve the coating thickness uniformity on a workpiece with flat and curved regions through an iterative post-optimization. The 3D coating height profile used in their simulation was analytically formulated, which did not consider some specific operating parameters, such as electric arc current in the plasma spraying or fuel gas flow rate in HVOF thermal spraying. Chen et al. [11] constructed a numerical coating profile model through a Gaussian distribution fitting. Chen's model is specific to the cold spray process and focuses on the kinematic parameters, such as spray angle, spray distance, and nozzle traverse speed. Zhang et al. [10] created a spherical surface coating thickness model. A spherical coordinate system was employed to describe the Gaussian coating profile which required a preliminary experiment to estimate the height of a single coating profile. Their model was validated by a specific set of parameters in the plasma spraying. More recently, Deng et al. [27] reported an approach predicting coating thickness for components with complex curved surfaces, especially in the case of shadow effects. They further developed the Gaussian coating profile model concerning the parameters of spray angle, spray distance, nozzle traverse speed, and deposition efficiency. However, since the data used for constructing the Gaussian coating profile was obtained in cold spraying, their model is limited to this technology. As to the aspect of HVOF thermal spraying, the coating thickness prediction is still in its relative infancy due to the complexity. Li and Christofides [28] presented a stochastic simulation method to predict the microstructure of HVOF coatings. Besides, based on experimental validation, Mostaghimi et al. [29] proposed a set of rules to simulate the buildup of a small area of the coating. However, the coating thickness simulation was not shown in their papers.
Based on the abovementioned findings, most of the coating thickness models mainly concentrated on the effect of the kinematic parameters on the development of the coating thickness. Some operating parameters and working conditions in specific thermal spray technologies, e.g., the gas flow rates and the nozzle geometry in HVOF processes, are not involved.
To sum up, even though the numerical studies of HVOF thermal spraying are becoming increasingly mature, most previous research works have only aimed to understand the fundamentals of the complex phenomenon. The applications of these models are still restricted. To the authors' best knowledge, there has been no published work on linking the physical phenomenon models with the coating thickness model so that the specific operating parameters in the HVOF and the kinematic parameters in a robotized thermal spraying process could be integrally studied in a comprehensive simulation.

Proposed model
The general procedure for constructing the model is presented in Fig. 1. The first step is to simulate the in-flight flame/ particle behavior by using a sophisticated CFD model. Due to the revolved structure of the nozzle chamber, a 2D axisymmetric CFD model is always used to reduce the computation cost. Once the steady behaviors of the combustion, flame, and in-flight particles are obtained, the particle landing distribution on the substrate surface can be extracted from the CFD model. To evaluate the influence of particle size on the deposition distribution, by suiting the Rosin-Rammler diameter distribution [30,31], the particle size range is discretized into a number of groups that are separately inputted to the CFD model and get the corresponding particle distribution at the impact on the substrate. Then, according to the particle landing distribution, a rule-based coating growth model captures an initial coating profile in the second step. The third step is to convert this initial coating profile built on the 2D CFD model to a 3D single coating profile and fit the converted profile by a Gaussian distribution. The last step is to apply the Gaussian distribution and simulate the coating thickness with a spray path through a proposed linking method. The details of each step are expressed in the following subsections.

Parameterized HVOF in-flight behavior model construction
As the coating formation process in the HVOF is determined by the in-flight particle behavior prior to the impact on the substrate surface, which, in turn, is affected by the combustion, flame, and in-flight particle behavior, these physicochemical phenomenon simulations are noticeable for the coating thickness modeling. With 20-year development, these simulations by using CFD tools, intuitively named "HVOF in-flight behavior modeling" in our work, have been mature enough. This subsection aims to develop an approach to systematically integrate these previous efforts into the coating thickness prediction. The explanation of the CFD models and governing equations for the in-flight behavior are not the core of this work; therefore, they are summarized in the Appendix and employed in this work to construct a parameterized HVOF in-flight behavior model so that the proposed idea is applicable to various HVOF applications. The models and governing equations in the Appendix can be solved by ANSYS FLUENT [2,32].
Based on the investigation of the previous works, the inflight behavior influencing the particle landing distribution is controlled by the operating parameters, such as the gas flow rates, spray distance, and other operational conditions, such as the nozzle configuration and particle properties. Thus, for predicting the coating thickness profile, the parameterized HVOF in-flight behavior modeling needs to respect two aspects: the specific operating parameters and operational conditions in HVOF processes, as shown in Fig. 2. It should also contain the process of computational domain determination, Fig. 1 The procedure for simulating the coating thickness in HVOF processes reaction formula calculation, solver setup, and convergence analysis. Figure 2 starts from the construction of the computational domain for the in-flight behavior. It can be extracted from the nozzle chamber geometry via Boolean operations and in combination with the spray distance. For the purpose of demonstration, a homemade Diamond Jet Hybrid gun as shown in Fig. 3 and a set of operating parameters indicated in Table 1 (recommended by the manufacturer) are used here to get the results of each step. For example, a computational domain as shown in Fig. 4a is built according to the building elements and geometry parameters of the Diamond Jet Hybrid gun in Fig. 3. The domain from the nozzle exit to the substrate surface is determined by a spray distance given in Table 1. To exactly capture the behavior, the mesh of the important local parts is refined as shown in Fig. 4b.
Typically, the combustion process is fueled by hydrocarbon gases, such as propylene, propane, and many others, or (This method is adapted on the basis of our previous research work [33] and specifically designed for the coating thickness simulation) hydrocarbon liquids as kerosene. Oxygen and compressed air are usually used as the oxidant. When the combustion temperature reaches above 2000 K, the combustion of hydrocarbon gases is complex because the reaction products will be dissociated into a number of low-molecular species [34]. To model the combustion model accurately, the reaction is considered as the following form: Then the combustion process is solved by the eddy dissipation model [2]. The physics models for calculating the flame are set up according to the parameters in the Appendix, and the operating parameters are inputted via the boundary conditions tagged in Fig. 4a. It is worth mentioning that the reaction formula is derived from an instantaneous equilibrium code [35] which needs the equivalence ratio and an assumed reaction chamber pressure. However, the chamber pressure from the simulation result is very likely different from the preliminary assumed pressure. Therefore, to improve the reaction model accuracy, several iterations of running this part of the algorithm are needed to ensure the consistency between these two parameters up to a certain tolerance range. After the convergence state check and grid independence analysis, the steady flame properties, such as velocity, pressure, temperature, and momentum, can be obtained. For example, Fig. 5 presents the velocity magnitude contour and pressure contour of the flame flow with the operating parameters given in Table 1.
Owing to the very low particle loading (less 4% usually) [32], a one-way coupling between the gas phase and the particle phase is assumed; in other words, the momentum and heat of the particle are solved by the Lagrangian approach after the gas flow fields are determined, and the particles have no influence on the gas phase [6,36]. It should be noted that this model for computing the in-flight particle behavior was widely used in previous studies and partial particle properties have been tested against experimental results [2,20,32]. For example, the simulated particle velocity and temperature at impact in Li and Christofides's work [8,32] are consistent with experimental observations of Zhao et al. [37]. In a more recent publication, this algorithm, computing particle behavior after solving the gas phase, was further tested in a more advanced HVOF technology, suspension HVOF [20]. In this publication, the simulated particle velocity and temperature at different standoff distances were in good agreement with the experiment, and the variation trend was similar. Collectively, this model is adequate for simulating particle behavior.
Therefore, in our work, once the flame behavior simulation reaches a certain stable state, powder particles can be inserted into the model through the corresponding boundary condition as shown in Fig. 4a and released from the grid nodes. The length of releasing particles along the feedstock inlet boundary shown in Fig. 4a and the initial particle velocity are carefully set up according to considerations of the real scenario and may need to be reasonably tuned to provide results for important characteristics of the coating profile (such as the width of the single coating profile) that are consistent with the available experimental results. For the homemade HVOF (This homemade Diamond Jet Hybrid gun was also used in our previous research work. The schematic representing the gun structure is adapted from Ren et al. [33]) system used in this work, this length is set up as 3.5 mm from the axis of the gun along the radial direction, and the particles are released from the feedstock inlet boundary with an initial velocity of 0.01 m/s. The powder particles should be evenly dispersed on the area of the feedstock inlet. Their thermophysical properties and shape equivalent factors can be correspondingly assigned to the parameters shown in the Appendix. In this work, WC-12Co powder particles are used to demonstrate the modeling process. Table 2 shows their thermophysical properties [38]. Figure 6 is a scanning electron microscope (SEM) image of the WC-12Co powder particles, which indicates the spherical particles should have a shape factor of 1. Apart from these parameters, particle size distribution is a critical factor affecting the coating profile [39]. The following subsection introduces the particle size discretization and gives a demonstration of the in-flight particle behavior simulation.

Particle size distribution and the in-flight particle behavior simulation
Because the size distribution of the powder particles used in HVOF thermal spraying is polydisperse, and particles of different sizes have different dynamic behavior during the flight due to different momentum, the particle size distribution [32] should be taken into account for accurately calculating the particle deposition distribution on the substrate. The Rosin-Rammler distribution of cumulative mass [40] widely exists in the literature and industry for describing the powder particle size [30,32]. To simulate this effect, our approach discretizes the particle size range into numerous groups of different sizes with the same number of particles firstly. Then these groups are separately inserted into the stable flame simulation for computing the in-flight trajectory of particles of different sizes. Meanwhile, the Rosin-Rammler distribution of cumulative mass is converted to the count distribution of particles [40], which will be introduced in the next subsection to simulate an initial coating thickness by combining a rule-based coating growth model. For further demonstrating this step, taking the WC-12Co powder particles as an example, the particle size is from 7.5 to 50 μm and follows a Rosin-Rammler distribution of cumulative mass as shown in Fig. 7. The model computes the trajectories of 30 groups of different particle diameters. The discretized particle sizes are represented by the orange circles in Fig. 7. The star signs are the measurement result from a laser diffraction analyzer, which validates this Rosin-Rammler distribution. Figure 8 shows the trajectories of the partial particles under the aforementioned flame behavior. Please note that, for including the effect of turbulent velocity fluctuations on the particle trajectories, the  particle motion in the turbulent flow is calculated in a stochastic tracking approach [30], and up to 15 independent stochastic particle tracking calculations are performed to obtain the statistical description of the particle trajectories in the gas flame. Therefore, in the simulation, 480 particles per group are released on the boundary condition and computed in the gas flame (this number also relevant to the node number at the particle-insert boundary condition). Apart from the particle trajectories, other in-flight particle parameters with the consideration of the particle size distribution, e.g., temperature, melting ratio, and flight time, can be computed according to the equations shown in the Appendix.

Splat morphology approximation and coating growth
Once the particle trajectory calculations are completed, the particle physics behavior and distribution at impact (on the substrate surface as shown in Fig. 4a) can be extracted. Then, based on the extracted data, an initial coating thickness on the area in the red rectangle shown in Fig. 9 can be computed by modeling the coating growth.
However, coating growth is a complicated process. Generally, thermally sprayed coatings mainly consist of lamellar splats interspersed with a small number of pores and unmelted particles as shown in Fig. 10b. The splats are the major building blocks of the coating, which are formed by the impact, deformation, spreading, and solidification of the individual powder particles [28,41]. Each particle is accelerated and heated in the flame during the process, after which it impinges on the substrate or previously deposited coating layer in a molten, partially molten, or solid state. At impact, the sudden stagnation of the particles leads to a pressure buildup and deformation occurs. Then the molten part of particles spreads outward from the point of impact and forms a splat [28,33]. Figure 10a gives a scanning electron microscope (SEM) image of a single splat on a flat substrate surface. The pores as shown in Fig. 10b are deformed by the interaction of the new deposited particles and the previously deposited coating layers. The coating growth process can be treated as the pileups of these compositions as depicted in Fig. 10c.
It is obvious that the complete virtual model with the consideration of all physics phenomena involved is hard to be realized. Moreover, modeling the impact and solidification of just a single drop even requires many hours of computing time. Simulating the pileups of a huge number of particles is impractical. Since this work concentrates on the macro coating thickness profile (the radius of a single coating profile is typically several millimeters), the microlevel splat morphology (the radius of a single splat typically ranges from several micrometers to several hundred micrometers) is simplified with reasonable assumptions, which are discussed in the following paragraphs. Here, the splat morphology [42] is the final shape of a single particle on the substrate or previously deposited coating layer after the impact, spreading, and solidification. In addition, according to the previous major characteristics of the coating formation, some rules for calculating the pileups of particles are adapted in this work on the basis of previous related research works [2,8,28].
Regarding the particle deformation at impact, there are two assumptions in HVOF processes. As shown in Fig. 11a, if a particle can reach a fully melted state just prior to impact, the droplet becomes a lamellar splat as a result of deformation  Cumulative mass distribution of the WC-12Co powder particles and particle size discretization which can be treated approximately as a cylinder shape [43,44]. The splat flattening degree ξ is defined by the following equation: where D is the approximated diameter of the splat and d p is the particle diameter prior to impact. The previous research indicated that for HVOF coatings, ξ predominantly depends on the Reynolds number which represents the viscous dissipation of the inertia forces [28,43]. The flattening ratio can be approximated by a function of this parameter as follows [28,43]: If a particle is partially melted, the unmelted particle will form a hemisphere with the equivalent volume, and the melted part will form a ring around the hemisphere as depicted in Fig.  11b [28,45]. The flattening ratio of the melted part can be calculated by using the previous equation.
To simulate the pileups of splats, some rules are adapted in this work on the basis of the related research works [2,8,28] to program different events that might take place when a particle hits the substrate and previously deposited coating layer. A brief summary of these rules is given here. For a more detailed interpretation, the reader could refer to Shi et al. [28]. Fig. 12a, if a particle at impact is partially melted, the unmelted part will form a hemisphere on the previously deposited layer, and the melted part will form a Fig. 8 Trajectories of the partial particles Fig. 9 Schematic depiction of the single coating profile calculation ring around this hemisphere. Besides, if the unmelted part of the particle locates at the previously deposited layer where cavities exist, pores will form under the hemisphere. 2. If the unmelted part of the particle hits at the previously deposited layer that is formed by an unmelted particle, it will bounce off, and only the melted part will be deposited (Fig. 12b). 3. If the particle is able to cover a whole step and the splat encounters a dead end, it will break at the corner of the step as shown in Fig. 12c. 4. If the particle impinges at the corner or boundary of a step and the ratio of the height of the step to the thickness of the splat is not great, the splat will flow over the corner and outer surface as shown in Fig. 12d. 5. If the splat comes to the previously deposited layer with small gaps as shown in Fig. 12e, the splat will cover the gaps, and pores will form. 6. In the fourth rule, if the ratio is large enough, the splat will flow over the corner and form a slope at the corner as shown in Fig. 12f.

Combination of the in-flight particle behavior and the splat morphology approximation
By combining the aforementioned rules and the particle parameters at impact from the in-flight behavior model (e.g., Reynolds number, melting ratio, particle distribution), the coating growth can be calculated. Finally, an initial coating height based on the 2D numerical model with the consideration of the particle size distribution is the summation of these groups by multiplying their corresponding count distribution, as shown in Fig. 13, derived from the abovementioned cumulative mass distribution. The result of this step is corresponding to the initial coating height on the area highlighted by the red rectangle in Fig. 9.

Initial single coating profile calculation
After getting the initial coating height on the area of the red rectangle in Fig. 9, the next step is to create a circular pattern Fig. 10 The major characteristics of coating formation: a SEM image of a single splat under a set of specific operating parameters, b microstructure of a thermal spray coating, and c schematic of the coating formation Fig. 11 Deformation of a fully melted particle (a) and partially melted particle (b) upon impact on the substrate (adapted from Shi et al. [28]) of the initial coating height model using the axis of the nozzle and calculate the overlapping parts, such as the overlapping area of the red and yellow rectangles as shown in Fig. 9. It is worth noting that, for getting a continuous coating distribution, the instances in the circular pattern could be evaluated by Eqs. 4-6.

number of instances
where Length of sample can be estimated by the maximum radial distance of the particle distribution at impact and Width of sample can be calculated by using the average diameter of particles as the following equations: where CMD is the count mean diameter, d i is the median value of the particle diameter in the groups in Fig. 13, n i is the count distribution, and N is the total count of the particles.

Width of sample
where h splat is the function of the equivalent height of the splat which can be derived from Eqs. 2 and 3. Figure 14a presents the initial coating model after the circular pattern construction. Through numerical and experimental validations, for most thermal spray processes, the Gaussian distribution can be used to accurately represent the single coating profile [11,13,25]. Hence, a symmetric Gaussian distribution is applied to characterize the cross-sectional profile of the coating model, as illustrated in Fig. 14b. Equation 7 is its mathematical expression.
where φ, called the single coating profile, is the coating height of the position r meters away from the position of the center of peak (b) where the nozzle stationarily stays at a position for 1 s. Thus, the unit of φ is m/(m 2 ·s). a is the height of the curve's peak, and c is the standard deviation. For the scenario described earlier, a equals 4.271 × 10 −5 , b is 4.292 × 10 −7 , and c is 3.743 × 10 −3 . δṁ; d m ; γ is a function relevant to the powder feed rate (ṁ ) and the deposition efficiency (γ) [46]. Their estimations will be introduced in the following subsection.
In the Cartesian coordinate system, the distance away from the position of the center of the peak, r, can be calculated by the following equation: where (u x , u y ) is the center coordinate of the coating profile on a substrate surface.

Deposition efficiency determination
So far, the mass rate of the powder particles and the deposition efficiency have not been considered in the modeling. In HVOF processes, in addition to the unmelted particles bouncing off the substrate, the deposition efficiency is affected by many other parameters and factors, such as the spray distance, gas flow rates, and robot kinematic parameters. Among these parameters, the fuel flow rate is a critical one. A low fuel flow cannot heat the particle to reach a molten state, which results in deposition efficiency reduction. On the other side, a high fuel flow may overheat particles and cause particle oxidation, which also leads to a low deposition efficiency. For simplification, in this work, only this critical parameter, fuel flow rate, was studied to quantify the variation of the deposition efficiency. A series of fuel flow rates changing from 80 to 380 SLPM were implemented for measuring their deposition efficiencies. The experimental measurement results and the regression model are given in Fig. 15.

Single coating profile model
After determining the deposition efficiency under a specific gas flow rate, the number of the deposited particles can be estimated by using Eqs. 9 and 10. The unmelted particles bouncing off the substrate, which is introduced in Section 3.2.1, should not be considered as the deposited particles.
whereṄ deposited is the number of deposited particles per second andṄ total is the total number of particles fed into the nozzle per second, which can be approximated by the following equation: where ρ p is the density of powder particles,ṁ is the powder feed rate, and d m is the diameter of mass average [40], which can be derived by using the count distribution.
δṁ; γ used in Eq. 7 is a function to include the effect of the powder feed rate and the deposition efficiency on the single coating profile. Its analytical derivation is given by Eqs. 11 and 12.
where n i is the count percentage of different size particles, for example, for the WC-12Co used in this work, the count percentage of different size particles is given in Fig. 13. N group is the number of particles in the groups of discretized particle sizes.
Then, substitute Eqs. 9 and 10 into Eq. 11. δṁ; γ is given as Finally, by substituting Eq. 12 into Eq. 7, a single coating profile model can be obtained. Figure 16 shows the 3D model of the single coating profile under the scenario described earlier. Equation 13 is its mathematical expression. The algorithm is performed in the software Matlab 2020b on a 64core CPU with a utilization of 70-80%.
3.4 Integration of the spray path

Mechanism of the spray path integration
By combining the numerical model of the single coating profile with the parameters relevant to the nozzle movement, e.g., scanning speed and scanning step, it is able to build the model of coating thickness distribution on a substrate surface deposited by a nozzle path. The schematic is illustrated in Fig. 17, where the substrate surface is meshed by a grid with constant grid size and a part of an exemplary nozzle path is represented by a blue line. Firstly, the path is dispersed into a series of target points separated by a constant time step, as shown in Fig. 17a. The target points can be obtained according to the constant time step and nozzle speed interpolation. An appropriate time step value is important for the thickness simulation result due to the fact that a large time step can lead to low accuracy and a small one needs excess computation. Then the target points on the nozzle path are transformed into mapped target points on the substrate surface along with the orientation of the spray nozzle, as depicted in Fig. 17a. After the target point transformation, for a moment that the nozzle passes through a target point, the coating height at each grid node as shown in Fig. 17b can be calculated by substituting the node coordinate (x, y) into the integral in Eq. 14. Lastly, by repeating this step with the progress of the nozzle movement, coating thickness distribution on the entire substrate surface can be simulated. In this way, it can simulate the effect of the nozzle movement parameters on the coating thickness distribution, such as scanning step and scanning velocity.
where ∅ i is the coating height at the node coordinate (x, y) when the nozzle passes through the ith mapped target point (u xi , u yi ).

Validation of the mechanism
To validate the mechanism of the spray path integration, the Gaussi an distri buti on ( Eq. 13) obtained by the abovementioned operating parameters was applied to simulate the coating distribution on a flat-plane substrate as depicted in Fig. 18. The substrate surface has meshed with a 5 × 5-μm grid cell. To indicate the effect of different spray paths on the coating distribution, two sets of the kinematic parameters in Table 3 for describing the nozzle movements along a zigzag path shown in red in Fig. 18 are separately carried out. In these two simulations, the nozzle passes the substrate surface in a vertical direction and maintains the constant spray distance given in Table 1 (160 mm) and the constant scanning velocity in Table 3. The time step for the target point generation is 0.1

Experimental setup
For verifying the proposed model, simulation and experimental investigations under different operating parameters were implemented. The experiments were carried out by a homemade Diamond Jet HVOF spray system and a YASKAWA DX200 MA2010 robot system. For keeping the information consistent between the experiments and simulations, the building elements and geometries of the nozzle as depicted in the modeling demonstration ( Fig. 3) were acquired from the direct measurement of the real spray gun. Besides, the WC-12Co powder particles used in Section 3 were kept on using in the experiments as feedstock to coat a flat-plane steel substrate (size: 200 mm × 200 mm × 20 mm). To reduce the effect due to the robot kinematics as much as possible, the spray path used for coating is designed as simple as possible, which is a straight line shown in red in Fig. 21. Besides, to avoid the deformation of the coating profile due to gravity during the coating formation, the substrate was placed horizontally. The impact of gravity on the in-flight particles can be neglected due to the high velocities and high-velocity gradients [18,32]. In order to obtain an obvious bead-like shape for observation and measurement, the nozzle reciprocated up and down along the red path with a scanning velocity of 300 mm/s for 12 passes. Furthermore, the lateral distance is 100 mm, which is far away from the substrate edge, so that the robot can reach a stable motion before passing through the substrate.
The spray distance was kept constant during each experiment. Table 4 shows the detailed spray distance and operating parameters for different experiments. Apart from the set of parameters (baseline) for the model demonstration in Section 3, three additional sets were conducted to further evaluate the predictive capability of the proposed model. Two of them were conducted respectively above and below the gas flow rates of fuel and oxygen in case 1, and one is a set of parameters with a shorter spray distance.

Results and discussion
The results of experiments and simulations under the scenarios in Table 4 are introduced in this subsection. The left side from top to bottom in Fig. 22 respectively shows the top views of the coatings after 12 passes in the experiments of case 1, case 2, case 3, and case 4. The right side in Fig. 22 is their corresponding simulation results. It is known that the simulated coating thickness along the path direction is almost continuous and the same due to the very short time step (0.1 ms). Thus, only partial lengths of the coatings are displayed in Fig.  22 as the representatives of the simulation results. To conveniently compare the width of the coating profile, the origin of coordinates lies at the entry point of the nozzle on the substrate surface. The positive and negative values of the width all mean the radial direction of the coating profile. From Fig.  22a, b, it can be observed that the width of the coating bead in the experiment is around 17 mm which is very closed to that of the simulation. Similarly, in Fig. 22c-h, we can observe the  For accurately evaluating the capability of the proposed model, in the experiments, the cross section at the middle of the length of the coating beads was measured by contour measuring equipment (HOMMEL-ETAMIC nanoscan 855). These measurement results were compared with the corresponding simulation results as shown in Fig. 23. The measurement results under the 4 sets of parameters are displayed in orange in Fig. 23a-d, respectively. The blue plots in Fig. 23a, b, are their simulated cross sections. Obviously, the simulated coating profiles basically fit well with the measured ones. It also can be found that the coating thickness of case 2 is higher than that of case 1, but the width is narrower than that of case 1. This was because the higher gas flow rates in case 2 generated a flame jet with a higher axial velocity, which, in turn, accelerated the powder particles along the axial  direction. Hence, the radial displacement of powder particles was relatively reduced, and the particles tended to impact around the center of the coating profile. A similar phenomenon was found in case 3. The short spray distance in case 3 also reduced the radial displacement of the powder particles. Case 4 has a flat coating profile with the lowest peak value due to a low deposition efficiency. When the gas flow rates came to a low level, most powder particles did not reach a molten state, so they could not stick on the substrate. Besides, the low gas flow rates generated a low-power flame which could not restrain the particle to deposit around the center of the coating, so the profile tended to be flat. Further, for quantitatively assessing the capacity of the model, the mean absolute error (MAE), root-mean-square error (RMSE), and normalized root-mean-square error (NRMSE) of the cross-sectional profile heights were calculated via more than 10,000 measurement points. Table 5 shows the calculation results. The MAE, RMSE, and NRMSE are around 0.04, 0.05, and 10% for the first three cases, which is acceptable. Case 4 has a relatively large error. In general, these results and analysis indicate that the proposed model enables the prediction of the coating thickness in HVOF processes with respect to the kinematic parameters for describing the spray path and HVOF specific operating parameters. However, due to the complexity of the technology, there exist differences between the experimental and numerical results. Multiple reasons induce the differences. From the experimental point of view, the continuously high-energy jet during the coating may lead to the substrate to deformation, which affects the measurement accuracy. It can be found that due to the thermal deformation, the experimental profiles, especially in Fig. 23a, b, are not symmetric and tend to be rightward. From the simulation point of view, the particle trajectories were solved by the stochastic tracking approach under the Reynold and Favre averaging flame. The particles may tend to land around the center of the profile. Thus, it can be observed that the simulated coating profiles are narrower and higher than the experimental results. In addition, the particle  deformation at impact and the coating growth are approximately modeled with the consideration of very limited factors. However, the process is very complicated and affected by a number of factors. In the future, more factors can be consolidated in the empirical function to enhance the accuracy of the model.

Summary
This article contributes to an integrated numerical model for simulating the coating thickness in HVOF processes by linking spray paths with the Gaussian coating profile. The kinematic parameters for describing spray paths are incorporated into our model; thus, the effect of spray paths on the coating thickness can be studied. Besides, the Gaussian coating profile is obtained via an HVOF in-flight behavior model. The HVOF specific parameters, such as fuel type and gas flow rates, and working conditions, such as nozzle chamber geometry and powder particle property, are used to construct the Gaussian coating profile. Hence, the proposed model can also analyze the influence of the HVOF specific parameters and conditions on the coating distribution. Through experimental investigations, the capacity of the model has been verified in this work. The results indicate that it can flexibly simulate the coating thickness while considering different HVOF scenarios. However, the current work is still in its infancy.
The modeling of the splats' morphology and coating growth is not sophisticated. In the future, it can be captured by applying a more advanced numerical model that encapsulates the more physics features of the coating formation, such as the deformation of the previously deposited coating due to the impinging of the flame and particles. In addition, the substrate surface in the HVOF in-flight behavior model will be changed to different geometries for exploring the coating distribution on different substrate shapes. A 3D CFD model for the in-flight behavior simulation can be developed so that the spray angle of the nozzle relative to the substrate surface can be adjusted. Moreover, the comparison of simulation results of two different spray paths in Section 3.4.2 implies that the model has the potential to control and optimize the coating distribution by implementing an optimization methodology.