A Fabri

M Teillaud
email: monique.teillaud@inria.fr

Presqu'île de Giens �Var) CGAL -The Computational Geometry Algorithms Library

Keywords: geometric data structures and algorithms, open source, C++, generic programming, exact geometric computing

We present fundamental geometric data structures and algorithms offered by CGAL, the Computational Geometry Algorithms Library. As geometry is ubiquitous this library is used by application developers in a variety fileds such as medical imaging, VLSI, CAD/CAM, geophysics, computer graphics, GIS. . . In the presention we will give a broad overview, and go in more detail on surface and volume mesh generation, as well as on geometry for periodic domains.

Introduction

CGAL, the Computational Geometry Algorithms Library, is a C++ library of geometric algorithms and data structures which is developed by the CGAL Open Source project [1]. This project started in 1996 as a joint effort of several research groups working in computational geometry, and it has been partially funded by the European and national science foundations. CGAL consists of 600K lines of code, organized in 80 packages, documented in a 3500 pages user and reference manual, with 1000 subscribers to the user forum, and about 100 commercial users, about 20 active developers, and two types of licenses : open source and commercial. The project is steered by an Editorial Board, which is also responsible for reviewing submissions of new packages.

Algorithms

As CGAL is not an end user application, but a collection of geometric software components, the demos given during the presentation serve as illustrations of the algorithms and they are not specific for a particular application domain. The goal of the demos is to present an overview of the functionality, to show that CGAL is robust and fast, and can deal with complex data sets.

2D Triangulations and Mesh Generation : CGAL offers triangulations of point sets [START_REF] Boissonnat | Triangulations in CGAL[END_REF] and it can deal with constraints which allow to delimit domains. The mesh generation algorithms add Steiner points on constraints or inside domains such that triangles fulfill user defined quality criteria as minimal angles, or triangle size. The triangle size can be constant, or it can be governed by a sizing field, for example local feature size.

3D Triangulations� Surface and Volume Mesh Generation : The generated meshes are three-dimensional isotropic simplicial meshes [START_REF] Boissonnat | Provably good sampling and meshing of surfaces[END_REF][START_REF] Rineau | A generic software design for Delaunay refinement meshing[END_REF]. The discretized region may be a monodomain, i.e. a single connected component, or a multidomain subdivided into different subdomains. The domain is input as an oracle able to answer queries, of a few different types, on the domain. In the current version, the mesh generator is able to discretize domains described through implicit functions, 3D images or polyhedral boundaries. The output is a 3D mesh of the domain volume and conformal surface meshes for all the boundary and subdividing surfaces.

The algorithm can also respect one-dimensional features formed by intersections of surfaces, or sharp features of objects. The package further provides mesh optimization algorithms to improve the aspect ratio of tetrahedra [START_REF] Tournois | Interleaving Delaunay refinement and optimization for practical isotropic tetrahedron mesh generation[END_REF][START_REF] Cheng | Sliver exudation[END_REF].

Periodic Triangulations : This package allows to build and handle triangulations of infinite periodic point sets in 3D �see Figure 4), which can also be seen as triangulations of the three-dimensional flat torus [START_REF] Caroli | Computing 3D periodic triangulations[END_REF][START_REF] Caroli | Delaunay triangulations of point sets in closed euclidean d-manifolds[END_REF]. Triangulations are built incrementally and can be modified by insertion or removal of vertices. The package provides Delaunay triangulations and offers nearest neighbor queries and primitives to build the dual Voronoi diagrams.

FIGURE 4 -3D periodic Delaunay triangulation

Some work is in progress regarding periodic surface and volume mesh generation �Figure 5) [START_REF] Caroli | Meshing of triply-periodic surfaces in CGAL[END_REF]. Surface Reconstruction : This package allows to take a point cloud coming for example from a laser range scanner. Given a set of 3D points with oriented normals sampled on the boundary of a 3D solid, the Poisson surface reconstruction method [START_REF] Kazhdan | Poisson Surface Reconstruction[END_REF] solves for an approximate indicator function of the inferred solid, whose gradient best matches the input normals. The output scalar function, represented in a refined Delaunay triangulation, is then iso-contoured using the Surface mesh generator �see Figure 6).

CGAL offers tools to clean the input which is often noisy, has outliers, and must often be locally simplified, and tools to estimate normals.

Polyhedral Surfaces : CGAL offers algorithms that operate on polyhedral surfaces, such as mesh simplification [START_REF] Lindstrom | Fast and memory efficient polygonal simplification[END_REF], estimation of local differential properties, principal component analysis, as well as various surface subdivision schemes. The polyhedral surface are represented by a halfedge data structure. As they are limited to manifold surfaces, efforts are under way to implement combinatorial maps [START_REF] Damiand | Contributions aux Cartes Combinatoires et Cartes Généralisées : Simplification, Modèles, Invariants Topologiques et Applications[END_REF].

Boolean Operations : The 2D Boolean operations package is based on overlays of arrangements of curves [START_REF] Wein | Advanced programming techniques applied to CGAL's arrangement package[END_REF]. Polygon edges can be line segments, as well as circular arcs or Bezier curves. The former appear often in 2D CAD data, the latter in descriptions of fonts.

The 3D Boolean operations package is based on the theory of the Nef polyhedra, which can represent the result of any finite sequence of Boolean operations applied on half spaces [START_REF] Bieri | Nef polyhedra : A brief introduction[END_REF]. This allows open faces, antennas, isolated points, or volumes sharing only an edge.

Polygon Offsets : CGAL offers two variants for computing inner or outer offsets of a polygon. The first one is based on Minkowski sums which results in round caps at corners, and the second one is based on the straight skeleton which is a polygon with only straight line segments �see Figure 7). Some users simplify polygons by computing the inner offset of the outer offset of the input polygon. Neighbor Search : CGAL offers data structures for computing k-nearest/furthest neighbors for points in arbitrary dimensions. There are variants for approximative queries, and for various metrics. These data structures are based on the kd-tree.

CGAL further offers data structures to perform efficient intersection and distance queries against sets of finite 3D geometric objects. The set of geometric objects stored in the data structure can be queried for intersection detection, intersection computation and distance. The intersection queries can be of any type, provided that the corresponding intersection predicates and constructors are implemented in the traits class. The distance queries are limited to point queries. Examples of intersection queries include line objects �rays, lines, segments) against sets of triangles, or plane objects �planes, triangles) against sets of segments. An example of a distance query consists of finding the closest point from a point query to a set of triangles. The underlying data structure is an AABB tree [START_REF] Terdiman | OPCODE 3D Collision Detection library[END_REF].

Bounding Areas : This package offers algorithms to compute the smallest enclosing circle, ellipse, rectangle, quadrilateral, or strip for a given point set. It further offers the convex hull algorithm.

Voronoi Diagrams : For a given set of sites, a Voronoi diagram decomposes space in cells, generally one cell per site. CGAL offers Voronoi diagrams for 2D points, segments, and circles. The Voronoi edges are then arcs of lines, parabolas, and hyperbolas. The latter are dealt with symbolically, that is they are not discretized polylines. The segment Voronoi diagram is often used for computing the medial axis of a polygon.

Polyline Simplification : This package is work in progress that allows to simplify several polylines simultaneously while maintaining the topology. The simplification of a coastline will not produce an island which is no longer surrounded by water, and isolines will not touch each other, if they did not touch before the simplification.

Software issues

The design of CGAL [START_REF] Fabri | On the design of CGAL, a computational geometry algorithms library[END_REF][START_REF] Fogel | Generic programming and the CGAL library[END_REF] follows the generic programming and the exact geometric computing paradigm.

Generic Programming

Just as a std::set of the Standard Template Library, can store objects of an arbitrary type as long as a comparison operator for this type is defined, a CGAL data structure that deals with points can operate on an arbitrary point type as long as several geometric predicates like lexicographical order, and incircle test are defined for this point type. The types and predicates that must be provided are part of the documentation for each CGAL data structure. In the same spirit CGAL provides the models of the Boost Graph Library �BGL) boost::graph concept, which allows to run any graph algorithm of the BGL directly on CGAL data structures.

Generic Programming [START_REF] Austern | Generic Programming and the STL[END_REF], with the notion of concepts defining an API and models implementing it, is key to reuse, and to the ease of integration of CGAL components into existing applications. This is important, as few users start from scratch and build entire applications only on CGAL.

Exact geometric computing

Geometric algorithms are notorious for their lack of robustness when faced with numerical instability issues. In order to deal with this problem, CGAL follows the Exact Geometric Computing paradigm [START_REF] Li | Recent progress in exact geometric computation. Special issue on the practical development of exact real number computation[END_REF]. This paradigm states that computing the low-level primitives exactly guarantees the correctness of the higher-level algorithms.

