Study of the high speed blanking process with the
C-Nem

Lounès Illoul, Philippe Lorong, Francisco Chinesta

► To cite this version:

Lounès Illoul, Philippe Lorong, Francisco Chinesta. Study of the high speed blanking process with the C-Nem. International Conference on High Speed Industrial Manufacturing Processes, 2007, Cenlis, France. 10.1007/s12289-008-0103-3 . hal-00283309

HAL Id: hal-00283309

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

Submitted on 28 Oct 2021

**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
# Study of the high speed blanking process with the C-NEM

L. Illoul¹, P. Lorong¹, F. Chinesta¹

¹L.M.S.P., U.M.R. C.N.R.S. n° 8106, Arts et Métiers ParisTech
151, boulevard de l'Hôpital, 75013 Paris

URL: www.paris.ensam.fr/lbsp

e-mail: philippe.lorong@paris.ensam.fr

**ABSTRACT:** The C-NEM is a new approach that we apply to solve mechanical problem in the context of finite transformations as encountered in forming process simulation. The C-NEM (Constrained Natural Element Method) is based on the Voronoi diagram associated to the nodes spread on the studied domain. The Voronoi diagram is constrain to respect the geometry of the domain boundaries. Its ability for constructing the interpolation when the domain becomes highly distorted makes of that approach an appealing choice for simulation forming processes as for example 3D cutting and blanking. In this communication we present some aspects related to the C-NEM approach and we show its ability to simulate kinematics encountered in high velocity blanking process.

**KEYWORDS:** C-NEM, Finite Transformations, Simulation of Forming Processes

## 1 INTRODUCTION

The C-NEM has still a number of common points with the finite elements approach: the interpolation is defined on a set of nodes, degrees of freedom (DOF) are associated to the nodes, a localised in space shape function is associated to each node, the interpolation is nodal and is exact for linear fields. The C-NEM interpolation uses the constrained Voronoi diagram (dual of the constrained Delaunay's tesselation) associated to the set of nodes and the boundaries of the domain. The quality of the produced interpolation depends primarily on the distribution of the nodes. The presence of flat tetrahedrons in the dual mesh of Delaunay does not affect the interpolation quality.

In next paragraphs we initially reconsider the NEM interpolation. We then focus on the non-convex domains for which the construction of a constrained Voronoi diagram is needed. The aspects related to the partitioning of the domain are then evoked. Finally we give the example of a 3D simulation of blanking.

## 2 THE C-NEM INTERPOLANT

### 2.1 Headlines of the natural element method

The NEM (Natural Element Method) interpolant is based on the Sibson's natural neighbor coordinates

(shape functions) [3, 10] and is constructed on the basis of the Voronoi diagram. For a set of nodes $S = \{n_1, n_2, ..., n_N\}$ in $\mathbb{R}^{\text{dim}}$, ($\text{dim} \in \{2, 3\}$), the Voronoi diagram is the subdivision of $\mathbb{R}^{\text{dim}}$ into regions $T_i$ (Voronoi cells) defined by:

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

Sibson's coordinates of x with respect to a natural neighbor $n_i$ (see figure 1) are 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 x:

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

<figure><img src="image_2.png" /><figcaption>Figure 1: Sibson shape functions construction in 2D</figcaption></figure>

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(\mathbf{x}_j) = \delta_{ij}$ ($\delta_{ij}$ being the Kronecker delta). The properties of positivity, interpolation, partition of unity [10] and local coordinate property [9] are verified:

$$ \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 \mathbf{x} = \sum_{i=1}^{n} \phi_i(\mathbf{x})\mathbf{x}_i \quad (3) $$

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 (tetrahedrons) connecting the node $n_i$ (in this case $n$ is the number of natural neighbors of node $n_i$). Another important property of the NEM interpolant is the ability to reproduce linear functions over the boundary of convex domains. The proof can be found in Sukumar et al. [10]. This is not true in the case of non convex boundaries, and the next section focuses on an approach to circumvent this difficulty.

## 2.2 The Constrained natural element method

In its original form [10], the NEM can only be applied to strictly convex domains. For strongly non-convex domains (cracks, auto-contact...) some spurious influences between nodes of the boundaries appear [11]. In order to avoid these drawbacks, we have proposed in a previous paper [11] 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 [5] (see [11] for further details).

The new cells $T_i^C$, called constrained Voronoi cells, 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_{x \to n_i} \cap \Gamma = \emptyset, S_{x \to n_j} \cap \Gamma = \emptyset \} \quad (4) $$

where $\Gamma$ is the domain boundary, composed by a set of segments in 2D (triangular plane facets in 3D) 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$. For both the trial and the test functions the new approximation is then:

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

where $V$ is the set of natural neighbors which are 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}$.

