Efficient Simulation of Nonlinear Heat Transfer during Thermal Spraying of Complex Workpieces

The quality of coatings, produced by thermal spraying processes, considerably decreases with the occurrence of higher residual stresses, which are especially pronounced for complex workpiece geometries. To understand the occurring effects and to aid in the planning of coating processes, simulations of the highly transient energy flux of the HVOF spray gun into the substrate are of great value. In this article, a software framework for the simulation of nonlinear heat transfer during (HVOF) thermal spraying is presented. One part of this framework employs an efficient GPU-based simulation algorithm to compute the time-dependent input boundary conditions for a spray gun that moves along a complex workpiece of arbitrary shape. The other part employs a finite-element model for a rigid heat conductor adhering to the computed boundary conditions. The model is derived from the fundamental equations of continuum thermodynamics where nonlinear temperature-depending heat conduction is assumed.


Introduction
Thermal spraying is a cost-efficient coating technique for the production of wear-resistant surfaces consisting of various materials tailored to particular applications.For the coating of forming tools, e.g., hard material coatings of tungsten carbide (WC) and cobalt (Co) are used, [1], because of their superior wear-resistance compared to chrome (Cr) and nickel (Ni) coatings, [2].As a disadvantage, the High Velocity Oxygen Fuel (HVOF) thermal spraying process induces a large amount of energy into the heterogeneous coating and the substrate which leads to a complex transient thermo-mechanical problem, as illustrated in Figure 1.For an overview about HVOF thermal spraying, the reader is referred to, e.g., [3] [4].While thermal-spraying of planar work pieces delivers coatings of quite satisfying quality, the quality considerably decreases with the rising degree of complexity of the work pieces' geometry, for instance due to radii or curvatures.Experimental analyses suggest significant temperature exaltations in the substrate in dependence of the component geometry, whereby the coating quality decreases due to residual stresses during the cooling process.Thus, for the computer aided planning of thermal spray processes, the ability to predict the temperature development for a given workpiece is of high value.
This article presents a simulation framework for the computation of the temperature development for a given workpiece during the thermal spraying process.Therein the "outer" part of the framework handles the time stepping process, the spray gun movement and computes the time-dependent input boundary conditions for the "inner" simulation module, which evaluates a thermodynamically consistent transient and nonlinear heat conduction formulation for a rigid heat conductor based on the governing equations of continuum thermodynamics, as introduced by Coleman and Noll [5].
For complex workpieces the setup of the input boundary conditions is a computationally expensive procedure by itself, since the (non-ambient) load distribution on the workpiece surface due to thermal radiation is not locally restricted and may affect a large fraction of the workpiece.The computation of the affected parts constitutes a computationally expensive visibility problem.For this, a novel approach is presented, which makes use of an elegant GPU-acceleration technique and is in fact an adaption of a method developed earlier in the context of mass deposition simulation [6].This two-level simulation approach enables the highly detailed thermodynamic model, which is specifically adjusted to the HVOF spraying process, to be efficiently applied even to larger workpieces such as forming tools in the automotive industry or complex geometries like turbine blades.Another contribution of this work is the "inner" finite element simulation module including a multi-threaded C++ based implementation which uses the solver library Eigen.
The following two sections describe the fundamental equations of continuum thermodynamics for a rigid heat conductor and the resulting finite element discretisation.Section 4 outlines the functionality of the robot guided spray simulation, [6], and provides important details regarding the GPU-accelerated computation of the input boundary conditions.Section 5 presents a demonstration of the software tool for the simulation of a real work piece.The material parameters of the underlying constitutive relations-the heat capacity c as well as the heat conduction coefficient λ -are represented by suitable functions which are fitted to experimental data.

