# Natural Neighbour Strategies for the Simulation of Laser Surface Coating Processes

David Gonzáles, David Bel, Elías Cueto, Francisco Chinesta, Manuel Doblaré

► To cite this version:

David Gonzáles, David Bel, Elías Cueto, Francisco Chinesta, Manuel Doblaré. Natural Neighbour Strategies for the Simulation of Laser Surface Coating Processes. International Journal of Forming Processes, 2007, 10 (1), pp.89 - 108. 10.3166/ijfp.10.89-108. hal-01007101

HAL Id: hal-01007101

https://hal.science/hal-01007101

Submitted on 29 Jan 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
Natural Neighbour strategies for the simulation of
laser surface coating processes

D. González<fnref n="1" /> — D. Bel<fnref n="1" /> — E. Cueto<fnref n="1" /> — F. Chinesta<fnref n="2" /> — M. Doblaré<fnref n="1" />

ABSTRACT. 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.

KEYWORDS: Meshless methods, Partition of Unity methods, Natural Neighbour interpolation, laser surface coating.

<footnote n="1">¹ Aragon Institute of Engineering Research. University of Zaragoza.
Edificio Betancourt. María de Luna, 5, E-50018. Zaragoza, Spain
URL: http://mmcyte7.cps.unizar.es
e-mail: ecueto@unizar.es</footnote>

<footnote n="2">² LMSP (Laboratoire de Mécanique des Systèmes et des Procédés)
UMR 8106 CNRS-ENSAM-ESEM
151 boulevard de l'Hôpital, F-75013 Paris, France
URL: http://www.paris.ensam.fr/LMSP
e-mail: francisco.chinesta@paris.ensam.fr</footnote>
# 1. 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 Fig. 1. These kind of problems can be viewed as a particular instance of Stefan problems.

<figure><img src="image_2.png" /><figcaption>**Figure 1.** *Typical apparatus to texture a superconductor.*</figcaption></figure>

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 neighbour-based 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
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Delaunay triangulation and Voronoi diagram of a cloud of points. On the right, an example of a degenerate distribution of nodes, with the two possible triangulations depicted. In this last case, four points lie in the same circumcircle and thus no single triangulation exists.</figcaption></figure>

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.

## 2. 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, $\mathcal{D}$, used to discretise the domain. The Delaunay triangulation of a set of points $N = \{n_1, n_2, \dots, 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,

$$A_I = \{ \mathbf{x} \in \mathbb{R}^n : d(\mathbf{x}, \mathbf{x}_I) < d(\mathbf{x}, \mathbf{x}_J) \ \forall J \neq I \}, \quad (1)$$
where $\mathcal{A}_I$ is the Voronoi cell and $d(\cdot, \cdot)$ 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 velocity-pressure 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:

$$ \mathcal{A}_{IJ} = \{ \mathbf{x} \in \mathbb{R}^n : d(\mathbf{x}, \mathbf{x}_I) < d(\mathbf{x}, \mathbf{x}_J) < d(\mathbf{x}, \mathbf{x}_K) \forall K \neq J; \forall K \neq I \}. \quad (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 $\mathbf{x}$ with respect to one of its neighbours $I$ as the ratio of the cell $\mathcal{A}_I$ that is transferred to $\mathcal{A}_x$ to the total area of $\mathcal{A}_x$. In other words, being $\kappa(\mathbf{x})$ and $\kappa_I(\mathbf{x})$ the Lebesgue measures of $\mathcal{A}_x$ and $\mathcal{A}_{xI}$ respectively, the natural neighbour coordinates of $\mathbf{x}$ with respect to the node $I$ is defined as

$$ \phi_I^{sib}(\mathbf{x}) = \frac{\kappa_I(\mathbf{x})}{\kappa(\mathbf{x})}. \quad (3) $$

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 Fig. 5. Natural neighbour (Sibson) shape functions possess linear completeness [SUK 98] and form a partition of unity. Therefore, it is possible to enrich natural neighbour interpolants in order to increase the order of the polynomial reproduced by the interpolation, as proposed in [BAB 96].

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} = \{\mathbf{x} \in T_I \cap T_J, J \neq I\}$ (note that $t_{IJ}$ may be an empty set) we can define the value

$$ \alpha_J(\mathbf{x}) = \frac{|t_{IJ}|}{d(\mathbf{x}, \mathbf{x}_J)}. \quad (4) $$
<figure><img src="image_4.png" /><figcaption>Figure 3. Definition of the Natural Neighbour coordinates of a point x.</figcaption></figure>

