List of Figures

2.9

Comparison between a linear mixture law (blue) and a quadratic mixture law (red)

2.10 The effect of using a variable time step on the mass conservation.. . . . . . . .

2.11

The effect of using P1mixture law instead of a P0 mixture law on the mass conservation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2.12

The effect of imposing the velocity on the inner of the dies . . . . . . . . . . .

2.13

Flat bunch simulation using a linear coefficient for mixture law. . . . . . . . . .

2.14

Representation of the mixture law for a flat bunch simulation using different coefficient weights: (a) linear -(b) inverse -(c) geometric -(d) logarithmic . . .

2.15

Flat bunch simulation using a inverse coefficient for mixture law. . . . . . . . .

2.16

Flat bunch simulation using a logarithmic coefficient for mixture law. . . . . . .

2.17

Top view of the zero iso-value of a flat bunch simulation in 3d showing the mixture law when using: a-linear coefficient, b-logarithmic coefficient. . . . . .

2.18

Instability of the levellerT solver. . . . . . . . . . . . . . . . . . . . . . . . .

2.19

Two different ways to compute the mesh size . . . . . . . . . . . . . . . . . . .

2.20

Comparison between two methods to compute the mesh size in CIMLib: (a) using the Edge min definition and (b) using the Altitude min definition. . . . .

2.21

A 3d flat bunch simulation after the ameliorations (logarithmic weigh, new mesh size definition ...). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2.22

Domain considered using a porous upper tool to show the Eulerian approach capacity to detect trapped air . . . . . . . . . . . . . . . . . . . . . . . . . . .

2.23

Initial mesh and mixture law for the problem in Figure 2.22. . . . . . . . . . .

2.24

The effect given by the imposition of the kinematic on the inner of the upper tool.

2.25

The evolution of the porous boundary condition in time. . . . . . . . . . . . . .

2.26

The evolution of the zero iso-value of the deformable body when using the porous boundary condition.

3.3

Illustration of a node containing bi-processors with three level of cache memory.

3.4

Top five super computers june 2013. [www.top500.org, ] . . . . . . . . . . . . .

3.5

Different types of speed up when running a parallel simulation : the common speed up (blue); the linear speed up (red); the super-linear speed up (green) . . .

3.6

Example of partitioning a global mesh in two sub-meshes . . . . . . . . . . . .

List of Figures

4.7

Geometry of the crankshaft immersed in the computational domain between the tools represented by the anisotropic mesh ( t = 0 s ). . . . . . . . . . . . . . . .

4.8

Mesh projected on the zero iso-value of the upper tool . . . . . . . . . . . . . .

4.9

Geometry of the crankshaft at t ≈ 0.9 s . . . . . . . . . . . . . . . . . . . . . . b) and (c) represent respectively the velocity components v x , v y and v z on the whole crankshaft geometry at t ≈ 0.9 s. . . . . . . . . . . . . . . . . . . . 4.11 (a) and (b) represent respectively the contours of the crankshaft between the tools at t ≈ 0.9 s and t ≈ 0.96 s. A growing penetration of the upper tool in the deformable piece is noticeable . . . . . . . . . . . . . . . . . . . . . . . . . .

List of Figures

4.12

The zero-isovalue of the connecting rod between the tools at t ≈ 1.6s . No penetration between the different bodies is noted. . . . . . . . . . . . . . . . . . . .

4.13

Initial geometry of the 10 parallelepiped immersed in the computational domain.

4.14

The uniform mesh of the domain presented on a section: Inside the virtual box, the mesh size is almost ten times smaller than the external one. . . . . . . . . .

4.15

The deformed geometries when using an isotropic mesh with a mixture law thickness equal to 0.9 mm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4.16

The initial mesh adapted anisotropically around the deformable parallelepiped (on the left) with a zoom (on the right).

. . . . . . . . . . . . . . . . . . . .

4.17

Introduction

This present work is devoted entirely to modeling multi-domain large deformation problems using a monolithic Eulerian approach in a massively parallel environment.

Generally, these type of problems are rather treated using Lagrangian or Arbitrary Lagrangian Eulerian (ALE) formulations (check chapter 1). Lagrangian approaches are spread widely and have almost perfect management for large deformation problems. But when it comes to multi-domain parallel deformation problems, they often show some shortcoming points. These imperfections are usually caused by the contact search algorithms affecting the scalability of the whole parallel simulation (see discussion in chapter 1). If one considers the contact problem between 10 deformable bodies treated on 16 cores (illustrated in Figure 1), the computation domain associated to each core is no longer a physical domain. It is rather a numerical domain containing parts of each deformable body. This means that every time contact search algorithm needs to identify a node in contact, communication between different cores is absolutely required. In other words, the communication cost is directly dependent of both the number of cores and interacting bodies.

In Cemef, we dispose of a library based on an Eulerian formulation fully functional in a parallel environment. In this context, our approach relies on a single meshed domain where different bodies are immersed and identified via the Level Set method. The different bodies are no longer individual domains but rather heterogeneity in the global computational one. We should note that one single mechanical set of equations is solved all over the domain. From the parallel point of view, it is expected to show great scalability even when using hundred cores since it depends of one unique mesh. In addition. this approach represents the advantage of not requiring any contact search algorithm. It is managed automatically via the Level Set method (see chapter 2). Now that all has been said, we propose using this Eulerian monolithic approach to try and prevent the previously mentioned problems which leads us to the motivations of this work.

Industrial Motivation and Aim of this study

Transvalor, house of Forge2 and Forge3, owns in partnership with Armines a C++ Library developed in Cemef called CIMLib developed by [Digonnet and Coupez, 2003], [START_REF] Digonnet | Adaptation de maillage en parallèle, application à la simulation de la mise en forme des matériaux[END_REF], [Mesri et al., 2009]. This library is based on an Eulerian description and is fully parallelized.

CIMLib is mostly used for computational fluid dynamics and it was never been used for large deformation problems (treated normally using Lagrangian formulations). To prevent difficulties confronted in Lagrangian approaches such as the ones mentioned earlier ( multi-domain large deformation problems, the complicated search contact algorithms, etc.), an Eulerian approach can offer some insight and eventually solutions for these problems.

Transvalor, interested in the capacity of this monolithic Eulerian approach considers utilizing this library for application similar to the ones carried-out in Forge. To do so, a proposition for the present thesis was conceived.

Though this monolithic Eulerian approach acquires impressing potential , it was never applied to these kind of problems before. Thus, it requires further developments and manipulations to be fully operational. So, the aim of this work can be summarized by the following headlines:

Study the feasibility of multi-domain deformation problem using an Eulerian approach in a massively parallel environment.

Identifying along this study several points such as where the approach excels and 2. The layout of this thesis 3 when it falls short.

Proposing improvements and new developments to enable applications to large deformation problems

In other words, the fundamental objective is to i) formulate a general idea if this monolithic approach is fit for large deformation problems, ii) to offer further needed developments and iii) to conclude if it shows promise to be industrialized on the long haul.

The layout of this thesis

This thesis is divided into five main chapters.

In chapter 1 a litterature survey for modeling large deformation problems is presented (Lagrangian, Eulerian and Arbitrary Lagrangian Eulerian formulations). It focuses mainly on previous works using Lagrangian formulations evoking both advantages and difficulties (contact search algorithms, mesh distorsions,..). The latter justify our choice leading to utilizing an Eulerian approach.

In chapter 2, the CIMLib library [Digonnet and Coupez, 2003] developed in a C++ parallel Eulerian framework is presented. First the Eulerian formulation, basis of this library, is detailed. Then, the different needed tools are introduced : from the level set method and the front convection, to the mixture laws and the anisotropic adapted mesh.

Last but not least, several tips and improvements are proposed such as introducing new mixture laws, new boudary conditions to treat air trapping and new definitions to handle some numerical instabilities. In addition, direct applications are offered to help readers understand how this existent approach can be adapted to modelling large deformation problems.

In chapter 3 the parallel enviroment in CIMLib is detailed such as the adopted parallel strategy, the solvers parallelization and the mesh partitionning. A full study dealing with the performance and scalabillity of different parts of the codes is presented. The overall performance was found satisfactory.

In chapter 4, advanced applications are displayed using all the previous developpment combined. Applications are divided into two halfs. The first models large deformations of complicated industrial pieces along with a comparison with Forge simulations. The second half is mainly dedicated to multi-domain deformation problem in a massively parallel enviroment. The approach performance is tested for upto 250 deformable bodies using couple hundreds of cores. We take advantage as well to point out some difficulties to be improved.

In chapter 5 the reasearch is pushed further by investigating the feasability of several concepts usefull to large deformation problems.

A first attempt to introduce friction notion using an Eulerian description is detailled. Direct applications and comparison with Forge simulation are presented. A new directional solver is developed as a first pursue to simulate contact resolution in an Eulerian formulation.

The different results and limitations of the model developments are discussed. This offers areas of further development along with ways to improve the proposed models for new research topics (not necessarly directly related to this present study)

Finally, general conclusions, perspectives and recommendations for futur works are provided.

Bibliography [Digonnet and Coupez, 2003] Digonnet, H. and Coupez, T. (2003). Object-oriented programming for fast and easy development of parallel applications in forming processes simulation. In 2nd MIT Conference on Computational Fluid and Solid Mechanics, page 1922.

[ [START_REF] Digonnet | Adaptation de maillage en parallèle, application à la simulation de la mise en forme des matériaux[END_REF] Digonnet, H., Silva, L., and [START_REF] Digonnet | Adaptation de maillage en parallèle, application à la simulation de la mise en forme des matériaux[END_REF]. Cimlib: a fully parallel application for numerical simulations based on components assembly. In MA-TERIALS PROCESSING AND DESIGN; Modeling, Simulation and Applications; NU-MIFORM'07; Proceedings of the 9th International Conference on Numerical Methods in Industrial Forming Processes, volume 908, pages 269-274. AIP Publishing.

[ [START_REF] Farhat | The two-level {FETI} method for static and dynamic plate problems part i: An optimal iterative solver for biharmonic systems[END_REF]] Farhat, C., Chen, P.-S., Mandel, J., and Roux, F. X. (1998). The two-level {FETI} method part ii: Extension to shell problems, parallel implementation and performance results. Computer Methods in Applied Mechanics and Engineering, 155(1-2):153 -179.

[ Mesri et al., 2009] Mesri, Y., [START_REF] Bruchon | Using a signed distance function for the simulation of metal forming processes: Formulation of the contact condition and mesh adaptation. from a lagrangian approach to an eulerian approach[END_REF]

Introduction

Modeling large deformation problems, such as forging and rolling processes, is one of the complicated problems in solid mechanics.

Their complexity requires a great care to details, specially the consideration of all interacting aspects in order to insure the precision of the simulation. We cite important aspects such as the material and geometrical non-linearities, the good manipulation of the contact and friction problems. For all these reasons, the use of numerous numerical techniques is a must. This subject has been widely discussed in the literature. A survey shows that the Lagrangian description is the most used dealing with such problems. The Arbitrary Lagrangian Eulerian (ALE) description comes in second. Finally, the Eulerian approach is rarely used for this type of problems (in solid mechanics). It is due to the lack of works presenting complete models dedicated to large solid deformation problems. Actually, the Eulerian approach is usually used in fluid mechanics.

In the present chapter, we will present • A quick description of the Lagrangian framework focusing on the contact algorithm.

• A brief overview of the Eulerian formulation. Note that this part will not be discussed in detail since chapter 2 will be entirely dedicated for it.

• A small review of the Arbitrary Lagrangian Eulerian description.

• A discussion summarizing the relevant points and justifying our choice selecting the Eulerian approach.

Lagrangian Description

The Lagrangian description of a deformable moving body consists on identifying the deformation of a material point (on an instant t) from its initial position (on instant t 0 ).

If a deformable body is occupying the domain Ω 0 in its reference configuration and Ω t in its current one (See Figure 1.1), a point deforming with the body is called:

• M 0 at t = 0 and occupies the position X = (X 1 , X 2, X 3 )

• M t at t and moves to the position x = (x 1 , x 2, x 3 ) Figure 1.1: Mesh deformation when using a Lagrangian approach.

The new position x is defined as a function of X and t:

φ : Ω 0 × T → Ω t (X, t) → φ(X, t) = x (1.1)
where φ is continuous and bijective. Its inverse Φ = φ -1 is continuous and bijective as well.

The function Φ exists and is bijective means that two distinct points initially remains distinct through the entire problem. Using Φ the position x can be computed as follows:

Φ : Ω t × T → Ω 0 (x, t) → Φ(x, t) = X (1.2)
In the Lagrangian description, all physical quantities (position, velocity, density, strain, stress, etc.) are expressed as functions of X and t:

The displacement or position vector is written as

OM t = x = φ(X, t) where φ(X,0)=X (1.3)
The velocity of M t takes the following form:

v(X, t) = ∂x ∂t = ∂φ(X, t) ∂t (1.4)
The same goes for the acceleration of M t :

γ(X, t) = ∂V ∂t = ∂ 2 φ(X, t) ∂t 2 (1.5)
To finish this section, one should at least present the constitutive equations to be solved in Lagrangian description:

dF dt = ∂F ∂t X = ∂F ∂t + ∂F ∂x ∂x ∂t = ∂F ∂t + (v.∇) F
This derivative will simply help transforming equations from Lagrangian description to the Eulerian description once needed in sec. 1.3.

Spatial discretization

To describe the deformation of a body, it is necessary to characterize the variations of lengths and angles. Since an analytical solution is quite impossible to find, a weak formulation is required. The transition to the weak formulation resort in introducing Sobolev spaces. Then the weak formulation is obtained by multiplying the equilibrium equation (1.6) by the test function and integrating the product over the entire area in the current configuration Ω t :

Find (v, p) ∈ V × P            ´Ωt ρ ∂v ∂t w + ´Ωt s(v) : ε(w) dV -´Ωt pdiv w dV + ´Γ τ.wdS = 0 ´Ωt qdiv vdV = 0 ∀(w, q) ∈ V 0 × P (1.7)
where s represents the deviatoric part of Cauchy stress tensor, τ the friction contribution and:

V = {v ∈ (H 1 (Ω)) 3 /(v -v r ).n = 0 over the contact surf ace} V 0 = {v ∈ (H 1 (Ω)) 3 /v.n = 0 over the contact surf ace} P = {p ∈ L 2 (Ω)}

Temporal discretization

When solving the system of equations (1.7), the choice of the time integration scheme varies mostly between an explicit scheme or an implicit scheme. By adopting an explicit scheme, the numerical system is the following at the increment i:

M i a i + F i int = F i ext (1.8)
where M, the mass matrix, is diagonal, F int is the internal vector forces and F ext is the external vector forces. Using the former scheme, the stiffness matrix is updated at the end of the incremental procedure. Furthermore the stiffness matrix obtained by adopting such scheme is diagonal, therefore the resolution of the numerical system is often given by the inversion of M. Moreover, no convergence loop is needed because the resolution is directly. This method is well suited for fast dynamics ( shocks, crash, explosion...) where the physical phenomenon is very violent and induces strong non linearities over a short time. However, this scheme needs many small increments for good accuracy and it is time consuming. Whereas, implicit schemes have interest in slow dynamics (such as vibration in structure problems...). The stiffness matrix is no longer diagonal.This scheme allows the use of bigger time steps. However in non-linear problems, Newton-Raphson (N-R) algorithm is often needed to insure equilibrium. This leads to one important drawback which is reconstructing the stiffness matrix at each N-R iteration.

To overcome the restriction on the time step when using an explicit scheme, some works as in [Bourel, 2006] propose to couple the multi-domain decomposition technique with a multi-scale time step. When adopting the latter technique, the time steps used are different on each sub-domain. Since the equilibrium equations are locally verified at each sub-domain at different times, the challenge is to develop a method capable of managing the exchange among sub-domains, even if they are never at equilibrium in the same time.

Figure 1.2: Illustration of the multi-scale time step adopted in explicit dynamic [Bourel, 2006].

Contact

In large deformation problems, contact algorithms have great impact on both solution quality and computational cost. In fact for such type of problems, the interacting geometries deform drastically leading to significant changes in the contact surface. The authors [Barboza, 2004] and [Yastrebov, 2011] recommend to update the contact surface at each time step. In [START_REF] Fourment | Numerical formulations and algorithms for solving contact problems in metal forming simulation[END_REF], the contact surface is updated at each Newton-Raphson iteration. Thus, the contact problem should be treated carefully. In other words, a poor treatment of the contact (e.g. non-smooth contact surface discretization) yields to an undesired numerical solution. In addition, as the number of nodes increases; the detection of nodes in contact affects the computational cost of the simulation specially when the contact surface should be updated frequently.

In what follows, we concentrate on frequently used numerical techniques dealing with contact problems. Theoretical aspects related to the existence and the uniqueness of the contact problem can be found in [Benson, 1992, Martins andOden, 1987] and [START_REF] Sofonea | Mathematical Models in Contact Mechanics[END_REF]. A large number of reviews on contact algorithms related to solid mechanics are presented in [START_REF] Aliabadi | Computational methods in contact mechanics[END_REF], Laursen, 2002, Hughes, 2012]. Additional publications devoted to investigate contact algorithms applied to metal forming can be found in [Chenot et al., 2002], [START_REF] Fourment | Numerical formulations and algorithms for solving contact problems in metal forming simulation[END_REF] and [START_REF] Hachani | A smoothing procedure based on quasi-c1 interpolation for 3d contact mechanics with applications to metal forming[END_REF].

As mentioned earlier, when using a Lagrangian formulation the boundaries of a body are represented by a set of boundary cells and nodes. In other words, the contact is highly dependent of the moving mesh which brings us to the main objective of this section.

In every contact problem, we can distinguish two major parts. The first one is devoted to contact detection in which a spatial search retrieving nodes in/or may possibly go in contact is established. Note that in most large deformation problems, the contact surface cannot be predicted and has to be detected at each time step. The second part is related to the resolution of contact problem in which several schemes can be used . It aims to eliminate the penetration between different bodies by applying repulsive forces.

In a large number of reviews, contact takes place between a rigid body and a deformable one. Note that in a Lagrangian framework, one should take into consideration self contact e.g. different parts of the same body are in contact, which is the case in a large deformation problems. Publications in this direction can be found in [Oliveira et al., 2008]. In addition, multiple contact problems can occur in many applications containing multi-body geometries. Work investigating the multiple contact problem can be found in [Chenot et al., 2002] , [Barboza, 2004] and [START_REF] Habraken | Contact between deformable solids: The fully coupled approach[END_REF]].

Contact search algorithms

Numerical treatment of contact problems confront many difficulties. The algorithms needed can consume a great part of the computational time and require in general a significant computer resources. These difficulties are related to three different topics: the geometry discretization, the contact detection and the spatial search.

The literature provides three different type of geometry discretization. The simplest method is called node to node discretization [START_REF] Francavilla | A note on numerical computation of elastic contact problems[END_REF]. It is recommended in case of small deformation problems. In fact, this technique consists on forming contact pair nodes. These pairs are formed by nodes belonging to a surface matched with nodes on the opponent surface. Once contact nodes are defined, it does not change during the time step. The construction of contact pairs and the assumption taken in consideration does not allow large deformation of the mesh. The second technique is called node to segment discretization. It demands more effort put into the implementation. Contact elements are formed by first searching for the closest nodes in two surfaces (the slave and the master). Then, the slave node forms with the segment containing the corresponding node of the master surface a contact element. The final technique is the segment to segment discretization. It requires even more effort in implementation since detection algorithms based on surface topologies are needed.

In large deformation problems, the node to segment discretization is frequently used since it is well adapted for a good representation of sliding problems. In other words, it is easier to represent sliding between a node and a segment then between two nodes or two segments [START_REF] Fourment | Master/slave algorithm for contact between deformable bodies and axial symmetries -application to 3D metal forging[END_REF].

Next, we will address the contact algorithms. As mentioned earlier, contact can occur at non-predectible surfaces. If an update is needed at each time step, which is always the case in our problem, then it is necessary to insure the efficiency of these algorithms. Generally, in every contact algorithm, two phases are distinguished : the global algorithm representing the spatial search and the local one representing the contact detection [Wriggers, 2002].

Among the spatial search algorithms, we cite the all-to-all algorithm, the bucket sort algorithm and the spherical sorting algorithm. 1.2.3.1.1 All-to-All algorithm is one of the first contact searching algorithm. It is easy to implement. However for problems with great numbers of nodes, this algorithm presents an essential drawback. The number of operations required to identify the contact pairs is proportional to the square of the number of boundary elements or nodes [Bourago and Kukudzhanov, 2005]. To reduce the computational time needed to identify the contact surface, several optimizations can be found in the literature. The main idea of these optimizations is to introduce a parameter d max -called maximal distance -. It is used to determine which nodes are relevant for identifying the contact surface e.g. all nodes lying at a distance greater than d max are not considered among the nodes involved in contact analysis. The choice of d max is discussed in [Yastrebov, 2011]. A poor choice of d max leads to a wrong contact surface definition as shown in Figure 1. 3 ( left part). According to [Yastrebov, 2011] work, the d max should be computed as:

d max = max i=Nm,j=N i n -1,k=N i n max i=1,j=1,k=j+1 r i j -r i k ; 2max Nc i=1 |∆r i | (1.9)
where N m is the number of the master segments, N i n is the number of master nodes forming the i -th master segment, N c is the total number of slave and master nodes and r i is the coordinate of the i -th node and ∆r i its displacement (Figure 1.3). In case of re-meshing, it is recommended to compute d max after each re-meshing step.

Figure 1.3: Improvement of the spatial search in the contact analysis [Yastrebov, 2011].

Even with the optimization proposed above, the number of nodes for the contact analysis remains significant. An additional improvement was brought by [Yastrebov, 2011]. The main idea is to determine a bounding box on each master and slave boundary surface. Then the nodes taken into consideration in the contact detection are reduced to the intersection of the two bounding boxes (see Figure 1.4). More details about the construction of these bounding boxes can be found in [Yastrebov, 2011].

Figure 1.4: Improvement of the spatial search suggested by [Yastrebov, 2011].

Bucket sort algorithm is an alternative technique used for spatial search.

It was proposed and applied for modeling self-contact surfaces in [START_REF] Benson | A single surface contact algorithm for the post-buckling analysis of shell structures[END_REF] and [START_REF] Belytschko | Contact-impact by the pinball algorithm with penalty and lagrangian methods[END_REF]. Several works relying on this algorithm can be found such as the works of [Heinstein et al., 2000] and [Paik et al., 2006]. Each of the previous cited works adapted and presented new variant of the algorithm responding to their own needs. In what follows, only the technique described in [Yastrebov, 2011] will be presented briefly.

Similar to the all-to-all algorithm, the preliminary phase is to define an appropriate d max . The author recommends to choose the biggest segment of the master surface. Then the bounding boxes are defined. Finally an internal bucket grid should be created, where the size w of each bucket should be chosen properly. For instance, w should be greater than √ 2d max when using linear elements. If the bucket size is smaller, contact spatial search can be mistreated. In this case, one should consider treating the contact detection on several bucket layers.

Spherical sorting algorithm

was used in [Papadopoulos and Taylor, 1993]. This algorithm is used as well in Forge ® for the contact spatial search. To identify areas where contact may occur, several spheres are created. The latter are defined for every boundary element. Its center is the barycenter of the element and its radius is the biggest distance connecting the barycenter with the nodes of the elements.

Once all the sphere are determined, the set of neighbor spheres is transformed into a new bigger sphere covering the whole set. This step can be repeated several times depending on the complexity of the geometry. The initial spheres are called the firstlevel-spheres, the bigger spheres englobing them are called the second-level-spheres and so on. Forge ® considers only three spheres levels ( see Figure 1.5). Once this hierarchy is established, the search algorithm consists on finding for every node (belonging to the slave surface) the closest third-level-sphere. The selected sphere englobes smaller ones. So another search is need to determine the closest level-two-sphere. The latter is formed by the surface elements of the master. This means that the set of the closest surface elements is known. The last step consists on defining the closest element belonging to this set.

Contact detection

After having found the areas where contact may occur, the local search is performed. Moreover, in contact detection algorithms, one should distinguish a contact between rigid and deformable body from contact between two deformable bodies. Several algorithms can be used to accomplish this detection such as the pinball algorithm [START_REF] Belytschko | Contact-impact by the pinball algorithm with penalty and lagrangian methods[END_REF] and the gap function [START_REF] Bathe | Stability of finite element mixed interpolations for contact problems[END_REF], [El-Abbasi and Bathe, 2001]. Other algorithms are inspired from computer graphics as the aircraft's shadow projection method. According to [Yastrebov, 2011], the latter will eliminate some drawbacks of the gap function.

Herein we present briefly a description of the gap function used in Forge ® . Recall that a master-slave formulation is put in place to describe the contact dynamics.

Contact between rigid and deformable bodies.

Let us denote by Ω r the rigid body and by Ω d the deformable body. Γ r and Γ d represent respectively the set of the boundary nodes of the rigid and deformable body. We begin by introducing the contact zone as an intersection between the boundary of two distinct bodies:

Γ c = Γ r ∩ Γ d (1.10)
Moreover, the contact should fulfill the no-penetration condition and the equilibrium condition on the interface. Those conditions are expressed by equations (1.11) and (1.12).

Ω r ∩ Ω d = 0

(1.11)

σ n1 + σ n2 = 0 with σ k = σ k n k ∀k ∈ {1, 2} (1.12)
The two previous conditions when combined provide Signorini law which determine the type of contact between the different domains. For the theoretical aspects, the readers can refer to [Wriggers, 2002]. The simplest way to express the no-penetration condition is by introducing the gap function :

g(x 2 , t) = min x 1 ∈Γ 1 ||x 2 -x 1 ||sign(n 1 .(x 1 -x 2 )) 0 ∀x 2 ∈ Γ 2 (1.13)
where n 1 is the unit outward normal at the surface Γ 1 (see Figure 1.6).

Contact between deformable bodies.

In many applications, specially in a multi-domain cases, multiple contact between deformable bodies occur. A broad range of algorithms allowing to fulfill the task can be found in the literature [START_REF] Flickinger | What's wrong with collision detection in multibody dynamics simulation?[END_REF], [START_REF] Hahn | An explicit multi-body contact algorithm[END_REF] and [Iwai et al., 1999]. Similar to the contact between rigid and deformable body, a contact graph is created. A review of different algorithms and implementations are found in [Erleben, 2004]. Note that in multi-domain algorithms, nodes do not always represent a physical entity. For instance, it can be a virtual entity used to trigger the contact. These nodes will form fictitious elements used for contact analysis. The latter concept is applied in Forge ® .

Herein a brief description is presented. For more details, the readers can refer to [Barboza, 2004]. The following algorithm can be considered as a generalization of the contact analysis between a rigid and a deformable body. The global search relies on the spherical sorting algorithm presented above.

As for the local detection algorithm, let Ω k and Ω k denote two different deformable bodies, where Γ k and Γ k represent respectively their boundaries. For each node x k of the slave surface, the contact element is constructed as follows:

π : ∂Ω k h -→ ∂Ω k h x k -→ x k = π(x k )
where π is the orthogonal projection. The projection determines the location of x k on the master surface e.g. x k belonging to which element. The position of x k is computed using the element's barycenter and nodes {A 1 , A 2 , A 3 } (see Figure 1.7). The distance between the two bodies is given by : Before proceeding into the parallel strategies, we would like to mention that the literature proposes alternatives methods to make solving large contact problems a possibility. We cite for instance the domain decomposition method presented in [START_REF] Sassi | Generalization of lions' nonoverlapping domain decomposition method for contact problems[END_REF]Garaud, 2009]. Many variants can be found also as in [Farhat andMandel, 1998, Farhat et al., 1998]. In general the computational domain is divided into non-overlapping sub-domains. A part of the solution is evaluated on a sub-domain then a gluing phase is needed to enforce field continuity [Roux et al., 2005]. These decomposition algorithms make solving large problem a possibility without resorting to parallel algorithm.

g n (x k ) = |x k x k | (1.14) n(x k ) =      - 1 g n (x k ) ---→ x k x k if g n (x k ) = 0 n F k if g n (x k ) =

Parallel strategies in contact detection

Nowadays, most problems are oriented towards industrial applications which sometimes require complex geometries. The latter will involve a large number of contact nodes.

Even with all simplifications introduced to replace a full scale contact detection, the execution time of such algorithms remains important. To optimize such algorithms, developers used parallelism as a mean to decrease the required computed time.

In what follows, the parallel strategy used in Forge ® is presented briefly. The focus will be on the generation of fictitious elements in multi-domain cases. For sake of clarity, two deformable bodies Ω 1 and Ω 2 are taken into account next. To begin with, we should mention that the parallelism of contact search algorithm relies on a S.P.M.D strategy ( see chapter 3 for more details ). Every body is partitioned into n p sub-domains :

1.2. Lagrangian Description 21 Ω 1 = np p=1 Ω 1 p , Γ 1 pl = Ω 1 p ∩ Ω 1 l (1.16) Ω 2 = np p=1 Ω 2 p , Γ 2 pl = Ω 2 p ∩ Ω 2 l (1.17)
where n p is the number of cores, Γ 1 pl and Γ 2 pl are respectively the interfaces between sub-domains of Ω 1 and Ω 2 . After partitioning Ω 1 and Ω 2 , the sub-domains are seen as the sub-meshes formed by each core :