Continuum Thermodynamics
The governing well-established equations of continuum thermodynamics for a rigid heat conductor are illustrated in the following section.The fundamental balance equation for a rigid heat conductor is the balance of energy, here in local form, , e q wherein e  denotes the rate of the internal energy, q denotes the heat flux vector and q denotes the external heat supply, while ρ represents the material density, [5].The balance of energy equation represents the first law of thermodynamics.By assuming the entropy flux to be proportional to the heat flux, as proposed by Coleman and Noll [5], the general entropy inequality results in the so called Clausius-Duhem inequality 0.
In the above inequality, η is the rate of the entropy and θ is the absolute temperature.The Clausius- Duhem inequality represents the second law of thermodynamics.A Legendre transform of the internal energy ( ) e η results in the free energy density ( ) e ψ θ θη = − so that η ψ θ = − ∂ ∂ and subsequently e ψ θ ψ θ = − ∂ ∂ .Following Liu [7] and references cited therein, the only remaining part of inequality ( 2) is the thermal dissipation therm 0 , also referred to as Fourier's inequality which is fulfilled if the heat flux points from a hot spot to a cold spot.To maintain generality, the specific heat is introduced temperature-dependent, i.e.
This equation is to be solved within a general thermodynamical initial boundary value problem.For this purpose, the boundary ∂ of the body  is decomposed into three disjoint parts, i.e., . On D ∂ Dirichlet boundary conditions are prescribed for the temperature θ , whereas Neumann and Robin boundary conditions are prescribed for the heat flux q on N ∂ respectively R ∂ : ( ) Fourier's law of heat conduction is applied throughout this work.For a review on non-Fourier heat conduction the reader is referred to other recent works, e.g., the works of Atefi and Talaee [8] [9].

Finite Element Discretisation
The different representations of the energy balance equation derived above are given in strong form.To calculate a solution for the desired field of the temperature by means of the finite element method in the context of inhomogeneous initial boundary value problems, the temperature-based balance of energy has to be reformulated in weak form.Therefore, we transfer Equation (3) to a residual form, ( ) where ω is a scalar-valued test function, see [10] for detailed background information.This relation also holds under integration over the volume of a body  , ( ) The divergence of the heat flux in Equation ( 6) can be reformulated by the application of Gauß's theorem and integration by parts, Here, n denotes the outward unit vector on the boundary of the body  .By interpreting the test function ω as the virtual temperature δθ , which is 0 at Dirichlet boundaries D ∂ , Equation (7) Equation ( 8) is the weak form of the initial boundary value problem of a rigid heat conductor which now has to be discretised in time and space, following the procedure outlined, e.g., in [11] or Kuhl et al. [12].For time discretisation, a differential-quotient-based Backward Euler integration scheme is applied [13], where [ ] • symbolises an arbitrary quantity of interest and t ∆ denotes the time increment .
In addition to the Backward Euler method applied in this work, recent works which address thermal problems use the Crank-Nicolson method [14] or energy-momentum consistent schemes in the context of thermoelastodynamics [15].At time 1 n t + , Equation ( 11) leads to the relation ( ) for the unknown temperature 1 n θ + .In view of the discretisation in space, it has to be taken into account that Robin boundary conditions require additional effort, since they represent one type of temperature dependent loads.For a review on the implementation and algorithmic treatment of deformation dependent loads, the reader is referred to the textbook by Bonet and Wood [16] and references cited therein.For that purpose, on the one hand the body  is approximated by a finite number of el n  volume elements Following the spirit of the isoparametric concept, the geometry of the body, in terms of position vectors X , as well as the temperature of the body θ and the virtual temperature δθ are interpolated element-wise by shape functions en en en en Hence, the gradients ∇ X X , θ ∇ X and δθ ∇ X are approximated as en en en en The insertion of the approximations given by Equations ( 14) and ( 15) into Equation ( 12) leads to the fully discretised balance of energy in terms of the unknown temperature where the discrete inertia, volume, internal and surface fluxes are expressed by ( ) In the above equations, A represents the assembly operator over all .Note, that q and N q are assumed to be temperatureindependent throughout this work.With the definitions in Equation ( 17), Equation ( 16) represents the discretised nonlinear temperature field equation.The solution is performed by an incremental Newton-Raphson scheme, see [17] for detailed background information in the context of the finite element method.For that purpose, the Jacobian of the residual with respect to the temperature θ has to be determined.Therefore, the iterative scheme can be expressed as , , wherein d I r is approximated by the linear term of a Taylor series expansion, The derivative .
The first term represents the dynamic contribution and therefore characterises the time-dependency of the problem, whereas the second term corresponds to the internal heat flux where the constitutive relation for the heat flux q still has to be defined.The third term reflects the surface element contribution by Robin boundary conditions