Thus, the point x shape function value with respect to node 4 in Fig. 4 is defined as

$$ \phi_4^{ns}(\mathbf{x}) = \frac{\alpha_4(\mathbf{x})}{\sum_{J=1}^{n} \alpha_J(\mathbf{x})} = \frac{s_4(\mathbf{x})/h_4(\mathbf{x})}{\sum_{J=1}^{n} [s_J(\mathbf{x})/h_J(\mathbf{x})]}, \quad (5) $$

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><img src="image_5.png" /><figcaption>Figure 4. Definition of non-Sibsonian coordinates.</figcaption></figure>

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><img src="image_6.png" /><figcaption>**Figure 5.** (a) Natural Element (Sibson) shape function and (b) non-Sibsonian or Laplace shape function (photo courtesy N. Sukumar).</figcaption></figure>

In the context of two- and three-dimensional approximations, the unknown variable is approximated in the form:

$$T^h(\mathbf{x}) = \sum_{I=1}^{n} \phi_I(\mathbf{x}) T_I \qquad (6)$$

where $T_I$ is the vector of nodal values and $n$ the number of natural neighbours of each point $\mathbf{x}$. This leads to a $C^0$ interpolation scheme.

## 3. Modelling the melting process as a Stefan problem involving an interphase tracking

### 3.1. 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, $\Omega_1$ and $\Omega_2$ (Figure 6).

In the boundary of the domain, $\Gamma = \partial\Omega$, either the temperature ($T = \bar{T}$ on $\Gamma_1$), or heat flux ($q = \bar{q}$ on $\Gamma_2$) are imposed. The problem is defined in a time interval $[0, t_{max}]$, being the initial temperature $T(\mathbf{x}, t_0) = T_0$ prescribed in the whole domain. At time $t = t_0$ the temperature on $\Gamma_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 $\Gamma_1$ can be seen as the area under the laser beam, that can move or change during the process.
<figure><img src="image_7.png" /><figcaption>**Figure 6.** A domain formed by a liquid phase and a solid one.</figcaption></figure>

A melting front $\Gamma_I$ is generated that moves at a velocity $\mathbf{v}$, dividing the domain $\Omega$ in the two before mentioned regions.

Each phase is assumed to be homogeneous and isotropic. Equilibrium equations in each phase can be stated as:

$$ \begin{cases} c_1 \frac{\partial T(\mathbf{x},t)}{\partial t} = \nabla \cdot (k_1 \nabla T) & \text{in } \Omega_1 \\ c_2 \frac{\partial T(\mathbf{x},t)}{\partial t} = \nabla \cdot (k_2 \nabla T) & \text{in } \Omega_2 \end{cases} \qquad (7) $$

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:

