M A Martínez

E Cueto

I Alfaro

Manuel Doblaré
email: mdoblare@posta.unizar.es

F Chinesta

Updated Lagrangian free surface flow simulations with natural neighbour Galerkin methods

Keywords: meshless methods, natural element method, -shapes, Newtonian fluid, Norton-Hoff plasticity, short fibre-reinforced thermoplastics

In this paper, a new method to simulate free surface fluid flows within an updated Lagrangian framework is described. It is based on the use of a meshless technique coined as natural element method (NEM) or, more recently, as natural neighbour Galerkin method. The position of the flow front or the geometry of the fluid domain is handled by invoking the geometrical concept of -shape of the cloud of points, thus avoiding the explicit definition of the boundary of the domain as it evolves. This method also avoids the traditional need of remeshing typical in finite element simulations of this kind of processes. Three types of fluid behaviour have been considered, namely a purely Newtonian fluid, a non-Newtonian short fibrereinforced thermoplastic, and finally a Norton-Hoff viscoplastic behaviour. Benchmark examples showing the performance of the technique are included in the paper.

INTRODUCTION

In the last decade a great effort of research has been paid by researchers in the development of a wide family of methods coined as meshless methods (see, among a wide possible bibliography, [START_REF] Belytschko | Element-free Galerkin methods[END_REF][START_REF] Babuška | The partition of unity method[END_REF][START_REF] Liu | Reproducing kernel particle methods[END_REF][START_REF] Duarte | An H-p adaptive method using clouds[END_REF][START_REF] Atluri | New concepts in meshless methods[END_REF][START_REF] De | Towards an efficient meshless computational technique: the method of finite spheres[END_REF]). These methods are quite different in nature, although share some common characteristics. Their main advantage is the ease to simulate finite strains due to the fact that the 'element' connectivity is computed in such methods in a process transparent to the user. The traditional burden associated to the construction of a mesh in a finite element (FE) 1 model is thus avoided. This class of methods are also not sensitive to 'mesh' distortions, thus constituting an attractive alternative to finite element techniques.

Although not fully competitive with FE techniques in some aspects such as the imposition of essential boundary conditions, or some others like numerical integration, the use of meshless methods appears very appealing in simulating physical phenomena like those involving high mesh distortion. In particular, free surface flows of Newtonian or non-Newtonian fluids is one of these fields. Numerical simulation of free surface flows has an outstanding importance in many industrial problems, such as mould filling, metal casting, extrusion and many others.

As remarked in Reference [START_REF] Lewis | A mixed lagrangian-eulerian approach to modelling fluid flow during mould filling[END_REF], the natural approach to deal with moving boundary problems is a Lagrangian (or Lagrangian-Eulerian) framework, since these methods enable us to track explicitly the free surface position at each time step. In an updated Lagrangian approach the position of the mesh is updated at each time step, thus being able to determine the position of the flow front by taking into account the velocity of the nodes at the front itself. However, in finite element simulations, this method has the important drawback of mesh distortion. In two-dimensional problems this is somehow solvable, but in three-dimensional simulations high-quality meshers are difficult to obtain.

In an arbitrary Lagrangian-Eulerian approach [START_REF] Donea | Arbitrary Lagrangian-Eulerian finite element methods[END_REF], the mesh is moving with an arbitrary velocity. This velocity is incorporated to the convective term of the equations, somehow alleviating mesh distortions.

A third possibility in the simulation of free surface flows is the use of fixed mesh (fully Eulerian) approaches [START_REF] Montón | A 3d mold filling simulation for short fibres reinforced thermoplastics[END_REF]. In this last case, it is necessary to introduce a new variable in the problem (often referred to as presence of fluid) that takes into account if the element is full, empty or in the filling process. This variable-that ranges from 0 to 1-is advected with the flow in order to determine the front position. This approach has the drawback of introducing some numerical diffusion in the position of the front, that becomes unclear, since it is located somewhere along a strip of elements partially full. Finally, the empty region of the mould is assigned fictitious material properties. Other techniques substitute the presence of fluid function by the employ of a level set approach to determine the front position within the fixed mesh [START_REF] Sussman | A level set approach for computing solutions to incompressible two-phase flows[END_REF].

Meshless methods represent an alternative to the techniques explained before. They allow for an updated Lagrangian approach of the problem without care for mesh distortions-in fact this is one of the main properties of meshless methods, the accuracy does not depend significantly on the 'quality' of the mesh in an FE sense-. In this paper, a natural neighbour Galerkin approach [START_REF] Sukumar | The natural element method in solid mechanics[END_REF][START_REF] Cueto | Imposing essential boundary conditions in the natural element method by means of density-scaled -shapes[END_REF][START_REF] Cueto | Modeling three-dimensional piece-wise homogeneous domains using the -shape based natural element method[END_REF] to free surface flows is presented and analysed.

