I Alfaro

D González

D Bel

E Cueto
email: ecueto@unizar.es

M Doblaré

F Chinesta

Recent Advances in the Meshless Simulation of Aluminium Extrusion and other Related Forming Processes

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.

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.

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 [START_REF] Zienkiewicz | Flow formulation for numerical solution of forming processes[END_REF] 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 [START_REF] Zienkiewicz | Flow of solids during forming and extrusion: some aspects of numerical solutions[END_REF] [57] [START_REF] Zienkiewicz | Plastic flow in metal forming. (I) Coupled thermal (II) Thin sheet forming[END_REF].

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 [START_REF] Zienkiewicz | Flow of plastic and visco-plastic solids with special reference to extrusion and forming processes[END_REF]. 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 σ = D(d, T ) • d, [START_REF] Alfaro | Three-dimensional simulation of Aluminium extrusion by the α-shape based Natural Element Method[END_REF] where d represents the strain rate tensor (symmetric part of the velocity gradient) and T the temperature. More recently, Bonet [START_REF] Bonet | The incremental flow formulation for the numerical analysis of plane stress and thin sheet forming processes[END_REF] 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 divergencefree condition is imposed on displacements theselves rather than in incremental form, typical of fluids.

Most researchers assumed (see [START_REF] Zienkiewicz | Flow formulation for numerical solution of forming processes[END_REF], 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 [START_REF] Zhou | 3D FEM simulation of the whole cycle of aluminium extrusion throughout the transient state and the steady state using the updated Lagrangian approach[END_REF]. 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, [START_REF] Lof | Elasto-viscoplastic fem simulation of the aluminium flow in the bearing area for extrusion of thin-walled sections[END_REF]).

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 [START_REF] Zienkiewicz | Flow formulation for numerical solution of forming processes[END_REF]. 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:

ρc p dT dt + ∇(k∇T ) + ṙ = 0 ( 2
)
where k denotes the thermal conductivity, ṙ 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

ρc p ∂T ∂t + v∇T + ∇(k∇T ) + ṙ = 0. (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 [START_REF] Alfaro | Meshless methods with application to metal forming[END_REF] now looks like

ρc p ∂T ∂t + (v -v)∇T + ∇(k∇T ) + ṙ = 0, (4)
where v represents this prescribed velocity. Clearly if v = v we are back to a Lagrangian formulation, whereas if 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 [START_REF] Atzema | Finite Element Analysis of Forward/Backward Extrusion using ALE Techniques[END_REF] [40] [35] [51].

Other recent contributions on this topic can be found in [START_REF] Lof | FEM simulations of the extrusion of complex thin-walled aluminium sections[END_REF] [33] [START_REF] Lof | Developments in finite element simulations of aluminium extrusion[END_REF]. Unfortunately, the prescribed velocity value v is not easy to determine and must be carefully chosen. In any case, the convective term, although is not so important, still remains.

Meshless Methods

In the last decade a new possibility has opened with the irruption of meshless methods [START_REF] Nayroles | Generalizing the finite element method: Diffuse approximation and diffuse elements[END_REF] [9] [START_REF] Babuška | The partition of unity finite element method: Basic theory and applications[END_REF] [48] [START_REF] De | Towards an efficient meshless computational technique: the method of finite spheres[END_REF]. 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 [START_REF] Belytschko | Element-Free Galerkin Methods[END_REF], modified variational principles or other techniques must be used [START_REF] Belytschko | Meshless methods: An overview and recent developments[END_REF]. 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) [START_REF] Sukumar | The Natural Element Method in Solid Mechanics[END_REF]. This method is based on the use of natural neighbor interpolation schemes in the framework of a Galerkin procedure [START_REF] Cueto | Overview and recent advances in Natural Neighbour Galerkin methods[END_REF]. In [START_REF] Sukumar | The Natural Element Method in Solid Mechanics[END_REF] 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 [START_REF] Cueto | Imposing essential boundary conditions in the Natural Element Method by means of density-scaled α-shapes[END_REF] 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 [START_REF] Chen | A stabilized conforming nodal integration for galerkin mesh-free methods[END_REF] and also to natural element techniques [START_REF] Gonzalez | Numerical integration in Natural Neighbour Galerkin methods[END_REF]. Also, an in-deep study of the behavior of these methods and some rigorous mathematical proofs of convergence can be found in [START_REF] Babuška | Survey of meshless and generalized finite element methods: A unified approach[END_REF].

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 [START_REF] Yoon | Accelerated meshfree method for metal forming simulation[END_REF] the application of Reproducing Kernel Particle Methods [START_REF] Liu | Reproducing kernel particle methods[END_REF] to the simulation of various forming processes, including isothermal extrusion and sheet forming is discussed. In [START_REF] Chen | A lagrangian reproducing kernel particle method for metal forming analysis[END_REF] an example of upsetting is also analyzed. Bonet [START_REF] Bonet | Correction and stabilization of smooth particle hydrodynamics methods with applications in metal forming simulations[END_REF] employed the incremental flow formulation to simulate upsetting problems with Corrected Particle Hydrodynamics methods [START_REF] Bonet | Variational and momentum preservation aspects of smooth particle hydrodynamic formulations[END_REF]. More recently, the authors have employed natural element techniques to the simulation of foams forming processes [START_REF] Chinesta | Induced anisotropy in foams forming processes: modelling and simulation[END_REF] and also to the simulation of aluminium extrusion, casting [START_REF] Alfaro | Meshless methods with application to metal forming[END_REF] and forging [START_REF] Sukumar | Meshless Methods and Partition of Unity Finite Elements[END_REF].

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.