Ω = Ω 1 ∪ Ω 2 = np p=1 Ω p , Ω p = Ω 1 p ∪ Ω 2 p (1.18) Γ pl = Γ 1 pl ∪ Γ 2 pl = Ω 1 p ∪ Ω 2 p ∩ Ω 1 l ∪ Ω 2 l (1.19)
After redefining the computational domain Ω as in equation (1.18), the simulation is launched on each core. The contact surface Γ c is duplicated on each core eliminating any needed communications. Using this strategy, the contact analysis is carried out between each sub-domain Ω 1 p and the global surface ∂Ω 2 . As for the parallelization of the fictitious element algorithm generation, the mesh of both the domains and the sub-domains are identified via the following numbering:

M(Ω) = {1, ..., nbnoe} where Ω = Ω 1 ∪ Ω 2 (1.20) M(Ω p ) = {1, ..., nbnoe(p)} where Ω p = Ω 1 p ∪ Ω 2 p (1.21)
in which nbnoe is the total number of nodes and nbnoe(p) the number of nodes assigned to a core p.

Recalling paragraph 1.2.3.2.2, the contact analysis requires to compute the distance between a node j ∈ M (Ω p ) and a surface A = {A 1 , A 2 , A 3 } ∈ Ω ( using a node to segment technique ). The set of node (j, A 1 , A 2 , A 3 ) form a contact element which does not belong necessary to Ω p . As a solution, [Barboza, 2004] proposes to add each node A i={1,2,3} to the core containing the node j. In fact, the interface Γ pl should be updated frequently by adding all nodes in contact (belonging to Ω l ) to the core p (e.g. to the set on nodes associated to the core p). This technique is not optimal in term of memory consumption. This strategy is tested in a triaxial compression simulation. Table 1.1 shows the efficiency of this strategy adopted by Forge ® . We can notice that when using 4 cores, the efficiency decrease to 0.65.

Contact resolution

Numerous methods can be found in the literature dedicated to the contact resolution problem. A survey on this subject confirms that the most used methods are the Lagrange multipliers, the penalty method and the augmented Lagrangian method.

Though other methods can be used as well, such as the invariant point algorithm [Chabrand et al., 1998, Duvaut andLions, 1972], the below descriptions will be only restricted to the previously cited methods.

Lagrange Multiplier Method

This method consists on introducing new terms representing contact forces [Kikuchi andOden, 1988] [Wriggers and[START_REF] Wriggers | [END_REF]. These terms are called Lagrange multipliers. Thus the set of equations (1.7) is replaced by (1.22) introducing λ.h the Lagrange multiplier.

       ´Ω s(v) : ε(w)dV -´Ω pdiv wdV -´Γc λh(w)dS = 0 ´Ω qdiv vdV = 0 -´Γc λ * h(v)dS ≤ 0 (1.22)
1.2. Lagrangian Description where (1.22.c) represents the non-penetration criterion (to be respected) and h(v) = (v -v r ).ng ∆t . This technique presents both drawbacks and advantages. Though it increases the number of unknowns (i.e. the computational time), it respects exactly the contact condition (or what we called non-penetration condition). In addition, it has no need for further parameter manipulation -normally confronted in penalty method-(see section below)

Penalty Technique

This is one of the firsts techniques treating contact problems. It relies on the fact that numerical solutions are the approximation of the actual solution. So the contact interface is identified with a certain tolerance ε prec . Respecting the contact criterion depends of this precision ε prec as well.

This method consists on introducing penalty coefficient ρ to approximate the contact law. The set of equations (1.7) is now transformed into:

           ´Ω s(v) : ε(w)dV -´Ω pdiv wdV -ρ x k ∈Γc [h (v (x k ))] + S (x k ) (n (x k ) .v) = 0 ´Ω qdiv vdV = 0 (1.23)
where ρ is the penalty coefficient, S(x k ) is the surface associated to the node k and verifies S(x k ) = ˆ∂Γ N k dS. Thus, σ n the normal contact stress is approximated by

σ n (x k ) = -ρh(v h (x k )).
Note that the precision ε prec ≈ O 1 ρ . So, to respect exactly the non-penetration conditions ε prec → 0 and ρ → +∞ must be verified.

Though the penalty method does not increase the unknowns number, it cannot respect exactly the non-penetration conditions unless ρ → +∞ . This generates further complications regarding the badly conditioned matrices (affecting the algorithm convergence).

Augmented Lagrangian Method

As presented previously for resolution methods, a choice between two alternatives has to be made. Using the "Lagrangian multiplier method" insures the respect of the contact criterion but increases the unknowns. Using the "penalty method" offers exactly the opposite. Which brings us to a new compromise: The "augmented Lagrangian method".

As proposed in [START_REF] Fortin | Augmented lagrangian methods: Applications to the numerical solution of boundary-value problems[END_REF], a hybrid method called augmented Lagrangian is introduced. The concept is to solve the problem using the Lagrange multiplier and adding a penalty contribution as well. This method combines the advantages of the latter ones. The Lagrange multiplier contribution helps verifying the non-penetration condition and the penalty contribution eases the algorithm convergence. For more details, the readers are invited to check [Kuss, 2008] [Chamoret, 2002].

Eulerian Description

Although the physical fields are the same, the Eulerian and Lagrangian descriptions are quite different. Rather than following each material point, the evolution of fields properties are recorded at every point in space as time varies. In other words, two reference frames are used to evaluate the physical fields. This approach is very common in fluid mechanics. Equation (1.1) represents the transformation function φ in the Lagrangian description. As mentioned earlier, an inverse transformation function Φ = φ -1 can be defined. The initial position of a material point M t can now be defined as follows:

X = Φ(x, t) (1.24)
where x is the spatial variable describing the current point position.

By using (1.24), the deformation is no longer expressed in terms of the Lagrangian variable (X, t). Instead it is expressed in terms of the Eulerien spatial variables (x, t). Let us take the velocity for instance:

v(x, t) = dx dt = ∂φ(X, t) ∂t = ∂φ(Φ(x, t)) ∂t (1.25)
For this reason, when applying a derivative in the Eulerian description two terms show up according to the remark on page 11:

dF dt = ∂F ∂t + (v.∇)F (1.26)            ρ ∂v ∂t + v.∇v = ∇.σ + ρb = ρg ∇.v = 0 σ = behavior law (1.27)
In fact, by adopting an Eulerian description the mesh is considered fixed in space. The fixed mesh does not represent a physical characteristic. Therefore in a multi-domain problem, tools representing each sub-domain are required; which brings us to the next point: Interface modeling.

Eulerian Description

Interface modeling

A literature survey shows the existence of several methods for dealing with the interface modeling in Eulerian formulation. Mainly, these methods can be divided into two categories: the interface tracking based on a Lagrangian approach and the interface capturing based on an Eulerian approach.

For interface tracking, two methods stand out the most: the volume tracking and the front tracking. the first was introduced by [START_REF] Harlow | Numerical Calculation of Time-Dependent Viscous Incompressible Flow of Fluid with Free Surface[END_REF]. It consists on adding particles to represent the volume. The efficiency of this method depends highly on the particles number and their distribution. As a result, it can be expensive in term of memory and computing time. The front tacking method was presented as an improvement of the latter in [Daly and Pracht, 1968]. Instead of distributing particles all over the volume, they are only distributed over the interface. This change plays a key-role in decreasing the memory cost. For these methods, the main drawbacks are the need to reposition particles at each time step and the lack of information regarding the topology changes.

As for interface capturing, two methods exist as well : The volume of Fluid [START_REF] Hirt | Volume of fluid (vof) method for the dynamics of free boundaries[END_REF] and the Level set method [START_REF] Osher | Fronts propagating with curvature dependent speed: algorithms based on hamilton-jacobi formulations[END_REF]. Both methods rely on a scalar parameter to describe the domains interface. The volume of fluid method uses a fraction function C determining if a mesh element is fully or partially filled. If the interface is beyond an element -i.e. element fully filled-, the fraction function is equal to one. If the interface does not reach the element, C takes a zero value. Otherwise, C is equal to the filling of the element volume ratio. Whereas, the level set method uses a signed distance function to describe the interface of a domain. The interface is captured for a time t as the zero of the level set function φ ( i.e., Γ(t) = {x|φ(x, t) = 0} . Then, φ is chosen to be positive inside the domain and negative outside. The interface motion in both methods is given by convecting the corresponding scalar parameter with the velocity v. We describe the level set method in detail in the next chapter of this manuscript.

Contact formulation

After reviewing the literature, rare are the works dedicated to contact problems in Eulerian framework. To our knowledge, the only existing work is presented in [START_REF] Bruchon | Using a signed distance function for the simulation of metal forming processes: Formulation of the contact condition and mesh adaptation. from a lagrangian approach to an eulerian approach[END_REF]. Herein, we present a brief description.

Γ c = {x ∈ Ω; |φ d (x, t)| + |φ s (x, t)| = 0}
(1.28) Thus, the no-penetration condition is expressed by:

φ d (x, t) + φ s (x, t) 0 ∀(x, t) ∈ Γ c × R + (1.29)
If v and v s denote respectively the velocity of the deformable body and the tool, the authors propose linearizing the previous inequality in the following manner:

g(φ d , φ s ) = φ d ∆t + φ s ∆t + (v -v s ) .n s 0 ∀(x, t) ∈ Γ c × R + (1.30)
where g is the gap function defined using the level set functions and the normal vector n s . Proceeding by analogy with a Lagrangian formulation, the authors wrote the contact term in the weak formulation using a penalty method.

ˆΓc r [g(φ d , φ s )] + w.n s dΓ (1.31)
where r is the penalty parameter.

However, this surface integral is not computable within an Eulerian framework. In fact, the contact surface is not defined explicitly (see equation (1.28) ). Thus, surface integral (1.31) has to be converted into a volume integral.

At first, the authors of [START_REF] Bruchon | Using a signed distance function for the simulation of metal forming processes: Formulation of the contact condition and mesh adaptation. from a lagrangian approach to an eulerian approach[END_REF] introduce the function ξ : R → R such as :

ξ(y) =      1 2 (1 + cos(πy)) if |y| < 1 0 if |y| 1 (1.32)
Using the function defined in equation (1.32), the authors approximate the Dirac mass by the right term in equation (1.33) when ε tends towards zero:

δ Γc ≈ 1 ε ξ |φ d | + |φ s | ε |∇ |φ d | + ∇ |φ s || (1.33)
Since the contact term is written in the vicinity of the zero isovalue of the level set function, the previous equation is simplified using the property of distance function.

δ Γc ≈ 2 ε ξ |φ d | + |φ s | ε (1.34)
Finally, the contact term added to the Eulerian weak formulation is given by: 1.

Arbitrary Lagrangian Eulerian Description

This section presents the theoretical aspects of Arbitrary Lagrangian Eulerian formulation (ALE) and the different resolution techniques used in solid mechanics. This approach was designed to address problems confronted by adopting a Lagrangian approach while maintaining its advantages. When the ALE approach was introduced, each subdomain was either purely Lagrangian or purely Eulerian. Then the approach was extended to the first arbitrary relative movements between the mesh and material. Nowadays, ALE formulations are widely used treating fluid structure interactions problems [Richter andWick, 2010, Kucharik et al., 2010]. Nevertheless, some works use such approach in solids mechanics [START_REF] Wang | Formulation and survey of {ALE} method in nonlinear solid mechanics[END_REF]. A review of material forming processes using the ALE formulation can be found in [Chenot and Bellet, 1995].

From the mathematical point of view, an area located in the space can be identified using three different frames : a material, a spatial and a reference one. For each frame, a system of coordinates is associated. Position vectors are respectively denoted X, x and χ. The reference frame can be chosen properly to facilitate solving the numerical problem. In a Lagrangian approach, the reference frame coincides with the material one. While in an Eulerian approach, it coincides with the spatial frame. For the ALE approach the reference frame is chosen arbitrary and the material moves with a velocity different from the mesh velocity. The separation of these two velocities requires a new bijective function ψ joining the material frame with the referential one :

χ = ψ(X, t)
(1.36)

The vector position χ can be expressed using the functions φ and φ * crossing from the Lagrangian to the Eulerian frame then to the ALE frame (Figure 1.9) :

χ = φ * • φ(X, t)
(1.37) Figure 1.9: Mesh deformation when using an Arbitrary Lagrangian Eulerian approach.

To transform equations from the Lagrangian to the ALE description the following relation is used for every characteristic

F (χ) = F (χ(x, t)) = F (χ(x(X, t), t)) : dF dt = ∂F ∂t X = ∂F ∂t + ∂F ∂χ ∂χ ∂t = ∂F ∂t + ∂F ∂x ∂x ∂χ ∂χ ∂t = ∂F ∂t + (c.∇) F (1.38)
in which c = v -v is the convective velocity. v is the velocity of the material and v is the velocity of the grid. It represents the relative velocity between the material frame and the grid (the spatial frame).

Remark:

• In a Lagrangian approach, since the material is represented by the mesh, the mesh velocity will be equal to the material velocity : v = v =⇒ c = 0.

• In an Eulerian approach, since the mesh is fixed, the mesh velocity is given by v = 0. Equation (1.38) will be expressed as follows:

∂ν ∂t = ∂ν ∂t X -v∇ X .ν (1.39)
Finally, the conservation laws are transformed using equation (1.38). This formulation is also called quasi-Eulerian formulation, since these equations are convected with the relative velocity c.

ρ   ∂v ∂t χ + (c.∇) v   = ∇.σ + ρb (1.41) ρ   ∂e ∂t χ + c.∇e   = σ : D + ρr + ∇.q (1.42)
where ρ is the density, v is the velocity, σ is the Cauchy stress tensor, D is the strain rate tensor, b is the body force, e is the internal energy.

Coupled Resolution

Solving the previous equations using a finite element discretization requires the introduction of a variational formulation. The principle of the latter formulation is similar to the one explained in section (sec. 1.2.1).

ˆρ ∂v i ∂t χ w i dΩ + ˆρc j ∂v i ∂x j w i dΩ + ˆ∂w i ∂x j σ ij dΩ = ˆρb i w i dΩ + ˆσij n j w i dΓ (1.43)
In a coupled resolution, a convective term appears in equation (1.43). This requires stabilization techniques. The literature proposes several techniques such as SUPG -SCPG are the most used. To solve this type of problem, two methods are mentioned in the literature. The first is to consider that the mesh movement is known a priori [Liu et al., 1986]. This technique makes the mesh velocity independent from the material velocity which remains the only unknown of the problem. However, in most cases the definition of the mesh motion is not obvious. Therefore another method to solve this problem is needed where the mesh motion is unknown. Under the latter assumption, the introduction of new equations describing the mesh motion is required [START_REF] Schreurs | Simulation of forming processes, using the arbitrary eulerian-lagrangian formulation[END_REF]. For a more elaborate bibliography on this matters the readers are invited to check [Boman, 2010, Philippe, 2009].

Split operator

In solids mechanics, the decoupled approach is more used. The following explained method is extracted from the work of [Benson, 1997, Benson, 1998a]. This method eliminates the convective term of the previous formulation (1.43). It consists on solving the equations in two consecutive steps. The system of equations (1.44) to (1.46) presented previously is written in the following conservative form:

∂ρ ∂t + ∇.(ρv) = 0 (1.44) ∂ρv ∂t + ∇.(ρv ⊗ v) = ∇.σ + ρb (1.45) ∂ρe ∂t + ∇.(ρe.v) = σ : D + ρb.v (1.46)
where ρ is the density, v is the velocity, σ is the Cauchy stress tensor, D is the strain rate tensor, b is the body force and e is the internal energy.

Denoting φ = (ρ, ρv, ρe), F = (ρv, ρv ⊗ v, e.v) the flux function and S = (0, ∇.σ + ρb, σ : D) the source term, equations (1.44-1.46) are rewritten as :

∂φ ∂t + ∇.F = S (1.47)
The splitting Operator divides equation (1.47) into two equations:

∂φ ∂t = S (1.48) ∂φ ∂t + ∇.F = 0 (1.49)
This method consists in solving the equations in two alternating steps. The first solves the equation using purely a Lagrangian approach where all the above definitions remains valid. The second step, called Eulerian, transfers all necessary fields from the deformed Lagrangian mesh to the fixed Eulerian mesh. Figure 1.10 depicts this two step with the mesh evolution. The decoupled method offers several advantages over the coupled method. First, this method is faster and more robust. In addition, the time step is independent of the convective velocity. Finally, the decoupled approach is very easy to implement in a 1.5. Discussion and Conclusion 31 Lagrangian code. It only requires the addition of a module for determining the speed of the mesh and a transport module. It is at this point that each increment is completed. Several algorithms are available to complete the task. The most intuitive transfer method is to perform an interpolation between the two meshes. This type of techniques is also used after a complete re-meshing for which the number of nodes and the number of elements changes [Philippe, 2009]. An alternative transport technique is represented by solving (1.49).

As mentioned earlier, the ALE approach helps overcoming the problem of mesh distortion encountered in a Lagrangian approach. The mesh motion is given by re-positioning algorithms. According to the nodes nature, these algorithms can be very different. For instance, a node on the mesh surface is generally constrained to remain thereon, while an internal node can be moved with more freedom. These algorithms are divided into two categories: i) mono-material algorithms known as Simple ALE [Benson, 1989] and ii) multi-material [START_REF] Benson | Stable time step estimation for multi-material eulerian hydrocodes[END_REF]. Multi-material algorithms require sophisticated techniques to determine the boundaries between the involved bodies and then manage the interactions between them [Benson, 1997, Benson andOkazawa, 2004].

Though this section does not cover all aspects in an ALE approach, we should mention that it exhibits the same drawbacks as the Lagrangian approach regarding the efficiency of contact detection algorithms in a parallel environment.

Discussion and Conclusion

For large deformation problems, the most widespread and most developed methods rely on Lagrangian formulations [Benson, 1997]. Though the mesh suffers a distortion, a re-meshing algorithm is enough to fix this problem and insure the solution accuracy. Some works choose to use an Arbitrary Lagrangian Eulerian to prevent the mesh distortion. With the increasing number of industrial applications including multiple interacting components and complex geometries, a parallel environment is a must. Lagrangian approaches present drawbacks mainly for multi-domain problems in a parallel framework. In particular, the contact remains the topical issue where the spatial search algorithm affects highly the scalability of the whole computation. In Forge ® , this defect is presented in Table 1.1. Almost half the efficiency is lost when using four cores. In addition, the authors of [Paik et al., 2006] present a study on this particular topic using the bucket sort algorithm. Figure 1.11 confirms that the loss of efficiency is mainly due to the contact calculation. A worth mentioning work dealing with the parallel contact detection is presented in [Li et al., 2001]. The authors introduce a mesh-free method presenting several improvement such as a simple implementation and a better scalability.

Facing all these needs, the Eulerian approach seems to be a good candidate. This formulation uses a single fixed mesh even in presence of different interacting bodies. It presents the advantage of automatically detection contact. All these advantages does not mean that the formulation will not present limitations. For instance, it opens doors on how to treat friction and contact in such problems. Due to the small number of works treating the large deformation problem in solid mechanics, it shows a lack of development and unclear insight on what will or will not work. This is where our work comes in handy. It studies the feasibility of the Eulerian approach to model solid mechanics applications.

Résumé en Français

Dans ce chapitre, une enquête citant différentes approches disponibles dans la littérature dédiée pour la modélisation des problèmes de grandes déformations multi-domaines a été présentée.

Pour des problèmes de grandes déformations, les approches les plus utilisées et les plus développées se reposent sur des formulations Lagrangiennes où le maillage peut subir des distorsions. Quoique ces distorsions puissent être résolues en utilisant des algorithmes de remaillage, plusieurs travaux se servent des approches Lagrangiennes Arbitrairement Eulériennes comme une alternative.

Avec la demande croissante de modéliser des applications industrielles incluant de multiples géométries (complexes) interagissant au cours du procédé, le travail dans un environnement parallèle est prouvé nécessaire. Les approches lagrangiennes confrontent principalement des difficultés pour les problèmes multi-domaines dans un cadre parallèle. Plus particulièrement, le problème de contact reste le sujet d'actualité où l'algorithme de recherche spatiale affecte fortement l'efficacité du calcul entier.

Face à ces besoins, l'approche Eulérienne semble être un bon candidat. Cette formulation utilise un maillage unique même en présence de différents corps interagissant durant la simulation. La détection de contact est gérée automatiquement et une très bonne scalabilité est attendue. Vis-à-vis à ces avantages, l'utilisation de cette approche mène à de nouvelles questions préalablement non traitées telles que la modélisation du frottement et du contact.

En raison du petit nombre de travaux dédiés à la modélisation des problèmes de grandes déformations de la mécanique solide, les idées sont peu claires sur la direction dont nous devons mener les développements.

Ceci est où notre travail arrive à point nommé. Il étudie la faisabilité de l'approche Eulérienne en proposant les atouts et les développements additionnels nécessaires pour modéliser des problèmes de grandes déformations multi-domaines.

[ Boman, 2010] [Chenot and Bellet, 1995] Chenot, J. and Bellet, M. (1995). The ale method for the numerical simulation of material forming processes. Proc. NUMIFORM, eds, S.-F. Shen and P.R. Dawson, Balkema A.A., Rotterdam, pages 39-48. [Chenot et al., 2002] Chenot, J.-L., Fourment, L., and Mocellin, K. (2002). Numerical treatment of contact and friction in {FE} simulation of forming processes. Journal of Materials Processing Technology, 125-126(0):45 -52.

[Daly and Pracht, 1968] Daly, B. J. and Pracht, W. E. ( 1968). Numerical study of density current surges. Physics of Fluids , 11(1):15-30. [Erleben, 2004] Erleben, K. ( 2004). Contact graphs in multibody dynamics simulation.

[ Heinstein et al., 2000] Heinstein, M. W., Mello, F. J., Attaway, S. W., and Laursen, T. A. (2000). Contact-impact modeling in explicit transient dynamics. Computer Methods in Applied Mechanics and Engineering,:621 -640.

[ Hirt and[START_REF] Hirt | Volume of fluid (vof) method for the dynamics of free boundaries[END_REF][START_REF] Hirt | Volume of fluid (vof) method for the dynamics of free boundaries[END_REF]. Volume of fluid (vof) method for the dynamics of free boundaries. Journal of Computational Physics, 39(1):201 -225.

[Hughes, 2012] Hughes, T. J. ( 2012). The finite element method: linear static and dynamic finite element analysis. Courier Dover Publications.

[ Iwai et al., 1999] Iwai, T., Hong, C.-W., and Greil, P. (1999). Fast particle pair detection algorithms for particle simulations. International Journal of Modern Physics C, 10(05):823-837.

[ Kikuchi and Oden, 1988] Kikuchi, N. and Oden, J. T. (1988). Contact problems in elasticity: a study of variational inequalities and finite element methods. Philadelphia :

SIAM.

[ Kucharik et al., 2010] Kucharik, M., Garimella, R. V., Schofield, S. P., and Shashkov, M. J. ( 2010). A comparative study of interface reconstruction methods for multimaterial ale simulations. Journal of Computational Physics, 229(7):2432 -2452.

[ Kuss, 2008] Kuss, F. (2008). Méthodes duales pour les problèmes de contact avec frottement. PhD thesis, Université de Provence.

[Laursen, 2002] Laursen, T. ( 2002). Computational Contact and Impact Mechanics: Fundamentals of Modeling Interfacial Phenomena in Nonlinear Finite Element Analysis.

Springer-Verlag, Heidelberg.

[ Li et al., 2001] Li, S., Qian, D., Liu, W. K., and Belytschko, T. (2001). A meshfree contact-detection algorithm. Computer Methods in Applied Mechanics and Engineering,:3271 -3292. Advances in Computational Methods for Fluid-Structure

Interaction.

[ Liu et al., 1986] Liu, W. K., Belytschko, T., and Chang, H. (1986). An arbitrary lagrangian-eulerian finite element method for path-dependent materials.

Introduction

This chapter is devoted to the numerical approach used in this work to model large deformations problems. This approach is based on an Eulerian description of the monolithic system governing a multi-domain deformation problem. The used approach proved to be a powerful tool modeling CFD and FSI [START_REF] Feghali | Stable/stabilized mixed formulation for fluid-structure interaction: theory and application[END_REF] problems. Despite its potential, this approach was rarely if never applied to large deformation problems such as metal forming processes (treated later on in this study).

Unlike Lagrangian formulations, in a Eulerian one the mesh itself does not represent the material neither its evolution. This fact leads us to introduce first interface capturing techniques such as distance function. The evolution of a certain domain is followed by solving the convective problem corresponding to the level set function.

Next, the interface capture method is coupled with an heterogeneous anisotropic mesh adaptation technique presented in [Coupez, 2011]. The main goal is to obtain a high precision solution while maintaining the lowest possible number of nodes. It is insured since the mesh is only adapted where the important physical phenomenons appear, while it remains coarse elsewhere. Thus, a large scale simulations can be handled at an affordable cost (i.e. computing time).

To simplify the different steps leading to the monolithic system describing a multidomain large deformation problem, we depict them using a simple flat simulation problem (containing two tools and one deformable sub-domains to generalize later on). One single mesh is used and the different sub-domains are represented as heterogeneities in the global domain. Since our aim is to solve one set of equations and since the different subdomains can easily have different physical characteristics, mixture laws come in handy at this stage to introduce global characteristics defined all over the computational domain using the local ones corresponding to each sub-domain. New additions to this monolithic approach are proposed to help adapt it to large deformation problems. For instance, a novel quadratic law is introduced that can be helpful to represent an additional phase ( detecting the presence of a lubricant for example ). Air trapping is studied as well. Two options are offered to treat this topic depending on the problem need. In addition, some ameliorations are brought to the existent Eulerian approach to insure better results. For every improvement, academic applications are presented to understand both capacities and limitations of this approach.

The chapter is ended with a quick discussion covering the confronted difficulties, what has to be done to prevent them and what remains unsolved to have a more-or-less complete approach treating every aspect of large deformation problems.

Interface Capturing

Aiming to present a "fully" Eulerian approach, the interface capturing should rely on an Eulerian description. Volume of fluids method (VOF) was the most often used. It is based on a mass conservation and consists on convecting the characteristic function over the whole mesh. However, it faces several difficulties such as reconstructing the interface, at each step, specially in 3 d geometries. Later on, Osher introduced the level set method in [START_REF] Osher | Fronts propagating with curvature dependent speed: algorithms based on hamilton-jacobi formulations[END_REF]. In addition to its fluency taking into account every slight topology changes; this method has proven not only its effectiveness to describe precisely the interface but also a quick and smooth transition from 2d to 3d without any further development and complications. This method and other variant ones will be presented next.

Distance function

Consider Ω ⊂ R n the fixed Eulerian computational domain, in which all different bodies are immersed. A standard level set method consists on defining a signed distance function φ verifying (2.1).

|φ( x, t)| = d( x) = min x 0 ∈Γ( x) (| x -x 0 |) ∀ x = (x 1 , . . . , x n ) ∈ Ω (2.1)
where Γ( x) = { x ∈ Ω; /φ( x, t) = 0} is the interface of the studied domain. In this work, φ is chosen to be positive inside the domain and negative otherwise.

If i different bodies (Ω i ) i are immersed in Ω, for each body Ω i (t) a signed distance function φ i (t) is associated and Ω is chosen to cover all bodies:

∪ Ω i (t) ⊆ Ω (2.2)
Both the immersed body and its distance function can depend of the time variable t, since they are supposed to evolve during simulations later on. Note that we do not assign a level set to the air domain. Instead it is deduced by complementarity as follows:

Ω a = Ω -∪Ω i (2.3)
To illustrate this concept, we take the example of one sphere, denote Ω s , immersed in the computing domain Ω. Figure 2.1 depicts the immersed sphere along with a horizontal slice. The isovalues of this slice -0.4 < iso < 0.285 are projected on the bottom of Ω. The latter is a mean to clarify how the distance function associated to the sphere is defined. The projected isovalues confirm that they are negative outside the sphere and positive inside. Thus, the isovalue zero corresponds to the interface of this sphere. Now that the distance function is defined, the integral on the sub-domain itself can be defined as well. Let f be an integrable function on Ω. The integral of f on Ω S can be expressed as the integral of the function f.H on Ω:

ˆΩS f dΩ S = ˆΩ f H(φ s )dΩ (2.4)
where H is the Heaviside function :

H(φ) =    1 if φ > 0 0 if φ < 0 (2.5)

Level Set Method

Assuming that the signed distance function is advected by a given velocity field v, the distance function convection is equivalent to solving the following system on Ω:

   φ t + v.∇φ = 0 φ(t = 0, x) = φ 0 ( x) (2.6)
In fact, the equation (2.6) describes the evolution of the level set function in time. To conserve its stability during the convection, the following feature should be always verified:

|∇φ| = 1 (2.7)
Despite all extensive research and improvements made in stabilization technique; the equation (2.7) was not insured by simply solving (2.6). Since maintaining this gradient property is a must to guarantee the stability of the numerical system, it is then necessary to restore it by solving a Hamilton-Jacobi equation (called also the reinitialization equation):

   φ τ + s(φ)(|∇φ| -1) = 0 φ(τ = 0, x) = φ(t, x) (2.8)
where τ is a virtual time and s(φ) is the sign function defined as follows:

s(φ) =        1 if φ > 0 0 if φ = 0 -1 if φ < 0 (2.9)
The system (2.8) does not affect the position of the zero-value of the level set function; it simply repositions the other isovalues to verify the equation (2.7). The previous reinitialization algorithm was introduced in [START_REF] Sussman | A level set approche fort computing solutions to incompressible two phase flow[END_REF]. The frequency of reinitialization after solving the system (2.6) was found a delicate matter, it can affect considerably the computing time. Authors like [START_REF] Osher | [END_REF]Sethian, 1988, Sussman et al., 1994] urge to reinitialize φ after each time step, others find it sufficient to reinitialize it every other step.

