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

HAL
open science

# Simulation of moving solid-phase boundaries using the partition of unity concept

Dominique Chamoret, Joris Remmers, Sergio Turteltaub, René de Borst

► To cite this version:

Dominique Chamoret, Joris Remmers, Sergio Turteltaub, René de Borst. Simulation of moving solid-phase boundaries using the partition of unity concept. 7e colloque national en calcul des structures, CSMA, May 2005, Giens, France. hal-01812896

HAL Id: hal-01812896

https://hal.science/hal-01812896

Submitted on 11 Jun 2018

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
# Simulation of moving solid-phase boundaries using the partition of unity concept

D. Chamoret<fnref n="1" /> — J.J.C. Remmers<fnref n="2" /> — S. Turteltaub<fnref n="2" /> — R. de Borst<fnref n="2" />

**ABSTRACT.** In this paper, numerical techniques are proposed to model phase transformation mechanisms. The evolution of interfaces is governed by field equations and jump conditions, and by a kinetic relation. Jumps in displacement gradients are incorporated by means of the partition of unity concept. The interface is described by a level set function.

**RÉSUMÉ.** Dans cet article, de nouvelles approches numériques pour modéliser les mécanismes de transformations de phases sont présentées. L'évolution de l'interface est prise en compte dans l'écriture des équations de champs associées à une relation cinématique. La discontinuité du champ de déformation à l'interface est modélisée par l'intermédiaire de la méthode de la partition de l'unité. Pour la description de l'interface, l'approche level set est mise œuvre.

**KEYWORDS:** Phase transformation mechanisms, moving boundaries, partition of unity, level set method.

**MOTS-CLÉS :** Transformation de phase, frontières mobiles, partition de l'unité, level set.

## 1. Introduction

Shape memory alloys constitute a special class of metallic alloys that have the ability to generate large work output per unit volume. They are particularly suited for applications in small devices such as in micro-electro-mechanical systems. Depending on the temperature, they have the possibility of recovering the initial shape of an object after a thermomechanical cycle (i.e., the shape memory effect) or upon a purely mechanical cycle (i.e., the pseudo-elasticity effect). These physical phenomena, which endow shape memory alloys of their unique characteristics, are related to crystallographically reversible, displacive phase transformations.

<footnote n="1">¹ Laboratoire M3M
Université de Technologie de Belfort Montbéliard, Site de Sévenans
90010 Belfort Cedex, France
dominique.chamoret@utbm.fr</footnote>

<footnote n="2">² Engineering Mechanics
Delft University of Technology, Faculty of Aerospace Engineering
Kluyverweg 1, 2629 HS Delft, The Netherlands
{J.J.C.Remmers, S.R.Turteltaub, R.deBorst}@lr.tudelft.nl</footnote>
Our goal is to obtain a comprehensive understanding of the phase transformation mechanisms at the length scale of individual interfaces between the parent and product phases, austenite and martensite, respectively. To this end, we propose a methodology suitable to analyze problems where the evolution of the interfaces throughout a transformation process is modelled explicitly. The numerical treatment of these problems raises several challenging issues, such as a mesh-independent description of sharp interfaces and the need of an efficient methodology to keep track of the interfaces' evolution. To address these requirements, we make use of the level set method [HOU 99] to monitor the evolution of interfaces (section 2) and we employ the partition of unity property [MEL 96] of finite element shape functions to model weak discontinuities (section 3).

## 2. Interface conditions and kinetic relation

### 2.1. Governing equations and jump conditions

Consider a deformable body $B$ that occupies the domain $\Omega$ in $\mathbb{R}^2$. Suppose that $\Omega$ contains phase 1 in region $\Omega^+$ and phase 2 in region $\Omega^-$, as shown in Figure 1(a) (e.g., austenite in $\Omega^+$ and martensite in $\Omega^-$). Interfaces between austenite and martensite in shape memory alloys are typically coherent, hence the displacement $\mathbf{u}$ is continuous across the interface but the displacement gradient $\nabla \mathbf{u}$ experiences a finite jump (weak discontinuity) [BHA 03]. These conditions can be shown to reduce to the following jump conditions:

$$ [[\nabla \mathbf{u}]] = \mathbf{c} \otimes \mathbf{n} \quad \text{on } \Gamma_t \qquad (1) $$

for some nonzero vector $\mathbf{c}$ and where $\mathbf{n}$ is the normal to the interface. Here, $[[w]] = w^+ - w^-$ is the jump of the function $w$. For a quasi-static process, where inertial effects are neglected, the balance of linear momentum at points where $\mathbf{u}$ is smooth takes the form

