# Recent advances in the meshless simulation of aluminium extrusion and other related forming processes

Icíar Alfaro, David González, David Bel, Elías Cueto, Manuel Doblaré, Francisco Chinesta

► To cite this version:

Icíar Alfaro, David González, David Bel, Elías Cueto, Manuel Doblaré, et al.. Recent advances in the meshless simulation of aluminium extrusion and other related forming processes. Archives of Computational Methods in Engineering, 2006, 13 (1), pp.3-43. 10.1007/BF02905930. hal-01007236

HAL Id: hal-01007236

https://hal.science/hal-01007236

Submitted on 21 Jan 2017

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 Attribution 4.0 International License
# Recent Advances in the Meshless Simulation of Aluminium Extrusion and other Related Forming Processes

I. Alfaro, D. González, D. Bel, E. Cueto, M. Doblaré

Group of Structural Mechanics and Material Modeling
Aragón Institute of Engineering Research (I3A) University of Zaragoza
Betancourt Building. María de Luna, 5. E-50018 Zaragoza, Spain
ecueto@unizar.es

F. Chinesta

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

In this paper we review some recent results in the field of numerical simulation of extrusion and other forming processes obtained by the authors by using a meshless approach, together with a wide review of the existing bibliography on the topic. Three main alternatives exist in the literature, namely (updated) Lagrangian, Eulerian and arbitrary Lagrangian-Eulerian (ALE) methods. A review of the most important characteristics of each of these three approaches is here presented and their possible advantages are pointed out. Finally, an updated Lagrangian approach over a meshless approximation, based on a class of methods globally coined as natural element methods (also as natural neighbour Galerkin methods) is analysed and its relative advantages studied. Some numerical examples are included that clearly show the potential capabilities of the proposed method.

## 1 INTRODUCTION

Aluminium extrusion is one of the most extended forming processes in industry. It consists, roughly speaking, in forcing an aluminium billet to pass through a die, thus producing a long profile whose cross section closely resembles that of the die. A simplified, schematic, representation of an extrusion process is shown in Figure 1, together with the typical evolution of the force in the ram during the whole process.

In this paper we briefly review the state of the art in the numerical simulation of this specific forming process. Some recent advances in the application of meshless methods to the simulation of forming processes are also reviewed. In particular, the application of the Natural Element Method (NEM) to the simulation of non-isothermal aluminium extrusion is analyzed in a higher detail.

### 1.1 Numerical Simulation of Extrusion and other Forming Processes

The numerical simulation of extrusion and other processes is of utmost importance in order to avoid the traditional trial-and-error approaches in the equipment design. The first attempts to perform extrusion simulations date back to the late sixties (see [56] and references therein) but the most extended approach during the last decades has been Finite Elements couples with the so-called *flow formulation* due to Zienkiewicz and co-workers [59] [57] [58].

Although, to some extent, an elastic recovery exists at the die outlet, this is often neglected. In addition, the Cauchy stress is usually related to the strain rate tensor. This leads to a formulation that closely resembles that of non-newtonian fluids, and hence the term *flow formulation* [57].
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Scheme of an extrusion process and associated evolution of the pressure on the ram</figcaption></figure>

Thus, the equations governing the metal deformation can be expressed in terms of velocities rather than displacements. Stresses produced in the forming process can be set in a simple form as

$$ \sigma = D(d, T) \cdot d, \qquad (1) $$

where *d* represents the strain rate tensor (symmetric part of the velocity gradient) and *T* the temperature.

More recently, Bonet [10] presented a formulation for this class of non-newtonian flow in terms of displacements. He called it the *incremental flow* formulation. Among its advantages is the possibility of direct imposition of incompressibility, since the divergence-free condition is imposed on displacements themselves rather than in incremental form, typical of fluids.

Most researchers assumed (see [56], for instance) a stationary character for the extrusion process. This assumption allows the use of a fixed mesh (Eulerian) frameworks. Although an intermediate stage exists in which this assumption can be considered as sufficiently accurate —at least as a first approximation, see Figure 1—, a pressure drop due to both the temperature increase during the process and the decreasing billet length can be observed [55]. Near the end of the process the pressure raises again as the material of the remaining billet tends to flow radially towards the die aperture. Nevertheless, it is possible to employ an Eulerian approach to simulate non-stationary processes, but some corrections must be added to adjust the geometry of the free surface of the flow, i.e., at the die outlet. These corrections are often performed iteratively (see, for instance, [33]).

It seems therefore interesting to develop numerical methods able to simulate the whole extrusion cycle. Similar comments can be made to some extent for other forming processes such as rolling. The use of (updated) lagrangian approaches appears to be appealing and this possibility has been studied since the first works on this field [56]. However, this approach has also some important drawbacks. Consider, for instance, a non-isothermal
process, in which the heat transfer equation in Lagrangian coordinates reads:

$$ \rho c_p \frac{dT}{dt} + \nabla(k\nabla T) + \dot{r} = 0 \qquad (2) $$

where k denotes the thermal conductivity, r the heat generation rate due to plastic dissipation during the forming process, ρ the specific density and c_p the specific heat of the metal. If Eulerian coordinates are chosen, the time derivative must be replaced by a material derivative, and this equation will read

$$ \rho c_p \left( \frac{\partial T}{\partial t} + v \nabla T \right) + \nabla(k\nabla T) + \dot{r} = 0. \qquad (3) $$

It is well known that in this last case some stabilization procedure must be used due to the presence of the additional convective term in Eq.(3). These procedures are not necessary when a Lagrangian procedure is followed.

On the contrary, the most obvious drawback of Lagrangian procedures is mesh distortion as the material flows. Good quality three-dimensional mesh generation algorithms are still a matter of research, although some commercial codes that employ tetrahedral meshes and extensive remeshing exist. Many researchers have preferred to employ *Arbitrary Lagrangian Eulerian* (ALE) procedures to reduce this problem. In these ALE methods the mesh is moved with a prescribed velocity, independent of the fluid motion. The heat equation (2) now looks like

$$ \rho c_p \left( \frac{\partial T}{\partial t} + (v - \bar{v}) \nabla T \right) + \nabla(k\nabla T) + \dot{r} = 0, \qquad (4) $$

where $\bar{v}$ represents this prescribed velocity. Clearly if $\bar{v} = v$ we are back to a Lagrangian formulation, whereas if $\bar{v} = 0$ we recover the Eulerian formulation (and hence the name of the method). The objective is to minimize mesh distortion and therefore the remeshing stages of the method. Some examples of the application of ALE formulations to the simulation of aluminium extrusion were presented by Huétink and co-workers [3] [40] [35] [51]. Other recent contributions on this topic can be found in [34] [33] [32]. Unfortunately, the prescribed velocity value $\bar{v}$ is not easy to determine and must be carefully chosen. In any case, the convective term, although is not so important, still remains.

## 1.2 Meshless Methods

In the last decade a new possibility has opened with the irruption of meshless methods [41] [9] [6] [48] [23]. These methods, although of very different nature, share the characteristic of no need of explicit definition of the nodal connectivity. Also —what is probably more important— their accuracy is not greatly affected by the quality of nodal distribution. In order words, they do not suffer from *mesh distortion*.

There exist meshless methods based on both Galerkin and collocation approaches. Probably, among the most popular ones are the so-called *moving least squared methods* [41] [9]. Due to the youth of these methods some drawbacks persist, although a considerable improvement has occurred in the last years.

One of the most important problems with meshless methods is related to the imposition of essential boundary conditions. Since moving least squares approximations are not strictly interpolant, it is not possible to directly prescribe nodal values in the discrete system of equations emerging from the Galerkin method, as in the FEM, by directly substituting nodal values by the prescribed ones. Instead, some especial technique like Lagrangian multipliers [9], modified variational principles or other techniques must be used [8]. Nevertheless, even if node-wise imposition of essential boundary conditions is achieved, this does not ensure
an “exact” interpolation, since it is necessary to preclude interior nodes to have influence on the boundary (this is straightforward in the FEM, since shape functions associated to interior nodes vanish along the boundary). This is not the case in one of these methods, the so-called Natural Element Method (NEM) [48]. This method is based on the use of natural neighbor interpolation schemes in the framework of a Galerkin procedure [22]. In [48] the interpolant character of this method along convex boundaries was demonstrated and therefore its ability to exactly interpolate (up to linear precision) essential boundary conditions. More recently, in [21] the ability of the method to exactly interpolate essential boundary condition also on non-convex boundaries was also demonstrated.

Another problem of meshless methods is related to the numerical integration errors due to the employ of non-polynomial approximation functions, together with Gauss-Legendre quadratures. However, considerable improvement has been done in the use of especial nodal integration techniques with application to moving least squared methods [16] and also to natural element techniques [27]. Also, an in-deep study of the behavior of these methods and some rigorous mathematical proofs of convergence can be found in [5].