$$ \left\{ \begin{array}{ll} T(\mathbf{x}, t = 0) = T_0, & \forall \mathbf{x} \in \Omega \\ T(\mathbf{x}, t) = \bar{T}(\mathbf{x}, t), & \forall \mathbf{x} \in \Gamma_1 \\ -k \nabla T(\mathbf{x}, t) \cdot \mathbf{n} = \bar{q}(\mathbf{x}, t), & \forall \mathbf{x} \in \Gamma_2 \end{array} \right. \qquad (8) $$

Evolution of the interphase $\Gamma_I$ is governed by the so-called *Stefan's condition*:

$$ \mathbf{v}(\mathbf{x} \in \Gamma_I(t)) = \frac{|[q]|}{L} \cdot \mathbf{n}_I \qquad (9) $$

where $\mathbf{v}$ represents the velocity of the interphase, $L$ is the latent heat, $\mathbf{n}_I(\mathbf{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:

$$ |[q]| = (k_1 \nabla T|_{\Gamma_I^+} - k_2 \nabla T|_{\Gamma_I^-}) \cdot \mathbf{n}_I \qquad (10) $$

Finally, one last restriction affects the melting temperature, that must remain constant along the interphase:

$$ T(\mathbf{x}, t) = T_m \quad \forall \mathbf{x} \in \Gamma_I(t) \qquad (11) $$
## 3.2. Natural neighbour Galerkin discretization

The weak form of the problem (7)-(8) will read:

Find $T \in H^1(\Omega)$ satisfying $T = \bar{T}$ on $\Gamma_I$ such that:

$$ \int_{\Omega} c \frac{\partial T}{\partial t} \delta T d\Omega = - \int_{\Omega} k \nabla T \cdot \delta T d\Omega + \int_{\Gamma_I(t)} |[q]| \delta T d\Gamma, \quad \forall \delta T \in H_0^1(\Omega) \quad (12) $$

being $H^1$ and $H_0^1$ 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:

$$ \mathbf{C}\dot{\mathbf{T}} + \mathbf{K}\mathbf{T} = \mathbf{F} \quad (13) $$

where $\mathbf{T}$ is a vector containing the nodal temperatures. Equation (12) is discretised through the generalised midpoint rule:

$$ \frac{\partial T^{n+1}}{\partial T} = \frac{T^{n+1} - T^n - (1 - \alpha)\Delta t \frac{\partial T^n}{\partial t}}{\alpha \Delta t} \quad (14) $$

thus obtaining:

$$ (\mathbf{C}^{n+1} + \alpha \Delta t \mathbf{K}^{n+1}) \mathbf{T}^{n+1} = \mathbf{F}^{n+1} \quad (15) $$

where

$$ \mathbf{C}^{n+1} = \int_{\Omega^{n+1}} \mathbf{N}^t c \mathbf{N} d\Omega \quad (16) $$

$$ \mathbf{K}^{n+1} = \int_{\Omega^{n+1}} \mathbf{B}^t k \mathbf{B} d\Omega \quad (17) $$

and

$$ \mathbf{F}^{n+1} = \mathbf{C}^{n+1} \mathbf{T}^n + (1 - \alpha) \Delta t \int_{\Omega^{n+1}} \mathbf{N}^t \frac{\partial \mathbf{T}^n}{\partial t} d\Omega + \alpha \Delta t \int_{\Gamma_I^{n+1}} \mathbf{N}^t |[q]|^{n+1} d\Gamma \quad (18) $$

where $\mathbf{N}$ is a vector containing the shape functions and $\mathbf{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.

## 3.3. 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<fnref n="1" /> to discretise the domain and a moving one that discretises the interphase.

<footnote n="1">1. If a deformation process is linked to the problem, this set of nodes can also move, but this movement is not related with the thermal problem.</footnote>
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 $\alpha$-shapes [EDE 94]. Briefly, an $\alpha$-shape is a polytope (a subset of the Delaunay triangulation of the cloud) whose geometry depends on a value $\alpha$ 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 circumcircle possesses a radius bigger than $\alpha$ is eliminated from the triangulation. The employ of $\alpha$-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.

<figure><img src="image_8.png" /><figcaption>**Figure 7.** *Algorithm to extract the shape of each phase.*</figcaption></figure>

## 4. 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:

$$H(T) = \begin{cases} \int_{T_{\text{ref}}}^{T} c_2 dT & \text{if } T < T_m \\ \int_{T_{\text{ref}}}^{T_m} c_2 dT + \Delta H_m + \int_{T_m}^{T} c_1 dT & \text{if } T > T_m \end{cases} \quad (19)$$

where $c_2$ and $c_1$ are used to identify the solid and liquid phases respectively. $\Delta H_m$ represents the fusion enthalpy and $T_{\text{ref}}$ is a reference temperature.

The energy balance, taking into account eventual phase changes, results:

$$\frac{\partial H}{\partial t} - \nabla \cdot (k \nabla T) = 0 \quad (20)$$

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

$$\frac{H_i^{n+1} - H_i^n}{\Delta t} - \frac{k_{i+1/2} \frac{\partial T}{\partial x}|_{i+1/2}^n - k_{i-1/2} \frac{\partial T}{\partial x}|_{i-1/2}^n}{h} = 0 \quad \forall i \in [2, \dots, N-1] \quad (21)$$

Being the nodes uniformly distributed, the heat flux derivatives at the cell boundaries can be approximated by:

$$\left. \frac{\partial T}{\partial x} \right|_{i-1/2}^n = \frac{T_i^n - T_{i-1}^n}{h} \quad (22)$$

and

$$\left. \frac{\partial T}{\partial x} \right|_{i+1/2}^n = \frac{T_{i+1}^n - T_i^n}{h} \quad (23)$$

Now, from the enthalpy just computed using Eqs. (21)-(23), the nodal temperatures $T_i^{n+1}$ can be updated using Eq. (19), i.e.:

$$\left\{ \begin{array}{ll} H_i^{n+1} = \int_{T_{\text{ref}}}^{T_m} c_2 dT + \Delta H_m + \int_{T_m}^{T_i^{n+1}} c_1 dT & \text{if } H_i^{n+1} > \int_{T_{\text{ref}}}^{T_m} c_2 dT + \Delta H_m \\ T_i^{n+1} = T_m & \text{if } H_i^{n+1} > \int_{T_{\text{ref}}}^{T_m} c_2 dT \text{ and} \\ & H_i^{n+1} < \int_{T_{\text{ref}}}^{T_m} c_2 dT + \Delta H_m \\ H_i^{n+1} = \int_{T_{\text{ref}}}^{T_i^{n+1}} c_2 dT & \text{if } H_i^{n+1} < \int_{T_{\text{ref}}}^{T_m} c_2 dT \end{array} \right. \quad (24)$$

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 $\Omega_i$:

$$ \frac{H_i^{n+1} - H_i^n}{\Delta t} |\Omega_i| - \int_{\partial \Omega_i} k \nabla T^n \cdot n \, dS = 0 \quad (25) $$

where $n$ represents the unit outwards vector, $|\Omega_i|$ the area or volume of cell $\Omega_i$ and $\partial \Omega_i$ its boundary. The integral term in Eq. (25) can be approximated as:

$$ \int_{\partial \Omega_i} k \nabla T \cdot n \, dS = \sum_{j=1}^{j=N_i} L_{ij} \nabla T^n|_{ij} \cdot n_{ij} \quad (26) $$

where $N_i$ is the number of neighbor cells related to cell $\Omega_i$, $L_{ij}$ the length of the common edge between $\Omega_i$ and $\Omega_j$, $\nabla 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.

# 5. Numerical results

## 5.1. 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].

<table><caption>**Table 1.** Properties of liquid and solid phases.</caption><thead><tr><th>PROPERTY</th><th>SOLID</th><th>LIQUID</th></tr></thead><tbody><tr><td>Heat capacity (cal · °C<sup>-1</sup>cm<sup>-3</sup>)</td><td>0.49</td><td>0.62</td></tr><tr><td>Conductivity (cal · cm<sup>-1</sup>s<sup>-1</sup> °C<sup>-1</sup>)</td><td>9.6 × 10<sup>-3</sup></td><td>6.9 × 10<sup>-3</sup></td></tr><tr><td>Melting temperature (°C)</td><td></td><td>0.0</td></tr><tr><td>Latent heat (cal · cm<sup>-3</sup>)</td><td></td><td>19.2</td></tr></tbody></table>

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 Fig. 8. The time step was chosen to be $\Delta t = 0.002s$ in an explicit scheme.

Solutions for different time steps are shown in Figs. 9 to 11.
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Initial cloud of points and associated triangulation. The second column of nodes on the left represents the moving set.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 9.** (a) Temperature distribution at $t = 0$. (b) Solid phase (blue) and liquid (red)</figcaption></figure>

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 shown in Fig. 12. Similar results can be obtained for three-dimensional settings, as shown in Fig. 13.
<figure><img src="image_11.png" /><figcaption>**Figure 10.** (a) *Temperature distribution after 200 time steps.* (b) *Solid phase (blue) and liquid (red)*</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 11.** (a) *Temperature distribution after 300 time steps.* (b) *Solid phase (blue) and liquid (red)*</figcaption></figure>

**5.2. 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
<figure><img src="image_13.png" /><figcaption>**Figure 12.** *Theoretical vs. numerical solution for the solidification problem described in § 5.1.*</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 13.** *Temperature distribution in the three-dimensional problem at time steps 1 (a) and 60 (b).*</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Figure 14.** Schematic representation of the problem described in § 5.2.</figcaption></figure>

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).

**5.3. 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.

**6. 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.
<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /></figure>

<figure><img src="image_21.png" /></figure>

<figure><img src="image_22.png" /></figure>

<figure><img src="image_23.png" /></figure>
<figure><img src="image_24.png" /><figcaption>**Figure 16.** *Transversal cut of a textured Bi-2212 sample. Photo courtesy Instituto de Ciencia de Materiales de Aragón, CSIC.*</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>**Figure 17.** *Evolution of the melted zone geometry at time steps 1 (a) and 60 (b). The background triangulation is represented in wireframe.*</figcaption></figure>

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.
<figure><img src="image_26.png" /><figcaption>**Figure 18.** *Domain geometry.*</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>**Figure 19.** *Temperature profiles (temperature versus time) at the three different points indicated in Fig. 18.*</figcaption></figure>

## Acknowledgements

The work here presented has been supported by the Spanish Ministry of Education and Science, through project CICYT-DPI2002-01986. This support is gratefully acknowledged. Experimental results provided by the *Instituto de Ciencia de Materiales de Aragón*, CSIC, are also acknowledged.
7. References

[BAB 96] BABUŠKA I., MELENK J. M., “The partition of Unity Finite Element Method: Basic Theory and Applications”, *Comp. Meth. in Appl. Mech. and Eng.*, vol. 4, 1996, p. 289-314.

[CUE 00] CUETO E., DOBLARÉ M., GRACIA L., “Imposing essential boundary conditions in the Natural Element Method by means of density-scaled α-shapes.”, *International Journal for Numerical Methods in Engineering*, vol. 49-4, 2000, p. 519-546.

[CUE 03] CUETO E., SUKUMAR N., CALVO B., MARTÍNEZ M. A., CEGOÑINO J., DOBLARÉ M., “Overview and recent advances in Natural Neighbour Galerkin methods.”, *Archives of Computational Methods in Engineering*, vol. 10(4), 2003, p. 307–384.

[DEL 34] DELAUNAY B., “Sur la Sphère Vide. A la memoire de Georges Voronoi”, *Izvestia Akademii Nauk SSSR, Otdelenie Matematicheskii i Estestvennyka Nauk*, vol. 7, 1934, p. 793-800.

[EDE 94] EDELSBRUNNER H., MÜCKE E., “Three Dimensional Alpha Shapes”, *ACM Transactions on Graphics*, vol. 13, 1994, p. 43-72.

[GON 04a] GONZALEZ D., CUETO E., MARTINEZ M. A., DOBLARE M., “Numerical integration in Natural Neighbour Galerkin methods.”, *International Journal for Numerical Methods in Engineering*, vol. 60(12), 2004, p. 2077-2104.

[GON 04b] GONZÁLEZ D., CUETO E., DOBLARÉ M., “Volumetric locking in Natural Neighbour Galerkin methods.”, *International Journal for Numerical Methods in Engineering*, vol. 60(12), 2004, p. 2077-2104.

[HIY 99] HIYOSHI H., SUGIHARA K., “Two Generalizations of an Interpolant Based on Voronoi Diagrams”, *International Journal of Shape Modeling*, vol. 5, num. 2, 1999, p. 219–231.

[SIB 80] SIBSON R., “A Vector Identity for the Dirichlet Tesselation”, *Mathematical Proceedings of the Cambridge Philosophical Society*, vol. 87, 1980, p. 151-155.

[SIB 81] SIBSON R., “A Brief Description of Natural Neighbour Interpolation”, *Interpreting Multivariate Data*. V. Barnett (Editor), p. 21-36, John Wiley, 1981.

[STE 91] STEFAN J., “Über die Theorie der Eisbildung, inbesondere über die Eisbildung im Polarmee”, *Annalen der Physik und Chemie*, vol. 42, 1891, p. 269–286.

[SUK 98] SUKUMAR N., MORAN B., BELYTSCHKO T., “The Natural Element Method in Solid Mechanics”, *International Journal for Numerical Methods in Engineering*, vol. 43(5), 1998, p. 839-887.

[SUK 01] SUKUMAR N., MORAN B., SEMENOV A. Y., BELIKOV V. V., “Natural Neighbour Galerkin Methods”, *International Journal for Numerical Methods in Engineering*, vol. 50(1), 2001, p. 1-27.

[THI 11] THIESSEN A. H., “Precipitation averages for large areas.”, *Monthly Weather Report*, vol. 39, 1911, p. 1082–1084.

[VOR 08] VORONOI G. M., “Nouvelles Applications des Paramètres Continus à la Théorie des Formes Quadratiques. Deuxième Memoire: Recherches sur les parallélloèdres Primifs”, *J. Reine Angew. Math.*, vol. 134, 1908, p. 198-287.

[YVO 05] YVONNET J., CHINESTA F., RYCKELYNCK D., LORONG P., “ The constrained natural element method (C-NEM) for treating thermal models involving moving interfaces”, *International Journal of Thermal Sciences*, vol. 44, 2005, p. 559-569.