Simulation of Forming Processes by the α -Shapes-Based Natural Element Method

In this paper we review some of the authors’ more recent developments in the field of Forming Processes simulation by means of meshless methods. In particular, all simulations are performed by employing the Natural Element Method (NEM), which has shown some particular characteristics that make it appear as an appealing tool for this kind of problems. Applications include forging, aluminum extrusion and other related forming processes. Particularly, the treatment of the free surface deserves some comments, since it is done in this work by means of alpha-shapes, a particular shape constructor.


Introduction
Forming processes constitute a field of engineering where very complex phenomena occur and where very often multi-physic and coupled problems arise. In addition, these usually involve large transformations (i.e., large strains or large displacements.) This has traditionally posed some problems to the Finite Element technology, related to mesh distortion, remeshing, etc. Thus, forming processes seem to be a natural candidate for the use of meshless methods.
In this article we review some of the authors' work on meshless simulation of forming processes, particularly those related to aluminium extrusion and other processes where coupled thermo-mechanical problems are of utmost importance.
Although many different meshless methods exist nowadays, our developments have been focused on the so-called Natural Element Method (NEM) or, equivalently, Natural Neighbour Galerkin methods [4,18]. This choice has been motivated mainly by the unique features of the NEM among meshless methods, that include, for instance, its ability to exactly reproduce essential boundary conditions.

1
The NEM also allows for an updated Lagrangian description of the movement, and this is especially interesting when dealing with free surface flows, for instance. This framework can thus be competitive to existing ALE formulations, avoiding advection terms in the equations of conservation and their associated numerical complexity.
This article is organized as follows. Firstly, we describe the NEM and its main properties. In Section 3 we describe the technique used to track the free-surface during extrusion, forging or related processes. We then include some numerical examples showing the performance of the proposed method. The article is finished with some conclusions.

The Natural Element Method
The NEM is basically a Galerkin procedure where the essential field is approximated by any of the existing natural neighbour interpolation schemes [9,15,20]. Prior to the description of these schemes, it is necessary to review some basic geometrical concepts, such as the Voronoi diagram of a cloud of points or the Delaunay triangulation.

Natural Neighbour Interpolation
Models are constructed in our version of the NEM based on nodes only. Consider a cloud of irregularly distributed nodes N = {n 1 , n 2 , . . . , n N }. The Delaunay triangulation [5] of the set N, D, is the only triangulation of the cloud that verifies the empty circumcircle criterion, i.e., no circle containing the three nodes of a triangle contains an additional node of the cloud, see Figure 1. The dual structure of a Delaunay triangulation is the so-called Voronoi diagram [21], which is the unique decomposition of the space into non-overlapping cells containing the locus of points closer to a given node than to any other. Formally, where T I is the Voronoi cell and d(·, ·) represents the Euclidean distance. n represents the space dimension (two or three in the examples presented in this article). Two nodes whose Voronoi cells share one edge are called natural neighbours and hence the name of these interpolation schemes. The first, and most obvious, interpolation scheme based on natural neighbours is the so-called nearest neighbour or Thiessen interpolation [20]. If we give the nodal value to the whole associated Voronoi cell, we obtain a C −1 interpolation scheme. This interpolation scheme is not suitable for solving second-order partial differential equations, but has been employed in [8] to construct mixed velocity-pressure approximations for the simulation of incompressible media.  The most extended form of natural neighbour-based interpolation is due to Sibson [16]. For the definition of Sibson interpolation it is necessary to previously introduce the concept of second-order Voronoi cell. It is defined as the locus of the points that have the node n I as the closest node and the node n J as the second closest node: If a new point is added to a given cloud of points, the Voronoi cells will be modified by its presence. Sibson [15] defined the natural neighbour coordinates of a point x with respect to one of its neighbours I as the ratio of the cell T I that is transferred to T x to the initial cloud of points to the total area of T x . In other words, being κ(x) and κ I (x) the Lebesgue measures of T x and T xI respectively, the natural neighbour coordinates of x with respect to the node I is defined as The resultant shape function depends obviously on the relative position of the nodes. An example of a node surrounded by other six is depicted in Figure 3.
Recently, some new interpolation schemes based on the concept of natural neighbors have been proposed [9]. One of them, coined as Laplace or non-Sibsonian interpolation, has received considerable attention, since it involves magnitudes of one order less of the space dimension (i.e., the calculus of areas in three-dimensional problems, for instance, instead of volumes). If we define the cell intersection t I J = {x ∈ T I T J , J = I } (note that t I J may be an empty set) we can define the value Thus, the point x shape function value with respect to node 4 in Figure 4 is defined as where s J represent the length of the Voronoi segment associated to node J and n represents the number of natural neighbours of the point under consideration, x.
Derivatives of the Laplace shape function are not defined along the edges of the Delaunay triangles that lie inside its support (see [19]). For the purposes of the work here presented, Sibson's interpolation has been considered.