These primitives are categorized in two types. On one hand, geometric predicates compute Booleanlike properties, like the orientation of three points in the plane, or the decision whether two segments intersect. On the other hand, geometric constructions build new geometric objects or numerical quantities, like the computation of a distance or the intersection point of two lines.

Language bindings

CGAL is written in C++, but there is work in progress on offering Java and Python bindings based on Swig. The 2D and 3D triangulations of MATLAB are wrapped version of the CGAL classes. 1

License

CGAL is distributed under an Open Source which makes that the library has a large user community. As the QPL license does not allow closed source development, the library is at the same time available under commercial licenses. Commercial users include companies as Agilent, Total, ESRI, The Mathworks, Schlumberger, Navteq, to name just a few.

FIGURE 1 -

 1 FIGURE 1 -The Delaunay surface mesh has nicely formed triangles and less triangles where the surface is almost flat, compared to the mesh generated with marching cubes.

FIGURE 2 -

 2 FIGURE 2 -Volume mesh for 3D segmented voxel data.

FIGURE 3 -

 3 FIGURE 3 -Volume mesh and remeshed surface mesh for a mechanical part.

FIGURE 5 -

 5 FIGURE 5 -Periodic surface �left) and volume �right) meshes. �left) Data courtesy of Maarten Moesen, Department of Metallurgy and Materials Engineering �MTM), K.U. Leuven, Belgium.

FIGURE 6 -

 6 FIGURE 6 -The Poisson surface recinstruction algorithm copes well with undersampled areas.

FIGURE 7 -

 7 FIGURE 7 -Offsets, based on : Minkowski sum �left) and straight skeleton �right)

video : http://www.mathworks.com/products/demos/shipping/matlab/NewMATLABMathematicsFeaturesinR2009a.html

http://www.cgal.org/people.html

Acknowledgments

Our thanks to the CGAL developers 2 who wrote the demos and made some of the pictures.


L Rhode-Barbarigos
email: giosef.rhode-barbarigos@epfl.ch

C Schulin

N Bel Hadj

R Motro
email: rene.motro@univ-montp2.fr

I F C Smith

9-13 mai 2011, Presqu'île de Giens (Var) Deployment analysis of a pentagonal tensegrity-ring module

Keywords: Tensegrity structures, deployable structures, dynamic relaxation

Ring modules are tensegrity systems that include a single strut circuit and recently, they have been shown to be viable systems for pedestrian bridges. Furthermore, their shape can be controlled using cable actuation. This paper focuses on the deployment of a pentagonal tensegrity-ring module. A geometric study is conducted to identify the deployment-solution space without strut contact. Deployment paths and actuation requirements are explored. The structural response of the module during deployment is analyzed using a modified dynamic relaxation method.

Introduction

Tensegrity structures are spatial reticulated structures that are composed of cables and struts. Therefore, they contain only axially loaded elements. The stability of tensegrity structures is based on the equilibrium among cables and struts. Since the introduction of the tensegrity concept in the 1950s, there have been many definitions for tensegrity systems [START_REF] Motro | Tenségrité: Hermes Science[END_REF]. A widely accepted definition was proposed by Motro [START_REF] Motro | Tenségrité: Hermes Science[END_REF] in 2003: "A tensegrity is a system in stable self-equilibrated state comprising a discontinuous set of compressed components inside a continuum of tensioned components". This definition includes systems with interconnected compressed elements into the broader class of tensegrity structures. Moreover, a class characterization was proposed to distinguished tensegrity systems that contain interconnected compressed elements: "class I" tensegrity systems have no connection among struts, while "class II" systems include strut-to-strut connections [START_REF] Skelton | An introduction to the mechanics of tensegrity structures[END_REF]. Throughout the four last decades, scientists and engineers have shown interest in the tensegrity concept [START_REF] Rovira | Control and simulation of a tensegrity-based mobile robot[END_REF][START_REF] St | Using Actuated Tensegrity Structures to Produce a Responsive Architecture[END_REF][START_REF] Ingber | The architecture of life[END_REF]. Although tensegrity systems can be combined to create lightweight modular structures with high strength to mass ratios, few examples of tensegrity structures have been used for civil engineering applications. The Rostock tower in Germany is an assembly of six "simplex" modules designed by Schlaich [START_REF] Schlaich | The messeturm in Rostock -A tensegrity tower[END_REF]. The structural system of the roof of the World Cycling Center Arena in Aigle (Switzerland) is composed of a tensegrity structure [START_REF] Paronesso | The 2002 World Cycling Center Arena Aigle[END_REF]. Furthermore, a tensegrity grid was proposed for the Swiss National Exhibition in 2002 [START_REF] Paronesso | The cloud of Yverdon[END_REF]. Research on double layer tensegrity grids resulted in the construction of an 81 m 2 double layer tensegrity grid [START_REF] Quirant | Designing tensegrity systems: the case of a double layer grid[END_REF]. Tensegrity modules are attractive systems not only for conventional structures but also for active/adaptive structures and deployable structures. Active structures can interact with complex environments with the use of embedded actuated elements. In tensegrity systems, actuated elements and structural elements can be combined. Consequently, both actuated cables and telescopic struts can be used for active control. Furthermore, tensegrity systems require small amounts of energy for actuation and control [START_REF] Skelton | An Introduction to the Mechanics of Tensegrity Structures[END_REF]. A modular full-scale active tensegrity structure with telescopic struts was designed and constructed by Fest [START_REF] Fest | Active Tensegrity Structure[END_REF]. Active control was enhanced with advanced informatics such as stochastic search and case-based reasoning in order to improve service performance and perform self-diagnosis as well as self-repair [START_REF] Domer | An Active Structure that Learns[END_REF][START_REF] Adam | Reinforcement Learning for Structural Control[END_REF]. Moreover, vibration control was also studied and verified experimentally [START_REF] Hadj | Dynamic behavior and vibration control of a tensegrity structure[END_REF]. Deployable structures change shape from a compact configuration to an expanded one. Actuation can thus have a dual action: control the shape of the structure and enhance service performance after deployment. Furuya [START_REF] Furuya | Concept of deployable tensegrity structures in space applications[END_REF] studied the deployment of a tensegrity mast. Tibert and Pellegrino [START_REF] Tibert | Deployable Tensegrity Masts[END_REF] compared deployable tensegrity masts to conventional ones. Lack of stiffness during deployment and low bending stiffness were identified as obstacles to practical applications for deployable tensegrity masts. Sultan and Skelton [START_REF] Sultan | Deployment of tensegrity structures[END_REF] studied the deployment of tensegrity systems and showed that cable control has the advantage of keeping the system close to an equilibrium manifold throughout deployment. The system can thus maintain stiffness during deployment. Bouderbala and Motro [START_REF] Bouderbala | Folding tensegrity systems[END_REF] studied folding of "expanded octahedron" assemblies with cable and strut actuation. Smaili and Motro [START_REF] Smaili | Folding/unfolding of tensegrity systems by removal of self-stress[END_REF] investigated folding of tensegrity systems using finite mechanisms. Similarly, Sultan [START_REF] Sultan | Tensegrity motion control using internal mechanisms[END_REF] studied shape control of tensegrity structures through infinitesimal mechanisms. Moreover, Schenk et al. [START_REF] Schenk | Zero stiffness tensegrity structures[END_REF] studied the application of statically balanced mechanisms in tensegrity systems. Although cable actuation has been used for deployment, no study explored the use of continuous cables. Furthermore, most of the studies do not focus on civil engineering applications. This paper focuses on the deployment of a pentagonal tensegrity-ring module using actuated cables. The deployment-solution space is identified through a geometrical study. Deployment paths are compared based on the number of actuators required. Furthermore, continuous cables and spring elements are used to decrease the number of actuators. The structural behavior of the modules during deployment is investigated using a deployment-analysis algorithm.

FIG. 1 -Illustration of a pentagonal tensegrity-ring module including struts, layer cables and x-cables 2 Pentagonal tensegrity-ring module topology

Tensegrity-ring modules were presented in 2006 by Motro [START_REF] Motro | Tensegrity Rings and the Hollow Rope[END_REF]. They are named after their topology that includes interconnected struts in a single circuit in the form of a ring (Figure 1). Motro studied their topology and their deployability [START_REF] Motro | Tensegrity Rings and the Hollow Rope[END_REF]. Ring-module topology is based on straight prism geometry. Moreover, straight prism geometry can be used to further characterize the geometry of the module. The authors showed that pentagonal ring modules are viable systems for a pedestrian bridge [START_REF] Rhode-Barbarigos | Designing tensegrity modules for pedestrian bridges[END_REF][START_REF] Hadj Ali | Design optimization and dynamic analysis of a tensegrity-based footbridge[END_REF]. A single pentagonal ring-module is shown in Figure 1. The module includes 15 struts and 30 cables. There are two sets of cables: layer cables and x-cables. Layer cables form a pentagon on the front and back side of the module, while x-cables form an "x" on each lateral side of the module. Ring-module topology guarantees deployability with cable actuation [START_REF] Motro | Tensegrity Rings and the Hollow Rope[END_REF]. If cables can change rest length, then modules can change shape (deploy) from a folded to an unfolded configuration. Dung studied the structural response of the deployed configuration of a pentagonal ring-module and its folding [START_REF] Nguyen | Etude de comportement mécanique et du pliage d'un anneau de tenségrité à base pentagonale[END_REF]. Cable stiffness was found to be more important for the overall stiffness compared with strut stiffness. Folding was simulated using FEM where the structure is compacted under nodal external forces leading to some cables to go slack. Therefore, this method is not suitable for studying deployment under cable actuation.