It was shown in [11] that the use of the constrained Voronoi diagram does not affect the properties of the NEM interpolation.

## 2.3 Constrained Delaunay tetraedrisation

The constrained Delaunay triangulation does not always exist in 3D without adding new nodes [4]. Some techniques for constructing 3D constrained Delaunay tessellations are available by adding Steiner points [6, 7]. These points are always set on triangles belonging to the boundary of the domain. It is possible to connect each of these points to at most 3 initial nodes (those defining triangles containing the new points). The new points are treated as Slave nodes, the initial nodes being the Master nodes. For each Slave node $n_i^S$, we set the following kinematic linear constraint:

$$ \mathbf{u}_i^S = \sum_{j=1}^{3} \eta_{ij} \cdot \mathbf{u}_j^{M(S^i)} \quad (6) $$

where $\mathbf{u}_j^{M(S^i)}$ are the nodal values at the 3 Master nodes of $n_i^S$, and $\eta_{ij}$ the reduced coordinates of the Slave node on its supporting triangle. To keep the NEM interpolation properties, the interpolation must be built using all the Master and Slave nodes values. With this partition, the interpolation becomes:

$$ \mathbf{u}^h(\mathbf{x}) = \sum_{i \in M_i} \phi_i^C(\mathbf{x}) \mathbf{u}_i^M + \sum_{i \in S_i} \phi_i^C(\mathbf{x}) \mathbf{u}_i^S \quad (7) $$

$M_i$ and $S_i$ being respectively the set of Master and Slave natural neighbors which are visible from $\mathbf{x}$. The software TETGEN (Quality Tetrahedral Mesh Generator and Three-Dimensional Delaunay Triangulator, [8]) is currently used to build the constrained Delaunay tetrahedrisation.

# 3 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 as follows:

$$ \begin{aligned} & \int_{\Omega^t} \rho(t) \ddot{\mathbf{u}} \cdot \boldsymbol{\eta} d\Omega^t + \int_{\Omega^t} \boldsymbol{\sigma}^t : \nabla_x \boldsymbol{\eta} d\Omega^t \\ &= \int_{\Omega^t} \rho(t) \mathbf{b} \cdot \boldsymbol{\eta} d\Omega^t + \int_{\partial \Omega^t} \tau \cdot \boldsymbol{\eta} d\Gamma^t \quad \forall \eta \in \mathcal{V} \end{aligned} \quad (8) $$
where $\rho$ is the density, b and $\tau$ represent the body forces and applied tractions respectively, and $\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:

$$ 
\begin{aligned}
& \int_{\Omega_0} \rho_0 \ddot{\mathbf{u}} \cdot \boldsymbol{\eta} d\Omega_0 + \int_{\Omega_0} \mathbf{P}^t : \nabla_x \boldsymbol{\eta} d\Omega_0 \\
& = \int_{\Omega_0} \rho_0 \mathbf{b} \cdot \boldsymbol{\eta} d\Omega_0 + \int_{\partial\Omega_\sigma^t} \boldsymbol{\tau} \cdot \boldsymbol{\eta} d\Gamma^t \quad \forall \boldsymbol{\eta} \in \vartheta
\end{aligned}
\quad (9) $$

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 (5) of the variational form (9) results in the discrete set of algebraic time dependent equations, 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 (10) $$

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.

As shown in the next section, the use of the SCNI quadrature [1] results in a $\mathbf{M}$ matrix diagonal, whose diagonal terms are given by $m_i = \rho_0 \Omega_i$, with $\Omega_i$ the volume 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:

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

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

## 4 NUMERICAL INTEGRATION

In the context of finite elements, a natural partition of the domain is carried out by the elements. With C-NEM the concept of element does not exist. A natural partitioning is given by the cells of the Voronoi diagram. When cells do intersect the boundary of the domain, only their part inside the domain is kept. However, for strongly non-convex domains in 3D, some boundary cells may have negative volume inside the domain (some faces may intersect themselves inside the domain). It is thus necessary to use another partitioning. For reasons of simplicity we have chosen a partitioning based on the Delaunay tetrahedrons: each elementary volume $\Omega_i$ is defined by the quarters of all tetrahedrons connected to the node $n_i$.

We then use the stabilized conforming nodal integration (SCNI) proposed by Chen et al. [1] to define a deformation gradient at node $n_i$:

$$ \tilde{\nabla} f_i^h = \frac{1}{|\Omega_i|} \int_{\Omega_i} \nabla f^h(\mathbf{x}) d\Omega \quad (13) $$

In order to avoid projections between two successive actualisations of the reference configuration all variables are computed with this gradient at the nodes.

## 5 NUMERICAL SIMULATION OF BLANKING