It seems therefore interesting to study the application of meshless methods to the numerical simulations of forming processes. Since their accuracy is not affected by the nodal distribution they seem to be especially well suited to be employed in updated Lagrangian approaches to these problems. Little has been done, however, in the application of these methods in forming processes simulation. In [53] the application of Reproducing Kernel Particle Methods [31] to the simulation of various forming processes, including isothermal extrusion and sheet forming is discussed. In [15] an example of upsetting is also analyzed. Bonet [11] employed the incremental flow formulation to simulate upsetting problems with Corrected Particle Hydrodynamics methods [12]. More recently, the authors have employed natural element techniques to the simulation of foams forming processes [17] and also to the simulation of aluminium extrusion, casting [2] and forging [47].

In the subsequent sections we will review the basics of the natural element method and study its application to the simulation of a non-isothermal aluminium extrusion process, among other forming processes.

## 1.3 Scope of the Paper

In section 2 we will present the flow formulation applied to a Sellars-Tegart model [44] for the aluminium. This type of behaviour will be applied to the simulation of aluminium extrusion. In section 3 we will review the basics of the natural element method and its good characteristics for the simulations here presented. The coupled non-isothermal procedure employed is introduced in section 4, while in section 5 some representative numerical simulations are presented. Finally, the paper is completed with some concluding remarks.

# 2 A VISCOPLASTIC MODEL FOR THE ALUMINIUM BEHAVIOR DURING EXTRUSION

As mentioned before, forming processes involve large deformation and the permanent change of shape that occurs make the elastic strains to be negligible. Although an elastic recovery or spring-back exists, it has been often neglected a priori by many researchers as a first approximation. The obvious advantage of this assumption is that the material can then be modelled as a non-newtonian (visco-plastic) fluid. This assumption is known as the flow formulation in the forming processes community [57]. Thus, the essential variables of the problem will be velocities and pressures, instead of displacements and pressures.

The deviatoric stresses will be, under this assumption,

$$s = 2\mu d, \qquad (5)$$
being, as usual,

$$ \sigma = s - pI \qquad (6) $$

where $p = -tr(\sigma)/3$ and $I$ stands for the second-order identity tensor. Obviously, in the most general case, the parameter $\mu$ will depend on both the level of strain (and hence the non-linear character of the behavior) and the temperature, as stated in Eq. (1). To derive the expression of the parameter $\mu$ it is a common practice to write the strain rate tensor as emerging from a visco-plastic potential. Following Perzyna [43]

$$ d^{vp} = \dot{\gamma} \frac{\partial Y(\sigma, q)}{\partial \sigma}, \qquad (7) $$

where $Y$ is the viscoplastic potential —usually coincident with the plastic criterion as has been considered here—, $\dot{\gamma}$ is a scalar function given by

$$ \dot{\gamma} = \frac{\langle g(Y(\sigma, q)) \rangle}{\eta} \quad \text{with} \quad \langle x \rangle = \frac{x + |x|}{2}, \qquad (8) $$

$\langle g \rangle$ is a monotonic function that takes zero value only if $Y(\sigma, q) \le 0$, $\eta$ 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:

$$ Y(\sigma, q) = \bar{\sigma}(s) - \sigma_y(\bar{d}, T), \qquad (9) $$

where $\bar{\sigma} = \sqrt{\frac{3}{2}s : s} = \sqrt{3J_2}$ represents the *effective* stress and $\sigma_y$ represents the uniaxial yield stress. $\bar{d}$ is the only internal variable in this model and is sometimes called *effective strain rate*:

$$ \bar{d} = \sqrt{\frac{2}{3} \bar{d} : \bar{d}}. \qquad (10) $$

Among the most extended laws to account for the aluminium yield stress is the so-called hyperbolic sine or Sellars-Tegart law [44]

$$ \sigma_y(\bar{d}) = S_m \operatorname{arcsinh} \left[ \left[ \left( \frac{\bar{d}}{A} \right) e^{\frac{Q}{RT}} \right]^{\frac{1}{m}} \right]. \qquad (11) $$

$S_m$, $m$ and $A$ are material parameters. $A$ is a factor that depends on the magnesium and silicium matrix solute content (see [34] 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 Eq.(11) a initial strain rate, $\bar{d}_0$ so as to render a *modified* Sellars-Tegart law:

$$ \sigma_y(\bar{d}) = S_m \operatorname{arcsinh} \left[ \left( \frac{\bar{d}_1}{A} \right) e^{\frac{Q}{RT}} \right]^{\frac{1}{m}} \quad \text{with} \quad \bar{d}_1 = \max\{\bar{d}, \bar{d}_0\}. \qquad (12) $$
If we combine now the general form of the strain rate tensor given in Eq.(7), with Eq.(9), we arrive to

$$d = \dot{\gamma} \frac{3s}{2\bar{\sigma}}. \qquad (13)$$

It is immediate now, by combining Eq.(10) and the definition of effective stress, $\bar{\sigma}$, to prove that $\dot{\gamma}$ is precisely the effective strain rate:

$$\bar{d} = \frac{\dot{\gamma}}{\bar{\sigma}} \sqrt{\frac{3}{2} s : s} = \dot{\gamma}. \qquad (14)$$

On the other hand, and by following the Perzyna-like model employed in Eqs.(7) and (8) and taking $g(f) = f$, we arrive to a relationship between equivalent stress and equivalent strain rate:

$$\bar{d} = \dot{\gamma} = \frac{\langle \bar{\sigma} - \sigma_y \rangle}{\eta} \Rightarrow \bar{\sigma} = \eta \bar{d} + \sigma_y \text{ if } \bar{\sigma} \ge \sigma_y \qquad (15)$$

that, introduced in Eq.(13), accounting Eq.(14), gives the following visco-plastic constitutive equation:

$$s = 2 \frac{\eta \bar{d} + \sigma_y(\bar{d})}{3\bar{d}} d. \qquad (16)$$

Note that, depending on the $\eta$ value, the return to the yield surface is done with different velocity. Since it is common to describe aluminium behaviour as rigid-plastic (rather than viscoplastic) we employ null viscosity, so as to enforce $Y = \bar{\sigma} - \sigma_y = 0$, leading to

$$s = \frac{2\sigma_y}{3\bar{d}} d. \qquad (17)$$

Finally, the constitutive equation, accounting the incompressibility of plastic flow results:

$$\sigma = 2\mu d - pI, \text{ with } \mu = \frac{\sigma_y}{3\bar{d}}. \qquad (18)$$

Of course, this simple model has important limitations. Undoubtedly, the lack of elastic behaviour is one of the most important. Thus, spring-back can not be predicted. However, as mentioned before, it has rendered good results and seems to be widely accepted in the forming processes community [57] [59] [34].

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):

$$\mu = \mu_0 \bar{d}^{n-1} \qquad (19)$$

where $n$ is the so-called sensitivity index. The interested reader can consult [47] for some applications on this class of problems.

In the following section we briefly review the basics of the natural element method, a necessary step prior to introducing the NE discretization of the governing equations.
<figure><img src="image_3.png" /><figcaption>Figure 2. Delaunay triangulation and Voronoi diagram of a cloud of points. On the right, an example of a degenerate distribution of nodes, with the two possible triangulations depicted. In this last case, four points lie in the same circumcircle and thus no single triangulation exists</figcaption></figure>

# 3 THE NATURAL ELEMENT METHOD FOR FLOW PROBLEMS

## 3.1 Natural Neighbor Interpolation

The natural element method is one of the youngest meshless methods and has received much less attention than others. Briefly speaking, it is based on the use of any natural neighbor-based interpolation scheme to construct the approximation spaces of the Galerkin method. The interested reader can consult a wide overview of the recent advances in this method in [22].

The interpolation scheme used in the NEM is based on the geometrical concepts of Delaunay triangulations or their dual structure, the Voronoi diagram. Consider a cloud of irregularly distributed nodes $N = \{n_1, n_2, \dots, n_N\}$. This cloud is, in general, the only information needed to discretize the domain, as will be seen in the following section.

The Delaunay triangulation [24] of the set $N$, $\mathcal{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 2. The dual structure of a Delaunay triangulation is the so-called Voronoi diagram [52], 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,

$$T_I = \{ \mathbf{x} \in \mathbb{R}^n : d(\mathbf{x}, \mathbf{x}_I) < d(\mathbf{x}, \mathbf{x}_J) \ \forall J \neq I \}, \quad (20)$$

where $T_I$ is the Voronoi cell and $d(\cdot, \cdot)$ represents the Euclidean distance. $n$ represents the space dimension (two or three in the examples presented in this paper). 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 [50]. 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 [28] 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 [46]. 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:

$$T_{IJ} = \{ \mathbf{x} \in \mathbb{R}^n : d(\mathbf{x}, \mathbf{x}_I) < d(\mathbf{x}, \mathbf{x}_J) < d(\mathbf{x}, \mathbf{x}_K) \ \forall K \neq J; \forall K \neq I \}. \quad (21)$$
<figure><img src="image_4.png" /><figcaption>Figure 3. Definition of the Natural Neighbour coordinates of a point $x$</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4. Natural Element (Sibson) shape function (courtesy N. Sukumar)</figcaption></figure>

If a new point is added to a given cloud of points, the Voronoi cells will be modified by its presence. Sibson [45] 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 $\kappa(x)$ and $\kappa_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

$$ \phi_I^{sib}(x) = \frac{\kappa_I(x)}{\kappa(x)}. \qquad (22) $$

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 4.

Recently, some new interpolation schemes based on the concept of natural neighbors have been proposed [29]. 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,
<figure><img src="image_6.png" /><figcaption>Figure 5. Definition of non-Sibsonian coordinates</figcaption></figure>

instead of volumes). If we define the cell intersection $t_{IJ} = \{x \in T_I \cap T_J, J \neq I\}$ (note that $t_{IJ}$ may be an empty set) we can define the value