Geometrical study of deployment and cable actuation

The motion of ring modules during deployment is composed of a translation, a rotation and a dilation of their pentagonal faces. The translation between front and back pentagonal faces defines the deployment length L. The transverse rotation θ of the two faces and their dilation given by the external radius R describe the repositioning of the strut circuit under the new deployment length. Both R and θ increase with folding and decrease with unfolding. Consequently, during deployment ring-module topology can be described based on the module length L, the external radius R and the transverse rotation θ between front and back pentagonal face of the module (Figure 2). However, there is no explicit relationship among the parameters since various configurations have the same length. If module topology is respected, R can be expressed in terms of L and θ. Consequently, the deployment space is reduced into a 2D-space described by L and θ. Figure 2 shows the deployment-solution space for the pentagonal module with a total length of 75 cm with an external radius of approximately 57 cm. Struts have a length of 100 cm. Cable length is set at approximately 68 cm and 51 cm for layer cables and x-cables respectively. The space in white describes the space where deployment is conducted without strut contact. Isometric curves in Figure 2 represent the closest distance between struts from 1 (inner curve) to 9 cm (outer curve).

FIG. 2 -Deployment-solution space and studied deployment paths

A deployment path is defined based on the tip of isometric curves along which the risk of strut contact is minimized (maximum strut-to-strut distance). The same path is followed for both unfolding and folding (Figure 2). This path requires the actuation of all cables in the pentagonal ring module (10 layer cables and 20 x-cables). Layer cable length is decreased for unfolding while x-cable length is

θ 2π/5 -θ z y y x z R R L
increased. Inversed actions take place for folding. Consequently, layer cable action controls unfolding while x-cable action controls folding. However, actions in both types of cables are required to obtain new stable configurations.

A deployment path that reduces the number of actuated cables form 30 to 20 was found. Along this path, x-cables that are coplanar with struts do not require actuation as their length remains constant. Layer cables remain actuated to control unfolding. This path is similar to the path with minimum risk of strut contact (Figure 2). A lower transverse rotation and a slightly larger folded length are the main differences.

In order to further reduce the number of actuators, continuous cables and spring elements are used in the module. Continuous cables affect the mechanics of tensegrity systems by changing the number of independent states of self-stress and the number of infinitesimal mechanisms in the module [START_REF] Moored | Investigation of clustered actuation in tensegrity structures[END_REF]. The use of continuous cables in all actuated cables (layer cables and x-cables) results in an unstable configuration. Therefore, continuous cables are used only in actuated x-cables found on each lateral side of the module. The total number of actuators decreases thus to 15: 10 actuators in layer cables and 5 actuators in x-cables. Furthermore, spring elements can be used to replace actuated layer cables. Spring elements allow length changes without requiring actuation. Both deployment phases are thus controlled by x-cables. The configuration including continuous cables in actuated x-cables and spring elements for layer cables requires only 5 actuators for deployment.

FIG. 3 -Snapshots of the deployment of small scale physical model of the pentagonal ring-module

A small scale physical model of the pentagonal ring module was constructed for experimental validation of the deployment path (Figure 3). Therefore, the module was designed according to the ring-module topology used in the geometrical study. The module is composed of aluminum struts and steel cables. Joints allow the required movement for deployment. The minimum actuation configuration is used for the control of the model. Continuous cables are thus used in actuated x-cables and spring elements are used to replace layer cables. Cable actuation is conducted manually using cranks placed on the support of the model. Both unfolding and folding were conducted successfully without strut contact. However, a systematic error in the transverse rotation is observed due to joint design.

Structural analysis of the deployment

A deployment-analysis algorithm was conceived for the structural analysis of tensegrity modules during deployment. The algorithm is based on a modified version of the dynamic relaxation method. Dynamic relaxation is a method suitable for the analysis of highly non-linear structures such as tensegrity systems. The modified version of the dynamic relaxation accommodates continuous cables [START_REF] Hadj | Analysis of clustered tensegrity structures using a modified dynamic relaxation algorithm[END_REF]. Furthermore, the method was adapted to take into account spring elements. Structural analysis during unfolding and folding is conducted with self-weight as unique loading. Friction and dynamic effects are not part of this study. Moreover, boundary nodes should allow deployment while blocking rigid body movements.

The deployment-analysis algorithm integrates the ring-module actuation scheme. The algorithm starts from a known folded or unfolded topology. Cable actuation is used to create a finite mechanism that will allow the module to change length. Increasing the length of x-cables leads to unfolding the module. The inversed action is used for folding the module. Length changes in spring elements (replacing actuated layer cables) are driven by cable actuation. Increasing x-cable length results in a contraction of spring elements corresponding to a decrement in the length of layer cables. After each actuation, a new equilibrium configuration is found based on the modified dynamic relaxation method. Equilibrium configurations are checked for strut contact. Deployment is thus based on a series of sequences of "mechanism -equilibrium -constraint". If strut contact is detected then the algorithm provides the details of the event. If deployment is completed successfully, then the deployment-path and the internal forces during unfolding and folding are given. This algorithm is summarized in Figure 4.

FIG. 4 -Illustration of the deployment-analysis algorithm

The deployment of a module made of aluminum struts and steel cables is analyzed using the deployment-analysis algorithm. Strut and cable cross-sections are set at 2.54 cm 2 and 0.03 cm 2 respectively. Spring elements have a constant of 20 N/cm and an initial length of approximately 66 cm. Figure 5 shows the evolution of internal forces during unfolding and folding for a single module. Structural elements are separated in actuated x-cables, spring elements (actuated layer cables), nonactuated x-cables and struts. Two types of struts are identified based on their internal forces. Strut types are in relation with ring-module topology. Figure 5 reveals that intermediate equilibrium configurations for unfolding and folding are the same. Internal forces in all structural elements increase with folding and decrease with unfolding. This trend is an effect of spring action as spring elements elongate with folding and contract with unfolding. Internal forces in actuated cables are lower compared with non-actuated elements. Finally, no strut contact is observed during unfolding or folding.

FIG. 5 -Internal force evolution in the pentagonal ring module during deployment

Discussion

The deployment of a pentagonal ring module is analyzed in this paper. Although deployment motion can be described using three geometrical parameters: L, θ and R (translation, rotation and dilation of the pentagonal faces), no explicit relationship between the parameters exists for deployment. The deployment-solution space is defined based on avoiding strut contact. A deployment path with minimum contact risk was identified. The path requires the actuation of all 30 cables for a single module. Furthermore, a deployment path with only 20 actuators was found based on observations on the path with minimum contact risk. The number of actuators is further reduced to 5 with the use of continuous cables and spring elements. Deployment was successfully validated without strut contact in a small scale physical model including both continuous cables and spring elements. The deployment of a pentagonal ring module controlled by the proposed actuation scheme is composed of a series of equilibrium configurations. Furthermore, deployment motion is entirely controlled by cable actuation. Large actuation steps in x-cables result in unstable configurations while small actuation steps are computationally expensive. Moreover, spring elements are important for maintaining ring-module topology during deployment. Spring characteristics affect the stability of the system. Low values of spring constant combined with large values of initial spring length do not provide sufficient stiffness to the system resulting in unstable configurations. On the contrary, large values of spring constant combined with low values of initial spring length induce large internal forces that may exceed element strength during deployment. Therefore, actuation-step size and spring characteristics are important parameters for deployment. Spring elements also affect the energy in the tensegrity system. Spring elements contract during unfolding and elongate during folding. It was found that the energy stored in spring elements during folding can be used for unfolding the module. Consequently, only folding requires an external energy supply. This finding was validated experimentally on the small scale physical model. Finally, a single actuation step and a single spring constant were applied in actuated elements for this study. Different actuation steps and different spring constants could be considered for each actuated element. This may result in better control of the structure during deployment.

Conclusion

The conclusions from this study are as follows:

• The deployment of pentagonal ring modules is feasible through cable control.

• Continuous cables and spring elements can be used to reduce the number of actuated cables.

• Actuation steps and spring constants are critical parameters for deployment.

• Energy stored in spring elements during folding can be used for unfolding.

Work in progress includes the experimental validation of structural behavior of the module during deployment. Furthermore, a prototype of the tensegrity bridge will be built and studied experimentally.