$$ \operatorname{div} \boldsymbol{\sigma} + \mathbf{f} = 0 \quad \text{in } \Omega - \Gamma_t \qquad (2) $$

where $\boldsymbol{\sigma}$ is the stress tensor and $\mathbf{f}$ is a body force. At points on the interface $\Gamma_t$, the balance of linear momentum reduces to the continuity of traction,

$$ [[\boldsymbol{\sigma}]] \mathbf{n} = 0 \quad \text{on } \Gamma_t. \qquad (3) $$

### 2.2. Constitutive model

Within the framework of small deformations, the total strain $\boldsymbol{\epsilon}$ is additively decomposed into an elastic strain $\boldsymbol{\epsilon}_e$ and a transformation strain $\boldsymbol{\epsilon}_{tr}$, i.e.,

$$ \boldsymbol{\epsilon} = \boldsymbol{\epsilon}_e + \boldsymbol{\epsilon}_{tr}. \qquad (4) $$

Neglecting rigid body rotations, the transformation strain is given as

$$ \boldsymbol{\epsilon}_{tr} = \mathbf{0} \quad \text{in } \Omega^+ \quad \text{and} \quad \boldsymbol{\epsilon}_{tr} = \frac{1}{2}(\mathbf{b} \otimes \mathbf{m} + \mathbf{m} \otimes \mathbf{b}) \quad \text{in } \Omega^- . \qquad (5) $$
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Description of the interface</figcaption></figure>

The shape strain vector **b** provides information about the unconstrained transformation shear and dilation in the martensite, measured with respect to the interface, and the vector **m** corresponds to the normal vector of an unstressed austenite/martensite interface. The vectors **b** and **m** can be obtained based on crystallographic information. Furthermore, the constitutive relation for each phase is given as

$$ \sigma = \mathbb{C}^A \epsilon_e \quad \text{in } \Omega^+ \quad \text{and} \quad \sigma = \mathbb{C}^M \epsilon_e \quad \text{in } \Omega^- \qquad (6) $$

where $\sigma$ is the stress tensor and $\mathbb{C}^A$ and $\mathbb{C}^M$ are the elasticity tensors of austenite and martensite, respectively. Since the elastic properties of the two phases are different, the elastic strains are discontinuous across an interface.

In order to characterize the evolution of the interface, it is necessary to give an additional piece of constitutive information, namely a kinetic relation that relates the driving force $f$ to the normal velocity $V_n$ of the interface [ABE 90]. The driving force for a slowly-propagating interface, across which the stress is continuous, is given as

$$ f = [[W]] - \sigma \cdot [[\epsilon]] \qquad (7) $$

where $W$ is the strain energy of the material, which, in view of (5) and (6), is given by

$$ W(\epsilon) = \frac{1}{2} \epsilon \cdot \mathbb{C}^A \epsilon \quad \text{in } \Omega^+ \quad \text{and} \quad W(\epsilon) = \frac{1}{2}(\epsilon - \epsilon_{tr}) \cdot \mathbb{C}^M (\epsilon - \epsilon_{tr}) \quad \text{in } \Omega^- . \quad (8) $$

The kinetic relation is formally given by a function $V$, defined for points on the interface, such that

$$ V_n = V \cdot n = V(f) \qquad (9) $$

where $V$ is the velocity of the interface and $n$ is the normal vector to the stressed interface (see Figure 1(a)). A specific form of kinetic relation (9) can be found in [HOU 99].
### 3. Description of the interface: level set

In order to represent the interface, we have chosen to use the level set method since it simplifies the mesh step and the geometrical data pre-processing [OSH 03]. The idea is to describe the interface as the zero level set of a smooth, scalar-valued function $\phi$. The level set function $\phi$ can be defined as the signed distance function:

$$ \Gamma_t = \{\boldsymbol{x} \in \Omega : \phi(\boldsymbol{x}, t) = 0\} \quad \text{with} \quad \phi(\boldsymbol{x}, t) = \pm \min_{\boldsymbol{x} \in \Gamma_t} \|\boldsymbol{x} - \tilde{\boldsymbol{x}}\|. \quad (10) $$

As shown in Figure 1(a), the function $\phi$ is chosen positive in $\Omega^+$ (phase 1), negative in $\Omega^-$ (phase 2). The interface motion is governed by a Hamilton-Jacobi evolution equation, i.e.,

$$ \dot{\phi}(\boldsymbol{x}, t) + V_n \|\nabla \phi(\boldsymbol{x}, t)\| = 0 \quad \text{in } \Omega. \quad (11) $$

