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

**HAL**
open science

# Simulating Dynamic Thermo-Elasto-Plasticity in large Transformations with Adaptive Refinement in the NEM. Application to Shear Banding

Julien Yvonnet, Philippe Lorong, David Ryckelynck, Francisco Chinesta

► To cite this version:

Julien Yvonnet, Philippe Lorong, David Ryckelynck, Francisco Chinesta. Simulating Dynamic Thermo-Elasto-Plasticity in large Transformations with Adaptive Refinement in the NEM. Application to Shear Banding. International Journal of Forming Processes, 2005, 8. hal-01508714

HAL Id: hal-01508714

https://hal.science/hal-01508714

Submitted on 25 Apr 2017

**HAL** is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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
Simulating dynamic thermo-elasto-plasticity
in large transformations with adaptive refine-
ment in the natural element method : appli-
cation to shear banding

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

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

ABSTRACT. Simulation of forming processes usually involves large transformations and/or stain localization. In the context of the finite element method (FEM), a common approach to deal with such phenomena is the use of continuous adaptive remeshing. Nevertheless, constructing a quality mesh is a delicate task, especially for complex 3D problems. The natural element meshfree method is a recent numerical technique which uses the features of the Delaunay triangulation of the set of nodes, providing an accuracy equivalent to the quadrilateral/hexahedral finite elements, even if the Delaunay triangles are very distorted. In this context, inserting or removing nodes is an easy task, as no special geometrical criterion is imposed on the relative position of the nodes. Furthermore, the use of a stabilized conforming nodal integration allows to define nodal internal variables, which simplifies transfer of these data through successive updates of the reference configuration in a Lagrangian procedure. In this paper, a framework for practical implementation of the NEM in the context of explicit thermo-elastoplasticity at finite strains is provided. An adaptive strategy, based on the attractive features of the NEM is developed, including simple error indicators for history-dependent nonlinear problems, a refinement procedure for unstructured scattered of nodes based on the underlying Voronoi diagram, and the use of a nodal integration to alleviate issues associated with internal variables transfer. This promising technique for forming processes simulation is illustrated in the context of dynamic shear bands propagation after impact.

KEYWORDS: Natural element method, dynamic thermo-elastoplasticity, large transformations, adaptive refinement, shear bands
# 1. Introduction

Simulation of forming processes usually involves large transformations and/or strain localization. In the context of the finite element method (FEM), it is necessary to construct a mesh satisfying some quality criteria for a given deformed geometry and which also conforms the size of its elements to the localized phenomena, i.e. strain localization or cracks. Triangular and tetrahedral meshes can be automatically generated based on the Delaunay criterion in the FEM whatever the complexity of geometry shapes of the analyzed domain. Nevertheless, there are no guarantee that no flat elements can be generated, and the triangular and tetrahedral elements are rarely adopted due to their poor accuracy. Quadrilateral/hexahedral elements are usually preferred for its better accuracy, resulting in the great difficulty of mesh generation, especially for complex geometries and in the context of adaptive refinement. In recent years, meshless approximations have become interesting and promising methods in solving partial differential equations due to their flexibility in practical applications. Several meshless methods, such as, the smooth particle hydrodynamics (SPH) [24], the diffuse element method (DEM) [27], the element-free Galerkin method (EFG) [6], the reproducing kernel particle method (RKPM) [22], the partition of unity (PUM) [2], and the natural element method (NEM) [39] have been proposed. The NEM offers outstanding properties for the simulation of problems involving large transformations and adaptive refinements due to the following properties: (a) essential boundary can be enforced directly as the NEM shape functions satisfies the Kronecker delta property, and are also strictly linear over the boundaries of the convex hull (it has been extended to non convex geometries in [46]); (b) the support of the shape functions, constructed on the basis of the underlying Voronoi diagram (dual of the Delaunay triangulation), automatically adapts its size and shape to the local neighborhood, whatever the complexity and density of the surrounding nodal distribution; (c) the accuracy does not depend significantly on the relative position of the nodes (shape of the Delaunay triangles), and it results equivalent to quadrangular/hexahedral finite elements [42] even when Delaunay triangles become very distorted. Is noteworthy that this property greatly simplifies the treatment of complex three-dimensional domains, as the Delaunay triangulation is unique for a given set of nodes.

In this paper, we analyze the application of the C-NEM to the simulation of dynamic thermo-elastoplastic behaviors, involving large transformations combined with an adaptive strategy based on the attractive features of the NEM. In section 2, the natural element method is briefly reviewed, as well as its extension for arbitrary (non convex) geometries that results in the constrained natural element method (C-NEM) [46]. In section 3, a Lagrangian formulation of the coupled thermo-mechanical problem in large transformations is established. In section 4, a strategy for refinement and adaptivity in the NEM (or C-NEM) context is proposed, based on the attractive features of the NEM. A simple error indicator using the Zienkiewicz-Zhu [47] paradigm in tandem with the stabilized conforming nodal integration is proposed. A refinement strategy based on the underlying Voronoi diagram is proposed, which allows refinement through arbitrary unstructured clouds of nodes. A nodal-NEM formulation [43]
<figure><img src="image_2.png" /><figcaption>Figure 1. Construction of the Sibson shape functions.</figcaption></figure>

is also used to define the internal variables, usually defined at the integration points in the context of FEM, at the nodes, which greatly simplifies the transfer of these variables through successive updates of the reference configuration. Finally, the technique is applied in section 5 in the context of dynamic shear bands propagation in metals.

## 2. The constrained natural element method

### 2.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 Braun and Sambridge [31], and Sukumar et al. [39]. 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}^{\text{dim}}$. The Voronoi diagram is the subdivision of $\mathbb{R}^{\text{dim}}$ into regions $T_i$ (Voronoi cells) defined by:

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

The Sibson coordinates of $\mathbf{x}$ with respect to a natural neighbor $n_i$ (see Fig. 1) 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}(afghe)}{\text{Area}(abcde)} \quad (2) $$
If the point x coincides with the node $n_i$, i.e. ($x = x_i$), $\phi_i(x_i) = 1$, and all other shape functions are zero, i.e. $\phi_i(x_j) = \delta_{ij}$ ($\delta_{ij}$ being the Kronecker delta). The properties of positivity, interpolation, and partition of unity are then verified [39]:

$$ \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. \qquad (3) $$

The natural neighbor shape functions also satisfy the local coordinate property [40], namely:

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

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