Local Distance function

The level set advection though required, it is not necessary on the whole domain. In computations, the important information are located in the vicinity of the zero isovalue i.e. on the body interface. Therefore, it could be handy to use a function that crops all isovalues in the interface surroundings. Several functions exist in the literature. For example, [Ville et al., 2011] used a sinusoidal function. In this work, a hyperbolic tangent function φ is chosen to narrow the distance function φ:

φ = ε.tanh φ ε (2.10)
where ε represents the thickness of the interface and depends on the mesh size. Figure 2.2 shows that both functions φ and φ match completely in the vicinity of zero. Thus, the use of the hyperbolic filter φ instead of the distance function φ in the interface proximity does not affect at all the results precision. On the contrary, it can be beneficial decreasing the computing domain and time.

Signed distance function Hyperbolic tangent function

Convected level set method

To summarize, a classical level set approach depends on solving both (2.6) and (2.8) systems. Recall that the system (2.8) is not linear. [START_REF] Peng | A pde-based fast local level set method[END_REF] introduced in his work a linearization based on evaluating the gradient at the previous virtual time:

|∇φ| = |∇φ| 2 |∇φ| ≈ ∇φ.∇φ - |∇φ -| (2.11)
Substituting the gradient norm in the first equation of the system (2.8) by the formula (2.11) , the reinitialized convection equation is transformed into:

∂φ ∂τ + s(φ -) ∇φ - |∇φ -| ∇φ = s(φ -) (2.12)
Notice that equation (2.12) is the convection equation of φ with a velocity field

u = s(φ -) ∇φ - |∇φ -|
. The convective reinitialization was introduced by Coupez in [Coupez, 2006]. In addition, he showed that equation (2.12) along with equation (2.6) leads to a convected level set method. In fact, the virtual time step can be decomposed as follows:

∂φ ∂t = ∂φ ∂τ × ∂τ ∂t (2.13
) ∂τ ∂t will be approximated by h ∆t and noted λ, where h is the mesh size and ∆t is the time step [Ville et al., 2011]. Using equation (2.13), equation (2.12) is rewritten to obtain:

∂φ ∂t + λs (φ) (|∇φ| -1) = 0 (2.14)
Adding equations (2.14) to (2.6), the resulting equation is called the convective level set equation:

∂φ ∂t + (v + λu)|∇φ| = λs (φ) (2.15)
Finally, the signed distance function is substituted by the local distance one defined in (2.10), but for the sake of clarity we maintain the same notation φ:

∂φ ∂t + (v + λu).∇φ = λs(φ)g(φ) φ(t = 0, x) = φ 0 ( x) (2.16)
where

g(φ) = 1 - φ ε 2 .

Weak form

Until now, we presented a convective level set method. The remaining part is discretizing the corresponding system. In this section, the finite element discretization is presented briefly.

In order to approximate the level set function as a linear finite element, a space approximation must be determined :

W h = w h ∈ C 0 (Ω h ), w h|K ∈ P 1 (K), ∀K ∈ T h (Ω) (2.17)
Multiplying the first equation of the system (2.16) by w h ∈ W h , the temporal-spatial discretization of the weak formulation is given by :

ˆΩ φ h ∆t w h + ˆΩ φ - h ∆t w h = θ ˆΩ λsg(φ h )w h -ˆΩ(v + λu).∇φ h w h + (1 -θ) ˆΩ λsg(φ - h )w h -ˆΩ(v + λu).∇φ - h w h (2.18)
where φ -= φ(x, t -1) and θ is a parameter used to choose the type of the temporal scheme. In this work, we use a semi-implicit scheme (i.e θ = 1 2 ).

As we are dealing with a convection problem, a stabilization method is needed. Streamline Upwind Petrov-Galerkin method (SUPG) is one of the most popular used methods. Therefore, the space W h is replaced by W h , where the shape function is modified as mentioned in (2.19)

W h = wh , wh|K = w h + τ SU P G v.∇w h , ∀K ∈ T h (Ω) (2.19)
In other words, replacing W h by W h means giving the shape functions an additional artificial weight in the flow direction (see Figure 2.3).

For the determination of the parameter stabilization τ SU P G the reader can refer to [START_REF] Brooks | Streamline upwind/petrov-galerkin formulations for convection dominated flows with particular emphasis on the incompressible navier-stokes equations[END_REF]. To maintain the highest precision possible in this eulerian approach, the anisotropic mesh adaptation is required. Next, we will present the technique implemented in CIMLib to generate and adapt the mesh.

Mesh Generation

A large effort has been devoted to develop mesh generation algorithms. A survey of the literature offers several algorithms capable of generating an unstructured triangle/tetrahedral mesh. The three major algorithms are: Quadtree/Octree described in [START_REF] Tchon | Constructing anisotropic geometric metrics using octrees and skeletons[END_REF], Shephard, 1988, Shephard and Georges, 1991], Frontal Method presented in [START_REF] Bonet | [END_REF]Peraire, 1991, Lo, 1985] and finally the Delaunay Triangulation detailed in [START_REF] Frey | Maillages -applications aux elements finis[END_REF]. The quality of the mesh is crucial for all computations. To ensure an optimal mesh quality, these three algorithms require an optimization step after the generation step.

However, this work lies on a topological mesh. It was first introduced in [START_REF] Coupez | 3-d finite element modelling of the forging process with automatic remeshing[END_REF] and proved in [Coupez, 2000]. Unlike the previous algorithms, the topological mesh optimizes the mesh automatically while its generation. It uses the below definition to construct the mesh and the minimal volume theorem to optimize its quality. The demonstration can be found in [Coupez, 2000, Gruau and[START_REF] Gruau | [END_REF] .

Anisotropic mesh adaptation

Working in an Eulerian framework, an anisotropic mesh adaptation is an efficient means to ensure an accurate prediction of the numerical solution. In fact, an uniform global mesh increases the number of degrees of freedom and may not lead to an improvement of the obtained solution. Starting from a coarse initial mesh, this technique consists on refining it locally in the regions where physical informations requires more nodes. Coupled with the Level set method, the anisotropic mesh adaptation is capable of handling complex geometries. The purpose of this section is to describe briefly the edge based error estimator and the metric construction algorithm implemented in CIMLib. Further details of the metric calculation and the error estimator used to control the generation of the anisotropic mesh are available in [Coupez, 2011, Coupez andHachem, 2013].

2.4. Anisotropic mesh adaptation 51

Edge based error estimator

In large deformations problems, the arise of a new arbitrary surface is common. As mentioned earlier, an isotropic mesh leads to a loss in the accuracy of the desired solution, due to the arbitrary evolution of the interfaces of the different sub-domains. Therefore, an anisotropic mesh adaptation is required. As a matter of fact, discretization errors occur during simulations. The goal of the error estimator is to convert the numerical error into an additional information used to refine locally the mesh where needed. In [Coupez, 2011], the author proposed an edge based error estimator.

Let us begin with recalling some useful notations and definitions.

Consider that v ∈ C 2 (Ω) = V and V h the P 1 finite element approximation such that :

V h = {v h ∈ C 0 (Ω), v h | k ∈ P 1 (K), K ∈ K}, where K represents the set of element of the mesh. Furthermore let us denote by X = {X i ∈ R d , i = 1, • • • , N} the set
of nodes, and by U i the projection of the solution U at the node X i given by the Lagrange interpolation. Finally let X ij = X j -X i and U ij = U j -U i be respectively the length between the nodes i and j and the variation between the solution evaluated on these nodes.

Since a P 1 approximation is used, the solution is continuous and can be written as:

U j = U i + ∇u h X ij (2.21) so, ∇u h X ij = U ij (2.22)
Notice that the gradient of u h is a piecewise constant vector i.e. constant vector on each element. Though, the gradient field is discontinuous from element to element, it is continuous along the element edges since it depends only on the variation of the field value at the extremities (see equation (2.22) ). Most error estimators are based on the Hessian approximation, nevertheless the author of [Coupez, 2011] proposes an error based on the gradient.

Using identities ( 2.21 -2.22) and the mean value theorem, the following inequality is determined :

∇u h .X ij -∇u(X i ).X ij ≤ max Y ∈[X i ,X j ] H(u)(Y )X ij .X ij (2.23)
where H(u) is the Hessian solution. Equation (2.23) is the main reason behind the choice of the following gradient based error estimator:

e ij =| G ij X ij | (2.24)
where G ij = G i -G j is the difference between the gradient defined respectively on the nodes i and j.

Since the gradient is only determined on the elements, an additional step is needed to recover its values on the nodes. It is accomplished by solving the following equation:

G i = argmin G   j∈Γ(i) (G -∇u h ) .X ij 2   (2.25)
where

Γ(i) = {j ∈ N , ∃K ∈ K, X ij ∈ K}.

Metric construction

To pilot such a technique, a metric is needed. A metric is a symmetric defined positive tensor. It can be transformed into a diagonal one using:

M = R         λ 1 0 • • • 0 0 . . . . . . . . . . . . 0 0 • • • 0 λ d         R T (2.26)
where R is a rotation tensor containing the eigenvectors of M and (λ i ) i=1,••• ,d represent the eigenvalues of the metric. The eigenvalues (λ i ) i=1,••• ,d are strongly related to the shape and the size of the obtained mesh. It will be inversely proportional to the square of the mesh size. The orientation of the mesh elements will be determined by both R and R T . During computations, a generated mesh formed by equilateral triangles is the most suitable. That is why the metric notion will introduce a distance between the set of the nodes with the aim of maintaining the norm close to one :

X ij M = MX ij , X ij = 1 (2.27)
in which X ij is the vector formed by the nodes i and j.

Following the work of [Coupez, 2011] , the proposed metric is constructed on each node of the mesh verifying (2.27). In fact, the author demonstrates that a unit metric can be built via an affine transformation A ij . The latter transformation will only adjust each edge length to be equal to 1. Therefore, the metric can be written as M = t A ij A ij (defined on the node i).

When summing (2.27) on Γ(i) = {j ∈ N , ∃K ∈ K, X ij ∈ K}, we obtain the following identity:

2.4. Anisotropic mesh adaptation 53 M :   j∈Γ(i) X ij ⊗ X ij   = |Γ(i)| (2.28)
The tensor in (2.28) on the node i is called the length distribution tensor defined (see Figure 2.5) :

X i = 1 | Γ(i) |   j∈Γ(i) X ij ⊗ X ij   (2.29) ij X i X j X Figure 2.5:
The length distribution tensor for the node i is determined from the collected edges Since X i is symmetric defined positive, when at least d non aligned edge vectors exist the metric is defined as :

M = 1 d X i -1 (2.30)
In order to introduce the error variations calculated in (2.24) to the metric defined in (2.30), a new edge length is defined based on a stretch factor s ij . The introduction of the latter factor leads to a new definition of the error given by e ij = s 2 ij e ij . More precisely, the stretch factor leads to generate a mesh under the constraint of a fixed number of edges. In fact, by scaling the edges by s ij , the number of edges is multiplied by s -p ij . Thus for a given number of edges A, the proposed stretching factor is computed as follows:

s ij = λ e ij 1 p (2.31) in which Chapter 2. Monolithic approach tools 54 λ =           i j∈Γ(i) e p p+2 ij   A         p+2 p
Finally the optimized metric is given by:

M i = 1 d   j∈Γ(i) s 2 ij X ij ⊗ X ij   -1
(2.32)

A full review of the metric construction, lemma, propositions and proofs can be found in [Coupez, 2011]. Figure 2.6 shows the mesh obtained when using the technique described above for two different geometries. In Figure 2.6 (a-b) we choose a cube, it will be used latter on to study the case of a flat bunch simulation. To prove that this technique can handle more complex geometries, Figure 2.6 (c) shows the adapted mesh of the upper tool of a connecting rod. We point out that in these cases the mesh is adapted on the level set φ.

Extension to multi-component field

Many simulations require to adapt the mesh over a multi-field. Even though, the most important aspect is to capture the evolution of the new arbitrary surface generated by the level set motion; it come in handy to adapt the mesh over the mixture law and/or the velocity field. Though metric intersection methods are commonly used to define one metric when several ones are present [Loseille, 2008]; the author of [START_REF] Coupez | Solution of highreynolds incompressible flow with stabilized finite element and adaptive anisotropic meshing[END_REF] chose to extend the above theory to suit a multi-component field adaptation. This approach was validated with a driven cavity problem. the mesh was adapted not only using the velocity norm, but also using all velocity vector components.

Mainly this approach requires a global vector containing all fields to be adapted.

Denote u = {u 1 , • • • , u m } the vector containing those fields. Using the error estimator described in (sec. 2.4.1), the computed error is now a vector given by:

e ij = e 1 ij , e 2 ij , . . . , e m ij (2.33)
Each component of the error vector is obtained by applying the already described theory on each component of u. To provide a scalar value to the error estimator; [START_REF] Coupez | Solution of highreynolds incompressible flow with stabilized finite element and adaptive anisotropic meshing[END_REF] proposes to evaluate the vector e ij by its L 2 or L ∞ norm.

Mechanical problem

This section is devoted to elaborate an appropriate mathematical model capable of describing the large deformation problem. The desired model is written based on an Eulerian formulation of the fundamental mechanics equations. First, the system is adapted to suit this problem needs. The monolithic system is established from the weak formulation written in a velocity/pressure variables. Finally, the Finite Element formulation is explained using a Variational Multi-Scale method for stabilization.

Governing equations

In the framework of large deformation problem, let us consider a simple forging process illustration for a better understanding of the mechanical problem (see Figure 2.7). Since an Eulerian formulation is adopted, the computational domain Ω ⊂ R d (d = 2 or 3) is chosen in a way to cover all the present sub-domains, even the air. The deformable body, denoted by Ω d , is deformed when the upper die moves uniformly in the -z direction. Usually, the lower die remains fixed but nothing prevent it to move in the z direction. In this work, the temperature is assumed to be constant during the simulation. As a result, the system of equations is governed by the mass conservation equation and the momentum equation, omitting the heat transfer equation. The mass conservation equation is expressed by: where ρ is the material density and v represents the velocity. Since the computational domain is assumed to be incompressible, the equation ( 2.34) is reduced to:

∇.v = 0 (2.35)
The momentum equation is given by:

ρ ∂v ∂t + v.∇v -∇.σ = ρg (2.36)
in which g is the gravity and σ = s -pI d is the Cauchy stress tensor, where p is the pressure and I d is the identity matrix. The deviatoric component s is defined in terms of each sub-domain Ω i . In Ω d , s follows a Norton-hoff law:

s = 2K d √ 3 ε m-1 ε(v) (2.37)
in which K d is the consistency of the deformable body and m (0 < m 1) is the strain rate sensitivity coefficient. The strain rate tensor

ε(v) = 1 2 (∇v + t ∇v), while ε = 2 3 ε(v) : ε(v).
In Ω a which denotes the air, s obeys to a Newtonian fluid law:

s = 2η a ε(v) (2.38)
where η a represents the air viscosity. In the rigid bodies the following equation should be verified: Chapter 2. Monolithic approach tools 58

ε(v) = 0 (2.39)
In fact, equation (2.39) is not considered within the system. Nonetheless, a different procedure is taken into account to verify this equation. In fact, to be consistent with what was presented above, the stress deviatoric part of the rigid bodies is written as

s = 2K s ε(v) (2.40) in which K s , the tools consistency, should respect K s K d .
In addition to the previous restriction, the signed distance function φ r is computed at each time step. Thus the upper die is placed into it right place. Under the latter conditions, equation (2.40) is equivalent to equation (2.39).

Weak formulation

In the previous section, the equations governing the physical problem were presented. Now the aim is to combine all these sub-systems in order to obtain the monolithic system defined all over the domain. To begin with, the computational domain Ω is defined to cover all heterogeneous environment including the air ( see sec. 2.2.1 ). To obtain the monolithic system, using the weak formulation is mandatory. Thus, we must first define Sobolev spaces in which the solution is approximated:

V = v, v ∈ H 1 (Ω) d |v = g on Γ (2.41) Q = p, p ∈ L 2 (Ω), ˆΩ pdΩ = 0 (2.42) W = w, w ∈ H 1 (Ω) d |w = 0 on Γ (2.43)
In the spaces defined above, the scalar product is given by: (v, w) Ω = ˆΩv : w (2.44)

A variational formulation lies mostly on the Green formula. The problem amounts to find (v, p)

∈ V × Q/ ∀w ∈ W, ∀Ω i ⊂ Ω :        ρ ∂v ∂t , w Ω i + (ρv.∇v, w) + (s : ε(w)) Ω i -(p, ∇.w) Ω i = (ρg, w) Ω i (∇.v, q) Ω i = 0
(2.45) 2.5. Mechanical problem 59

Monolithic system

At this stage, the monolithic system is defined all over the computational domain Ω by adding the previous weak formulations corresponding to each sub-domain Ω i .

Find (v, p) ∈ V × Q/ ∀w ∈ W:     ρ s H(φ r ) + ρ d H(φ d ) + ρ a H(φ a ) ρ     ∂v ∂t + v.∇.v, w +     2     K s H(φ r ) + K d √ 3 ε m-1 H(φ d ) + η a H(φ a ) η     ε(v) : ε(w)     Ω -(p, ∇.w) Ω = (ρg, w) Ω (2.46)
where the introduction of the Heaviside function allows to preserve the physical characteristic of each heterogeneity. In fact the terms ρ and η in equation ( 2.46) denotes the mixture law, it will be detailed and discussed in (sec. 2.6) and (sec. 2.7.1.2).

Finite Element Formulation

Since solving system (2.46) by finding an analytical solution is quite impossible, it is necessary to approximate its solution by a space and time discretization.

Spatial discretization

A finite element discretization of the monolithic system above begins by approximating the functional infinite spaces (2.41 -2.43) by finite size sub-spaces. Let V h ⊂ V, Q h ⊂ Q, W h ⊂ W a family of sub-spaces defined as follows:

V h = v h | v h ∈ C 0 (Ω) d , v h|k ∈ P 1 (K) d , ∀K ∈ T h (2.47) Q h = p h | p h ∈ C 0 (Ω) d , p h|k ∈ P 1 (K), ∀K ∈ T h (2.48) W h = w h , w h ∈ C 0 (Ω) d |w h = 0 on Γ h (2.49)
In the theory of finite elements, the approximation of the solution v (for instance) on an element is determined using only its nodes values (called degrees of freedom). Then, the solution is approximated by

∀K ∈ T h v h|k = D i=1 v i N i p h|k = D i=1 p i N i (2.50)
in which i represents the node and N i the shape function on this node. The numerical system is given by inserting equations (2.50) in the system (2.46). However the solution obtained by just solving this system is not accurate and represents many oscillations due to the convective terms.

Variational Multi-scale Stabilization

The accurate simulation of such problem remains a major challenge. Indeed, the convective term in (2.46) disrupts highly the solution. A variational multi-scale approach is chosen to stabilize the solution. This stabilization technique was originally proposed by [START_REF] Hughes | The variational multiscale method-a paradigm for computational mechanics[END_REF]] . The following description of this technique is for a finite element formulation using a tetrahedral mesh. Several development and publications of VMS methods adapt this technique to other type of formulation such as the finite volume for example [START_REF] Dahmen | Adaptive Multiresolution Finite Volume Discretization of the Variational Multiscale Method: General Framework[END_REF]. Inspired by the work of [Codina et al., 2007], [START_REF] Hachem | Stabilized finite element method for incompressible flows with high reynolds number[END_REF] developed a rigorous derivation of the VMS method in 2009. Here, a brief description of the VMS method implemented in CIMLib is presented.

The anatomy of the method relies on a two-scale decomposition of the velocity/pressure fields. Thus the following sub-spaces are introduced:

V k = v h | v h|k ∈ P 1 (K) d ∩ H 1 0 (K) d , ∀K ∈ T h (2.51) Q = p h | p h|k ∈ P 1 (K) ∩ H 1 0 (K), ∀K ∈ T h (2.52)
The solution is decomposed into a coarse scale and a fine scale:

v h + v and p h + p (2.53)
Due to the orthogonality of spaces, the system of equations is split into two parts called coarse/fine scale. After Integrating by parts and considering Dirichlet boundary conditions, the coarse scale is written as :

ρ ∂ (v h + v) ∂t , w Ω + ρ ((v h + v) .∇ (v h + v) , w) Ω + (2ηε (v h ) : ε (w)) Ω -((p h + p) , ∇.w) Ω = (ρg, w) Ω (2.54) 2.5. Mechanical problem 61 (∇. (v h + v) , q) = 0 (2.55)
while the fine scale is given by:

ρ ∂ (v h + v) ∂t , w Ω + ρ ((v h + v) .∇ (v h + v) , w) Ω + (2ηε ( v) : ε ( w)) Ω -(p, ∇.(w h + w)) Ω = (ρg, (w h + w)) Ω (2.56) (∇. (v h + v) , q) = 0 (2.57)
Solving the previous system of equations (2.54) and ( 2.56) is similar to the one used by a Mini Element formulation [Basset, 2006]. It normally begins by solving the fine scale on the sub-grid without being tracked in time. The VMS method in CIMLib adopts many assumptions to simplify the problem.

The non-linear term in equation ( 2.56) is approximated using only the coarse scale part:

(v h + v) .∇ (v h + v) ≈ v h .∇v h + v h .∇ v (2.58)
Tough the resolution of the fine scale equation using both variables (v, p) provides a better stabilization, [Hachem, 2009] follows the work of [START_REF] Franca | Pressure bubbles stabilization features in the stokes problem[END_REF]] and approximates the pressure in equation (2.56) by:

p = -τ c ∇.v (2.59)
For the definition of τ c , he adopts the one proposed by [Codina and Principe, 2007]

τ c =   ν 2 + c 2 |v| k c 1 h 2   1/2 (2.60)
where c 1 and c 2 are two constants, independent of h. h is the characteristic length of an element and ν represents the kinematic viscosity. Now equation (2.59) can be integrated directly into the coarse scale equation (2.54). Thus, the fine scale equation is reduced to the following:

ρ (v h .∇ v, w) Ω + (2ηε ( v) : ε ( w)) Ω = (R m , w) Ω ∀ w ∈ W (∇. v, q) = (R c , q) Ω ∀q ∈ Q (2.61)
where R c designs the term -∇.v h and R m gathers all terms considered in the second Chapter 2. Monolithic approach tools 62 member of the equation :

R m = ρg -ρ ∂v h ∂t -ρ (v h + v) .∇v h -∇p h , v (2.62)
Many methods were developed capable of solving the latter system. A full review can be found in [Hachem, 2009]. In CIMLib, the used approach consists on employing the bubble functions as the shape functions on each element. Moreover since the problem is a convection problem, the shape functions are modified to an upwind shape functions. Thus, the solution of (2.61) is given by :

v = 1 ρ (v i-1 .∇b, b * ) + (2ηε (b) : ε (b)) (R m , b) (2.63)
where b represents the bubble shape functions.

As for the coarse scale equation, the non linear term is approximated by keeping the terms of the first order at the i th iteration as:

(v.∇v, w) = v i .∇v i-1 , w + v i-1 .∇v i , w -v i-1 .∇v i-1 , w (2.64)
where u i-1 is the previous know Newton-Raphson's iteration. After extracting ṽ the solution of the fine scale problem, equation (2.54) is integrated by parts and ṽ is substituted by it value, the coarse scale is expressed by equations (2.65 -2.66).

ρ ∂v h ∂t , w h Ω + ρ (v h .∇v h , w h ) Ω + (2ηε (v h ) : ε (w h )) Ω -(p h , ∇.w h ) Ω + (∇.v h , q h ) -(ρg, w h ) A + K∈T h τ k ρ ∂v h ∂t -ρv h .∇v h + ∇p h -ρg, ρv i-1 h ∇.w h Ω B + K∈T h τ k ρ ∂v h ∂t -ρv h .∇v h + ∇p h -ρg, ∇.q h Ω C + K∈T h (τ c ∇v h , ∇.w h ) Ω =0 D (2.65) (∇.v h , q h ) - K∈T h (τ k R m , ∇.q h ) Ω = 0 (2.66)
Comparing with a standard Galerkin formulation, the terms denoted A in equation 2.5. Mechanical problem 63 (2.65) refers to a standard Galerkin formulation, the terms denoted B represents the upwind stabilization terms and the terms denoted C refers to the pressure stabilization terms. Where the last term denoted by D is induced by the enrichment of the pressure spaces. The appearance of this new term will not interfere with the existence and unicity of the solution i.e. equation ( 2.65) satisfies the inf-sup condition for the velocity and pressure interpolations.

Stabilization parameter

Following the same analogy of [Hachem, 2009], this section presents the different evaluation of the stabilization parameter τ k . A comparison will be presented between a condensed Mini Element formulation and the VMS stabilization technique. To eliminate all ambiguities matrix notations shown are similar to those in a Mini Element formulation. Using a the latter formulation, the local matrix formulation of the problem presented above is given by:

     A vv A t vb A t vp A vb A bb A t bp A vp A bp 0           v h v h p h      =      B v B b B p     
(2.67) when using a regular bubble function, (i.e. the bubble function is centered in the element) the terms in velocity/bubble will be vanished. Solving (2.67) consists into expressing v h in function of v h and p h . Then by injecting it in the third and first line, the final matrix system is obtained. Therefore the inversion of the matrix A bb is quite necessary when adopting such formulation. The genuinely of the VMS method is it capacity to tune automatically the stabilization parameter. In fact, the matrix formulation of the problem presented above is given by:

  A vv A t vp A vp A pp     v h p h   =   B v B p   (2.68)
where the additional terms appearing in equations (2.65 -2.66) are explicitly present in the matrix according to their correspondence. Reasoning by analogy, the A pp contribution in the VMS formulation is equivalent to the A bp A -1 bb A t bp contribution after the condensation of the system (2.67). In what follows, the choice of τ k is discussed. This choice has been the subject of several studies. Many works evaluate τ k by computing the local Reynolds number. In addition, it depends highly of the problem type: convection or diffusion [Medic and Mohammadi, 1999]. In CIMLib, the adopted parameter is computed from Fourier analysis proposed by [Codina, 2002]. It combines both regimes and evaluated as follows:

τ k =   c 1 ν h 2 2 + c 2 | v | k h 2   -1 2 (2.69)
where c 1 and c 2 are two algorithmic constants chosen equal to 4 and 2 respectively for linear elements. Notice that the previous stabilization parameter is set to be a static formulation. For a time dependent problem, [Tezduyar and Osawa, 2000a] proposed the following stabilization term

τ k =   2 ∆t 2 + c 1 ν h 2 2 + c 2 | v | k h 2   -1 2 (2.70)
For anisotropic meshes with highly stretched elements, the definition of the mesh size within each element is still an open problem [START_REF] Harari | What are c and h?: Inequalities for the analysis and design of finite element methods[END_REF]. Different definitions are proposed in [START_REF] Principe | On the stabilization parameter in the subgrid scale approximation of scalar convection-diffusion-reaction equations on distorted meshes[END_REF] and [Tezduyar and Osawa, 2000a] for instance. In CIMLib, the mesh size element is computed in terms of the velocity v, its norm and the derivative of the shape function N as follows:

h = ne i v v ∂N ∂x i (2.71)
where n e denotes the number of nodes per element.

Mixture laws

After defining all sub-domains geometrically with a distance function, the physical characteristics of each sub-domain are determined by mixture laws. In this section, different types of mixture laws are described briefly. The impact of using each of these mixture laws is discussed in the next section. In fact, to represent density and viscosity discontinuities over the interface a Heaviside function is needed. It is equal to 1 if the mesh element is fully filled and equal to 0 if it is not. If the material interface passes through an element (i.e. mesh element partially filled), the Heaviside function H(φ) is computed as follows:

H(φ) |k = 1 2 1 + φ + k |φ| k then H ∈ [0, 1] (2.72)
where φ + k is the sum of all signed distance functions evaluated at each node of the element k and |φ| k is the sum of the absolute value of φ.

By choosing to calculate the Heaviside function in this manner, the mixture law will be similar to the Volume of Fluid Method (VOF) where the fraction of the volume will 2.6. Mixture laws 65 be H. Then, the mixture law is given by:

η = n i=1 H(φ) |k η i (2.73)

Linear mixture law

In computations, to achieve numerical robustness, a smoothed Heaviside is often used:

H ε (φ) =            1 if φ > ε 1 2 1 + φ ε if | φ | ε 0 if φ < -ε (2.74)
where generally ε is taken as the interface thickness interface ( see sec. 2.2.2.1 ). The regularized Heaviside will ensure a lighter transition between the different values of η corresponding to the different heterogeneity. Notice that, in Figure 2.8, the difference in the transition is very clear to the naked eye.