$$ \alpha_J(\mathbf{x}) = \frac{|t_{IJ}|}{d(\mathbf{x}, \mathbf{x}_J)}. \qquad (23) $$

Thus, the point $\mathbf{x}$ shape function value with respect to node 4 in Figure 5 is defined as

$$ \phi_4^{ns}(\mathbf{x}) = \frac{\alpha_4(\mathbf{x})}{\sum_{J=1}^{n} \alpha_J(\mathbf{x})} = \frac{s_4(\mathbf{x})/h_4(\mathbf{x})}{\sum_{J=1}^{n} [s_J(\mathbf{x})/h_J(\mathbf{x})]}, \qquad (24) $$

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, $\mathbf{x}$.

Derivatives of the Laplace shape function are not defined along the edges of the Delaunay triangles that lie inside its support (see [49]). For the purposes of the work here presented, Sibson's interpolation has been considered.

## 3.2 Properties of Natural Neighbour Interpolation

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

$$ \phi_I(\mathbf{x}_J) = \delta_{IJ} \qquad (25) $$

with $\delta_{IJ}$ 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 [46]):

$$ \sum_{I=1}^{n} \phi_{I}(\mathbf{x}) \mathbf{x}_{I} = \mathbf{x} \qquad (26) $$

One important property that derives from Eq.(25) 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 [21] by using the theory of α-shapes. This proof was later generalised in [20]. Idelsohn and co-workers later adopted this same approach based on α-shapes in the context of the *Meshless Finite Element Method* [30]. Another method to impose essential boundary conditions in the context of the NEM by employing visibility criteria was developed by Chinesta and co-workers [54].

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 [49] [19].

# 4 GOVERNING EQUATIONS AND NATURAL ELEMENT DISCRETIZATION

In this section we describe the numerical model of the extrusion process, as explained before, in the context of a natural neighbour Galerkin method. We firstly review the governing equations and the mixed approximation used to deal with the incompressibility condition. The technique used to extract the geometry of the free surface as it evolves is also explained.

## 4.1 Governing Equations

The equilibrium or momentum balance equations, considering that the material is subjected to inertial effects and gravity loads can be written as:

$$ \nabla \cdot \sigma + \mathbf{b} = \rho a \quad \text{in } \Omega \qquad (27) $$

where **b** stands for the gravity loads and **a** represent the acceleration. Together with these equations, appropriate boundary conditions are considered:

$$ \sigma \cdot n = \bar{t} \quad \text{on } \Gamma_t \qquad (28) $$

$$ \mathbf{v} = \bar{\mathbf{v}} \quad \text{on } \Gamma_v \qquad (29) $$

with $\Gamma_t \cup \Gamma_v = \partial\Omega$ and $\Gamma_t \cap \Gamma_v = \emptyset$.

In an updated Lagrangian framework, the acceleration can be obtained after differentiation of the velocity terms, i.e.

$$ a = \frac{dv}{dt} = \dot{v} \qquad (30) $$

whereas if eulerian coordinates were chosen, a material (total) derivative must be employed:

$$ a = \frac{Dv}{Dt} = \dot{v} + v \cdot \nabla v \qquad (31) $$

Again, the convective terms in Eq. (31) pose additional difficulties, since any stabilization technique must be used. In this work, as widely discussed before, an updated Lagrangian approach has been chosen, so this problem is avoided.

The incompressibility conditions of the flow, mentioned before, can be expressed in rate form as

$$ \nabla \cdot \mathbf{v} = 0, \qquad (32) $$
Aluminium extrusion processes are slow enough so as to take $\mathbf{a} = \mathbf{0}$ and, in practice, the vast majority of the simulations available in the literature do so. However, it might be interesting to take into account these terms when simulating other forming processes, such as the phenomena occurring at the injection chamber of a die casting machine, for instance [36]. If this is the case, natural element methods provide a simple means to discretize in time the resulting equations by employing the method of characteristics. By taking into account that the characteristic lines of Eq.(27) are precisely the nodal pathlines, the acceleration can be approximated as

$$ \frac{dv}{dt} \approx \frac{v^t(\mathbf{x}) - v^{t-\Delta t}(\mathbf{X})}{\Delta t}, \qquad (33) $$

where $t$ and $t - \Delta t$ represent the current and previous time steps, respectively. Similarly, $\mathbf{X}$ represents the position occupied by point $\mathbf{x}$ at time $t - \Delta t$.

The weak form of the problem can thus be expressed as

$$ \int_{\Omega} 2\mu \mathbf{d} : \mathbf{d}^* d\Omega - \int_{\Omega} p \mathbf{I} : \mathbf{d}^* d\Omega = - \int_{\Omega} \rho \mathbf{b} \mathbf{v}^* d\Omega + \int_{\Omega} \rho \frac{dv}{dt} \mathbf{v}^* d\Omega \qquad (34) $$

$$ \int_{\Omega} \nabla \cdot \mathbf{v} \, p^* d\Omega = 0 \qquad (35) $$

An example of the performance of such simple technique is presented in section 6 below.

In spite of the simple formulation presented before, for the numerical simulation of extrusion, inertial and gravity terms are often neglected, as mentioned before, and we will do so in the examples that follow. If this is the case, the equations presented before will read

$$ \nabla \cdot \sigma = 0, \qquad (36) $$

and

$$ \nabla \cdot \mathbf{v} = 0. \qquad (37) $$

This mechanical problem is coupled with the following heat transfer equations:

$$ \nabla(k\nabla T) + \dot{r} - (\rho c_p \dot{T}) = 0 \qquad (38) $$

where $k$ denotes the thermal conductivity, $\dot{r}$ the heat generation rate, $\rho$ the specific density and $c_p$ the specific heat of the metal. The rate of heat generation in the aluminium billet due to plastic deformation is given by

$$ \dot{r} = \beta \sigma : \mathbf{d} \qquad (39) $$

where $\beta$ represents the fraction of mechanical energy transformed to heat that has been assumed to be 0.9 [55].

If we write the incremental variational equation at time $t + \Delta t$ we arrive to:

$$ \int_{\Omega(t+\Delta t)} \left( -(p^t + \Delta p)\mathbf{I} + 2\mu(\mathbf{d}^t + \Delta \mathbf{d})(\mathbf{d}^t + \Delta \mathbf{d}) \right) : \mathbf{d}^* d\Omega = 0. \qquad (40) $$
Due to the non-linear character of the constitutive equations, an iterative approach has been applied, using the Newton-Raphson scheme, thus leading to

$$ \int_{\Omega(t+\Delta t)} \left( -\Delta\Delta p\mathbf{I} + 2\mu \left( \frac{\partial\mu(d_k^{t+\Delta t})}{\partial\mathbf{d}} : \Delta\Delta\mathbf{d} \right) d_k^{t+\Delta t} + \right. \\ \left. \qquad +2\mu(d_k^{t+\Delta t})\Delta\Delta\mathbf{d} \right) : \mathbf{d}^* d\Omega = \\ = -\int_{\Omega(t+\Delta t)} (-p_k^{t+\Delta t}\mathbf{I} + 2\mu(d_k^{t+\Delta t})d_k^{t+\Delta t}) : \mathbf{d}^* d\Omega, \qquad (41) $$

where the subscript k indicates the iteration within a time increment. The incremental form of the incompressibility condition results

$$ \int_{\Omega(t+\Delta t)} \nabla \cdot (\Delta\Delta\mathbf{v}) p^* d\Omega = - \int_{\Omega(t+\Delta t)} \nabla \cdot (\mathbf{v}_k^{t+\Delta t}) p^* d\Omega. \qquad (42) $$

If we approximate velocities and pressures, as well as their variations, by employing a finite-dimensional set of basis functions, we arrive to a discrete form of the previous equations (Bubnov-Galerkin method)

$$ \Delta\Delta\mathbf{v}^h(\mathbf{x}) = \sum_{I=1}^{nn} \phi_I(\mathbf{x})\Delta\Delta\mathbf{v}_I; \quad \mathbf{v}^*(\mathbf{x}) = \sum_{I=1}^{nn} \phi_I(\mathbf{x})\mathbf{v}_I^* \qquad (43) $$

$$ \Delta\Delta p^h(\mathbf{x}) = \sum_{I=1}^{nn} \psi_I(\mathbf{x})\Delta\Delta p_I; \quad p^*(\mathbf{x}) = \sum_{I=1}^{nn} \psi_I(\mathbf{x})p_I^*. \qquad (44) $$