It turns out that the support of $\phi_i(\mathbf{x})$ is the union of the $n$ circles (spheres in 3D) passing through the vertices of the $n$ Delaunay triangles (tetrahedra) containing the node $n_i$ (in this case $n$ is the number of natural neighbors of node $n_i$). The support of a node $n_i$ in a particular nodal distribution is depicted in figure 2.

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 [14] have shown that the Sibson interpolant belongs to a more general class of Voronoi-based interpolants, called k-th order standard coordinates, proving that the interpolant generated by the k-th order standard coordinates have $C^k$ continuity on the Delaunay circles (spheres) boundaries. In this context, the Sibsonian and non-Sibsonian (Laplace) coordinates [4] 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. [39]. An illustration is depicted in Fig. 1 (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 function (vectorial, scalar or tensorial) $\mathbf{u}(\mathbf{x})$, in the form:

$$ \mathbf{u}^h(\mathbf{x}) = \sum_{i=1}^{n} \phi_i(\mathbf{x}) \mathbf{u}_i \qquad (5) $$
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Support of the natural shape function related to node $n_i$.</figcaption></figure>

where $u_i$ are the nodal values 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. (5) defines a local interpolation scheme. Thus, the trial and test functions used in the discretization of the variational formulation describing both the mechanical and thermal problems treated in this paper take the form of Eq. (5).

## 2.2. The Constrained natural element method

### 2.2.1. Constrained Voronoi diagram

In its original form [39], the NEM can only be applied to strictly convex domains. For non-convex domains, two main issues occur: (a) It was proved in [39, 46] and [9] that a loss of linearity in the interpolation along boundaries of non convex domains appear. Thus essential boundary conditions can only be imposed directly over convex boundaries; (b) For strongly non-convex domains (cracks, auto-contact...) some spurious influences between nodes of the boundaries appear [46]. An additional treatment is thus required to maintain all the properties of the NEM for any geometry.

In order to avoid these drawbacks, we have proposed in a previous paper [46] an extension of the NEM in which a visibility criterion is 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), which is the strict dual to the constrained Delaunay triangulation, introduced by Seidel in [35], instead of the Voronoi diagram (see [46] for further details). The intersection between the CVD and the domain results into new cells $T_i^C$, called constrained Voronoi cells, 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_{x \to n_i} \cap \Gamma = \emptyset, S_{x \to n_j} \cap \Gamma = \emptyset \} \quad (6) $$

where $\Gamma$ is the domain boundary, composed by a set of segments $l_i \in L$ and $S_{a \to b}$ denotes the segment between the points a and b. In this framework, a point located inside a cell $T_i^C$ is closer to the node $n_i$ than to any other visible node $n_j$.
The constrained Delaunay triangulation does not always exist in 3D without adding new nodes, as shown in [32]. Nevertheless, some techniques for constructing 3D constrained Delaunay tessellations are available and provided in [33, 34] by adding Steiner points.

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

$$ \mathbf{u}^h(\mathbf{x}) = \sum_{i=1}^{V} \phi_i^C(\mathbf{x})\mathbf{u}_i \quad (7) $$

