Natural Neighbour Strategies for the Simulation of Laser Surface Coating Processes

In this paper two different formulations for the numerical simulation of laser surface coating processes are presented and analysed. Both are based on the use of natural neighbour interpolation, but one employs a Galerkin approach and takes temperature as the primary variable while the second one is based on the use of finite differences and enthalpy as primary variable. The main practical difference is thus the description of the interphase of the melted zone during the process. While in the first method the interphase is described by a set of nodes that evolve in time, in the second one it is located somewhere between a string of nodes. Both formulations are described and compared, showing their potential benefits for the simulation of the before mentioned process.


Introduction
In this work we analyse the behaviour of methods based on natural neighbour interpolation applied to the simulation of laser surface coating problems arising, for instance, in superconductor texturing processes or ceramic tiles processing, see The NEM (also known as Natural Neighbour Galerkin Method) belongs to the wide family of meshless methods and is based on the use of any natural neighbourbased interpolation scheme in the framework of a Galerkin scheme. NEM allows the use of a mobile set of nodes that describe the boundary between phases and moves over a cloud formed by fixed points or formed by mobile points which movement is independent of the phase's interaction.
Using natural neighbour interpolants it is also possible to construct finite differences schemes over irregular clouds of points (See [CUE 03] and references therein). A method of finite differences based on the enthalpy formulation of the Stefan problem is also presented and compared with the Galerkin formulation.
Methods based on FE usually conform element boundaries to the interfaces as it evolves, but these methods are limited by the severe mesh distortion. Thus, remeshing is needed, which is a delicate task in the 3D case. Meshless methods, however, discretize a continuum body by a finite number of nodes and the field on study is interpolated without the aid of an explicit mesh. In this paper, a technique to detect the interface location using Natural Neighbour-based Galerkin methods is presented and it is applied to several problems to show its behaviour. The results obtained are not affected by the relative position of the nodes; this is an important characteristic of meshless methods. Results are also compared with those obtained by employing natural neighbour finite difference schemes.
The final aim of this work would be, therefore, to analyse the influence of the support temperature on the solidification profile and its correlation with the microstructure properties of the samples textured using the Laser Zone Melting (LZM) processes by means of numerical simulation of this behaviour.

Natural Neighbour interpolation
Natural neighbour Galerkin or finite differences methods can be considered members of the vast family of meshless methods. They rely on any of the various natural neighbour-based interpolation schemes to construct the approximation. So, it is necessary to introduce these kind of interpolation prior to describe the characteristics of the proposed methods.
There exist various types of natural neighbour-based interpolations, but they all are based on the construction of the Delaunay triangulation [DEL 34] of the cloud of points, D, used to discretise the domain. The Delaunay triangulation of a set of points N = {n 1 , n 2 , . . . , n N } is the unique triangulation of the set that verifies the so-called empty circumcircle criterion. This means that no point of the set lies in the interior of a circle that passes through the three vertices of each triangle, see Fig. 2.
The dual structure of the Delaunay triangulation is the Voronoi diagram of the cloud [VOR 08]. For a given node n I , its associated Voronoi cell is composed by all of the points which are closer to the node n I than to any other node. Formally, where A I is the Voronoi cell and d(·, ·) represents the Euclidean distance. In the problems considered in this paper, n = 3. Two nodes whose Voronoi cells share one edge are called natural neighbours and hence the name of these interpolation schemes.
The first, and most obvious, interpolation scheme based on natural neighbours is the so-called nearest neighbour or Thiessen interpolation [THI 11]. If we give the nodal value to the whole associated Voronoi cell, we obtain a C −1 interpolation scheme. This interpolation scheme is not suitable for solving second-order partial differential equations, but has been employed in [GON 04b] to construct mixed velocitypressure approximations for the simulation of incompressible media.
The most extended form of natural neighbour-based interpolation is due to Sibson [SIB 81]. For the definition of Sibson interpolation it is necessary to previously introduce the concept of second order Voronoi cell. It is defined as the locus of the points that have the node n I as the closest node and the node n J as the second closest node: (2) If a new point is added to a given cloud of points, the Voronoi cells will be modified by its presence. Sibson [SIB 80] defined the natural neighbour coordinates of a point x with respect to one of its neighbours I as the ratio of the cell A I that is transferred to A x to the total area of A x . In other words, being κ(x) and κ I (x) the Lebesgue measures of A x and A xI respectively, the natural neighbour coordinates of x with respect to the node I is defined as The resultant shape function depends obviously on the relative position of the nodes. An example for a node surrounded by other six is depicted in Fig. 5.
The resultant shape function has some remarkable properties (see [SUK 98] or [CUE 03] for more in-deep explanations and rigorous proofs of this behaviour). Firstly, it is smooth (C 1 ) everywhere except from the nodes, as can be seen from Recently, Hiyoshi and co-workers [HIY 99] have generalised the form of natural neighbour interpolants. One different type of interpolation has attracted the interest of researchers, since it is slightly faster to compute, although gives less smooth interpolations. It has received the name of Laplace interpolant. The Laplace interpolant is defined by using geometrical entities of one dimension less than the original space under consideration. If we define the cell intersection t IJ = {x ∈ T I T J , J = I} (note that t IJ may be an empty set) we can define the value .
(4) Thus, the point x shape function value with respect to node 4 in Fig. 4 is defined as where s J represent the length of the Voronoi segment associated to node J and n represents the number of natural neighbours of the point under consideration, x.