Main Properties of the NEM
Firstly, unlike most approximation techniques used in meshless methods, Sibson's interpolation scheme is strictly interpolant, i.e., the approximated surface pass through the data. This can be expressed as with δ I J the Kronecker delta tensor. The natural neighbour (Sibson) interpolant possesses linear consistency, and thus it is amenable to be used in the approximation of second-order partial differential equations. This can be demonstrated from the before mentioned partition of unity property and its ability to exactly reproduce linear fields (also known as local coordinate property, see [16]): One important property that derives from Equation (6) is that the imposition of essential boundary conditions in NE methods is straightforward. Like in the FEM, it is sufficient to prescribe nodal values to reproduce essential boundary conditions. In addition, natural neighbour interpolants can be strictly interpolant (up to linear precision) on the boundaries (convex or not), under very weak conditions, as demonstrated in [3] by using the theory of α-shapes. This proof was later generalised in [2]. Idelsohn and co-workers later adopted this same approach based on α-shapes in the context of the Meshless Finite Element Method [10]. Another method to impose essential boundary conditions in the context of the NEM by employing visibility criteria was developed by Chinesta and co-workers [22].
The exact imposition of essential boundary conditions means that no interior point of the domain takes influence on the boundary. This is of first importance in the simulation of processes where friction or other phenomena occurring on the boundary are significant.
In its application to two-and three-dimensional linear elasticity, the natural element method has proved to render more accurate results than linear triangular and tetrahedral finite elements respectively [1,19].
But, maybe the most important characteristic of the Galerkin scheme obtained by using natural neighbour interpolation is that the accuracy of the approximation is not affected by the distortion of the triangulation [18]. This confers the method the "meshless" character that allows us to implement updated Lagrangian descriptions of the processes. It is also important to note that natural neighbour interpolation is continuous when the data sites move continuously, while other mesh-based methods (i.e., finite elements based on a Delaunay triangulation) are not. This is especially important when dealing when updated Lagrangian descriptions of the movement.
In the following section we describe the technique used to track the free surface appearing in the description of the different forming processes.

The α-Shapes-Based Natural Element Method
When dealing with processes involving a free surface -which is commonly the case in classic forming processes simulations, such as forging or extrusion -it is extremely important to accurately track its position. In aluminium extrusion, for instance, non-uniform fluxes of the hot billet through the die can lead to non-straight profiles. Extrusion of hollow profiles, such as tubes, is particularly interesting from a geometrical point of view, since the flow of the hot metal should divide in the interior of the die and re-join just before exiting it.
In many Lagrangian or Lagrangian-Eulerian codes tracking of the free-surface is done by means of markers. The free surface segments (in 2d) or facets (in 3d) are detected and stored at each time step, see [11] for an elegant description of such algorithms in mould filling simulations. Previously, checking of inter-penetration of free-surfaces or development of holes should also be done. This is usually a cumbersome process.
In our approach a different approach has been developed, based on the use of shape constructors. Shape constructors are geometrical entities that give a continuous shape to a discrete cloud of nodes. One of such constructors is the family of α-shapes of the cloud [6].