where $V$ is the number of natural neighbors visible from point $\mathbf{x}$ and $\phi_i^C$ is the constrained natural neighbor shape function related to the $i$-th node at point $\mathbf{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 [46] and [45] 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 to any geometry, convex or not.

The ability of the C-NEM for treating problems involving cracks has been illustrated in [46] and in the context of moving interfaces in [45], where we have shown how the C-NEM simplifies the treatment of material discontinuities in meshfree methods, due to the continuity of the approximation across interfaces (consequence from the interpolant character and linearity of the C-NEM shape functions over any external boundary or internal interface). It is noteworthy that no size parameter is involved in the definition of the influent nodes (neighbors), which is essential for the robustness of refinement procedures, in which the nodal density varies in some parts of the domain.

## 3. Formulation of the coupled thermo-mechanical problem

### 3.1. Preliminaries

In this section, we summarize briefly some fundamental kinematic relations and introduce the appropriate notations. We consider two configurations of a body $B$: the first one, the reference configuration $\Omega_0 \subset \mathbb{R}^{dim}$ with external boundary $\partial\Omega_0$ (not necessarily the initial configuration), where $\mathbf{X}$ denotes the coordinates of a point in this configuration. The second, called current configuration, is denoted by $\Omega_x \subset \mathbb{R}^{dim}$ with external boundary $\partial\Omega_x$ at time $t$, with $\mathbf{x}$ the coordinates of a point in the current configuration. $\mathbf{x}$ is related to $\mathbf{X}$ by:

$$ \mathbf{x} = \mathbf{X} + \mathbf{u}(\mathbf{X}, t) \quad (8) $$
the deformation gradient is defined by:

$$ \mathbf{F} = \frac{\partial \mathbf{x}}{\partial \mathbf{X}} = \mathbf{1} + \nabla_{\mathbf{x}} \mathbf{u}, \quad J = \det(\mathbf{F}) > 0 \qquad (9) $$

From now on, we use a hyperelastic formulation based on the multiplicative decomposition of the deformation gradient, with elastic response described by a hyper-elastic stored energy function. In this context, elastic predictor becomes exact and the need for incrementally objective algorithms is entirely avoided.

The basic hypothesis underlying this approach to finite strain elastoplasticity is the multiplicative split of the deformation gradient, F, into elastic and plastic parts:

$$ \mathbf{F} = \mathbf{F}^e \mathbf{F}^p \qquad (10) $$

This assumption, firstly proposed by Lee [18], admits the existence of a local unstressed intermediate configuration. Following the multiplicative split of F, the velocity gradient, L = ḞF⁻¹ can be decomposed additively as

$$ \mathbf{L} = \mathbf{L}^e + \mathbf{L}^p \qquad (11) $$

where L^e and L^p are, respectively, the elastic and plastic contributions defined by:

$$ \mathbf{L}^e = \dot{\mathbf{F}}^e [\mathbf{F}^e]^{-1}, \quad \mathbf{L}^p = \mathbf{F}^e \dot{\mathbf{F}}^p [\mathbf{F}^p]^{-1} [\mathbf{F}^e]^{-1} \qquad (12) $$

Similarly, the stretching tensor, D := sym[L], can be decomposed as:

$$ \mathbf{D} = \mathbf{D}^e + \mathbf{D}^p \qquad (13) $$

with the elastic and plastic stretching tensors given by

$$ \mathbf{D}^e = \mathrm{sym}[\mathbf{L}^e], \quad \mathbf{D}^p = \mathrm{sym}[\mathbf{L}^p] \qquad (14) $$

## 3.2. Hyperelastic constitutive law

From polar decomposition, F^e is given by

$$ \mathbf{F}^e = \mathbf{R}^e \mathbf{U}^e \qquad (15) $$

where U^e and R^e are, respectively, the elastic right stretch tensor and the elastic rotation.

Let ε^e denote the Eulerian (or spatial) logarithmic strain tensor [3]

$$ \epsilon^e = \ln[\mathbf{U}^e] \qquad (16) $$
where $ln[\cdot]$ above denotes the tensor logarithm of (.) which involves spectral decomposition of $U^e$. Following Periò et al. [30], we assume the existence of a quadratic strain energy function $\psi^e(\epsilon^e)$ in the form of a scalar symmetric function of its stretches $\lambda_i$ ($i = 1, 2, 3$) given by:

$$ \psi^e(\lambda_1^e, \lambda_2^e, \lambda_3^e) = \mu [\ln(\lambda_1^e)^2 + \ln(\lambda_2^e)^2 + \ln(\lambda_3^e)^2] + \frac{1}{2}\lambda(J^e)^2 \quad (17) $$

where $\mu$ and $\lambda$ are Lamé's parameters and $(J^e) = \lambda_1^e \lambda_2^e \lambda_3^e$ is the Jacobian. After applying standard procedure, the following hyperelastic constitutive equation is obtained:

$$ T = \frac{\partial \psi^e}{\partial \epsilon^e} = C^e : \epsilon^e \quad (18) $$

where $T$ is the *rotated* stress tensor. Assuming incompressibility of the plastic flow, it is expressed as:

$$ T = [R^e]^T \tau R^e \quad (19) $$

where $\tau = J\sigma$ is the *Kirchhoff* stress tensor, $C^e$ is the fourth-order isotropic elastic tensor. Further details about the thermo mechanical foundations can be found in [12].

### 3.3. Numerical integration of the constitutive equations

Typically, within an incremental numerical procedure for solving history dependent problems, a numerical approximation to the material constitutive law is needed to update stresses $\tau$ as well as the internal variables $\alpha \equiv \bar{\epsilon}^p$ (equivalent plastic strain) within each time (load) increment. In the present context, given the values of the variables $\{\tau_n, F_n^p, \alpha_n\}$ at the beginning of a generic increment $[t_n, t_{n+1}]$, an algorithm is required to update $\{\tau_{n+1}, F_{n+1}^p, \alpha_{n+1}\}$ at the end of the increment.

Under the assumption of elastic isotropy, the elastic Eulerian logarithmic strains tensor is updated according to:

$$ \epsilon_{n+1}^e = \epsilon_{n+1}^{\text{trial}} - \Delta\gamma \mathbf{N}_{n+1} \quad (20) $$

which has the same format that the standard return mapping used in the infinitesimal theory [12, 36, 11, 38].

$$ \mathbf{N}_{n+1} = \frac{3 \mathbf{dev}(\mathbf{T}_{n+1})}{2 J_2(\mathbf{T}_{n+1})} \quad (21) $$

where $J_2(.)$ is the second invariant of (.) The trial elastic logarithmic strain $\epsilon_{n+1}^{\text{trial}}$ is given by:

$$ \epsilon_{n+1}^{\text{trial}} = \ln[\mathbf{U}_{n+1}^{\text{trial}}] \quad (22) $$
where $\mathbf{U}_{n+1}^{\text{trial}}$ results from polar decomposition of $\mathbf{F}_{n+1}^{\text{trial}}$ defined by:

$$ \mathbf{F}_{n+1}^{\text{trial}} = \mathbf{F}_{n+1}(\mathbf{F}_n^p)^{-1} \qquad (23) $$

From Eq. (18) we have:

$$ \mathbf{dev}(\mathbf{T}_{n+1}) = 2\mu\mathbf{dev}(\boldsymbol{\epsilon}_{n+1}^{e}) \qquad (24) $$

Using Eq. (20), Eq. (21) and Eq. (24) we obtain:

$$ \mathbf{dev}(\mathbf{T}_{n+1}) = \mathbf{dev}(\mathbf{T}_{n+1}^{\text{trial}}) - 3\mu\Delta\gamma \frac{\mathbf{dev}(\mathbf{T}_{n+1})}{J_2(\mathbf{T}_{n+1})} \qquad (25) $$

with:

$$ \mathbf{T}_{n+1}^{\text{trial}} = \mathbf{C}^{e} : \boldsymbol{\epsilon}_{n+1}^{e \text{ trial}} \qquad (26) $$

which leads to:

$$ J_2(\mathbf{T}_{n+1}) = J_2(\mathbf{T}_{n+1}^{\text{trial}}) - 3\mu\Delta\gamma \qquad (27) $$

We assume that $f_{y_{n+1}}^{\text{trial}} = J_2(\mathbf{T}_{n+1}^{\text{trial}}) - \sigma_y(\Delta\gamma)$ is a non-linear scalar function. In the case of a plastic increment ($f_{y_{n+1}}^{\text{trial}} > 0$), we must solve for $\Delta\gamma$:

$$ J_2(\mathbf{T}_{n+1}^{\text{trial}}) - 3\mu\Delta\gamma - \sigma_y(\Delta\gamma) = 0 \leftrightarrow J_2(\mathbf{T}_{n+1}) - \sigma_y(\Delta\gamma) = 0 \qquad (28) $$

A classical Newton-Raphson procedure has been used in this work to solve the above equation.

Using Eq. (25) we obtain:

$$ \mathbf{dev}(\mathbf{T}_{n+1}) = \frac{\mathbf{dev}(\mathbf{T}_{n+1}^{\text{trial}})}{1 + \frac{3\mu\Delta\gamma}{J_2(\mathbf{T}_{n+1})}} \qquad (29) $$

with $J_2(\mathbf{T}_{n+1})$ defined according to Eq. (27) and:

$$ \mathbf{T}_{n+1} = \mathbf{dev}(\mathbf{T}_{n+1}) + \frac{1}{3} \mathrm{Tr}(\mathbf{T}_{n+1}^{\text{trial}}) \qquad (30) $$

The Cauchy stress tensor may thus be obtained as:

$$ \boldsymbol{\sigma}_{n+1} = \frac{1}{\det(\mathbf{F}_{n+1})} \left[ \mathbf{R}_{n+1}^{\text{trial}} \right]^{-T} \mathbf{T}_{n+1} \left[ \mathbf{R}_{n+1}^{\text{trial}} \right]^{-1} \qquad (31) $$

Finally using Eq. (20) and Eq. (22), and assuming that:

$$ \mathbf{R}_{n+1}^{e} = \mathbf{R}_{n+1}^{e \text{ trial}} \qquad (32) $$
Box 1. Algorithm for integration of constitutive equations.

(i) For given displacement $\mathbf{u}^{n+1}$, evaluate total deformation gradient

$$ \mathbf{F}^{n+1} = \mathbf{1} + \nabla_X \mathbf{u}^{n+1} = \mathbf{1} + \mathbf{B}(\mathbf{X})\mathbf{u}^{n+1} $$

with $\mathbf{B}(\mathbf{X})$ a matrix containing the shape functions derivatives in the reference configuration.

(ii) Evaluate elastic trial deformation gradient

$$ (\mathbf{F}^e)_{n+1}^{trial} = (\mathbf{F}_{n+1})(\mathbf{F}_n^p)^{-1} $$

(iii) Perform polar decomposition of $(\mathbf{F}^e)_{n+1}^{trial}$

$$ (\mathbf{F}^e)_{n+1}^{trial} = \mathbf{R}_{n+1}^e \mathbf{U}_{n+1}^{trial} $$

(iv) Evaluate elastic trial logarithmic strain tensor

$$ (\boldsymbol{\epsilon}^e)_{n+1}^{trial} = \ln[\mathbf{U}_{n+1}^{trial}] $$

(v) Evaluate trial stress tensor :

$$ \mathbf{T}_{n+1}^{trial} := \mathbf{C}^e : (\boldsymbol{\epsilon}_{n+1}^e)^{trial} $$

(vi) Check plastic consistency condition

IF $J_2(\mathbf{T}_{n+1}^{trial}) - \sigma_y^n \le 0$ THEN :

(.)$_{n+1} = (.)_{n+1}^{trial}$ and RETURN

ELSE go to (vii)

(vii) Plastic corrector (solve for $\Delta\gamma$)

$$ J_2(\mathbf{T}_{n+1}) - \sigma_y^{n+1}(\Delta\gamma) = 0, \text{ with } J_2(\mathbf{T}_{n+1}) = J_2(\mathbf{T}_{n+1}^{trial}) - 3\mu\Delta\gamma $$

(viii) Update Cauchy stress tensor

$$ \boldsymbol{\sigma}_{n+1} = J_{n+1}^{-1} [\mathbf{R}_{n+1}^e]^{-T} \mathbf{T}_{n+1} \mathbf{R}_{n+1}^{e-1}, \text{ with } \mathbf{T}_{n+1} = \mathbf{dev}(\mathbf{T}_{n+1}) + \frac{1}{3} \mathrm{Tr}(\mathbf{T}_{n+1}^{trial}) \mathbf{1} $$

$$ \mathbf{dev}[\mathbf{T}_{n+1}] = \frac{\mathbf{dev}[\mathbf{T}_{n+1}^{trial}]}{1+\Xi}, \Xi = \frac{3\mu\Delta\gamma}{J_2(\mathbf{T}_{n+1})} $$

(ix) Update plastic part of deformation gradient

$$ \mathbf{F}_{n+1}^p = \exp[\Delta\gamma\mathbf{N}_{n+1}] \mathbf{F}_n^p, \quad \mathbf{N}_{n+1} = \frac{3}{2} \frac{\mathbf{dev}(\mathbf{T}_{n+1})}{J_2(\mathbf{T}_{n+1})} $$

we obtain, after some calculations, the following incremental law for $\mathbf{F}^p$ (see [12, 36, 11]):

$$ \mathbf{F}_{n+1}^p = \exp[\Delta\gamma\mathbf{N}_{n+1}] \mathbf{F}_n^p \quad (33) $$

As noticed in [30], as a consequence of the exponential mapping in the implicit integration of the plastic flow rule, the incompressibility of the plastic flow for pressure insensitive yield criteria is carried over exactly to the incremental rule (33). The algorithm, therefore, generalizes the standard *return mapping algorithms* [37] of the infinitesimal theory. The overall algorithm for the incremental stress update is outlined in Box 1.
### 3.4. Explicit Lagrangian procedure

With the principle of virtual work as a basis of kinematically based C-NEM solution scheme, the corresponding continuum incremental boundary value problem is formulated in the spatial configuration as follows.

$$ \int_{\Omega^t} \rho(t) \ddot{\mathbf{u}} \cdot \boldsymbol{\eta} d\Omega^t + \int_{\Omega^t} \boldsymbol{\sigma}^t : \nabla_{\mathbf{x}} \boldsymbol{\eta} d\Omega^t = \int_{\Omega^t} \rho(t) \mathbf{b} \cdot \boldsymbol{\eta} d\Omega^t + \int_{\partial \Omega^t_{\sigma}} \mathbf{t} \cdot \boldsymbol{\eta} d\Gamma^t \quad \forall \eta \in \vartheta \quad (34) $$

where $\vartheta$ is the space of virtual displacements. The properties $d\Omega^t = J^t d\Omega_0$ and $\rho_0 d\Omega_0 = \rho(t) d\Omega^t$ are used, which leads to:

$$ \int_{\Omega_0} \rho_0 \ddot{\mathbf{u}} \cdot \boldsymbol{\eta} d\Omega_0 + \int_{\Omega_0} \mathbf{P}^t : \nabla_{\mathbf{x}} \boldsymbol{\eta} d\Omega_0 = \int_{\Omega_0} \rho_0 \mathbf{b} \cdot \boldsymbol{\eta} d\Omega_0 + \int_{\partial \Omega^t_{\sigma}} \mathbf{t} \cdot \boldsymbol{\eta} d\Gamma^t \quad \forall \eta \in \vartheta \quad (35) $$

where $\mathbf{P}$ denotes the first Piola-Kirschhoff stress tensor related to $\boldsymbol{\sigma}$ by $\mathbf{P} = J\mathbf{F}^{-1}\boldsymbol{\sigma}$.

The C-NEM discretization (7) of the variational form (35) results in the discrete set of algebraic time dependent equations which may be expressed, in matrix form, as:

$$ \mathbf{M}\ddot{\mathbf{u}}_{n+1}(t) = \mathbf{F}_{n}^{\text{ext}}(t) - \mathbf{F}_{n}^{\text{int}}(\mathbf{u}_n, t) \quad (36) $$

where $t$ is the time, $\mathbf{M}$ denotes the mass matrix, $\mathbf{F}_{n}^{\text{int}}(\mathbf{u}, t)$ the internal force vector, while $\mathbf{F}_{n}^{\text{ext}}(t)$ is the external force vector, expressed, respectively, by:

$$ \mathbf{M} = \int_{\Omega_0} \rho_0 \phi^T(\mathbf{X}) \phi(\mathbf{X}) d\Omega_0 \quad (37) $$

$$ \mathbf{F}_{n}^{\text{int}}(\mathbf{u}_n, t) = \int_{\Omega_0} J_n \mathbf{F}_{n}^{-1} \sigma_n \mathbf{B}(\mathbf{X}) d\Omega_0 \quad (38) $$

$$ \mathbf{F}_{n}^{\text{ext}} = \int_{\partial \Omega_{\sigma}^{t}} \phi^{T}(\mathbf{x}) \mathbf{t} d\Gamma^{t} \quad (39) $$

with $\phi(\mathbf{X})$ a matrix containing the shape functions in the reference configuration and $\mathbf{B}(\mathbf{X})$ a matrix containing the shape functions derivatives also in the reference configuration. As shown in the next section, the use of the SCNI quadrature [10] results in a $\mathbf{M}$ matrix diagonal, whose diagonal terms are given by $m_i = \rho_0 \Omega_i$, with $\Omega_i$ the area (volume in 3D) of the Voronoi cell related to node $n_i$.

The velocity $\mathbf{v} = \dot{\mathbf{u}}$ and acceleration $\ddot{\mathbf{u}} = \dot{\mathbf{v}}$ are approximated by using central differences with variable time steps. Thus, we have:
$$ \mathbf{v}_{n+1/2} = \mathbf{v}_{n-1/2} + \frac{\Delta t_1 + \Delta t_2}{2} \ddot{\mathbf{u}}_n \quad (40) $$

$$ \mathbf{u}_{n+1} = \mathbf{u}_n + \Delta t_2 \mathbf{v}_{n+1/2} \quad (41) $$

Finally, the general explicit algorithm is outlined as follows. Being known the initial conditions or the computed solution at time: $t_n$: $\mathbf{u}_n$, $\mathbf{v}_{n-1/2}$, $\ddot{\mathbf{u}}_n$, $\mathbf{F}_n$, $\mathbf{F}_n^p$:

(i) Update displacements and velocity

$$ \mathbf{v}_{n+1/2} = \mathbf{v}_{n-1/2} + \frac{\Delta t_1 + \Delta t_2}{2} \ddot{\mathbf{u}}_n $$

$$ \mathbf{u}_{n+1} = \mathbf{u}_n + \Delta t_2 \mathbf{v}_{n+1/2} $$

(ii) Evaluate $\sigma_{n+1}$, $\mathbf{F}_{n+1}^p$ by using Box 1.

(iii) Update accelerations $\ddot{\mathbf{u}}_{n+1}$

$$ (\ddot{u}_i)_{n+1} = \frac{1}{m_i} \left\{ (F_i^{ext})_{n+1} - [F_i^{int}]_{n+1} \right\} $$

Remark: In case of contact, $\ddot{u}_{n+1}$ is only used as a predictor phase, which has to be corrected according to the prescribed displacements or traction.

### 3.5. Thermo-mechanical coupling

The weak form of the heat balance can be expressed as

$$ \int_{\Omega^t} \rho(t)c(t)\dot{T}\eta d\Omega^t + \int_{\Omega^t} k(t)\nabla_x T \cdot \nabla_x \eta d\Omega^t = \int_{\Omega^t} r^t \eta d\Omega^t + \int_{\partial_2 \Omega} \bar{q}\eta d\Gamma^t \quad \forall \eta \in V_T \quad (42) $$

where $c(t)$ is the specific heat, $k(t)$ is the thermal conductivity for isotropic conduction, and $r^t$ a heat source related to the inelastic deformations, defined at time $t$. In the following, we assume $c$ and $k$ constant in time. $V_T$ is the space of virtual temperatures and $\bar{q}$ represents the heat transfer at the boundary $\partial_2 \Omega$. Following similar arguments as in the mechanical problem (34), the weak form of the heat balance can be expressed in the reference configuration by

$$ \int_{\Omega_0} \rho_0 c \dot{T} \eta d\Omega_0 + \int_{\Omega_0} k \nabla_x T F^{-T} F^{-1} \nabla_x \eta d\Omega_0 = \int_{\Omega_0} \det(F) r^t \eta d\Omega_0 + \\ + \int_{\partial_2 \Omega} \bar{q} \eta d\Gamma^t \quad \forall \eta \in V_T \quad (43) $$

The C-NEM discretization of Eq. (43) results in the ODE system:

$$ C\dot{T} + KT = Q \quad (44) $$
which in the context of an explicit scheme can be written as

$$CT^{n+1} = [C + \Delta t K] T^n + \Delta t Q^n$$

with

$$C = \int_{\Omega_0} \rho_0 c \phi^T(\mathbf{X}) \phi(\mathbf{X}) d\Omega_0 \quad (45)$$

$$K = \int_{\Omega_0} k B^T(\mathbf{X}) F^{-T} F^{-1} B(\mathbf{X}) d\Omega_0 \quad (46)$$

$$Q^n = \int_{\partial \Omega_2^t} \phi^T(\mathbf{x}) \bar{q} d\Gamma^t \quad (47)$$

The matrix C becomes diagonal in the context of C-NEM when stabilized conforming nodal integration is used, being the diagonal terms $c_i = \rho_0 c \Omega_i$. The heat source resultig from the inelastic deformations is given by

$$r^t = \chi \sigma^t : \hat{D}_t^p \quad (48)$$

with

$$\hat{D}_t^p = \text{sym} \left[ \left( \dot{\mathbf{F}}_t^p \right) \left( \mathbf{F}_t^p \right)^{-1} \right] \quad (49)$$

where $\sigma^t$ is the Cauchy stress tensor et time $t$, and $\chi$ is the Taylor-Quinney parameter [41] representing the fraction of plastic work converted into heat. In the present work, we have used $\chi = 0.9$ [26]. The thermo-mechanical coupling is carried out by the effects of inelastic deformations (48) and by the softening effects due to temperature in the hardening law. The hardening law adopted in the present simulation is given by the classical Johnson-Cook law [15]

$$\sigma_y(\bar{\epsilon}^p) = [A + B(\bar{\epsilon}^p)^n] \left[1 + C \ln\left(\frac{\dot{\bar{\epsilon}}^p}{\bar{\epsilon}_0^p}\right)\right] \left[1 - \left(\frac{T-T_0}{T_f-T_0}\right)^m\right] \quad (50)$$

where A, B, C, m, n and $\dot{\bar{\epsilon}}_0^p$ are material parameters, $\bar{\epsilon}^p$ represents the equivalent plastic strain and $\dot{\bar{\epsilon}}^p$ represents the rate of plastic strain. T and $T_0$ correspond to the current and initial temperature respectively.

When mechanical and thermal problems are coupled, a staggered solution approach is usually adopted, which solves the mechanical and the thermal problems in an uncoupled manner with data exchange performed at the end of each time step or increment. In particular, the nodal temperatures are transferred to the mechanical procedure, while plastic work is communicated to the thermal solver.
<figure><img src="image_4.png" /><figcaption>Figure 3. Representative domain around a node for the evaluation of $\tilde{\mathbf{F}}_i$.</figcaption></figure>

## 4. Refinement and adaptivity in the C-NEM

In this section, an adaptive strategy is proposed for non-linear problems with history-dependent internal variables in the C-NEM context, including: (a) a Zienkiewicz-Zhu error indicator [47] based on equivalent plastic strain; (b) a transfer technique based on the stabilized conforming nodal integration [10] and (c) a refinement strategy for domains evolving in time based on the Voronoi cells.

### 4.1. Transfer of internal variables with nodal integration

In the context of the finite element method, the transfer of internal variables between successive remeshing operations is an important issue (see i.e. [19, 7, 30]). In this work, we use the stabilized conforming nodal integration proposed by Chen et al. [10] to define all variables at the nodes, in order to avoid projection between two successive actualization of the reference configuration (update of the Voronoi diagram and shape functions).

We thus consider a deformation gradient at node $n_i$:

$$ \tilde{\mathbf{F}}_i = \frac{1}{|\Omega_i|} \int_{\Omega_i} (\nabla_X u^h(\mathbf{x}) + \mathbf{1}) d\Omega = \mathbf{1} + \tilde{\nabla}_X u^h(\mathbf{x}_i) = \mathbf{1} + \tilde{\mathbf{B}}_i \mathbf{u}_i \quad (51) $$

where $\Omega_i$ is a representative domain around the node (typically a Voronoi cell) (see fig. 3). Introducing C-NEM discretization scheme into (51) we obtain:

$$ \tilde{\nabla}_X u^h(\mathbf{x}_i) = \begin{bmatrix} u_{,X}^h(\mathbf{x}_i) \\ u_{,Y}^h(\mathbf{x}_i) \\ v_{,X}^h(\mathbf{x}_i) \\ v_{,Y}^h(\mathbf{x}_i) \end{bmatrix} = \tilde{\mathbf{B}}_i \mathbf{u}_i \quad (52) $$
where $\tilde{\mathbf{B}}_i$ is expressed by:

$$ \tilde{\mathbf{B}}_i = \begin{bmatrix} \tilde{\phi}_{1,X} & 0 & \tilde{\phi}_{2,X} & 0 & \dots & \tilde{\phi}_{N,X} & 0 \\ \tilde{\phi}_{1,Y} & 0 & \tilde{\phi}_{2,Y} & 0 & \dots & \tilde{\phi}_{N,Y} & 0 \\ 0 & \tilde{\phi}_{1,X} & 0 & \tilde{\phi}_{2,X} & \dots & 0 & \tilde{\phi}_{N,X} \\ 0 & \tilde{\phi}_{1,Y} & 0 & \tilde{\phi}_{2,Y} & \dots & 0 & \tilde{\phi}_{N,Y} \end{bmatrix} \quad (53) $$

with

$$ \tilde{\phi}_{i,X} = \int_{\Omega} \phi_{i,X}(\mathbf{X}) d\Omega \quad (54) $$

$$ \tilde{\phi}_{i,Y} = \int_{\Omega} \phi_{i,Y}(\mathbf{X}) d\Omega \quad (55) $$

where $\mathbf{u}_i$ are the nodal displacements. In the context of a Lagrangian procedure, different options can be considered:

1) A *total* Lagrangian procedure, where the Voronoi diagram as well as the shape functions are computed only once at the beginning of the simulation;