In this case, a formulation based on the use of -shapes [START_REF] Edelsbrunner | On the shape of a set of points in the plane[END_REF][START_REF] Edelsbrunner | Three dimensional alpha shapes[END_REF] is introduced. This technique, firstly employed in a natural neighbour (often referred to as natural element) approach in Reference [START_REF] Cueto | Imposing essential boundary conditions in the natural element method by means of density-scaled -shapes[END_REF], allows to extract the shape of the domain being simulated by employing a cloud of nodes only. This alleviates the burden associated to geometrical checks in mesh inter-penetrations, etc., during the flow. The use of natural neighbour Galerkin methods is also motivated by the fact that, in contrast to most meshless methods, the interpolation scheme used in this technique is interpolant-instead of only approximant-and thus essential boundary conditions are imposed exactly in FE-like way.

This paper is organized as follows: in Section 2, the basics about natural neighbour interpolation are addressed. In Section 3, the formulations of the problems under consideration-namely, purely Newtonian fluids, short fibre-reinforced thermoplastics and Norton-Hoff plasticity-are described. In Section 4, some aspects related to mixed natural neighbour interpolation are addressed. Issues concerning volume tracking in an -shape based natural neighbour Galerkin approach are investigated in Section 5. In Section 6, some benchmark examples are analysed and, finally, conclusions are enumerated in Section 7.

NATURAL NEIGHBOUR INTERPOLATION

Most meshless methods employ interpolation techniques that do not rely in a background mesh to construct the approximation. Thus, element free Galerkin methods [START_REF] Belytschko | Element-free Galerkin methods[END_REF], hp-clouds [START_REF] Duarte | An H-p adaptive method using clouds[END_REF] or other methods [START_REF] Atluri | New concepts in meshless methods[END_REF][START_REF] De | Towards an efficient meshless computational technique: the method of finite spheres[END_REF] use moving least-squares approximation to build the discrete system of equations of the Galerkin method. By the contrary, natural neighbour Galerkin methods use different natural neighbour-based (and hence the name) interpolations. This kind of interpolation was firstly introduced by Sibson [START_REF] Sibson | A brief description of natural neighbour interpolation[END_REF] and is usually referred to as Sibsonian interpolation, in contrast to the non-Sibsonian or Laplace interpolation [START_REF] Belikov | Non-Sibsonian interpolation on arbitrary system of points in Euclidean space and adaptive isolines generation[END_REF][START_REF] Hiyoshi | Two generalizations of an interpolant based on Voronoi diagrams[END_REF]. In this section, both interpolation schemes are reviewed, although only Sibson interpolation will be used in the examples that are included in Section 6. For an in-deep review of the method, the interested reader can consult [START_REF] Cueto | Overview and recent advances in natural neighbour Galerkin methods[END_REF].

Voronoi diagrams of a cloud of points and Delaunay triangulations

Given a cloud of points N = {n 1 , n 2 , . . . , n M } ⊂ R d there is a unique decomposition of the space into regions such that each point within these regions is closer to the node to which the region is associated than to any other of the cloud. This kind of space decomposition is called Voronoi diagram of the cloud of points and each Voronoi cell is formally defined as (see Figure 1):

T I = {x ∈ R d : d(x, x I ) < d(x, x J ) ∀J = I } (1)
being d(• , •) the euclidean distance function.

The dual structure of the Voronoi diagram is the Delaunay triangulation, obtained by connecting nodes that share a common (d -1)-dimensional facet. While the Voronoi structure is unique, the Delaunay triangulation is not, there being the some so-called degenerate cases in which there are two or more possible Delaunay triangulations (consider, for example, the case of triangulating a square in 2D). Other way to define the Delaunay triangulation of a set of nodes is by invoking the empty circumcircle property, which means that no node of the cloud lies within the circle covering a Delaunay triangle. Two nodes sharing a facet of their Voronoi cell are called natural neighbours. The Delaunay triangulation is very useful in 2D mesh generation, since provides the optimal triangulation of set of points from the FE viewpoint. This is because of its property of maximizing the minimum interior angle of the triangulation. This property does not hold in three dimensions.

Sibson interpolation

Interpolation by natural neighbours was introduced by Sibson [START_REF] Sibson | A vector identity for the Dirichlet tesselation[END_REF] by defining the second-order Voronoi cell 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 I J = {x ∈ R 3 : d(x, x I ) < d(x, x J ) < d(x, x K ) ∀J = I = K} (2)
Thus, if a new node is added to a given cloud of points the Voronoi cells will be altered.