On the left, the P0-mixture law (2.73) is computed relatively per element. Thus depends highly on the refinement of the mesh. On the right, the mixture law is computed on each node. It does not relay on the mesh refinement as the P0-mixture law does. It is called the Linear mixture law and is given by: In this work, a new mixture law is implemented in CIMLib. This section will just describe the method, while the utility of this law will be discussed in chapter 5. It will be used to manage the sliding between two different heterogeneities. The main idea is to introduce a third viscosity between the two other heterogeneities. The newly added viscosity represents the viscosity of the lubricant. The advantage of this method is that lubricant is added implicitly. Thus no need to deal with it displacement. It is given by the mixture law.

η(H ε ) = n i=1 η i × H ε (φ i ) (2.75) (a) (b)
To obtain the desired mixture law, one must choose a function admitting a minimum. Therefore, the general form of a quadratic function is considered (see equation (2.76)). This minimum will represent the viscosity of the lubricant.

Φ : [0, 1] → R H ε → η(H ε ) = a × H 2 ε + b × .H ε + c (2.76)
where

a = -2(η L -η 0 ) + (η L -η 0 )(η L -η 1 ) + (η 1 -η 0 ) b = 2(η L -η 0 ) -(η L -η 0 )(η L -η 1 ) c = η 0 (2.77)
In equation (2.77), η 0 and η 1 represents respectively the consistence of the deformable body and the tool, while η L represents the lubricant viscosity. Figure 2.9 depicts the difference between a linear and a quadratic mixture law. This simulation consists on compressing a workpiece between two flat dies, as described in (Figure 2.7). The upper die presses the workpiece with a constant velocity v imp =0,25 m/s, while the lower die stays fixed. As mentioned earlier, the material obeys the Norton-Hoff behavior law with a strain rate sensitivity coefficient fixed equal to 1. Table 1 shows the different parameters used in this simulation.

The initial height of the workpiece is 0.5 m .

After several tests, some deficiencies were noticed such as mass loss. In order to solve this deficiency, a short parametric analysis is presented next. Several parameters are proved to have a great impact on the mass loss appearance: the time step, the mixture law choice, and the boundary condition imposition.

Several improvements are introduced and presented next highlighting their positive effect on minimizing the mass loss:

1. Since the use of a really small time step is out of question do to the time of the simulation, we introduced a variable time step in CIMLib.

The heterogeneous time step is introduced using a formula in terms of the workpiece initial height h 0 , the workpiece height h i at the current increment i and the initial time step Δt 0 :

∆t i = h i h 0 ∆t 0 Figure 2
.10 shows better mass conservation around 14% using the variable time step. Figure 2.12: The effect of imposing the velocity on the inner of the dies 2.7. New additions and Numerical Applications 69 2. we changed the mixture law from P0 to P1 interpolation. Figure 2.11 shows the result of this change. An improvement of 8% is noted.

3. Last but not least, we treated the imposition of the dies boundary condition differently. Since a P1 mixture law is used, a critical transition zone is present when the contact is established. That is why, an inner part of the die is defined taking into account the transition zone. For every die defined on L × l, its inner part is defined at an appropriate distance of the interface e.g. (Lε) × (lε) where ε is the same parameter used for smoothing the Heaviside function in the mixture law. Now, the velocity is no longer imposed on the whole die, it is only imposed on its inner part.

All these improvements combined, minimize remarkably the mass loss. From now on and unless mentioned otherwise, they will be used in all applications presented in this work.

Weight of mixture laws

The issue remaining is the shape of the forged piece. As Figure 2.13 shows, using equation ( 2.75) to compute the mixture law affects the form obtained. In fact using such law allows higher consistency to outweighs the smaller one. It leads to preserving the initial trace of the deformable body on the tools. In addition, it prevents the deformable body from entering perfectly in contact with the tools during the simulation as shown in Figure 2.13. Therefore, we need to consider new mixture laws by changing the weight of the coefficients η i (associated to the sub-domain i). .14 shows a comparison between the four mixture law presented above. As mentioned previously when using linear coefficients, the body with the higher consistency outweighs the body with the smaller one. But when using 1 η i , the smaller consistency outweighs the higher consistency. The geometric mean and the logarithmic mixture laws are the best reducing the gap between the different consistencies. Still, the logarithmic weight is the only one capable of centering the transition between two different values on the interface.

1 η(H ε ) = n i=1 1 η i × H ε (φ i ) (2.78) n η(H ε ) = n i=1 n √ η i × H ε (φ i ) (2.79) log 10 η(H ε ) = n i=1 log 10 η i × H ε (φ i ) (2.80)
Next, we decide to launch the flat bunch simulation using only three different mixture laws defined in (2.75), (2.78) and (2.80). The geometric mean coefficient is not used since it is very similar to the logarithmic coefficient.

Figure 2.15 shows clearly that the inverse coefficient mixture law is poorly adapted to our case. When the deformable body does not perfectly touch the tool ( dashed lines ) i.e. the gap is filled with air. Knowing that the air has the smallest viscosity the inverse mixture law will give it advantage on the other sub-domains. Numerically, this means that these gaps will be preserved during the whole simulation and will affect the form of the final piece ( which is not the case in reality ).

Figure 2.16 confirms that the logarithmic weighs of the mixture law is the one with the best results. It offers a smooth transition without any remarkable imperfections. The transition between different consistencies are perfectly centered on the interface between the sub-domains. In addition, Figure 2.17 depicts the top view of the deformable body obtained in a 3d flat bunch simulation. A comparison between the classical linear mixture law and the logarithmic one displayed on the zero isovalue confirms that the contact problem confronted using the linear mixture law (Figure 2.13) faded away when using the logarithmic one : At the end of the simulation the deformable body is perfectly flat.

For all these reasons, the logarithmic weighs mixture law will be used in all our simulations.

Transport stabilization

In more than one occasion, stabilization problems were confronted in the convection solver levellerT. They were more noticeable in cases displaying highly anisotropic meshes. These numerical oscillations are mainly represented by positive level set values, indicating new material appearance in unexpected places. Figure 2.18 shows one example exhibiting instabilities in the computational domain. The physical characteristics are not mentioned since the problem is purely numerical. After lengthy analysis, we found the origin of this instability. It is nothing other than the mesh size defined in the solver. As shown in Figure 2.19, in CIMLib, the size of every stretched element is determined as the smaller edge of the circumscribed rectangular. We propose another definition for the mesh size noted h previously. It can be defined as the element altitude (Figure 2.19). Knowing that it is not the most optimal definition, it remains logical since it is determined in the material flow direction. To clarify the impact of the new definition, we use a 2D-unity domain containing 6 stretched elements illustrated in Figure 2.20. The different mesh elements were computed and presented in both Figure 2.20 and Table 2.2. Notice that the new value Altitude min is always smaller or equal to the older Edge min value. It is particularly smaller when the element is stretched (for the elements numbered 2 and 5 for instance). The new smaller mesh size Altitude min stabilize the solver by respecting better the relationship binding h and the virtual time step τ (defined earlier in the reinitialized convected level set equation in sec. 2.2.2.2). After this adjustment, the same 3D case (Table 2.1) presenting instabilities in Figure 2.18, is launched again to check if the problem persists. A logarithmic weight is used for the mixture law coefficients. Notice that, the mesh is anisotropically adapted/refined on the interface of the body and is coarse elsewhere. The body is deformed uniformly everywhere; we should mention though that at t = 2.975 s, the deformable body shows an impression and a sort of out of plane displacement. This caused by the small tools along with the boundary conditions imposed : v. n = 0. The deformable body has no way to move but in this manner. As approved in Figure 2.21 and even though the anisotropic mesh presents highly stretched elements, no numerical oscillations were detected. Note that this new mesh size definition will be applied in all following simulations to prevent any instabilities.

Air trapping

During a multi-domain problem, the lubricant or the air can be trapped between different interacting bodies. This section addresses the trapping issue between a rigid body and a deformable one. Figure 2.22 illustrates a problem in which the upper tool presents a porosity. Notice that the air ( blue domain ) fills the porosity between the upper tool ( red domain ) and the deformable body ( orange domain ) as presented in Figure 2.23 .

To treat this topic in a Lagrangian approach, an additional function must be introduced to the code, allowing to locate pockets of air or lubricant. The purpose of such a function is to take into account the impact exerted by the volume variation of the lubricant or the trapped air. In fact, the volume variation introduces a pressure variation exerted to the deformable body preventing it to fill the porous tools.

However in our Eulerian approach, their is no need to implement a function defining the air presence. The computational domain is defined to cover all heterogeneity including Similar to the previously presented case, the kinematic of the tools is imposed on their inner defined as follows :

Γ inner t (x) = x ∈ Ω, φ s ( x, t) > ε 2 (2.81)
Figure 2.24 shows the evolution of the isovalue zero of φ d when considering the boundary conditions (2.81). This imposition generates a great pressure in the trapped air. Combined with the incompressible assumption, they prevent the deformable body from taking the tool shape. In addition, the air pockets are clearly detected in Figure 2.24(c-d). When simulating industrial cases with complex geometries, these problems (trapped air and pockets) are faced frequently.

Following this discussion, two paths can be proposed depending on what is more interesting for the study. The first focuses on the pockets formation due to the trapped air. The second, discussed in the next section, consists on creating vents to ensure the filling of the tools.

Porous boundary conditions

To overcome the trapping effect, a new boundary condition called porous boundary condition, is defined. The kinematic of the tools is imposed accordingly to the new boundary condition allowing the escape of the air.

To define the porous boundary condition, we define the following set of nodes denoted Γ φ t :

Γ φ t (x) = { x ∈ Ω, φ d ( x, t) > -0.9 × ε} (2.82)
Γ φ t represents a fictitious domain slightly bigger than the actual deformable body. The next step is to define the porous boundary condition as the intersection of Γ φ t and the classical boundary condition:

Γ p = Γ inner t ∩ Γ φ t (2.

Conclusion 81

Figure 2.25 shows that the set Γ p evolves in time. The evolution is expected since the contact zone changes as well. This porous boundary condition shows a great capacity to adapt with the changing geometry of the deformable body.

To complete this discussion, the evolution of the deformable body is followed in Figure 2.26. The new boundary condition has proven to be a good fit to our problem. The air escapes allowing the body to deform properly filling the tool porosity.

Conclusion

The present chapter is divided into two main parts. The first details the used Eulerian monolithic approach. The different tools are presented beginning by the convected level set method insuring the tracking of the interface. The mesh generator is described as well resulting in a single anisotropic mesh all over the global domain guaranteeing better precision. Since one mesh is applied per computational domain, the different interacting bodies are represented as heterogeneities via mixture laws. This part is ended by the mechanical problem coupling the Navier-Stokes equations with the level set convection equations.

The second part is mostly dedicated to our contribution and improvements to this existent approach. A new quadratic mixture law is introduced. It plays a key role in large deformation problems exhibiting additional phases like the presence of a lubricant for instance. Several technical ameliorations are proposed also such as i) the choice of the mixture laws weight, ii) the choice of interpolation (P0 vs. P1), iii) the choice of an adaptive time step and iv) a new choice of the mesh size. These enhancements are proven to be of a great importance when it comes to stabilization and volume conservation. A last addition deals with the way to treat the existent air in our domain. The method by default traps the air inside. We offer the option to let it clear out using a new manner to impose boundary conditions. It is equivalent to creating vents in the tools for instance. The air trapping condition remains available if one wants to follow bubbles formation in the domain.

To conclude, the existent approach needed several additions to be adapted to large deformation problems involving multiple interacting bodies. Moreover, a great care is needed to choose the optimal parameters since it shows a high sensitivity. Combined with the new improvements, it shows a great capacity to respond to our needs. From now and unless specified otherwise, all ameliorations proposed in this chapter will be used for every application to ensure the optimal results. Note that several development remains necessary. Even though the contact detection is managed automatically, the contact resolution is not taken into account. Sticking friction and sliding friction should be studied likewise. The latter developments will be treated in details in chapter 5.

Since one of the objectives is to work in a highly scalable parallel environment, an exhaustive study is presented in the next chapter. It is entirely dedicated for the parallelism of the CIMLib library where the performance of this approach is inspect closely.

Résumé en Français

Le présent chapitre est divisé en deux parties principales. La première détaille l'approche Eulerienne monolithique utilisée. Les différents outils sont présentés commençant par la méthode de convection de Level Set (assurant le suivi de l'interface). La génération d'un seul maillage anisotrope défini sur le domaine de calcul garantissant une meilleure précision est décrite aussi. Puisqu'un maillage est appliqué sur tout le domaine, les différentes pièces interagissant sont représentées comme des hétérogénéités par les lois de mélange. Cette partie prend fin avec la description du problème global couplant les équations de Navier-Stokes avec les équations de convection de la Level Set. La deuxième partie est principalement dédiée à nos contributions et améliorations proposées sur cette approche existante. Une nouvelle loi de mélange quadratique est introduite. Elle joue un rôle clé dans les problèmes de grandes déformations présentant des phases supplémentaires, comme la présence d'un lubrifiant par exemple. Plusieurs améliorations techniques sont proposées également telles que i) le choix du poids des lois de mélange, ii) le choix de l'interpolation (P0 vs P1), iii) le choix d'un pas de temps adaptatif et iv) un nouveau choix de la taille de maille. Ces améliorations sont prouvées être d'une grande importance en ce qui concerne la stabilisation et la conservation du volume. Un dernier ajout traite la modélisation de l'air dans notre domaine de calcul. La méthode par défaut emprisonne l'air à l'intérieur. Nous offrons la possibilité de le laisser évacuer le domaine à l'aide de nouvelles conditions aux limites. C'est équivalent à la création des évents dans les outils par exemple. La condition de piégeage de l'air reste disponible si l'on veut suivre la formation de bulles au cours du temps. Pour conclure, l'approche existante avait besoin de plusieurs ajouts pour être adaptées aux problèmes de grandes déformations. En outre, un grand soin est nécessaire pour choisir les paramètres optimaux. Combinée avec les nouvelles améliorations, l'approche montre une grande capacité à répondre à nos besoins. Notez que plusieurs développements restent nécessaires. Bien que la détection de contact soit gérée automatiquement, la résolution de contact n'est pas prise en compte. La modélisation du frottement glissant et collant doivent être étudiés de même. Ces derniers développements seront traités en détails dans le chapitre 5. Puisque l'un des objectifs est de travailler dans un environnement massivement parallèle une étude exhaustive est présentée dans le chapitre suivant. Il est entièrement dédié au parallélisme de la bibliothèque CimLib où la performance de cette approche est inspectée de près.

[ Ville et al., 2011]

Introduction

The last few years has experienced a significant progress in computer technology. This progress has influenced mostly the processor power. Modern processors are formed by multi-cores. Hence, the execution of multiple instructions is available.

Simultaneously, the complexity of mathematical models has increased as well. These models can describe complex physical problem such as weather forecast, aerospace applications, crash test vehicles, material forming and much more applications. The resulting high computational demands inspired scientific programmers to introduce parallel computing concept.

The use of parallelism in scientific computing reduces considerably computational time when comparing it to the sequential time. Thus, handling larger scale physical problems are now affordable. For instance, the number of freedom can be increased to obtain a more accurate solution without a drastic increase of the computing time.

To be beneficial, parallelizing algorithms has to be treated with a great care. In fact, the efficiency of parallelism depends highly on the computer architecture and other software aspects such as the load balancing.

This chapter provides a quick overview of the parallelism in general and focuses on CIMLib's parallel environment in particular. We begin by explaining some important concepts and terminologies to simplify the topic for non-initiative readers. For more details, the readers are invited to consult [START_REF] Magoulès | Calcul scientifique parallèle -Cours, exemples avec Openmp et Mpi, exercices corrigés[END_REF], [START_REF] Eijkhout | Programming on Parallel Machines: GPU, Multicore, Clusters and More[END_REF] and [START_REF] Rauber | Parallel Programming-For Multicore and Cluster Systems[END_REF]. Next, a brief guide on how to judge a parallel performance is presented highlighting hardware and software limitations. Then, the parallelism in CIM-Lib is described restricting to solvers and dynamic load balancing (consuming the bigger part of the computing time). Finally, the efficiency of the parallel library is inspected via numerical applications followed by a conclusive discussion.

Concepts and Terminology

A first step into understanding parallel programming begins by introducing some concepts and terminologies. Even though, parallel algorithms play a key role to obtain a high performance computing application, it should be adapted to our computer/cluster architecture. Since many different possibilities of algorithms decomposition exist for the same application, choices should be made carefully. For computer architectures, several classifications exists. Herein, we present in sec. 3.2.1 the most popular taxonomy defined in [Flynn, 1972]. Flynn's classification scheme is based on the notion of a stream of information and instructions. Memory is another important factor affecting performance. Not only the memory size but also the speed in which the information is accessible impact To clarify the previous notations, one must mention the following points:

• A typical example of the SISD is the conventional sequential computer.

• Using the MISD represents a big challenge due to its nature executing different instructions on a single data memory. Thus it can be very restrictive and tricky.

• The SIMD and MIMD are the most used approaches. For instance, the SIMD is widely spread in multimedia/games applications especially to generate 3d prospects. While the MIMD is mostly used in multiprocessors systems or clusters. Between the SIMD and MIMD, it should be mentioned that the parallel algorithm for the former is considerably easier to implement.

Memory classification and access

As for the memory, the conventional classification consists in dividing it into two major types: a shared memory (Figure 3.2a) and a distributed memory (Figure 3.2b).

-For a shared memory system (called also global memory system), the memory is considered as a pool in which all the problem data are stored. All processors access the data via this global memory. In term of parallel computing, a shared memory machine is associated with threads based on a Compiler directives (OpenMP, TBB) -Distributed memory machines contain a number of processing elements called nodes and interconnection network. A node is formed by a processor and a local memory. Each processing element has its own local memory. The latter is private and can be accessed by one and only one processing element. The interconnections between processing elements -not sharing memory-are ensured via a network. In term of parallel computing, a distributed memory machine is associated with a process using a message passing interface to communicate between nodes ( MPI or Parallel Virtual Machine). Note that a third category can be found as well: a hybrid memory combining shared and distributed memory in one configuration (see Figure 3.2c).

Another important factor affecting a simulation computing time is the time needed to access the memory. In a shared memory machine, the access time to the data storage depends highly on CPU characteristics in a NUMA systems ( No Uniform Memory Access). Further details dealing with the memory management can be found in [START_REF] Hager | Introduction to High Performance Computing for Scientists and Engineers[END_REF]. Nowadays all processors has a cache memory. A cache is a small memory formed at most by three levels : L1, L2 and L3 (Figure 3.3). The L2 cache is the most important one. Actually, when requesting data, the processor will check L1, followed by L2 and L3. If the requested data is found in one of these levels, a cache hit is obtained, if not ( i.e. the data is found in the memory ) a cache miss is obtained. The more a cache hits faster the simulation will go. Thus the choices made can affect significantly the computing time.

Performance of a Parallel Code

Hardware performance

Processors are the most important components in a computer. They affect significantly the hardware performance. Since 2005, manufacturers like Intel and AMD deliver processors containing two or more cores (multicore processors). The processor performance is measured using several parameters. According to [www.top500.org, ], the most important ones are :

• #P roc: Number of processors (Cores)

• R peak : Theoretical peak performance provided by the manufacturer or the seller.

This performance indicates the number of instructions treated by all cores per second ( FLOPS ) .

• R max : The maximum value of R peak which cannot be exceeded by the cores. The real performance of the cores is given by R max . It is measured by the LINPACK Benchmark. The LINPACK Benchmark was introduced by Jack Dongarra. It consists to solve a dense system of Linear equations. A detailed description is available in [START_REF] Dongarray | The linpack benchmark: Past, present, and future[END_REF]. A parallel implementation of the Linpack benchmark can be found in [START_REF] Petitet | A portable implementation of the high-performance linpack benchmark for distributedmemory comput[END_REF] .

• N max : Problem size for achieving R max

Software performance

Despite all progress on recent multiprocessors and their performance, parallel algorithms remain essential to obtain an efficient parallel code. The performance is judged mainly by analyzing three conventional parameters: the execution time, the speed up and the efficiency. The execution time of any code is a primary concept in measuring the parallel implementation efficiency. Consider a size N problem running on p cores, the execution time adopts the following notation:

T p = t(N, p) (3.1)
The main criterion to evaluate a parallel algorithm is "Speed up". It refers to how much a parallel algorithm is faster than a sequential one; it is defined by the following formula:

S p = T 1 T p (3.2)
where p is the number of used cores, T 1 is the execution time of the sequential algorithm and T p is the execution time of the parallel algorithm using p cores.

Theoretically, the linear or ideal Speed up is obtained if S p = p, i.e. when executing an algorithm using p cores, the parallel computing time should be equal to the sequential time divided by p. Speed up could be classified into three categories: sub-linear, linear, super-linear (see Figure 3.5). Even though we aim for a linear speed up, the commonly reached is a sub-linear one. Nevertheless sometimes if the problem is small enough to be contained fully in the cache memory, a super linear speed up is obtained. In this case, the simulation runs the fastest.

The efficiency of a parallel algorithm is given by the ratio:

E = S p p (3.3)
Equation ( 3.3) will precise the scalability of the parallel algorithm: if the speed up is linear (S p = p), the efficiency is equal to 1 i.e. algorithm is as efficient running on 10 cores as on 1000 and does not require any further improvements. We should note that T p , the execution time of the parallel algorithm, plays a key-role impacting the algorithm efficiency.

Limits and Costs of Parallel Programming

Summarizing what has been discussed previously, the ideal speed up and efficiency verify:

S p = p and E = 1 (3.4)
In practice, it can not be easily reached. Many reasons affect the scalability of the code (such as the serial fraction, inter-cores communication and the load balancing) and will be presented next.

Serial fraction

In practice, several parts of the code cannot be parallelized and remain sequential. These algorithms are called the serial fraction of the code. [Amdahl, 1967] and [Gustafson, 1988] laws provide a theoretical approximation of speed up in codes containing a serial fraction.

Consider a problem of size N . The execution time needed to solve the problem using one core is noted by t(N, 1). The algorithm contains both sequential and parallel parts. Thus, the execution time is divided into two parts as well: t(N, 1) = t seq + t par where t seq and t par represent respectively the execution time of the sequential and parallel algorithms.

The serial fraction is defined as follows:

f seq = t seq t seq + t par = t seq t(N, 1) (3.5)

Amdahl's law

This law evaluates the computational time of a fixed size problem according to the number of cores p and the serial fraction f seq [Amdahl, 1967]:

t(N, p) = t seq + t par p (3.6)
Using (3.6) , Amdahl's law offers an approximation of the maximum speed-up that can be reached by the code:

S p(max) = t(N, 1) t(N, p) = t seq + t par t seq + t par p ≤ t seq + t par t seq = 1 f seq (3.7)
Equation (3.7) confirms that the presence of a serial fraction affects the scalability of a parallelized code. The speed-up is limited by a threshold inversely proportional to f seq . If one aims to improve the code scalability (i.e. speed up closer to p), the serial fraction must be decreased as much as possible (i.e. replacing sequential algorithms by parallel ones) .

Gustafson law

The Gustafson law consists on evaluating the speed up for a particular problems where the problem size varies linearly with the number of cores p (N is replaced by pN for p ≥ 1) [Gustafson, 1988]. The Amdahl law is replaced by the Gustafson law in which the the parallel execution time is taken as the reference : t(pN, p) = t seq + t par (3.8) If t par is the execution time on p processors (size problem equals pN ), the sequential processing time (size problem equals N ) is : t(N, 1) = t seq + pt par (3.9) Then the corresponding Speed up is : Gustafson's law offers a linear speed, a better estimation of the one proposed by Amdahl. The main difference is that the total execution time is now dependent of p.

With these kind of problems, the smartest step is to simulate the largest cases possible. In other word, if the data size increases the sequential part decreases relatively. According to equation (3.10) , the speed up increases with the number of cores and tends towards it (if f seq is small enough then S p → p).

Communications between cores

Inter-cores communication limitations are confronted when using a distributed memory machines: the communication between cores is required and the data exchange is done via an interconnection network. These operation lead to a loss of efficiency mainly related to:

1. Task synchronization: It is compulsory to synchronize the task management when the simulation runs using multiprocessors. If a processor needs information residing in the memory of another processor, it is necessary that the latter is available for communication. If the processor requesting information must wait, we are missusing these resources leading to loss of efficiency and computing time. Thus, it is necessary to have a good coordination between processors to insure the best synchronized communication possible.

2. Transfer time: Another important reason affecting the algorithm scalability is the time needed to transmit information from one node to another. The cost is assessed as follows: (3.11) where :

t comm = t latency + t pack + N v transf ert + (d -1)t transit
• t latency represents the time needed to establish a connection. It is equivalent to the time needed for sending a zero length message.

• t pack : represents the time needed to load the information into a buffer.

• v transf ert : represents the speed of data transfer in the interconnection network.

• t tansit : represents the time needed to transfer a message via a network if using distributed memory machines. The time of transfer, using a network of d internodes connections, is assessed equal to (d -1)t transit . For instance, in a shared memory machine, using a UMA, d is equal to 1 i.e. one direct connection. Thus, no transit time what so ever is needed.

Notice that the speed of data transfer depends highly on the hardware (used interconnection) as well as on the message transfer protocol. In order to minimize the communication time between nodes, it is recommended to send one long message rather than several short ones.

Load balancing

In scientific computing, load balancing algorithms are used to distribute the amount of work ( i.e. the workload ) as equally as possible on the number of cores in a parallel environment. These algorithms can be roughly divided into two different categories : static and dynamic ones. The static algorithms distribute the load evenly between cores at the beginning of the simulation. Such algorithms can affect the scalability since the current state of the simulation is not taken into account. Thus, a dynamic load balancing algorithm is recommended to fix this limitation. For instance, a finite element problem using a remeshing or mesh adaptation technique will present eventually a load imbalance. Their use increases locally the number of nodes on a particular core leading to an imbalance in distributing the workloads. An other example of imbalance workload occurs in contact search algorithms in Lagrangian approaches ( see chapter 1 for more details ). An introduction and discussion related to this topic can be found in [START_REF] Hendrickson | Dynamic load balancing in computational mechanics[END_REF]]. An overview on load balancing algorithms in parallel environment can be found in [Jimack, 1999].

For dynamic algorithms, a literature survey splits the existent methods into two major families. The first consists on partitioning the domain into sub-domains recursively using a geometric or topological criteria [START_REF] Williams | Performance of dynamic load balancing algorithms for unstructured mesh calculations[END_REF] and [ Van Driessche and Roose, 1995]. These methods are known as the direct methods. The others, called iterative, uses a diffusion algorithm. They start from an initial partition and exchange nodes between neighbor cores in order to improve the balance of the workload. The dynamic load balancing used in CIMLib, implemented by [Digonnet, 2001], will be described in sec. 3.4.2

Parallel computing in CIMLib

Though CIMLib is a fully parallelized library, we only highlight two important parts: the mechanical solver and the Dynamic load balancing. The first represents a large part of the simulation time. Whereas the second is considered as a crucial component to insure the work equilibrium since we use a mesh adaptation technique. Therefore their

S.P.M.D parallelization of the Navier-Stokes Solver

Let M = (N , T ) be the mesh discretizing Ω where N and T are respectively the set of nodes and simplexes forming the mesh. Note P the set of cores used for running the simulation. If M i is a partition of M then it can be expressed as follows:

M = card(P) i=1 M i (3.12) Every sub-mesh M i = (N i , T i
) is associated to a particular core p i ; where T i is the set of simplexes assigned to a core p i and N i is the set of nodes forming those simplexes. Note that a simplex can only belong to one and only one sub-domain. Nevertheless the set of nodes N i can be shared between several cores. The shared nodes form an interface between the different cores :

Γ ij = {M i ∩ M j /p i = p j } (3.13)
As illustrated in (Figure 3.6), the interface nodes are not replicated. It means that if two cores share the same nodes, they are stocked on one core and ghost nodes are created and stocked on the other one. Once the partition is created, each core p i reads the data related to the associated sub-mesh. Note that for now, no communications between cores is needed. Parallelizing the solver is summarized by two steps: The assembly of the stiffness matrix and the iterative resolution of the problem:

Ax = b (3.14)
where A is the global stiffness matrix, x and b are two global vectors (i.e. defined on the whole domain Ω).

To determine the global matrix A, each core p i must first assemble its own local matrix A i . We should mention that the dimension of each A i varies between cores depending on the number of nodes in the sub-mesh M i . The local rigidity matrix is assembled on each core. Then, the global stiffness matrix A is written as follows:

A = card(P) i=1 Âi (3.15)
where Âi is the extension of the local matrix A i . This extension is determined by resizing the matrix A i to fit the global dimension of A. It is accomplished by adding zeros for nodes belonging to only one core respecting the relations between the local and global positions of the variable nodes via a binding matrix ( see Figure 3 Now that the stiffness matrix A is determined, equation (3.14) can be solved. CIM-Lib uses the Portable Extension Toolkit for Scientific Computation ( PETSc ) [START_REF] Balay | Efficient management of parallelism in object oriented numerical software libraries[END_REF], [Balay et al., 2014a] and [START_REF] Balay | PETSc Web page[END_REF] for an iterative resolution method. Such methods are based on a matrix-vector product: Let us define a local vector rp representing the restriction of r on the core p. So the product matrix-vector is written as follows:

y = Ar =   p∈P Âp   r = p∈P (A p rp ) = p∈P (Ar) p (3.17)
With locally y p = (Ar) p = A p r p .

Nevertheless, the product y p does not represent the restriction of Ār p onto the core p. It represents the contribution of the latter product to the p -th core . Thus, the global vector y is given by: y = p ŷp (3.18) where ŷp represents the local vector resized as we mentioned before. So to use an iterative method we should define the restriction ȳp as : ȳp = y p + q∈Γpq y qp (3.19)

Dynamic load balancing

In CIMLib, the dynamic load balancing is treated as a renumbering of the mesh (nodes and simplexes ) under constraints [Digonnet, 2001]. In fact, the nodes/simplexes are numbered using a numbering couple (n 1 , n 2 ) :

• n 1 ∈ N is the number of the nodes/simplexes.

• n 2 ∈ N is the number of the domain p.

The above technique does not change the definition of the mesh, it simply introduce a new numbering of its nodes and simplexes. It is obvious that the numbering is not unique. Therefore a cost function is needed to compare the different numbering sets. The cost function assessing the computational time spent on a core p can be decomposed into:

φ = max p∈P t p + max p∈P c p (3.20)
where t p is the computational time required by the core p and c p the sum off all the communication related to the core p including the synchronization time. The previous cost function uses the maximum criterion. It generates difficulties when searching for the optimal numbering. Therefore the author introduces locally an order relation leading to optimize the cost function ( i.e. minimizing the cost function globally ).

In a parallel framework, the procedure above is based on an iterative algorithm coupling two processors at a time. As mentioned earlier, optimizing a partition can be achieved by optimizing the corresponding numbering. To manage this task, the author introduces three numbering operators R pp , R pq and I p . The operator R pp = (r n pp , r s pp ) acts locally (on a processor p) renumbering both the nodes and simplexes. The operator R pq = (r n pq , r s pq ) acts between processors p and q renumbering the nodes and simplexes. The operator I p is an identity operator leaving the numbering of a processor p intact. Repartitioning can be conceived as a succession of the renumbering operators R pp , R pq , I p . This renumbering strategy can be repeated as long as the re-meshing evolves in order to maintain both an optimal partitioning and cost function.

Applications

This section is dedicated to examine the performance of the monolithic Eulerian approach implemented in CIMLib. Through all this study, a 3d generic flat bunch simulation is used to inspect closely the parallel scalability. In fact, neither the type of the chosen application nor the complexity of the geometry should affect the obtained speed up. It should only increase the computational time of the simulation. This is true since the approach uses one mechanical solver all over the domain and each deformable body is tracked via the level set technique. In other words, the parallel scalability (i.e. speed up) is not affected regardless of the number of deformable bodies and their geometry (even if complex).

Case presentation

To begin with, we define the domain Ω = 40.5 mm × 33.5 mm × 51 mm where the deformable body of dimensions 27 mm × 20 mm × 15 mm is completely immersed in Ω (see Figure 3.8). The work domain Ω is meshed using one of the following meshes:

• a coarse isotropic mesh containing 65 017 nodes and 377 522 elements where the mesh size is equal to 1.25 mm;

• a refined isotropic mesh containing 491 055 nodes, 2 879 257 elements and a mesh size of 0.625 mm;

• an anisotropic mesh of approximately the same size of the coarse isotropic one.

The forging process is stopped when the forged piece reaches a deformation of 10% of its initial height. Thus, it is equivalent to 32 increments on the coarse mesh using a time step of 5.10 -3 s and 64 increments on the refined mesh using a time step of 25.10 -4 s. The scalability of the parallel code is tested versus several variables. Once the results are established, the refined mesh is used to highlight its effects if any on the scalability. The anisotropic mesh is used, at the end, to inspect the adaptation scalability used in CIMLib.

Navier-Stokes vs Stokes

We begin by a quick comparison between two different solvers in CIMLib : i) Navier-Stokes and ii) Stokes using the isotropic coarse mesh. The comparison is done using 2 n where n ∈ {0, 1, • • • , 5} cores. Table 3 In terms of resolution both solvers consume approximately the same amount of time even though they use completely different types of resolution : i) Generalized minimal residual algorithm ( GMRES) for Navier-Stokes and ii) Conjugate Gradient for the Stokes solver. As for the assembly time, Stokes solver takes approximately half the time needed 3.5. Applications 107 by Navier-Stokes solver. This behavior is expected since less terms are in play. For sake of convenience, Stokes solver will be used through all the study presented in this chapter.

Coarse mesh: MPI Version and Bind to Core Option

Using the 1.2.6 version of MPI and the coarse mesh, the simulations were launched on 2 n cores where n ∈ {0, 1, . . . , 5}. As mentioned in sec. 3.4.1, solvers are the ones consuming most of the computational time. For example, the Stokes solver represents 62.41% of the computational time when the simulation runs using one core. As for the convection solver, levellerT, it represents 15.5%. Thus, the time consumed by both solvers takes up to ≈ 78% of the computational time.

Note that the assembly time required by the LevellerT is very comparable to the mechanical solver assembly. Even though the local stiffness matrix to be assembled ( size 4×4 ) is smaller then Stokes ( size 16×16). Thus, we recommend optimizing the levellerT assembly to improve the global computational time.

For sake of clarity, we separate the assembly time from the resolution time. Table 3.3 and Table 3.4 summarize respectively the time used by the solvers ( Stokes solver and the convection solver -LevellerT-) and by other relevant algorithms of the simulation (generating output files, computing generalized strain, etc.). These numbers confirm that the same comment (for the sequential computing) is valid for parallel computing as well. To begin, the evolution of speed up and efficiency with the increasing number of cores are presented in Figure 3.9. Figures 3.9a and 3.9b show that, for the assembly part, a good scalability of 0.8 is obtained when using 32 cores. For the resolution part, a good scalability is obtained for a small problem such as the convection problem. However, we notice a loss of scalability for the Stokes resolution. Recent processors allow each core to communicate with its neighbors trying to optimize the attached process. Analyzing the results, we suspect that these communications can lead to an important loss of efficiency. Thus, we decided to migrate towards the 1.4.3 version of MPI. This newer version of MPI offers an option called "bind to core". This option is used to attach a process to a core throughout the whole simulation and hopefully guarantees a better performance. The same tests are relaunched but this time using the "bind to core" option. Again, the time consumed by the different parts of the code are presented in Table 3.5 and Table 3.6. The "bind to core" option does not change noticeably the final computational time. When examining the speed-up and efficiency presented in Figure 3.10, we notice that the efficiency of the Stokes resolution is improved from 0.43 to 0.58 (for 32 cores). Nevertheless, the loss of efficiency is mainly visible for less than 16 cores. We find that the latter is reasonable when eliminating the communication between cores (i.e. "bind to core" option) and when utilizing less than the full cores of a computing node.

To summarize, we can say that the 'bind to core' option has proven to be essential to insure a good scalability. The overall efficiency is satisfactory, it decreases linearly when using more than 4 cores. However we can notice a loss of efficiency when using less than 4 cores. This behavior is clear in Figure 3.9. This local loss can be attribute to the architecture of the multi-core processor. The next section is the key to prove this point. When iterating between different computing nodes (i.e using Round Robin option), we gain respectively 34 s and 15 s when using 2 and 4 cores (check Table 3.8). So we consider this loss (for p ≤ 4) as a direct outcome of the hardware limitation. Accordingly, these algorithms especially the resolution algorithm are scalable. If leaving aside the scalability of the output generating algorithm, we prove that the scalability of our code, is mostly limited by the scalability of the Stokes system.

The scalability loss of the output files generation is investigated next. The aim is to determine how the hardware can be manipulated to obtain the best results our machines can offer.

Coarse mesh: Hardware limitation and Round Robin Option

When analyzing the efficiency of our code in the previous sections, it was noticeable that the scalability of files generation algorithms was not satisfactory (Figure 3.9c). This result was expected since it is directly related to the used hardware. To understand better, we should explain that the used hard disc is formed by 3 platters and 3 heads which enable it to write in parallel. The speed up was obtained until reaching a specific number of cores (equal to 4 cores in Figure 3.9c) which means until all cores of a node are used. Afterward, the hard disc fails to manage properly writing in parallel.

This point is confirmed in Table 3.7: Though the size of the results files are equally divided on the number of used cores (i.e. each core writes the same amount of data), the hardware is limited and fails once more than 4 cores are used.

So the reasons behind the loss of scalability for output files generation is i) simulations are launched locally on nodes while results files are written in a parallel way (using vtk file format ".pvtu") and ii) the hard disc head's cannot manage data when different cores are in the process of writing output files. To try and extract the best of what this hardware can offer, the use of "Round Robin" option is recommended. This option can be handy to promote the scalability of the written output files.

The "Round Robin" option works in the following manner: it iterates on the available nodes attempting to distribute all processes equally on each node. This means that as soon as the code is launched on more than two cores, it uses two different nodes and accordingly two Hard Discs. Next, the same simulations are repeated using a "Round Robin" option.

The results presented in Table 3.8 and Table 3.9 confirms that nothing is affected but the output algorithms. In addition, Table 3.9 and Figure 3.11 show that both time and speed-up of the output algorithms are improved. By redistributing processes equally on each nodes, the "Round Robin" option increased the limit number of cores from 4 to 8 cores (i.e. for which the hard disc stop writing properly afterward). Still, after 8 cores the time used to write output files increases and the speed up drops again. To summarize, the output generation scalability is bound by the used hardware limitation. "Round Robin" option proved to be a great way to boost the scalability but still narrowed by the hard disc capacities.

Refined mesh

Using a refined mesh is expected to give better scalability. In fact, the workload of a core will increase. Thus, the ratio workload/communication using a refined mesh will increase as well leading to a better scalability. The same simulation is launched using up-to 128 cores. The results for time consuming are summarized in Table 3.10 and Table 3.11. The simulation time is considerably greater since more nodes are used in the new refined mesh. The overall behavior of the solver remains the same: the solvers consume the bigger part of the total computing time. 3.12 completes this discussion. The Stokes resolution remains the least scalable with an efficiency of 0.8. It is interesting to mention that after 32 cores LevellerT changes behavior and becomes less scalable. It is totally expected for small problems such as the convection problem (i.e LevellerT). Increasing the number of cores is efficient to a certain limit; afterward the scalability decreases due to an increasing communication between cores.

version of MPI with the bind to core option

To conclude, the study presented hereby proves the scalability of the code using up to 128 cores. The global behavior is very consistent with the results obtained on the coarse mesh. Nevertheless an undeniable improvement was detected: Increasing the number of nodes improves enormously the obtained speed up (check the speed up using 32 cores for instance). The behavior using refined mesh was proven to be satisfactory and consequently, only hardware limitations can be held accountable for any loss of efficiency.

Anisotropic mesh: Adaptation Scalability

Anisotropic mesh adaptation, described in chapter 2, is an important part of our approach. Its scalability highly matters in our framework. Its performance will be tested in this section in order to know if any improvements are needed.

Recall that the anisotropic mesh adaptation used is a h-adaptation. Thus at each adaptation step, three algorithms should be executed : i) the re-meshing algorithm where the mesh generator "MTC" generates the mesh according to the metric constructed at each node. ii) The transport algorithm where the solution and all essential fields are transfered from the previous mesh to the newly generated one in order to continue the computation. In other words, this step relies simply on interpolating the solution between the old mesh and the new one. iii) The dynamic load balancing algorithm to ensure a good workload distribution between different cores. Table 3.12 and Figure 3.13 show explicitly the results ( computation time, speed up, efficiency ) for both re-meshing and transport algorithms. The total adaptation time contains in addition the time spent executing the dynamic load balancing algorithm. Figure 3.13b shows that the scalability is steady over 4 cores. Even the transport caught up with the remeshing part when using 32 cores. The loss of scalability between 1 and 2 cores is due to the number of time the adaptation scheme is used depending of the computation environment (whether parallel or sequential). In a parallel execution the adaptation scheme is repeated two times to ensure a good quality of the mesh. 3.12: Times (in s) spent to accomplish the anisotropic mesh adaptation and its different parts. The loss of efficiency visible for the total adaptation time is due to the iterative method used in CIMLib. The algorithm starts with a given initial partitioned mesh. At first, the nodes shared between cores are fixed; and the re-meshing is done on nodes belonging to one and only one core. Then the dynamic load balancing algorithm is executed. This process is repeated until the couple mesh/solution has converged. This iterative scheme is illustrated in Figure 3.14. For more details, we encourage the readers to refer to [START_REF] Digonnet | Adaptation de maillage en parallèle, application à la simulation de la mise en forme des matériaux[END_REF], Digonnet et al., 2013]. Finally, It is important to mention that the time spent to adapt the mesh is equal to the time needed for Navier Stokes solver when adapting every 10 time steps (adaptation frequency = 10 ). In other words, the time is doubled. Despite this fact, we recommend using the adaptation technique since it improves drastically the obtained precision (see chapter 4).

Conclusion

In this chapter, we presented a general description of parallelism in modern computers. We focused essentially on the parallel environment in CIMLib. The main objective was to test the response of the library in a massively parallel environment to determine if any improvements are needed. In particular, important solver such as the Stokes and the convection solver LevellerT were inspected closely. The study has shown the need to optimize the latter solver in term of assembly. An MPI migration to newer version was recommended to insure better performance when exceeding 16 cores. As a result, the solvers performance was improved and found satisfactory using up to 128 cores. The efficiency of the mesh adaptation was tested as well. A loss of efficiency was noticed in the resulting dynamic load balancing. Though not crucial, we should note that it is mainly attributed to the used iterative partitioning method. Output generation exhibited a loss of efficiency also, it is for the most part associated the our hardware limitation and does not reflect programing deficiency.

To conclude and despite the casual loss, the overall scalability is found satisfactory. The performance of the major parts including the solvers and mesh adaptation was judged solid.

Advanced applications combining information and recommendations from both chapter 2 and chapter 3 are the focal point of the next chapter. The results are assessed and confronted with Forge ® simulations.

Résumé en Français

Dans ce chapitre, nous avons présenté une description générale du parallélisme dans les ordinateurs modernes. Nous nous sommes concentrés essentiellement sur l'environnement parallèle au CIMLib. L'objectif principal était de tester la réponse de la bibliothèque dans un environnement massivement parallèle pour déterminer si des améliorations sont nécessaires. En particulier, le solveur mécanique de Stokes et le solveur de convection LevellerT ont été inspectés de près. L'étude a montré que c'est recommandé d'optimiser le solveur de convection en termes d'assemblage. Une migration de MPI vers une nouvelle version a été recommandée aussi. Ainsi, meilleures performances sont assurées lorsque nous utilisons plus que 16 coeurs au cours d'une simulation. Par conséquent, la performance des solveurs a été améliorée et jugée satisfaisante en utilisant jusqu'à 128 coeurs. L'efficacité de l'adaptation de maillage a été testée aussi bien. Une perte d'efficacité est remarquée dans l'équilibrage de charge dynamique. Bien que non essentiel, il faut noter qu'il est principalement attribuable à la méthode de partitionnement itérative utilisée. La génération de sorties présentait une perte d'efficacité aussi. C'est associé à la limitation du matérielle utilisé et ne reflète pas de carence en programmation. Pour conclure et malgré la perte occasionnelle, la performance globale est jugée satisfaisante. Dans ce qui suit, des applications avancées combinant les recommandations des chapitre 2 et 3 sont le point focal du prochain chapitre. Les résultats sont évalués et confrontés à des simulations de Forge ®.

Introduction

This chapter is devoted to the different applications of the Eulerian monolithic approach described in chapter 2. Two sets of applications are presented. In the first, industrial applications such as the forging processes of rods and crankshafts are highlighted. The second is dedicated to the deformation of more simple geometries (such as spheres and Parallelepiped) but in a multi-domain context. In other words, the numerous deformable bodies (up to 250 bodies) will be stacked in the computational domains where there evolution will be followed during the simulation.

The aim of these cases will be testing both the capacity and the reliability of the method for complex geometries (in the industrial applications) in multi-domain situations.

Both the advantages and the limitations will be detailed and discussed. Consequent improvements will be proposed as well leading to a decisive conclusion if the approach is fit to these kind of applications. Whether casted, forged or machined the production of the connecting rod and crankshaft remains critical. Connecting rods should sustain complex loads such as the cycle of high compressive loads followed by high tensile ones; while cranks should be hardened and strengthen to endure bending and twisting during full functioning.

Industrial applications

For these reasons alone, such pieces are the topics of numerous researches including ours. In the aim of providing the industrial partners a simulation tool responding to their needs, different cases are explored next, using the monolithic Eulerian approach presented in chapter 2. This attempt is first of a kind since these problems are typically treated using fully Lagrangian approach

Connecting rod

In this section, the method capacity to model the deformation of a complex geometries is tested. For instance, the forging of a connecting rod is illustrated next. Figure 4.2 shows the position of the initial geometry ( orange ) between the two tools for t = 0s. Notice that, the details of the initial geometry of the deformable body and the tools are perfectly captured due to the anisotropic mesh. The used mesh is adapted on the gradient of both deformable body and tools level set. It contains 119 000 nodes and 687 784 elements. The simulation is launched using 16 cores. For a better visualization, a longitudinal section is represented in Figure 4.3. Both the contour of the deformable body ( orange ) and the tools ( blue ) are depicted and are perfectly smooth. The lower tools is fixed (v = 0 mm/s ), whereas the upper tools moves with a velocity equal to -10 mm/s along the Z axis. The tools consistency and the density are respectively equal to 1000 M P a.s and 8 g/cm 3 . The deformable body consistency and density are respectively equal to 352 M P a.s and 2.8 g/cm 3 . As for the air, the viscosity is equal to 1 × 10 -7 M P a.s and the density is equal to 1.2 g/cm 3 . 4.4 represents the zero iso-value of the deformable body for t = 0 s ( on the left ) and t = 2.6 s ( on the right ). We mention that, based on the same analogy presented in chapter 2, the air is allowed to evacuate by imposing a porous boundary conditions. The latter is equivalent to defining vents in the tools. Qualitatively, our simulation gives a good results. To validate quantitatively our approach, the same case is simulated using Forge ® . The Lagrangian mesh contains 37 483 nodes 42 438 elements. A comparison between the different results is presented in The chosen cut corresponds to the mid-plane of the deformable body. Results are very similar and represent the same features. The main difference, thought not significantly important, is mostly located in the center. For instance, the velocity along the x-axis represents a gap of 1.7 mm/s. This is directly related to the meshes illustrated in Figures 4.5a and 4.5b. In the CIMLib case, lesser nodes are used in the center leading to this difference. The same reason is behind the difference for the pressure and the velocity along the y-axis. Though, the Eulerian cut contains approximately half the nodes of the Lagrangian one ( 2900 nodes vs. 6000 nodes ), the results remains similar and satisfactory.

Using CIMLib

Using Forge ® To conclude, the Eulerian monolithic approach gives good results when compared with the Lagrangian approach. The slight differences shown in Figure 4.5 and Figure 4.6 are the outcome of the different used re-meshing techniques (an anisotropic adaptation technique in the Eulerian approach versus a re-meshing technique in the Lagrangian approach). Figure 4.5a shows that the nodes concentration in the Eulerian approach is located on the border of the deformable body. Whereas the re-meshing in the Lagrangian approach is based on the deformation and thus distributes the nodes appropriately in the discretized deformable body. Therefore the maps in the right part of Figure 4.5 are smoother then the left one.

(a) (b) (c) (d) (e) (f) (g) (h)

Crankshaft

The aim of the next case is to test the program performance for even more complex geometries such as the crankshaft. To push the evaluation even further, the same program used in sec. 4.2.1 for the connecting rod is used again. Only the tools and piece geometries are replaced. The same characteristics are used as well . The mesh remains intact. In other words, the same number of nodes and elements is used.

At t = 0 s, the deformable piece is placed between the two tools (illustrated in Figure 4.7). The upper tool moves with a velocity equal to v z = -80 mm/s while the lower tool remains fixed (v = 0 mm/s). Though the tools show more intricate details than the last ones, the mesh seems to captures every aspect. In Figure 4.8, the perfectly meshed upper tool confirms that the current mesh is versatile. At t ≈ 0.9 s, the isovalue of the deformable piece is illustrated in Figure 4.9. The crankshaft seems to take well the tools form. To complete the discussion, the velocities maps all over the deformable piece are presented in Figure 4.10. The velocities components v x and v y are bigger on the edges depicting the direction of the material flow. On the top, v z takes the upper tool velocity value 80 mm/s which is expected. The overall behavior seems satisfactory. Nevertheless, further inspection shows some lack of precision. Figure 4.11 plots the zero-isovalue of the tools and the deformable piece on a longitudinal section in the thickness. We have noticed an increasing penetration between the upper tool and the deformable piece. This suggests that even though the current mesh is sufficient to perfectly represent the initial geometries; it needs more nodes to capture the complex shape of the deformable piece during time. Notice that the mesh generator tried to salvage the situation by migrating nodes from the edge to the contact zone, but fell short. Figure 4.12 confirms this argument. The mesh was able to capture less complicated forms (such as the connecting rod) and thus no penetration was noticeable. In addition, the upper tool is translated analytically along the z-axis while the level set of the deformable body is convected using the velocity ( solution of Navier-Stokes ). In other words, the movement of the upper tool does not depend of the mesh since it is analytically determined at every time step. This can result in a clash of the bodies for lack of precision. This is why no penetration occurs with the lower tool. The latter remains fix while the deformable piece takes the corresponding form accordingly. This can be another reason behind this anomaly.

To reenforce our argument, we can compare the penetration of the deformed material into the two contact surfaces in Figure 4.11 ( the upper tool surface and the lower tool surface ). More penetration is visible with the upper tool. This is not surprising since the lower tool is fixed during the whole simulation, the mesh is more stable resulting in less (to no) penetration.

To summarize, the overall behavior is acceptable but the mesh size was proven to have a crucial impact on the precision maintain during the simulation progress.

Multi-domain applications

In the following section, the method capacity is tested when using a great number of deformable bodies at once. Two major applications are proposed. The first handles 10 parallelepiped deformed simultaneously using both anisotropic and uniform meshes. The second stacks up to 256 deformable spheres in the computational domain. In both, computations are treated in a highly parallel environment.

Ten Parallelepiped

A first test handling more than two deformable bodies at once is presented hereby. In a computational domain Ω = [0; 100] × [0; 100] × [0; 51] mm 3 , we immerse 10 identical parallelepiped of dimensions 5 mm × 5 mm × 10mm (see Figure 4.13). The deformable pieces are placed between two flat rigid tools. The lower tool remains fixed and the upper tool moves with a velocity equal to -10 mm/s. The deformable pieces viscosity is equal to 353 M P a.s while the rigid tools consistency is equal to 1000 M P a.s. Note that it is the same case mentioned in the General Introduction page 1 . We chose to use a uniform mesh. Since the computational domain is relatively big in comparison with the parallelepiped dimension, we created a new virtual box en-globing the deformable bodies but is still smaller than the domain. Though the mesh is uniform, we impose two different mesh sizes h int and h ext . Inside the virtual box, the size mesh h int is equal to 0.3 mm whereas outside the box h ext is equal to 2.5 mm (see Figure 4.14). This choice enable us to use less mesh nodes and elements. The overall mesh contains 2 282 109 nodes and 13 445 963 elements. The simulation is launched on 256 cores. We should indicate as well that 6 different level set functions were defined and associated manually to these 10 deformable pieces. problem what so ever. Nonetheless, the results were found unsatisfactory. The pieces did not deform properly: their top remained intact during the whole simulation. When examined closely, this irregularity is associated to the mesh size h int . Even though it is small, it remains significantly large when compared with the mixture law thickness ≈ 0.9 mm . If one decreases even more the mesh size, the problem will be prevented but it will be costly in computing time and the number of cores needed to run the simulation.

Although the results were inadequate for this particular case, we learned interesting information regarding the solvers performance. Considering the big number of unknowns (> 8 000 000) to be treated in this particular case, we were expecting an unrealistic number of iteration to insure the convergence. In fact, an average of 1 500 iterations were enough.

Figure 4.15: The deformed geometries when using an isotropic mesh with a mixture law thickness equal to 0.9 mm.

We were brought once again to use the anisotropic mesh generator. Leaving all the other components intact, only the mesh is changed. It is an anisotropic adaptive mesh formed by 147 500 nodes and 856 442 elements (illustrated in Figure 4.16). The simulation is now launched on only 16 cores. Remark that both the elements/nodes and the cores numbers were divided by 16. Not only less nodes/elements and cores were needed for the exact simulation, but also the results were much better. Figure 4.17 completes the discussion. The deformable geometries are acceptable. In addition, the velocity components and pressure predict a good behavior. v x and v y are bigger on the edges whereas v z is the highest on top. In addition, |v z | takes the value of the upper tool velocity on top (10 mm/s) and the lower tool on the bottom (0 mm/s). The pressure anticipates an adequate performance: it takes the maximal value on top where the upper tool compresses the most. Until now, we proved more than once that even though the overall performance of the approach is impressing, it remains extremely dependent of the mesh. It is an approach involving a big spectrum of parameters dependent one of the other (such as the mesh size, the thickness of the mixture law, ......).

In conclusion, important points should be retained from this section. In multidomain problems and up to 2 000 000 nodes and 256 cores, no suspicious behavior in the solver was detected. However, the use of an adapted anisotropic mesh were found crucial to insure the precision of the simulation.

256 spheres

To increase even more the number of the deformable bodies, we begin by simulating a problematic case proposed by Transvalor. In both Forge and CIMLib, 148 spheres are stacked to form a cubic domain partitioned on 8 cores. Before imposing any mechanical problem, the partition is inspected closely. In Forge ® , the partitioning algorithm is unable to converge and the tasks are not equally divided between the 8 cores. Figure 4.18 shows that the workload is almost handled by only one core which affects the scalability of the parallel computation (see chapter 3). Whereas, CIMLib confronted no difficulties partitioning the computational domain on 8 cores. As shown in Figure 4.20 , The workloads seem to be almost equally divided between the different cores. Note that the lack of smoothness is only due to the unrefined mesh used for this test. As mentioned earlier, the dynamic load algorithm does not allow to exchange nodes between unneighbored domains ( i.e. does not share a node ).

In order to generalize the algorithm, an exception should be added for the initial partitioning when dealing with a sequential mesh. Thus, the exchange between two domains -not sharing an interface-is allowed only for the first step in the dynamic load algorithm. A node is chosen and is affected to an empty core then the classical algorithm is executed.

However, when considering a mesh formed by multiple connected geometries, the classical algorithm fails to partition properly the domain. A random node is chosen to begin with. Once the connected geometry containing this node is handled by the algorithm, it is unable to handle another one. This is logic since the core is no longer empty and the domain has no neighbors. Thus, the dynamic load algorithm stops even though the imbalance remains.

To overcome this limitation, a simple yet effective modification is introduced. the core is allowed to handle new nodes when it has no neighbor and even if it is not empty. The latter modification allows to improve the imbalance shown in Figure 4.18 (see Figure 4.19 ) Note that this is the best partitioning we can offer considering the fact that no geometrical notions are present in the algorithm. Since CIMLib had no problems what so ever partitioning a domain containing 148 spheres on 8 cores, why not push this study further to test the limits of this approach.

Multi-domain applications 139

To do so, 256 identical spheres of a radius equal to 0.1 m are placed in a computational domain of dimension [0; 4]×[0; 4]×[0; 1.63] m 3 to form a cube. Then, the spheres are placed between two rigid flat tools with a consistency of 10 6 P a.s. The lower tool is fixed while the upper one moves with a velocity equal to 1 m/s. The whole computational domain is anisotropically meshed using 252 204 nodes and 1 440 718 elements. Since the deformable spheres will interact at one stage of the simulation, we cannot associate one level set for all the 256 spheres. Instead we should associate one level set for each deformable body. This can be hard to accomplish due to the high number of deformable bodies. To pull-off this step, we utilize an existent technique based on a coloration algorithm developed by [Hitti, 2011] . This coloring algorithm enable us to associate one level set to one or more deformable spheres with only one restriction. Two different spheres can share the same level set if and only if they never interact during the simulation. In other words, the same level set cannot be associated to neighbor spheres. Applying this technique, we were able to define only 40 level set functions to all 256 spheres. The spheres consistencies are chosen arbitrary between 10 5 P a.s and 7 × 10 5 P a.s. The simulation is launched on 64 cores. Note that we could of used less level set functions. In-spite of this fact, we chose to define 40 functions to test the limits of the program . We expected to confront a memory clash but we did not. For a great number of level sets, the resulting data will be normally overwhelming and could lead to a lack of memory space. In this case, we recommend to compress the unused data to increase the remaining space. Once the data is needed, it can be uncompressed. This can be a solution when significant number of deformable pieces is needed in future works.