Amdouni Saber

Khalil Mansouri
email: khalil.mansouri@insa-lyon.fr

Yves Renard
email: yves.renard@insa-lyon.fr

Makrem Arfaoui
email: makremarfaoui@yahoo.fr

Maher Moakher
email: maher.moakher@enit.rnu.tn

S Amdouni
email: amdouni@insa-lyon.fr

Numerical convergence and stability of mixed

come

Introduction

The Finite Element Methods (F.E.M.) are widely used in engineering analysis, and has deep mathematical foundations. However, the treatment of incompressible or nearly incompressible problems, with only one unknown field variable, leads to numerical locking of this approximation. For the analysis of such problems, it is necessary to use a mixed formulation. However, mixed finite element methods are not stable in all cases, some of them show spurious pressure oscillations if displacement and pressure spaces are not properly chosen [1]. For a mixed formulation to be stable it must verify consistency, ellipticity and the so called inf-sup (or LBB) condition. The later is a severe condition which depends on the connection between the approximation spaces for the displacement and pressure [START_REF] Chapelle | The inf-sup test[END_REF]. The presence of a crack in a structure reveals two types of discontinuities : a strong discontinuity that requires an adapted mesh to the shape of the crack, hence the domain is meshed at each time step ; and a weak discontinuity that requires refinement at the crack tip. These two operations lead to a huge computational cost. In order to overcome these difficulties and to make the finite element method more flexible, many alternative methods had been developed. Apart from mixing meshless methods and finite elements, another alternative consists in using the concept of "partition of unity" introduced by Babuška, to creates an independence between the mesh and the crack. Among the class of partition of unity finite element methods, the Generalised Finite Element Method (GFEM) and the eXtended Finite Element Method (X-FEM) are the most advanced. The X-FEM allows to model cracks, material inclusions and holes on nonconforming meshes. It was introduced by Moës and al [START_REF] Moës | A finite element method for crack growth without remeshing[END_REF]. It consists in enriching the basis of the classical finite element method by a step function along the crack line and by some non-smooth functions representing the asymptotic displacement around the crack tip. In spite of the singular enrichment of the finite element basis, the obtained convergence error of X-FEM remains of order if linear finites elements are used. To improve this convergent rate and to obtain an optimal accuracy, Chahine and al introduced a new enrichment strategy [START_REF] Chahine | Crack tip enrichment in the X-FEM method using a cut-off function[END_REF] : the so called X-FEM cut-off. This enrichment strategy uses a cut-off function to locate the crack tip surface. In their work, Chahine et al have shown that the X-FEM cut-off has an optimal convergence rate of order h and that the conditioning of the stiffness matrix does not deteriorate. In this work, we extend the numerical results given by Chahine et al. [START_REF] Chahine | Crack tip enrichment in the X-FEM method using a cut-off function[END_REF] to an incompressible isotropic linear plane elasticity problem in fracture mechanics. In particular, this formulation must satisfy the famous inf-sup condition (or LBB). Let Ω be a two-dimensional cracked domain, Γ c denotes the crack and Γ the boundary of Ω. We assume that Γ \ Γ c is partitioned in two parts : Γ N where a Neumann surface force t is applied and Γ D where a Dirichlet condition u = 0 is imposed (see Fig. 1). We assume that we have a traction-free condition on Γ c . Let f be the body force applied on Ω. The equilibrium equation and boundary conditions are given by

-div σ(u) = f, in Ω, (1) σ(u) = C : ε(u), in Ω, (2)
u = 0, on Γ D , (3) σ
(u) • n = t, on Γ N , (4) σ
(u) • n = 0, on Γ c . ( 5
)
with n is the outside normal to the domain Ω. Let p be the hydrostatic pressure defined in two dimensions by : p = -tr(σ) 2 .

Now we decompose the stress tensor σ in two parts : the spherical part and the deviatoric part σ d given by :

σ d (u) = σ(u) + p I = 2µε d (u), where ε d (u) = ε(u) - div(u) 2 I.
Then, the strong mixed formulation is written as follows

-div[σ d -p I] = f in Ω, (6)
div u = 0 in Ω. ( 7
) Let V = v ∈ H 1 (Ω) with u = 0 on Γ D and Q = L 2 (Ω).
Taking the inner poduct of (6) with a test function v ∈ V , and multiplying (6) by a test function q ∈ Q. On applying Green's formula for elasticity, we find the weak mixed formulation

   Find (u, p) ∈ (V, Q) such that : Ω σ d (u) : ε(v) dΩ -Ω p div v dΩ = Ω f • v dΩ + Γ N t • v dΓ, ∀v ∈ V, Ω q div u dΩ = 0, ∀q ∈ Q.
Accordingly, the weak mixed formulation of the isotropic incompressible linear elastic problem is written as :

   Find (u, p) ∈ (V, Q) such that : a(u, v) -b(v, p) = L(v), ∀v ∈ V, b(u, q) = 0, ∀q ∈ Q, (8) with a(u, v) = Ω σ d (u) : ε(v) dΩ, b(v, p) = Ω p div v dΩ, L(v) = Ω f • v dΩ + Γ N t • v dΓ.
Discretization of the elasticity problem follows the usual steps. We approximate (u, p) by (u h , p h ) ∈ V h × Q h . The subspaces V h and Q h are finite dimensional spaces that will be defined later. The discretized problem is then :

   Find (u h , p h ) ∈ (V h , Q h ) such that a(u h , v h ) -b(v h , p h ) = L(v h ), ∀v h ∈ V h , b(u h , q h ) = 0, ∀q h ∈ Q h . (9)
REMARK : The existence of a stable finite element approximate solution (u h , p h ) depends on choosing a pair of spaces V h and Q h such that the following condition holds :

inf q∈Q h sup v∈V h b(v h , q h ) q h 0,Ω v h 1,Ω ≥ β 0 .
where β 0 > 0 independent of h. This is the "Ladyzhenskaya-Brezzi-Babuška condition" (LBB) also called the "inf-sup condition" [1]. The verification of this condition for a couple (V h , Q h ) is very difficult to prove in practical situations. Therefore, the numerical evaluation of the inf-sup has been widely studied (see Chapelle and Bathe [START_REF] Chapelle | The inf-sup test[END_REF]). The numerical evaluation gives an indication of the verification of the LBB condition for a given finite element discretization. The numerical inf-sup test is based on the following proposition.

Proposition 1 ( [START_REF] Chapelle | The inf-sup test[END_REF]) Let [M] uu and [M] pp be the matrices associated with the H 1 -inner product in V h and the L 2 -inner product in Q h , respectively, and let µ min be the smallest nonzero eigenvalue of the eigenvalue problem :

[B] T [M] -1 uu [B] V = µ 2 [M] pp V.
Then the value of β 0 in the LBB condition is simply µ min .

The numerical test proposed in [START_REF] Chapelle | The inf-sup test[END_REF] is to test the stability of the mixed formulation by calculating β h with increasingly refined meshes. Indeed, if log(β h ) continuously decreases as h goes to 0, Chapelle and Bathe [START_REF] Chapelle | The inf-sup test[END_REF] predicted that the finite element violates the LBB condition. Otherwise, if β h stabilizes when the number of elements increases, then the numerical inf-sup test is verified.

X-FEM cut off

The idea of XFEM is to use a classical finite element space enriched by some additional functions. These functions result from the product of global enrichment functions and some classical finite element functions. This variant of X-FEM uses a cut-off function to find the singular enrichment surface. Indeed, for this method, the enrichment is introduced on all crack tip surface and not only on the nodes. The classical enrichment strategy for this problem is to use the asymptotic expansion of the displacement and pressure fields at the crack tip area. Indeed, the displacement is enriched by the Westergaard functions.

F = F j (x), 1 ≤ j ≤ 4 = √ r sin θ 2 , √ r cos θ 2 , √ r sin θ 2 sin θ, √ r cos θ 2 sin θ .
For the pressure, the asymptotic expansion at the crack tip is given by [START_REF] Legrain | Stability of incompressible formulations enriched with X-FEM[END_REF] p(r, θ)

= 2K I 3 √ 2πr cos θ 2 + 2K II 3 √ 2πr sin θ 2 .
This expression is used to obtain the basis of enrichment of the pressure in the area of the crack tip [START_REF] Legrain | Stability of incompressible formulations enriched with X-FEM[END_REF] :

     F p 1 = 1 √ r cos θ 2 , F p 2 = 1 √ r sin θ 2 .
We note that the displacement and pressure are also enriched with a Heaviside function at the nodes whose shape functions support is totally cut by the crack. Using this enrichment strategy, the discretisation spaces V h and Q h take the following forms :