Implementation
Two aspects are vital for the computer aided planning of thermal spray processes: the prediction of the coating thickness distribution on the surface of a workpiece and the prediction of the temperatures reached in the workpiece during the process.For the simulation of the coating thickness, an efficient GPU-accelerated, C++ based simulation was developed in [6].In the present work, this framework is modified so as to compute a time dependent temperature distribution ∞ θ on the surface of the discretised workpiece mesh which serves as an input for the finite element framework developed in the previous section.In order to model the energy input of the unloaded flame of the spray gun into the workpiece, the surface heat flux is assumed to be applied by convection and radiation.An unloaded flame means that the problem at hand is restricted to heat transfer only without accounting for mass transport phenomena.One simple approach to mathematically cover the convective heat contribution is the introduction of a so-called film condition , wherein c h denotes the convective heat transfer coefficient which is assumed to be a constant, and where θ ∞ denotes the environmental temperature which either represents the HVOF spray gun or the room temperature.It should be noted that c h represents a rather complex fluid mechanical process, not captured by the finite element approach as this work proceeds.The heat radiation is modelled by 2 2 with , where ε is the emissivity and σ represents the Stefan-Boltzmann constant, cf.Comini et al. [18].Both, c q and r q , are Robin boundary conditions in manner of Equation ( 4), i.e.

( ) ( ) ( )
where ( ) θ k has to be positive semi-definite for any θ , see e.g. [7],and is thus represented by where I characterises the second order identity tensor and where ( ) 0 λ θ > is the heat conduction coefficient.
The assumption of ( ) θ k being proportional to I leads to an isotropic heat conduction model.The nonlinear heat conduction model specified via Equations ( 21)-( 25) is implemented in C++ as part of the already mentioned simulation framework.The remaining part of this section describes the computation of the film condition parameters i θ ∞ for all boundary triangles of the tetrahedron volume mesh representing the workpiece.Apart from the triangulated workpiece mesh, the simulation approach takes a robot guided spray gun movement path as an input, which is represented as a sequence of 3-tuples

{ }
, , x q for 1, , i n =  discrete robot path positions.Therein i x is the spray gun position in global coordinates and i q is a quaternion describing the orientation of the spray gun at time i t .For detailed background information on quaternions the reader is referred to Argyris [19], Betsch et al. [20], Altmann [21] and the references cited therein.
During the simulation a virtual spray gun is moved along the given path in discrete time steps t ∆ , where positions are interpolated linearly and where quaternion slerp interpolation is used for the orientations [22].For every gun position i p the following computations, where steps 1-3 are entirely GPU-based and implemented in OpenGL and the OpenGL Shading Language (GLSL), are performed: 1) Compute the subset loaded N of nodes of the workpiece mesh that receive a heat flux from the current spray gun position.This is defined to be the subset of nodes, the coordinates of which are inside the spray plume, approximated by a (spray) cone which expands from the gun center position ( ) t x , see Figure 2.Only those nodes visible from the spray gun center are considered; nodes on the back side of the workpiece or nodes shadowed by other parts of the workpiece do not receive a high temperature load.
This visibility computation for arbitrary meshes with potentially hundreds of thousands of nodes and faces can be efficiently implemented by exploiting the computational power of modern graphics cards, tailored by design to perform visibility computations on large triangle meshes in real-time.For the rendering of virtual scenes consisting of triangle meshes, the graphics card projects the three dimensional nodes of the mesh in a perspective correct manner onto a two-dimensional image plane.In the first visibility computation step, projected nodes outside the rectangular screen area are discarded, limiting the remaining nodes to the ones inside a pyramidal viewing frustum as depicted in Figure 2. In an additional step, the graphics hardware uses the depth-buffer algorithm [23] to remove any nodes that are shadowed by the geometry for the particular viewpoint considered.For the current spray gun position, first, the opening angle of the viewing frustum of the virtual camera is adjusted to match the desired opening angle of the spray cone.Subsequently, the workpiece mesh is rendered and the desired nodes that receive a heat flux from the gun are determined.The case that the viewing frustum is rectangular, but the spray cone should be circular is corrected in the next step.
2) For every node in the subset loaded N determine the thermal load node θ ∞ in order to set up Robin type boundary conditions.
The nodes are projected onto a plane orthogonal to the central axis of the spray cone at the stand off distance SOD d from the gun center.The thermal load for every node is then computed based on its distance no r from the cone center axis.Due to the lack of precise measurement data, the load function is modelled as a rotationally symmetric Gaussian with standard deviation s and amplitude max θ which is offset by in r in order to cover a larger circular area of magnitude max . The threshold max r is used to limit the function to a circular radius thereby transforming the rectangular viewing frustum into the desired circular spray plume approximation.For the simulations presented in this article an HVOF spraying process is assumed and the values SOD 170 mm d = , max gun θ θ ∞ = , in 10 mm r = , max 66 mm r = and 2 mm s = are used.The shape of the load function is shown in Figure 3.
3) The remaining nodes are set to receive an ambient thermal load amb θ ∞ in order to set up Robin type boundary conditions.
4) The thermodynamic model is evaluated.For each Newton iteration, the C++ implementation of the thermodynamic model assembles the global iterative residual vector I r and the global tangent matrix IJ K according to Equations ( 19) and ( 20), and then employs a conjugate gradient iterative solver with Jacobi (diagonal) preconditioning (pCG) to solve the system.