Sibson [START_REF] Sibson | A vector identity for the Dirichlet tesselation[END_REF] defined the natural neighbour co-ordinates 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 when adding x to the initial cloud of points to the total volume of T x . In other words, being (x) and I (x) the Lebesgue measures of T x and T xI , respectively, the natural neighbour co-ordinates of x with respect to the node I is defined as

I (x) = I (x) (x) (3)
In Figure 2, the shape function associated to node 1 may be expressed as being straightforward to prove that NE shape functions form a partition of unity [START_REF] Babuška | The partition of unity method[END_REF], as well as some other properties like positivity (i.e. I ∈ [0, 1]) and interpolation:

1 (x) = A abf e A abcd (4)
I (x J ) = I J ( 5
)
Sibson interpolation also satisfies the so-called co-ordinate property (see Reference [START_REF] Sukumar | The natural element method in solid mechanics[END_REF] and references therein) or, equivalently, linear completeness:

x = n I =1 I (x)x I (6)
being x I the nodal co-ordinates and n the number of natural neighbours of the point. This allows the use of Sibson co-ordinates to solve second-order partial differential equations in a Galerkin framework (Figure 3).

Laplace interpolation

The non-Sibsonian interpolant was introduced independently by Belikov and Semenov [START_REF] Belikov | Non-Sibsonian interpolation on arbitrary system of points in Euclidean space and adaptive isolines generation[END_REF] and by Hiyoshi and Sugihara [START_REF] Hiyoshi | Two generalizations of an interpolant based on Voronoi diagrams[END_REF] and is defined by using geometrical entities of one dimension less than the original space under consideration. If we define the cell intersection t I J = {x ∈ T I ∩ T J , J = I } (note that t I J may be an empty set) we can define the value Thus, the value of the shape function associated to node 4 at the point x (see Figure 4) is defined as

J (x) = |t I J | d(x, x J ) (7)
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)] (8)
This kind of interpolation shares many of its properties with Sibson's co-ordinates (like positivity, interpolation and linear completeness), being less computationally expensive, as noted in Reference [START_REF] Sukumar | Natural neighbor Galerkin methods[END_REF]. However, it is less smooth than the previous one since the derivatives of the non-Sibsonian-shape function are not defined along the edges of the Delaunay triangles that lie within its support. For the computations performed in this work, Sibson's interpolation has been preferred by the authors.

Properties of Sibson's interpolation

It is easy to prove (see Reference [START_REF] Sukumar | The natural element method in solid mechanics[END_REF]) that in two dimensions, Sibson's interpolation leads to constant strain finite element shape functions if the considered point has three neighbours and to bilinear finite element shape functions if it has four neighbours placed on a regular grid. In other case, Sibson's functions are rational quartic functions.

The interpolation thus constructed is smooth everywhere except from the nodes, where it is C 0 . A proof of this assertion can be found in Reference [START_REF] Sukumar | The natural element method in solid mechanics[END_REF]. It is also possible to construct a C 1 interpolation based on natural neighbours that can be applied to discretize the variational formulation of fourth-order partial differential equations, but these will not be discussed here. For more details, see Reference [START_REF] Sukumar | C 1 natural neighbour interpolant for partial differential equations[END_REF] and references therein.

One of the main properties of Sibson's co-ordinates, as mentioned before, lies in its ability to strictly interpolate nodal values [START_REF] Sukumar | The natural element method in solid mechanics[END_REF]. This allows to enforce essential boundary condition in a straightforward manner if we deal with convex domains. If not, in Reference [START_REF] Cueto | Imposing essential boundary conditions in the natural element method by means of density-scaled -shapes[END_REF] and also in Reference [START_REF] Cueto | On the imposition of essential boundary conditions in natural neighbour Galerkin methods[END_REF] two distinct proofs of the ability of the method to interpolate linear fields along the boundary of any kind (convex or not) of domains are developed. The first of these proofs is based on the use of the concept of -shapes of the cloud of points [START_REF] Edelsbrunner | On the shape of a set of points in the plane[END_REF], that will be revisited in Section 5.

This makes natural neighbour Galerkin methods an attractive choice among the wide field of meshless methods because of their simplicity in handling essential boundary conditions.

Computational costs

The key difference between the formulation here proposed and one based upon finite element techniques (like in Reference [START_REF] Lewis | A mixed lagrangian-eulerian approach to modelling fluid flow during mould filling[END_REF]) is the computation of shape functions. Prior to this, it is necessary to compute the Delaunay triangulation of the cloud of nodes. However, this is a well-known procedure for computational geometers and can be done in a very fast way (see Reference [START_REF] Sukumar | The natural element method in solid mechanics[END_REF] and references therein). The search of natural neighbours for a given point, usually an integration point, is by far the most time-consuming task. An obvious choice is to perform a search through all the triangles for each integration point (O(n 2 ) algorithm), but for a big number of nodes other more sophisticated algorithms (like the walking triangle [START_REF] Lawson | Software for C 1 surface interpolation[END_REF]) are preferable.