Scope of the Paper

In section 2 we will present the flow formulation applied to a Sellars-Tegart model [START_REF] Sellars | Hot workability[END_REF] 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.

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 [START_REF] Zienkiewicz | Flow of plastic and visco-plastic solids with special reference to extrusion and forming processes[END_REF]. 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μd, (5)
being, as usual,

σ = s -pI (6)
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, as stated in Eq. (1). 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 [START_REF] Perzyna | Fundamental problems in visco-plasticity[END_REF]

d vp = γ ∂Y (σ, q) ∂σ , ( 7
)
where Y is the viscoplastic potential -usually coincident with the plastic criterion as has been considered here-, γ is a scalar function given by

γ = g(Y (σ, q)) η with x = x + |x| 2 , ( 8
)
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:

Y (σ, q) = σ(s) -σ y (d, T ), (9)
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:

d = 2 3 d : d. ( 10
)
Among the most extended laws to account for the aluminium yield stress is the so-called hyperbolic sine or Sellars-Tegart law [START_REF] Sellars | Hot workability[END_REF]

σ y (d) = S m arcsinh d A e Q RT 1 m . ( 11
)
S m , m and A are material parameters. A is a factor that depends on the magnesium and silicium matrix solute content (see [START_REF] Lof | FEM simulations of the extrusion of complex thin-walled aluminium sections[END_REF] 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, d 0 so as to render a modified Sellars-Tegart law:

σ y (d) = S m arcsinh d 1 A e Q RT 1 m with d 1 = max{d, d 0 }. ( 12
)
If we combine now the general form of the strain rate tensor given in Eq.( 7), with Eq.( 9), we arrive to

d = γ 3s 2σ . ( 13
)
It is immediate now, by combining Eq.( 10) and the definition of effective stress, σ, to prove that γ is precisely the effective strain rate:

d = γ σ 3 2 s : s = γ. ( 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:

d = γ = σ -σ y η ⇒ σ = ηd + σ y if σ ≥ σ y ( 15
)
that, introduced in Eq.( 13), accounting Eq.( 14), gives the following visco-plastic constitutive equation:

s = 2 ηd + σ y (d) 3d d. ( 16
)
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 rigid-plastic (rather than viscoplastic) we employ null viscosity, so as to enforce Y = σσ y = 0, leading to

s = 2σ y 3d d. ( 17
)
Finally, the constitutive equation, accounting the incompressibility of plastic flow results:

σ = 2μd -pI, with μ = σ y 3d . ( 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 [START_REF] Zienkiewicz | Flow of plastic and visco-plastic solids with special reference to extrusion and forming processes[END_REF] [59] [START_REF] Lof | FEM simulations of the extrusion of complex thin-walled aluminium sections[END_REF].

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

μ = μ 0 d n-1 ( 19
)
where n is the so-called sensitivity index. The interested reader can consult [START_REF] Sukumar | Meshless Methods and Partition of Unity Finite Elements[END_REF] 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.

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 [START_REF] Cueto | Overview and recent advances in Natural Neighbour Galerkin methods[END_REF].

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 , . . . , 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 [START_REF] Delaunay | Sur la Sphère Vide. A la memoire de Georges Voronoi[END_REF] 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 2. The dual structure of a Delaunay triangulation is the so-called Voronoi diagram [START_REF] Voronoi | Nouvelles Applications des Paramètres Continus à la Théorie des Formes Quadratiques. Deuxième Memoire: Recherches sur les parallélloèdres Primitifs[END_REF], 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 = {x ∈ R n : d(x, x I ) < d(x, x J ) ∀ J = I}, ( 20
)
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 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 [START_REF] Thiessen | Precipitation averages for large areas[END_REF]. 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 [START_REF] González | Volumetric locking in Natural Neighbour Galerkin methods[END_REF] 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 [START_REF] Sibson | A brief description of natural neighbour interpolation[END_REF]. 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 = {x ∈ R n : d(x, x I ) < d(x, x J ) < d(x, x K ) ∀ K = J; ∀ K = I}. ( 21
)
x 1 If a new point is added to a given cloud of points, the Voronoi cells will be modified by its presence. Sibson [START_REF] Sibson | A Vector Identity for the Dirichlet Tesselation[END_REF] 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

φ sib I (x) = κ I (x) κ(x) . ( 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 [START_REF] Hiyoshi | Two generalizations of an interpolant based on Voronoi diagrams[END_REF]. 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 IJ = {x ∈ T I T J , J = I} (note that t IJ may be an empty set) we can define the value

α J (x) = |t IJ | d(x, x J ) . ( 23
)
Thus, the point x shape function value with respect to node 4 in Figure 5 is defined as

φ ns 4 (x) = α 4 (x) n J=1 α J (x) = s 4 (x)/h 4 (x) n J=1 s J (x)/h J (x) , ( 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, x. Derivatives of the Laplace shape function are not defined along the edges of the Delaunay triangles that lie inside its support (see [START_REF] Sukumar | Natural Neighbor Galerkin Methods[END_REF]). For the purposes of the work here presented, Sibson's interpolation has been considered.

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

φ I (x J ) = δ IJ ( 25
)
with δ 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 [START_REF] Sibson | A brief description of natural neighbour interpolation[END_REF]):

n I=1 φ I (x)x I = x ( 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 [START_REF] Cueto | Imposing essential boundary conditions in the Natural Element Method by means of density-scaled α-shapes[END_REF] by using the theory of α-shapes. This proof was later generalised in [START_REF] Cueto | On the imposition of essential boundary conditions in Natural Neighbour Galerkin methods[END_REF]. Idelsohn and co-workers later adopted this same approach based on α-shapes in the context of the Meshless Finite Element Method [START_REF] Idelsohn | The meshless finite element method[END_REF]. Another method to impose essential boundary conditions in the context of the NEM by employing visibility criteria was developed by Chinesta and co-workers [START_REF] Yvonnet | A new extension of the Natural Element method for non-convex and discontnuous problems: the Constrained Natural Element method[END_REF].

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 [START_REF] Sukumar | Natural Neighbor Galerkin Methods[END_REF] [START_REF] Cueto | Modeling three-dimensional piece-wise homogeneous domains using the α-shape based Natural Element Method[END_REF].

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.

Governing Equations

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

∇ • σ + b = ρa in Ω ( 27
)
where b stands for the gravity loads and a represent the acceleration. Together with these equations, appropriate boundary conditions are considered:

σ • n = t on Γ t (28) v = v on Γ v (29) with Γ t ∪ Γ v = ∂Ω and Γ t ∩ Γ v = ∅.
In an updated Lagrangian framework, the acceleration can be obtained after differentiation of the velocity terms, i.e.

a = dv dt = v ( 30
)
whereas if eulerian coordinates were chosen, a material (total) derivative must be employed:

a = Dv Dt = v + v • ∇v (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

∇ • v = 0, (32)
Aluminium extrusion processes are slow enough so as to take a = 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 [START_REF] López | On the critical plunger speed and threedimensional effects in high-pressure die casting injection chambers[END_REF]. 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

dv dt v t (x) -v t-Δt (X) Δt , ( 33
)
where t and t -Δt represent the current and previous time steps, respectively. Similarly, X represents the position occupied by point x at time t -Δt. The weak form of the problem can thus be expressed as

Ω 2μd : d * dΩ - Ω pI : d * dΩ = - Ω ρbv * dΩ + Ω ρ dv dt v * dΩ (34) Ω ∇ • v p * dΩ = 0 (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

∇ • σ = 0, (36)
and

∇ • v = 0. ( 37
)
This mechanical problem is coupled with the following heat transfer equations:

∇(k∇T ) + ṙ -(ρc p Ṫ ) = 0 ( 38
)
where k denotes the thermal conductivity, ṙ the heat generation rate, ρ 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 ṙ = βσ : d [START_REF] Martínez | Fixed mesh and meshfree techniques in the numerical simulation of injection processes involving short fiber suspensions[END_REF] where β represents the fraction of mechanical energy transformed to heat that has been assumed to be 0.9 [START_REF] Zhou | 3D FEM simulation of the whole cycle of aluminium extrusion throughout the transient state and the steady state using the updated Lagrangian approach[END_REF].

If we write the incremental variational equation at time t + Δt we arrive to:

Ω(t+Δt) -(p t + Δp)I + 2μ(d t + Δd)(d t + Δd) : d * dΩ = 0. ( 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

Ω(t+Δt) -ΔΔpI + 2μ ∂μ(d t+Δt k ) ∂d : ΔΔd d t+Δt k + +2μ(d t+Δt k )ΔΔd : d * dΩ = = - Ω(t+Δt) (-p t+Δt k I + 2μ(d t+Δt k )d t+Δt k ) : d * dΩ, ( 41
)
where the subscript k indicates the iteration within a time increment. The incremental form of the incompressibility condition results

Ω(t+Δt) ∇ • (ΔΔv) p * dΩ = - Ω(t+Δt) ∇ • (v t+Δt k )p * dΩ. ( 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)

ΔΔv h (x) = nn I=1 φ I (x)ΔΔv I ; v * (x) = nn I=1 φ I (x)v * I (43) ΔΔp h (x) = nn I=1 ψ I (x)ΔΔp I ; p * (x) = nn I=1 ψ I (x)p * I . ( 44
)
where nn represents the number of nodes considered in the approximation. The functions ψ I (x) and φ I (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:

K t+Δt k G G T 0 ΔΔv ΔΔp = f t+Δt k 0 (45)
being the update given by

v t+Δt k+1 = v t+Δt k + ΔΔv ( 46
)
p t+Δt k+1 = p t+Δt k + ΔΔp (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 ψ and φ 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.

Mixed Natural Neighbor Approximations

In the problem defined by Eqs.( 36) and [START_REF] Martin | An experimental study of the collapse of liquid columns on a rigid horizontal plane[END_REF] there are two essential variables, namely velocity and pressure. Thus, we deal with a Hellinger-Reissner-like variational principle. It is well know 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 [START_REF] Ladyzhenskaya | The Mathematical Theory of Viscous Incompressible Flows[END_REF] [4] [START_REF] Brezzi | On the existence, uniqueness and approximation of saddle-point problems arising from lagrange multipliers[END_REF], together with the ellipticity condition.

The discrete form of the LBB condition may be written as:

inf p h ∈P h sup v h ∈U h Ω p h ∇ • v h dΩ ||p h || 0 ||v h || 1 = γ h ≥ γ > 0 ( 48
)
where γ is a positive constant independent of the mesh size, h. U h and P h represent, respectively, the displacement and pressure approximation spaces whose norms are defined as

|| • || 2 0 = Ω (•) 2 dΩ ( 49
)
|| • || 2 1 = Ω 3 i,j=1 ∂(•) i ∂x j 2 dΩ ( 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 [START_REF] Chapelle | The inf-sup test[END_REF] or [START_REF] Bathe | Finite Element procedures[END_REF], 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 [START_REF] Sukumar | The Natural Element Method in Solid Mechanics[END_REF], an equivalent form is developed:

inf w h ∈P h sup u h ∈U h W T h G h U h W T h G h W h • U T h S h U h = γ h ≥ γ > 0 (51)
where W h and U h are vectors of the nodal displacements corresponding to w h and u h . G h and S h are the matrices associated with the norms

||p h || 2 0 = W T h G h W h ( 52
)
||u h || 2 1 = U T h S h U h ( 53
)
G h and S h are positive semidefinite and positive definite [START_REF] Chapelle | The inf-sup test[END_REF], respectively. It is then demonstrated that the first non-zero eigenvalue, λ k , of the problem

G h φ h = λS h φ h ( 54
)
is related to the searched value γ h through the expression

γ h = λ k . ( 55
)
The obvious problem when developing mixed natural neighbor approximations is that, as mentioned before, natural neighbor interpolants are, at most, linearly complete. In [START_REF] González | Volumetric locking in Natural Neighbour Galerkin methods[END_REF] a solution is proposed to this problem based on the use of the partition of unity property (see [START_REF] Sukumar | The Natural Element Method in Solid Mechanics[END_REF]) 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 [START_REF] Babuška | The partition of unity finite element method: Basic theory and applications[END_REF].

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 [START_REF] Martínez | Fixed mesh and meshfree techniques in the numerical simulation of injection processes involving short fiber suspensions[END_REF] [START_REF] Martinez | Updated Lagrangian free surface flow simulations with Natural Neighbour Galerkin methods[END_REF] and rarely presents locking. Nevertheless, some enrichments have been studied in [START_REF] González | Volumetric locking in Natural Neighbour Galerkin methods[END_REF] 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, 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 [START_REF] Sukumar | The Natural Element Method in Solid Mechanics[END_REF]. 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 [START_REF] Bathe | Finite Element procedures[END_REF].

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 [START_REF] Martinez | Updated Lagrangian free surface flow simulations with Natural Neighbour Galerkin methods[END_REF] 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.

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 [START_REF] Edelsbrunner | On the shape of a set of points in the plane[END_REF] [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 ≤ α < ∞. 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:

(i) If σ T ∈ C then σ T ∈ C for every T ⊆ T .
(ii) 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 [START_REF] Edelsbrunner | Three dimensional alpha shapes[END_REF] 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:

S α = |C α | ∀0 ≤ α < ∞ ( 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 twodimensional 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 Equivalent plastic strain rate is depicted 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 [START_REF] Cueto | Imposing essential boundary conditions in the Natural Element Method by means of density-scaled α-shapes[END_REF] and [START_REF] Cueto | On the imposition of essential boundary conditions in Natural Neighbour Galerkin methods[END_REF] for an in-deep demonstration of this property). We have coined this method the α-shape based NEM, or simply α-NEM [START_REF] Cueto | Imposing essential boundary conditions in the Natural Element Method by means of density-scaled α-shapes[END_REF].

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.

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 [START_REF] Alfaro | Three-dimensional simulation of Aluminium extrusion by the α-shape based Natural Element Method[END_REF].

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: 2. Velocity updating. Compute v t+Δt 0 using Eq.( 41) and ( 42) assuming the temperature T t unchanged (the temperature field affects the pseudo-viscosity parameter μ (Eq.( 18)) through the expression of the yield stress (Eq.( 12)).

x t+Δt I = x t I + v t I Δt
3. Temperature updating. Update T t+Δt 0 solving Eq. [START_REF] Martinez | Updated Lagrangian free surface flow simulations with Natural Neighbour Galerkin methods[END_REF]. The material derivative of the temperature field is discretised in a Lagrangian form, i. e.,

dT dt I = T (x t+Δt I ) -T (x t I ) Δt
The resulting thermal diffusion problem is integrated by a semi-implicit Cranck-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 t+Δt Note that, by using a fully explicit time discretisation, the previous algorithm becomes uncoupled at each time step. Also if, after convergence in step [START_REF] Babuška | The partition of unity finite element method: Basic theory and applications[END_REF], 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 [START_REF] Alfaro | Three-dimensional simulation of Aluminium extrusion by the α-shape based Natural Element Method[END_REF].

CHECKING THE ACCURACY OF THE α-NEM APPROACH

The α-NEM accuracy and capabilities have been analyzed in some of our former works [START_REF] Cueto | Imposing essential boundary conditions in the Natural Element Method by means of density-scaled α-shapes[END_REF] [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.

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 [START_REF] Martin | An experimental study of the collapse of liquid columns on a rigid horizontal plane[END_REF]. 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, 81m/s 2 and the density of the water is 1000kg/m 3 . The fluid is assumed to be purely newtonian. Its viscosity was set to μ = 0.1P a • 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 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 and15 the absolute and relative errors in mass conservation are depicted. This error is below 3 % during the whole simulation.

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 [START_REF] Chinesta | α-NEM and model reduction: two new and powerful numerical techniques to describe flows involving short fibers suspensions[END_REF] 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 [START_REF] Chinesta | α-NEM and model reduction: two new and powerful numerical techniques to describe flows involving short fibers suspensions[END_REF] 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.

Flow kinematics

The governing equations for a short fiber suspension include the equilibrium equations and mass conservation, already described in Eqs.( 36) and [START_REF] Martin | An experimental study of the collapse of liquid columns on a rigid horizontal plane[END_REF]. In addition, the constitutive equation for such a suspension is described by:

σ = -pI + 2μ{D + N p (a 4 : D)} (57)
where N p is a parameter depending on the fiber concentration and the fiber aspect ratio. a 4 is the fourth-order fiber orientation tensor, defined by:

a 4 = ρ ⊗ ρ ⊗ ρ ⊗ ρψ(ρ)dρ ( 58
)
where ρ is the unit vector aligned with the fiber axis and ψ(ρ) is the orientation distribution function.

It is also possible to define a second-order orientation tensor, a, such that a = ρ ⊗ ρψ(ρ)dρ [START_REF] Zienkiewicz | Flow of solids during forming and extrusion: some aspects of numerical solutions[END_REF] Considering spheroidal fibres immersed in a dilute suspension, the orientation evolution is described by the so-called Jeffery equation:

dρ dt = Ωρ + k Dρ -(D : (ρ ⊗ ρ))ρ ( 60
)
where Ω 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

dψ dt + ∂ ∂ρ ψ dρ dt = 0 (61)
Taking into account Eqs.( 59), ( 60) and (61) we arrive to a general expression for the evolution of the second-order orientation tensor:

da dt = Ωa -aΩ + k Da + aD -2(a 4 : D) (62)

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:

θ i = i × 2π N , i∈ [0, . . . , N -1] (63)
In two dimensions, the orientation vector can be defined by

ρ = cos θ sin θ (64)
and, introducing it in the Jeffery's equation and after some simple algebraic manipulation we get: and thus the dependence of the rotation velocity on: (i) the fiber orientation θ ; (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 P a.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 [START_REF] Chinesta | α-NEM and model reduction: two new and powerful numerical techniques to describe flows involving short fibers suspensions[END_REF], 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 ability of the FEM to accurately simulate this class of problems can be argued. In threedimensional 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.

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.

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. [START_REF] Cueto | Modeling three-dimensional piece-wise homogeneous domains using the α-shape based Natural Element Method[END_REF]. We have considered μ 0 = 1.0 and n = 0.3. These values are obviously nonphysical, 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 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.

Two-Dimensional Simulation of Aluminium Extrusion

In this section we explore the ability of the α-NEM to simulate a simple, academic, twodimensional 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. 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 [START_REF] Alfaro | Meshless methods with application to metal forming[END_REF]. To obtain a profile velocity of about 0.5m/min a velocity of 1 mm/s is imposed. 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 α = 0.02 mm was taken. Different α-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.

in Figure 27 the equivalent strain rate is depicted. These values are in agreement with those of Lof et al [START_REF] Lof | Developments in finite element simulations of aluminium extrusion[END_REF] by employing ALE techniques.

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 [START_REF] Zhou | 3D FEM simulation of the whole cycle of aluminium extrusion throughout the transient state and the steady state using the updated Lagrangian approach[END_REF]. In Figures 28 and29 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 mm s , in order to obtain an exit velocity of 1 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 and33 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 and35. 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 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 Equivalent strain rate (s -1 ) for the extrusion simulation 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).

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. Meshless methods allow the use of an updated Lagrangian approach to this kind of nonstationary 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 α-shapes is used that allows for a very accurate extraction of the geometry of the domain, without the need of any geometrical check 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.

Figure 1 .

 1 Figure 1. Scheme of an extrusion process and associated evolution of the pressure on the ram

Figure 2 .

 2 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

Figure 3 .Figure 4 .

 34 Figure 3. Definition of the Natural Neighbour coordinates of a point x

Figure 5 .

 5 Figure 5. Definition of non-Sibsonian coordinates

Figure 6 .Figure 7 .

 67 Figure 6. Cloud of 9 nodes on a regular lattice

Figure 8 .Figure 9 .

 89 Figure 8. Geometry of the problem used to evaluate the numerical inf-sup test

Figure 10 .

 10 Figure 10. Two snapshots of a two-dimensional simulation of an extrusion process.Equivalent plastic strain rate is depicted

Figure 11 .

 11 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 ∞ (the convex hull of the set)

0 . 4 . 1 . 5 . 1 . 6 .

 041516 Velocity correction. We recompute step (2) assuming T t+Δt 0 in the yield stress expression, leading to the new velocity field v t+Δt Temperature correction. We recompute step (3) using v t+Δt 1 in the rate of generation expression, leading to the field T t+Δt Convergence check. then proceed to the next time step. Otherwise return to (4).

Figure 12 .

 12 Figure 12. Results for the two-dimensional simulation of a water collum collapse

Figure 13 .Figure 14 .

 1314 Figure 13. Comparison between experimental and numerical results for the water column collapse problem

Figure 15 .

 15 Figure 15. Relative error vs. time in the mass conservation for the water column collapse

Figure 16 .

 16 Figure 16. Cavity geometry

θ = -sin θ cos θ

Figure 17 .

 17 Figure 17. Evolution of the orientation field

Figure 18 .

 18 Figure 18. Simulation of an extrusion flow. Velocity contour (a) and arrow (b) map and α-shape at an intermediate time step

Figure 19 .Figure 20 .

 1920 Figure 19. The piece to be forged (a) and schematic representation of the forging process (b)

  d

  (d) 168th time increment.

Figure 21 .

 21 Figure 21. Equivalent strain rate distribution (s -1 ) Throughout the forging process

P

  (d) 168th time increment.

Figure 22 .

 22 Figure 22. Pressure ( N mm 2 ) distribution throughout the forging process

Figure 23 .Figure 24 .

 2324 Figure 23. Extrusion of a flat profile

Figure 25 .

 25 Figure 25. α-shapes at different time steps for the two-dimensional extrusion simulation

  (a) 101th time step. (b) 150th time step. (c) 200th time step. (d) 275th time step.

Figure 27 .

 27 Figure 27. Equivalent strain rate (s -1 ) at different time steps

Figure 28 .Figure 29 .Figure 30 .

 282930 Figure 28. Geometry of the container and the die for the cross-shaped profile problem

Figure 31 .

 31 Figure 31. Equivalent strain rate (s -1 ) for the extrusion simulation

(

Figure 33 .

 33 Figure 33. Presure ( N mm 2 ) at the 40th time step

  (d) 40th time step.

Figure 36 .

 36 Figure 36. Yield stress ( N mm 2 ) through the extrusion process

Figure

  Figure 37. Comparison of the obtained velocities with and without slipping boundary conditions ( mm s )

  (b) Sticking boundary conditions.

Figure 38 .

 38 Figure 38. Comparison of the temperature (K) profiles for the extrusion simulation

  3 two-and threedimensional examples of aluminium extrusion are studied.

		0	0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.018 0.02
	0	0.01	0.02	0.03	0.04	0.05	0.06	0.07	0.08	0.09	0.1

Table 1 .

 1 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 Material parameters for AA6063 Aluminium alloy (values taken[START_REF] Lof | FEM simulations of the extrusion of complex thin-walled aluminium sections[END_REF])

	Parameter Units	Value
	d 0	(s -1 )	0.005
	S m	(MPa)	25
	m		5.4
	A	(s -1 )	6× 10 9
	Q R	( J mol ) 1.4 × 10 5 ( J molK ) 8.314

  a)

		P
		82.9894
		77.0616
		71.1337
		65.2059
		59.2781
		53.3503
		47.4225
		41.4947
		35.5669
		29.6391
		23.7112
		17.7834
		11.8556
		5.92781
		0
		(b)
	Figure 32. Presure (	N mm 2 ) at the 20th time step
		P
		94.7663
		87.9972
		81.2282
		74.4592
		67.6902
		60.9212
		54.1521
		47.3831
		40.6141
		33.8451
		27.0761
		20.3071
		13.538
		6.76902
		0
	(a)	(b)

ACKNOWLEDGEMENT

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

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.