Figure 4. Definition of non-Sibsonian coordinates.
Derivatives of the Laplace shape function are not defined along the edges of the Delaunay triangles that lie inside its support (see [SUK 01]).

Figure 5. (a) Natural Element (Sibson) shape function and (b) non-Sibsonian or Laplace shape function (photo courtesy N. Sukumar).
In the context of two-and three-dimensional approximations, the unknown variable is approximated in the form: where T I is the vector of nodal values and n the number of natural neighbours of each point x. This leads to a C 0 interpolation scheme.

Governing equations
Mathematically, laser surface coating processes can be viewed as a particular instance of Stefan problems [STE 91]. Consider a domain formed by two phases, liquid and solid, respectively, Ω 1 and Ω 2 ( Figure 6).
In the boundary of the domain, Γ = ∂Ω, either the temperature (T = T on Γ 1 ), or heat flux (q = q on Γ 2 ) are imposed. The problem is defined in a time interval [0, t max ], being the initial temperature T (x, t 0 ) = T 0 prescribed in the whole domain. At time t = t 0 the temperature on Γ 1 is set to a value T 1 > T m , being T m the material's melting temperature. In laser surface coating processes, the region Γ 1 can be seen as the area under the laser beam, that can move or change during the process. Figure 6. A domain formed by a liquid phase and a solid one.
A melting front Γ I is generated that moves at a velocity v, dividing the domain Ω in the two before mentioned regions.
Each phase is assumed to be homogeneous and isotropic. Equilibrium equations in each phase can be stated as: where c 1 and c 2 represent, respectively, the heat capacities of solid and liquid phases, while k 1 and k 2 represent the thermal conductivity of each phase. Initial and boundary conditions of the problem are: Evolution of the interphase Γ I is governed by the so-called Stefan's condition: where v represents the velocity of the interphase, L is the latent heat, n I (x) is the outward unit vector of the interphase, directed towards the liquid phase and, finally, |[q]| represents the discontinuity of the flux across the interphase: Finally, one last restriction affects the melting temperature, that must remain constant along the interphase:

Natural neighbour Galerkin discretization
The weak form of the problem (7)-(8) will read: being H 1 and H 1 0 the usual Sobolev spaces of order one. Substituting in this expression the unknown and admissible variations fields by a suitable approximation obtained through natural neighbour interpolation, we arrive to: where T is a vector containing the nodal temperatures. Equation (12) is discretised through the generalised midpoint rule: thus obtaining: where and where N is a vector containing the shape functions and B la matrix with shape function derivatives.
Integration of the before presented equations is achieved by employing three quadrature points within each Delaunay triangle. In [GON 04a] it has been demonstrated that this technique leads to errors in the quadrature that can be avoided, for example, by employing stabilised conforming nodal quadrature schemes.

Interphase tracking
Interphase tracking in FE-like discretisation of Stefan problems deserves some comments. As mentioned before, two different sets of nodes are considered, namely a fixed one 1 to discretise the domain and a moving one that discretises the interphase. This last set of nodes will create a planar straight-line graph in two-dimensional problems, while in three-dimensional ones it will discretise a surface.
The problem is thus to determine the geometry of the evolving phases as the moving set of nodes is governed by the Stefan condition, Eq. (9). The planar straight line graph is equipped with a normal vector at each segment, pointing to the liquid phase, in order to determine if a node in the first set belongs to the melted or solid phases. Once each node is located within its sub-domain (see Fig. 7) the shape of the sub-domain is extracted by employing a shape constructor. A shape constructor is a geometrical entity that gives shape to a set of nodes. In our case we employ α-shapes [EDE 94]. Briefly, an α-shape is a polytope (a subset of the Delaunay triangulation of the cloud) whose geometry depends on a value α that represents the level of detail up to which we wish to represent the domain. Of course, this level of detail can not be finer than the nodal spacing parameter, h. Thus, every triangle whose cimcurcircle possesses a radius bigger than α is eliminated from the triangulation. The employ of α-shapes has another advantage, since it allows to obtain linear precision along the boundary of the domain, whenever it is convex or not [CUE 00]. Thus, a fully conforming method is obtained.