Globally, the computational cost in NEM has been estimated in about two times more than the traditional finite element matrix computation [START_REF] Sukumar | The natural element method in solid mechanics[END_REF]. Note, however, that this does not include the time to perform geometrical checks-interpenetration of boundaries, waves or holes formation, for instance-, as discussed in Section 5, necessary in a mesh-based approach.

GOVERNING EQUATIONS

In this paper, we consider three types of material behaviour in order to show the advantages of using natural neighbour Galerkin methods in the simulation of free surface flows. They correspond to the benchmark problems presented in Section 6. Firstly, a purely Newtonian fluid behaviour is considered in Section 3.1. In Section 3.2, a model that takes into account the influence of short fibre reinforcement in Newtonian fluids is presented. The equations that govern the behaviour of such a fluid have a hyperbolic character. However, its characteristic lines coincide with the nodal pathlines, thus making straightforward the use of the method of characteristics [START_REF] Chinesta | Numerical simulation of the steady recirculating flows of fiber suspensions[END_REF] to integrate the constitutive equation. Finally, in Section 3.3 a Norton-Hoff plasticity model is presented. In such a model, the viscosity of the fluid varies with the strain rate tensor, thus leading to a fully non-linear behaviour that has been linearized through the Newton-Raphson method.

In essence, the three models can be viewed within a unified framework, being their governing equations the following:

1. Equilibrium equations (balance of linear momentum in the absence of inertial and body forces):

∇ • = 0 (9)
2. Incompressibility of the fluid:

∇ • v = 0 ( 10
)
The difference between the three models lies in the constitutive equation, which is described in detail in the following sections.

Newtonian fluid behaviour

The simplest behaviour for the fluid is described by the purely Newtonian constitutive equations, namely

= -pI + 2 D (11)
being , the viscosity, assumed constant, I the second-order unit tensor and D the strain rate tensor.

The variational (weak) formulation associated to Equations ( 9)-( 11) results

(t)

:

D * d = t t • v * d ∀v * ∈ V ( 12
) (t) (-∇v + p)p * d = 0 ∀p * ∈ L 2 ( (t)) ( 13
)
where

U = {v|v ∈ (H 1 ( (t))) 2 , v| v = v}, V = {v * |v * ∈ (H 1 ( (t))) 2 , v| v = 0},

Short fibre-reinforced thermoplastics

In this section, we consider the model for a short fibre-reinforced thermoplastic. This model has been previously presented and tested in a number of previous works (see Reference [START_REF] Azaiez | State-of-the-art on numerical simulation of fiber-reinforced thermoplastic forming processes[END_REF] and references therein) from a Eulerian point of view, but its main features make it appealing for its use within an updated Lagrangian framework. In essence, we consider the previous constitutive equations (Equation ( 11)) with an additional term which describes the contribution of the fibres. This leads to = -pI + 2 {D + N p tr(aD)a} [START_REF] Edelsbrunner | On the shape of a set of points in the plane[END_REF] where p denotes the pressure, I the unit tensor, the equivalent suspension viscosity, D the strain rate tensor, N p a scalar parameter depending on both the fibre concentration and the fibre aspect ratio and a the second-order orientation tensor defined by

a = ⊗ ( ) d ( 15
)
where is the unit vector defining the fibre axis direction, and ( ) is the orientation distribution function, verifying the normality condition

( ) d = 1 (16)
If ( ) = ( -ˆ ), with () the Dirac's distribution, all the orientation probability is concentrated in the direction defined by ˆ , and the corresponding orientation tensor is â = ˆ ⊗ ˆ , where ⊗ denotes the dyadic product. In a planar case, the isotropic orientation is defined by the uniform distribution function

( ) = 1 2 (17)
The orientation tensor related to a planar isotropic orientation distribution is then

a = I 2 (18)
From a physical point of view, we can consider that the eigenvalues of the second-order orientation tensor (a) represent the probability of finding the fibres in the direction of the corresponding eigenvectors.

With a quadratic closure relation [START_REF] Azaiez | State-of-the-art on numerical simulation of fiber-reinforced thermoplastic forming processes[END_REF] the orientation equation is expressed as

da dt = a -a + k(Da + aD -2 tr(aD)a) + D r a - I 2 ( 19
)
where a satisfies

a = a t tr(a) = 1 (20)
is the vorticity tensor (skew-symmetric component of ∇v), k is a constant that depends on the fibre aspect ratio r (fibre length to fibre diameter ratio): k = (r 2 -1)/(r 2 + 1), and D r a diffusion coefficient to take into account the fibre-fibre interactions.