Finally, the deformed domain and the different velocity components are illustrated in Figure 4.22 at t = 0.64 s. We notice that every sphere deforms differently inside the domain (due to their different characteristics) and interacts with neighbor spheres to fill properly the gaps. v x and v z are represented on a longitudinal section in the thickness, while v y is represented on a longitudinal section in the width. These maps portray a good behavior. The Eulerian monolithic approach was proven more resilient than Forge ® when undertaking multi-domain problems containing significant number of deformable bodies. In particular, the partitioning algorithm was found more efficient. In addition, the results issued by the mechanical problem give a good description of the deformation evolution with time.

Conclusions and Discussions

This chapter is entirely dedicated to numerical applications using the monolithic Eulerian approach along with the proposed improvements (described in chapter 2). All 4.4. Conclusions and Discussions 141 applications were carried out in a scalable parallel environment on at least 16 cores. The efficiency of the method was studied in more than one context using complex geometries (such as crankshaft and connecting rod) or simple deformable geometries stacked in great numbers (up to 256 spheres).

The results were judged via comparison with simulations launched in Forge and found satisfactory. The overall performance of this approach was solid. It was even considered more resilient than Forge ® in some cases; in particular when undertaking multidomain problems containing significant number of deformable bodies.

Although the performance is generally impressive, some points can be improved even more. The approach itself involve a big spectrum of parameters dependent one of the other that need to be mastered (such as the mesh size, the thickness of the mixture law, ......). In more than one occasion, we proved that the method precision remains extremely dependent of the mesh. When dealing with evolving geometries exhibiting details more visible with time, the mesh felt short and problems such penetration were detected. In addition, when comparing with Forge, we noticed that our results exhibited less precision in the center than on the edge. It is mainly because more nodes are distributed on the edge respecting the level set gradient.

Following these observations several upgrades/enhancement can be proposed:

• An increasing number of nodes can be used as the simulation progress. It can be determine using a relationship between the geometry surface and the mesh size. This can easily solve the problem of penetration when using complex deformable geometries.

• The mesh adaptation can be based not only on the level set gradient but on additional fields such as the velocity. This can improve the precision in the center (for instance) to match the outer precision depending of the application type.

To this stage, this approach was proven effective. Furthermore, the superposition with Forge ® was reassuring: two different softwares based on different approaches provided not only promising results but extremely comparable ones as well. The fact that a recent approach like the Eulerian one can compete with Forge -a well known software perfected for at least 20 years-drives us to enrich the method and search for new answers. Hence, new ways to model friction and contact in an Eulerian context are proposed and tested in the next chapter.

Résumé en Français

Ce chapitre est entièrement dédié aux applications numériques utilisant l'approche Eulérienne monolithique ainsi que les améliorations proposées dans le chapitre 2. Toutes les applications ont été réalisées dans un environnement parallèle utilisant au moins 16 coeurs. L'efficacité de la méthode a été testée dans plusieurs contextes en utilisant des géométries complexes (telles que le vilebrequin et la bielle) ou des géométries simples empilées en grand nombre (jusqu'à 256 sphères par exemple). Les résultats ont été jugés par comparaison avec des simulations lancées en Forge et trouvés satisfaisants. La performance globale de cette approche est solide. Elle est même considérée plus pertinente que Forge® dans certains cas; en particulier dans les cas multi-domaines avec un grand nombre de corps déformables. Bien que la performance soit généralement impressionnante, certains points peuvent être améliorés. L'approche elle-même implique un grand spectre de paramètres dépendants l'un de l'autre, qui doivent être maîtrisés (comme la taille de la maille, l'épaisseur de la loi de mélange, ......). A plusieurs reprises, nous avons prouvé que la précision de la méthode reste extrêmement dépendante de la taille de maille. En traitant l'évolution d'une géométrie présentant de plus en plus de détails au cours du temps, le maillage ne semble pas entièrement satisfaisant et entraine des problèmes de pénétration par exemple. De plus, suite à des comparaisons avec Forge pour ces caslà, nous avons remarqué que nos résultats présentent moins de précision au centre que sur les bords. C'est principalement dû à la distribution des noeuds plutôt localisés au bord de la pièce déformable. Suite à ces observations plusieurs améliorations peuvent être proposées: • Un nombre croissant de noeuds peut être utilisé vis-à-vis l'avancement de la simulation. Ça peut être déterminé en utilisant une relation entre la surface de la géométrie et le maillage. Cela peut facilement résoudre le problème de pénétration entre des pièces déformables complexes et les outils. • L'adaptation du maillage peut se baser non seulement sur le gradient de la Level Set mais sur des champs supplémentaires tels que la vitesse. Cela certainement améliorera la précision des solutions au centre de la pièce.

Pour le moment, cette approche a été prouvée efficace. En outre, la superposition avec Forge® est rassurante: deux logiciels différents basés sur différentes approches ne fournissent pas seulement des résultats prometteurs, mais très comparables aussi. Le fait que l'approche Eulérienne, une approche relativement nouvelle pour telles applications, puisse rivaliser avec Forge nous motive à enrichir la méthode même plus. Par conséquent, de nouvelles façons à modéliser le frottement et le contact dans un contexte Eulérien sont proposées et testées dans le chapitre suivant.

Introduction

Friction and contact are two important features in every large deformation problem. As pointed several times through out our study, Eulerian approaches are rarely used in this context. Thus, a lack of information handling this topic is detected in the literature. This chapter offers first steps towards friction and contact modeling in an Eulerian environment.

The chapter is mainly divided into two parts. The first part is dedicated to friction modeling. We propose different manners to numerically imitate friction depending on the simulation. A linear mixture law emulates a perfectly sticking friction, whereas a quadratic mixture law introduces, via a boundary layer, a lubricant to simulate sliding friction. The latter requires additional manipulation to determine the appropriate lubricant consistency for each case. An identification technique superposing both Lagrangian and Eulerian formulations is used. For each friction coefficient in the Lagrangian approach corresponds a lubricant consistency spread along a certain thickness. Followed the technique description, 3D test cases are launched and results are validated with Forge ® simulations.

Recalling the literature survey in chapter 1, the contact detection algorithm affect enormously the parallel scalability of the code. In Eulerian approaches, contact detection is managed automatically. Still, no way to model different types of contacts is available. To our knowledge, there are very few works in the literature dealing with this topic. The only work found is proposed by [START_REF] Bruchon | Using a signed distance function for the simulation of metal forming processes: Formulation of the contact condition and mesh adaptation. from a lagrangian approach to an eulerian approach[END_REF]. In the second part of this chapter, a daring proposition is presented to model the contact problem. Searching for inspiration in other domains, we decided to apply the same reasoning used in geology and oceanography. A directional solver is a great candidate and may hold the answers to our needs. The formulation of the anisotropic solver is detailed then followed by different applications. Its capacity is assessed ending in a small discussion and conclusions.

Friction: from Lagrangian to Eulerian description

It is important to model friction in large deformation problems. For instance in forging process, friction leads to heat and tools wear. Thus, it can affect the final shape and quality of the product.

In reality, to control the effects of friction, lubrication techniques are used. Generally in an Eulerian context, modeling surface forces ( such as friction ) consists on crossing from surface to volume integral using Dirac function [START_REF] Bruchon | Using a signed distance function for the simulation of metal forming processes: Formulation of the contact condition and mesh adaptation. from a lagrangian approach to an eulerian approach[END_REF] (see chapter 1).

In this work, a simpler approach based on boundary layer is used. A new viscosity corresponding to the lubricant is introduced. To determine this appropriate viscosity, a similar approach to [Foudrinier, 2007] and [Mondalek, 2012] is applied hereby. The first worked on reactive molding processes. Whereas the work presented in [Mondalek, 2012] treats modeling of spark plasma sintering process.

Based on a Lagrangian description of friction, this approach relies on finding a new formula describing friction in the Eulerian context. The detailed methodology is presented next.

Lagrangian Friction law

For every multi-domain problem, the physical phenomena in the contact zone is modeled using friction laws. Friction can be described as the force acting tangentially along surfaces in contact and resisting the motion of solid surfaces.

For instance, in a case containing a tool and a deformable body in contact. Let us denote respectively v t and v the velocities of the tool and the deformable body. The sliding velocity v s is computed in terms of v and v tool :

v s = ( v -v tool ) -[( v -v tool ) . n] n (5.1)
where n is the normal to the surface.

In addition, the shear stress is given by:

τ = σ n -(σ n. n) n (5.2)
In literature, numerous behavior laws are introduced. These laws describe the relation between the sliding velocity v s and τ and are divided into two families. The first family is based on a threshold concept such as in Tresca and Coulomb friction laws (Figure 5.1a); where below this threshold no sliding can occurs (equation (5.3)). For Coulomb law, the threshold τ c is proportional to the normal stress σ n .

   τ < τ c no slip τ = τ c slip (5.
3)

The second family includes regularized law such as Norton-Hoff law. In our study, Norton-Hoff friction law is used due to its continuity and regularized form (Figure 5.1b). Using the latter law, the relation between the shear stress τ and the sliding velocity v s is expressed by:

τ = -α f η.|v s | p-1 v s (5.4)
where η is the material consistency and α f the friction coefficient ( chosen depending on the material ). p can be considered equal to m the strain rate sensitivity coefficient since the lubricant is not considered as a third body.

Eulerian friction law

To apply the boundary layer approach, we begin by explaining its concept. We consider the existence of a boundary layer located between the two bodies in contact (tool/deformable body). This new layer represents the lubricant with its own viscosity (see Figure 5.2).

Unlike other works [Mondalek, 2012], where the boundary layer is explicitly represented, this viscosity once determined helps representing implicitly the boundary layer via quadratic mixture law presented in chapter 2. To proceed in determining the layer viscosity, we must first express the shear stress in the Eulerian context. Then by identification, the viscosity is determined from the Lagrangian law (equation 5.4).

Friction: from Lagrangian to Eulerian description 149

In the boundary layer, the shear rate is written as :

γ = | v s | ε (5.5)
where v s is the tangential velocity defined in equation ( 5.1) and ε is the thickness of the boundary layer (see Figure 5.2). Evaluating the strain rate tensor in the latter layer, it can expressed as :

ε(v) =       0 v s 2ε 0 v s 2ε 0 0 0 0 0       (5.6)
Thus the equivalent strain rate tensor is computed in the following manner :

ε(v) = 2 3 ε(v) : ε(v) = | v s | ε √ 3 (5.7)
Using the definition of the shear stress given in equation ( 5.2) and substituting ε and ε by their new form (5.6) and (5.7) computed in the boundary layer, the shear stress is expressed as :

τ = -2η lub .   √ 3 v s ε √ 3   m-1 v s 2ε = - η lub ε m .| v s | m-1 v s (5.8)
where η lub represents the lubricant viscosity.

Both equations (5.4) and (5.8), describe the form of the shear stress. One in a Lagrangian formulation and the other in an Eulerian one. These formulations must be equivalent and it is ensured via an identification technique :

η lub = α f ηε m
(5.9)

The layer viscosity η lub depend highly on the chosen thickness and the friction coefficient. η lub will be finally transferred into equation (2.76) (see chapter 2) to represent this heterogeneity in the monolithic approach.

Friction validation

In large deformation problems, especially the ones involving multiple materials in contact, friction is considered an important component. It interferes in the deformation evolution and affects the final shape of the piece (in this case the ring). Therefore this section illustrates the numerical results of friction in an Eulerian formulation as presented in the previous description. In the computational domain [0; 75] × [-75; 75] × [0; 55] mm 3 , we immersed a perfectly symmetrical ring illustrated in Figure 5.3 between two rigid tools. The inner and outer radius denoted respectively r int and r ext are equal to 40 mm and 50 mm. The half-height h is equal to 40 mm. Since the ring is symmetrical, only one fourth of the piece is deformed to decrease the computational time. Imposing a symmetry is equivalent to verifying the condition σ. n = 0. The mesh of a quarter ring contains 39 000 nodes. Note that the lower tool is fixed while the upper tool moves with a velocity of -5 mm/s. The rigid tools consistency η r and density ρ r are respectively equal to 1 000 M P a.s and 8 000 kg/m 3 . As for the deformable body, η d and ρ d are equal to 66 M P a.s and 2 800 kg/m 3 . The air viscosity and density are respectively fixed to 10 -7 M P a.s and 1.2 kg/m 3 . Different scenarios demonstrating how to represent friction in this Eulerian description are considered. We start with a perfectly sticking friction and finish by a sliding friction. The results are superposed with Forge ® simulations. Both advantages and limitations are highlighted before ending with a small discussion.

To begin with, we impose a perfectly sticking friction. In our approach this can be accomplished by using a linear mixture law illustrating the different domains. In fact, since the ratio η r /η d exceeds 100, the interface consistency is closer to the rigid body consistency. Therefore, it can be considered as an excessive numerical friction or "sticking friction". The latter is illustrated on a cutting plane in the thickness in Figure 5.5a.

The evolution of the ring in time is followed. Figure 5.4 illustrates the deformed ring at t = 3.467 s. The half-height of the deformable ring is equal to 22.7 mm. The inner and outer radius for h = 0 are approximately equal to 54.89 mm and 66.34 mm. The velocities and pressure maps are presented as well for t = 3.467 s in Figure 5.5. The overall behavior seems satisfactory. The velocity component v z is equal to -5 mm/s on the top conform with the upper tools velocity. The velocity components v x and v y take a maximum of approximately 4.77 mm/s on the lower edge, which suggests that the deformation is symmetric as expected. To validate the results, the same simulation is launched in Forge taking a sticking friction. Since the deformable geometry is hollow inside, we did not plot the curve on a vertical section. Instead, we took a cutting plane in the thickness for y = 0 in both interfaces (Paraview for CIMLib and Glview for Forge ® ). Important information can be extracted. First, by analyzing Figures 5.7a and 5.7b we notice that the velocities v x plotted on the outer outline in both softwares are globally very similar. Still, decreasing the mixture law thickness from ε = 0.5 mm to ε = 0.25 mm affects locally our results to fit better the Forge simulation on the top (z = h). Since this contact region is of a great importance, the parameter ε should be mastered to find the optimal results. We should mention that decreasing the mixture law thickness is not that simple. Smaller the chosen value is, more difficulties will be confronted to mesh this region. In other words, more stretched elements will be generated and smaller time steps will be needed. The end game is to find a thickness value ε small enough to increase the precision but big enough to prevent meshing difficulties. We should mention as well that small thickness values ε can not represent the boundary layer when treating sliding friction. This topic will be treated in details in the following cases.

Using the smaller mixture law thickness ε = 0.25 mm, v z and v x are respectively plotted on the outer and inner outlines. Once again the curves are similar and the average error does not exceed 2 %. This slight difference can be attributed to the different used meshes as well as the different interpolations used by Glview and Paraview. The second part of this study is dedicated to sliding friction. As mentioned earlier, "sliding friction" is equivalent to using the quadratic mixture law introduced in equation (2.76) (chapter 2). In other words an additional phase is added. It is called the boundary layer and it represents the presence of a lubricant with η lub (Eulerian formulation). The consistency of this lubricant is chosen to verify equation (5.9). To be more precise, the couple (η lub , ε) should be chosen to give the same effects as the coefficient α f in the Norton-Hoff friction Law (Lagrangian formulation). For instance, if we need to simulate a case with sticking friction using α f = 0.3 in a Lagrangian approach, for a fixed mixture law thickness ε = 0.25 mm we should use η lub = 4.95 M P a.s in our monolithic Eulerian approach. This means that the couple (η lub , ε) = (4.95, 0.25) in CIMLib should give the same results as when using α f = 0.3 in Forge ® .

The curves are plotted in the same manners as before. The results are compared in Figure 5.8. Mainly, the curves show good agreement in both softwares. The overall average error does not exceed 4%. To assess even more the capacity of this Eulerian method to impose sliding friction, we attempt to model the same case by decreasing drastically α f to 0.03. For the same thickness ε = 0.25 mm, the corresponding η lub is now equal to 0.495 M P a.s. On the same cutting plane, the results were inspected. In Figure 5.9(a), the superposition of curves plotted in both softwares reveals this time a relatively big gap between the results. They seem almost parallel one to another. The velocity v x obtained in CIMLib is almost exactly the same found for sticking friction. When examining in details the cause of this difference, we noticed that for this particular α f , η lub =0.495 M P a.s is very small in comparison with η d = 66 M P a.s. And since we are using a small thickness mixture equal to 0.25 mm, the exact value of η lub was not reached in the boundary layer. Actually, the latter is almost non-existent as illustrated in Figure 5.9(b). This can be solved by increasing the mixture thickness to ε = 0.5 mm for instance. The curves are compared again in Figure 5.9(c). The curves fit almost perfectly except on the edge where an odd hill is visible in the vicinity of z = h. Figure 5.9(d) presents the reason behind this bump: even though the boundary layer is perfectly represented, a miss-functioning is detected on the upper corner. Since η lub =0.99 M P a.s (corresponding to the new mixture thickness 0.5 mm) is very small in comparison with η d = 66 M P a.s. The resulting viscosity is much smaller than η d as well and is translated in bigger displacement on the top (for z = h). Following the last analysis, a major limitation was noticed. The proposed method gives altogether great results for the chosen α f values . Nevertheless, it cannot be generalized for the whole range of α f (respectively (η lub , ε)) especially for smaller values. The mixture thickness ε was proven to be a key-parameter affecting significantly the results.

It should be big enough to represent properly the boundary layer but remains insufficient for very small values of α f .

To conclude, we can say that the boundary layer should be a good candidate to simulate fiction in an Eulerian framework. However, a proper rpresentation of the boundary layer was proven to be a delicate matter as shown in Figure 5.9.

Next, we present a directional solver attempting to model bilateral sliding friction.

Directional Solver

To model the contact problem, we propose in this section developing an anisotropic solver. A long reasoning was behind this choice. Two types of contact are usually available: bilateral and unilateral contact. The bilateral contact can be expressed as an extension of the formulation proposed in the previous section. We need to impose a linear mixture law in the normal direction and a quadratic mixture law in the tangential direction. The first insures a sticking contact in the normal direction whereas the second allows sliding in the tangential direction. If we successfully represent a directional tool, the same results obtained in Figure 5.9 are expected but with a slight improvement in the critical zones ( in the vicinity of the corners ). This is the main focus of this section. We should note that an anisotropic solver can be exploited in many ways.

After breaking these facts, developing an anisotropic solver seems to be the right candidate to our approach. We began searching in the literature for similar works or at least for any directional notion that we can exploit to our favor. Two particular methods stand out. The first category is used frequently in oceanography and geology modeling [START_REF] Bellec | Sur des modèles asymptotiques en océanographie[END_REF]][Iftimie and Planas, 2006][Zhang and Fang, 2008]. They generally proposes to represent the consistencies anisotropy using a diagonal tensor. In other words, the scalar consistency is replaced by a tensor, where the different consistencies are directionally dependent. They frequently suppose that the consistency in the z-direction is a small value tending toward zero:

η =     η x 0 0 0 η y 0 0 0 η ε    
(5.10) where η x , η y and η ε are respectively the consistencies in x, y, and z directions with η ε → 0.

Using this tensor, the Navier-Stokes equations are rewritten: The second category based on the fiber theory proposes a more general form to portray the anisotropic viscosity by introducing an orientation tensor. Based on the work of [START_REF] Folgar | Orientation behavior of fibers in concentrated suspensions[END_REF], [Redjeb, 2007] proposes a new decomposition using the strain rate tensor ε(v). As a result a new behavior law taken into account the directional viscosity is presented.

∂v ∂t -η x ∂ 2 ∂x 2 + η y ∂ 2 ∂y 2 + ε ∂ 2 ∂z 2 v + v.∇v = -∇p div v = 0 (
Note that several work were dedicated to prove the existence and unicity of the solution in the anisotropic equations [Paicu, 2005] [Iftimie and[START_REF] Iftimie | [END_REF]. We highlight in particular, the work presented in [START_REF] Besson | Un modèle asymptotique en océanographie[END_REF] proving the existence and unicity of the anisotropic solution in the same Sobolev spaces used in our standard formulation (described in sec. 2.5.4.1).

In this work, we rather chose the second approach to express the anisotropic consistency. Inspired by their decomposition, the stain rate tensor ε(v) is multiplied by an orientation tensor a to express the anisotropy. The orientation tensor will determine mainly in which direction the deformation/velocity is greater.

The deviatoric stress tensor can be expressed as follows:

s = [ε(v).a + a.ε(v)]
(5.12) in which a is the orientation tensor depending of the anisotropic consistency η.

If η =     η x 0 0 0 η y 0 0 0 η z   
 denotes the diagonal consistency tensor with η x , η y and η z respectively the consistencies in the x, y and z directions, the orientation tensor a can take the following form:

a = t RηR (5.13)
where R is a rotation matrix .

Equation ( 5.13) guaranties that the couple (η, R) offers a general description of the anisotropy (i.e. in all directions).

Combining equations (5.12) and (5.13), the deviatoric tensor is now

s = ε(v). t RηR + t RηR.ε(v) (5.14)
Since R is antisymmetric and η is diagonal then the orientation tensor a is symmetric ( t a = t ( t RηR) = t RηR = a). By that, we insure that the deviatoric part of the Cauchy tensor remains symmetric as well.

Note that if R is equal to the identity tensor I and all the directional consistencies are equal (i.e. η x = η y = η z ), we recover the classic form of the deviatoric tensor s = 2ηε(v).

To conclude, the system to be solved is summarized by: (5.15) in which the inertia term are dropped for sake of simplicity.

       ∇.σ = 0 ∇.v = 0 σ = [ε(v).a + a.ε(v)] -pI

Weak formulation

To retrieve the weak formulation of (5.15), the same steps as in chapter 2 are used. The constitutive equations are multiplied by the test function couple (w, q). The only noticeable difference is replacing s by ε(v).a + a.ε(v). This substitution leads to new terms in the variational formulation. The integral (2ηε(v) : ε(w)) Ω i will be replaced by the following one (ε(v).a : ε(w)

) Ω i + (a.ε(v) : ε(w)) Ω i :    (ε(v).a : ε(w)) Ω i + (a.ε(v) : ε(w)) Ω i -(p, ∇.w) Ω i = 0 (∇.v, q) Ω i = 0 (5.16)

Spatial discretization

In this section, we recall some basics about the spatial discretization of the incompressible Stokes equations by the finite element method. The emphasis will be on the new terms introduced in system (5.15). Then we describe the implementation of our code in CIMLib. In order to simplify the problem into a practical for computation, we will be working locally on each simplex K.

Let (N i ) i=1,...,D denote the nodal shape functions. Then the solution v and the test function are written as follows :

v K = D i=1 v i N i and w K = D j=1 w j N j
As mentioned earlier, the only modifications will concern the local matrix denoted A loc vv (see equation (2.68) in chapter 2 ). It is the matrix form of the following term found in the classic solver :

K ε N i k : ε N j l dK =                    1 2 ˆK ∂N i ∂x l ∂N j ∂x k dK if k = l 1 2 K d m ∂N i ∂x m ∂N j ∂x m dK + 1 2 K ∂N i ∂x l ∂N j ∂x l dK if k = l (5.17)
where N i k is the k th component of the shape function on the i th node. However computing the new A loc vv in system (5.16) proves to be more expensive than the classical one. For instance,

ε i k .a =                          ∂N i ∂x 1 a k1 • • • ∂N i ∂x 1 a kk • • • ∂N i ∂x 1 a kd ∂N i ∂x 2 a k1 • • • ∂N i ∂x 2 a kk • • • ∂N i ∂x 2 a kd . . . . . . . . . . . . . . . d m ∂N i ∂x m a m1 + ∂N i ∂x k a k1 • • • d m ∂N i ∂x m a m1 + ∂N i ∂x k a kk • • • d m ∂N i ∂x m a m1 + ∂N i ∂x k a kd . . . . . . . . . . . . . . . ∂N i ∂x d a k1 • • • ∂N i ∂x d a kk • • • ∂N i ∂x d a kd                         
(5.18)

Instead, we write (ε i k .a + a.ε i k ) : ε i k using the gradient decomposition. The new A loc vv is now computed using the following form:

A loc vv = 1 4 ˆK ∇N i k .a + t ∇N i k .a + a.∇N i k + a. t ∇N i k : ∇N j l + t ∇N j l dK (5.19)
Equation (5.19) gives eight contraction terms leading to eight local matrices. This increases the computational time. To overcome this increase, we compute only four matrices defined by equations (5.20) to (5.23) . The other four are deducted by computing their transpose (5.24) :

A loc 1 = ˆK ∇N i k .a : ∇N j l dK ( 5
A loc 5 = t A loc 1 , A loc 6 = t A loc 2 , A loc 7 = t A loc 3 , A loc 8 = t A loc 4
(5.24)

The local matrix A loc vv is now the sum of all the individual matrices defined above:

A loc vv = 1 4 A loc 1 + A loc 2 + • • • + A loc 7 + A loc 8 (5.25)
Term by term decomposition:

∇N i k .a n,m =        0 if n = k m p=1 ∂N i ∂x p a pm if n = k ∀m = 1, . . . , d (5.26) t ∇N i k .a n,m =                 ∂N i ∂x 1 a k1 • • • • • • • • • ∂N i ∂x
a k1 • • • • • • • • • ∂N i ∂x d a kd                 (5.27) ∇N i k .a : ∇N j l =        0 if k = l m n ∂N i ∂x n a nm ∂N j ∂x m if k = l (5.28) ∇N i k .a : t ∇N j l = p ∂N i ∂x p a pl ∂N j ∂x k (5.29) t ∇N i k .a : ∇N j l = d m ∂N i ∂x l a km ∂N j ∂x m (5.30) t ∇N i k .a : t ∇N j l = d n ∂N i ∂x n a kl ∂N j ∂x n (5.31)

Variational Multi-scale Method

Similar to the classic Navier-Stokes Solver, the anisotropic Stokes solver is unstable as well. A stabilization technique should be put in place. A decision had to be made between the P1+/P1 and the VMS technique. An equivalence between these two methods was established in [START_REF] Brezzi | b = ´g[END_REF]. The choice fell on the latter. It is more general (it can be used in many cases) and offers a simpler yet efficient stabilization to our problem. The stabilization is insured via one parameter called τ c . In other words, the condensation needed in P1+/P1 is avoided and instead automatically managed using τ c . In what follows, an extension of the VMS stabilization technique presented in chapter 2 is adapted to our system of equation (5.12). The upside of this choice is our ability to use the same VMS technique by only replacing the stabilization parameter with the new one.

Stabilization techniques are used for numerous problems nowadays, the literature offers many eligible stabilization parameters [Tezduyar and Osawa, 2000a]. We mainly were interested in the work presented in [Principe, 2008] where a new stabilization parameter taking into account the anisotropy is proposed.

The same reasoning is used as before. The solution spaces V and Q are respectively decomposed into two spaces V h and Q h (resolvable space) and Ṽk and Q k (sub-grid scale space

): (V, Q) = (V h ⊕ V k , Q h ⊕ Q k ) .
Thus the unknowns of the problem (v, p) are replaced by (v h + ṽ, p h + p). The weak formulation is written using the same decomposition. Rather than solving it directly, the author of [Principe, 2008] transforms the equation to the reference domain, then a Fourier analysis of the sub-scale problem is used as proposed in [Codina, 2002].

After a lengthy calculation, a 2D stabilization parameter is computed in terms of the velocity, viscosity and the mesh size.

We adapted it to our needs by supposing that the viscosity and velocity take respectively the following forms η =

    η 1 η 12 η 13 η 12 η 2 η 23 η 13 η 23 η 3     , v =     v 1 v 2 v 3     .
The new stabilization term, customized for anisotropic consistency and mesh orientation can be expressed as:

τ c = (4M η + c) -1/2 where c = 4 v 2 1 h 2 1 + v 2 2 h 2 2 + v 2 3 h 2 3
( 5.32) in which M is the mesh metric and h 1 , h 2 and h 3 are the mesh sizes in the different directions. Note that since we are solving anisotropic Stokes equations, the convection term c is dropped in the following applications.

Applications

As explained previously, the main focus behind the anisotropic solver was to recreate a bilateral sliding contact. Nevertheless, this solver can be of a great importance for other applications in particular for anisotropic flows of viscoplastic materials. We begin by applying the anisotropic solver on a free flow case for two reasons: First, to insure that the solver is perfectly running and second to give a general idea on how it can be used in other applications. Afterward, the solver is used to model bilateral sliding contact in large deformation problems. The center of attention is to evaluate this solver e.g if it is general enough to offer the needed information or if it requires additional manipulation.

Free flow

In the computational domain [0; 1] × [0; 1] × [0; 1] m 3 , an anisotropic material is filled initially in the back corner of the domain to take the form of a quarter of a cylinder (symmetries are applied on both sides). The radius R and height H are respectively equal to 0.4 m and 0.85 m (see Figure 5.10). Since the viscosity in the x-direction is bigger than the one in the y-direction, the material is thicker and should move slower in this direction. .12 completes this discussion. To visualize even more the progress of the anisotropic material, the zero-isovalue is presented at t ≈ 2.1 s. The evolution is indeed more visible in the y-direction. When inspected closely, we noticed an uneven top (Figure 5.12a). This is not considered as an unusual behavior. In fact, since the material advances more in one direction under nothing other than gravity, less material is located in the center which means less support. This can generate a slack region in the center (the corner vicinity). Though this case is simple and maybe considered academic, it is the perfect example to illustrate the performance of the anisotropic solver. Its behavior is overall satisfactory.

Towards contact modeling

An intermediate step before applying the anisotropic solver to contact problem is presented hereby. The computational domain [0; 1] × [0; 1] × [0; 1] m 3 is formed by two sub-domains, the deformable material (in red) and the lubricant (in blue) covering the top and the bottom (Figure 5.13). Note that the lubricant in this case is no longer represented implicitly via a quadratic mixture law. Instead it is explicitly defined as a sub-domain of thickness 0.1 m. On the domain's upper and lower boundaries, we impose velocities of respectively 5 m/s and -5 m/s. These boundary conditions emulate in the simplest way a contact problem without actual contact between tools and the deformable body.

Note that the chosen lubricant thickness is somewhat larger than what is expected in reality, but this case is purely academic to test the behavior of the anisotropic solver. Note that we fixed η lub 3 = 10 7 P a.s, and we control the anisotropy by changing η lub 1 and η lub 2 . For instance, η lub 1 and η lub 2 takes three different values: 10 7 P a.s, 10 5 P a.s and 10 3 P a.s. For the first value, the lubricant consistency is isotropic. The resulting velocity profile should be the same as for a sticking friction. The latter values suggest that the anisotropic ratios are respectively 10 2 and 10 4 . We should expect more sliding impacting the velocity profiles. Figure 5.14 illustrates the velocities plotted on an outer section in the thickness positioned on the middle of the cube side. Notice that for η lub 1 = η lub 2 = 10 7 P a.s, the velocity profile is perfectly parabolic. This is the kind of response we get when sticking contact is modeled. When decreasing the lubricant consistency, the anisotropy ratio increases and more sliding is allowed on the edges of the deformable body. More sliding results in bigger velocities on the edges and the velocity profile (in the material) is more and more flat with the decreasing consistency. The velocity curves in the deformable body confirm our analysis. Note that the velocity picks are located in the lubricant domain which is perfectly normal. Smaller the lubricant consistency, lesser it resists to displacement (e.g bigger the velocity) . To conclude, the behavior of the anisotropic solver is once again proven logic and satisfactory.

Contact

Now that the well-behavior of the anisotropic solver is approved, we attempt using it to model a 2D-contact problem. The lubricant phase is presented implicitly using the quadratic mixture law (contrary to the previous case). Combined together, we await an imitation of a bilateral sliding contact. The monolithic anisotropic consistency tensor η defined all over the domain is obtained in a delicate manner. In the normal direction, it is recovered using a linear mixture law. Whereas in the tangential directions, a quadratic mixture law enabling η lub is used: 5.36) where η quad and η lin are the quadratic and linear mixture laws defined all over the domain.

η =   η quad 0 0 η lin   (
To fix the ideas, we should mention that the anisotropy in η is a direct outcome from the lubricant anisotropy: In this 2D contact problem, the upper tool moves downward with a velocity of 0.3 m/s. The lubricant, deformable body and tools consistencies are respectively equal to 700 P a.s, 10 4 P a.s and 10 5 P a.s.

The results were inspected and did not meet our expectations. The overall results issued from this case were similar to the ones obtained using a linear mixture law e.g. sticking contact.

The cause was immediately identified. Figure 5.15 illustrates the consistency computed all over the domain using the anisotropic solver and a quadratic mixture law. A zoom focused on the corner uncovers an irregular zone in which the consistency is unexpectedly large. To complete the discussion, the flow material presented using arrows is illustrated in Figure 5.16 at t = 0.014 s. In the corner, when using the anisotropic solver the deformable body is unable to slide properly due to the barrier created by the large consistency. The flow should be more like the one issued when using the quadratic law.

Conclusion 169

Though the behavior of the anisotropic solver was proven satisfactory for other applications, it remains limited when it comes to contact modeling. This limitation arises due to the level set gradient. The computed gradient shows singularities and leads to an ill-defined consistency tensor on the corners. It leads to rotated consistency in the x direction forming a barrier preventing the material from sliding freely in the tangential direction.

Conclusion

This chapter attempts to model friction and contact in an Eulerian framework. We presented a quadratic mixture law to introduce a lubricant implicitly between the tool and the deformable body. The lubricant consistency and thickness are deduced by an identification technique with the friction parameter in a Norton-Hoff friction law. We carry out a ring test case and the results were compared with Forge®. They were found comparable but highly dependent of the boundary layer thickness. For a small thickness, the lubricant consistency is overpowered by the tools/deformable body consistencies and consequently the boundary layer is not properly presented. However for a large thickness, an odd behavior is detected near the corners.

The development of the directional solver was meant to model bilateral sliding contact. A VMS stabilization technique was used to insure the stability. The new solver was validated with two different cases: i) in a free flow simulation and ii) an explicit contact case. However when applying the directional solver to a 2d case, where the friction is represented implicitly using a quadratic mixture law, some difficulties were confronted. The gradient has proven to be the problem, it introduces a rotated consistency values on the corner of the deformable body. It acts like a barrier preventing the sliding of the deformable body. An alternative method, non dependent of the mesh ( here we mean the interface thickness of the level set ), is the use of Navier boundary conditions. The latter boundary condition is a combination between a Dirichlet conditions and the Newman ones. However to apply it, we should compute explicitly the contact surface between the tool and the deformable body.

Résumé en Français

Ce chapitre est dédié à la modélisation du frottement et du contact dans un cadre Eulérien. Nous avons présentés une loi de mélange quadratique pour introduire implicitement un lubrifiant entre l'outil et le corps déformable. La viscosité et l'épaisseur du lubrifiant sont déduites par une technique d'identification avec le paramètre de frottement dans une loi Norton-Hoff. En utilisant le cas de déformation d'un anneau, les résultats de notre approche Eulérienne utilisant la loi quadratique ont été comparés ave Forge. Les résultats sont satisfaisants mais fortement dépendant de l'épaisseur de la couche limite. Pour une petite épaisseur, les consistances de l'outil et du corps déformable emportent sur la viscosité du lubrifiant. Par conséquent, la couche limite n'est pas proprement représentée. Pour surmonter ce problème, le maillage est la réponse-clé pour une propre représentation. Dans le but de modéliser un problème de contact glissant bilatéral, nous avons proposé un solveur directionnel stabilisé par la technique VMS. Le nouveau solveur a été validé avec deux cas différents : i) dans une simulation d'un écoulement libre d'un matériau anisotrope et ii) un cas de contact explicite. Cependant en appliquant le solveur directionnel à cas 2D, où le frottement est représenté implicitement en utilisant une loi de mélange quadratique, quelques difficultés ont été confrontées. Le gradient s'avérant être l'origine, il agit comme une barrière empêchant le glissement du corps déformable. Pour surmonter cette difficulté, nous devons revisiter la définition des normales dans notre formulation directionnelle.

[ Zhang andFang, 2008] Zhang, T. andFang, D. (2008). Global wellposed problem for the 3-d incompressible anisotropic navier-stokes equations. Journal de Mathématiques Pures et Appliquées, 90(5):413 -449.

Conclusions 1 Conclusions

This work is one of the first attempts studying the feasibility of multi-domain large deformation problems using an unusual approach, "the monolithic Eulerian technique". The resulting model is desired to be able of i) following complex geometries evolution during time, ii) describing multi-domain problems, iii) modeling friction and contact and iv) giving results more or less comparable to real applications (using Forge ® as the base of our superpositions). All the mentioned objectives has to be accomplished in a highly parallel scalable environment.