2) An *updated* Lagrangian procedure where the Voronoi diagram and the shape functions are adapted every time steps.

3) A *periodically updated* Lagrangian procedure where the Voronoi diagram and the shape functions are computed after several total Lagrangian steps.

In the context of a total Lagrangian procedure, no internal variable transfer is necessary. Nevertheless, in updated Lagrangian procedures (2) and (3), it is necessary to update the reference configuration, which implies to reconstruct the new Voronoi diagram associated with the new nodal distribution, as well as the new shape functions. In that case, the use of the assumed gradient defined in (51) allows the definition of all internal variables at the nodes. In the case of an explicit procedure, all quantities defined in the former time step are defined at the nodes in the same way. If the same cloud of nodes is used through the whole simulation, the transfer is thus direct as integration points coincide with the nodes. If a different cloud of node is used during the simulation (due to the refinements or nodal repositioning, a direct C-NEM interpolation (7) can be performed to define the internal variable the new nodes.

## 4.2. Error indicator based on equivalent plastic strain

In the context of the finite element method, adaptive strategies have been extensively applied and continuously developed for linear [17] [1] [47] and some class of nonlinear problems and history-dependent nonlinear problems over the last two decades or so (see i.e. Ladevèze et al. [16], Belytschko et al. [5], Ortiz et Quigley [28], Gallimard et al. [13], among many others).

The present section aims to propose a simple error indicator for adaptive solutions of large elasto-plastic transformations in the C-NEM framework. In the context of
the finite element method, an elementary procedure for the error estimation may be defined by substituting the exact solution by some post-processed field obtained from the available FEM solution. When the finite element solution is accurate enough, the post-processed solution is expected to be more accurate than the original FEM solution. In particular, the a posteriori error estimation procedure originally proposed and used by Zienkiewicz and Zhu [47] for linear elliptic problems is based on the observation that exact stress σ may be represented accurately by smoothed stress σ* obtained by a suitable projection of stresses σ^h.

In the context of meshless methods, the field σ^h is smooth in the general case, unlike in the FEM. Many options have been considered to construct error indicators in meshfree methods, mostly using some recovery Zienkiewicz-Zhu fields (see the work by Liu et al. [23], You et al. [44], Chung et Belytschko [8], Lee et Zhou [21, 20] and Lu et Chen [25] for an overview of the recent proposed techniques). In most of the referred papers in which a MLS approximation is used, an important issue is the influence of the shape function support size on the efficiency of error estimates [8]. Furthermore, it has been shown in [20] that adaptive refinement in MLS meshless techniques is a delicate task due to the necessity to adapt locally the shape function support according to the local nodal density. In the following, an error indicator based on the NEM shape functions is proposed to circumvent the just referred difficulties, as the NEM shape functions support automatically adapts its shape and size to the surrounding neighborhood.

We propose a simple error indicator using the Zienkiewicz-Zhu idea in tandem with the stabilized conforming nodal integration. Let α(X, t) a variable either associated with the spatial derivatives (i.e. F(X, t)), or with internal variables (i.e. ϵ^p(X, t)). In the C-NEM context, the stabilized conforming nodal integration scheme proposed by Chen et al. [43] produces constant piece-wise fields associated with the different derivatives, discontinuous across the Voronoi cells. This is a consequence of the stabilization scheme used in (51). In addition, the different internal variables can be considered at the nodes, as we have nodal integration. Constant fields α^h(X, t) associated with variables α(X, t) can thus be considered in each Voronoi cell Ω_i. A simple solution for recovery fields α^*(X, t) with assumed better accuracy is to interpolate the nodal values of α^h(X, t) with the C-NEM interpolation scheme:

$$ \alpha(\mathbf{X}, t)^* = \sum_{i=1}^{V} \phi_i^C(\mathbf{X}) \alpha^h(\mathbf{X}_i, t) \quad (56) $$

Then, the error based on the equivalent plastic strain for any cell Ω_i may be defined in the reference configuration, as:

$$ |e|_{\Omega_i}^2 = \int_{\Omega_{i0}} ([\bar{\epsilon}^p(\mathbf{X}, t)]^* - [\bar{\epsilon}^p(\mathbf{X}_i, t)]^h)^2 J_i d\Omega_0 \quad (57) $$

with:
$$[\bar{\epsilon}^p(\mathbf{X}, t)]^* = \sum_{j=1}^{V} \phi_j^c(\mathbf{X}) [\bar{\epsilon}^p(\mathbf{X}_j, t)]^h, \quad J_i = \det(\tilde{\mathbf{F}}_i) \qquad (58)$$

Where $[\bar{\epsilon}^p(\mathbf{X}_j, t)]^h$ is the equivalent plastic strain associated with the natural neighbor $n_j$ of point $\mathbf{x}$, assumed constant in the constrained Voronoi cell $\Omega_j$. In order to evaluate (57), the constrained Voronoi cells are triangulated and standard Gauss quadrature is applied on the triangles.

Any kind of error indicators can be constructed on the above framework, i.e. error indicators based on rate of plastic work or based on damage [30]. For the sake of simplicity, we only focus in this study on the error indicator described below. The global error is obtained, in a standard way, by : $|e|_{\Omega}^2 = \sum_{i}^{N} |e|_{\Omega_i}^2$. In addition, the relative error $\eta$ is defined as $\eta_i := |e|_{\Omega_i} / |e|_{\Omega}$

### 4.3. Refinement strategy based on the Voronoi cells

To adapt the nodal density to the evolution of plastic deformation, we must define a limit global error $\bar{\eta}$ and a local (in each cell) limit error $\bar{\eta}_i$ satisfying:

$$\bar{\eta}^2 = \sum_{i=1}^{N} (\bar{\eta}_i)^2 \qquad (59)$$

If we assure a uniform limit error $\bar{\eta}_i = \bar{\eta}_*$, $\forall i$ then:

$$\bar{\eta}_* = \frac{\bar{\eta}}{\sqrt{N}} \qquad (60)$$

Following Zienckiewicz et. al [48], we propose to determine the characteristic length of the new cell compared to its parent cell (which has been divided to generate new cells). For this purpose, we assume a rate convergence in energy norm of $p$ ($O((h^p))$). We thus have:

$$\frac{h_i^{new}}{h_i^{old}} = \left[ \frac{\bar{\eta}_*}{\eta_i} \right]^{\frac{1}{p}} = \left[ \frac{\bar{\eta}|e|_{\Omega}}{\sqrt{N}|e|_{\Omega_i}} \right]^{\frac{1}{p}} \qquad (61)$$

As noticed in former studies [46, 39], the convergence rate in the natural element method is of the same order as that obtained in the FEM using linear finite elements (due to the linear consistency of the Sibson interpolant), that is $p = 1$. To determine when a cell has to be subdivided, we assume that the characteristic length of the offspring cells are twice lower (see fig. (4)), for sake of simplicity.

$$\frac{h_i^{new}}{h_i^{old}} \le \frac{1}{2} \qquad (62)$$
<figure><img src="image_5.png" /><figcaption>Figure 4. Refinement procedure based on the Voronoi cells. (o): additional nodes.</figcaption></figure>

The criterion for Voronoi cell subdivision is then defined by:

$$ \left[ \frac{\bar{\eta}|e|_{\Omega}}{\sqrt{N}|e|_{\Omega_i}} \right]^{\frac{1}{p}} \leq \frac{1}{2} \qquad (63) $$

Figure 4 shows different refinement possibilities based on the Voronoi diagram. On the left is depicted the Voronoi diagram before insertion of new nodes. On the right is depicted the updated Voronoi diagram. In the strategy depicted in figure 4 (a), new nodes are added on the vertex of the Voronoi cell when the criterion (63) is satisfied. This strategy has been used by several authors [44] [25] for refinement in unstructured scattered of nodes. In the strategy depicted in figure 4 (b), new nodes are added between the node associated with the cell and the neighbors, when the criterion (63) is satisfied. We recommend to use the strategy (b), as in the case of strategy (a), for regular grids, several additional nodes can be inserted at the same position. These two strategies allow refinement in unstructured set of nodes and are consistent with the error indicators based on the Voronoi cells described below. In addition, no parameter is associated with the shape functions support size (see Fig. 2), which simplifies the adaptive analysis when the nodal density is significantly different in some parts of the domain.
### 4.4. Dynamic refinement strategy

We propose the following strategy for dynamic refinement: starting for the reference configuration (at time $t_0$), several Lagrangian increments are performed until a chosen time $t_a$. At time $t_a$, errors are evaluated using (57). If maximal accepted error is reached in some cells, new nodes are added according to the refinement strategy described in fig 4 (b) in the reference configuration. The computation then restart from time $t_0$ with the new nodes until time $t_a$. The procedure continues while the error is lower than the prescribed limit in all the cells of the domain. Then the configuration is updated at time $t_a$, which will constitute the new reference configuration.

## 5. Numerical example

High-speed shearing is a new metal cutting operation which offers significant advantages compared to traditional shearing. In this process, the high velocity of the punch induces adiabatic shear bands which combined with the damage phenomena produces high-quality of the cutting zone, without any burr formation. This process allows to perform shearing in high strength material such as titanium or stainless steels, or to produce holes in plates with relatively large thickness (up to 20mm) in traditional metals such as aluminum.

The example shows a first insight into the simulation of a plane-strain, high-speed shearing operation using the error indicator and the refinement procedure described previously. The geometry of the problem is depicted in Fig. 5. The following material data for stainless steel 304L have been used: $A = 253$ MPa, $B = 685$ MPa, $C = 0.0973$ $n = 0.312$, $m = 2.044$, $\bar{\epsilon}_0^p = 1$; $T_0 = 296$ K, $T_f = 1698$ K, $E = 210$ GPa, $\nu = 0.33$, $c = 500$ J/kgK, $k = 5.86$ J/JgmK, $\rho = 7850$ kg/m$^3$. The initial velocity of the punch is taken as $V = 11$ m.s$^{-1}$.

The main aim was to illustrate the capability of the refinement procedure to capture the details of the shear bands when different geometries of the tools are used. The high speed of the punch, combined to the material characteristics, induces high strain localization. The high rate of energy generation induced by the plastic deformation and the low thermal diffusivity of the material causes localized heating, which in turn, causes thermal softening and consequently localized plastic deformation. Two test have been performed using the following geometrical parameters: $R_p = 5$ mm, $L = 20$ mm, $H = 10$ mm and $j = 0.2$ mm.

In a first test, zero corner radii $R_1$ and $R_2$ have been used. Figure 6 shows the refined nodal distribution during the shear band formation. Fig. 7 illustrates the evolution of the underlying Voronoi diagram used for the refinement procedure in half of the workpiece. Fig. 8 depicts the equivalent plastic strain. In this first example, the direction of the shear band is observed parallel to the punch speed direction.

In a second test, corner radii $R_1 = 0.5$ mm and $R_2 = 0.5$ mm have been used. We can notice that the increase in the corner radius induces inclined shear bands, which
<figure><img src="image_6.png" /><figcaption>**Figure 5.** *High-speed shearing: geometry of the problem.*</figcaption></figure>

causes conical deformation in the workpiece and lower quality of the manufactured workpiece. The present procedure was able to capture the phenomena, as illustrated from Fig. 9 and Fig 10.

**6. Conclusion**

In this paper, a framework for practical implementation of the natural element method in the context of explicit thermo-elastoplasticity involving a strategy for dynamic adaptive refinement in this framework has been presented. The present procedure offers some outstanding advantages in large transformations including localization: (a) the C-NEM provides equivalent accuracy than the quadrilateral/hexahedral finite elements by only using the Delaunay triangulation of the current set of nodes; (b) the meshfree features of the technique allows frequent updates of the reference configuration without special attention on the relative position of the nodes; (c) the use of a stabilized conforming nodal integration simplifies the transfer of internal variables between successive updates as integration points coincide with the nodes; (d) the underlying Voronoi diagram allows refinement through unstructured scattered of nodes; (e) no size parameter is involved in the definition of the neighbors, which is essential in refinement procedure, where the nodal density is different in some parts of the domain. An example involving dynamic shear bands progression in metal has been presented to illustrate the potentiality of the proposed technique.
<figure><img src="image_7.png" /><figcaption>**Figure 6.** High-speed shearing ($R_1$ = 0 mm and $R_2$ = 0 mm): adaptive refinement during shear band formation</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 7.** High-speed shearing ($R_1$ = 0 mm and $R_2$ = 0 mm): adaptive refinement during shear band formation and constrained Voronoi cells.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 8.** High-speed shearing ($R_1$ = 0 mm and $R_2$ = 0 mm): shear band progression (equivalent plastic strain).</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 9.** High-speed shearing ($R_1$ = 0.5 mm and $R_2$ = 0.5 mm): adaptive refinement during shear band formation.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>**Figure 10.** *High-speed shearing (R₁ = 0.5 mm and R₂ = 0.5 mm): shear band progression (equivalent plastic strain).*</figcaption></figure>

**7. References**

[1] Babuška I., Rheinboldt W.C. Adaptive approaches and reliability estimations in finite element analysis. *Computer Methods in Applied Mechanics and Engineering* 1979; **17/18**:519–40;

[2] Babuška I., Melenk J.M., The partition of unity finite element method: basic theory and applications. *Computer Methods in Applied Mechanics and Engineering* 1996; **4**:289–314.

[3] Bathe K.J., *Finite Element Procedures* Prentice Hall, 1986.

[4] Belikov V.V., Ivanov, V.D., Kontorovich V.K., Korytnik S.A., Semenov A.Y., The non-Sibsonian interpolation : a new method of interpolation of the values of a function on an arbitrary set of points. *Computational Mathematics and Mathematical Physics* 1997; **37**(1):9–15.

[5] Belytschko T., Wong B.L., Plakacz E.J., Fission-fusion adaptivity in finite elements for nonlinear dynamics of shells. *Computers and structures* 1989; **33**:1307–1323.

[6] Belytschko T., Lu Y.Y., Gu L. Element-free Galerkin methods. *International Journal for Numerical Methods in Engineering* 1994; **37**:229–256.

[7] Camacho G.T., Ortiz M., Computational modeling of impact damage in brittle materials, *International Journal of Solids and Structures* 1996; **33**:2899–2938.
[8] Chung H-J., Belytschko T., An error estimate in the EFG method, *Computational Mechanics* 1998; **21**:91–100.

[9] Cueto E, Doblaré M, Gracia L., Imposing essential boundary conditions in the natural elements method by means of density-scaled alpha-shapes. *International Journal for Numerical Methods in Engineering* 2000; **49**:519–546.

[10] Chen J.S., Yoon Y., Wu C.T., Non-linear version of stabilized conforming nodal integration for Galerkin mesh-free methods. *International Journal for Numerical Methods in Engineering* 2002; **53**:2587–2615.

[11] Cuitiño A., Ortiz M., A material-independent method for extending stress-update algorithms from small-strain plasticity to finite plasticity with multiplicative kinematics. *Engineering computations* 1992; **9**:437–451.

[12] Eterovic A.L., Bathe K.-J., A hyperelastic-based large strain elasto-plasstic constitutive formulation with combined isotropic-kiematic hardening using the logarithmic stress and strain measures. *International Journal for Numerical Methods in Engineering* 1990; **30**:1099–1114.

[13] Gallimard L., Ladevèze P., Pelle J.P., Error estimation and adaptivity in elasto-plasticity. *International Journal for Numerical Methods in Engineering*. **39**:189–217, 1996.

[14] Hiyoshi H, Sugihara K, Improving continuity of Voronoi-based interpolation over Delaunay spheres. *Computational Geometry*. 2002; **22**:167–183.

[15] Johnson G.R., Cook W.H. A constitutive model and data for metals subjected to large strains, high strain rates and high temperatures. In *Proceedings of the Seventh International Symposium on Ballistic, The Hague, The Netherlands*, page 54, 1983.

[16] Ladavèze P., Coffignal G., Pelle J.P., Accuracy of elastoplastic and dynamic analysis, In *Babuška et al. eds.*, Accuracy estimates and adaptive refinements in finite element computations, John Wiley, Ney York, june 1998; pp. 86–95, association for computing machinery

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

[18] Lee E.H., Elastic-plastic deformation at finite strains. *J, Appl. Mech.* 1969; **36**:1–6.

[19] Lee N.-S., Bathe K.-J. Error indicators and adaptive remeshing in large deformation analysis. *Finite Elem. Anal. Des.* 1991; **16**:99–139.

[20] Lee C.K., Zhou C.E. , On error estimation and adaptive refinement for element free galerkin method: Part I: stress recovery and a posteriori error estimation. *Computer and structures* 2003; **82**(4-5): 4293–443;

[21] Lee C.K., Zhou C.E. , On error estimation and adaptive refinement for element free galerkin method: Part II:adaptive refinement. *Computer and structures* 2003; **82**(4-5): 413–428.

[22] Liu W.K., Jun S., Zhang Y.F., Reproducing Kernel Particle Methods. *Int. J. Numer. Methods Fluids*. 1995; **21**:1081–1106.

[23] Liu W.K., Uras R.A., Chen Y. Enrichment of the finite element method with the reproducing kernel particle method. *Journal of Applied Mechanics*, ASME 1997; **64**:861–870.

[24] Lucy L.B., A numerical approach to the testing of fusion process. *The astronomic journal* 1977; **88**:1013–1024.

[25] Lu H., Chen J.S., Adaptive meshfree particle method. *Lecture notes in Computational Science and Engineering* 2002; **26**:251–267.
[26] Marusich T., Ortiz M., Modelling and simulation of high-speed machining. *International Journal for Numerical Methods in Engineering* 1995; **38**:3675–3694.

[27] Nayroles B., Touzot G., Villon P. Generalizing the finite element method: diffuse approximation and diffuse elements. *Computational mechanics* 1992; **10**:307–318.

[28] Ortiz M., Quigley J.J., Adaptive mesh refinement in strain localization problems. *Comput. Methods Appl. Mech. Engrg.* 1991; **90**: 781–804.

[29] Perić D., Owen D.R.J., A model for finite strain elastoplasticity based on logarithmic strains : computational issues. *Comput. Methods Appl. Mech. Engrg.* 1992; **94**:35–61.

[30] Perić D., Vaz M., Owen D.R.J., On adaptive strategies for large deformations of elasto-plastic solids at finite strains: computational issues and industrial applications. *Comput. Methods Appl. Mech. Engrg.* 1999; **176**:279–312.

[31] Sambridge M., Braun J, McQueen M. Geophysical parameterization and interpolation of irregular data using natural neighbors. *Geophysi. J. Int.* 1995; **122**:837–857.

[32] Schönhardt E., Über die Zerlegung von dreieckspolyedern in tetraeder. *Math. Annalen* 1928; **98**.

[33] Shewchuck JR, Tetrahedral mesh generation by delaunay refinement *Proceedings of the fourteenth annual symposium on computational geometry*, Minneapolis, Minnesota, june 1998; pp. 86–95, association for computing machinery

[34] Shewchuck JR, Sweep algorithms for constructing higher-dimensional constrained Delaunay triangulations *Proceedings of the sixteenth annual symposium on computational geometry*, Hong-Kong, june 2000; pp. 350–359, association for computing machinery

[35] Seidel R., Constrained Delaunay triangulations and Voronoi diagrams with obstacles In "1978-1988 Ten Years IIG" 1988; 178–191.

[36] Simo J.C., Algorithms for static and dynamic multiplicative plasticity that preserve the classical return-mapping algorithm schemes of the infinitesimal theory. *Computer methods in Applied Mechanics and Engineering* 1992; **99**: 61–112.

[37] Simo J.C., Hughes T.J.R., *Computational inelasticity* Springer-Verlag, New York, 1998.

[38] Simo J.C., Miehe C., Associative coupled thermoplasticity at finite strains : formulation numerical analysis and implementation. *Computer methods in Applied Mechanics and Engineering* 1998; **98**:41–104.

[39] Sukumar N., Moran B., Belytschko T. The natural elements method in solid mechanics. *International Journal for Numerical Methods in Engineering* 1998; **43**:839–887.

[40] Sibson R., A vector Identity for the Dirichlet tessellations. *Math. Proc. Camb. Phil. Soc.* 1980; **87**:151–155.

[41] Taylor G., Quinney H.. The latent energy remaining in a metal after cold working. *Proceedings of the Royal Society* 1934; **143**:307–326.

[42] YongChang C., Hehua Z., A meshless local natural neighbor interpolation method for stress analysis of solids. *Engineering analysis with boundary elements* 2004; **28**:607–613.

[43] Yoo J., Moran B., Chen J.-S., Stabilized conforming nodal integration in the natural-element method. *International Journal for Numerical Methods in Engineering* 1998; **60**:861–890.

[44] You Y., Chen J.S., Lu H., Filters, Reproducing Kernel, and Adaptive Meshfree Method. *Computational mechanics* 2003; **31**:316–326.
[45] Yvonnet J., Chinesta F., Lorong P., Rynckelynck D., The constrained natural element method (C-NEM) for treating thermal models involving moving interfaces, *International Journal of Thermal Sciences*, accepted.

[46] Yvonnet J., Ryckelynck D., Lorong P., Chinesta F., A new extension of the natural element method for non convex and discontinuous domains : the constrained natural element method (C-NEM). *International Journal for Numerical Methods in Engineering* 1998; **60**:1451–1474.

[47] Zienkiewicz O., Zhu J., A simple error estimator and adaptive procedure for practical engineering analysis. *International Journal for Numerical Methods in Engineering* 1987; **24**:337–357.

[48] Zienkiewicz O.C., Boroomand B, Zhu J.Z., Recovery procedures in error estimation and adaptivity Part I: Adaptivity in linear problems. *Comput. Methods Appl. Mech. Engrg.* 1999; **176**:111–125.