In this expression $\dot{\phi}(\boldsymbol{x}, t)$ represents the total time derivative of $\phi$ and $V_n$ is the normal component of $\mathbf{V}$ the interface velocity.

#### 3.1. Weak formulation of the evolution equation

Equation (11) must be solved over the entire domain $\Omega$. However, the velocity $V_n$ is known only at the interface $\Gamma_t$. To circumvent this problem, the velocity field $V_n$ can be extended to the domain $\Omega$ by using an interpolation process [OSH 03]. To this end, we formally extend $V_n$ to the entire domain by introducing a function $\hat{V}_n$ defined as

$$ \hat{V}_n = V_n \delta(\Gamma_t), \quad (12) $$

where $\delta$ is the Dirac-delta distribution. In the weak formulation of the problem we transform the volume integral over the entire domain into a surface integral over the interface. So, the weak form of (11) is

$$ \int_{\Omega} \eta \dot{\phi}(\boldsymbol{x}, t) d\Omega + \int_{\Omega} \eta \hat{V}_n \|\nabla \phi(\boldsymbol{x}, t)\| d\Gamma = 0, \quad \forall \eta, \quad (13) $$

where $\eta$ is a test function. Equation (12) and the properties of the Dirac-delta distribution leads to the final form of the weak formulation, i.e.,

$$ \int_{\Omega} \eta \dot{\phi}(\boldsymbol{x}, t) d\Omega + \int_{\Gamma_t} \eta V_n \|\nabla \phi(\boldsymbol{x}, t)\| d\Gamma = 0, \quad \forall \eta. \quad (14) $$

#### 3.2. Space and time discretization

The level sets are approximated by finite elements $\phi(\boldsymbol{x}, t) = \mathbf{N}\Phi_t$ where $\mathbf{N}$ is matrix of the classical shape functions. Discrete expression of equation (14) can be written as

$$ \mathbb{M} \dot{\Phi}_t + f_{\phi_t} = 0 \quad \text{with} \quad \mathbb{M} = \int_{\Omega} \mathbf{N}^T \mathbf{N} \, d\Omega \quad \text{and} \quad f_{\phi_t} = \int_{\Gamma_t} \mathbf{N}^T V_n \|\nabla \phi_t\| \, d\Gamma. \quad (15) $$
The time-dependent equation (15) is discretized in time using the (explicit) forward-Euler integration scheme to obtain

$$ \mathbf{\Phi}_{t+\Delta t} = \mathbf{\Phi}_t - \Delta t \mathbf{M}^{-1} \mathbf{f}_{\phi_t} \quad (16) $$

where the subscript t denotes the time and Δt is the time step.

# 4. FEM, partition of unity and weak discontinuities

## 4.1. FEM and partition of unity

To describe weak discontinuities in the interior of elements, an “enhanced basis” is added to the classical finite element (polynomial) approximation by exploiting the partition of unity property [MEL 96]. Consequently, the approximation of the displacement field can be decomposed in two parts: the standard displacement field, composed of the classical shape functions $N_i$ and degrees of freedom $a_i$ connected to nodes i, and the enrichment part, composed of enrichment functions $N_j \psi_j$ and additional degrees of freedom $b_j$ introduced at enriched nodes j, i.e.,

$$ \mathbf{u}^h(\mathbf{x}) = \underbrace{\sum_{i \in \mathcal{I}} N_i(\mathbf{x}) a_i}_{\text{standard part}} + \underbrace{\sum_{j \in \mathcal{I}_e} N_j(\mathbf{x}) \psi_j(\mathbf{x}) b_j}_{\text{enriched part}} \quad (17) $$

Here, $\mathcal{I}$ represents the set of all nodes i while $\mathcal{I}_e$ is the set of enriched nodes j, which are chosen such that the support of $N_j$ is intersected by the interface $\Gamma_t$ (Figure 1(b)).

## 4.2. Weak discontinuities

In view of the strain decomposition (4), the discontinuity in the total strain across the interface $\Gamma_t$ is associated with (i) a discontinuity [[$\epsilon_e$]] in the elastic strain and (ii) a discontinuity [[$\epsilon_{tr}$]] in the transformation strain. From (5), the discontinuity on the transformation strain can be explicitly computed from the crystallography of the austenite and the martensite. The discontinuity on the elastic strain needs to be computed using a finite element formulation that utilizes enrichment functions to describe weak discontinuities. The derivatives of such functions must be discontinuous across the interface. A possible choice for $\psi_j$ is to take the absolute value of the level set [BEL 01], although other possibilities can be found in [MO03]. This function has indeed discontinuous first derivatives at the interface. The approximation (17) can be written as