V h = v h = ∑ i∈I α k ψ u,k + ∑ i∈I H β k Hψ u,k + 4 ∑ j=1 γ j F u j χ; α k , β k , γ j ∈ R ∀i, j , Q h = p h = ∑ i∈I p i ϕ p,i + ∑ i∈I H b p i Hϕ p,i + 2 ∑ j=1 c p j F p j χ; p i , b p i , c p j ∈ R ∀i, j ,
with ψ u,k is the vector shape function defined by :

ψ u,k =            ϕ u,i 0 if i = k + 1 2 , 0 ϕ u,i if i = k 2 , (10)
and ϕ u,i (resp.ϕ p,i ) is a scalar schape function for displacement (resp. for pressure). χ is a is a polynom of order 5 between r 0 and R 1 such that :

χ(r) = 1 if r < r 0 , χ(r) = 0 if r > r 1 . ( 11
)
4 Numerical study

In this section we numerically study the inf-sup condition and a comparison between the convergence rates of the X-FEM fixed area and X-FEM cut off. The numerical test is made on a non-cracked domain The Von Mises stress for this test is presented in Fig. 3(b). As expected the Von Mises stress is concentrated at the crack tip.

Numerical inf-sup test

In this section we numerically study the inf-sup condition and its dependence on the position of the crack. First, the inf-sup condition is approximated using gradually refined structured triangulation meshes. The evolution of the numerical inf-sup value is plotted in Fig. 4 with respect to the element size. From this figure we can conclude that the numerical inf-sup value is stable for all studied formulations. Let δ be the crack position as shown in Fig. 5. To test the influence of the position of the crack on the inf-sup condition, we calculate the LBB condition by decreasing δ. The tests are made, on a P + 1 /P 1 formulation, with h = 1/60 (see Fig. 6(a)) and h = 1/40 (see Fig. 6(b)). The results presented on Figs. 6(a) and 6(b) show that the inf-sup condition remains bounded regardless of the position of the crack. Hence, one can conclude that our formulation is stable independently of the position of the crack. Figure 7(b) confirms that the convergence rate for the energy norm is of order h for both variants, of the X-FEM : with fixed area and cut-off. Figure 7(a) shows that the convergence rates for the L 2 -norm in displacement is of order h 2 for both variants. Figure 7(c) shows that the convergence rates for the L 2 -norm in pressure is h for both variants. Furthermore, the rate of convergence for the L 2 -norm is better than the one for the H 1 -norm, which is quite usual. Compared to the X-FEM method with a fixed enrichment area, the convergence rate is very close but the In order to test the computational cost of X-FEM cut-off, Table (1) shows a comparison between the number of degrees of freedom for different refinements of the classical method X-FEM with fixed enrichment area and the cut off method. This latter enrichment leads to a significant decrease in the number of degrees of freedom. Figure 8 shows that the conditioning of the linear system associated to the cut-off enrichment is much better than the one associated with the X-FEM with a fixed enrichment area. We can conclude that, similarly to the X-FEM with fixed enrichment area, the X-FEM cut off leads FIGURE 8 -Conditioning number of the stiffness matrix for the mixed problem. to an optimal convergence rate but it reduces the approximation errors without significant additional costs.

Convergence of the High-order X-FEM

The numerical tests of the higher order X-FEM method (P + 2 /P 1 disc, P + 2 /P 1 , P 2 /P 1 and P 2 /P 0 ) do not give an optimal order of convergence (see Figs. 9(a), 9(b), 9(c) and 9(d)). This means that the enrichment function does not capture the behavior of the solution at the crack tip. This result was expected as the main singularity belongs to H 3/2 (Ω). Then, for the X-FEM cut-off, the convergence rate remains limited to h 3/2 with high order polynomials.

To have an optimal convergence rate, we must make an asymptotic expansion of order 2 to find the correct expression of the enrichment basis for the displacement and pressure.

conclusion

By this study we can conclude that the X-FEM cut-off mixed formulation is stable, regardless of the position of the crack. Similar to the X-FEM with fixed enrichment area, the X-FEM cut-off gives an optimal convergence rate but without significant additional costs. For shape functions of higher order, the convergence rate is limited to h 3/2 . This result was expected as the main singularity belongs to H 3/2 (Ω).

Références

[1] F. Brezzi and M. Fortin. Mixed and Hybrid Finite Element Methods, volume 15. Springer Series in Computational Mathematics, New York, 1991.

FIGURE 1 -Cracked domain 2

 2 FIGURE 1 -Cracked domain

FIGURE 2 -

 2 FIGURE 2 -Cracked specimen

FIGURE 3 -FIGURE 4 -

 34 FIGURE 3 -Structured triangulation and Von Mises stress

FIGURE 5 -

 5 FIGURE 5 -Position of the crack.

  FIGURE 6 -Evolution of the inf-sup condition as, a function of the position of the crack

(a) Enriched P + 2 /

 2 FIGURE 9 -Error for the mixed problem with enriched P + 1 /P 1 elements (logarithmic scales)


M Papadrakakis

Mastering Computational Demanding Problems in Mechanics with Neural Network Predictions

The solution of Computational Mechanics problems were traditionally based on rigorous scientific procedures that are formed on mathematical methods and the principles of theoretical mechanics and are solved with the implementation of numerical simulation methods based on discretized continua. However, three decades ago new families of computational methods, denoted as soft computing (SC) methods, have emerged. These methods are based on heuristic approaches rather than on rigorous mathematics. Despite the fact that these methods were initially received with suspicion, their use in various areas of engineering science is continuously growing. Artificial Neural Networks (ANN), Genetic Algorithms and Fuzzy Logic are the most popular approaches of SC. An ANN can store experimental knowledge and make it available for later use. Its features are adaptive learning, self-organizing capability during training and fault imprecision during applications. The main advantage of using ANN is that they can deal with problems that do not have an algorithmic solution or for which an algorithmic solution is too complex to be found. Among other fields of applications, ANN have been used for the identification of nonlinear dynamic systems and damage assessment. Predictions of the structural behaviour by neural networks have been employed in the context of design optimization and structural reliability analysis. Furthermore, ANN have been applied for the prediction of systems response under static or dynamic loading conditions. Performance-Based Design (PBD) is the current trend for the design of structural systems where the structural performance is assessed for multiple hazard levels. PBD framework is performed by means of either a static or a dynamic analysis. In the latter case multiple dynamic analyses have to be performed for every hazard level requiring significant computational effort. Structural reliability analysis can be performed either with simulation methods, such as the Monte Carlo Simulation (MCS), or with other approximation methods. Although the major advantage of MCS is that accurate solutions can be obtained for almost every problem, yet it requires excessive computational effort in many cases. Variance reduction techniques, such as Importance Sampling, Directional Simulation, Antithetic Variates or Adaptive Sampling, have been proposed in order to reduce the computational effort of MCS. However, recent studies reveal that variance reduction techniques still require significant number of system's response evaluations in order to estimate failure probabilities of low order. In this work the capabilities and limitations of ANN are investigated in a number of applications of Computational Mechanics with particular emphasis on Structural Mechanics, such as design optimization, reliability analysis and nonlinear structural dynamics. Adaptive neurocomputing schemes are proposed in order to predict the nonlinear structural behavior and drastically reduce the excessive computational cost. This approach has a broad field of applications making ANN a valuable tool for solving computational demanding problems in Mechanics.


Abdelaziz Timesli
email: timesli@univ-metz.fr

Abdelhadi Moufki
email: moufki@univ-metz.fr

Hamid Zahrouni
email: zahrouni@univ-metz.fr

Bouaaza Braikat
email: b.braikat@gmail.com

Hassan Lahmam
email: lahmam.hassane@gmail.com

Hassan Lahmam Numer

V Favier

Numerical model based on SPH method to simulate Friction Stir Welding

Keywords: FSW, Simulation, sans maillage, SPH

Introduction

