# Enriched Reproducing Kernel Particle Approximation for Simulating Problems Involving Moving Interfaces: Application to Solidification Problems

Jean Trunzler, P. Joyot, Francisco Chinesta, Patrick Reuter

► To cite this version:

Jean Trunzler, P. Joyot, Francisco Chinesta, Patrick Reuter. Enriched Reproducing Kernel Particle Approximation for Simulating Problems Involving Moving Interfaces: Application to Solidification Problems. 8th ESAFORM Conference on Material Forming (ESAFORM 2005), Apr 2005, Cluj-Napoca, Romania. pp.97-100. [hal-00020869](https://hal.science/hal-00020869v1)

HAL Id: hal-00020869

https://hal.science/hal-00020869v1

Submitted on 12 Aug 2020

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>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Enriched Reproducing Kernel Particle Approximation For Simulating Problems Involving Moving Interfaces: Application To Solidification Problems

J. Trunzler¹, P. Joyot¹, F. Chinesta², P. Reuter¹

¹LIPSI-ESTIA, Technopole Izarbel 64210 Bidart

URL: www.estia.fr

e-mail: {j.trunzler,p.joyot,p.reuter}@estia.fr

²LMSP - UMR CNRS 8106, ENSAM-ESEM, 151 Bd. de l'Hôpital 75013 Paris

URL: www.paris.ensam.fr/LMSP

e-mail: francisco.chinesta@paris.ensam.fr

**ABSTRACT:** This paper presents a new method to simulate problems with moving interfaces. We use a meshless method where the approximation functions can reproduce discontinuities in the field derivatives. The position of the discontinuity is defined independent of the nodes that define the approximation. We apply our method to the propagation of a melting boundary.

**Key words:** ERKPM, Stefan problem, Numerical simulation, Meshless methods, Collocation techniques

## 1 INTRODUCTION

Finite element methods are commonly used to model melting problems where the problem domain is defined using a mesh support. When the domain evolves the mesh becomes very distorted, and remeshing as well as field projection between successive meshes are then required [5]. Moreover, to represent the thermal transmission conditions through the interface (melting or solidification fronts) the interface must be represented by facets (in 3D) or segments (2D) of the mesh skeleton, which involves continuous remeshing in the context of the FE approximations.

Recently, several authors proposed to use meshless methods in order to get rid of the mesh-related problems. In fact, the term meshless methods covers several different methods, and we classify these methods by their ability to deal with the solidification problem according to the three following criteria:

The first criterion is the way how the local discontinuous approximation is set up. In [2], Ji uses the Moving Least Squares (MLS) method associated with a Partition of Unity (PUM) enrichment, and in [8], Yvonnet uses the constrained natural method (CNEM). In a prior work, we proposed to use the enriched reproducing kernel particle approximation (ERKPA)[4, 3].

The second criterion is the way in which the problem equations are discretized. There are two possible ways. First, the weak formulation of the problem is used as in the EFG [1] or XFEM method. But these two cases

require the use of an integration mesh and thus reduce the meshless character of the methods. Second, the strong formulation of the problem is used to discretize the problem as in the point collocation methods. In this paper, we use this last approach.

The third criterion is the way how the solidification front is defined. The boundary can either be defined by using a mesh as in the CNEM method proposed by Yvonnet [8], or the front is defined totally independent from the nodes that define the approximation. For this latter case, Ventura [7] or Ji [2] proposed to use Level Set Methods [6].

## 2 PROBLEM STATEMENT

Consider the biphasic Stefan problem in 1D. The domain $\Omega$ with boundary $\Gamma$ consists of a pure solid or liquid material. The domain $\Omega$ can thus be subdivided into two disjoint subdomains $\Omega_S$ and $\Omega_L$, where $S$ and $L$ denote whether the material is in liquid or solid state. The position of the interface between liquid and solid is given by the coordinate $x_d$. Moreover, we consider the material to be homogeneous and isotropic, and it does not contain any internal source. The heat transfer problem results:

$$C_S \frac{\partial T}{\partial t} = \nabla \cdot (k_S \nabla T), \quad x \in \Omega_S \quad (1)$$

$$C_L \frac{\partial T}{\partial t} = \nabla \cdot (k_L \nabla T), \quad x \in \Omega_L \quad (2)$$
where $C_S$ and $C_L$ denote the volumetric heat capacities, and $k_S$ (resp. $k_L$) the thermal conductivity of the material in $\Omega_S$ (resp. $\Omega_L$).

The evolution of the discontinuity is given by Stefan's law:

$$LV(x_d) = - \left[ k_L \frac{\partial T}{\partial x} \bigg|_{x_d^-} - k_S \frac{\partial T}{\partial x} \bigg|_{x_d^+} \right] \quad (3)$$

The temperature on the interface is defined by

$$T = T_m, \quad x = x_d, \quad (4)$$

where $T_m$ is the melting temperature, $L$ the volumetric latent heat of fusion and $V$ the interface velocity.

For solving the thermal problem we need to prescribe the initial material temperature:

$$T(x, t = 0) = T^0(x), \quad x \in \Omega \quad (5)$$

In this work, we only impose one Dirichlet condition on the boundary $\Gamma$.

$$T(x, t) = T_{\Gamma}(x, t), \quad x \in \Gamma, \quad t \in [0, t_{\max}] \quad (6)$$

Thus, the thermal problem can be summarized as follows:

Find $T: (x, t) \to T(x, t)$ such that

$$\left\{
\begin{array}{l}
C_S \frac{\partial T}{\partial t} = \vec{\nabla} \cdot (k_S \vec{\nabla} T), \quad x \in \Omega_S \\
C_L \frac{\partial T}{\partial t} = \vec{\nabla} \cdot (k_L \vec{\nabla} T), \quad x \in \Omega_L \\
T = T_m, \quad x = x_d \\
T(x, t) = T_{\Gamma}(x, t), \quad x \in \Gamma, \quad t \in [0, t_{\max}]
\end{array}
\right.
\quad (7)$$

The Stefan's law (Eq. (3)) allows the interface location updating.

## 3 TIME DISCRETIZATION

We approximate the time derivative according to

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

that inserted into Eq. (7), leads to the following implicit scheme:

$$\left\{
\begin{array}{l}
C_S \frac{T^{n+1}-T^n}{\Delta t} = \vec{\nabla} \cdot (k_S \vec{\nabla} T^{n+1}), \quad x \in \Omega_S \\
C_L \frac{T^{n+1}-T^n}{\Delta t} = \vec{\nabla} \cdot (k_L \vec{\nabla} T^{n+1}), \quad x \in \Omega_L \\
T^{n+1} = T_m, \quad x = x_d \\
T^{n+1}(x) = T_{\Gamma}^{n+1}(x), \quad x \in \Gamma
\end{array}
\right.
\quad (9)$$

The temperature at the time $t = t^n = n\Delta t$ is approximated by the ERKPA defined by

$$T^n(x) = \sum_{i=1}^{N} \Psi_i^n(x) T_i^n. \quad (10)$$

The shape functions $\Psi_i^n(x)$ are defined with respect to the interface position $x_d^n$ at time $t^n = n\Delta t$, and $x_i$ ($i = 1 \cdots N$) denotes the coordinates of the approximation nodes that are considered to be fixed.

## 4 SPATIAL DISCRETIZATION

The coefficients $T_i^{n+1}$ have to be determined at every time step. For this purpose, we use a collocation method. We will see that there are two ways to impose the melting temperature $T_m$ at the liquid-solid interface:

* In the first one we impose the condition $T = T_m$ by penalization,

* In the second one we use an iteration procedure for solving the resulting non linear problem as described later.

### 4.1 Initialization

At the time $t = 0$, the coefficients $T_i^0$ have to satisfy the following conditions:

$$T^0(x_j) = \sum_{i=1}^{N} \Psi_i^0(x_j) T_i^0, \quad \text{for } j = 1 \cdots N \quad (11)$$

resulting in the following linear system:

$$\left[
\begin{array}{c}
\Psi_1^0(x_1) \quad \cdots \quad \Psi_N^0(x_1) \\
\vdots \quad \vdots \quad \vdots \\
\Psi_1^0(x_N) \quad \cdots \quad \Psi_N^0(x_N)
\end{array}
\right]
\left[
\begin{array}{c}
T_1^0 \\
\vdots \\
T_N^0
\end{array}
\right]
=
\left[
\begin{array}{c}
T^0(x_1) \\
\vdots \\
T^0(x_N)
\end{array}
\right]
\quad (12)$$

By solving this system, we obtain the coefficients $T_i^0$.
## 4.2 Incremental procedure

Given the coefficients $T_i^n$ and the interface position at time $t^n x_d^n$, we proceed as follows:

* The interface velocity $V^n$ is determined by

$$V^n(x_d^n) = -\frac{1}{L} \left[ k_L \frac{\partial T^n}{\partial x} \Big|_{x_d^-} - k_S \frac{\partial T^n}{\partial x} \Big|_{x_d^+} \right]. \quad (13)$$

* We determine the interface position $x_d^{n+1}$ at time $t^{n+1}$ using the explicit formula:

$$x_d^{n+1} = x_d^n + V^n \Delta t \quad (14)$$

* From the positions $x_d^n$ and $x_d^{n+1}$, we construct the shape functions $\Psi_i^n(x)$ and $\Psi_i^{n+1}(x)$, $\forall i$.

* Then, we determine the coefficients $T_i^{n+1}$ satisfying Eq. (9). To solve it, there are two different approaches that we describe in the following paragraphs.

## 4.3 Penalization method

Let $\Sigma_n$ be the set of the $N-1$ nodes that are the farthest away from the interface and let $x_{i_d}$ the coordinate of the node closest to the interface. By using a collocation discretization of Eq. (9), we obtain the following linear equations:

If $x_j \in \Gamma$:

$$\sum_{i=1}^{N} \Psi_{i}^{n+1}(x_j) T_{i}^{n+1} = T_{\Gamma}(x_j, t_{n+1}) \quad (15)$$

for $x_j \in \Sigma_n, x_j \notin \Gamma$:

$$\sum_{i=1}^{N} C(x_j) \Psi_{i}^{n+1}(x_j) T_{i}^{n+1} - k(x_j) \Delta t \Psi_{i,xx}^{n+1}(x_j) T_{i}^{n+1} = C(x_j) \sum_{i=1}^{N} \Psi_{i}^{n}(x_j) T_{i}^{n} \quad (16)$$

for the node $x_{i_d}$:

$$\sum_{i=1}^{N} C(x_{i_d}) \Psi_{i}^{n+1}(x_{i_d}) T_{i}^{n+1} - k(x_{i_d}) \Delta t \Psi_{i,xx}^{n+1}(x_{i_d}) T_{i}^{n+1} + \alpha \left( \sum_{i=1}^{N} \Psi_{i}^{n+1}(x_{d}^{n+1}) T_{i}^{n+1} - T_{m} \right) = C(x_{i_d}) \sum_{i=1}^{N} \Psi_{i}^{n}(x_{i_d}) T_{i}^{n} \quad (17)$$

We impose the condition $T(x_d) = T_m$ by penalization, being $\alpha$ the penalization coefficient. By solving this system, we obtain the coefficients $T_i^{n+1}$.

## 4.4 Iterative method

Another technique for prescribing the condition $T(x_d) = T_m$ lies in changing Eq. (17) by:

$$\int_{\Omega_{SD}} C(x) \left( \sum_{i=1}^{N} \Psi_{i}^{n+1}(x) T_{i}^{n+1} - \Psi_{i}^{n}(x) T_{i}^{n} \right) dx - \Delta t \left( \sum_{i=1}^{N} (k_L \Psi_{i,x}^{n+1}(x_{SDf}) - k_S \Psi_{i,x}^{n+1}(x_{SDi})) T_{i}^{n+1} \right) + |[q]|^{n+1} = 0$$

where $x_{SDf}$ and $x_{SDi}$ are the coordinates of the boundaries of the subdomain $\Omega_{SD}$, and $|[q]|$ denotes the jump of heat flux through the interface defined by:

$$|[q]|^{n+1} = \left( k_S \left. \frac{\partial T^{n+1}}{\partial x} \right|_{x_d^{n+1}-} - k_L \left. \frac{\partial T^{n+1}}{\partial x} \right|_{x_d^{n+1}+} \right) \quad (18)$$

In order to solve the resulting non-linear problem, we apply an iterative Newton technique, where the residual $R$ is defined:

$$R = T^{n+1}(x_d^{n+1}) - T_m \quad (19)$$

# 5 NUMERICAL EXAMPLE

Consider the domain $\Omega = [0, 1]$, with an initial temperature $T^0$ greater than the melting temperature $T_m$. At time $t_0$, we impose a temperature $T_1$ lower than $T_m$ at the boundary $x = 0$. This leads to a generation and propagation of a solidification front. See [8] for the analytic solution of the interface evolution.

In our numerical example, we consider the domain $\Omega$ filled with water at the initial temperature $T^0 = 4^\circ C$. The thermal characteristics of water can be found in [8]. At time $t_0$, we impose the temperature $T_1 = -10^\circ C$ on the boundary $x = 0$. In order to solve the numerical problem, we will impose the analytical solution at the boundary $x = 1$. The approximation of the temperature $T$ is built from 40 nodes that are uniformly distributed, and we are using a time step of $\Delta t = 2s$.

We can notice that the results obtained by using the penalization and the iterative method are very close.
<figure><img src="image_2.png" /><figcaption>Fig. 1. Temperature field $T(x)$ at $t = 20 \text{ s (a)}, t = 32 \text{ s (b)}, t = 64 \text{ s (c)}$.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 2. Evolution of the interface position $x_f(t)$.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 3. Error on the interface position predicted by using the penalization (left) and the iteration (right) techniques.</figcaption></figure>

The error related to the predicted interface position is lower than 1% (see Figs. 3a and 3b). Figure 1 proves that we reproduce the temperature field accurately. The enrichment of the approximation function allows then capturing the change of the temperature slope across the interface.

## 6 CONCLUSION

In this work we have solved a 1D thermal problem involving a moving interface by using a collocation method in combination with an ERKPA approximation. This technique can also be applied in 2D by using for example a Fast Marching method for the interface updating.

## REFERENCES

1. T. Belytschko, Y. Kronggauz, D. Organ, and M. Fleming. Meshless methods: an overview and recent developments. *Comput. Meths. Appl. Mech. Engrg.*, 139:3–47, 1996.

2. H. Ji, D. Chopp, and J.E. Dolbow. A Hybrid Extended Finite Element / Level Set Method for Modeling Phase Transformations. *Int. J. Numer. Meths. Eng.*, 54(8):1209–1233, 2002.

3. P. Joyot, J. Trunzler, and F. Chinesta. Discontinuous derivative enrichment in reproducing kernel particle approximations (rkpa). In the 7th International ESAFORM Conference on Material Forming, pages 69–72. Norwegian University of Science and Technology, Sigurd Stören, ISBN 82-92499-02-04, 2004.

4. P. Joyot, J. Trunzler, and F. Chinesta. Enriched reproducing kernel approximation: Reproducing functions with discontinuous derivatives. In Marc Alexander Griebel, Michael; Schweitzer, editor, *Meshfree Methods for Partial Differential Equations II*, volume 43 of *Lecture Notes in Computational Science and Engineering*, pages 93–108, 2004. ISBN: 3-540-23026-2.

5. R. Lewis and K. Ravindran. Finite element simulation of metal casting. *Int. J. Numer. Meths. Eng.*, 47:29–59, 2000.

6. J.A. Sethian. Evolution, implementation and application of level set and fast marching methods for advancing fronts. *J. Comput. Phys.*, 169:503–555, 2001.

7. G. Ventura, J.X. Xu, and T. Belytschko. A vector level set method and new discontinuity approximations for crack growth by efg. *Int. J. Numer. Meths. Eng.*, 54:923–944, 2002.

8. J. Yvonnet, D. Ryckelynck, P. Lorong, and F. Chinesta. The meshless constrained natural element method (C-NEM) for treating thermal involving moving interfaces. 2005.