The quadratic closure relation is exact as long as the fibres are locally perfectly aligned. If we consider other closure relations, the orientation equation is changed, but its hyperbolic character remains unchanged. Thus, all the numerical procedures proposed in this paper can be also applied when other closure relations are considered.

The flow model is defined in the volume occupied by the fluid at time t, (t). On its boundary, (t) ≡ * (t) either the velocity or the traction are imposed:

v(x ∈ v ) = v (21) n(x ∈ t ) = t ( 22
)
where v is the part of where non-null velocity is prescribed (inflow boundary). On the other hand, the orientation equation ( 19) is hyperbolic, thus needing an initial value for its integration a(x ∈ -) = a 0 [START_REF] Cueto | On the imposition of essential boundary conditions in natural neighbour Galerkin methods[END_REF] With the flow kinematics known at time t, v t (x), the position of nodes can be updated at the same time that the fibres orientation evolution, given by Equation [START_REF] Cueto | Overview and recent advances in natural neighbour Galerkin methods[END_REF]. This is computed by using the method of characteristics. The simplest explicit updating consists in writing

x t+ t I = x t I + v t I t, ∀I (24)
and

a
where D t I and t I are the symmetric and skew-symmetric components of the gradient of velocity tensor, at time t in the node x I , respectively. It is seen how, by virtue of the meshless character of the natural neighbour Galerkin methods, the employ of the method of characteristics can be afforded in a straightforward manner. Other higher-order explicit strategies can also be used.

The weak form of Equation ( 9) with the substitution of material behaviour (Equation ( 14)) can be formulated as

(t+ t) {-pI + 2 (D t+ t + N p tr(a t+ t D t+ t )a t+ t )} : D * d = 0 ∀v * ∈ V ( 26
)
where a null traction is assumed on the fluid flow front t , and D t+ t is the strain rate tensor related to the velocity field v (t+ t) (x). p also refers to the pressure at time t + t. Here, the dependence of the weak form of the problem on the fibre orientation tensor a is clear. At the same time, this orientation depends on the velocity field, as expressed by Equation (25).

Norton-Hoff viscoplasticity

The third type of constitutive equation considered here is the so-called generalized Newtonian behaviour, whose description can be made as follows:

= -pI + 2 (D)D ( 27
)
If we assume a Norton-Hoff plasticity model (viscoplastic flow with null yield stress), viscosity is a function of the second invariant of the strain rate tensor, namely

(D) = 0 ( √ 2D : D) n-1 (28)
being 0 the so-called consistency coefficient and n the pseudo-plasticity coefficient.

If we write the incremental variational equations about time t we arrive to

(t+ t) (-(p t + p)I + 2 (D t + D)) : D * d = 0 ( 29
)
Owing to the non-linear character of the constitutive equations, an iterative approach has been accomplished, through 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 ( 30
)
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 (31)

MIXED NATURAL NEIGHBOUR APPROXIMATION

We review here how velocities and pressures are approximated in the formulations before presented. The velocity vectors v h (x) : ⊂ R 2 → R 2 are interpolated using the Sibson shape functions, and can be written in the form

v h (x) = n I =1 I (x)v I ( 32
)
where v I (I = 1, 2, . . . , n) are the vectors of nodal velocities at the n natural neighbours of point x, and I (x) are the shape functions associated with each node. The pressure field is interpolated by an interpolation scheme similar to that in Equation ( 32):

p h (x) = n I =1 I (x) pI (33)
where I (x) is the shape function associated with node I and pI are the nodal pressures. In Reference [START_REF] Sukumar | The natural element method in solid mechanics[END_REF] I (x) is proposed to be of the general form

I (x) = ( I (x)) k n J =1 ( J (x)) k (34)
where k is a non-negative integer and n is the number of natural neighbours of the point x. The case k = 1 reduces to natural neighbour interpolation, and all other k render shape functions that form a partition of unity, with the resulting interpolant satisfying only constant completeness. The case k = 1 showed spurious pressure modes in the tests performed during the development of this work. However, the case k = 0 has been employed in a number of previous works [START_REF] Martínez | Natural element meshless simulation of injection processes involving short fiber suspensions[END_REF], showing no spurious pressure modes:

p h (x) = n I =1 1 n pI (35)
Since the shape functions I (x) have compact support [START_REF] Farin | Surfaces over Dirichlet tesselations[END_REF] and have C 0 continuity, the shape functions I (x) = 1/n have C -1 continuity. Other C -1 pressure interpolants can be constructed using nearest neighbour interpolation [START_REF] Preparata | Computational Geometry: An Introduction[END_REF]. In this approach, if x ∈ T I , where T I is the Voronoi polygon related to node I , then we assign the nodal pressure pI to p(x). In essence, the pressure field is assumed to be a constant over the Voronoi polygon T I . Thus, it results