Le soudage par friction et malaxage (FSW) a été inventé par l'Institut de Soudure Britannique TWI depuis les années 1990 pour les alliages d'aluminium [1]. Le principal avantage de cette technique est sa capacité à souder à l'état solide les alliages de métaux qui sont généralement difficiles à souder par des procédés de soudage conventionnels. L'assemblage de deux pièces par FSW consiste en la production de chaleur due principalement à l'épaulement et le malaxage du matériau grâce au pion. L'épaulement et le pion constituent les principaux éléments de l'outil de soudage. La chaleur est générée essentiellement par le frottement entre l'outil et les plaques à souder et par la dissipation plastique [START_REF] Colligan | Material Flow Behavior during Friction Stir Welding of Aluminium‖[END_REF]. Malgré les nombreuses applications industrielles du procédé FSW, son développement reste largement empirique et basé sur une grande connaissance expérimentale du procédé. Les travaux expérimentaux sont nombreux ; les simulations numériques aussi mais le traitement du malaxage reste une difficulté majeure pour les numériciens. La modélisation numérique du procédé de soudage par Friction et malaxage FSW (friction stir welding) a été étudiée par plusieurs auteurs considérant l'aspect thermique ou thermomécanique. Différentes formulations ont été proposées dans ces contributions concernant des formulations eulériennes, lagrangiennes ou ALE [START_REF] Guerdoux | A 3D numerical simulation of different phases of friction stir welding[END_REF][START_REF] Schmidt | Thermal modelling of friction stir welding‖[END_REF][START_REF] Lorrain | A contribution to a critical review of FSW numerical simulation‖[END_REF][START_REF] Bastier | Modelling of the residual state of friction stir welded plates[END_REF][START_REF] Feulvarch | 3D Modelling of Thermofluid Flow in Friction Stir Welding‖[END_REF]. Le choix d'une formulation spécifique dépend principalement du phénomène étudié. Le malaxage est difficile à réaliser à l'aide de la méthode des éléments finis (en lagrangien) puisque la zone proche de l'outil de soudage est le siège de grandes déformations. Quelques techniques sans maillage ont été proposées dans la littérature pour la simulation du soudage par friction et malaxage [START_REF] Alfaro | Numerical Simulation of Friction Stir Welding by Natural Element Methods‖[END_REF][START_REF] Tartakovsky | Modeling of Friction Stir Welding (FSW) Process with Smooth Particle Hydrodynamics (SPH)[END_REF]. Nous proposons dans ce travail d'utiliser la méthode sans maillage smoothed particle hydrodynamics (SPH) [START_REF] Monaghan | An introduction to SPH[END_REF][START_REF] Monaghan | Smoothed Particle Hydrodynamics‖[END_REF][START_REF] Liu | Smoothed Particle Hydrodynamics -A Meshfree Particle Method‖[END_REF][START_REF] Monaghan | On the problem of penetration in particle methods‖[END_REF] pour la simulation de ce procédé. Cette technique est basée sur un formalisme lagrangien, elle permet de prendre en compte les grandes déformations ainsi que l'effet de l'histoire du chargement. Ce modèle considère la zone de soudure comme un fluide non newtonien faiblement compressible où la viscosité dépend de la température. On se limite dans cette étude à une modélisation bidimensionnelle mais nous travaillons sur la modélisation tridimensionnelle avec une géométrie complexe de l'outil et une loi de comportement plus réaliste qui prend en compte l'effet de la vitesse de déformation ainsi que de l'écrouissage du matériau.

A notre connaissance, le seul travail qui propose la simulation du procédé FSW à l'aide de SPH est présenté dans la référence [START_REF] Tartakovsky | Modeling of Friction Stir Welding (FSW) Process with Smooth Particle Hydrodynamics (SPH)[END_REF] ; dans ce travail, le contact est géré par une force répulsive pour empêcher la pénétration entre les particules de l'outil et celles des plaques à souder. Dans notre modèle, le contact entre l'outil et les plaques à souder est géré par une force due au terme de pression ce qui permet d'éviter les complexités supplémentaires concernant l'estimation de la force de répulsion.

Les résultats du modèle proposé sont comparés avec ceux obtenus en utilisant un code industriel Fluent basé sur une formulation eulérienne de l'écoulement de fluide et une discrétisation par volumes finis du domaine étudié [START_REF]FLUENT 6.1 User's Guide[END_REF].

Equations de problème

On se propose dans cette section de présenter les différentes équations nécessaires à la modélisation du soudage par friction et malaxage. Ce procédé s'apparente aux procédés de fabrication par enlèvement de matière de type usinage. En effet, ce procédé utilise l'outil de soudage pour chauffer et malaxer le matériau. La zone proche de l'outil est évidemment le siège de grandes déformations. L'apport de la chaleur est dû à la fois au frottement de l'outil de soudage sur les plaques à souder et à la déformation importante de la matière dans le cordon de soudure. Du point de vu de la modélisation, nous avons un couplage thermo-mécanique comme dans la plupart des procédés de fabrication. Nous utilisons une loi de comportement simplifiée et qui ne dépend que de la température. Une loi plus complexe sera proposée dans des travaux futurs. Dans le présent travail, nous avons choisi de mettre l'accent sur l'implantation de la technique SPH pour la modélisation du malaxage du matériau. Les conditions de contact entre l'outil et la tôle sont considérées de type collant.

Nous présentons maintenant les différentes équations nécessaires à la modélisation du procédé de soudage FSW. Comme tous les procédés de fabrication, ce procédé peut être décrit par les lois de conservations suivantes:

Conservation de masse:

v dt d .      (1)
Conservation de quantité de mouvement:

    1 dt dv (2)
Conservation d'énergie:

v p q T k t T c       ) .(  (3) Source de chaleur:   v q v    :   (4)
Où v est le vecteur vitesse,  désigne la densité du matériau, σ est le tenseur des contraintes, c p désigne la chaleur spécifique, T est le champ de température, k désigne la conductivité thermique du matériau et  représente la fraction de l'énergie mécanique transformée en chaleur (ici 0.9   ) [15]. Le tenseur déviatorique  est lié au tenseur des contraintes  par la relation suivante:

     pI
(5) Où p est la pression hydrostatique et I le tenseur identité d'ordre 2. La relation liant le tenseur déviatorique des contraintes au taux de déformation peut s'écrire sous la forme suivante :

    ] . 3 2 [ I v v v T         (6)
où la viscosité des particules qui diminue avec l'augmentation de la température est considérée dans cette étude sous la forme suivante:

1 m ref melt ref TT A TT             (7)
A et m sont des paramètres du matériau. T ref est la température de référence et T melt est la température de fusion. La pression est obtenue en utilisant l'équation d'état suivante :

 2 c P  (8)
c étant la vitesse du son.

Principe de la méthode SPH

La méthode SPH est basée sur une approximation intégrale [START_REF] Monaghan | Smoothed Particle Hydrodynamics‖[END_REF]:

( ) ( ') ( ', ) ' f r f r W r r h dr    (9)
Où r est le vecteur position d'une particule, h définit la taille du domaine de lissage et W est la fonction noyau. Les formes approximatives de l'intégrale et du gradient de la fonction f par SPH s'écrivent sous la forme suivante:

  1 ( ) , Np a a b a b b b m f r f W r r h      (10)   1 ( ) ( ) , Np a a a b a b b b m f r f f W r r h        (11)
Où m b , ρ b et N p représentent respectivement la masse, la densité volumique de la particule b et le nombre des particules voisins. Plusieurs choix pour la fonction noyau sont possibles, en voici un exemple de fonction de type cubique spline [START_REF] Liu | Smoothed Particle Hydrodynamics -A Meshfree Particle Method‖[END_REF]:

2 3 3 3 3 1 si q 1 2 4 1 ( , ) (2
) si 1 q 2 4 0 si q 2 a b d q q W r r r h q                     (12) avec a b r r q h   et     3D en ) h 1/(
P P m v v r r m W W dt r r                            ( 14
) ) T - h(T - q ) , ( ) ( ) ( 4 ref a v ,               b a b a b b a b a b a b a b a a p r r h r r W T T k k k k m dt dT C   (15) 2 2 2 [( ).( )] ( ) ( ) ( ) b a b a b a b a b v a ab b a b a b a b m v v r r r r q W rr               (16)                     m ref melt réf a a T T T T A 1  (17) a a c p  2  (18)

Mouvement des particules

Les particules se déplacent en utilisant la variante XSPH [START_REF] Monaghan | On the problem of penetration in particle methods‖[END_REF] [START_REF] Monaghan | On the problem of penetration in particle methods‖[END_REF]:

() a b a a b ab b ab dr m v v v W dt       (19) Où 2 / ) ( ab b a     
et  est un paramètre donné dépendant du problème à résoudre. Cette méthode permet le déplacement de la particule ‗a' avec une vitesse moyenne par rapport aux particules voisines.

Conditions aux limites

Pour éviter la pénétration des particules représentant le matériau de la plaque à souder à travers l'outil de soudage, les particules de l'outil sont forcées de satisfaire les mêmes équations que les particules du matériau. Cependant, elles ne se déplacent pas en fonction de la variante XSPH, leurs vitesses sont imposées de la manière suivante:

, , , ( , , 0
) i i c t i c i c i c v r v r x x z z       (20)
Où  est la vitesse de rotation et t v est la vitesse de translation de l'outil. c x et c z sont les coordonnées de l'axe de rotation de l'outil. Quand une particule du matériau s'approche de l'outil, la densité des particules de l'outil augmente selon l'équation de continuité (Eq. 13) ; ce qui entraîne l'augmentation de pression suivant l'équation d'état (Eq. 18). Par conséquent, la force appliquée sur la particule du matériau augmente à cause du terme ) / ( 2  P dans l'équation de quantité de mouvement (Eq. 14). Finalement, si la distance entre une particule de l'outil et une particule de la plaque à souder devient inférieure à 2h, la densité, la pression et la force appliquée sur la particule augmentent et génèrent un mécanisme répulsif.

Schéma de temps utilisé

Le schéma de temps utilisé pour traiter les dérivées temporelles est un schéma explicit basé sur l'algorithme -velocity Verlet‖ [START_REF] Allen | Computer simulation of liquids[END_REF] comme suite:

1 n 1 1 n 1 a a 1 2 1 1 n 1 a 0.5 ( F ) ; 0.5 ( D ) 0.5 ; 0.5 ( E ) n n n n n n a a a a a a n n n n n n n a a a a a a a v v t F t D r r tV t F T T t E                           (21) avec a a dv F dt  , a a D dt d   , a a dr V dt  , a a E dt dT 

5

Résultats numériques

On considère une plaque formée d'un alliage d'aluminium de longueur L=50 mm et de largeur l=20 mm. Les paramètres de la loi (eq. 17 6). On peut observer qu'une erreur relative de moins de 5% est obtenue ; ce résultat montre la validité de l'algorithme proposé.

Conclusion

Dans ce travail, nous avons proposé un premier modèle basé sur la méthode SPH pour simuler le procédé de soudage FSW. L'avantage principal de cette technique concerne le malaxage du matériau autour de l'outil qui est difficilement réalisable avec d'autres méthodes numériques comme la méthode des éléments finis. Nous travaillons actuellement sur l'extension du modèle proposé à trois dimensions. Une loi de comportement plus réaliste dépendante à la fois de la température, de la vitesse de déformation et de l'écrouissage du matériau sera proposée. Par ailleurs, nous travaillons sur le couplage de modèles qui permettra la simulation du soudage pour des plaques de grandes dimensions. Un modèle local qui simulera le malaxage au voisinage de l'outil sera couplé à un modèle thermique qui prendra en compte l'ensemble de la structure et son environnement. Dans ce cadre, des géométries complexes de l'outil proches des géométries réelles des outils de soudage seront étudiées.

  et le schéma de temps utilisé L'utilisation des approximations SPH des fonctions et des dérivées spatiales permet d'écrire les équations du problème sous la forme suivante:

  ) sont définis comme suit : A = 200 Mpa.s -1 , m=1, T ref = 25°C, T melt = 502°C. L'outil est considéré rigide et fait en un acier avec une géométrie cylindrique (r p =2.5 mm rayon). Les propriétés physiques (conductivité matérielle k, chaleur spécifique Cp, et densité volumique ρ) utilisées pour la plaque d'aluminium sont les suivantes: physiques de la plaque à souder Dans le travail proposé, on a utilisé un coefficient de transfert thermique h lat =30 W/m2.K, une temperature ambiante T0=25°C, une vitesse de soudage de 2 mm/s et une vitesse de rotation de 20 rad/s.

