# The constrained natural element method (C-NEM) for treating thermal models involving moving interfaces

Julien Yvonnet, Francisco Chinesta, Philippe Lorong, David Ryckelynck

► To cite this version:

Julien Yvonnet, Francisco Chinesta, Philippe Lorong, David Ryckelynck. The constrained natural element method (C-NEM) for treating thermal models involving moving interfaces. International Journal of Thermal Sciences, 2005, 44 (6), pp.559 - 569. 10.1016/j.ijthermalsci.2004.12.007. hal-01508705

HAL Id: hal-01508705

https://hal.science/hal-01508705

Submitted on 15 Apr 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.

Public Domain
# The constrained natural element method (C-NEM) for treating thermal models involving moving interfaces

J. Yvonnet, F. Chinesta, Ph. Lorong, D. Ryckelynck

LMSP, UMR 8106 CNRS ENSAM-ESEM, 151 boulevard de l'Hôpital, F-75013 Paris, France

Received 29 May 2004; received in revised form 29 November 2004; accepted 13 December 2004

In this paper, the features of the natural neighbor (Sibson) interpolant are used within the context of a constrained Voronoi diagram, dual to the constrained Delaunay triangulation, for treating moving interface (Stefan) problems. The constrained natural element method (C-NEM) uses the Voronoi cells instead of the Delaunay triangles for both interpolation and integration, and thus permits the use of very distorted triangles in the Delaunay triangulation without loss of accuracy. The resulting interpolation is $C^1$ ($C^\infty$ in the most part of the domain) everywhere except at the nodes and across the interface, where it is $C^0$, which allows an accurate treatment of models involving moving interfaces. Numerical experiments are performed for proving the accuracy of the method, as well as its promise.

**Keywords:** Natural neighbor interpolation; C-NEM; Constrained Voronoi diagram; Moving interface; Stefan problem

## 1. Introduction

Phase boundaries represent material interfaces across which several fields may exhibit sharp gradients, and even discontinuities. A wide range of numerical methods have been developed for treating these problems according to the pertinent physics and assumptions about the interface [12]. When a sharp interface is considered, its motion is governed by the jump in the temperature gradient normal to the phase boundary and is accompanied by latent heat effects (Stefan condition). In order to satisfy these conditions the most common approach lies in explicitly tracking the interface motion. Within the interface tracking approach two main alternatives exist: the moving mesh methods and the mixed Eulerian–Lagrangian methods. Moving finite element mesh methods conform element boundaries to the interface as it evolves. Although these methods are very accurate, they are limited by severe mesh distortion. Thus, frequent remeshing is needed, with the associated field projections between

successive meshes. Moreover, remeshing is, even today, a delicate task in three dimensions. To alleviate remeshing efforts a number of Eulerian–Lagrangian methods have been developed recently that track the interface while solving the equations on a fixed grid [28]. Many of these methods effectively smear the discontinuity over a few grid cells, and are therefore not capable of representing the true discontinuity across the interface.

A new approach for representing localized behaviours has recently emerged in the field of the finite element method, known as the partition of unity method [16]. The main idea is to extend the classical approximation considering the product of the standard shape functions and local enrichment functions. The extended finite element method (X-FEM) is a variation on this technique. Recently, the X-FEM has been coupled to the Level Set Method [27] to represent interface [25]. In this way, the discontinuity evolution can be properly represented on a fixed background mesh, just by adding an appropriate enrichment in the functional approximation in the elements that are intersected by the moving discontinuity [8]. However, when the material in which the interface is moving, is subjected to large displacements, an
Nomenclature

<table><tr><td>T</td><td>temperature</td></tr><tr><td>c</td><td>volumetric heat capacity</td></tr><tr><td>k</td><td>thermal conductivity</td></tr><tr><td>V</td><td>interface velocity</td></tr><tr><td>|[q]|</td><td>thermal flux jump across the interface</td></tr><tr><td>L</td><td>volumetric latent heat of fusion</td></tr><tr><td>T<sub>m</sub></td><td>melting temperature</td></tr></table>

<table><tr><td>β</td><td>thermal diffusivity</td></tr><tr><td>η</td><td>ratio of thermal diffusivities</td></tr><tr><td>Q</td><td>heat sink intensity</td></tr><tr><td>erf</td><td>error function</td></tr><tr><td>erfc</td><td>complementary error function</td></tr><tr><td>Ei</td><td>exponential integral</td></tr></table>

updated Lagrangian description could be a better choice. When the background mesh evolves, remeshing will also be required to avoid large mesh distortions. In the context of moving phase boundary problems, it has been shown in [8] that due to the $C^0$ finite element shape functions, errors are introduced when the moving interface crosses the edge of the elements.

To alleviate dependence on the mesh and to provide smoother shape functions, the use of meshfree or meshless methods has been investigated. Meshless methods discretize a continuum body by a finite number of particles (or nodes) and the field of interest is interpolated under these nodes without the aid of an explicit mesh. Many meshless methods have been proposed, including the smooth particle hydrodynamics (SPH) [14], the Radial Basis Function methods [9, 10] the diffuse element method (DEM) [18], the element free Galerkin (EFG) [2], the reproducing kernel particle method (RKPM) [13], the HP clouds [6] and the partition of unity method (PUM) [16].