5) The simulation time index is increased by t
∆ .
This scheme is iterated until the end of the robot guided path is reached.A key aspect for the efficient computation of the boundary conditions is the similarity of the geometric relationships of the spray plume expanding from the gun nozzle towards the workpiece on the one hand, and the geometry-projection process towards a virtual camera center in the visibility computation on the other.Due to this similarity, the first three steps can easily be performed by exploiting GPU acceleration techniques and are implemented in OpenGL and the OpenGL Shading Language (GLSL) within our framework.

Examples
The implemented framework is now applied to carry out simulations of a realistic deep drawing tool in order to demonstrate the capabilities of the novel software.The starting temperature 0 t θ of the material is considered to coincide with the room temperature amb 293 K θ ∞ = whereas the temperature of the spray gun is set to gun 3073 K θ ∞ = .The steel material is St 35.8 (1.0305), the physical properties of which have been investigated by Richter [24] [25] up to a temperature of 873 K.The parameters for Equation (23), the polynomial of the heat capacity, are taken from [24] [25] and summarised in Table 1.Note, that the polynomial is parameterised in ˚C not in K.For the sake of completeness, the polynomial is depicted with a K-scale in Figure 4(a).It can be seen that the parameter set is not valid for temperatures beyond the interval , the implemented framework remains thermally stable even when the physical valid temperature range is excelled.This turns out to be helpful for testing the algorithmic framework.Furthermore, the parameters for the heat conduction coefficient ( ) λ θ , given by Equation (25), are also summarised in Table 1.Because the parameters are the result of fitting a polynomial of Richter, [25], both functions are depicted in Figure 4(b).The solid line represents the fitted function and, in contrast, the dashed line represents the polynomial taken from Richter, [25].Obviously, the fit is very good within the physical valid temperature range.Beyond this range, the fitted curve asymptotically tends to a positive limit while the Richter polynomial has a root at 1451.58 K θ = . As an advantage, the fitted function guarantees ( ) 0 λ θ > for any θ, thereby avoiding numerical problems while testing the algorithmic framework.To complete the set of parameters, values for the convective heat transfer coefficients air c h and gun c h as well as for the emissivity ε and the Stefan-Boltzmann constant σ are summarised in Table 2.We assume the surface surrounded by air to be exposed to natural convection, as the surface loaded by the HVOF spray gun is assumed to be charged with forced convection.
The integrated coating and thermodynamic simulation is employed to predict the workpiece temperature for the manufacturing of thermal sprayed sheet metal forming tools.In this regard, the heat transport into a deep drawing tool of a moving and unloaded HVOF thermal spray gun is simulated using a complex robot tool path at three different movement speeds.The fastest simulation run completes the tool path within 40 s, the medium speed simulation in 80 s and the slowest run within 160 s.The tool path of the spray gun movement is depicted ; see [25] for further details.
( )      in Figure 5.The colors of the robot path represent the temperature load resulting from the respective gun position.For the robot path, the blue regions on the right and left parts have not been processed yet, since the gun started at the lower end of the path proceeding upwards along the circular opening of the workpiece onto the plain surface on top.For the path already processed, the coloring represents the surface load temperatures induced by the gun, where grey areas represent zero thermal load.In other words, these grey parts represent positions of the spray gun where the spray cone does not intersect with the workpiece.The surface colors represent the surface temperature at different time steps for the medium gun speed simulation run.
The simulation keeps track of the volume-averaged temperature and the peak workpiece temperature both with respect to time to identify critical sections of the path, which may cause an overheating of the workpiece.Figure 6 shows the simulated temperatures averaged over the entire workpiece volume and the peak temperature occurring in the workpiece-both plotted versus time-for the three gun speed simulation runs.The spraying process can be divided into four sections-compare Figures 5(b)-(d): First, the spray gun aims at the front radius working its way up to the plain surface, which is then coated in circular arc movements from left to right and back.During this process, peak temperatures are reached when the gun approaches the border of the workpiece and turns around for the next arc in opposite direction.The fluctuations of the peak temperature curves are due to the gun leaving and entering the workpiece at the borders.The third and fourth distinguishable sections are the right and left border areas of the workpiece including the screw holes for mounting the part.At the beginning of these sections, the peak and average temperatures increase considerably, because the large screw holes provide a significantly larger surface area that is heated by the spray gun.Afterwards the temperatures drop again when the gun gradually leaves the workpiece area.It can be seen that, for the considered robot tool path, the slowest simulated gun movement speed leads to the highest mean temperature of the workpiece as well as to the highest surface peak temperatures.The heat output to the environment through the surface is lower than the heat input by the spray gun.