In other words, this work is the first step toward an Eulerian Forge ® : It has to assess if the current Eulerian approach -available in CIMLib-is fit to large deformation modeling.

It has to propose and test feasible development and recommend answers to unsolved aspects for future works.

It has to conclude if the resulting approach has industrializing potential on the long run.

After an overview of the different approaches fit for large deformation modeling, the literature confirmed that the Eulerian approach is rarely used for these kind of problems (chapter 1). Therefore, adopting the Eulerian technique will require an outside of the box thinking to find answers to our problems. However, the undeniable potential of this method makes it a great candidate (e.g the capacity of managing automatically contact detection). The monolithic approach -base of the CIMLib-is described in details in the first part of chapter 2.

To adapt the standard version to large deformation problems, additional developments and manipulations are needed . The second part of chapter 2 is fully dedicated to these improvements:

A new quadratic mixture law was implemented. It offers the option of adding a new phase between two interacting bodies to simulate the presence of a lubricant for instance. Note that this is a great deal when attempting large deformation modeling.

Especially when there is a need to model friction between different domains.

A new manner to impose boundary conditions was presented. It enables the users to let the air escape from the computational domain emulating actual vents.

Technical ameliorations were proposed as well. We recommend an optimal choice of parameters such as i) mixture law weight, order of interpolation and a variable time step to improve volume conservation and ii) mesh size choice for the stabilization of the convection solver.

chapter 3 and chapter 4 were entirely dedicated to numerical applications. The first is devoted for the parallel environment in CIMLib and its efficiency. The overall performance was found satisfactory. The second presents different applications to evaluate the capacity of the approach in several scenarios: i) deforming complex geometries such as crankshafts and connecting rod, ii) deforming simple geometries but stacked in great numbers (up to 256 deformable geometries). The results were validated via Forge ® simulations and found very satisfactory. In some cases (with high numbers of deformable bodies), the performance was even more resilient than Forge ® . After these comparisons, we were reassured that this approach exhibits real potential to compete with an industrial software. To push even further, attempts to model friction and contact problems in an Eulerian approach were the subject matter of chapter 5. Simple yet effective, propositions were formed:

To model a sticking friction, a linear mixture law is the answer. The transition between a high consistency (for the solid) and a lower consistency (for the deformable body) results with a high resistance emulating numerically friction. The results were validated by deforming a ring between two rigid tools in both CIMLib and Forge ® . The results were a perfect match.

To control friction, a quadratic mixture law is used enabling a fictitious boundary layer imitating the presence of a lubricant between the different bodies. The lubricant viscosity is determined in terms of the deformable body consistency and the boundary layer thickness. This relationship was established via a superposition with Lagrangian friction to make sure the same cases are modeled. Again, a comparison 2. Recommendation and Perspectives 175 with simulations carried out in Forge ® were established. The results showed the same behavior but exhibited some limitations as well for certain cases. The results were highly dependent of the boundary layer thickness. If the needed thickness is very small, the mesh may not detect the boundary layer and the lubricant consistency is not reached. In this case, the same results as in sticking friction were found.

If the thickness is increased to prevent this problem but a big gap is present between the different consistencies, an unrealistic behavior is detected in the vicinity were the different values meet (on the corners).

Even though the quadratic law was proven effective for some cases in controlling friction, the overall performance cannot be generalized.

The directional solver is introduced to model bilateral sliding contact. The consistency is expressed as a tensor. Although our intentions were to model contact in a general way, the solver can be used for other high-demand applications such as anisotropic polymer flow. To test the well-functioning of the solver, it was utilized for a simple anisotropic free flow case. Its behavior was found impressive. High hopes were expected for contact applications. When testing it in a simple 2D case a draw-back was noticed. The directional formulation does not take into account the critical zone (on the corners) where points may have an ill-defined normals. A rotated consistencies resulted near the corner and limited the proper progress and sliding of the material.

To avoid this matter, this crucial zone should be properly managed. Several steps are proposed such as a new formulation using Navier slipping conditions.

To conclude, the Eulerian approach express high potential to compete, in the long run, with a perfectly mastered software such as Forge ® . Additional developments are needed. Still, its capacity is undeniable especially when taking into consideration that few works deal with these large deformation problems.

Note that during this study only materials with Newtonian behavior law were considered. Thermal effects were neglected as well. Which brings us to the next point: the recommendations and perspectives.

Recommendation and Perspectives

To recapitulate, future works should be directed towards :

Studying the feasibility of multi-domain contact problem in this Eulerian approach:

• The directional solver should be revisited and explored to asses if it will be ever able to describe properly contact problems. Defining the normal using the level set of the deformable body leads to a rotated consistency on the corners preventing the sliding. Redefining the normal in a new manner could be the solution to this limitation.

• In term of friction, we were able to represent sliding. However, representing the consistency in the boundary layer was challenging. The adopted method relies heavily on the lubricant consistency and thickness (η lub , ε). An alternative method is to introduce Navier slip/no-slip conditions. It does not depend of the mesh (especially the thickness interface ε). However the use of such method requires an explicit representation of the contact zone. In other words, the zero isovalue should be reconstructed to be defined on nodes as in Lagrangian formulations.

Varying Norton-Hoff laws with a sensibility parameter m different then 1.

We should note that in order to pursue this point, the Norton-Hoff law may need need regularization. It can be accomplished using a hyperbolic tangent function for instance.

Taking into account thermal effects by coupling diffusion-convection equations with the current mechanical system.

In an Eulerian framework, coupling the thermal and mechanical problem can be more beneficial than in a Lagrangian formulation. It takes into account different interacting heterogeneities and enables us to study/follow the thermal transfers easier.

Optimizing different algorithm affecting the computing time and the parallel efficiency such as the convection solver assembly and the balancing load algorithm.

• The assembly time of the LevellerT solver is very comparable to a Stokes problem. In fact, the convection problem is relatively small when compared to the mechanical problem (Stoke or Navier-Stokes) and should not consume a large amount of time. We recommend thus, rewriting the solver in order to optimize the assembly time.

• The parallel strategy of the adaptation algorithm should be improved as well.

In a parallel framework, the adaptation is executed twice. In addition, the mesh generator acts locally on nodes and therefore needs an excessive interface displacement (interfaces between cores). The latter strategy should be revisited to reduce its high cost. Modelling multi-domain large deformation problems using an Eulerian monolithic approach in a massively parallel environment ABSTRACT : Modeling of multi-domain problems is addressed in a Purely Eulerian framework. A single mesh is used all over the domain. The evolution of the different interacting bodies is described using numerical tools such as the Level Set method. The characteristics of the subdomains, considered as heterogeneities in the mesh, are determined using mixture laws. This work is one of the first attempts applying fully Eulerian Approach to Model large deformation problems. Therefore, the capacity of this approach is tested to determine necessary developments. The friction between the different objects is managed by adding a boundary layer implying the presence of a lubricant. Combined with an identification technique, a new quadratic mixture Law is introduced to determine the lubricant viscosity. Comparisons have been performed with Forge ® and results were found satisfactory. To treat the contact problem between the different objects, a directional solver was developed. Despite the interesting results, it remains the topic of further improvements. The scalability of the approach in a massively parallel environment is tested as well. Several recommendations were proposed to ensure an optimal performance. The technique of a single mesh guarantees a very good scalability since the efficiency of parallelism depends of the partition of a single mesh (unlike the Lagrangian Methods). The proposed method presents undeniable capacities but remains far from being complete. Ideas for future Improvements are proposed accordingly.

Modélisation des problèmes de grandes déformations multi-domaines par une approche

Keywords : monolithic approach, Level Set, anisotropic mesh adaptation, quadratic mixture law, directional solver, parallel computing.

x 2 . 7

 27 Illustration of a 2d flat bunch simulation involving two rigid bodies (tools) and a deformable body. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.8 Different Types of mixture laws: (a) P0-mixture law and (b) P1-mixture law . .

3 . 2

 32 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 The Flynn classification of computer architectures: (a) SISD, (b) MISD, (c) SIMD and (d)MIMD. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . The memory classification of parallel computers: (a) shared memory, (b) distributed memory and (c) hybrid memory. . . . . . . . . . . . . . . . . . . . .

xi 3 . 7 4 . 2 4 . 3 4 . 4 4 . 5

 3742434445 Local matrix assembly of the two sub-domains: red and blue (a-b) and their projection into the global matrix (c) . . . . . . . . . . . . . . . . . . . . . . . 3.8 Physical domain used to study the parallel performance of the monolithical approach. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.9 Speed up (a-c) and efficiency (b-d) for the mechanical and transport solvers and different important part of the code obtained on the coarse mesh using the 1.2.6 version of MPI. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.10 Speed up (a-c) and efficiency (b-d) for the mechanical and transport solvers and different important part of the code obtained on the coarse mesh using the 1.4.3 version of MPI with the bind to core option. . . . . . . . . . . . . . . . . . . . 3.11 Speed up (a) and efficiency (b) for different part of the code obtained on the coarse mesh using the openmpi round robin option. . . . . . . . . . . . . . . . 3.12 Using the refined mesh: Speed up (a-c) and efficiency (b-d) for the mechanical and transport solvers and different important part of the code using the 1.4.3 version of MPI with the bind to core option . . . . . . . . . . . . . . . . . . . 3.13 Speed up (a) and efficiency (b) for the anisotropic mesh adaptation technique used in CIMLib. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.14 Dynamic load balancing algorithm used in CIMLib. . . . . . . . . . . . . . . . 4.1 Examples of a connecting rod (a) and a crankshaft (b). . . . . . . . . . . . . . Geometry of the connecting rod immersed in the computational domain and adapted using an anisotropic mesh ( t = 0 s ). . . . . . . . . . . . . . . . . . . The contours of the deformable body and tools presented on a longitudinal section in the thickness at t=0 s . . . . . . . . . . . . . . . . . . . . . . . . . . . . . From left to right: the isovalue zero of the connecting rod respectively for t = 0 s and t = 2.6 s. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Comparison of the results obtained with CIMLib and a Forge ® simulation at t ≈ 1.6 s. The same longitudinal section in the width is used in both softwares: (a) and (b) illustrate the mesh, (c) and (d) illustrate the pressure, (e) and (f) illustrate the velocity component v x and (g) and (h) illustrate the velocity component v y . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.6 Comparison of the the velocity component v z obtained with CIMLib (a) and Forge ® (b) at t ≈ 1.6 s. The same longitudinal section in the thickness is used in both softwares. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

xii 4 .

 4 10 (a), (

  The deformed parallelepiped using an anisotropic adapted mesh. (a), (b), (c) and (d) represent respectively the maps of the velocity components v x , v y , v z and the pressure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.18 Forge ® failed attempt to partition the domain containing 148 spheres on 8 cores.

4. 19 5 . 2 5 . 3 5 . 4

 19525354 From left to right: The partition of a domain containing 148 spheres in Forge ® taking into account the above improvment using respectively 4 and 8 cores. . . . . . . 4.20 From left to right: The partition of a domain containing 148 spheres in CIMLib using respectively 4 and 8 cores. . . . . . . . . . . . . . . . . . . . . . . . . . 4.21 Initial partition of a domain containing 256 spheres on 64 cores. . . . . . . . . 4.22 Velocity components v x and v z are respectively represented on a longitudinal section in the thickness in (a) and (c). Velocity component v y represented on a longitudinal section in the width in (b) (t = 0.64 s). . . . . . . . . . . . . . . . 5.1 Graphic representation of the different friction laws. . . . . . . . . . . . . . . . Schematic representing the boundary layer. . . . . . . . . . . . . . . . . . . . An illustration of the initial ring geometry along with the symmetry planes . . . The deformed ring geometry at t = 3.467 s . . . . . . . . . . . . . . . . . . . . 5.5 (a) illustrates the mixture law on a cutting plane in the thickness. (b) to (e) represent the different velocity components and magnitude. (f) illustrates the pressure repartition. The different maps are presented on a forth of the ring at t = 3.46 s. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.6 The deformed ring outline for y = 0 at t = 3.467 s . . . . . . . . . . . . . . . . Comparison of the velocities obtained with sticking friction in both Forge ® and CIMLib on a cutting plane in the thickness for y = 0 at t = 3.46 s : (a) and (b) compare the velocity v x plotted on the outer outline using two different mixture thickness ε = 0.5 mm (on the left) and ε = 0.25 mm (on the right). (c) compares the velocity v z plotted on the outer outline using a mixture thickness ε = 0.25 mm. (d) superposes the velocity v x plotted on the inner outline of the ring using ε = 0.25 mm. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.8 Comparison of the velocities obtained with sliding friction in both Forge ® and CIMLib. Curves are plotted at t = 3.46 s for α f = 0.3 on a cutting plane in the thickness for y = 0 : (a) and (b) compare the velocity v x respectively plotted on the inner and outer outline using a mixture thickness ε = 0.25 mm. (c) and (d) compares the velocity v z respectively plotted on the inner and outer outline using a mixture thickness ε = 0.25 mm. . . . . . . . . . . . . . . . . . . . . . . . . 5.9 (a) and (c) present a comparison of the velocities v x obtained with sliding friction in both Forge ® and CIMLib using respectively ε = 0.25 mm and ε = 0.5 mm . Curves are plotted at t = 3.46 s for α f = 0.03 on the outer outline of a cutting plane in the thickness for y = 0. (b) and (d) present a zoom of the quadratic mixture law using respectively ε = 0.25 mm and ε = 0.5 mm . . . . . . . . . . . 5.10 Initial position of the anisotropic material. . . . . . . . . . . . . . . . . . . . . 5.11 Maps of the velocity components v x and v y at t ≈ 1 s . . . . . . . . . . . . . . 5.12 Side and upper views of the final position of the material corresponding to t ≈ 2.1 s. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.13 The computational domain formed by two phases: the deformable body and the lubricant. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.14 The impact of a variant lubricant anisotropic consistency on the velocity profile v x . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.15 At t = 0 s, a superposition between the anisotropic consistency (on the left) and the consistency using the quadratic mixture law (on the right). . . . . . . . . . 5.16 At t = 0.014 s, a superposition between the material flow (illustrated by arrows) using an anisotropic consistency (on the left) and a quadratic mixture law (on the right). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Efficiency of the contact analysis algorithm in Forge ® . . . . . . . . . . . . . . 22 2.1 Parameters used in the flat bunch simulation . . . . . . . . . . . . . . . . . . 67 2.2 Comparison between the different definitions to compute the simplex size. . 763.1 The simulations parameters used for speed up tests . . . . . . . . . . . . . . . 106 3.2 Time (in s) spent to assemble, solve the mechanical problem with two different solvers ( Navier-Stokes vs Stokes ) using the 1.2.6 version of MPI . . . . . . . 106 3.3 Time (in s) spent to assemble, solve the mechanical and transport solver and to accomplish the simulation using the 1.2.6 version of MPI . . . . . . . . . . . . 107 3.4 Time (in s) spent to accomplish other relevant part of the code during simulation using the 1.2.6 version of MPI . . . . . . . . . . . . . . . . . . . . . . . . . . 108 3.5 Time (in s) spent to assemble, solve the mechanical and transport solver and to accomplish the simulation using the 1.4.3 version of MPI with the bind to core option . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 3.6 Time (in s) spent to accomplish other relevant part of the code during simulation using the 1.4.3 version of MPI with the bind to core option . . . . . . . . . . . 109 3.7 Size of the Output files in parallel computing . . . . . . . . . . . . . . . . . . . 111 3.8 Time (in s) spent to assemble, solve the mechanical and transport solver and to accomplish the simulation using the 1.4.3 version of MPI with the opemmpi round robin option . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112 3.9 Time (in s) spent to accomplish other relevant part of the code during simulation using the 1.4.3 version of MPI with the openmpi round robin option . . . . . . 112 3.10 Using the refined mesh: Time ( in s) spent to assemble, solve the mechanical and transport solver and to accomplish the simulation using the 1.2.6 version of MPI 114

Figure 1 :

 1 Figure 1: Partition of a domain formed by 10 deformable bodies using a Lagrangian approach.

Figure 1 . 5 :

 15 Figure 1.5: Contact spatial search analysis in Forge ® .

Figure 1 . 6 :

 16 Figure 1.6: Illustration of the gap function.

  Figure 1.7: Construction of contact fictitious elements

Figure 1 . 8 :

 18 Figure 1.8: SPMD strategie for parallel computing used in Forge ® .

Figure 1 .

 1 Figure 1.10: Decoupled resolution using an ALE approach .

Figure 1 .

 1 Figure 1.11: Time gain for different algorithm .

[

  Duvaut and Lions, 1972] Duvaut, G. and Lions, J.-L. (1972). Les inéquations en mécanique et en physique. [El-Abbasi and Bathe, 2001] El-Abbasi, N. and Bathe, K.-J. (2001). Stability and patch test performance of contact discretizations and a new solution algorithm. Computers & Structures, 79(16):1473 -1486.

Figure 2 .

 2 Figure 2.1: A sphere immersed in a cubic computational domain Ω. The signed distance function isovalues of the central slice are projected onto the bottom of the domain.

Figure 2 . 2 :

 22 Figure 2.2: Comparison between a signed distance function (blue) and a hyperbolic tangent function (red)

Figure 2 . 3 :

 23 Figure 2.3: Different weight of the shape functions using a Galerkin approximation and a Streamline Upwind Petrov-Galerkin .

Definition 2 .

 2 3.0.1. Let (N , T ) denotes respectively the nodes and the set of simplex whose vertices belong to N . Considering F the set of faces of the simplex, T is a mesh topology of Ω if and only if:(i) Each face of F shares at most two elements.

Figure 2 . 4 :

 24 Figure 2.4: Mesh generation by local optimization using "MTC"

Figure 2 . 6 :

 26 Figure 2.6: The anisotropic mesh obtained when using the edge based error estimator metric: (a-b) for a cube, (c) for the upper tool of the connecting rod. .

Figure 2 . 7 :

 27 Figure 2.7: Illustration of a 2d flat bunch simulation involving two rigid bodies (tools) and a deformable body.

Figure 2 . 8 :

 28 Figure 2.8: Different Types of mixture laws: (a) P0-mixture law and (b) P1-mixture law

Figure 2 . 9 :

 29 Figure 2.9: Comparison between a linear mixture law (blue) and a quadratic mixture law(red)

Figure 2 .

 2 Figure 2.10: The effect of using a variable time step on the mass conservation..

Figure 2 .

 2 Figure 2.11: The effect of using P1mixture law instead of a P0 mixture law on the mass conservation.

Figure 2 .

 2 Figure 2.13: Flat bunch simulation using a linear coefficient for mixture law.

Figure 2 . 14 :

 214 Figure 2.14: Representation of the mixture law for a flat bunch simulation using different coefficient weights: (a) linear -(b) inverse -(c) geometric -(d) logarithmic .

Figure 2

 2 Figure2.14 shows a comparison between the four mixture law presented above. As mentioned previously when using linear coefficients, the body with the higher consistency

Figure 2 .

 2 Figure 2.15: Flat bunch simulation using a inverse coefficient for mixture law.

Figure 2 .

 2 Figure 2.16: Flat bunch simulation using a logarithmic coefficient for mixture law.

Figure 2 .

 2 Figure 2.17: Top view of the zero iso-value of a flat bunch simulation in 3d showing the mixture law when using: a-linear coefficient, b-logarithmic coefficient.

Figure 2 .

 2 Figure 2.18: Instability of the levellerT solver.

Figure 2 .

 2 Figure 2.19: Two different ways to compute the mesh size .

Figure 2 . 20 :

 220 Figure 2.20: Comparison between two methods to compute the mesh size in CIMLib: (a) using the Edge min definition and (b) using the Altitude min definition.

Figure 2 .

 2 Figure 2.21 illustrates the evolution of the zero iso-value of the deformable body at different time steps.

Figure 2 .

 2 Figure 2.21: A 3d flat bunch simulation after the ameliorations (logarithmic weigh, new mesh size definition ...).

Figure 2 .

 2 Figure 2.22: Domain considered using a porous upper tool to show the Eulerian approach capacity to detect trapped air .

Figure 2 .

 2 Figure 2.23: Initial mesh and mixture law for the problem in Figure 2.22.

Figure 2 .

 2 Figure 2.24: The effect given by the imposition of the kinematic on the inner of the upper tool.

  83) is equivalent to impose the kinematic only on the vicinity of the contact zone and letting the air escape otherwise.To test the efficiency of the proposed porous boundary condition, the same case treated earlier is used. This time the kinematic are imposed on Γ p instead of Γ inner .

Figure 2 .

 2 Figure 2.25: The evolution of the porous boundary condition in time.

Figure 2 .

 2 Figure 2.26: The evolution of the zero iso-value of the deformable body when using the porous boundary condition.

3. 4

 4 Parallel computing in CIMLib . . . . . . . . . . . . . . . . . . . 100 3.4.1 S.P.M.D parallelization of the Navier-Stokes Solver . . . . . . . 101 3.4.2 Dynamic load balancing . . . . . . . . . . . . . . . . . . . . . . 103 3.5 Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104

Figure 3 . 1 :

 31 Figure 3.1: The Flynn classification of computer architectures: (a) SISD, (b) MISD, (c) SIMD and (d)MIMD.

3. 2 . 3 . 4 .

 234 Concepts and Terminology93 Single-Instruction, Multiple-Data (SIMD) : It uses multiprocessing element in order to execute one single instruction using multiple data storage. In fact, each processing element obtains the same instruction to execute, but loads a different and separate data. This type of architecture is compatible with a shared or distributed memory. Multiple-Instruction, Multiple-Data (MIMD) : It uses multiprocessing element capable of executing multiple instruction for separate data memory units. Each processing element loads a different instruction, executes it to a different data memory and stores the result back into the corresponding storage.

Figure 3 . 2 :Figure 3 . 3 :

 3233 Figure 3.2: The memory classification of parallel computers: (a) shared memory, (b) distributed memory and (c) hybrid memory.

Figure 3 . 5 :

 35 Figure 3.5: Different types of speed up when running a parallel simulation : the common speed up (blue); the linear speed up (red); the super-linear speed up (green)