All the simulations given below are dealing with the same example. It is a simulation of blanking where the gap $e$ between the punch and the die is equal to 0.6mm and the radius $r$ of the cutting edges are equal to 0.2mm (Figure 2). The sample has a thickness (in the direction of cutting) equal to 5mm, and its behavior is modeled via a Johnson-Cook law [2] (Eq.14) with the following coefficients (SI): $A = 253.10^6$, $B = 685.10^6$, $C = 0.097$, $m = 2.044$, $n = 0.3128$, $\dot{\epsilon}_0^p = 1$, $T_0 = 296$. These coefficients correspond to a 306L stainless steel.

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

$\bar{\epsilon}^p$ represents the equivalent plastic strain, $\dot{\epsilon}^p$ the rate of plastic strain and $T$ the current temperature. The punch and the die are considered as rigid and the contact with the sample is without friction. The two symmetries of the experiment are taken into account and only a quarter of the sample is studied.

<figure><img src="image_3.png" /><figcaption>Figure 2: The blanking device</figcaption></figure>

Figure 3 gives the evolution of the Mises equivalent stress. The velocity of the punch is equal to $10m.s^{-1}$. The elastic behavior takes a major place in the localization of the shearing band. The elastic waves take $0.910^{-6}s$ to pass across the sample and figure 4 shows that the shearing band is clearly localized after
only 2.0 10⁻⁶s when elastic waves are still governing
the stress distribution in the sample. On figure 4 we
compare the equivalent plastic strain for 2 different
punch velocities.

<figure><img src="image_4.png" /><figcaption>Figure 3: Stress field in blanking - Punch velocity : 10ms⁻¹</figcaption></figure>

This figure exhibits the asymmetry of the shearing
band which is more important close to the die, this
effect being more important when the velocity of the
punch goes from 10 to 30m.s⁻¹.

<figure><img src="image_5.png" /><figcaption>Figure 4: Plastic strain localisation after 2.00µs</figcaption></figure>

On both previous simulations the punch and the die
are taken as rigid although their young modulus and
mass density are close to those of the sample. This
increases the effect of plastic strain near the contact
and inside the sample. This aspect will be taken into
account for future simulations.

6 CONCLUSIONS

With the possibility of dealing with three-dimensional
problems involving domains of any type (convex or
not), C-NEM is from now a general tool. Its develop-
ment is currently oriented towards the adaptation of
the node distribution (refinement in particular) spe-
cially in zones where high gradients occur, and also
to the adaptation of the boundary description when
high deformations have occurred. For instance, a cor-
rect description of the domain boundary is necessary

to update the constrained Voronoi diagram in the con-
text of an updated Lagrangian approach.
First results on blanking are nevertheless interesting.

ACKNOWLEDGEMENT

his work has been done under the technical and financial support
of CETIM Senlis (Centre Technique des Industries Mécanique).

REFERENCES

[1] J.S. Chen, Y. Yoon, and C.T. Wu. Non-linear version of stabilized conforming nodal integration for galerkin mesh-free methods. *Int. J. Numer. Meth. in Engng*, 53:2587–2615, 2002.

[2] G. R. Johnson and W. H. Cook. A constitutive model and data for metals subjected to large strains, high strain rates and high temperatures. In *Proceeding of the 7th International Symposium on Ballistics*, pages 541–547, The Hague, Holland, April 1983.

[3] M. Sambridge, J. Braun, and H. McQueen. Geophysical parameterization and interpolation of irregular data using natural neighbours. *Geophys. J. Int.*, 122:837–857, 1995.

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

[5] R. Seidel. Constrained delaunay triangulations and voronoi diagrams with obstacles. In *1978-1988 Ten Years IIG*, pages 178–191, 1988.

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

[7] JR Shewchuck. Sweep algorithms for constructing higher-dimensional constrained delaunay triangulations. In *Proceedings of the sixteenth annual symposium on computational geometry*, pages 350–359, June 2000. association for computing machinery.

[8] H. Si and K. Gärtner. Meshing piecewise linear complexes by constrained delaunay tetrahedralizations. In *Berlin/Heidelberg B.W. HANKS, Springer, editor, 14th International Meshing Roundtable*, pages 147–163, 2005.

[9] R. Sibson. A vector identity for the dirichlet tesselations. *Math. Proc. Camb. Phil. Soc.*, 87:151–155, 1980.

[10] N. Sukumar, B. Moran, and T. Belytschko. The natural elements method in solid mechanics. *Int. J. Numer. Meth. in Engng*, 43:839–887, 1998.

[11] J. Yvonnet, D. Ryckelynck, P. Lorong, and F. Chinesta. A new extension of the natural element method for non convex and discontinuous domains : the constrained natural element method (c-nem). *Int. J. Numer. Meth. in Engng*, 60:1451–1474, 2004.