p(x) = pI if x ∈ T I

(36)

TRACKING THE BOUNDARY OF THE DOMAIN

The issue of tracking the boundary of the domain as it evolves is crucial in free surface simulations. In Lagrangian FE simulations the nodes lying at the boundary are usually detected by marking coincident edges of elements (see Reference [START_REF] Lewis | A mixed lagrangian-eulerian approach to modelling fluid flow during mould filling[END_REF], for instance). Thus, non-coincident edges must lie at the boundary of the domain. In most meshless methods this approach is not useful, since the concept of element edge disappears. The formation of bubbles as the fluid evolves is difficult to handle (especially in three-dimensional cases), since all boundary edges must be checked not to intersect with other boundary edges. A similar approach must be accomplished with internal obstacles that may appear in the domain. These obstacles must be discretized and the intersection of the boundary with them must be checked in each time step.

The approach here presented is essentially different. It relies in the well-known geometric concept of -shapes [START_REF] Edelsbrunner | On the shape of a set of points in the plane[END_REF][START_REF] Edelsbrunner | Three dimensional alpha shapes[END_REF]. The problem can be formulated briefly as follows: can the cloud of points itself contain information enough about the geometry of the domain over which it is defined? If the density of the cloud is enough, the answer is affirmative. Thus, there is no need to take into account boundary segments and all the geometrical checks they involve.shapes have also been employed in References [START_REF] Cueto | Imposing essential boundary conditions in the natural element method by means of density-scaled -shapes[END_REF][START_REF] Cueto | Modeling three-dimensional piece-wise homogeneous domains using the -shape based natural element method[END_REF] to enforce essential boundary conditions, leading to what we have called the -shape based natural element method ( -NEM in what follows).

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

Following this, we can define the family of sets F k, as the sets of -exposed k-simplexes 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, . As remarked before, an -shape is a polytope that can be triangulated by a subset of the Delaunay triangulation or tetrahedrization, that is, by an -complex.

In the case of non-uniform nodal distributions, the possibility of using density-scaled -shapes [START_REF] Teichmann | Surface reconstruction with anisotropic density-scaled alpha shapes[END_REF] has also been studied in Reference [START_REF] Cueto | Imposing essential boundary conditions in the natural element method by means of density-scaled -shapes[END_REF]. An example of three-dimensional family of -shapes is shown in Figure 5. As shown, the resultant geometry of the domain ranges from the cloud of points itself for = 0 to the convex hull of the cloud for = ∞. If there exists a sufficiently dense nodal sampling it would be easy to find an value that gives an accurate geometry definition. See References [START_REF] Cueto | Imposing essential boundary conditions in the natural element method by means of density-scaled -shapes[END_REF][START_REF] Martínez | Natural element meshless simulation of injection processes involving short fiber suspensions[END_REF] for a deeper discussion on how the value of can affect the results.

In the simulations presented in Section 6 one has to deal also with the conservation of mass in the sense that the method used to extract the geometry of the domain must prevent from loss in volume of fluid. This is straightforward if the geometry of the domain is convex (any value strictly greater than the nodal spacing h would suffice) but somewhat complex if the geometry becomes non-convex. In the case of two or more fluid fronts that get in contact, care must be taken in avoiding fluid inter-penetration (and thus loss of volume). In this case, the -NEM, assume that the flow front interaction occurs when they approach to a distance lower than . Thus, the error in the mass conservation is of order , which at its turn is of the same order than the mesh density. In this form, we can conclude that the -NEM approximation converges towards a null error in the mass conservation as the mesh density is reduced.

On the other hand, in order to avoid material interpenetration, the time step must verify the condition

t = min x (x, t) v(x, t) (37)
in order to ensure that the path covered by a node would not be greater than the nodal distance h. These ideas will be further investigated in Section 6.1 below.

NUMERICAL RESULTS

Newtonian fluid passing over a cylinder

Consider the flow of a Newtonian fluid in a two parallel plates channel where a cylindrical obstacle is placed, as depicted in Figure 6. The initial cloud of nodes is composed of 384 nodes. Additional nodes are introduced at each time step through the inflow boundary placed at x = 0 with prescribed velocity v x (x = 0, y) = (4V /h 2 )(hy)y (V = 0.1 m s -1 ). Viscosity was set to = 1.0 Pa s and the time step t to 0.05 s. No-slipping boundary conditions were assumed between the fluid and the walls.

It is clear that, as explained in Section 5, for the actual volume of the domain to be represented, the parameter must be chosen such that the minimum radius of curvature of the boundary is properly reproduced. In this case, is bounded at least by the radius of the cylindrical obstacle, R = 0.0015. In the case in which > R the domain is not well reproduced, which is illustrated in Figure 7(a). In Figure 8, the horizontal velocity plot at the tenth time increment is depicted.