Summary
This work presents a novel integrated simulation approach for the computation of the workpiece temperature development during thermal spraying which is suited for large and arbitrarily complex shaped components.As its core, a thermodynamically consistent nonlinear heat transfer model based on the fundamental equations of continuum thermodynamics is employed.The equations obtained are transformed into their weak forms that represent the underlying equations of the presented non-linear finite-element-framework, implemented in C++ and the Eigen library.The material parameters are adjusted to experimental data by fitting suitable functions, and the parameters of the Robin boundary conditions are chosen in order to model realistic heat transfer.The presented simulations of a complex shaped deep drawing tool demonstrate a possible field of application for the novel software tool.Future work will include the consideration of the mass flow of the HVOF gun in order to model the coating process of the substrate with hot particles.Furthermore, the developed framework shall be used to extend the optimisation software described in [26] in order to be able to optimise temperature quantities.Another aspect of future research is the extension towards a thermo-mechanically coupled framework, [27], which will enable the analysis of residual stresses and damage due to thermal spraying as discussed in, e.g., [28] [29].

Figure 1 .
Figure 1.Photograph of the HVOF spraying process, kindly provided by LWT, TU Dortmund.
balance (1) is now restated in the strong form of the temperature field equation valid for a general rigid heat conductor, ( )

Figure 2 .
Figure 2. Visibility computation to determine the nodes of the workpiece to be coated and to receive a heat flux from the gun (white nodes).Nodes outside the spray cone or nodes that are shadowed from the gun by parts of the workpiece are not coated.

Figure 3 .
Figure 3. Shape of the radially symmetric load function node θ ∞ used in the emissivity ε and the Stefan-Boltzmann constant σ . air

Figure 4 .
Figure 4. (a) Richter polynomial of the heat conduction ( ) c θ ; (b) Fitted curve and the Richter polynomial of the heat

Figure 5 .
Figure 5. Simulated deep drawing workpiece and the employed robot movement path.The surface colors represent the surface temperatures and the path colors the temperature load after (a) 6 s; (b) 20 s; (c) 30 s; and (d) 40 s of the medium gun speed simulation run.Blue path segments have not been processed yet.

Figure 6 .
Figure 6.Simulated development of peak and mean workpiece temperature of the workpiece shown in Figure 5.The temperatures are depicted for the three gun speed simulation runs.
represents the virtual temperature problem and can be written in terms of dynamic, volume, internal and surface terms + , whereas subscript n denotes a quantity at the previous time step n t .Application of Equation (10) to the temperature θ results in 1

Table 2 .
Values for the convective heat transfer coefficients air