3. 4 .

 4 Parallel computing in CIMLib 101 parallelization has a great importance.

( a )Figure 3 . 6 :

 a36 Figure 3.6: Example of partitioning a global mesh in two sub-meshes

  .6 and Figure 3.7).

( a )

 a Matrix associated to the red sub-domain (b) Matrix associated to the blue sub-domain (c) Global Matrix

Figure 3 . 7 :

 37 Figure 3.7: Local matrix assembly of the two sub-domains: red and blue (a-b) and their projection into the global matrix (c)

y

  and r are two global vectors defined on the whole Ω.

Figure 3 . 8 :

 38 Figure 3.8: Physical domain used to study the parallel performance of the monolithical approach.

  Figures 3.9c and 3.9d show the performance of the remaining part of the simulation. The speed up of the function implemented in CIMLib to compute the generalized strain rate is acceptable. But, when it comes to the rest, a big loss of scalability is shown. The sources of this loss should be investigated in order to improve the overall scalability of the simulations.

Figure 3 . 9 :

 39 Figure 3.9: Speed up (a-c) and efficiency (b-d) for the mechanical and transport solvers and different important part of the code obtained on the coarse mesh using the 1.2.6 version of MPI.

Figure 3 . 10 :

 310 Figure 3.10: Speed up (a-c) and efficiency (b-d) for the mechanical and transport solvers and different important part of the code obtained on the coarse mesh using the 1.4.3 version of MPI with the bind to core option.

Figure 3 . 11 :

 311 Figure 3.11: Speed up (a) and efficiency (b) for different part of the code obtained on the coarse mesh using the openmpi round robin option.

Figure

  Figure 3.12 completes this discussion. The Stokes resolution remains the least scalable with an efficiency of 0.8. It is interesting to mention that after 32 cores LevellerT changes behavior and becomes less scalable. It is totally expected for small problems such as the convection problem (i.e LevellerT). Increasing the number of cores is efficient to a certain limit; afterward the scalability decreases due to an increasing communication between cores.

Figure 3 . 12 :

 312 Figure 3.12: Using the refined mesh: Speed up (a-c) and efficiency (b-d) for the mechanical and transport solvers and different important part of the code using the 1.4.3 version of MPI with the bind to core option

Figure 3 . 13 :

 313 Figure 3.13: Speed up (a) and efficiency (b) for the anisotropic mesh adaptation technique used in CIMLib.

Figure 3 .

 3 Figure 3.14: Dynamic load balancing algorithm used in CIMLib.

  Besides pistons, crankshafts and connecting rods are the main components of the reciprocating piston engine used in cars for instance.The connecting rod ( Figure4.1a ) is the element joining the piston to the rotating crankshaft ( Figure4.1b ).

Figure 4 . 1 :

 41 Figure 4.1: Examples of a connecting rod (a) and a crankshaft (b).

Figure 4 . 2 :

 42 Figure 4.2: Geometry of the connecting rod immersed in the computational domain and adapted using an anisotropic mesh ( t = 0 s ).

Figure 4 . 3 :

 43 Figure 4.3: The contours of the deformable body and tools presented on a longitudinal section in the thickness at t=0 s

Figure

  Figure 4.4 represents the zero iso-value of the deformable body for t = 0 s ( on the left ) and t = 2.6 s ( on the right ). We mention that, based on the same analogy presented in chapter 2, the air is allowed to evacuate by imposing a porous boundary conditions. The latter is equivalent to defining vents in the tools.

Figure 4 . 4 :

 44 Figure 4.4: From left to right: the isovalue zero of the connecting rod respectively for t = 0 s and t = 2.6 s.

Figure 4 .

 4 5 and Figure 4.6.

Figure 4 .

 4 Figure 4.5 illustrates the pressure and velocities maps for both software for t ≈ 1.6 s.The chosen cut corresponds to the mid-plane of the deformable body. Results are very similar and represent the same features. The main difference, thought not significantly important, is mostly located in the center. For instance, the velocity along the x-axis represents a gap of 1.7 mm/s. This is directly related to the meshes illustrated in Figures 4.5a and 4.5b. In the CIMLib case, lesser nodes are used in the center leading to this difference. The same reason is behind the difference for the pressure and the velocity

Figure 4 .

 4 Figure 4.6 completes the comparison. The velocity component v z is represented on a longitudinal section in the thickness. Once again, the results are found very comparable.

Figure 4 . 5 :

 45 Figure 4.5: Comparison of the results obtained with CIMLib and a Forge ® simulation at t ≈ 1.6 s. The same longitudinal section in the width is used in both softwares: (a) and (b) illustrate the mesh, (c) and (d) illustrate the pressure, (e) and (f) illustrate the velocity component v x and (g) and (h) illustrate the velocity component v y .

Figure 4 . 6 :

 46 Figure 4.6: Comparison of the the velocity component v z obtained with CIMLib (a) and Forge ® (b) at t ≈ 1.6 s. The same longitudinal section in the thickness is used in both softwares.

Figure 4 . 7 :

 47 Figure 4.7: Geometry of the crankshaft immersed in the computational domain between the tools represented by the anisotropic mesh ( t = 0 s ).

Figure 4 . 8 :

 48 Figure 4.8: Mesh projected on the zero iso-value of the upper tool

Figure 4 . 9 :Figure 4 .

 494 Figure 4.9: Geometry of the crankshaft at t ≈ 0.9 s

Figure 4 .

 4 Figure 4.11: (a) and (b) represent respectively the contours of the crankshaft between the tools at t ≈ 0.9 s and t ≈ 0.96 s. A growing penetration of the upper tool in the deformable piece is noticeable

Figure 4 .

 4 Figure 4.12: The zero-isovalue of the connecting rod between the tools at t ≈ 1.6s . No penetration between the different bodies is noted.

Figure 4 .

 4 Figure 4.13: Initial geometry of the 10 parallelepiped immersed in the computational domain.

Figure 4 .

 4 Figure 4.14: The uniform mesh of the domain presented on a section: Inside the virtual box, the mesh size is almost ten times smaller than the external one.

Figure 4 .

 4 Figure 4.15 represents the deformed geometries at t = .0.7 s The simulation was carried out on 256 cores with no problem what so ever. Nonetheless, the results were found unsatisfactory. The pieces did not deform properly: their top remained intact during the whole simulation. When examined closely, this irregularity is associated to the mesh size h int . Even though it is small, it remains significantly large when compared with the mixture law thickness ≈ 0.9 mm . If one decreases even more the mesh size, the problem will be prevented but it will be costly in computing time and the number of cores needed to run the simulation.

Figure 4 .

 4 Figure 4.16: The initial mesh adapted anisotropically around the deformable parallelepiped (on the left) with a zoom (on the right).

Figure 4 .

 4 Figure 4.17: The deformed parallelepiped using an anisotropic adapted mesh. (a), (b), (c) and (d) represent respectively the maps of the velocity components v x , v y , v z and the pressure.

Figure 4 .

 4 Figure 4.18: Forge ® failed attempt to partition the domain containing 148 spheres on 8 cores.

Figure 4 .

 4 Figure 4.19: From left to right: The partition of a domain containing 148 spheres in Forge ® taking into account the above improvment using respectively 4 and 8 cores.

Figure 4 .

 4 Figure 4.20: From left to right: The partition of a domain containing 148 spheres in CIMLib using respectively 4 and 8 cores.

Figure 4 .

 4 21 illustrates the whole domain partition. Again, CIMLib confronted no difficulties dividing the workloads properly on 64 cores.

Figure 4 .

 4 Figure 4.21: Initial partition of a domain containing 256 spheres on 64 cores.

  |v x | and |v y | are bigger on the edges and are almost symmetrical suggesting the direction of the material flow. |v z | is bigger on the top and takes the value of the moving upper tool.

Figure 4 .

 4 Figure 4.22: Velocity components v x and v z are respectively represented on a longitudinal section in the thickness in (a) and (c). Velocity component v y represented on a longitudinal section in the width in (b) (t = 0.64 s).

  (a) Tresca friction law (b) Norton-Hoff friction law

Figure 5 . 1 :

 51 Figure 5.1: Graphic representation of the different friction laws.

Figure 5 . 2 :

 52 Figure 5.2: Schematic representing the boundary layer.

Figure 5 . 3 :

 53 Figure 5.3: An illustration of the initial ring geometry along with the symmetry planes .

Figure 5 . 4 :

 54 Figure 5.4: The deformed ring geometry at t = 3.467 s .

Figure 5 . 5 :

 55 Figure 5.5: (a) illustrates the mixture law on a cutting plane in the thickness. (b) to (e) represent the different velocity components and magnitude. (f) illustrates the pressure repartition. The different maps are presented on a forth of the ring at t = 3.46 s.

5. 2 .Figure 5 . 6 :

 256 Figure 5.6: The deformed ring outline for y = 0 at t = 3.467 s .

Figure 5 . 7 :

 57 Figure 5.7: Comparison of the velocities obtained with sticking friction in both Forge ® and CIMLib on a cutting plane in the thickness for y = 0 at t = 3.46 s : (a) and (b) compare the velocity v x plotted on the outer outline using two different mixture thickness ε = 0.5 mm (on the left) and ε = 0.25 mm (on the right). (c) compares the velocity v z plotted on the outer outline using a mixture thickness ε = 0.25 mm. (d) superposes the velocity v x plotted on the inner outline of the ring using ε = 0.25 mm.

Figure 5 . 8 :

 58 Figure 5.8: Comparison of the velocities obtained with sliding friction in both Forge ® and CIMLib. Curves are plotted at t = 3.46 s for α f = 0.3 on a cutting plane in the thickness for y = 0 : (a) and (b) compare the velocity v x respectively plotted on the inner and outer outline using a mixture thickness ε = 0.25 mm. (c) and (d) compares the velocity v z respectively plotted on the inner and outer outline using a mixture thickness ε = 0.25 mm.

  (a) using a mixture thickness ε = 0.25 mm.(b) using a mixture thickness ε = 0.25 mm.(c) using a mixture thickness ε = 0.5 mm. (d) using a mixture thickness ε = 0.5 mm.

Figure 5 .

 5 Figure 5.9: (a) and (c) present a comparison of the velocities v x obtained with sliding friction in both Forge ® and CIMLib using respectively ε = 0.25 mm and ε = 0.5 mm . Curves are plotted at t = 3.46 s for α f = 0.03 on the outer outline of a cutting plane in the thickness for y = 0. (b) and (d) present a zoom of the quadratic mixture law using respectively ε = 0.25 mm and ε = 0.5 mm .

  usually used to reduce the equations into Saint-Venant equations (known as shallow water equations as well).

Figure 5 .

 5 Figure 5.10: Initial position of the anisotropic material.

Figure 5 .

 5 11 confirms this Chapter 5. Friction and contact 164 fact. The velocities maps in both directions are presented at t ≈ 1 s. It is clear that the velocity v y is bigger than the velocity v x , the material moves forward faster in the y-direction .

Figure 5 .

 5 Figure 5.11: Maps of the velocity components v x and v y at t ≈ 1 s

Figure 5 . 12 :

 512 Figure 5.12: Side and upper views of the final position of the material corresponding to t ≈ 2.1 s.

Figure 5 .

 5 Figure 5.13: The computational domain formed by two phases: the deformable body and the lubricant.

Figure 5 .

 5 Figure 5.14: The impact of a variant lubricant anisotropic consistency on the velocity profile v x .

  36) expresses the tensor in the Cartesian frame. To generalize its form, we apply a rotation matrix R.The rotation matrix R is determined using n, the normal vector to the level set function φ: R = (n, n ⊥ )

  Zoom focused on the corner (d) Zoom focused on the corner

Figure 5 .

 5 Figure 5.15: At t = 0 s, a superposition between the anisotropic consistency (on the left) and the consistency using the quadratic mixture law (on the right).

  Zoom focused on the corner (d) Zoom focused on the corner

Figure 5 .

 5 Figure 5.16: At t = 0.014 s, a superposition between the material flow (illustrated by arrows) using an anisotropic consistency (on the left) and a quadratic mixture law (on the right).

  Eulérienne monolithique massivement parallèle RESUME : La modélisation des problèmes multi-domaine est abordée dans un cadre purement Eulérien. Un maillage unique, ne représentant plus la matière, est utilisé. Les différentes frontières et leur évolution sont décrites via des outils numériques tels que la méthode Level Set. Les caractéristiques locales de chaque sous domaines sont déterminées par des lois de mélange. Ce travail est une des premières tentations appliquant une approche Eulérienne pour modéliser de problèmes de grandes déformations. Dans un premier temps, la capacité de l'approche est testée afin de déterminer les développements nécessaires. Le frottement entre les différents objets est géré par un lubrifiant ajouté dans une couche limite. Combinée avec une technique d'identification, une nouvelle loi de mélange quadratique est introduite pour décrire la viscosité du lubrifiant. Des comparaisons ont été effectuées avec Forge ® et les résultats sont trouvés satisfaisants. Pour traiter le contact entre les différents objets, un solveur directionnel a été développé. Malgré que les résultats soient intéressants, il reste le sujet de nouvelles améliorations. La scalabilité de l'approche dans un environnement massivement parallèle est testée aussi. Plusieurs recommandations ont été proposées pour s'assurer d'une performance optimale. La technique du maillage unique permet d'obtenir une très bonne scalabilité. L'efficacité du parallélisme ne dépend que de la partition d'un seul maillage (contrairement aux méthodes Lagrangiennes). La méthode proposée présente des capacités indéniables mais reste loin d'être complète. Des pistes d'amélioration sont proposées en conséquence. Mots clés : approche Eulérienne monolithique, Level Set, adaptation de maillage anisotrope, loi de mélange quadratique, solveur directionnel, calcul parallèle.

  Interface modeling . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.2 Contact formulation . . . . . . . . . . . . . . . . . . . . . . . .

	Chapter 1. State of the art	8
	1.3.1

. Advanced parallel computing in material forming with cimlib. European Journal of Computational Mechanics/Revue Européenne de Mécanique Numérique, 18(7-8):669-694.

CHAPTER 1

State of the art

Contents 1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Lagrangian Description . . . . . . . . . . . . . . . . . . . . . . . 1.2.1 Spatial discretization . . . . . . . . . . . . . . . . . . . . . . . . 1.2.2 Temporal discretization . . . . . . . . . . . . . . . . . . . . . . 1.2.3 Contact . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.3.1 Contact search algorithms . . . . . . . . . . . . . . . 1.2.3.1.1 All-to-All algorithm . . . . . . . . . . . . . . 1.2.3.1.2 Bucket sort algorithm . . . . . . . . . . . . . 1.2.3.1.3 Spherical sorting algorithm . . . . . . . . . 1.2.3.2 Contact detection . . . . . . . . . . . . . . . . . . . . 1.2.3.2.1 Contact between rigid and deformable bodies. 1.2.3.2.2 Contact between deformable bodies. . . . . . 1.2.3.3 Parallel strategies in contact detection . . . . . . . . . 1.2.3.4 Contact resolution . . . . . . . . . . . . . . . . . . . . 1.2.3.4.1 Lagrange Multiplier Method . . . . . . . . . . 1.2.3.4.2 Penalty Technique . . . . . . . . . . . . . . . 1.2.3.4.3 Augmented Lagrangian Method . . . . . . . . 1.3 Eulerian Description . . . . . . . . . . . . . . . . . . . . . . . . 1.4 Arbitrary Lagrangian Eulerian Description . . . . . . . . . . . 1.4.1 Coupled Resolution . . . . . . . . . . . . . . . . . . . . . . . . 1.4.2 Split operator . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.5 Discussion and Conclusion . . . . . . . . . . . . . . . . . . . . . Résumé en Français . . . . . . . . . . . . . . . . . . . . . . . . . . . . Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  X, called the material variable, is used in Lagrangian description.x, called the spatial variable, is used in Eulerian description.

	1.2. Lagrangian Description		11
	ρ div(v) = 0 ∂v + divσ = 0 ∂t	dynamic equilibrium equation continuity equation	(1.6)
	behaviour law	+ boundary conditions
	Remarks and notations:

All quantities, for instance F , can be expressed as functions of X and t in Lagrangian description or as function of x and t in Eulerian description: F = F (X, t) value of F experienced by the initial particle of position X at instance t (Lagrangian description). F = F (x(X, t), t) value of F experienced by the particular of position x instantaneously (Eulerian description).

  Boman, R. (2010). Développement d'un formalisme Arbitraire Lagrangien et Eulérien tridimensionnel en dynamique implicit.Application aux opérations de mise en forme. PhD thesis, Université de Liège. [Bourago and Kukudzhanov, 2005] Bourago, N. G. and Kukudzhanov, V. N. (2005). A review of contact algorithms. Technical report, The Institute for problems in mechanics of RAS. [Bourel, 2006] Bourel, B. (2006). Calcul Multi-domaines et approches Multi-Échelles pour la simulation numérique de crashs automobiles. PhD thesis, LaMCoS -INSA de Lyon.

[START_REF] Bruchon | Using a signed distance function for the simulation of metal forming processes: Formulation of the contact condition and mesh adaptation. from a lagrangian approach to an eulerian approach[END_REF]

] Bruchon, J.,

[START_REF] Bruchon | Using a signed distance function for the simulation of metal forming processes: Formulation of the contact condition and mesh adaptation. from a lagrangian approach to an eulerian approach[END_REF]

. Using a signed distance function for the simulation of metal forming processes: Formulation of the contact condition and mesh adaptation. from a lagrangian approach to an eulerian approach. International Journal for Numerical Methods in Engineering, 78(8):980-1008.

[

Chabrand et al., 1998]

Chabrand, P., Dubois, F., and Raous, M. (1998). Various numerical methods for solving unilateral contact problems with friction. Mathematical and Computer Modelling, 28(4-8):97 -108. Recent Advances in Contact Mechanics. [Chamoret, 2002] Chamoret, D. (2002). Modélisation du contact: nouvelles approches numériques. PhD thesis, École Centrale de Lyon.

  Oliveira et al., 2008] Oliveira, M.,Alves, J., and Menezes, L. (2008). Algorithms and strategies for treatment of large deformation frictional contact in the numerical simulation of deep drawing process. Archives of Computational Methods in Engineering, 15(2):113-162.[START_REF] Osher | [END_REF][START_REF] Osher | Fronts propagating with curvature dependent speed: algorithms based on hamilton-jacobi formulations[END_REF][START_REF] Osher | Fronts propagating with curvature dependent speed: algorithms based on hamilton-jacobi formulations[END_REF]. Fronts propagating with curvature dependent speed: algorithms based on hamilton-jacobi formulations. Journal of Computational Physics, 79:12-49.

	Computer
	Methods in Applied Mechanics and Engineering, 58(2):227 -245.

[Martins and Oden, 1987

] Martins, J. and Oden, J.

(1987)

. Existence and uniqueness results for dynamic contact problems with nonlinear normal and friction interface laws. Nonlinear Analysis: Theory, Methods and Applications, 11(3):407 -428.

[

[START_REF] Fourment | Numerical formulations and algorithms for solving contact problems in metal forming simulation[END_REF]

Mocellin, K. (1999). Contribution à la simulation numérique tridimensionnelle du forgeage à chaud : Étude du contact et calcul multigrille. PhD thesis, École nationale supérieure des Mines de Paris. [[Paik et al., 2006] Paik, S. H., Moon, J. J., Kim, S. J., and Lee, M. (2006). Parallel performance of large scale impact simulations on linux cluster super computer. Computers & Structures, 84(10-11):732 -741. [Papadopoulos and Taylor, 1993] Papadopoulos, P. and Taylor, R. (1993). A simple algorithm for three-dimensional finite element analysis of contact problems. Computers & Structures, 46(6):1107 -1118. [Philippe, 2009] Philippe, S. (2009). Dévéloppement dd'un formulation Arbitraire LaLagrangien Eulérienne pour la simulation tridimensionnelle du laminage de produits plats. PhD thesis, École Nationale Supérieure des Mines de Paris. [Richter and Wick, 2010] Richter, T. and Wick, T. (2010). Finite elements for fluidstructure interaction in ale and fully eulerian coordinates. Computer Methods in Applied Mechanics and Engineering, 199(41-44):2633 -2642. [Roux and Garaud, 2009] Roux, F.-X. and Garaud, J.-D. (2009). Domain decomposition methodology with robin interface matching conditions for solving strongly coupled fluidstructure problems. International Journal for Multiscale Computational Engineering, 7(1). [Roux et al., 2005] Roux, F.-X., Magoulés, F., Series, L., and Boubendir, Y. (2005). Approximation of optimal interface boundary conditions for two-lagrange multiplier feti method. In Barth, T., Griebel, M., Keyes, D., Nieminen, R., Roose, D., Schlick, T., Kornhuber, R., Hoppe, R., Périaux, J., Pironneau, O., Widlund, O., and Xu, J., editors, Domain Decomposition Methods in Science and Engineering, volume 40 of Lecture Notes in Computational Science and Engineering, pages 283-290. Springer Berlin Heidelberg.

Table 2 .

 2 1: Parameters used in the flat bunch simulation

2.7 New additions and Numerical Applications 2.7.1 Test cases and proposed ameliorations 2.7.1.1 Optimal parameter choices

Table 2 .

 2 2: Comparison between the different definitions to compute the simplex size.

	76

  Ville, L., Silva, L., and Coupez, T. (2011). Convected level set method for the numerical simulation of fluid buckling. International Journal for Numerical Methods in Fluids, 66(3):324-344.

	CHAPTER 3
	Parallel Computing
	. 97
	3.3.3.1.1 Amdahl's law
	. . . . . . . . . . . . . . . . . . . . . . . . . 98
	3.3.3.1.2 Gustafson law
	. . . . . . . . . . . . . . . . . . . . . . . . . 98
	3.3.3.2 Communications between cores . . . . . . . . . . . . . 99
	3.3.3.3 Load balancing . . . . . . . . . . . . . . . . . . . . . . 100

Contents 3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91 3.2 Concepts and Terminology . . . . . . . . . . . . . . . . . . . . . 91 3.2.1 Flynn's Classical Taxonomy . . . . . . . . . . . . . . . . . . . . 92 3.2.2 Memory classification and access . . . . . . . . . . . . . . . . . 93 3.3 Performance of a Parallel Code . . . . . . . . . . . . . . . . . . 95 3.3.1 Hardware performance . . . . . . . . . . . . . . . . . . . . . . . 95 3.3.2 Software performance . . . . . . . . . . . . . . . . . . . . . . . 96 3.3.3 Limits and Costs of Parallel Programming . . . . . . . . . . . . 97 3.3.3.1 Serial fraction . . . . . . . . . . . . . . . . . . . . .

Table 3

 3 Ghz. Each core is characterized by a L1 memory of size 128 kb; L2 memory of size 512 kb with a speed of 2.3 Ghz and a L3 memory of size 12288 kb.

	.1

Table 3 .

 3 1:The simulations parameters used for speed up tests

  .2 shows the time needed to assemble and solve both Navier-Stokes and Stokes solvers.

	number of	Navier-Stokes	Stokes
	cores	Assembly Resolution Assembly Resolution
	1	352.9625	455.8391 177.2075	438.8229
	2	184.3498	235.1301	93.0619	241.1997
	4	94.4416	152.1518	48.5017	126.6885
	8	47.9751	86.8955	25.0818	72.7976
	16	24.5011	50.7569	13.0989	47.4226
	32	12.8818	30.1226	7.3847	30.2918
	Table 3.2: Time (in s) spent to assemble, solve the mechanical problem with two different solvers
	( Navier-Stokes vs Stokes ) using the 1.2.6 version of MPI

  Stokes and LevellerT solvers consume roughly 75% of the computational time whether in sequential or parallel simulations.

	number of	Stokes	LevellerT		Total Comput-
	cores	Assembly Resolution Assembly Resolution	ational Time
	1	177.2075	438.8229 143.9945	9.0872	987.0089
	2	93.0619	241.1997	77.2412	5.8230	544.9218
	4	48.5017	126.6885	39.5058	2.4788	305.4120
	8	25.0818	72.7976	19.7558	1.3304	221.6395
	16	13.0989	47.4226	10.0588	0.7003	218.0213
	32	7.3847	30.2918	5.6219	0.3675	126.0812
	Table 3.3: Time (in s) spent to assemble, solve the mechanical and transport solver and to
	accomplish the simulation using the 1.2.6 version of MPI

Table 3 .

 3

	number of	Stokes		LevellerT	Total Comput-
	cores	Assembly Resolution Assembly Resolution	ational Time
	1	177.5655	439.8508 142.0285	9.0568	984.0464
	2	94.7523	270.4347	75.7627	6.5459	568.8039
	4	50.7460	179.8544	39.9512	3.4741	351.5527
	8	25.4583	90.4282	19.8727	1.6952	232.6729
	16	12.9823	46.1396	10.0351	0.6594	220.7195
	32	6.9151	24.7715		5.3670	0.3563	119.2165
	number of	Written Generalized Total without Output Remaining
	cores	Output Files	strain rate	& initial partition	Time
	1	82.9321	37.1710		896.9749	128.4733
	2	49.9486	19.6609		512.8915	65.3959
	4	35.0572	9.7928		308.7822	34.7565
	8	63.6903	5.2728		156.9266	19.4722
	16	125.9468	2.6567		84.1534	14.3370
	32	59.3700	1.3556		47.8224	10.4125
	Table 3.6: Time (in s) spent to accomplish other relevant part of the code during simulation
	using the 1.4.3 version of MPI with the bind to core option

5: Time (in s) spent to assemble, solve the mechanical and transport solver and to accomplish the simulation using the 1.4.3 version of MPI with the bind to core option

Table 3 . 7 :

 37 Size of the Output files in parallel computing

	number of	Size of Output files
	cores	pvtu (Kb) vtu (Mb)
	1	1.5	442
	2	1.5	219
	4	1.6	108
	8	1.7	53
	16	2.1	27
	32	2.7	14

Table 3

 3

	number of	Stokes		LevellerT	Total Comput-
	cores	Assembly Resolution Assembly Resolution	ational Time
	1	178.1416	440.9533 146.1969	9.0786	904.0554
	2	93.0375	236.6903	78.0407	5.842	480.9263
	4	49.3545	143.6863	39.6819	2.9108	269.7911
	8	25.355	90.8500	20.0234	1.6955	156.3066
	16	13.3082	46.3393	10.1175	0.6729	81.5517
	32	6.9328	24.6166		5.3414	0.3682	47.9026
	number of	Written Generalized Total without Output Remaining
	cores	Output Files	strain rate	& initial partition	Time
	1	82.7806	37.6434		904.0554	129.6850
	2	43.2563	20.4194		480.9263	67.3158
	4	27.9493	9.9888		269.7911	34.1576
	8	19.3887	4.9998		156.3066	18.3827
	16	28.6615	2.6957		81.5517	11.1138
	32	61.0352	1.3785		47.9026	10.6436
	Table 3.9: Time (in s) spent to accomplish other relevant part of the code during simulation
	using the 1.4.3 version of MPI with the openmpi round robin option

.8: Time (in s) spent to assemble, solve the mechanical and transport solver and to accomplish the simulation using the 1.

4

.3 version of MPI with the opemmpi round robin option

Table 3 .

 3 11: Using the refined mesh: Time ( in s) spent to accomplish other relevant part of the code during simulation using the 1.4.3 version of MPI with the bind to core option

	number of	Stokes	LevellerT	Total Comput-
	cores	Assembly Resolution Assembly Resolution		ational Time
	1	3162.2112 25497.1773 2492.3585	427.7585		34844.7487
	2	1679.9621 13875.3077 1389.4449	275.8599		19218.6489
	4	850.2076 9572.5892 679.7318	129.1397		12424.3732
	8	411.8736 4127.7119 329.8567	57.4261		6291.4682
	16	211.6616 2100.1478 166.1422	25.2834		4556.1730
	32	106.0528 1021.6278	81.7341	10.4183		2597.8566
	64	57.3246	508.7427	42.2413	6.9786		1242.0377
	128	33.0700	287.6470	22.8955	5.8589		828.2126
	Table 3.10: Using the refined mesh: Time ( in s) spent to assemble, solve the mechanical and
	transport solver and to accomplish the simulation using the 1.2.6 version of MPI
	number of	Written Generalized Total without Output Remaining
	cores	Output Files	strain rate	& initial partition	Time
	1	1193.1272	629.5313		33616.8978	2037.3923
	2	833.5299	334.6115		18324.3467	1103.7721
	4	516.2075	165.2059		11812.8972	581.2289
	8	982.5606	80.9494		5217.2592	290.3909
	16	1712.9052	40.1335		2706.2488	203.0138
	32	1160.9286	21.0772		1325.9360	106.1030
	64	444.2012	10.3300		704.9010	89.6138
	128	291.3921	5.5548		417.4245	67.9531

Considering a computational domain Ω formed by a deformable body Ω d , a rigid tool Ω s and the air Ω a . The deformable body and the tool are represented by two distinct level set φ d and φ s . However the air is deduced by complementarity. We should point out that in a Eulerian framework, the contact is not defined explicitly. Therefore, the authors define the contact zone as follows:

Acknowledgements

This work has been a very special project. I am deeply greatful to my family and to all whose enthusiasm and energy transformed my dream into reality. I would like to express my sincere thanks :