Figure 1 :Figure 2 :Figure 3 : 2

 1232 Figure 1: géométrie et conditions aux limites de la configuration FSW traitée par SPH

Figure 4 :Figure 5 :Figure 6 :Figure 7 :Figure 8 :

 45678 Figure 4: géométrie et conditions aux limites de la configuration FSW traitée par Fluent


F Salvatore

T Mabrouki

H Hamdi

Presqu'île de Giens (Var) The use of numerical simulations to improve analytical chip formation model

Keywords: Johnson-Cook, analytical, cutting model, chip formation

In this paper, an analytical approach is proposed to model chip formation in the case of turning process. Numerical simulations of chip genesis are performed in order to fit efficiently the proposed analytical model. In particular, the variables like temperature and internal stress distribution are studied using finite element modeling. Numerical model setting is made with experimental and literature data using forces and the thickness of the primary shear zone.

Introduction

Turning and abrasion processes are widely used in different industries to cut different engineering parts. Usually the optimization of these processes is made by experimental methods often expensive and not able to be extrapolated to other machining configurations. To overcome these drawbacks some numerical simulations were carried out by many researchers but the major inconvenience of those methods are the long computing time, the high cost of numerical software etc. For all these reasons, in manufacturing industry, a highly interest in analytical methods like that of Merchant [START_REF] Merchant | Mechanics of the metal cutting process. I. Orthogonal cutting and a type 2 chip[END_REF] is usually researched because it is very practice and simple to use. Recently, Gilormini, Molinari and Oxley [START_REF] Gilormini | Modelisation thermomécanique de la formation du copeau en usinage à grande vitesse[END_REF][START_REF] Molinari | A new thermomecanical model of cutting applied to turning operations " Part 1 Theory -int[END_REF][START_REF] Oxley | Mechanics of machining an analytical approach to assessing machinability[END_REF] have the computation of the temperatures and shear zones thicknesses but the methodology to obtain those values is in general long and complicated for industrial community. In fact the differential equations presents in those thermo-mechanical models need time to be solved and also, in the same way, some parameters not easy to find in bibliography.

For all this reasons, the aim of the present paper is to present at the scientific communities a methodology to calculate the major numbers of variables of the material removing (cutting and ploughing forces, temperatures, spring back…). Here the equations can be directly used by engineers.

The proposed analytical modeling of the chip formation is performed using universal mechanical formulation and balance, and is improved by both numerical simulations and experimentations data (figure 1). In this case, first the experimental tests, the machines, materials employed and the measured data like forces and chip thickness are presented. Then numerical simulation using ABAQUS Explicit is developed in order to identify more accurately influences physical phenomena. Afterwards, an analytical modeling is proposed and numerical improvements are presented.

Experimental tests

In this section the experimental cutting tests are presented. In order to replicate orthogonal cutting conditions, the machining operation was done on disc with a diameter of 70 mm and a thickness (a p ) of 3 mm, which presents in the same time the cutting depth (Fig. 1).

The cutting tool is in a carbide grade (ref. TPKN 16 03 PP R SM30) with cutting edge radius R of 30 µm and the machined material is a steel alloy AISI 4140.

The values of the thickness of the material to cut (a p ) and the cutting radius of the tool (R) were chosen to minimize ploughing effects in order to have a real cutting process. Measures of forces in (x, y z) are done with dynamometer Kistler 9257 A with natural signal frequency 2000 Hz.

Details of the averages of forces, angle of the primary shear zone Φ and the thickness of the chip e c are presented in the table 1. The primary shear angle Φ is computed using equation 1

φ = Arc tan f e c       (1)

Numerical model

Numerical simulations based on Abaqus\Explicit were performed in order to both understand the physic of cutting and improve the analytic method. Johnson and Cook law is used for the workpiece material behavior [START_REF] Barge | Numerical modelling of orthogonal cutting: influence of numerical parameters[END_REF].

The workpiece is numerically modeled in 3 parts [START_REF] Mabrouki | A contribution to a qualitative understanding of thermo-mechanical effects during chip formation in hard turning[END_REF][START_REF] Salvatore | Analytical model of removal material in case of cutting and abrasion processes[END_REF]. In the part "final piece" only the behavior law of Johnson and Cook is defined and for the part "chip" and "transition" also the Johnson and Cook damage law is taken into consideration.

In order to find the good mesh dimensions for the numerical simulations it is decided to study, in a first time, the influence of the mesh size versus the primary shear zone thickness and the cutting force as showed in figure 2. In this specified study frictionless assumption is done for the chip-tool contact. In figure 2, it is possible to underline that for small mesh size values the primary shear zone thickness is small too. This is due to localization phenomena, which induce high temperature and afterwards material softening leading to lower forces and then lower thickness.

The computed primary shear zone thickness for mesh size from 5 to 10 µm is closed to experimental values (25 µm) given by Shaw [START_REF] Bourdet | La coupe des métaux[END_REF].

For all these reasons, the mesh size 10x10 µm is adopted in this study and seems in agreement with literature [START_REF] Barge | Numerical modelling of orthogonal cutting: influence of numerical parameters[END_REF].

The "surface to surface" interaction option and the penalty contact method are chosen in ABAQUS\EXPLICIT.

In particular, the friction adhesive coefficient between tool and workpiece is determined by Zemzemi method [START_REF] Zemzemi | Caracterisation de modèles de frottement aux interfaces piece-outils-copeau en usinage : application au cas de l'usinage des aciers et de l'inconel 718[END_REF].

Using mesh dimension 10x10 µm and the previous adhesive friction coefficient [START_REF] Zemzemi | Caracterisation de modèles de frottement aux interfaces piece-outils-copeau en usinage : application au cas de l'usinage des aciers et de l'inconel 718[END_REF], it is then possible to perform numerical simulation and compare results to experimental cutting forces primary shear zone angle Φ, chip thickness e c (table 2) It is important to underline that, in order to have a good description of the vertical forces (Y) during the cutting process, the ALE formulation is also adopted. In particular with ALE the cutting edge influence is described with more accuracy than lagrangian one.

Analytical model and numerical improvements

Orthogonal cutting represented by a 2D model is studied. The cutting tool removes a specific layer of work material. f is the theoretical uncut thickness to remove with the tool, R is the cutting edge radius considered zero in this paper.