From Section 5, one can readily infer that a bad -value leads to the lack of fulfilment of the conservation of mass (or volume) throughout the simulation. In Figure 9, a comparison is shown between the analytical volume and the computed one, for three different values. It is seen that, if is too big a spurious increase in the volume of the domain is obtained (see Figure 7(a)). In Figure 7(b) the resultant -shape is depicted for an appropriate value of . On the other hand, if is too small a loss of volume is produced as the triangles become distorted by the fluid motion, since the distortion leads to an increase of the radius of the circumcircle. This can be circumvented, of course, by using refined clouds of points. Thus, is bounded between the nodal spacing h and the radius of the cylinder R to obtain a suitable geometric description. Here,

h < < R (38)
As commented before another source of error appears if the time increment t used is too large. In this case, if |v| t > the particle motion in a single time step is larger than the nodal spacing h, which would lead to an undesired mixture in the fluid downstream the cylinder. Thus, if = 0.001 and the velocity is in the order 0.1-0.01 m s -1 , then, as explained in the previous section, the appropriate time step results in the interval 0.01-0.1 s, which justify the value t = 0.05 s considered in these simulations. The resulting velocity field once the cylinder is overcome is depicted in Figure 10. It can be seen how a welding line is developed between the two colliding flow fronts, as predicted, for instance, in Reference [START_REF] Pichelin | Finite element solution of the 3D mold filling problem for viscous incompressible fluid[END_REF]

Mould filling with short fibre-reinforced thermoplastic

In this section, we analyse the filling of a benchmark mould with a short fibre-reinforced thermoplastic. Fluid behaviour is governed by the equations described in Section 3.2. The geometry of the mould is shown in Figure 11.

The suspension viscosity is fixed to 10 3 Pa s, the fibres are assumed to have a quasi-infinite aspect ratio k = 1, with a particle number of one or ten, N p = 1 or 10, and all diffusion effects neglected (i.e. D r = 0). At time t = 0 the suspension is assumed filling the inlet cavity (as shown in Figure 11) with an initial isotropic fibre orientation distribution (a(t = 0) = I/2). No-slipping boundary conditions were assumed between the fluid and the mould walls. In this simulation (depicted in Figure 12) we can notice the perfect flow symmetry and the good description of the flow front location. Figure 13 depicts the flow kinematics obtained using a fibre suspension characterized by the particle number N p = 10. As expected, for the first iterations of the filling process, when the flow is actually an extrusion flow, the die swelling increases with the particle number due to the normal stress differences increasing. In order to check the ability of the -shape technique to accurately extract the suspension volume at each time, we compare in Figure 14 the numerical volume of the fluid domain automatically extracted from the cloud of nodes with the volume of fluid really injected. We can notice the excellent accuracy in spite of the reduced number of nodes used in the simulation. Only some discrepancies appear between the theoretical volume and the predicted one in the final time steps, since the flow must fill the sharp corners. In this case, these errors are of the order of 5% at the final stages of the process, being under 1% at the most part of the simulation. Errors of about 3% are cited in Reference [START_REF] Lewis | A mixed lagrangian-eulerian approach to modelling fluid flow during mould filling[END_REF] by using arbitrary Lagrangian-Eulerian techniques and similar number of nodes, but applied to much simpler geometries.

Note the bad quality of the triangulation depicted in Figure 15 from an FE viewpoint. However, as emphasized before, NEM results are not significantly affected by the shape of the triangulation.

Finally, the fibre orientation field at an intermediate time step is depicted in Figure 16. In the graphical representation, the largest ellipse axis represent the most probable fibre orientation direction, and its length is proportional to that probability.

Simulation of extrusion processes

In this section we deal with the simulation of extrusion processes. We consider a benchmark example of an extrusion die whose geometry is shown in Figure 17. For the simulation, only the region of the metal nearest to the die was considered by enforcing appropriate displacements. The cloud of points is composed by 2989 nodes and remained unchanged throughout the simulation. The extruded metal was simulated as viscoplastic, with constitutive equations as described in Section 3.3, with consistency coefficient 0 = 1.0 MPa s and pseudo-plasticity n = 0.3. Note the highly non-linear character of the resultant behaviour. Similar material behaviour was employed in Reference [START_REF] Bonet | Correction and stabilization of smooth particle hydrodynamics methods with applications in metal forming simulations[END_REF] in the simulation of metal forging processes with corrected smooth particle hydrodynamics methods (CSPH). The mould is considered as perfectly rigid and slipping contact between metal and die was assumed. Of course, this model is a first attempt to validate the ability of natural neighbour Galerkin methods to handle such type of processes. The springback of the extruded metal can not be predicted, for instance, by assuming a Norton-Hoff plasticity model like the one exposed in Section 3.3.