The Family of α-Shapes of a Cloud of Points
A cloud of points defines a finite set of shapes that can be parameterized by the level of detail up to which we want to represent the geometry. The idea behind α-shapes is simple: the Delaunay triangulation will provide a straightforward connectivity of the nodes. We will eliminate from this triangulation those triangles (or tetrahedra) whose circumradius is greater than the desired level of detail, say α.
Formally, an α-shape is a polytope that is not necessarily convex nor connected, being triangulated by a subset of the Delaunay triangulation of the points. Let N be a finite set of points in R 3 and α a real number, with 0 ≤ α < ∞. A k-simplex σ T with 0 ≤ k ≤ 3 is defined as the convex hull of a subset T ⊆ N of size | T |= k + 1. Let b be an α-ball, that is, an open ball of radius α. A k-simplex σ T is said to be α-exposed if there exist an empty α-ball b with T = ∂b N where ∂ means the boundary of the ball. In other words, a k-simplex is said to be α-exposed if an α-ball that passes through its defining points contains no other point of the set N.
Thus, we can define the family of sets F k,α as the sets of α-exposed k-simplices for the given set N. This allows us to define an α-shape of the set N as the polytope whose boundary consists on the triangles in F 2,α , the edges in F 1,α and the vertices or nodes in F 0,α .
A three-dimensional simplicial complex is a collection, C, of closed k-simplices (0 ≤ k ≤ 3) that satisfies: The intersection of two simplexes in C is empty or is a face of both.
Each k-simplex σ T included in the Delaunay triangulation, D, defines an open ball b T whose bounding spherical surface (in the general case) ∂b T passes through the k + 1 points of the simplex. Let T be the radius of that bounding sphere, then, the family G k,α , is formed by all the k-simplexes σ T ∈ D whose ball b T is empty and T < α . The family G k,α does not necessarily form simplicial complexes, so Edelsbrunner and Mücke [6] defined the α-complex, C α , as the simplicial complex whose k-simplexes are either in G k,α , or else they bound (k + 1)-simplexes of C α . If we define the underlying space of C α , |C α |, as the union of all simplexes in C α , the following relationship between α-shapes and α-complexes is found: In order to clarify the before presented concepts, consider some examples of αshapes computed from a cloud of points corresponding to one particular simulation of a two-dimensional extrusion process. We restrict ourselves to geometrical concepts only.
Consider the extrusion example shown in Figure 5, where the contour plot of equivalent plastic strain rate is depicted. The key idea of the method here proposed is to extract the shape of the domain at each time step by invoking the concept of α-shape of the cloud of points. The α parameter will be obtained by geometrical considerations. In this case the radius at the inlet of the die, for instance, seems to be the smallest level of detail up to which the domain (i.e., the billet) must be represented. In order to appropriately represent this value, the nodal distance h must be accordingly chosen. In Figure 6 some members of the family of α-shapes of the cloud of points in its final configuration (corresponding to Figure 5(b)) are depicted. In Figure 6(a) the member for α = 0, i.e., the cloud of points itself, is shown. Note how, as α is increased, the number and size of the simplexes (in this case, triangles) that belong to the shape is increasing. For α = 1.0 we obtain an appropriate shape for the cloud. Note, however, that this is not an exact value to be determined at each time step. Since the number of α-shapes is finite, there generally exists an interval of valid α values for a single shape. Finally, by increasing the α value, we arrive to the convex hull of the cloud of points (Figure 6(f)).

Properties of the α-Shapes-Based Natural Element Method
Constructing the model by taking into account the actual shape of the domain has consequences not only in geometrical aspects of the method, but also in the quality of the approximation. As demonstrated in [3], the use of a proper α-shape for the definition of the domain -which means that we reproduce the geometry to the desired level of detail -ensures the linear interpolation of the essential field along the boundary. This is extremely important, since many of the meshless methods lack of the so-called Kronecker delta property, i.e., do not interpolate the field along the boundary, thus complicating the imposition of essential boundary conditions.
In practice, this requirement means that an appropriate cloud of points should be used in order to achieve a correct interpolation on the boundary. Here, the term "appropriate" should be understood in the sense that the cloud of points is dense enough so as to properly reproduce the geometry up to the desired level of detail.
Other methods exist to impose essential boundary conditions in the NEM, such as the so-called C-NEM [22], that uses visibility criteria and renders exactly the same approximation. Both methods can be considered equivalent.

Governing Equations
Many of the existing codes consider rigid-viscoplastic models for the constitutive equations of metals during forming processes. The obvious advantage of this method is that the material can be modelled as a non-Newtonian fluid. This approach is usually known as the flow formulation of the problem [23]. Thus a mixed Sibson-Thiessen approximation is used for the essential variables of the problem, namely, velocity and pressure. Details of the stability of this formulation can be found in [8].
The deviatoric stresses will be, under this assumption, being, as usual, σ = s − pI (10) where p = −tr(σ )/3 and I stands for the second-order identity tensor. Obviously, in the most general case, the parameter μ will depend on both the level of strain (and hence the non-linear character of the behavior) and the temperature. To derive the expression of the parameter μ it is a common practice to write the strain rate tensor as emerging from a visco-plastic potential. Following Perzyna [13] where Y is the viscoplastic potential -usually coincident with the plastic criterion as has been considered here -γ is a scalar function given bẏ g is a monotonic function that takes zero value only if Y (σ , q) ≤ 0, η is a positive parameter often called viscosity and q represents the hardening parameters. In what follows we will avoid the use of the vp superscript to indicate viscoplastic if there is no risk of confusion. For metals, there exist well-defined plastic yield rules and for aluminium it is a common practice to employ a von Mises criterion: where σ = 3 2 s : s = 3J 2 represents the effective stress and σ y represents the uniaxial yield stress. d is the only internal variable in this model and is sometimes called effective strain rate: Among the most extended laws to account for the aluminium yield stress is the so-called hyperbolic sine or Sellars-Tegart law [14] S m , m and A are material parameters. A is a factor that depends on the magnesium and silicium matrix solute content (see [12] and references therein), Q represents the activation energy of the deformation process, R is the universal gas constant and, finally, T is the absolute temperature. Note that, as the temperature increases, the yield stress decreases, as expected. Following this model, a state of null strain rate will give a null yield stress, and this is not in accordance to the well-known behavior of metals in general, and aluminium in particular. So it is a common practice to correct this effect by adding to Equation (15) a initial strain rate, d 0 so as to render a modified Sellars-Tegart law: If we combine now the general form of the strain rate tensor given in Equation (11), with Equation (13), we arrive to It is immediate now, by combining Equation (14) and the definition of effective stress, σ , to prove thatγ is precisely the effective strain rate: On the other hand, and by following the Perzyna-like model employed in Equations (11) and (12) and taking g(f ) = f , we arrive to a relationship between equivalent stress and equivalent strain rate: that, introduced in Equation (17), accounting Equation (18), gives the following visco-plastic constitutive equation: Note that, depending on the η value, the return to the yield surface is done with different velocity. Since it is common to describe aluminium behaviour as rigidplastic (rather than viscoplastic) we employ null viscosity, so as to enforce Y = σ − σ y = 0, leading to Finally, the constitutive equation, accounting the incompressibility of plastic flow results: Of course, this simple model has important limitations. Undoubtedly, the lack of elastic behaviour is one of the most important. Thus, spring-back cannot be pre- dicted. However, as mentioned before, it has rendered good results and seems to be widely accepted in the forming processes community [12,23,24].
This same model can be employed to model a wide variety of polymers, for instance, governed by a power law (also known as Norton-Hoff model): (23) where n is the so-called sensitivity index. The interested reader can consult [17] for some applications on this class of problems.