The introduction of moving discontinuities in these meshless methods can present difficulties for the following reasons:

(i) The quality of the approximation as well as the conditioning of the global system is pathologically dependent on the size of the support of the shape functions;

(ii) Imposition of essential boundary conditions needs special treatment;

(iii) Integration is not accurate enough; and

(iv) The physical discontinuity across the interfaces must be introduced accurately.

In order to overcome these different problems, we propose the use of the constrained natural element method (C-NEM) [30,31] for treating thermal models involving moving interfaces. This approach is an extension of the natural element method [26] in which both trial and test functions are constructed on the basis of the Voronoi based interpolants [7,24]. These interpolants satisfy the Kronecker delta property and their support is defined by the union of the Delaunay spheres passing through the visible nodes. In the C-NEM, the introduction of a visibility criterion and its related constrained Voronoi diagram preserves the appealing

properties of the NEM in any geometry (convex or not) and allows the introduction of material discontinuities.

Furthermore, the C-NEM is strongly related to the finite element method through the following points: (a) an underlying structure is needed (the constrained Delaunay triangulation/constrained Voronoi diagram); (b) the external boundary conditions can be imposed directly, due to the interpolant character and the strict linearity of the shape functions along any kind of boundary or interface. Nevertheless, in the C-NEM, both interpolation and integration are not constructed on the basis of the Delaunay triangles, but on its dual, the Voronoi diagram (which can be extended to constrained Delaunay triangulation [21] for non-convex domains or when they involve fixed or moving discontinuity), which provides smoother approximation and allows to proceed on the basis of very distorted Delaunay triangles. In a recent paper [31] the C-NEM has been successfully applied in some problems involving non-convex domains and cracks discontinuities [31].

The layout of this paper is as follow: in Section 2, we introduce a simple mathematical model of a thermal problem involving a moving interface (Stefan problem [29]). In Section 3, the C-NEM is summarized, and it will be applied in Section 4 to discretize the weak formulation of the nonlinear Stefan problem. A numerical benchmark is presented in Section 5, which allows us to test the accuracy of the proposed technique by comparing the numerical results with the exact solution of the problem.

## 2. Problem formulation

Let $\Omega \in \mathbb{R}^2$ be a bounded domain and $T$ the temperature field. On the domain boundary $\Gamma = \partial\Omega$ the temperature or the thermal flux is prescribed. We will denote by $\Gamma_1$ the domain boundary where the temperature is known $T(x \in \Gamma_1, t) = \bar{T}(x, t)$ and by $\Gamma_2$ the domain boundary where the heat flux $\bar{q}$ is imposed. The thermal model is defined in the time interval $[0, t_{max}]$. The initial temperature $T(x, t=0) = T_0$, where $T_0$ is assumed to be higher than the material melting temperature $T_m$. At time, $t=0$, a part of the domain boundary $\Gamma_1$ is suddenly exposed to a temperature $T_1 < T_m$. A moving solidification front $\Gamma_I$ is then generated, whose position evolves in time, i.e., $\Gamma_I(t)$, dividing the domain $\Omega$ in two regions $\Omega_1(t)$ (containing the solid
<figure><img src="image_1.png" /><figcaption>Fig. 1. Two phases problem.</figcaption></figure>

phase at time $t$) and $\Omega_2(t)$ (which contains the liquid phase) as shown in Fig. 1. For the sake of simplicity we will consider, from now on, a homogeneous and isotropic thermal model in both phases.

The heat transfer model is defined in each phase, neglecting volumetric source terms as well as the motion in the liquid phase induced by thermal gradients, by:

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

where $c_1$ and $c_2$ are the volumetric heat capacities of both phases, and $k_1$ and $k_2$ are the respective thermal conductivities. The associated initial and boundary conditions 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, \forall t \in [0, t_{\max}] \\ -k \nabla T(\mathbf{x}, t) \cdot \mathbf{n} = \bar{q}(\mathbf{x}, t) & \forall \mathbf{x} \in \Gamma_2, \forall t \in [0, t_{\max}] \end{array} \right. \quad (2) $$

The evolution of the interface $\Gamma_I(t)$ is described by a Stefan condition:

$$ \mathbf{V}(\mathbf{x} \in \Gamma_I(t)) = \frac{|[q]|}{L} \mathbf{n}_I(\mathbf{x}) \quad (3) $$

where $\mathbf{V}$ is the interface velocity, $L$ is the volumetric latent heat of fusion, $\mathbf{n}_I(\mathbf{x})$ is the normal vector to the interface at point $\mathbf{x}$ which is assumed to point into the liquid phase, and $|[q]|$ the thermal flux jump across the interface $\Gamma_I(t)$, i.e.,

$$ |[q]| = (k_1 \nabla T|_{\Gamma_I^-(t)} - k_2 \nabla T|_{\Gamma_I^+(t)}) \mathbf{n}_I \quad (4) $$

The additional constraint prescribed on the interface $\Gamma_I(t)$ is:

$$ T(\mathbf{x}, t) = T_m; \quad \forall \mathbf{x} \in \Gamma_I(t) \quad (5) $$

where $T_m$ is the melting temperature.

### 3. The constrained natural element method (C-NEM)