Equivalent plastic strain (second invariant of the strain rate tensor, √ 2D : D) for time steps 1, 100 and 200 is depicted in Figure 18. It can be noticed the accuracy obtained in the volume prediction (see Figure 19). In this case, total volume error is under 0.25%, much lower than those predicted in references like [START_REF] Lewis | A mixed lagrangian-eulerian approach to modelling fluid flow during mould filling[END_REF] by using ALE techniques and a similar number of nodes. This is in spite of the high distortion of the Delaunay triangles achieved throughout the In this example, it can be seen how natural neighbour Galerkin methods constitute an appealing choice among numerical methods to simulate forming processes in general, and among meshless methods in particular.

CONCLUDING REMARKS

In this paper, some key aspects of the simulation of free surface flows from an updated Lagrangian point of view by using natural neighbour Galerkin methods have been addressed. Meshless methods have became an interesting alternative to more traditional Eulerian or arbitrary Lagrangian-Eulerian (ALE) approaches in the simulation of such processes. In particular, natural neighbour Galerkin methods offer a sound geometrical basis for the treatment of free boundaries through the use of -shapes. It has been pointed out how if the domain is modelled by an appropriated -shape-which in turn requires a sufficiently dense nodal sampling of the domain-not only very accurate volume predictions can be obtained, but also a direct imposition of essential boundary conditions with linear interpolation on the boundary. See References [START_REF] Cueto | Imposing essential boundary conditions in the natural element method by means of density-scaled -shapes[END_REF][START_REF] Cueto | Modeling three-dimensional piece-wise homogeneous domains using the -shape based natural element method[END_REF] for a deeper explain on this last topic. The treatment of evolving boundaries such as those forming bubbles is straightforward by employing the proposed technique and without recursive checks on the geometry of the free surface. Three different material behaviours have been considered, namely a purely Newtonian fluid, a short fibre-reinforced thermoplastic model that takes into account the influence of fibre orientation in the resultant velocity field and, finally, Norton-Hoff plasticity. These behaviours have been studied with the aid of three respective examples that showed the capabilities of the method. In particular, it has been shown how the treatment of complex materials, whose behaviour is governed by hyperbolic equations, can be greatly simplified by noting that the nodal pathlines in a Lagrangian description coincide with the characteristic lines of the equation. Thus, it is possible to directly integrate these equations explicitly (forward or Euler schemes) with great accuracy.

Figure 1 .

 1 Figure 1. Delaunay triangulation and Voronoi diagram of a cloud of points.

Figure 2 .

 2 Figure 2. Definition of the natural neighbour co-ordinates of a point x.

Figure 3 .

 3 Figure 3. Function 1 (x) for the cloud of points in Figure 2.

Figure 4 .

 4 Figure 4. Definition of non-Sibsonian co-ordinates.

Figure 5 .

 5 Figure 5. Evolution of the family of -shapes of a cloud of points representing a mandible. Shapes: (a) S 0 ; (b) S 1.0 ; (c) S 1.5 ; (d) S 3.5 ; and (e) S ∞ are depicted.

3 Figure 6 .

 36 Figure 6. Initial geometry and cloud of points for the flow between two plates problem.

Figure 7 .

 7 Figure 7. Comparison of the effect of a badly chosen -value: (a) too high value ( = 0.01); and (b) appropriate ( = 0.001).

Figure 8 .

 8 Figure 8. Velocity map at the tenth time step.

Figure 9 . 3 Figure 10 .

 9310 Figure 9. Comparison between the theoretically computed volume and the actually computed for three different values.

Figure 11 .

 11 Figure 11. Geometry of the benchmark mould.

Figure 12 .

 12 Figure 12. Filling process simulation: flow kinematics (N p = 1).

Figure 13 .Figure 14

 1314 Figure 13. Comparison of the flow kinematics with N p = 10 and 1: (a) twenty-fifth time step; (b) forty-fifth time step; and (c) sixty-fifth time step.

Figure 15 .

 15 Figure 15. Horizontal velocity contour plot over the -shape of the fluid at the 60th time step.

Figure 16 .

 16 Figure 16. Fibre orientation field at the 60th time step (N p = 1).

Figure 17 .

 17 Figure 17. Geometry of the extrusion die and detail of the simulated region.

Figure 18 .Figure 19 .Figure 20 .

 181920 Figure 18. Equivalent plastic strain for time steps 1, 100 and 200.

Figure 21 .

 21 Figure 21. Velocity field at the 200th time step for the extrusion problem.

Contract/grant sponsor: Spanish Ministry of Science of Technology; contract/grant number: CICYT-DPI2002-