where nn represents the number of nodes considered in the approximation. The functions $\psi_I(\mathbf{x})$ and $\phi_I(\mathbf{x})$ in this work represent some form of natural neighbour interpolation, which will be studied in the following section. This leads to the linear system:

$$ \left( \begin{array}{cc} \mathbf{K}_k^{t+\Delta t} & \mathbf{G} \\ \mathbf{G}^T & \mathbf{0} \end{array} \right) \left( \begin{array}{c} \Delta\Delta\mathbf{v} \\ \Delta\Delta p \end{array} \right) = \left( \begin{array}{c} \mathbf{f}_k^{t+\Delta t} \\ \mathbf{0} \end{array} \right) \qquad (45) $$

being the update given by

$$ \mathbf{v}_{k+1}^{t+\Delta t} = \mathbf{v}_k^{t+\Delta t} + \Delta\Delta\mathbf{v} \qquad (46) $$

$$ p_{k+1}^{t+\Delta t} = p_k^{t+\Delta t} + \Delta\Delta p \qquad (47) $$

Obviously, Eqs.(41) and (42) will be solved iteratively until convergence.

As it is well known, not all of the possible combinations of functions $\psi$ and $\phi$ give rise to stable formulations and some of them tend to lock. The development of stable, mixed, natural neighbor Galerkin formulations has been a matter of research until recently. This topic is reviewed in the following section.

## 4.2 Mixed Natural Neighbor Approximations

In the problem defined by Eqs.(36) and (37) there are two essential variables, namely velocity and pressure. Thus, we deal with a Hellinger-Reissner-like variational principle. It is well known that not all the possible combinations of velocity and pressure spaces render stable formulations. These must verify the so called inf-sup or Ladyzhenskaya-Babuška-Brezzi condition [42] [4] [13], together with the ellipticity condition.
The discrete form of the LBB condition may be written as:

$$ \inf_{p^h \in \mathcal{P}_h} \sup_{v^h \in \mathcal{U}_h} \frac{\int_{\Omega} p^h \nabla \cdot v^h d\Omega}{\|p^h\|_0 \|v^h\|_1} = \gamma_h \geq \gamma > 0 \quad (48) $$

where $\gamma$ is a positive constant independent of the mesh size, $h$. $\mathcal{U}_h$ and $\mathcal{P}_h$ represent, respectively, the displacement and pressure approximation spaces whose norms are defined as

$$ \| \cdot \|_0^2 = \int_{\Omega} (\cdot)^2 d\Omega \quad (49) $$

$$ \| \cdot \|_1^2 = \int_{\Omega} \sum_{i,j=1}^{3} \left( \frac{\partial(\cdot)_i}{\partial x_j} \right)^2 d\Omega \quad (50) $$

This condition have been proved analytically only for some of the simplest finite element approximations, but there exist some numerical test in the literature that give valuable information on the stability of the studied approximation (see [14] or [7], for instance). This test is only valid for the given geometry and problem considered, but the approximation will not verify the LBB condition if the test is not passed.

In order to evaluate expression (48), an equivalent form is developed:

$$ \inf_{w_h \in \mathcal{P}_h} \sup_{u_h \in \mathcal{U}_h} \frac{\mathbf{W}_h^T \mathbf{G}_h \mathbf{U}_h}{\sqrt{\mathbf{W}_h^T \mathbf{G}_h \mathbf{W}_h} \cdot \sqrt{\mathbf{U}_h^T \mathbf{S}_h \mathbf{U}_h}} = \gamma_h \geq \gamma > 0 \quad (51) $$

where $\mathbf{W}_h$ and $\mathbf{U}_h$ are vectors of the nodal displacements corresponding to $w_h$ and $u_h$. $\mathbf{G}_h$ and $\mathbf{S}_h$ are the matrices associated with the norms

$$ \|p_h\|_0^2 = \mathbf{W}_h^T \mathbf{G}_h \mathbf{W}_h \quad (52) $$

$$ \|u_h\|_1^2 = \mathbf{U}_h^T \mathbf{S}_h \mathbf{U}_h \quad (53) $$

$\mathbf{G}_h$ and $\mathbf{S}_h$ are positive semidefinite and positive definite [14], respectively. It is then demonstrated that the first non-zero eigenvalue, $\lambda_k$, of the problem

$$ \mathbf{G}_h \phi_h = \lambda \mathbf{S}_h \phi_h \quad (54) $$

is related to the searched value $\gamma_h$ through the expression

$$ \gamma_h = \sqrt{\lambda_k}. \quad (55) $$

The obvious problem when developing mixed natural neighbor approximations is that, as mentioned before, natural neighbor interpolants are, at most, linearly complete. In [28] a solution is proposed to this problem based on the use of the partition of unity property (see [48]) This means that natural neighbor interpolants can be enriched so as to be able to reproduce higher order polynomial fields by employing the method firstly proposed by Babuška [6].

Although Sibson-Thiessen interpolation for velocities and pressures, respectively, has been found not to verify the inf-sup condition (some spurious pressure modes can be found for some especial boundary conditions), it has rendered excellent results in previous works in the field of Fluid Mechanics [39] [38] and rarely presents locking.
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Cloud of 9 nodes on a regular lattice</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Plot of the shape functions $\phi x$ (a), $\phi y$ (b) and $\phi xy$ (c)</figcaption></figure>

Nevertheless, some enrichments have been studied in [28] in order to develop a stable formulation. It was found that the enrichment with the monomial $xy$ is necessary in order to reproduce bi-linear fields and also to render stable formulations. The resultant shape functions, for the central node in Figure 6, can be seen in Figures 7(a), (b) and (c).

In order to test whether the proposed formulations comply with the inf-sup condition,
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Geometry of the problem used to evaluate the numerical inf-sup test</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Inf-sup test for the proposed approximations</figcaption></figure>

we checked a sequence of three meshes composed of 3 × 3, 4 × 4 and 5 × 5 nodes, regularly and irregularly distributed over a square of side unity (see Figure 8). It is accepted that, if the numerically computed inf-sup value, $γ_h$, is bounded away from zero as the mesh size is increased, the test is passed.

Results for these discretisations are shown in Figure 9, including when the test was applied over the 3/1 triangular Finite Element. These last were obtained by restricting the natural neighbours at a given point to the three nodes that make up the Delaunay triangle [48]. The Sibson approximation, enriched with $\{1, x, y, xy\}$ for the displacements, gives good results using both Thiessen ($C^{-1}$) and Sibson interpolations for the pressure. However, the enrichment with monomials $\{1, x^2, y^2\}$, initially proposed to avoid rank deficiency in the resultant stiffness matrix, gives rise to spurious pressure modes. These same spurious modes can be obtained for certain configurations in the Sibson-Thiessen element, as in the bilinear displacement–constant pressure quadrilateral Finite Element [7].

Although the results for a unique problem cannot be extrapolated, the proposed formulation for the $\{1, x, y, xy\}$ enrichment seems to be stable and adequate for a wide variety of problems. However, its use forces to employ some special technique (see [38] and references
therein) to solve the resultant system of equations, since it leads to rank deficiency in the system of equations. The {1, xy} enrichment avoids this drawback, while maintaining the stability properties of the formulation.

In the numerical examples that follow, a Sibson-Thiessen formulation has been employed, which has one nodal degree of freedom less than the previous one and, as mentioned before, rarely presents locking.

## 4.3 The α-Shape Based Natural Element Method

For the complete description of the method, the free-surface tracking algorithm deserves some comments. It is based on the use of natural neighbour approximations constructed over α-shapes. The main idea is simple: can the cloud of points itself contain geometrical information about the shape of the domain? The concept of shape has no associated formal meaning, but Edelsbrunner [25] [26] developed the theory necessary to extract a finite set of shapes from a cloud of points based on the use of a scalar parameter α.

An α-shape is a polytope that is not necessarily convex nor connected, being triangulated by a subset of the Delaunay triangulation of the points. Thus, the empty circumcircle criterion holds. Let N be a finite set of points in $R^3$ and α a real number, with $0 \le \alpha < \infty$. A k-simplex $\sigma_T$ with $0 \le k \le 3$ is defined as the convex hull of a subset $T \subseteq N$ of size $|T| = k + 1$. Let b be an α-ball, that is, an open ball of radius α. A k-simplex $\sigma_T$ is said to be α-exposed if there exist an empty α-ball b with $T = \partial b \cap 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,\alpha}$ 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,\alpha}$, the edges in $F_{1,\alpha}$ and the vertices or nodes in $F_{0,\alpha}$.

A three-dimensional simplicial complex is a collection, $C$, of closed k-simplices ($0 \le k \le 3$) that satisfies:

(i) If $\sigma_T \in C$ then $\sigma_{T'} \in C$ for every $T' \subseteq T$.

(ii) The intersection of two simplexes in $C$ is empty or is a face of both.

Each k-simplex $\sigma_T$ included in the Delaunay triangulation, $D$, defines an open ball $b_T$ whose bounding spherical surface (in the general case) $\partial b_T$ passes through the $k+1$ points of the simplex. Let $\rho_T$ be the radius of that bounding sphere, then, the family $G_{k,\alpha}$, is formed by all the k-simplexes $\sigma_T \in D$ whose ball $b_T$ is empty and $\rho_T < \alpha$. The family $G_{k,\alpha}$ does not necessarily form simplicial complexes, so Edelsbrunner and Mücke [26] defined the α-complex, $C_\alpha$, as the simplicial complex whose k-simplexes are either in $G_{k,\alpha}$, or else they bound $(k+1)$-simplexes of $C_\alpha$. If we define the underlying space of $C_\alpha$, $|C_\alpha|$, as the union of all simplexes in $C_\alpha$, the following relationship between α-shapes and α-complexes is found:

$$S_{\alpha} = |C_{\alpha}| \quad \forall 0 \le \alpha < \infty \qquad (56)$$

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 10, 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
<figure><img src="image_11.png" /><figcaption>Figure 10. Two snapshots of a two-dimensional simulation of an extrusion process. Equivalent plastic strain rate is depicted</figcaption></figure>

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 11 some members of the family of α-shapes of the cloud of points in its final configuration (corresponding to Figure 10(b)) are depicted. In Figure 11(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 11(f)).

The use of α-shapes has also a direct outcome in the imposition of essential boundary conditions. As mentioned before, if natural neighbour (Sibson or Laplace) interpolants are computed over an α-shape of the cloud, linear interpolation of essential boundary conditions on convex as well as non-convex boundaries is guaranteed (see [21] and [20] for an in-deep demonstration of this property). We have coined this method the α-shape based NEM, or simply α-NEM [21].

This simple algorithm provides a means to easily avoid complex geometrical checks when dealing with free surface flows that can develop holes or merge as a consequence of the movement. This is the case, for instance, when simulating the extrusion of hollow profiles, where the flow must be divided before passing through the die.

## 4.4 Coupled Thermo-Mechanical Simulations

When a coupled thermo-mechanical simulation has been carried out, we have employed a semi-implicit algorithm combined with a fixed-point iteration that allows to de-couple both problems. This algorithm was firstly employed in [1].

Assuming velocities and temperatures at time t, $v^t$ and $T^t$ respectively to be known, we proceed as follows:

1. *Explicit updating of the material domain.* This process is performed moving all the nodes with the associated velocity:

$$x_I^{t+\Delta t} = x_I^t + v_I^t \Delta t$$
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Some members of the family of α-shapes of the cloud of points used in the extrusion example. (a) $S_0$ (the cloud of points) (b) $S_{0.3}$ (c) $S_{0.5}$ (d) $S_{1.0}$ (e) $S_{1.5}$ and (f) $S_\infty$ (the convex hull of the set)</figcaption></figure>
2. Velocity updating. Compute $v_0^{t+\Delta t}$ using Eq.(41) and (42) assuming the temperature $T^t$ unchanged (the temperature field affects the pseudo-viscosity parameter $\mu$ (Eq.(18)) through the expression of the yield stress (Eq.(12)).

3. Temperature updating. Update $T_0^{t+\Delta t}$ solving Eq.(38). The material derivative of the temperature field is discretised in a Lagrangian form, i. e.,

$$ \left. \frac{dT}{dt} \right|_I = \frac{T(\boldsymbol{x}_I^{t+\Delta t}) - T(\boldsymbol{x}_I^t)}{\Delta t} $$

The resulting thermal diffusion problem is integrated by a semi-implicit Crank-Nicolson scheme. At each step, the velocity field that defines the rate of heat generation (Eq.(39)) is assumed to remain unchanged and given by the computed kinematics $v_0^{t+\Delta t}$.

4. Velocity correction. We recompute step (2) assuming $T_0^{t+\Delta t}$ in the yield stress expression, leading to the new velocity field $v_1^{t+\Delta t}$.

5. Temperature correction. We recompute step (3) using $v_1^{t+\Delta t}$ in the rate of generation expression, leading to the field $T_1^{t+\Delta t}$.

6. Convergence check.

$$ \text{if } ||v_1^{t+\Delta t} - v_0^{t+\Delta t}|| < \varepsilon $$

and

$$ |T_1^{t+\Delta t} - T_0^{t+\Delta t}| < \varepsilon $$

then proceed to the next time step. Otherwise return to (4).

Note that, by using a fully explicit time discretisation, the previous algorithm becomes uncoupled at each time step. Also if, after convergence in step (6), the previous nodal update is recomputed, that is, if we return to step (1), then the fixed point algorithm is applied to solve the fully implicit coupled model. This simple algorithm has proven to converge in the examples here presented in at most three iterations per time step. Additional details can be found in [1].

## 5 CHECKING THE ACCURACY OF THE α-NEM APPROACH

The α-NEM accuracy and capabilities have been analyzed in some of our former works [21] [19] [27] [38]. However, the use of a Lagrangian expression of the material derivative in the Galerkin discretisation or within an operator splitting technique in the framework of the NEM have not been analyzed enough. By this reason, and prior to introduce the numerical results concerning aluminium forming processes (making special emphasis on forging and extrusion processes, see section 6) we focus in this section in the description of two problems involving, in the first case, inertia terms and, in the second case, a microstructure that evolves during the process. Both aspects might be introduced if a more accurate mechanical model of aluminium processes, involving an elastoviscoplastic behaviour and in which inertial effects are retained, is developed.
<figure><img src="image_13.png" /><figcaption>Figure 12. Results for the two-dimensional simulation of a water column collapse</figcaption></figure>

## 5.1 Breaking Dam Problem

In order to show the capabilities of the α-NEM to accurately simulate dynamic problems we study here the well-known problem of a breaking dam (or water column collapse) (Figure 12). This problem was firstly experimentally studied by Martin and Moyce [37]. Consider a rectangular column of water with 0.05715 m wide and 0.05715 m high. The left side consists of an initially solid wall. In the simulation, the upper and the right boundary of the water column is considered as the free surface boundary. Initially, the particles are in hydrostatic equilibrium and when the right wall (dam) is removed, the water column collapses under the influence of the gravity, acting downwards.

The model is composed by 3364 nodes. The acceleration of gravity is $9,81 \text{ m/s}^2$ and the density of the water is $1000 \text{ kg/m}^3$. The fluid is assumed to be purely newtonian. Its viscosity was set to $\mu = 0.1 \text{ Pa} \cdot \text{s}$.

In Figure 13 a comparison is shown between the obtained, numerical, results and the experimental ones. The horizontal position of the waterfront is depicted in a dimensionless graph. Taking the initial width of the collum, $a$, time is depicted as $t\sqrt{g/a}$, where $g$ represents the acceleration of gravity. The waterfront position is depicted as $z/a$, where $z$ represents the actual position of the flow. The great accuracy of the results is noteworthy.

In order to show the accuracy of the proposed free-surface tracking algorithm based on the use of α-shapes, in Figures 14 and 15 the absolute and relative errors in mass conservation are depicted. This error is below 3 % during the whole simulation.
<figure><img src="image_14.png" /><figcaption>**Figure 13.** Comparison between experimental and numerical results for the water column collapse problem</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 14.** Absolute error vs. time in the mass conservation for the water column collapse</figcaption></figure>

## 5.2 Extrusion of Short Fiber-Reinforced Thermoplastics

The updated Lagrangian approach presented before presents some interesting capabilities, not only for the simulation of aluminium extrusion, but also for other materials. In this section we describe its use for the simulation of extrusion processes using short fiber suspensions. In [18] a formulation is presented that takes into account both the dependence of the kinematics with the fiber orientation and the orientation induced by the flow kinematics. In [18] model reduction techniques are used in order to alleviate the high number of degrees of freedom involved in usual discretisations. In this paper, we will restrict ourselves to the description of the flow kinematics and its numerical simulation by employing α-NEM techniques only. The interested reader is submitted to the before mentioned reference for more details in the model reduction techniques.
<figure><img src="image_16.png" /><figcaption>Figure 15. Relative error vs. time in the mass conservation for the water column collapse</figcaption></figure>

### 5.2.1 Flow kinematics

The governing equations for a short fiber suspension include the equilibrium equations and mass conservation, already described in Eqs.(36) and (37). In addition, the constitutive equation for such a suspension is described by:

$$ \sigma = -p\boldsymbol{I} + 2\mu\{\boldsymbol{D} + N_p(\boldsymbol{a}_4 : \boldsymbol{D})\} \quad (57) $$

where $N_p$ is a parameter depending on the fiber concentration and the fiber aspect ratio. $\boldsymbol{a}_4$ is the fourth-order fiber orientation tensor, defined by:

$$ \boldsymbol{a}_4 = \int \rho \otimes \rho \otimes \rho \otimes \rho \psi(\rho) d\rho \quad (58) $$

where $\rho$ is the unit vector aligned with the fiber axis and $\psi(\rho)$ is the orientation distribution function.

It is also possible to define a second-order orientation tensor, $\boldsymbol{a}$, such that

$$ \boldsymbol{a} = \int \rho \otimes \rho \psi(\rho) d\rho \quad (59) $$

Considering spheroidal fibres immersed in a dilute suspension, the orientation evolution is described by the so-called Jeffery equation:

$$ \frac{d\rho}{dt} = \Omega\rho + k(D\rho - (D:(\rho\otimes\rho))\rho) \quad (60) $$

where $\Omega$ represents the vorticity tensor and $k$ is a constant depending on the fiber aspect ratio. The fiber orientation distribution is governed by the Fokker-Plank equation that in absence of diffusion terms results

$$ \frac{d\psi}{dt} + \frac{\partial}{\partial \rho} \left[ \psi \frac{d\rho}{dt} \right] = 0 \quad (61) $$

Taking into account Eqs.(59), (60) and (61) we arrive to a general expression for the evolution of the second-order orientation tensor:

$$ \frac{da}{dt} = \Omega a - a\Omega + k(Da + aD - 2(a_4 : D)) \quad (62) $$
<figure><img src="image_17.png" /><figcaption>Figure 16. Cavity geometry</figcaption></figure>

### 5.2.2 Coupling the α-NEM and a particle approach for the Fokker-Plank equation

We assume that at each node at the initial time a set of N particles or virtual fibers with homogeneous orientation exist. These orientations are defined by:

$$ \theta_i = i \times \frac{2\pi}{N}, \quad i \in [0, \dots, N-1] \qquad (63) $$

In two dimensions, the orientation vector can be defined by

$$ \rho = \begin{pmatrix} \cos \theta \\ \sin \theta \end{pmatrix} \qquad (64) $$

and, introducing it in the Jeffery's equation and after some simple algebraic manipulation we get:

$$ \begin{aligned} \dot{\theta} = & (-\sin \theta \cos \theta) \begin{pmatrix} 0 & \frac{1}{2} \left( \frac{\partial u}{\partial y} + \frac{\partial v}{\partial x} \right) \\ -\frac{1}{2} \left( \frac{\partial u}{\partial y} + \frac{\partial v}{\partial x} \right) & 0 \end{pmatrix} \left( \begin{pmatrix} \cos \theta \\ \sin \theta \end{pmatrix} + \right. \\ & \left. + k (-\sin \theta \cos \theta) \begin{pmatrix} \frac{\partial u}{\partial x} & \frac{1}{2} \left( \frac{\partial u}{\partial y} + \frac{\partial v}{\partial x} \right) \\ \frac{1}{2} \left( \frac{\partial u}{\partial y} + \frac{\partial v}{\partial x} \right) & \frac{\partial v}{\partial y} \end{pmatrix} \right) \left( \begin{pmatrix} \cos \theta \\ \sin \theta \end{pmatrix} \right) \end{aligned} \qquad (65) $$

and thus the dependence of the rotation velocity on: (i) the fiber orientation $\theta$; (ii) the gradient of velocities and (iii) the fiber aspect ratio, becomes explicit in Eq. (65).

Figure 16 shows a cavity whose filling process was simulated using the natural element method just described. The suspension viscosity is fixed to $10^3$ Pa.s; fibers were assumed to have an aspect ratio of $k = 0.8$, with a particle number $N_p = 10$. At time $t = 0$ the suspension is assumed to fill the inlet cavity (as shown in Figure 16) with an initial isotropic fiber orientation distribution and 36 virtual fibers per node. The use of higher numbers of virtual particles per node (more than 360, for instance) considerably increases the computational cost of this model and hence the use of model reduction techniques, as explained in [18], is justified.

In Figure 17 the influence of the flow kinematics on the fiber orientation field can be observed, while in Figure 18 the orientation and velocity fields and α-shape for the nineteenth time step of the simulation are depicted.

Again the ability of α-NEM techniques to accurately reproduce large transformations from an updated Lagrangian perspective is shown. We can also notice the accurate ability of this technique for accounting the advection problems (Eq.61) in a Lagrangian framework. However, for two-dimensional geometries good-quality mesh generators exist and the
<figure><img src="image_18.png" /><figcaption>Figure 17. Evolution of the orientation field</figcaption></figure>

ability of the FEM to accurately simulate this class of problems can be argued. In three-
dimensional problems, however, mesh generation is far from being totally automatized and
the advantages of such an updated Lagrangian approach could be more evident. This will
be analyzed in the following section.

6 ALUMINIUM FORMING PROCESSES: NUMERICAL EXAMPLES

In this section we analyze the application of the α-NEM to the simulation of some aluminium forming processes. Although the section is devoted mainly to aluminium extrusion, we introduce next an application to forging as well. In sections 6.2 and 6.3 two- and three-dimensional examples of aluminium extrusion are studied.
<figure><img src="image_19.png" /></figure>

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

<figure><img src="image_21.png" /><figcaption>**Figure 18.** Simulation of an extrusion flow. Velocity contour (a) and arrow (b) map and α-shape at an intermediate time step</figcaption></figure>

## 6.1 Forging of a Workpiece

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 Eq. (19). We have considered $\mu_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.
<figure><img src="image_22.png" /><figcaption>**Figure 19.** The piece to be forged (a) and schematic representation of the forging process (b)</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Figure 20.** Geometry and dimensions of the forged piece</figcaption></figure>

The geometry of the piece is shown in Figure 19(a). The simulation deals with the forging of the central region of the piece, justifying the assumption of plane strain (see Figure 19(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 20.

The equivalent plastic strain at time steps 96, 120, 150 and 168 are shown in Figure 21. Very accurate results were obtained in spite of the large strains and displacements involved in the simulation.

In Figure 22 the distribution of pressure at intermediate time steps is shown. It can be seen that some slight pressure oscillation exists. This is due to the lack of conformity with the LBB condition of the proposed formulation. This can be eliminated, as explained before, by employing an enriched NE formulation. It can be seen, however, how the proposed method is a valuable tool to predict the formation of the flash.
<figure><img src="image_24.png" /><figcaption>Figure 21. Equivalent strain rate distribution (s⁻¹) Throughout the forging process</figcaption></figure>

## 6.2 Two-Dimensional Simulation of Aluminium Extrusion

In this section we explore the ability of the α-NEM to simulate a simple, academic, two-dimensional example of aluminium extrusion. In this case, we simulate the extrusion of two rectangular slabs, 1mm thick and wide enough so as to accept the hypothesis of plane strain. The geometry of the device is shown in Figure 23.
<figure><img src="image_25.png" /><figcaption>**Figure 22.** Pressure ($\frac{N}{mm^2}$) distribution throughout the forging process</figcaption></figure>

Due to the large size of the billet, only a small region of about $15 \ mm^2$ was simulated (see detail A in Figure 23) Approximate velocities were imposed on the boundary of this region in order to simulate the flow of the whole billet. If the ram moves with constant velocity, also constant velocities can be imposed on the boundary of this small region with great accuracy [2]. To obtain a profile velocity of about 0.5m/min a velocity of 1 mm/s is imposed.
<figure><img src="image_26.png" /><figcaption>**Figure 23.** Extrusion of a flat profile</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>**Figure 24.** Boundary conditions for the two-dimensional aluminium extrusion problem</figcaption></figure>

A Sellars-Tegart law was assumed for an AA6063 aluminium alloy, whose properties are summarized in Table 1. In this first example a constant temperature of 773 K was assumed.

The cloud of points was formed by 2869 nodes, that remained the same throughout the simulation (i.e., no addition nor deletion of nodes was used). After 275 time steps of 0.0025s a profile of about 8mm was obtained. Finally, for the proper reproduction of the geometry, a value $\alpha = 0.02$ mm was taken. Different $\alpha$-shapes for the billet at different time steps are shown in Figure 25.

The velocity profile at the die outlet is shown in Figure 26. It is well-known that, for a proper extrusion, the velocity profile must remain as symmetric as possible. Otherwise, curvatures may appear on the resulting extrudate.

Finally, in Figure 27 the equivalent strain rate is depicted. These values are in agreement with those of Lof *et al* [32] by employing ALE techniques.
<figure><img src="image_28.png" /><figcaption>Figure 25. α-shapes at different time steps for the two-dimensional extrusion simulation</figcaption></figure>

## 6.3 Three-Dimensional Simulation of the Extrusion of a Cross-Shaped Profile

In order to show the capabilities of the technique presented before, we analyze now the simulation of a cross-shaped aluminium profile. The geometry of the domain is the same analyzed in [55]. In Figures 28 and 29 the geometry of the die and the profile are shown.

The diameter of the billet is 50 mm and it has been assumed that a length of 20 mm
<table><caption>**Table 1.** Material parameters for AA6063 Aluminium alloy (values taken from [34])</caption><thead><tr><th>Parameter</th><th>Units</th><th>Value</th></tr></thead><tbody><tr><td>d&#773;<sub>0</sub></td><td>(s<sup>-1</sup>)</td><td>0.005</td></tr><tr><td>S<sub>m</sub></td><td>(MPa)</td><td>25</td></tr><tr><td>m</td><td></td><td>5.4</td></tr><tr><td>A</td><td>(s<sup>-1</sup>)</td><td>6 &times; 10<sup>9</sup></td></tr><tr><td>Q</td><td>(<sup>J</sup>&frasl;<sub>mpt</sub>)</td><td>1.4 &times; 10<sup>5</sup></td></tr><tr><td>R</td><td>(<sup>J</sup>&frasl;<sub>molK</sub>)</td><td>8.314</td></tr></tbody></table>

<figure><img src="image_29.png" /><figcaption>**Figure 26.** Velocity profile at the die outlet ($\frac{mm}{s}$)</figcaption></figure>

will provide the necessary information so as to arrive to the stationary stage. The model consisted of 1313 nodes distributed over one height of the geometry, thus imposing symmetry boundary conditions at planes $X = 0$ and $X = Y$. The initial configuration of nodes is shown in Figure 30.

Nodes located on the upper side of the billet were forced to move with a speed of $2\frac{mm}{s}$, in order to obtain an exit velocity of $1\frac{m}{min}$ approximately. Initially, slipping boundary conditions were considered between the billet and the die and the container. The initial temperature was set to $723K$. The whole model was considered adiabatic, including the profile surface in contact with air. The aluminium characteristics are enumerated in Table 1. In this case the rigid-plastic Sellars-Tegart material model described in section 2 was used.

The simulation ran over 42 time steps of 0.025s. The obtained evolution for the equivalent strain rate is depicted in Figure 31.

In Figures 32 and 33 the evolution of the pressure is depicted. Again, despite the lack of conformity with the LBB condition of the Sibson-Thiessen formulation employed, good pressure distributions were obtained.

The third essential variable of the problem, temperature, is shown in Figures 34 and 35. As expected, the highest values were obtained where the strain, and thus the heat generation, was higher.

Finally, it is interesting to observe the evolution of the yield stress through the extrusion
<figure><img src="image_30.png" /><figcaption>Figure 27. Equivalent strain rate (s⁻¹) at different time steps</figcaption></figure>

process. This is done in Figure 36. As expected, the yield stress is lower where the temperature is higher, according to the Sellars-Tegart model presented in section 2.

It must be noted that the obtained temperature values are lower than those expected in such an extrusion process. This is due to the slipping boundary conditions imposed. To quantify the relative influence of these boundary conditions the same simulation can be done by imposing sticking boundary conditions at the die and the container. The influence on
<figure><img src="image_31.png" /><figcaption>**Figure 28.** Geometry of the container and the die for the cross-shaped profile problem</figcaption></figure>

<figure><img src="image_32.png" /><figcaption>**Figure 29.** Dimensions of the die (mm). The simulated region is highlighted in dashed line</figcaption></figure>

<figure><img src="image_33.png" /><figcaption>**Figure 30.** Initial configuration of the cloud of points for the cross-shaped profile problem</figcaption></figure>
<figure><img src="image_34.png" /><figcaption>**Figure 31.** Equivalent strain rate ($s^{-1}$) for the extrusion simulation</figcaption></figure>

the velocity profile is shown in Figure 37. The appearance of a “dead” or recirculating zone is now evident. This has also a direct outcome, as mentioned before, in the temperature profile (see Figure 38).

# 7 CONCLUDING REMARKS

In this paper, a review of some of the most recent advances in the numerical simulation of extrusion and other related forming processes by meshless methods has been made.
<figure><img src="image_35.png" /><figcaption>Figure 32. Pressure ($\frac{N}{mm^2}$) at the 20th time step</figcaption></figure>

<figure><img src="image_36.png" /><figcaption>Figure 33. Pressure ($\frac{N}{mm^2}$) at the 40th time step</figcaption></figure>

Meshless methods allow the use of an updated Lagrangian approach to this kind of non-stationary problems and seem to be an appealing alternative to other formulations, like ALE finite element techniques or the use of continuous remeshing in Lagrangian FEM.

A study is made on the use of the natural element method for this purpose. This meshless method seems to be very well suited to problems where extremely large deformations occur. An algorithm based on the use of $\alpha$-shapes is used that allows for a very accurate extraction of the geometry of the domain, without the need of any geometrical check
<figure><img src="image_37.png" /><figcaption>Figure 34. Temperature distribution (K) at the 20th time step</figcaption></figure>

<figure><img src="image_38.png" /><figcaption>Figure 35. Temperature distribution (K) at the 40th time step</figcaption></figure>

regarding inter-penetration of free surfaces or development of holes in the interior of the domain.

A coupled algorithm has been presented that allows to take into account the influence of the temperature field during extrusion on the yield stress. Other forming processes, like forging, or the extrusion of other materials, like short fiber reinforced suspensions, have also been studied. Also, when dynamic effects are important, the above presented formulation allows for a simple and efficient treatment of these terms, through the use of the method of characteristics.
<figure><img src="image_39.png" /><figcaption>**Figure 36.** Yield stress ($\frac{N}{mm^2}$) through the extrusion process</figcaption></figure>

<figure><img src="image_40.png" /><figcaption>**Figure 37.** Comparison of the obtained velocities with and without slipping boundary conditions ($\frac{mm}{s}$)</figcaption></figure>
<figure><img src="image_41.png" /><figcaption>**Figure 38.** Comparison of the temperature (K) profiles for the extrusion simulation</figcaption></figure>

Of course, these methods need further development. For instance, the possibility of
developing a LBB-compliant formulation without an increase in the number of unknowns
seems to be an interesting field of research.

ACKNOWLEDGEMENT

This work has been partially supported by the Spanish Ministry of Science and Technology,
through project number CICYT-DPI2002-01986.

REFERENCES

1. I. Alfaro, D. Bel, E. Cueto, M. Doblaré and F. Chinesta (2004). Three-dimensional simulation of Aluminium extrusion by the α-shape based Natural Element Method. *Computer Methods in Applied Mechanics and Engineering*, submitted.

2. I. Alfaro, J. Yvonnet, E. Cueto, F. Chinesta and M. Doblaré (2004). Meshless methods with application to metal forming. *Computer Methods in Applied Mechanics and Engineering*, in press.

3. E.H. Atzema and J. Huétink (1995). Finite Element Analysis of Forward/Backward Extrusion using ALE Techniques. In S.F. Shen and P.R. Dawson, editors, *proceedings of the Numiform 95 conference*, pages 383–388.

4. I. Babuška (1973). The Finite Element Method with Lagrange multipliers. *Numerische Mathematik*, **20**, 179–192.

5. I. Babuška, U. Banerjee and J.E. Osborn (2003). Survey of meshless and generalized finite element methods: A unified approach. *Acta Numerica*, **12**, 1–125.

6. I. Babuška and J.M. Melenk (1996). The partition of unity finite element method: Basic theory and applications. *Comp. Meth. in Appl. Mech. and Eng.*, **4**, 289–314.

7. K.J. Bathe (1996). *Finite Element procedures*. Prentice-Hall.
8 T. Belytschko, Y. Krongauz, D. Organ, M. Fleming and P. Krysl (1998). Meshless methods: An overview and recent developments. *Computer Methods in Applied Mechanics and Engineering*, **139**, 3–47.

9 T. Belytschko, Y.Y. Lu and L. Gu (1994). Element-Free Galerkin Methods. *International Journal for Numerical Methods in Engineering*, **37**, 229–256.

10 J. Bonet (1994). The incremental flow formulation for the numerical analysis of plane stress and thin sheet forming processes. *Computer Methods in Applied Mechanics and Engineering*, **114**, 103–122.

11 J. Bonet and S. Kulasegaram (2000). Correction and stabilization of smooth particle hydro-dynamics methods with applications in metal forming simulations. *International Journal for Numerical Methods in Engineering*, **47**, 1189–1214.

12 J. Bonet and T.S.L. Lok (1999). Variational and momentum preservation aspects of smooth particle hydrodynamic formulations. *Computer Methods in Applied Mechanics and Engineering*, **180**, 97–115.

13 F. Brezzi (1974). On the existence, uniqueness and approximation of saddle-point problems arising from lagrange multipliers. *Revue Française d'Automatique Informatique Recherche Opérationelle, Analyse Numérique*, **8**, 129–151.

14 D. Chapelle and K.J. Bathe (1993). The inf-sup test. *Computers and Structures*, **47**(4-5), 537–545.

15 J.-S. Chen, C. Pan, C.M. O.L. Roque and H.-P. Wang (1998). A lagrangian reproducing kernel particle method for metal forming analysis. *Computational Mechanics*, **22**, 289–307.

16 Jiun-Shyan Chen, Cheng-Tang Wu, Sangpil Yoon and Yang You (2001). A stabilized conforming nodal integration for galerkin mesh-free methods. *International Journal for Numerical Methods in Engineering*, **50**, 435–466.

17 F. Chinesta, E. Cueto, P. Quintela and J. Paredes (2004). Induced anisotropy in foams forming processes: modelling and simulation. *Journal of Materials Processing Technology*, in press.

18 F. Chinesta, E. Cueto, D. Ryckelynck and A. Ammar (2004). α-NEM and model reduction: two new and powerful numerical techniques to describe flows involving short fibers suspensions. *Revue Europeenne de Elements Finis*, submitted.

19 E. Cueto, B. Calvo and M. Doblaré (2002). Modeling three-dimensional piece-wise homogeneous domains using the α-shape based Natural Element Method. *International Journal for Numerical Methods in Engineering*, **54**, 871–897.

20 E. Cueto, J. Cegoñino, B. Calvo and M. Doblaré (2003). On the imposition of essential boundary conditions in Natural Neighbour Galerkin methods. *Communications in Numerical Methods in Engineering*, **19**(5), 361–376.

21 E. Cueto, M. Doblaré and L. Gracia (2000). Imposing essential boundary conditions in the Natural Element Method by means of density-scaled α-shapes. *International Journal for Numerical Methods in Engineering*, **49-4**, 519–546.

22 E. Cueto, N. Sukumar, B. Calvo, M.A. Martínez, J. Cegoñino and M. Doblaré (2003). Overview and recent advances in Natural Neighbour Galerkin methods. *Archives of Computational Methods in Engineering*, **10**(4), 307–384.

23 S. De and K. J. Bathe (2001). Towards an efficient meshless computational technique: the method of finite spheres. *Engineering Computations*, **18**, 170–192.

24 B. Delaunay (1934). Sur la Sphère Vide. A la mémoire de Georges Voronoi. *Izvestia Akademii Nauk SSSR, Otdelenie Matematicheskii i Estestvennyka Nauk*, **7**, 793–800.
25 H. Edelsbrunner, D.G. Kirkpatrick and R. Seidel (1983). On the shape of a set of points in the plane. *IEEE Transactions on Information Theory*, **IT-29**(4), 551–559.

26 H. Edelsbrunner and E. Mücke (1994). Three dimensional alpha shapes. *ACM Transactions on Graphics*, **13**, 43–72.

27 D. Gonzalez, E. Cueto, M.A. Martinez and M. Doblare (2004). Numerical integration in Natural Neighbour Galerkin methods. *International Journal for Numerical Methods in Engineering*, **60**(12), 2077–2104.

28 D. González, E. Cueto and M. Doblaré (2004). Volumetric locking in Natural Neighbour Galerkin methods. *International Journal for Numerical Methods in Engineering*, in press.

29 H. Hiyoshi and K. Sugihara (1999). Two generalizations of an interpolant based on Voronoi diagrams. *International Journal of Shape Modeling*, **5**(2), 219–231.

30 S.R. Idelsohn, E. Oñate, N. Calvo and F. Del Pin (2003). The meshless finite element method. *International Journal for Numerical Methods in Engineering*, **58**, 893–912.

31 W.K. Liu, S. Jun, S. Li, J. Adee and T. Belytschko (1995). Reproducing kernel particle methods. *International Journal for Numerical Methods in Engineering*, **38**, 1655–1679.

32 J. Lof (2000). *Developments in finite element simulations of aluminium extrusion*. PhD thesis, Department of Mechanical Engineering, Applied Mechanics Section. University of Twente, P.O. Box 217, 7500 AE Enschede, The Netherlands.

33 J. Lof (2001). Elasto-viscoplastic fem simulation of the aluminium flow in the bearing area for extrusion of thin-walled sections. *Journal of Materials Processing Technology*, **114**, 174–183.

34 J. Lof and Y. Blokhuis (2002). FEM simulations of the extrusion of complex thin-walled aluminium sections. *Journal of Materials Processing Technology*, **122**, 344–354.

35 J. Lof and J. Huétink (1999). Numerical simulation of the aluminium extrusion process. In J.A. Covas, editor, *Proceedings of the 2nd ESAFORM Conference on Material Forming*, pages 29–32.

36 J. López, F. Faura, J. Hernández and P. Gómez (2003). On the critical plunger speed and three-dimensional effects in high-pressure die casting injection chambers. *Journal of Manufacturing Science and Engineering. Transactions of the ASME*, **125**, 529–537.

37 J.C. Martin and W.J. Moyce (1952). An experimental study of the collapse of liquid columns on a rigid horizontal plane. *Philosophical Transactions of the Royal Society of London*, **244**, 312.

38 M.A. Martinez, E. Cueto, I. Alfaro, M. Doblare and F. Chinesta (2004). Updated Lagrangian free surface flow simulations with Natural Neighbour Galerkin methods. *International Journal for Numerical Methods in Engineering*, **60**(13), 2105–2129.

39 M.A. Martínez, E. Cueto, M. Doblaré and F. Chinesta (2003). Fixed mesh and meshfree techniques in the numerical simulation of injection processes involving short fiber suspensions. *Journal of Non-Newtonian Fluid Mechanics*, **115**, 51–78.

40 H.G. Mooi and J. Huétink (1995). Simulation of complex aluminium extrusion using an arbitrary eulerian lagrangian formulation. In S.F. Shen and P.R. Dawson, editors, *proceedings of the Numiform 95 conference*, pages 869–874.

41 B. Nayroles, G. Touzot and P. Villon (1992). Generalizing the finite element method: Diffuse approximation and diffuse elements. *Computational Mechanics*, **10**, 307–318.

42 O.A. Ladyzhenskaya (1969). *The Mathematical Theory of Viscous Incompressible Flows*. Gordon and Breach, London.

43 P. Perzyna (1966). Fundamental problems in visco-plasticity. In *Recent advances in Applied Mechanics*. Academic Press, New York.
44 C. M. Sellars and W.J.M. Tegart (1972). Hot workability. *International Metallurgical Review*, **17**, 1-24.

45 R. Sibson (1980). A Vector Identity for the Dirichlet Tesselation. *Mathematical Proceedings of the Cambridge Philosophical Society*, **87**, 151-155.

46 R. Sibson (1981). A brief description of natural neighbour interpolation. In *Interpreting Multi-variate Data*. V. Barnett (Editor), pages 21-36. John Wiley.

47 N. Sukumar, J. Dolbow, A. Devan, J.Yvonnet, F. Chinesta, D. Ryckelynck, P. Lorong, I. Alfaro, M.A. Martínez, E. Cueto and M. Doblaré (2004). Meshless Methods and Partition of Unity Finite Elements. *International Journal of Forming Processes*, to appear.

48 N. Sukumar, B. Moran and T. Belytschko (1998). The Natural Element Method in Solid Mechanics. *International Journal for Numerical Methods in Engineering*, **43**(5), 839-887.

49 N. Sukumar, B. Moran, A. Yu Semenov and V.V. Belikov (2001). Natural Neighbor Galerkin Methods. *International Journal for Numerical Methods in Engineering*, **50**(1), 1-27.

50 A. H. Thiessen (1911). Precipitation averages for large areas. *Monthly Weather Report*, **39**, 1082-1084.

51 J. van de Langkruis, J. Lof, W.H. Kool, S. van der Zwaag and J. Huétink (2000). Comparison of experimental AA6063 extrusion trials to 3D numerical simulations, using a general solute-dependent constitutive model. *Computational Materials Science*, pages 381-392.

52 G. M. Voronoi (1908). Nouvelles Applications des Paramètres Continus à la Théorie des Formes Quadratiques. Deuxième Memoire: Recherches sur les parallélloèdres Primitifs. *J. Reine Angew. Math.*, **134**, 198-287.

53 Sangpil Yoon and Jiun-Shyan Chen (2002). Accelerated meshfree method for metal forming simulation. *Finite Elements in Analysis and Design*, **38**, 937-948.

54 J. Yvonnet, D. Ryckelynck, P. Lorong and F. Chinesta (2004). A new extension of the Natural Element method for non-convex and discontinuous problems: the Constrained Natural Element method. *International Journal for Numerical Methods in Engineering*, **60**(8), 1452-1474.

55 J. Zhou, L. Li and J. Duszczyk (2003). 3D FEM simulation of the whole cycle of aluminium extrusion throughout the transient state and the steady state using the updated Lagrangian approach. *Journal of Materials Processing Technology*, **134**, 383-397.

56 O.C. Zienkiewicz (1984). Flow formulation for numerical solution of forming processes. In R.D. Wood J.F. T. Pittman, O.C. Zienkiewicz and J.M. Alexander (Eds.), *Numerical analysis of forming processes*, pages 1-44. John Wiley & Sons.

57 O.C. Zienkiewicz and P.N. Godbolet (1974). Flow of plastic and visco-plastic solids with special reference to extrusion and forming processes. *International Journal for Numerical Methods in Engineering*, **8**, 3-16.

58 O.C. Zienkiewicz, E. Oñate and J.C. Heinrich (1978). Plastic flow in metal forming. (I) Coupled thermal (II) Thin sheet forming. In *Applications of Numerical Methods to Forming Processes*. *AMD-vol 28*, pages 107-120.

59 O.C. Zienkiewicz, P.C. Pain and E. Oñate (1978). Flow of solids during forming and extrusion: some aspects of numerical solutions. *International Journal of Solids and Structures*, **14**, 15-38.