In this section, the utility of the C-NEM to describe moving interfaces and discontinuities in a fixed cloud of nodes is discussed. After a brief review of the Voronoi-based interpolants, we introduce the constrained Voronoi diagram which is used for computing the shape functions in any domain.

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

<figure><img src="image_3.png" /><figcaption>Fig. 2. Construction of the Sibson shape functions.</figcaption></figure>

#### 3.1. Natural neighbor interpolation

We briefly touch upon the foundation of Sibson's natural neighbor coordinates (shape functions) that are used in the natural element method. For a more in-depth discussion on the Sibson interpolant and its application for solving second-order partial differential equations, the interested reader can refer to Sambridge et al. [19], and Sukumar et al. [26]. The NEM interpolant is constructed on the basis of the Voronoi diagram. The Delaunay tessellation is the topological dual of the Voronoi diagram.

Consider a set of nodes $S = \{n_1, n_2, \dots, n_N\}$ in $\mathbb{R}^2$. The Voronoi diagram is the subdivision of $\mathbb{R}^2$ into regions $T_i$ (Voronoi cells) defined by:

$$ T_i = \{ \mathbf{x} \in \mathbb{R}^2: d(\mathbf{x}, \mathbf{x}_i) < d(\mathbf{x}, \mathbf{x}_j), \forall j \neq i \}, \quad \forall i \quad (6) $$

The Sibson coordinates of $\mathbf{x}$ with respect to a natural neighbor $n_i$ (see Fig. 2) is defined as the ratio of the overlap area (volume in 3D) of their Voronoi cells to the total area (volume in 3D) of the Voronoi cell related to point $\mathbf{x}$:

$$ \phi_i(\mathbf{x}) = \frac{\text{Area}(\text{afghe})}{\text{Area}(abcde)} \quad (7) $$

If the point $\mathbf{x}$ coincides with the node $n_i$, i.e., $\mathbf{x} = \mathbf{x}_i$, $\phi_i(\mathbf{x}_i) = 1$, and all other shape functions are zero, i.e., $\phi_j(\mathbf{x}_j) = \delta_{ij}$ ($\delta_{ij}$ being the Kronecker delta). The properties of positivity, interpolation, and partition of unity are then verified [26]:

$$ \left\{ \begin{array}{l} 0 \le \phi_i(\mathbf{x}) \le 1 \\ \phi_i(\mathbf{x}_j) = \delta_{ij} \\ \sum_{i=1}^{n} \phi_i(\mathbf{x}) = 1 \end{array} \right. \quad (8) $$
The natural neighbor shape functions also satisfy the local coordinate property [24], namely:

$$ \mathbf{x} = \sum_{i=1}^{n} \phi_i(\mathbf{x})\mathbf{x}_i \qquad (9) $$

which combined with Eq. (8) implies that the natural neighbor interpolant spans the space of linear polynomials (linear completeness).

Natural neighbor shape functions are $C^\infty$ at any point except at the nodes, where they are only $C^0$, and on the boundary of the Delaunay circles (spheres in 3D) where they are only $C^1$, because of the discontinuity in the neighbors nodes across these boundaries. Hiyoshi and Sugihara [7] have shown that the Sibson interpolant belongs to a more general class of Voronoi-based interpolants, called kth order standard coordinates, proving that the interpolant generated by the kth order standard coordinates have $C^k$ continuity on the Delaunay circles (spheres) boundaries. In this context, the Sibsonian and non-Sibsonian (Laplace) coordinates [1] results to be the standard coordinates of order 1 and 0, respectively.

Another important property of this interpolant is the ability to reproduce linear functions over the boundary of convex domains. The proof can be found in Sukumar et al. [26]. An illustration is depicted in Fig. 2(b): as the areas associated to points on the boundary become infinite, the contribution of internal points vanish in the limit when the point approaches the convex boundary, and the shape functions associated with nodes $n_1$ and $n_2$ become linear on the segment ($n_1-n_2$). This is not true in the case of non-convex boundaries, and the next section focuses on an approach to circumvent this difficulty.

Consider an interpolation scheme for a scalar function $T(\mathbf{x}) : \Omega \subset \mathbb{R}^2 \to \mathbb{R}$, in the form:

$$ T^h(\mathbf{x}) = \sum_{i=1}^{n} \phi_i(\mathbf{x})T_i \qquad (10) $$

where $T_i$ are the nodal temperatures at the n natural neighbor nodes, and $\phi_i(\mathbf{x})$ are the shape functions associated with each neighbor node. It is noted that Eq. (10) defines a local interpolation scheme. Thus, the trial and test functions used in the discretization of the variational formulation describing the thermal problem treated in this paper take the form of Eq. (10).

## 3.2. The constrained natural element method

### 3.2.1. Constrained Voronoi diagram

It was proved in [26,31] and [5] that spurious influences between “non-visible” nodes and loss of linearity in the interpolation along boundaries of non-convex domains appear in the framework of the NEM. In order to avoid this drawback and to recover all properties of the method for any geometry (including non-convex domains containing cracks or involving field discontinuities), a visibility criterion is

<figure><img src="image_4.png" /><figcaption>Fig. 3. Interpretation of the constrained Voronoi diagram. (a) Schematic view of the CVD. (b) Constrained Delaunay triangulation with respect to segment $l_i$ and intersection between the dual CVD and the domain closure.</figcaption></figure>

introduced in order to restrict influent nodes among natural neighbors. The computation of the shape functions is done on the basis of the so-called constrained (or extended) Voronoi diagram (CVD), introduced by Seidel in [21]. Imagine $\mathbb{R}^2$ is a sheet of paper $\Sigma_0$, with the points of the set of nodes $S$ and the line segments defining the boundary in a set $L$, drawn on it. For each $l_i \in L$, we cut $\Sigma_0$ open along $l_i$ and glue another sheet $\Sigma_i$, which also cut open along $l_i$. The gluing is done around $l_i$ such that every traveler who crosses $l_i$ switches from $\Sigma_0$ to $\Sigma_i$ and vice versa. A schematic view of the particular gluing necessary to achieve that effect is illustrated in Fig. 3(a). We know what it means for two points on the primary sheet to be visible from each other. For other pairs we need a more general definition. For $i \neq 0$, points $x_0 \in \Sigma_0$ and $y_i \in \Sigma_i$ are visible if $xy$ crosses $l_i$, and $l_i$ is the first constraining line segment crossed if we traverse $xy$ in the direction from $x$ to $y$.

In Fig. 3(b), the intersection between the CVD and the domain closure is depicted. The resulting diagram is composed of cells $T_i^C$, one for each node $n_i$, such that any point $\mathbf{x}$ inside $T_i^C$ is closer to $n_i$ than to any other node $n_j$ visible from $\mathbf{x}$. We call this the constrained Voronoi cells, which are defined formally by:

$$ T_i^C = \{ \mathbf{x} \in \mathbb{R}^n : d(\mathbf{x}, \mathbf{x}_i) < d(\mathbf{x}, \mathbf{x}_j), \forall j \neq i, S_{\mathbf{x} \to n_i} \cap \Gamma = \emptyset, \\ S_{\mathbf{x} \to n_j} \cap \Gamma = \emptyset \} \qquad (11) $$
<figure><img src="image_5.png" /><figcaption>Fig. 4. Reproducing discontinuous derivatives using the constrained Voronoi diagram.</figcaption></figure>

where $\Gamma$ is the domain boundary, composed with segments $l_i \in L$ and $S_{a \to b}$ denotes the segment between the points a and b.

The constrained Delaunay triangulation does not always exist in 3D without adding new nodes, as shown in [20]. Nevertheless, some techniques for constructing 3D constrained Delaunay tessellations are available and provided in [22,23] by addition of Steiner points.

### 3.2.2. The constrained natural element approximation

In order to solve partial differential equations defined on non convex domains, or to reproduce functional discontinuities, we consider the following approximation for both the trial and the test functions:

$$T^h(\mathbf{x}) = \sum_{i=1}^{V} \phi_i^C(\mathbf{x}) T_i \quad (12)$$

where V is the number of natural neighbors visible from point x and $\phi_i^C$ is the constrained natural neighbor shape function related to the i-th node at point x. The computation of the C-n-n (constrained natural neighbor) shape functions is similar to the natural neighbor shape function, when one proceeds using the constrained Voronoi diagram introduced previously. It was shown in [30] and [31] that the use of the constrained Voronoi diagram does not affect the properties of the NEM interpolation, allowing the extension of the linearity of the shape functions on the convex domains boundaries, to any geometry, convex or not.

The ability of the C-NEM for treating problems involving cracks has been illustrated in [31]. In the present paper, we focus on its application in the context of a moving interface defining two domains with different thermal properties. Thus, defining at time t two CVD (constrained Voronoi diagrams) of $\Omega_1(t)$ and $\Omega_2(t)$, both with respect to the interface $\Gamma_I(t)$, it can be proved that the interpolated temperature field is $C^1$ everywhere, except at the nodes and on the interface $\Gamma_I(t)$ where it is only $C^0$. Thus, this interpolation seems to be appropriate to simulate the Stefan problem considered in this paper.

To illustrate this behavior, we consider the situation depicted in Fig. 4, where the point x moves from $\Omega_1$ to $\Omega_2$. If x is in $\Omega_1$, the interpolated field is constructed from Eq. (12) using the neighbor visible nodes from point x ($\Gamma_I$ is assumed opaque). If x is on $\Gamma_I$, according to the previous discussion,

the interpolated field is strictly linear because it only depends on the two neighbor nodes located on $\Gamma_I$. Finally, when x is in $\Omega_2$, the interpolated field is defined using the visible neighbor and visible nodes from point x ($\Gamma_I$ being opaque). The continuity of the interpolated field is then guaranteed, but a discontinuity appears in the field derivatives, because of a sudden change in the neighbor nodes across the interface. We can then reproduce the temperature field continuity, as well as the expected flux discontinuity on the interface.

## 4. C-NEM discretization

The weak formulation associated with Eq. (1) results:

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

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

where $H^1(\Omega)$ and $H_0^1(\Omega)$ are the usual Sobolev functional spaces. Substituting the trial and test functions (both approximated in the C-NEM framework) in the above equation and using the arbitrariness of the field $\delta T$, the following system of equations is obtained:

$$C\bar{T} + KT = F \quad (14)$$

where T is the vector containing the unknown nodal temperatures. We consider the solution on the time interval $[0, t_{max}]$, partitioned into steps as $[t^n, t^{n+1}]$ and the generalized trapezoidal time stepping algorithm characterized by the parameter $\alpha$:

$$\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 (15)$$

which leads to:

$$(C^{n+1} + \alpha \Delta t K^{n+1}) T^{n+1} = F^{n+1}(T^n, |[q]|^{n+1}) \quad (16)$$

with:

$$C^{n+1} = \int_{\Omega^{n+1}} N^t c N d\Omega \quad (17)$$
where $\Omega^{n+1} = \Omega_1^{n+1} \cup \Omega_2^{n+1}$;

$$ \mathbf{K}^{n+1} = \int_{\Omega^{n+1}} \mathbf{B}'k\mathbf{B} d\Omega \qquad (18) $$

and

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

where $\mathbf{N}$ is the vector containing the nodal shape functions:

$$ \mathbf{N} = \{ \phi_1 \ \phi_2 \ \dots \ \phi_N \} $$

and $\mathbf{B}$ is the matrix containing the shape functions derivatives:

$$ \mathbf{B} = \begin{Bmatrix} \phi_{1,x} & \phi_{2,x} & \dots & \phi_{N,x} \\ \phi_{1,y} & \phi_{2,y} & \dots & \phi_{N,y} \end{Bmatrix} $$

The stabilized conforming nodal integration proposed by Chen et al. in [4] is employed for the numerical integration of $\mathbf{K}$ (see our former work [31] for more details). A lumped mass matrix $\tilde{\mathbf{C}}$ is computed making use of the constrained Voronoi cells areas as nodal weights.

The iteration procedure is defined as:

Knowing $\mathbf{T}^n$ and $|[q]|^n$ at time $t^n$, the non-linear problem associated with Eq. (16) results in finding $\mathbf{T}^{n+1}$ and $|[q]|^{n+1}$ such that Eqs. (5) and (16) are satisfied. For this purpose we proceed as follows:

(1) Compute the interface velocity $\mathbf{V}^n(\mathbf{x})$ using Eq. (3) and update the interface position at time $t^{n+1}$ using the forward Euler formula:

$$ \mathbf{x}_J^{n+1} = \mathbf{x}_J^n + \Delta t \mathbf{V}^n(\mathbf{x}_J^n) \qquad (20) $$

where $\mathbf{x}_J$ are the nodes defining the interface.

(2) Update locally the constrained Voronoi diagram and the shape functions associated with integration points in the interface neighborhood. Then, we compute $\tilde{\mathbf{C}}^{n+1}$ and $\mathbf{K}^{n+1}$.

(3) Solve Eq. (16) using a Newton-Raphson procedure where the tangent matrix is computed numerically.

(4) Repeat while $t^{n+1} < t_{\max}$.

An alternative scheme using the Latin method [11] in the extended finite element framework can be found in Merle and Dolbow [17].

# 5. Numerical examples

## 5.1. Unidirectional solidification of a semi-infinite solid

In this section, we illustrate the potential of the proposed technique in simulating a two-phase Stefan problem. The

problem is essentially one-dimensional, but we solve it here in two dimensions to underscore the appealing features of the method.

The Stefan problem models the one-dimensional freezing of a semi-infinite domain ($x \ge 0$). The initial temperature $T_0$ is assumed constant in the whole domain, being higher than the melting temperature $T_m$. At time $t = 0$ the temperature at the left boundary $x = 0$ is suddenly prescribed to a value $T_1$ lower than the melting point, originating a solidification front that progresses from the boundary $x = 0$ in the $x$ direction. The exact flow front position $x_f(t)$ is given by:

$$ x_f(t) = 2\lambda\sqrt{\beta_s t} \qquad (21) $$

where $\beta_s = k_s/c_s$ is the thermal diffusivity of the solid phase, and the constant $\lambda$ satisfies the following relationship:

$$ \frac{e^{-\lambda^2}}{\operatorname{erf}(\lambda)} = \frac{k_l \sqrt{\eta} (T_0 - T_m) e^{-\eta \lambda^2}}{k_s (T_m - T_1) \operatorname{erfc}(\lambda \sqrt{\eta})} + \frac{\lambda L \sqrt{\pi}}{c_s (T_m - T_1)} \qquad (22) $$

with $\eta = \beta_s/\beta_l$ being the ratio of the thermal diffusivities and where $k_l$ represents the liquid phase conductivity. The temperature field in the solid phase $0 \le x \le x_f(t)$ is then:

$$ T(x, t) = T_1 + \frac{T_m - T_1}{\operatorname{erf}(\lambda)} \operatorname{erf}\left(\frac{x}{2\sqrt{\beta_s t}}\right) \qquad (23) $$

and in the liquid phase $x \ge x_f(t)$:

$$ T(x, t) = T_0 - \frac{T_0 - T_m}{\operatorname{erfc}(\lambda \sqrt{\eta})} \operatorname{erfc}\left(\frac{x}{2\sqrt{\beta_l t}}\right) \qquad (24) $$

where $\operatorname{erf}$ and $\operatorname{erfc}$ are the error function and complementary error function, respectively.

In the present investigation, we use the water-saturated sand thermal properties provided in [15] that are listed in Table 1. $T_1$ and $T_0$ were set to 263.15 K and 277.15 K, respectively ($\lambda = 0.3073$). We simulate the evolution of the temperature field in $\Omega = [0, 0.01] \times [0, 0.005]$ m. In order to use the infinite domain solution as reference solution, the temperature at $x = 0.01$ m is prescribed to its expected value according to Eq. (24).

In a first test, we consider in the domain $\Omega$ a $20 \times 10$ uniform grid and $\Delta t = 2$ s. In the following, a comparison between C-NEM and FEM has been performed. For this purpose, the FE shape functions have been computed on the basis of the constrained Delaunay triangulation, that can be obtained by connecting the natural neighbors in the constrained Voronoi diagram. Fig. 5 compares the computed interface position and the exact one. Excellent accuracy can be

<table><caption>Table 1
Thermal properties of the water saturated sand</caption><thead><tr><th>Properties</th><th>Solid</th><th>Liquid</th></tr></thead><tbody><tr><td>Volumetric heat capacity<br>(kg·m<sup>-1</sup>·s<sup>-2</sup>·K<sup>-1</sup>)</td><td>2.051 × 10<sup>6</sup></td><td>2.595 × 10<sup>6</sup></td></tr><tr><td>Thermal conductivity<br>(kg·m<sup>-1</sup>·s<sup>-3</sup>·K<sup>-1</sup>)</td><td>4.019</td><td>2.888</td></tr><tr><td>Melting temperature (K)</td><td></td><td>273.15</td></tr><tr><td>Volumetric latent heat of fusion<br>(kg·m<sup>-1</sup>·s<sup>-2</sup>)</td><td></td><td>8.038 × 10<sup>7</sup></td></tr></tbody></table>
<figure><img src="image_6.png" /><figcaption>Fig. 5. Computed FEM and C-NEM front position versus the exact solution using a 20 × 10 regular grid.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 6. Error in the front position.</figcaption></figure>

noticed for the C-NEM solution, as depicted in Fig. 6, where the error in the front position is represented. Fig. 7 shows the temperature profile at different times. The C-NEM accuracy is found greater than the FEM solution related to the same Delaunay triangles. In order to assess the thermal flux jump across the interface, the computed solution $|[q^h(t)]|$ is compared to the analytical solution $|[q^{ex}(t)]|$ given by:

$$|[q^{ex}(t)]| = \frac{L\lambda\beta_s}{\sqrt{\beta_s t}} \qquad (25)$$

Fig. 8 shows a comparison between the FEM and C-NEM solution for $\Omega = [0, 0.01] \times [0, 0.005]$ m containing 200 nodes, uniformly or randomly distributed. We can notice from Fig. 8 that C-NEM solution provides a higher accuracy compared to the FEM solution. Fig. 9 depicts the error between the C-NEM computed flux jump and the exact solution (25) for different refined meshes, proving the convergence of the approach. We point out that the discontinuity

<figure><img src="image_8.png" /><figcaption>Fig. 7. Temperature profiles along the line $y = 0.0025$ m using a regular grid.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Fig. 8. Error in flux: comparison FEM/C-NEM.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 9. Error in flux (C-NEM).</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Fig. 10. Computed interface position using an irregular cloud of nodes: (a) Cloud of nodes and interface position; (b) Constrained Voronoi cells.</figcaption></figure>

in the temperature gradient is accurately obtained, being in excellent agreement with the exact solution.

In a second test, we consider the domain $\Omega$ containing 200 nodes distributed at random. The purpose of the present test is to investigate the meshless feature of the technique, in which due to its meshless character no geometrical restrictions concerning the relative nodal positions are involved. Thus, neither the background nodal distribution nor the relative position of the nodes defining the moving interface with respect to the background nodes, induce a lack of accuracy when high distortions in the Delaunay mesh, used to compute the Voronoi diagram, takes place. This is the main difference between the proposed strategy and the standard finite element method whose accuracy depends significantly on the geometrical quality of the mesh. Moreover, this test approaches the situations encountered when the material is also moving, inducing highly irregular nodal densities and high background mesh distortions. Fig. 10 depicts the cloud of nodes and the interface position as well as the associated constrained Voronoi cells. Remarkably, despite the very irregular nodal distribution and density, we can notice that the interface computed in the C-NEM simulation is not distorted as it moves through the domain, unlike the FEM solution, from which we can conclude that the accuracy in the C-NEM is not significantly affected by the regularity of the nodal distribution. In Fig. 7 some temperature profiles along the line $y = 0.0025$ m are depicted, from which excellent accuracy can be inferred for the C-NEM solution.

In order to examine the convergence behaviour of the present approach, we conduct a series of calculations on increasingly refined meshes. The results for the relative error in the average front position for $t = 50$ s are provided in Fig. 11, which proves its convergence.

## 5.2. Solidification from a line heat sink

In this problem, a domain $\Omega = [-0.01, 0.01] \times [-0.01, 0.01]$ m initially entirely in a liquid state ($T_0 > T_m$) is exposed to a continuous line heat sink located at the point $(0, 0)$, originating the nucleation of an axisymmetric phase boundary separating the liquid and the solid phases whose
radius increases in time. This problem was investigated by Ji et al. in [8] in the context of the X-FEM method. The analytical solution for this problem defined in the unbounded $\Re^2$ can be found in [3]. The radius of the front is given by:

$$R_f = 2\lambda\sqrt{\beta_s t} \qquad (26)$$

<figure><img src="image_12.png" /><figcaption>Fig. 11. Relative error in the average front position for $t = 50$ s.</figcaption></figure>

The temperature in the solid region $r < R_f$ is given by:

$$T(r, t) = T_m + \frac{Q}{4\pi k_s} \left[ Ei\left(-\frac{r^2}{4\beta_s t}\right) - Ei(-\lambda^2) \right] \qquad (27)$$

and in the liquid region $r > R_f$ by:

$$T(r, t) = T_0 - \frac{T_0 - T_m}{Ei(\lambda^2 \eta)} Ei\left(-\frac{r^2}{4\beta_l t}\right) \qquad (28)$$

In the above, $\lambda$ is the root of the equation:

$$\frac{Q}{4\pi} e^{\lambda^2} = \lambda^2 \beta_s L - \frac{k_l (T_0 - T_m)}{Ei(-\lambda^2 \eta)} \qquad (29)$$

with $Q$ the heat sink intensity, $Ei$ the exponential integral and $\eta$ the ratio of thermal diffusivities:

$$\eta = \frac{\beta_s}{\beta_l} \qquad (30)$$

In our problem, only the upper quadrant $\Omega = [0, 0.01] \times [0, 0.01]$ of the domain is modeled, due to the symmetry. By imposing the exact solution on the external boundaries, we may model an axisymmetric problem in a square domain. It is noteworthy that the above solution is singular for $r = 0$. To overcome this difficulty, we move the node located at the origin to $(+h/2, +h/2)$, $h$ being the nodal distance between

<figure><img src="image_13.png" /><figcaption>Fig. 12. Evolution of the constrained Delaunay triangulation at times $t = 1.15$ s, $t = 2.65$ s and $t = 5.5$ s, and constrained Voronoi cells at time $t = 9.25$ s, using an unstructured nodal distribution.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Fig. 13. Comparison between computed and exact analytical interface position: (a) $t = 1.15$ s, (b) $t = 2.8$ s, (c) $t = 4.3$ s, (d) $t = 5.65$ s, (e) $t = 7.15$ s, (f) $t = 8.65$ s.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 14. Interface radius position.</figcaption></figure>

two nodes on the boundary. We consider in the numerical applications $Q = 4.186 \times 10^7$ kg·m$^{-1}$·s$^{-1}$ and $\lambda = 0.3513$. The thermal properties of the first example have been used here.

Simulations with regular and randomly distributed nodes are performed. In the first simulation, the domain contains 144 nodes on a uniform grid, and 8 additional nodes on the interface, which match the exact solution for $t_0 = 1$ s. Fig. 12 shows the evolution of the Delaunay triangles depicted at times $t = 1.15$ s, $t = 2.65$ s and $t = 5.5$ s. The triangulation is constrained by the moving interface. The dual constrained Voronoi cells are depicted at time $t = 9.25$ s.

In Fig. 13, the computed front is compared with the analytical solution. Higher accuracy is obtained using C-NEM in comparison with the FEM solution. Thus, an excel-

<figure><img src="image_16.png" /><figcaption>Fig. 15. Error in the interface location prediction.</figcaption></figure>

lent agreement with the analytical solution is obtained (see Figs. 14 and 15).

## 6. Conclusion

In this paper, the salient features of the C-NEM method are used for treating thermal problems involving moving interfaces. In the C-NEM framework, the Delaunay triangulation is used for its convenience in the frequent connectivity update during the motion of the interface. The use of a constrained Delaunay triangulation allows to conform the boundary with the interface and to accurately reproduce the jump in the heat flux across the interface. The traditional $C^0$ finite element shape functions are replaced in the Galerkin scheme by the smoother natural neighbor (Sibson) shape functions, which are constructed on the basis of the
dual to the Constrained Voronoi diagram (dual to the con-
strained Delaunay triangulation). A nodal integration is also
performed on the basis of the Voronoi cells by using this
dual structure instead of the triangles structure. This tech-
nique allows the use of very distorted Delaunay triangles
without loss of accuracy, and avoids the numerical problems
associated with the discontinuity across the interface. The
enrichment of the natural neighbor interpolation with dis-
continuous shape functions related to a level-set description
in the context of the partition of unity is a work in progress.

References

[1] V.V. Belikov, V.D. Ivanov, V.K. Kontorovich, S.A. Korytnik, A.Y. Semenov, The non-Sibsonian interpolation: a new method of interpolation of the values of a function on an arbitrary set of points, Comput. Math. Math. Phys. 37 (1) (1997) 9–15.

[2] T. Belytschko, Y.Y. Lu, L. Gu, Element-free Galerkin methods, Internat. J. Numer. Methods Engrg. 37 (1994) 229–256.

[3] H. Carslaw, J. Jaeger, Conduction of Heat in Solids, Clarendon Press, Oxford, 1959.

[4] J.S. Chen, C.T. Wu, Y. Yoon, A stabilized conforming nodal integration for Galerkin mesh-free methods, Internat. J. Numer. Methods Engrg. 50 (2001) 435–466.

[5] E. Cueto, J. Cegoñino, B. Calvo, M. Doblaré, On the imposition of essential boundary conditions in natural neighbor Galerkin methods, Comm. Numer. Methods Engrg. 19 (2003) 361–376.

[6] C.A. Duarte, J.T. Oden, An H-p adaptative method using clouds, Comput. Methods Appl. Mech. Engrg. 139 (1996) 237–262.

[7] H. Hiyoshi, K. Sugihara, Improving continuity of Voronoi-based inter-
polation over Delaunay spheres, Comput. Geom. 22 (2002) 167–183.

[8] H. Ji, D. Chopp, J.E. Dolbow, A hybrid finite element/level set method for modeling phase transformations, Internat. J. Numer. Methods 54 (2002) 1209–1233.

[9] E.J. Kansa, Multiquadratics—a scattered data approximation scheme with applications to computational fluid dynamics. I—Surface approximations and partial derivative estimates, Comput. Math. Appl. 19 (1990) 127–145.

[10] E.J. Kansa, Multiquadratics—a scattered data approximation scheme with applications to computational fluid dynamics. II—Solutions to hyperbolic, parabolic and elliptic partial differential equation, Comput. Math. Appl. 19 (1990) 147–161.

[11] P. Ladevèze, Non-Linear Computational Structural Mechanics, Springer, New York, 1998.

[12] R. Lewis, K. Ravindran, Finite element simulation of metal casting, Internat. J. Numer. Methods Engrg. 47 (2000) 29–59.

[13] W.K. Liu, S. Jun, Y.F. Zhang, Reproducing kernel particle methods, Int. J. Numer. Methods Fluids 21 (1995) 1081–1106.

[14] L.B. Lucy, A numerical approach to the testing of fusion process, As-
tron. J. 88 (1977) 1013–1024.

[15] D. Lynch, K. O'Neill, Continuously deforming finite elements for the solution of parabolic problems, with and without phase change, Internat. J. Numer. Methods Engrg. 17 (1981) 81–96.

[16] J.M. Melenk, I. Babuška, The partition of unity finite element method: Basic theory and applications, Comput. Methods Appl. Mech. Engrg. 139 (1996) 289–314.

[17] R. Merle, J.E. Dolbow, Solving thermal and phase change with the extended finite element method, Comput. Mech. 28 (5) (2002) 339–350.

[18] B. Nayroles, G. Touzot, P. Villon, Generalizing the finite ele-
ment method: Diffuse approximation and diffuse elements, Comput.
Mech. 10 (1992) 307–318.

[19] M. Sambridge, J. Braun, M. McQueen, Geophysical parameterization and interpolation of irregular data using natural neighbors, Geophys. J. Int. 122 (1995) 837–857.

[20] E. Schönhardt, Über die Zerlegung von Dreieckspolyedern in Tetraeder, Math. Ann. 98 (1928).

[21] R. Seidel, Constrained Delaunay triangulations and Voronoi diagrams with obstacles, in: 1978–1988 Ten Years IIG, 1988, pp. 178–191.

[22] J.R. Shewchuck, Tetrahedral mesh generation by Delaunay refinement, in: Proc. Fourteenth Annual Symposium on Computational Geometry, Minneapolis, Minnesota, ACM, 1998, pp. 86–95.

[23] J.R. Shewchuck, Sweep algorithms for constructing higher-
dimensional constrained Delaunay triangulations, in: Proc. Sixteenth
Annual Symposium on Computational Geometry, Hong-Kong, ACM,
2000, pp. 350–359.

[24] R. Sibson, A vector identity for the Dirichlet tessellations, Math. Proc. Cambr. Phil. Soc. 87 (1980) 151–155.

[25] N. Sukumar, D. Chopp, N. Mo's, T. Belytschko, Modeling holes and inclusions by level sets in the extended finite element method, Comput. Methods Appl. Mech. Engrg. 190 (2001) 6183–6200.

[26] N. Sukumar, B. Moran, T. Belytschko, The natural elements method in solid mechanics, Internat. J. Numer. Methods Engrg. 43 (1998) 839–887.

[27] M. Sussman, P. Smereka, S. Osher, A level set approach for comput-
ing solutions to incompressible two-phase flows, J. Comput. Phys. 114
(1994) 146–159.

[28] H. Udaykumar, R. Mittal, W. Shyy, Computation of solid–liquid phase fronts in the sharp interface limit on fixed grids, J. Comput. Phys. 153 (1999) 535–574.

[29] C. Vuik, Some Historical Notes about the Stefan Problem, CWI Tract 90, CWI, Amsterdam, 1993.

[30] J. Yvonnet, D. Ryckelynck, P. Lorong, P. Chinesta, Interpolation naturelle sur les domaines non convexes par l'utilisation du diagramme de Voronoi constraint—Méthode des éléments C-Naturels, Revue Européenne des éléments finis 12 (4) (2003) 487–509.

[31] J. Yvonnet, D. Ryckelynck, P. Lorong, F. Chinesta, A new extension of the natural element method for non-convex and discontinuous domains: the constrained natural element method (C-NEM), Internat. J. Numer. Methods Engrg. 60 (2004) 1451–1474.