$$ \mathbf{u}^h(\mathbf{x}) = \tilde{\mathbf{u}}(\mathbf{x}) + |\phi(\mathbf{x})|\hat{\mathbf{u}}(\mathbf{x}), \quad (18) $$

where, $\tilde{\mathbf{u}}(\mathbf{x}) = \sum_{i \in \mathcal{I}} N_i(\mathbf{x}) a_i$ and $\hat{\mathbf{u}}(\mathbf{x}) = \sum_{j \in \mathcal{I}_e} N_j(\mathbf{x}) b_j$. The gradient of this expression is clearly discontinuous and is given by

$$ \nabla \mathbf{u}^h(\mathbf{x}) = \nabla \tilde{\mathbf{u}}(\mathbf{x}) + |\phi(\mathbf{x})| \nabla \hat{\mathbf{u}}(\mathbf{x}) + \text{sign}(\phi) \hat{\mathbf{u}}(\mathbf{x}) \otimes \nabla \phi(\mathbf{x}) \quad (19) $$

Therefore the jump gradient at the interface is

$$ [[\nabla \mathbf{u}^h(\mathbf{x})]] = 2\hat{\mathbf{u}}(\mathbf{x}) \otimes \nabla \phi(\mathbf{x}) = 2\|\nabla \phi(\mathbf{x}, t)\| \hat{\mathbf{u}}(\mathbf{x}) \otimes \mathbf{n}. \quad (20) $$
This relation introduces the jump in the gradient of the displacement as specified by the jump Hadamard condition (1).

### 4.3. Discrete expression

The enriched finite element approximation for the enriched displacement $u^h$ and the enriched strain $\epsilon^h$ in an element $\Omega_e$ are defined as follows:

$$u^h = N_e a_e \quad \text{and} \quad \epsilon^h = B_e a_e \qquad (21)$$

where $a_e$ is the vector of degrees of freedom of the element. $N_e$ is the matrix of the shape functions (classical and enriched) and $B_e$ is the matrix of the derivatives of the shape functions. We then obtain the following finite element expression of the problem:

$$K a = f_{ext} + f_{tr} \qquad (22)$$

with $K = \int_{\Omega^+} \mathbf{B}^T \mathcal{C}^A \mathbf{B} d\Omega + \int_{\Omega^-} \mathbf{B}^T \mathcal{C}^M \mathbf{B} d\Omega$; $f_{ext} = \int_{\Omega} \mathbf{N}^T \mathbf{f} d\Omega + \int_{\Gamma_p} \mathbf{N}^T \mathbf{t} d\Gamma$

and $f_{tr} = \int_{\Omega^-} \mathbf{B}^T \mathcal{C}^M \epsilon_{tr} d\Gamma$.

The use of the present algorithm will be illustrated in a forthcoming publication.

### 5. References

[ABE 90] ABEYARATNE R., KNOWLES J., “On the driving traction acting on a surface of strain discontinuity in a continuum”, *Journal of the Mechanics and Physics of Solids*, vol. 38, num. 3, 1990, p. 345–360.

[BEL 01] BELYTSCHKO T., MOËS N., USUI S., PARIMI C., “Arbitrary discontinuities in finite element”, *International Journal for Numerical Methods in Engineering*, vol. 50, 2001, p. 993–1013.

[BHA 03] BHATTACHARYA K., *Microstructure of martensite. Why it forms and how it gives rise to the shape-memory effect*, Oxford Series on Materials Modelling, Oxford University Press, 2003.

[HOU 99] HOU T., ROSAKIS P., LEFLOCH P., “A Level-Set Approach to the Computation of Twinning and Phase-Transition Dynamics”, *Journal of Computational Physics*, vol. 150, num. 2, 1999, p. 302–331.

[MEL 96] MELENK J., BABUSKA I., “The partition of unity finite element method: Basic theory and applications”, *Computer Methods in Applied Mechanics and Engineering*, vol. 139, num. 1-4, 1996, p. 289–314.

[MOÕ3] MOËS N., CLOIREC M., CARTRAUD P., REMACLE J., “A computational approach to handle complex microstructure geometries”, *Computer Methods in Applied Mechanics and Engineering*, vol. 192, 2003, p. 3163-3177.

[OSH 03] OSHER S., PARAGIOS N., *Geometric Level Set Methods in Imaging, Vision, and Graphics*, Springer-Verlag Telos, 2003.