Simulation of a Forging Process
In this example we consider the application of the proposed α-NEM method to the simulation of the forge of a complex piece. In this case a viscoplastic behaviour is assumed, given by Equation (23). We have considered μ 0 = 1.0 and n = 0.3. These values are obviously non-physical, but can help us to show the behavior of the proposed technique when extremely high deformations are present.
The geometry of the piece is shown in Figure 7(a). The simulation deals with the forging of the central region of the piece, justifying the assumption of plane strain (see Figure 7(b)).
Accounting for the symmetry of the geometry, only one half of the domain was simulated, by imposing appropriate boundary conditions. Slip boundary conditions were assumed at the billet-punch-die contact region. The upper part of the punch, assumed perfectly rigid, moves towards the lower part, fixed throughout the simulation. Geometry and dimensions of the simulated region are shown in Figure 8.
The equivalent plastic strain at time steps 96, 120, 150 and 168 are shown in Figure 9. Very accurate results were obtained in spite of the large strains and displacements involved in the simulation.
It can be seen how the proposed method is a valuable tool to accurately predict the formation of the flash.

Simulation of the Extrusion of a Cross-Shaped Profile
In [7] an analysis is made of the process of extrusion of a cross-shaped profile, with particular interest on the effect of the misalignment between the die and the workpiece. Geometry and dimensions of the die are shown in Figure 10. The hole of the die was misaligned 10 mm in order to study the influence of this defects in the resulting quality of the extruded profile. The initial mesh is shown in Figure 11. Only one half of the geometry was analyzed, by applying appropriate symmetry boundary conditions.
The material is assumed to be lead, which is able to flow at room temperature, and is therefore easily characterizable by simple experiments with an universal testing machine. In [7], the flow rule of lead at room temperature was adjusted to a Norton-Hoff law, giving σ = 60ε 0.05 [MPa] (24) in a uniaxial test. Various snapshots of the mesh at different time steps are shown in Figure 12. The evolution of the equivalent strain rate is depicted in Figure 13.
Results are in good agreement with the experimental results performed by Filice and colleagues [7] and the resulting deviation of the profile is only slightly underestimated.
This approach notably simplifies those of the ALE methods, especially when large motion is produced in the direction perpendicular to the main flow of the extrudate, which is traditionally treated as purely Eulerian, as in [24]. This is the case when the flow is distorted by the misalignment of the die hole and the axis of the extruder.
15 Fig. 12. Evolution of the geometry of the extrudate throughout the process.

Conclusions
Forming processes usually involve large deformations and are therefore specially interesting as an objective for meshless methods. We have reviewed some of the more important characteristics of the α-shape-based Natural Element Method (α-NEM, see [3]) and studied the more interesting features of the method applied to such problems.
We believe that meshless methods and, in particular, the α-NEM are particularly well-suited to simulate that class of processes. As is well-known for all meshless methods, distortion of the cloud of points does not greatly affect the quality of the results, and therefore remeshing is avoided, in the sense that the cloud of nodes isor can be -the same throughout the whole simulation. The connectivity is updated by the method in a process transparent to the user, thus greatly alleviating the user effort. Other aspects, such as an improvement of the speed of calculation are currently under investigation, since the computation of Sibson coordinates is known to be considerably harder than the computation of FE shape functions.
In sum, we think that meshless methods, and particularly the NEM, are a good choice for the simulation of problems involving large deformations in a Lagrangian framework.