Modelling the melting process with an enthalpy formulation involving an interphase capturing
In this section we propose an alternative approach for the melt front treatment. In order to simulate the change of phase process, an explicit finite volumes scheme will be applied, which, as proved later, allows to proceed with a fixed mesh of the whole domain. In order to avoid the difficulties related to the change of phase accounting, an enthalpy formulation will be used. To describe easily the proposed strategy, which will considered in this work, first we consider the one-dimensional heat transfer problem defined in the interval containing N nodes uniformly distributed. A representative volume is associated to each node, where both the temperature and the enthalpy are assumed constant. The temperature is prescribed at both boundaries as well as the initial temperature. The enthalpy is defined by: where c 2 and c 1 are used to identify the solid and liquid phases respectively. ΔH m represents the fusion enthalpy and T ref is a reference temperature.
The energy balance, taking into account eventual phase changes, results: If we consider N nodes uniformly distributed (h being the distance between two consecutive nodes), and assuming that the initial enthalpy can be computed using equation (19), the discretisation of Eq. (20) results Being the nodes uniformly distributed, the heat flux derivatives at the cell boundaries can be approximated by: In the general 2D or 3D unstructured case the Voronoi cells can be considered as finite volumes. The temperature and enthalpy are assumed constant within each Voronoi cell. Eqs. (19), (20) and (24) are still valid. However, the discretized form related to Eq. (20) results after integration in each Voronoi cell Ω i : where n represents the unit outwards vector, |Ω i | the area or volume of cell Ω i and ∂Ω i its boundary. The integral term in Eq. (25) can be approximated as: where N i is the number of neighbor cells related to cell Ω i , L ij the length of the common edge between Ω i and Ω j , ∇T n | ij the temperature gradient evaluated at the middle point of that edge, and n ij the unit outwards vector related to that edge. In the scheme here presented the temperature gradient has been computed by using the natural neighbor shape functions derivatives.

Validation of the proposed Galerkin scheme
One of the few Stefan problems with analytical solution is employed here in order to validate the accuracy of the proposed technique. It consists of a rectangular plate subjected to a sudden fall of temperature along one of its sides, thus leading to frozen of the material. This problem was also analysed by employing a fully implicit method in [YVO 05].
Initial temperature T 0 is assumed constant and at time t = 0, one of the plate's sides is subjected to a temperature T 1 , under the temperature of solidification. For a semi-infinite plate, the analytical solution to this example can be found in [YVO 05].
The cloud of points for the discretisation of the problem consisted of 66 background nodes, discretising the geometry of the domain and 6 moving nodes, describing the interphase location, see

b) Solid phase (blue) and liquid (red)
Evolution of the interphase is near the vertical expected position, with some alterations near the boundary of the domain, where the number of neighbours is lower, and consequently the quality of the approximation decreases.
The position of the interphase (as a mean value of the position of its discretising nodes) compared with the analytical results is show in Fig. 12. Similar results can be obtained for three-dimensional settings, as shown in Fig. 13.

Simulation of a plate under a moving heat source
In this example we study the case of a moving heat source, representing the laser beam, that runs on the top of a plate, as represented in Fig. 14. The plate is isolated at the bottom. The heat source moves at a speed such that it remains at a node location at each time step. Material characteristics are the same as in Table 1.
This problem thus represents a simplified process of surface treatment by laser. Unfortunately, there is no analytical solution to this problem, but qualitative agreement with melted zone patterns obtained experimentally can be observed. Different   snapshots of the process are shown in Fig. 15. See Fig. 16 for a transversal cut of a laser textured superconductor. The shape of the melted zone can be observed in the center of the image.
In three dimensional cases this same qualitative agreement can be observed (see Fig. 17).

Voronoi finite volumes simulation
In this section we consider the application of a laser beam on a metallic surface. When the temperature reaches the vaporization temperature the volume is removed and the heating process can continue in the new domain geometry. Of course, at present we have not considered the effects of the vaporized metal because they require the plasma modelling.
In Fig. 19 we can appreciate the temperature evolution at the three positions indicated in Fig. 18. We can notice that the temperature remains constant whereas the considered volumes are changing of phase. Obviously, in this kind of approach it is not possible to define accurately the position of the interphase.

Conclusion
In this paper two different schemes are proposed and compared for the simulation of laser surface coating processes. Both are based in the use of natural neighbour interpolation. The first one is based in a Galerkin implementation, that allows for an accurate description of the interphase. By employing a moving set of nodes that discretize the interphase it is possible to accurately describe its position, while the resulting "elements" are not affected by the distortion, as in FE techniques.
The second technique is based on a natural neighbour finite difference scheme and a formulation of the problem based in the enthalpy, instead of temperature. This alternative formulation proves to be also very accurate, but the exact position of the interphase is lost.   Both techniques render accurate results and represent alternative ways of solving the same problem. However, a validation of the proposed techniques with some experiments (say, for instance, for the problem of superconductor texturing processes) is lacking. This validation constitutes nowadays the main work of the authors and will be published hopefully elsewhere.   Fig. 18.