The chip is assumed rigid except in zone 1 and 2 (figure 1) where all the deformations are concentrated. The thicknesses of these zones are d 1 and d 2 . In both zone, the assumption of thermoviscous-plastic-hardening material behavior is made. N 2 is the orthogonal force and T 2 is the tangential or friction force on the flank face of the tool. T 1 and N 1 are the tangential and normal forces applied on the primary shear zone.

It is supposed that in zone 2 Coulomb's law (equation 2) can be applied with friction coefficient µ ad .

2 2 2 ) ( N V T g ad µ = (2)
Different parameters, related to chip formation, can be expressed when balances of forces and moments applied to chip in static mechanical are written (like Merchant approach [START_REF] Merchant | Mechanics of the metal cutting process. I. Orthogonal cutting and a type 2 chip[END_REF]). Consequently, it is possible to write equations 3, 4, 5.

)) sin( ) cos( )( tan( ) sin( ) cos( ) sin( ) ( ) (cos( 1 2 γ γ µ γ µ γ - Φ - + Φ Φ + Φ = Tan T N (3) ) ( ) cos( )) sin( ) cos( ( 1 2 1 Φ + Φ - = Tan T N N γ γ µ (4) 2 1 1 N N l a l c = (5)
In equation 5, "a" depends on the distributions given to N 1 and N 2 . Φ is obtained by equation 6 like in Merchant's method where the minimum cutting power P C is considered.

{ } { }

)) cos( ) sin( (

2 2 γ γ N T V Min P Min C C + = (6)
It is important to underline that, if T 1 is constant like Merchant supposed, the values of Φ, computed with equation 6 is over evaluated and forces are under evaluated.

For this reason, it is decided to apply, in a first step, Pijpanen formulation [START_REF] Pijpanen | Theory of formation of metal chips[END_REF] to calibrate the model (equation 7). Where τ 10 is the simple shear stress characterizing the material and l 1 is the length of the primary shear zone.

It is important to note that equation 6 is the generalized Coulomb friction law where the first term is the "cohesive" contribution and the second one the "adhesive" term. The evolution of the coefficient c versus the cutting speed is presented in table 3. To compute analytically the coefficient c the followings steps are applied:

• Φ values based on experimental measurements are used in the analytical model as an input parameter.

• Analytical calculation is executed using different values of coefficient c. The computation is considered terminated when analytical Φ value is equal to experimental one.

Like showed in table 3 the Pijpanen coefficient computed with the analytical method is function of the cutting speed. Typically the bigger values of c are in correspondence of small cutting speed. This trend is similar to the friction coefficient trend between 2 surfaces.

In order to verify the values of c coefficient obtained analytically with equation 7, it is proposed a FEM validation (Fig. 3). In this last case T 10 is constant, T 1 and N 1 are calculated with numerical simulations and, than, it is possible to compute c numerically. The numerical and analytical results are similar. It is important to note that the method used to compute T 1 using c coefficient needs experimental or numerical approach for every configuration of parameters, in particular the values of this coefficient in equation 7. To overcome this drawback it is proposed in this paper also a different method using numerical simulations. In fact the ratio l c /l 1 remains constant when varying V C where lc is the contact length between chip and tool and l 1 is the length of the primary shear zone (table 4).

Afterwards, like showed in figure 6 and 7 it is possible to capture that normal stresses in zone 1 and 2 have similar values close to Tabor calculation (3σ y ). The region defined between zone 1 and 2 is totally plastically deformed and the triaxiality is important.

From those considerations it is possible to extract equation 9 where the normal stress n 1 is function of the tangential one τ 1 and le normal force N 2 in the second shear zone is a function of τ 1 .

1 2 1 3 τ = = n n . and c c l l n N 1 2 2 3 τ = = ( 9
)
Equation 8 is verified for different behavior law like showed in figure 4. For those reasons it is considered, in a first time, that it can be applied for different types of materials. It is important to underline that equations 8 and 9 replace equation 6 and 7, respectively. The total numbers of variables and equations are the same but now the model is predictive. Equation 8 and 9 in the proposed analytical model gives the calculation of cutting force and Φ according to experimental data (table 5), and this calculation don't need to minimize the cutting energy like Merchant. Many researchers think today that the employment of the equation of the minimum of the cutting energy to solve the problem is wrong; Molinari [START_REF] Molinari | The Merchant's model of orthogonal cutting revisited: a new insight into the modeling of chip formation[END_REF] suggest a correction of the Merchant formulation using the instability theory of the primary shear zone.

The primary shear force and stress is computed with Johnson and Cook law (equation 10 and 11).

                - - -                         +         + = • • + m f emp n n d T T T T d d CLn B A 0 0 1 1 1 1 1 10 1 1 3 3 γ γ γ τ ( 11
)
Now let study the calculation of the temperature in zone 1 to solve equation 11. For this reason the unitary volume of the zone where the deformation is concentrated is evaluated by equation 12, HH' is the displacement caused by the force T 1 , Q is the plastic energy converted into heat. The average value of temperature in zone 1 T emp1 is computed with equation 12 (definition of specific heat). In equation 12 δ 1 is the partition heat coefficient considered zero in this paper. In fact for modern cutting technology (Vc>100 m/min) the adiabatic hypothesis in zone 1 can be adopted [START_REF] Molinari | A new thermomecanical model of cutting applied to turning operations " Part 1 Theory -int[END_REF] like it is possible to see in figures 5 and 6. This result can be applied also to abrasion process where the local speed is significantly higher.

Conclusions

In this paper numerical simulation to improve a chip formation analytical model [START_REF] Salvatore | Analytical model of removal material in case of cutting and abrasion processes[END_REF] is presented. This model assumes that the chip can be considered rigid but all the deformations are concentrated in the shear zones. This model is based on numerical simulations and experimental results.

The main aim was to find simple model, which can be easily used in industry. In fact the output variables of the study (Forces, Temperatures, Φ, l c ) are computed with simple equations.

It is complete, taking into account the shear zones, the temperatures, and the contact length.

The present model can be qualified as predictive. In fact the calculation of the variables of the model are only function of the input parameters (working and workpiece parameters), and it is not required any experimental setting like chip thickness. In particular, the assumption that the tangential force in the primary shear zone is the sum of two terms (T 01 and cN 1 ) is rejected. Afterwards, the minimization of the cutting power proposed by Merchant is also discarded. Two new equations (8 and 9) are found using numerical simulations and are very important to solve the problem.

As a perspective the present analytical model can be exploited to predict residual stresses after machining and study a complex case of industrial cutting.

Fig. 1 -

 1 Fig. 1 -Description of the analytical model and the experimental and numerical methods used

Fig. 2 -

 2 Fig. 2 -Numerical simulations in case of AISI 4140 steel, f=0.15 mm, γ=0, µ=0, thickness of the primary shear zone vs cutting speed for different cutting speeds. For Vc=378 m/min the cutting forces are represented

Fig. 3 - 3 .

 33 Fig. 3 -The employed methodology to compute c in case of Vc=42 m/min

TABLE 4 .

 4 Numerical values of the ratio l c /l 1 in case of different cutting speed and coefficient m in Johnson and Cook law

Fig. 4 -

 4 Fig. 4 -Normal stress distribution in zone 1 and 2 in different cutting conditions and different sections

Fig. 5 :

 5 Fig. 5: Numerical simulations. Conduction and adiabatic case are represented.

Fig. 6 :

 6 Fig. 6: Numerical temperature distribution in zone 1 for conduction and adiabatic case.

TABLE 1 .

 1 Experimental values of X, Y forces and the chip thickness e c in the case of machining of AISI 4140, tool in carbide material (R=30µm), f=0.15mm, γ=0.

	Vc [m/min]	X [N]	Y [N]	e c [mm]	Φ Φ Φ Φ [Rad]
	42	1120	490	0.31	0.45
	126	1100	560	0.22	0.58
	378	1060		0.18	0.65

TABLE 2 .

 2 Numerical results: X and Y forces, Φ and the chip thickness e c in case of AISI 4140 steel, carbide tool f=0.15 mm, γ=0

		X [N]	Y [N]	e c [mm]	Φ [ Φ [ Φ [ Φ [Rad] ] ] ]
	Vc=42m/min	1150	520	0.32	0.46
	Vc=126m/min	1090	550	0.25	0.54
	Vc=378m/min	1070	590	0.20	0.63

TABLE 5 .

 5 Analytical section. Values of X force, Y force, angle Φ, temperature in zone 1 and the contact length tool-chip l c in case of 42CD4, tool in carbide material (R=30µm), f=0.15mm.

	Vc [m/min]	X [N]	Y [N]	Φ Φ Φ Φ [Rad]	l c x f	T emp1 [°C]
	42	1100	470	0.48	1.9	420
	126	1050	520	0.57	1.5	450
	378	1030		0.62	1.3	550
	5