Je Tiens Spécialement À Remercier

Vincent Robin

Alexandre Brosse

Rémi Lacroix

Maxime Dumas

François Josse

Kévin Chenegrin

Bastien Agard

Rémi Bertrand

Chaouki Tahri

Mathieu Girinon

Loïc Villié

Jian Cao

Hongfeng Ma

Wenxi Wang

Xuan Qi

Meng Wang

Kou

Keywords: Méthode des éléments finis, intégration nodale, thermomécanique, mouvement de matière, état stationnaire, repère mobile, soudage, laminage, usinage Finite element method, nodal integration technique, thermal-mechanical, material motion, steady-state, moving reference frame, welding, rolling, machining method . . . . . . . . . . . . . . .

J'ai envie d'adresser mes sincères remerciements à toutes les personnes qui ont contribué au succès de ma thèse. J'adresse mes premiers remerciements aux membres du jury pour avoir accepté d'évaluer mes travaux de thèse. Je remercie particulièrement à Jean-Baptiste LEBLOND pour m'avoir fait l'immense honneur de présider ce jury. Un grande merci à Alain RASSINEUX pour les échanges très fructueuse durant le congrè CSMA.

Je tiens tout particulièrement à remercier aux mes encadrants, Jean-Michel BERGHEAU et Jean-Christophe ROUX, qui m'ont soutenu, encouragés durant ma thèse. Je remercie infiniment pour leur patiences, leur disponibilités, et leur judicieux conseils. Je vous tiens ma plus profonde gratitude aux mes encadrants pour l'ensemble des discussions scientifiques et aussi pour les conseils que j'ai eu pour la vie et la vie professionnelle. Je vous remercie avec tous ma coeur. J'adresse mes sincères remerciements à tous les professeurs, intervenants et toutes les personnes qui par leurs paroles, leurs écrits, leurs conseils et leurs critiques ont guidé mes réflexions et ont accepté de me rencontrer et de répondre à mes questions durant mes recherches. Je tiens à remercier chaleureusement les équipes du LTDS. Je remercie particulièrement à Eric Feulvarch, Hédi Hamdi et Julien Sijobert pour leurs aides scientifiques. Je n'ai pas oublié de adresser mes sincères remerciements à Isabelle PLETTO, pour son côté maternel et ses soutiens qui ont été d'une grande aide. Merci Isabelle! J'en profite pour adresser un dernier mot de remerciement à Jérôme et Nicolas pour ses service informatique.

Résumé Simulation numérique des états stationnaires associés aux processus thermomécaniques

De nombreux procédés de fabrication thermomécanique comme le laminage, le soudage ou encore l'usinage mettent en jeu soit des sollicitations mobiles par rapport à la matière fixe, soit de la matière mobile par rapport à des sollicitations fixes. Dans tous les cas, après un régime transitoire en général assez court, les champs thermiques, métallurgiques et mécaniques associés à ces procédés atteignent un état stationnaire. La recherche de ces états stationnaires à l'aide de la méthode des éléments finis classique nécessite de mettre en oeuvre des modèles complexes et couteux où les sollicitations se déplacent par rapport à la matière (ou l'inverse).

La recherche directe des états stationnaires a fait l'objet de nombreux travaux de recherche ces trente dernières années. Des méthodes sont aujourd'hui disponibles et pour certaines sont proposées dans des codes de calcul du commerce. Ainsi, une option de calcul dite repère mobile proposée par différents auteurs est disponible dans le logiciel SY SW ELD T M . Cette méthode permet de calculer les états thermique, métallurgique et mécanique stationnaires associés à un procédé de soudage, en résolvant un problème de diffusion-convection en thermique et en intégrant, en mécanique, les équations constitutives du comportement du matériau le long des lignes de courant. Si cette méthode a été utilisée avec succès dans de nombreuses applications, elle présente néanmoins quelques limitations. Ainsi le maillage doit être structuré et la convergence des calculs est en général assez lente.

Nous proposons dans cette thèse de résoudre le problème mécanique dans un repère liè aux sollicitations, en nous appuyant sur une méthode de calcul par éléments finis reposant sur l'intégration nodale et la technique SCNI (Stabilized Conforming Numerical Integration). Cette méthode permet l'utilisation de maillages en tétraèdres (ou triangles en 2D) sans rencontrer de problème de verrouillage volumique résultant de l'incompressibilité plastique associée au critère de plasticité de von Mises. Plutôt que de rechercher directement l'état stationnaire, l'idée générale est ici de construire l'état stationnaire à vii partir d'une analyse transitoire en faisant entrer pas à pas la matière par l'amont et en la faisant sortir par l'aval d'un maillage fixe par rapport aux sollicitations et de taille limitée. L'état (quasi-)stationnaire n'est donc atteint qu'au bout d'un certain temps d'analyse. Les avantages de cette méthode résident dans:

1. L'utilisation d'un maillage libre en tétraèdres ou en triangles (au lieu des maillages structurés), 2. Un maillage raffiné uniquement dans la zone située au voisinage des sollicitations, 3. Une grande robustesse et notamment une bonne convergence des calculs mécaniques non linéaires liée à la résolution des états transitoires.

Après une introduction générale (Chapitre 1) et un état de l'art sur les méthodes existantes (Chapitre 2), nous présentons une approche de simulation du mouvement de matière dans le cadre de la méthode des éléments finis classique sur un problème de soudage (Chapitre 3). Nous y proposons également des conditions aux limites thermiques pertinentes pour calculer directement la distribution de températures en régime stationnaire.

La méthode des éléments finis reposant sur l'intégration nodale est ensuite décrite au Chapitre 4. Les avantages et inconvénients de la méthode sont discutés. La méthode est validée sur une application en grandes déformations élastoplastiques, un problème de flexion et une simulation thermomécanique de soudage.

La méthode des éléments finis reposant sur l'intégration nodale est alors développée pour prendre en compte un mouvement de matière (Chapitre 5). Trois types de mouvement sont considérés : en translation, circulaire et en hélice. Différentes méthodes de transport de champ sont abordées et discutées ainsi que le couplage thermomécanique. Des exemples d'application dans le domaine du laminage et du soudage pour différents mouvements de matière montrent l'efficacité de la méthode développée Des perspectives à ce travail sont proposées au Chapitre 6. Les perspectives envisagées visent d'une part à améliorer la méthode proposée et d'autre part, à développer la méthode pour simuler d'autres procédés. Une première application de la méthode à la simulation de la coupe orthogonale y est présentée.

The steady-state simulation in one increment has been the subject of much researches over the past thirty years. Methods are now available and some are integrated into calculation codes commercial. Thus, a so-called Moving Reference Frame method proposed by various authors is available in the SY SW ELD T M software. This method makes it possible to calculate the steady-state of thermal, metallurgical, and mechanical states associated with a welding process, by solving a thermal diffusion-convection problem in thermal-metallurgy and by integrating, in mechanics, the constitutive equations of the material along the streamline. Moreover, this method has been used successfully in many applications, it nevertheless has some limitations. Thus the mesh must be structured and the convergence of computations is generally quite slow.

In this thesis, we propose to solve the mechanical problem in a frame linked to the solicitations, by relying on a finite element calculation method based on nodal integration and the SCNI (Stabilized Conforming Numerical Integration) technique. This method allows the use of tetrahedron meshes (or 2D triangles) without encountering a locking problem resulting from the plastic incompressibility associated with the von Mises plasticity criterion. Rather than directly calculating the steady-state, the general idea here is to construct the steady-state from a transient analysis by bringing material step by step upstream and by making it exit downstream of a fixed mesh related to the solicitations ix and of the limited mesh size. The (pseudo-) steady-state is therefore only achieved after certain steps of analysis. The advantages of this method lie in:

1. The use of a tetrahedral or triangles mesh (instead of structured meshes), 2. A refined mesh is only needed in the area where the solicitation locates, 3. The method proposed shows the robustness and good convergence of nonlinear mechanical calculations because of the resolution of the transient states.

Apart from a general introduction (Chapter 1) and a state of the art on the existing methods (Chapter 2), we present an approach of simulation of the movement of material within the framework of the classical finite element method on a welding problem (Chapter 3). We also provide relevant thermal boundary conditions for directly calculating the steady-state of temperature distribution.

The finite element method based on the nodal integration technique is then described in Chapter 4. The advantages and disadvantages of the method are discussed. The nodalintegration-based finite element is validated by comparing its simulation results with classical finite element methods in large elastoplastic strains, a bending problem, and a thermomechanical simulation of welding.

The nodal-integration-based finite element is then developed and applied to simulate material motion (Chapter 5). Three types of movement are considered: translational, circular, and helical. Different methods of field transport are approached and discussed as well as thermomechanical coupling. Examples of applications in the field of rolling and welding for different movements of material show the efficiency of the developed method.

Perspectives for this work are presented in Chapter 6. The envisaged perspectives aim, on the one hand, to improve the proposed method and on the other hand, to develop the method to simulate other processes. A first application of the material motion method to the simulation of the orthogonal cut is presented there. The metal forming processes started from the Bronze Age in human history. These processes are mainly used to transform or change the metal parts to the desired shape by permanently mechanical deformation. Today, the high economic development and improvement of daily comfort level cannot be achieved without those processes.

Nowadays, competition among industrial companies demands developing products with increasingly reduced time and costs. The control of manufacturing industrial processes often determines the quality of products and thus constitutes an essential factor of success in this context.

Thermo-mechanical manufacturing processes, such as welding, rolling, machining, or additive manufacturing, are the most widely used metal forming processes. The numerical simulations of these processes are usually very time consuming due to the moving load, then a fine mesh and minor time step are used to capture the evolution of strain/stress. In the next subsections, since the numerical simulation of rolling, welding, and machining processes constitute the aims of thesis work, a short introduction of each process is presented.

Rolling processes

Rolling is a typical metal forming process, for example, the metal stock passes through one or more sets of pairs of rolls to reduce and uniformize the thickness of the sheet (figure-1.1). Rolling can be classified as either hot rolling or cold rolling, according to the temperature of the metal during the rolling process. Hot rolling induces coupled phenomena of recrystallization, phase transformation, and cooling that change the properties of the material. With the microstructural's modifications, metals could become stronger with improved weldability and formability, and these properties have a significant impact on its applications. According to the final form, the types of rolling processes can also be classified as ring rolling, roll bending, roll forming, profile rolling, and so on.

Rolling is a general operation for the metal produced. The rolling always takes place very early in the production chain, as soon as the metal becomes consistent enough to undergo plastic deformations. The rolling usually corresponds at an intermediate stage leading to a half-product that is subsequently intended for machining, cutting, or drawing, which means the quality of rolling will directly affect the final product's quality.

Even though the rolling process has been used for a very long time (the earliest rolling mills were slitting mills introduced in 1590). The paper of Buchmayr et al. [START_REF] Buchmayr | Future challenges in the steel industry and consequences for rolling plant technologies[END_REF] has provided an overview of recent steel market developments. It mentioned that the improvement of steel properties remains the main research field for many applications. For instance, an increase in the demand is expected regarding products with enhanced strength for light-weight design, improved toughness to ensure safety in case of earthquakes and fires, improved formability to overcome geometrical limits, etc. Therefore, developing a new efficient numerical model to model these physical problems is an interesting subject. In the rolling processes, the tools are usually fixed in space while the material will pass through between the tools.

Welding processes

Welding is one of the most common assembling processes to join materials by using high energy to melt the parts and welding material together, which can ensure good continuity of the assembly parts.

The welding processes often involve several physical coupled phenomena and interactions. The arc welding is chosen as an example figure-1.2, So first of all, the heat source providing the energy to melt the parts is an electric arc; the temperature will increase very quickly and create the melt pool. The Heat-affected zone (HAZ) is a volume surrounding the weld in which the temperature during the welding process exceeds the austenitization temperature. Combined with the stresses of uneven heating and cooling, the HAZ is generally where the high residual stress and strain appears after the welding process.

The mastering of welding processes is very important so that we can ensure the quality of the products. The liberty ship accident underlines the problems related to the welding process. To construct the ships in a more rapid and cheaper way, the liberty ships are constructed by using the electric arc welding technique which is a new efficient method but not a fully mastered technique at the time. Some ships suffered significant brittle fracture even broke in half without warning due to structural defects related to welding processes (figure-1.3).

The welding process is widely used in different domains, such as the nuclear industry, automobile, aeronautic. Welding processes lead to microstructural modifications, residual stresses, and distortions, the knowledge of which is of considerable interest. One thus wishes to be able to accurately predict the residual stress distribution in the weld region (molten zone and heat-affected zone), to estimate the structural integrity of a welded structure, or the distortions induced by welding, to control the feasibility of a welding process. Stresses and distortions mainly arise from temperature gradients and possible phase transformations occurring during cooling. Nowadays, mechanical engineers need to study and analyze the weldability of different materials and should also be able to optimize the welding processes parameters to make distortions and residual stresses within reasonable limits.

One should note that the temperature gradients induced by heat sources are very localized. The mechanical states are supposed to be stable once the heat source is far away. Therefore, Studying the welding process in a moving reference frame related to the solicitation could be useful and more efficient. Machining is a controlled material removal process, which can be used to transform the piece into the desired final shape. The principal machining processes can be classified as turning, drilling, and milling.

To respond to the demand of high integrity parts, industrials tend to improve the quality of their products. In the past decades, the experience mainly comes from the practice, therefore, the progress of machining technique is made slowly while the cost of development is considerable. Nowadays, the computation capacity and numerical algorithms are well developed, many researchers can model the machining process and study various factors such as the cutting velocity, thickness of the chip, lubrification condition, and so on by the numerical method. If we take an example of chip formation (figure-1.4), The main difficulty to control the machining process comes from the interactions that occur between cutting tools and workpiece. Although experimental set-up, with the use of a high-speed camera or other dispositions, can give us some information on the process, while some other information cannot be obtained, such as the strain/stress evolution or temperature at the separation zone. These limitations are the origin of the use of numerical simulation to fully understand machining processes.

The machining processes are very complicated processes because there is strong coupling between several physical phenomena. Moreover, in some fields (aeronautic, railways...), machined pieces require high precision control in terms of geometry, residual stresses, and distortions. The development of an efficient numerical method to model this complex phenomenon is one of the core subjects.

Introduction

Numerical simulation of thermo-mechanical processes

Numerical simulation is a great help in understanding and developing the thermomechanical processes. It allows us to analyze the influence of the various parameters and to access informations difficult or impossible to obtain by experimentation, such as temperature or stress evaluations in the core of the material during the process. The numerical simulation can reduce time, thus investigation costs in process development and mastering.

Nowadays, numerical simulation has become a decision-making tool. With the strongly increasing capabilities of computers, computational mechanics has been used more and more in engineering studies. Numerical simulation appears as a privileged tool to access the physical parameters characterizing the manufacturing process, not only the metal forming process. On the one hand, we can thus better understand the physical phenomena involved and optimize the operating conditions to achieve a better quality of the manufactured component. On the other hand, the material and mechanical modifications induced by the process must be taken into account during the design stage to reach the requirements in terms of mechanical strength and service life.

Most of the manufacturing processes, such as rolling, welding, and machining, involve the interaction between several physical phenomena such as heat transfer, metallurgy, and mechanics. As shown in Figure 1.5, the solid arrows represent couplings that are taken into account, and dotted arrows stand for interactions neglected to simplify computation models. Strong couplings are often performed to solve the thermo-metallurgical part and weak couplings are generally sufficient to predict mechanical state (displacement, residual stresses, distortions...). The strong coupling means that temperatures and phase proportions are solved in the same system of equations, while the mechanical simulation will take thermo-metallurgical results as input.

The numerical simulation of processes such as rolling, welding, or machining is 1.1. INTRODUCTION challenging because:

1. These processes involve moving loads (or heat source) in relation to solicitation leading to very high temperatures gradients or stress gradients.

2. The problem to solve is strongly nonlinear.

Therefore, the finite element simulation of this type of process generally involves meshes very fine along the solicitation path and very small time steps to capture temperature gradients and stress gradients. Computation models are often two-dimensional and limited to very localized areas near the interactions between the solicitations and the material in welding (or additive manufacturing), or between the heat source and the material in machining. Sometimes even though a reduced model is applied, it is still impossible to calculate the thermo-mechanical fields for a duration necessary to analyze some phenomena, such as the wear of the cutting tools in the machining process. Indeed, the steady-state can be reached after several seconds of cutting, while the explicit algorithms used to simulate the machining process can only simulate a few milliseconds.

It is obvious that solving the thermomechanical problems associated with this kind of process in a reference frame linked to the solicitations could be very interesting. Indeed, these simulations would involve smaller meshes that could be refined only near the solicitation zone thus enabling important cost reduction of the simulation. This is the aim of the present thesis work.

The thesis work is organized as follows:

1. In chapter 2, the state of art is presented. Different methods related to steady-state simulation are introduced.

2. The finite element method with material motion method is proposed in chapter 3. This material motion method is validated by comparing its simulation results with those of the Moving Reference Method.

3. In chapter4, the nodal-integration-based finite element method is developed. The numerical results obtained by the nodal-integration-based finite element are compared to those generated by classical simulations. The comparisons make the evidence of the efficiency of this new element proposed for various applications (volumetric locking test, bending test, thermal-elastic-plastic problems).

4. The nodal-integration-based finite element method with material motion simulation is introduced in chapter 5. As the nodal-integration-based finite element can calculate all the internal variables at nodes, the variables transfer procedure related to material motion simulation can perform in a more straight forward and easier way.

1. Introduction

State of the art

In the literature, different methods for steady-state simulations have been developed, such as (Updated) Lagrangian method, Eulerian formulation, Moving Reference Frame method, or Arbitrary Lagrangian-Eulerian (ALE) method, and other methods. In this section, we will try to discuss different methods reported in the literature, their advantages and limitations. Finally, a short conclusion will be addressed.

Lagrangian / updated Lagrangian formulation

The Lagrangian formulation is frequently used to solve the solid mechanic problem. The advantage of this formulation is that the material points follow the mesh nodes, and the history-dependent material properties related to the mechanical behavior of the solid can be solved in a natural and simple way. Therefore, the Lagrangian formulation is the primary choice for mechanical solid simulation, especially for thermal-mechanical simulations, which concern the couplings of metallurgy, thermal, and mechanical [START_REF] Bergheau | Editorial -numerical simulation of welding[END_REF][START_REF] Duranton | 3d modelling of multipass welding of a 316l stainless steel pipe[END_REF][START_REF] Leblond | Metallurgical and mechanical consequences of phase transformations in numerical simulations of welding processes[END_REF][START_REF] Siva Shanmugam | A transient finite element simulation of the temperature and bead profiles of t-joint laser welds[END_REF]. In order to take into account the large displacement/strain problem, the update Lagrange formulation can be applied by updating the coordinates of nodes at the beginning of each time step.

Despite its advantages, the Lagrangian formulation is not free of drawbacks. Firstly, as the loads move in space, the mesh should be fine enough along the trajectory of loads and the small-time step should be used to capture the rapid evolution of internal variables. This refine mesh sometimes leads to computational inefficiency. Secondly, oscillations appear in the Lagrangian description because of the interactions between time and space discretizations. In the case of large transformation, the mesh tangling may occur frequently, which could result in poor precision even divergence.

In order to make computation more efficient, Lindgren et al. [START_REF] Lindgren | Automatic remeshing for three-dimensional finite element simulation of welding[END_REF] and Bergheau et al. [START_REF] Boitout-F. Bergheau | Finite element simulation of processes involving moving heat sources. application to welding and surface treatment[END_REF] have developed and used the adaptive meshing technique for welding simulation(see . The local refine mesh will move with the load, as a high gradient of temperature or strain could appear in these areas. With this adaptive re-meshing technique, the mesh refinement is moving with the heat source, which could greatly decrease the size of the discrete problem to be solved. However, this technique requires special re-meshing tools and also physical quantities transfer procedures. And the re-meshing technique is usually limited to split the hexahedron into smaller hexahedrons. The penalty method or a Lagrange multiplier method should be applied to ensure the continuity for discontinuous mesh between refined and coarse areas. At the same time, a special algorithm should be developed to detect the refine mesh zone needed during the simulation. This adaptive re-meshing technique seems to be useful to reduce computational resources (CPU time, RAM) required for each time step while it is somehow complicated to use in the practice. Moreover, some strain and stress information at the fine mesh zone will be lost if the refined mesh disappears. The re-meshing procedure could also be used in simulations to avoid mesh tangling problems [START_REF] Im | Hexahedral mesh generation for remeshing in three-dimensional metal forming analyses[END_REF][START_REF] Léger | An updated lagrangian method with error estimation and adaptive remeshing for very large deformation elasticity problems: The threedimensional case[END_REF][START_REF] Rodrigues | Finite element remeshing in metal forming using hexahedral elements[END_REF]. Firstly, an efficient re-meshing tool is needed and one should note that the hexahedral element is not always possible for meshing all geometry, especially for geometry complex. Secondly, the data transfer from the old mesh to the new mesh should be performed after the re-meshing procedure. In reality, the efficiency and accuracy are often opposites to each other for variables transfer procedure, and choosing the suitable variables transfer algorithm requires lots of experience. The variable transfer method has been discussed a lot in the literature, especially in damage and fracture simulation [START_REF] Yang | Simulation of continuous damage and fracture in metal-forming processes with 3D mesh adaptive methodology[END_REF].

Eulerian method

The Lagrangian formulation needs the mesh refinement along the trajectory or adaptive meshing strategies with physical quantities transfer procedure, all these make Lagrangian description sometimes computational inefficient for steady-state computation. Then, Eulerian formulation with a spatial fixed grid seems to be more appropriate for steadystate simulation.

In the literature, the velocity-based Eulerian formulation has been widely applied for the steady-state simulations of the metal forming processes by the finite element method. An Eulerian formulation for steady-state simulations benefits from several advantages: the material flows through a frame of reference which is fixed in space, as we can see in figure 2.2. The mesh will not move with the material. Therefore, problems with grid distortion due to large deformation do not exist (in the case of rolling or machining). In the point of view of Eulerian formulation, gap entry and exit lines can thus be tracked precisely, avoiding the oscillations of forces, stresses, and strains from incremental models as nodes come into contact and finish contact (in case of rolling). Such spurious oscillations are particularly troublesome for simulation aiming to couple the thermo-elastic deformation via the contact stresses with the incremental Lagrangian formulation. For Eulerian formulation, the mesh refinement is only required in the areas where there is a high strain/stress gradient or contact zone, which could largely reduce the size of the discrete problem to be solved. As the mesh is fixed in space, Eulerian steady-state models sometimes require a free surface correction algorithm. In the literature, several Eulerian methods have been successfully implemented to model rolling processes in the past forty years. Nowadays, the steady-state of rolling simulation is a well-studied domain.

A three-dimensional steady-state deformation method is proposed by Mori et al. [START_REF] Mori | Simulation of three-dimensional deformation in rolling by the finite-element method[END_REF] and applied for plate rolling and edge rolling. This method is based on the rigid-plastic finite element method. However, the method presented works only when the deformation is uniform in the thickness direction. Then, this idea has been adopted by Yamada et al. [START_REF] Hamauzu | Three-dimensional analysis of mandrel rolling using rigid-plastic finite element method[END_REF], and the free surface is corrected on the basis of the newly traced streamlines.

Then, Bertand-Corsini et al. [START_REF] Chenot | A three-dimensional analysis of hot rolling with a steady-state thermo-mechanical approach[END_REF] and Lee et al. [START_REF] Dewhurst | Bulge prediction in steady-state bar rolling processes[END_REF] proposed a scheme for prediction of the free surface by adding the condition assuming that no mass flow can cross the free surface. The derived matrix equations are used to correct the free surface.

Kim et al. [START_REF] Kim | A new free surface scheme for analysis of plastic deformation in shape rolling[END_REF][START_REF] Kim | Steady state finite element simulation of bar rolling processes based on rigid-viscoplastic approach[END_REF] developed a steady-state simulation scheme for updating the mesh, by using a sectional sweeping technique based on the streamline tracing and three-dimensional contact algorithm. While there are no clear descriptions about how to precisely determine the contact area between the roll and workpiece, which is critical for the precise prediction of the deformation behavior.

In the Riper et al.'s [START_REF] Ripert | Méthode itérative de recherche de l'état stationnaire des procédés de mise en forme : application au laminage[END_REF] work, the viscoplastic model with surface correction and contact algorithm has been developed to model different steady-state of metal forming processes, which seems efficient. The thermal steady-state has also been solved in the Eulerian frame with streamline upwind Petrov-Galerkin (SUPG) technique or the leastsquares method without suffering thermal oscillations. However, this method cannot be used to solve the elastoplastic problems efficiently like cold rolling. One should know that the shrinkage of elastic deformation plays an important role in the final forme of the workpiece, especially for a thin plate.

In some works which consider steady-state simulations of metal forming processes [START_REF] Chenot | A three-dimensional analysis of hot rolling with a steady-state thermo-mechanical approach[END_REF][START_REF] Hamauzu | Three-dimensional analysis of mandrel rolling using rigid-plastic finite element method[END_REF][START_REF] Kim | A new free surface scheme for analysis of plastic deformation in shape rolling[END_REF][START_REF] Kim | Steady state finite element simulation of bar rolling processes based on rigid-viscoplastic approach[END_REF][START_REF] Ripert | Méthode itérative de recherche de l'état stationnaire des procédés de mise en forme : application au laminage[END_REF], the rigid-plastic or rigid-viscoplastic model without elasticity is commonly applied, because the elasticity can not only introduce numerical complicity of elastic computation based on velocity but also lead to convergence difficulties related to contact conditions (see. Other researchers studied the steady-state of a thermal-mechanical simulation with an Eulerian formulation. Different from the rolling simulation, the thermo-mechanical simulations, like laser surface treatment, require more information about residual stresses to estimate the impact on fatigue life of the structure. Besides, the behavior of the material used for thermo-mechanical processes is usually thermo-elastoplastic or thermo-elastoviscoplastic.

As the Eulerian formulation is a velocity-based formulation, calculating the stress from the velocity field is problematic. A solution is to directly establish the relationship between stress and strain rate by assuming that the material is incompressible. Since the strain rate can be derived directly from the velocity field, this method avoids complex deformation integration procedures by neglecting elastic deformations. However, this approach can not predict residual stress.

In order to take the elastic deformation and history-dependent variables into account, the elastic deformation can be computed by an 'elastic reanalysis' [START_REF] Zienkiewicz | Flow of solids during forming and extrusion: Some aspects of numerical solutions[END_REF] or by the 'initial stress-rate method' [START_REF] Dawson | Finite element analysis of steady-state elasto-visco-plastic flow by the initial stress-rate method[END_REF], and the history-dependent variables can be computed by integrating the evolution equations along their streamline [START_REF] Dawson | Finite element analysis of steady-state elasto-visco-plastic flow by the initial stress-rate method[END_REF] or the material response evolution equations can be solved simultaneously with the momentum balance equation via the Galerkin method [START_REF] Thompson | Some integration techniques for the analysis of viscoelastic flows[END_REF]. In the article of Agrawal et al [4], They note that, while both methods yield accurate results, the numerical implementation of the latter method is simpler than the streamline integration method. However, results obtained from the latter method exhibit oscillations due to the hyperbolic nature of the evolution equations. These numerical oscillations can be solved now by streamline upwind Petrov-Galerkin method.

Thompson and Yu [START_REF] Thompson | A flow formulation for rate equilibrium equations[END_REF] have proposed an Eulerian method using a rate equilibrium equation, which can solve pure elastic cases. However, for the elastoplastic condition, they simply reduced the maximum effective stress to certain yield stress instead of considering the plastic material flow.

Manitty et al. [START_REF] Maniatty | An eulerian elasto-viscoplastic formulation for steady-state forming processes[END_REF] presented an elasto-viscoplastic Eulerian method able to predict residual stresses. However, this method cannot handle low strain rate conditions. Qin and Michaleris [START_REF] Qin | Themo-elasto-viscoplastic modelling of friction stir welding[END_REF] developed a Galerkin Eulerian formulation with four unknown fields (velocity, stress, deformation gradient, and internal variable) to predict residual stresses of elasto-viscoplastic materials using Anand's model. A steady-state of welding and rolling simulations have been presented. The Eulerian elasto-viscoplastic model can accurately predict the residual stresses of a quasi-steady-state process as the lagrangian analysis. Due to too many degrees of freedom at each node (Figure 2.3), it mays leads to a significant increase in the size of the discrete problem to be solved, if the mesh model becomes bigger. Abdelkhalek et al. [START_REF] Abdelkhalek | Coupled and uncoupled approaches for thin cold rolled strip buckling prediction[END_REF][START_REF] Abdelkhalek | Strip flatness modeling including buckling phenomena during thin strip cold rolling[END_REF] has succeeded to apply the Asymptotic Numerical Method (ANM) shell element model to simulate the steady-state of bite bucking during clod rolling. The computed stress profiles simulated by ANM shell element consist well with experimental measurements.

Displacement-based reference frame formulation

In the the classical Eulerian formulation, we discuss the difficulty of residual stress computation and the two existing solutions for Eulerian formulation. Both solutions have advantages and limitations [4]. We can conclude that the Eulerian formulation seems to be more suitable for steady-state simulation. However, the main issue comes from the fact that Eulerian formulation use velocity as the primary field, which can not be used directly for residual stress computation.

Balagangadhar et al. [START_REF] Balagangadhar | and University of Illinois at Urbana-Champaign. A displacement-based reference frame formulation for steady-state thermo-elasto-plastic material processes[END_REF] introduced a displacement-based reference frame formulation for steady-state thermo-elastoplastic processes (see Figure 2.4). They adopt a Reference Frame kinematic description for simulating steady-state manufacturing processes. Displacement was chosen as the primary field instead of velocity, thereby obviating the need for free surface corrections and streamline integrations. As a result, the residual stress computation for history-dependent material becomes efficient, and this method was named mixed formulation. The hyperbolic material evolution equations are solved via the streamline upwind Petrov-Galerkin method [START_REF] Brooks | Streamline upwind/petrov-galerkin formulations for convection dominated flows with particular emphasis on the incompressible navier-stokes equations[END_REF] to eliminate the oscillations. Finally, a laser surfacing problem in 2 dimensions with temperature-independent material properties was studied. The results of displacement-based reference frame formulation give good agreement with those of the Lagrangian formulation.

The steady-state simulations of drawing and rolling processes have also been presented in a displacement-based reference frame by Balagangadhar et al. [START_REF] Balagangadhar | Design of large-deformation steady elastoplastic manufacturing processes. part i: a displacement-based reference frame formulation[END_REF]. These simulations show a good congruence with the Lagrangian formulation.

As we know, the temperature would change strongly in the welding process, and the temperature-dependant material properties should be used, which is not the case in Balagangadhar's simulation [START_REF] Balagangadhar | and University of Illinois at Urbana-Champaign. A displacement-based reference frame formulation for steady-state thermo-elasto-plastic material processes[END_REF]. The mixed formulation is reported to be difficult to converge for temperature-dependent material properties and significant relaxation was used according to Shanghvi et al. [START_REF] Shanghvi | Thermo-elasto-plastic finite element analysis of quasi-state processes in eulerian reference frames[END_REF].

In order to solve the problems related to temperature-dependent material properties, Shanghvi et al. [START_REF] Shanghvi | Thermo-elasto-plastic finite element analysis of quasi-state processes in eulerian reference frames[END_REF] has modified the governing equations of Reference [START_REF] Balagangadhar | and University of Illinois at Urbana-Champaign. A displacement-based reference frame formulation for steady-state thermo-elasto-plastic material processes[END_REF]. Temperaturedependent material properties have been implemented in these equations, and the convergence problems have been alleviated by applying smoothing functions proposed by Rajadhyaksha et al. [START_REF] Rajadhyaksha | Optimization of thermal processes using an eulerian formulation and application in laser surface hardening[END_REF]. The stress analysis has been realized in an Eulerian formulation with the mixed element [START_REF] Shanghvi | Thermo-elasto-plastic finite element analysis of quasi-state processes in eulerian reference frames[END_REF]. Moreover, for a laser surfacing steady-state simulation, they observe that the Eulerian formulation has a gain of 15 times faster than the Lagrangian formulation.

Moving Reference Frame method with integration of constitutive equations along the streamlines

This idea was originally proposed by Nguyen et al. [START_REF] Nguyen | Mouvement permanent d'une fissure en milieu elastoplastique[END_REF]. And it has been applied to study a moving crack in an elastoplastic medium; Dang Van et al. [38] and Maitournam [START_REF] Maitournam | Résolution numérique des problemes elastoplastiques stationnaires[END_REF] have applied this method for the evaluation of stresses induced in rails by the roll of wheels and other applications [START_REF] Inglebert | Sur un nouvel algorithme de calcul des structures elastoplastiques en regime stationnaire[END_REF][START_REF] Habibou | Modélisation par équations intégrales du frottement sur un demi-espace élasto-plastique[END_REF][START_REF] Maitournam | Résolution numérique des problemes elastoplastiques stationnaires[END_REF][START_REF] Hanus | Pressure rolling contact: Steady state flow analysis and comparison with experimental data[END_REF][START_REF] Nguyen | Mouvement permanent d'une fissure en milieu elastoplastique[END_REF][START_REF] Van | On a new method of evaluation of the inelastic state due to moving contacts[END_REF]. Hacquin et al. [START_REF] Hacquin | A steady state thermoelastoviscoplastic finite element model of rolling with coupled thermo-elastic roll deformation[END_REF] has formulated a steady-state thermo-elastoviscoplastic model for modeling rolling simulation, while this method becomes computational cost if the elastic is considered.

As a fully three-dimensional thermo-mechanical simulation is usually costly and timeconsuming, this method has been applied for 3D welding steady-state simulation. In order to circumvent this difficulty, the simulation presented in the article of Bergheau et al. [START_REF] Bergheau | Three-dimensional simulation of a laser surface treatment through steady state computation in the heat source's comoving frame[END_REF] is established on the direct evaluation of the stationary state in a comoving frame without calculating the transient intermediary states, which could significantly shorten the computation time because of reducing the number of time-steps to only one. Then a complete, thermal, metallurgical, and mechanical steady-state simulation of heat treatment has been presented, using the SY SW ELD T M software.

The thermo-mechanical steady-state simulation of the 3D welding process is then performed in a uncoupled manner: First, a steady-state thermal simulation in the heat source co-moving frame gives the temperature distributions corresponding to the end of the quasi-steady-state conditions. We have to deal with a diffusion-convection problem represented by equation (2.1) with appropriate boundary conditions.

ρW • gradH -divλgradT -Q = 0 (2.1)
In the above equation, ρ, H, and Q represent the density, the enthalpy of the material, and the heat source distribution respectively, and W is the speed of material points. For a translation movement, W is equal to -v which is a constant velocity of the heat source. The temperature-dependent material properties can be used in equation 2.1. Moreover, all the classical boundary conditions are available, such as prescribed temperature or heat flux.

The Galerkin method can be applied for such problems, while the finite element formulated could lead to spurious spatial oscillations due to high Peclet numbers P e = ρC∥W∥∆x λ , where the ∆x denotes the mesh size, and typically P e are greater than 2 with first-order elements. A Petrov-Galerkin variational formulation with special discontinuous test functions has been suggested by Hughes et al. [START_REF] Brooks | Streamline upwind/petrov-galerkin formulations for convection dominated flows with particular emphasis on the incompressible navier-stokes equations[END_REF] to avoid spurious spatial oscillations, and it works very well. As we can see in figure 2.5, the 3D mesh of the Moving Reference Frame (MRF) method should be generated by the translation of a cross-section in the welding direction. Therefore, the streamlines are parallel to the trajectory of the welding source. Before starting the simulation, a pre-processor will determine the sequences of Gauss points, the number of the preceding Gauss point on the streamline, and the element that contains these Gauss points. The pseudo time steps ∆t denotes the time needed for a material 2. State of the art particle to move from the previous Gauss point to the Gauss point studied, which depends on the distance between two successive Gauss points and the speed of the welding source.

In order to calculate the phase proportion, the step by step computation procedure starts from the first Gauss point (see figure 2.5) as the phase proportion corresponds to the initial state. The phase proportion of the second Gauss point will be calculated by using the phase proportion and temperature of the first Gauss point, the temperature of the second Gauss point, and the pseudo time steps given by the pre-processor. The phase proportion of 3rd, 4th, ... Gauss point can be computed in order. For the Gauss point studied, all the information is stored in the Gauss points at upstream, and we can call these points as the preceding points.

As the element has a mixture of phases, the latent heat effects can still be solved by the classical subroutines to compute the finite element vectors and matrices. The density and the enthalpy are calculated in the following manner:

ρ = N p i ρ i (2.
2)

H = N p i H i (2.3)
where N , p i , ρ i and H i are the number of phase, the proportion, the density and the enthalpy of phase i, respectively. These parameters could be the function of temperature. The equation (2.1) then becomes:

ρCW • gradθ -divλgradθ -Q + ρW • phases H i gradp i = 0 (2.4) with: C = N p i dH i dθ and W • gradp i = ∆p i ∆t
, in which ∆p i is the increment of the proportion of phase i between two successive Gauss points on a streamline and ∆t represents the pseudo time step. Therefore:

ρW • N H i gradp i = ρ N H i ∆p i ∆t (2.5)
The steady-state mechanical simulation is solved in a moving reference frame associated with the welding source. The mechanical analysis is expressed in classical form:

div σ + ρf = 0 (2.6)
To integrate constitutive equations of material, a similar procedure as the metallurgical transformation presented above is used. At each time step, the state variables at time t + ∆t is calculated from the state variables at time t (i.e. stress components, plastic strains, hardening parameters, etc.). At the beginning of each time step, The state variables of each Gauss point will be replaced by the state variables of his preceding Gauss point. Typically, the strain and stress of a Gauss point depend on all the preceding Gauss points in the streamline. In order to avoid a full non-zeros stiffness matrix, only the nearest preceding gauss point is used.

In conclusion, this procedure has been successfully applied for the simulation of surface treatment or welding processes. The results presented by Leblond et al. [START_REF] Leblond | Metallurgical and mechanical consequences of phase transformations in numerical simulations of welding processes[END_REF] and Bergheau et al. [START_REF] Bergheau | Three-dimensional simulation of a laser surface treatment through steady state computation in the heat source's comoving frame[END_REF] show good agreements with experimental results. Moreover, this method works well for classical elastoplastic and elastic-viscoplastic behavior.

The MRF method requires a regular mesh that should be obtained by translation of a cross-section as mentioned so that integration of the constitutive equation along streamline can be achieved directly with maximal efficacy. This method is dedicated to the small deformations and small displacements hypothesis.

Arbitary Lagrangian-Eulerian method

Noh et al. [START_REF] Noh | Cel: A time-dependent, two-space dimensional,coupled eulerianlagrangian code[END_REF] and Trulio et al. [START_REF] Trulio | Theory and structure of the afton codes[END_REF] originally developed the Arbitrary Lagrangian-Eulerian method. The Arbitrary Lagrangian-Eulerian formulation has been applied for fluid-structure interaction problems by Donea et al. [START_REF] Donea | An arbitrary lagrangian-eulerian finite element method for transient dynamic fluid-structure interactions[END_REF]. Then, ALE formulation was further extended to solve solid mechanics by various researchers Liu et al. [START_REF] Wing | An arbitrary lagrangianeulerian finite element method for path-dependent materials[END_REF], Benson [START_REF] Benson | Computational methods in lagrangian and eulerian hydrocodes[END_REF][START_REF] Benson | An efficient, accurate, simple ale method for nonlinear finite element programs[END_REF], Huétink et al. [START_REF] Han Huetink | Progress in mixed eulerian-lagrangian finite element simulation of forming processes[END_REF] etc. Huetink et al. [START_REF] Han Huetink | Progress in mixed eulerian-lagrangian finite element simulation of forming processes[END_REF] gives a review of the arbitrary Eulerian-Lagrangian finite element method for simulation of forming processes. The ALE formulation allows mesh points to move independently from material particles to reduce the possibilities of numerical difficulties due to large element distortions (see ). The free surface can be taken into account by adapting nodes in the Lagrangian way. The ALE formulation has been widely applied for metal forming simulations [START_REF] Jose Aymone | Simulation of 3d metal-forming using an arbitrary lagrangian-eulerian finite element method[END_REF][START_REF] Bui | Application of the arbitrary lagrangian eulerian formulation to the numerical simulation of cold roll forming process[END_REF][START_REF] Ponthot | Continuous roll forming simulation using arbitrary lagrangian eulerian formalism[END_REF][START_REF] Crutzen | Lagrangian and arbitrary lagrangian eulerian simulations of complex roll-forming processes[END_REF][START_REF] Aymone | Mesh motion techniques for the ale formulation in 3d large deformation problems[END_REF][START_REF] Ponthot | Enhanced ale data transfer strategy for explicit and implicit thermomechanical simulations of high-speed processes[END_REF][START_REF] Ponthot | Finite element simulation of lubricated contact in rolling using the arbitrary lagrangian-eulerian formulation[END_REF] etc. Each time step of ALE formulation can be generally split into two steps:

1. First step is a classical Lagrangian step, which means the mesh points are moving with material particles. The equilibrated state will be calculated in the Lagrangian configuration.

2. Second step is an Eulerian step. The mesh tool will define a new mesh by moving the mesh points of the old mesh. Then the data transfer from the old mesh to the new one will be performed. Compared to the re-meshing step in Update Lagrangian formulation, the ALE formulation keeps the same topology of the mesh, such as the number of elements and the element connectivity, during the entire simulation. Therefore, the transfer of physical quantities is generally more accurate because the convection technique can be applied in ALE simulations.

Wisselink et al. [START_REF] Wisselink | 3d fem simulation of stationary metal forming processes with applications to slitting and rolling[END_REF] presents the 3D FEM simulation of modeling the steady-state of slitting and shape rolling processes. Crutzen et al. [START_REF] Crutzen | Lagrangian and arbitrary lagrangian eulerian simulations of complex roll-forming processes[END_REF] apply ALE formalism to compute the stationary state for complex roll-forming simulation. A classical Lagrangian approach has also been used as a reference. They found that ALE and Lagrangian results are generally quite similar.

Mixed Eulerian-Lagrangian method

Unlike the ALE formulation, the Mixed Eulerian-Lagrangian coordinate system would employ an Eulerian coordinate in the longitudinal direction and Lagrangian coordinates in the remaining two directions to account for local deformation effects directly. This model requires a systematic transformation of the weak formulation of the governing equations and of its quantities involved in the intermediary reference frame. The velocity field, mean stress, and the contact pressure is determined by the mixed variational formulation [START_REF] Bathe | Finite element procedures in engineering analysis[END_REF][START_REF] Altan T Kobayashi | Metal forming and the finite-element methods[END_REF][START_REF] Taylor | The finite element method[END_REF].

According to Vetyukov et al. [START_REF] Yu | Mixed eulerian-lagrangian description in materials processing: deformation of a metal sheet in a rolling mill[END_REF], compared with ALE formulation, one important distinctive feature is the use of intrinsic strains, when the undeformed configuration of the body is incompatible because of the variable velocity of generation of material. Another difference is that the mixed Eulerian-Lagrangian method involves a two-stage mapping from the reference to the actual state, as shown in figure 2.7. Synka et al. [START_REF] Synka | A novel mixed eulerian-lagrangian finite-element method for steady-state hot rolling processes[END_REF] has applied the Mixed formulation for steady-state hot rolling simulation with a rigid-viscoplastic material behavior. As we know, the elasticity can change the final forme for thin plate processes. Therefore, it's still necessary to develop new methods that can take elasticity into consideration.

Conclusion

To conclude the state of art, the thermal-mechanical processes are typically simulated by Lagrangian formulation. Lagrangian formalism's benefits are that there is no free surface problem, and this formulation can take history-dependent material into account in a most straightforward manner. The price for the Lagrangian formulation is that the mesh should be fine along the load trajectory in space. However, in the case of large deformation, the mesh distortion can lead to less accurate results or even divergence in the numerical simulation.

In order to overcome the mesh distortion problem in large deformation problems, the updated Lagrangian formulation with the re-mesh technique has been proposed and applied for 3D forming processes. The updated Lagrangian formulation can take large displacement and deformation into consideration, and the re-mesh procedure can avoid mesh distortion. However, the re-mesh technique for the 3D hexahedrons for all geometry is still not always possible, and there is no robust mesh-generator for hexahedrons at the moment. At the same time, a data transfer strategy is also needed. Using an updated Lagrange formulation with re-mesh technique is less attractive if a large number of re-mesh steps are needed because of localized deformations.

In the case of welding simulation, Several authors have proposed to use adaptive re-meshing technique to reduce the size of meshes. The adaptive re-meshing technique consists of cutting a mesh locally according to the position and nature of the source. Since the thermal loads move on the structure, it is necessary to refine an initial coarse mesh at several steps corresponding to the position of heat source. This technique works well for the structure mesh (such as mesh is generated by translation or rotation of a section), in which the data transfer and refinement would not lead to the particular difficulties.

Among these simulations carried out by Lagrangian formulation, one can observe that the temperature, phase proportion, and residual stresses could reach steady-state after specific steps of transient analysis. The Eulerian formulation could be suitable for steady-state of rolling as the steady-state boundarys are knwon. The material flows through a frame fixed in space. Therefore, problems with mesh distortion due to large deformations do not exist in such a formulation. What's more, mesh refinement is only required around the load. However, in general, material boundaries are not equal to the element edges. Thus special procedures are needed to solve free surfaces problem or contact conditions and also to handle history-dependent material properties.

Mori et al. [START_REF] Mori | Simulation of three-dimensional deformation in rolling by the finite-element method[END_REF] solved free surface problems using an iterative surface correction method. This method has been modified by Kim et al. [START_REF] Kim | A new free surface scheme for analysis of plastic deformation in shape rolling[END_REF][START_REF] Kim | Steady state finite element simulation of bar rolling processes based on rigid-viscoplastic approach[END_REF] to take contact condition into consideration. Galerkin method [START_REF] Thompson | Some integration techniques for the analysis of viscoelastic flows[END_REF] and streamlines integration technique [START_REF] Dawson | Finite element analysis of steady-state elasto-visco-plastic flow by the initial stress-rate method[END_REF] has been applied for strain and stress histories by using the visco-plastic model. And streamline upwind Petrov-Galerkin (SUPG) technique [Balagangadhar et al. [8]] can remove numerical oscillations. Shanghvi et al. [START_REF] Shanghvi | Thermo-elasto-plastic finite element analysis of quasi-state processes in eulerian reference frames[END_REF] has noted that the streamline of Balagangadhar has difficulty in convergence, and significant relaxation was used. Therefore, they try to use the SUPG method with the mixed element to remedy that difficulty, and convergence problems have been alleviated by using smoothing functions. Their method has been applied for welding simulation, and they obtain a gain of 15 for computation time than classical Lagrange formulation. Bergheau et al. [START_REF] Bergheau | Three-dimensional simulation of a laser surface treatment through steady state computation in the heat source's comoving frame[END_REF] has also proposed the MRF method for steady-state simulation. This method gives good prediction results and is computationally efficient, while the mesh should be obtained by translation of a section in the flow direction and the MRF method doesn't converge for large transformation problems.

In this thesis work, we would like to propose a general moving reference frame method related to the solicitation for steady-state simulation. It would present several advantages:

1. A Nodal-integration-based finite element was developed, and triangle/tetrahedral mesh can be used without suffering volumetric locking.

2. A new mixed Eulerian-Lagrangian formulation will be proposed. An eulerian coordinate system is applied in solicitation's moving direction, while the local deformation can be taken into account directly in the Lagrangian manner. Therefore, the free surface problem doesn't exist.

3. All classical material behaviors are considered, such as viscoplastic, elastoplastic, elastic-viscoplastic, thermo-elastoplastic, and so on.

4. The mesh refinement would be required only nearby solicitations, which can greatly reduce the size of the discrete problem.

5. Compared with Lagrangian step by step simulation, the solicitation is fixed in space.

The contact conditions between tools and workpiece could be largely simplified and become more stable, which presents advantages in terms of convergence and a larger time step can be used.

Chapter 3

Finite element simulation with material motion

Introduction

In the first chapter, different methods in the literature have been presented, their advantages, and their shortcomings. The MRF (Moving Reference Frame [START_REF] Bergheau | Three-dimensional simulation of a laser surface treatment through steady state computation in the heat source's comoving frame[END_REF]) method seems very efficient and accurate for steady-state welding simulation.

The MRF method has been developed for welding simulations and solves the steadystate problem in a moving frame associated with the heat source. This method has been implemented in SYSWELD T M and can deal with a uniform translatory, circular, or helical motion. The latter option is interesting to simulate coating processes in cylindrical parts. However, for a translatory or helical motion, steady-state conditions would only exist for components presenting an "infinite" length in the welding direction. In practice, experience shows that quasi-steady-state conditions are fulfilled for a large part of the heat source path even if a circular motion in only one turn. Therefore, in this chapter, a new Moving Reference Frame model related to the heat source will be presented. With this model, we will calculate the steady-state by simulating the motion of material in a reference frame linked to the heat source (see figure 3.1). Our objective is to propose a new model that can give good results as the MRF method; at the same time, as the steady-state can be achieved by simulating transient states, this new model can remedy the difficulty in convergence that may occur in the MRF method. Finally, we will present a validation test at the end of this chapter. In order to clarify this idea, we consider a case of welding or rolling, in which the material moves in a uniform translation and is submitted to the local solicitations fixed in this moving frame. The mesh is fixed in this reference frame. Upstream and downstream of this area, mechanical fields are assumed to remain constant. From a physical point of view, the material flow is supposed to go through the computation frame as indicated in figure 3.2. This proposed method seems quite similar to the arbitrary Lagrangian-Eulerian (ALE) methods while the departing is quite different. There is no re-mesh procedure and the material point will move with the mesh. Rather than directly calculating the steadystate by integrating the materials' constitutive law along the streamline, we propose to simulate material motion in the reference frame, which has more advantages in terms of convergency.

Similar to the MRF method, the simulation is performed in two steps: first, a steadystate thermal simulation; Secondly, the steady-state mechanical simulation. The second step uses the thermal results of the first step, but the metallurgical phenomena as well as the heat generated by plastic deformation are disregarded.

Thermal steady-state computation

The steady-state exists with the condition that the model should be long enough in welding direction. Then the temperature distribution in the moving frame related to the heat source becomes time-independent. Therefore, in order to validate the steady-state thermal model, a transient thermal simulation has been prepared in advance, which is presented as the reference.

We start with the heat transfers in a solid medium occupying a domain Ω. The governing equation for transient heat transfer analysis is presented by the following:

ρ dH dt -div(λ grad T ) -Q = 0, in Ω (3.1)
λ grad T • n = q(T , t), with prescribed heat fluxes ∂Ω q (3.2)

T = T p (t), with prescribed temperatures ∂Ω T (3.3) ∂Ω = ∂Ω q ∪ ∂Ω T (3.4)
where ρ is density,

H mass enthalpy, λ thermal conductivity, T temperature, Q internal heat source,
n normal to the surface of ∂Ω, q surface heat flux which may depend on temperature and time to represent convection or radiation phenomena on the surface, ∂Ω q a part surface of ∂Ω on which a surface density of heat flux is imposed, T p temperature imposed, ∂Ω T a part surface of ∂Ω where the temperature is imposed.

The heat source can be represented by an internal volume heat source Q(x, y, z) or by a surface heat flux density q(x, y), which depends on the manufacturing processes. The thermal steady-state problem can be solved by a convection-diffusion problem with appropriate boundary conditions, at the same time, Petrov-Galerkin variational formulation with special discontinuous test functions has been applied to avoid spurious spatial oscillations of the solution for high Peclet numbers. Finally, the heat equation of steady-state in a moving frame associated with the heat source, takes the following form as we have presented in chapter 1, and W is a steady-state velocity (simply equal to -v s ):

ρW • gradH -divλgradT -Q = 0 (3.5)

How to choose boundary conditions

For the boundary conditions on the upstream and downstream faces of the mesh, one could think to prescribe the thermal conditions existing at infinity upstream and downstream. These conditions are:

1. Temperature imposed to room temperature which is supposed to be the temperature of the metal before welding:T = T 0 .

2. Conduction heat flux equal to zero: λ ∂T ∂n = 0.

But sometimes imposing both conditions is not possible. We, therefore, have to choose which condition to impose on each face. We suggest:

1. Upstream : Temperature imposed to room temperature: T = T 0 .

2. Downstream : Conduction heat flux equal to zero: λ ∂T ∂n = 0.

But to strongly enforce the upstream condition on the upstream face, special care should be paid to determine the distance between the heat source and the upstream face of the mesh. The distance chosen should be long enough to satisfy also λ ∂T ∂n ≈ 0 on the face and so, ensure that the temperature imposed on the upstream face will not affect the steady-state temperature distribution given by the model.

To impose the downstream condition on the downstream face, one should consider a mesh long enough so that the body can cool down to room temperature. One can notice that if only the condition of null conduction heat flux is imposed, nothing ensures that the temperature will reach the room temperature on the downstream face.

In order to limit the difficulties discussed above, we propose to define new boundary conditions to the mesh. These new conditions are based on a semi-analytical solution of the heat equation outside the mesh. We suppose that outside the mesh, the temperature distribution is uniaxial. The temperature is therefore supposed to be uniform in any section orthogonal to the weld path. We also supposed that the thermal properties are temperature independent.

Let's so consider the following model for the area beyond the upstream face (as we see in figure 3.6). One can note that for this situation, one has:

W = -W e x with W > 0 n = e x dT dn = dT dx (x = 0) (3.6)
The temperature distribution T (x) for x ∈ [0; +∞[ is governed by the following equation:

λ d 2 T dx 2 + ρCW dT dx + K (T 0 -T ) = 0 (3.7)
In the equation above, the term K(T 0 -T ) model the cooling condition with the ambient air through the outside surface of the model. The volumic heat exchange coefficient K is defined from the surface heat exchange coefficient H applied on the surface of the part by:

K = H * P S (3.8)
with P , the perimeter and S, the surface of a section.

The solution of equation 3.7 above is given by:

T (x) = T 0 + Ae (α-β)x + Be (α+β)x (3.9) With α = - ρCW 2λ and β = (ρCW ) 2 + 4Kλ 2λ = α 2 + K λ
and where A and B are 2 constants.

One can note that αβ < 0 and α + β > 0.

Because we should have lim x→+∞ T (x) = T 0 , it follows B = 0 and so

T (x) = T 0 + Ae (a-β)x (3.10)
From equation 3.10, we deduce:

dT dx = A(α -β)e (α-β)x (3.11)
And then

λ dT dn = λ dT dx (x = 0) = λA(α -β) (3.12)
Finally, we obtain the boundary condition to apply on the upstream face as following:

λ dτ dn = λ(β -α) (T 0 -T (x = 0)) (3.13)
Let's now consider the following model for the volume beyond the downstream face (figure 3.7). One can note that for this situation, one has:

W = -W e x with W > 0 n = -e x dT dn = - dT dx (x = 0) (3.14)
The temperature distribution T (x) for x ∈ [-∞; 0[ is still governed by the following equation 3.7:

The solution of equation 3.7 is again given by equation 3.9.

Because we should now have lim

x→-∞ T (x) = T 0 , it follows A = 0 and so

T (x) = T 0 + Be (α+β)x (3.15)
From equation 3.15, we deduce :

dT dx = B(α + β)e (α+β)x (3.16)
And then

λ dT dn = -λ dT dx (x = 0) = -λB(α + β) (3.17) λ dT dn = λ(α + β) (T 0 -T (x = 0)) (3.18)
Which is the boundary condition to apply on the downstream face.

One can notice that if the mesh used is long enough to satisfy λ ∂T ∂n ≈ 0 on the upstream or downstream face, we concurrently have T ≈ T 0 .

One should also notice that :

upstream: λ(β -α) = 1 2 ρCW       1 + 1 + K λα 2       = 1 2 ρCW        1 + 1 + 4λK (ρCW ) 2        downstream: λ(α + β) = 1 2 ρCW       1 + K λα 2 -1       = 1 2 ρCW         1 + 4λK (ρCW ) 2 -1 (3.19)
And so that λ(βα) > λ(α + β).

For air-cooling conditions, both convection and radiation losses are taken into account. The same coefficient as the transient simulation is set up. This study has identified that boundary conditions can make a significant difference in temperature computation. The steady-state temperature distribution will be simulated with the new boundary condition for the following tests.

Mechanical steady-state computation

Simulation of material motion

Once the thermal steady-state is validated, two mechanical steady-state simulations will be carried out based on the same thermal result. The first simulation is calculated by the MRF method [START_REF] Bergheau | Three-dimensional simulation of a laser surface treatment through steady state computation in the heat source's comoving frame[END_REF], which is the reference simulation. Another simulation is calculated with the material motion method that we proposed.

Figure 3.12 shows the diagram of the steady-state welding simulation.

Here, the material motion simulation diagram is presented in figure 3.13. We firstly have a computation frame in space, and we can imagine the material flow passes through this computation frame.

The mechanical steady-state simulation is carried out as follows:

1 With the thermal steady-state simulation results, the mechanical equilibrium is supposed to be solved at an instant t.

2 In order to calculate the mechanical equilibrium at time t+∆t, the mesh is translated for a distance of d = W * δt in space, which is supposed to be the length of a layer of the mesh. This is a new model A.

3 Internal variables (stresses, strains, displacement, and so on) transfer procedure are performed between two models: one is the model A of instant t + δt, another is the model B with zero strain and stress. The mechanical fields of model A will be For simplicity of the projection procedure, one should note that the length of each layer is identical and the new material layer arrives with zero strain/stress. These steps will be repeated until the mechanical steady-state is observed.

Application to 3D welding simulation

In order to validate this new method, a 3D welding simulation test has been carried out. The previous thermal steady-state result (B) (figure 3.10) has been adopted for both methods.

For the mechanical boundary condition (figure 3.14), there is a symmetry plane as mentioned in thermal simulation. A front red point has been fixed to avoid X, Y , Z translation. A behind violet point can prevent rotation around X axis.

A standard ASCII 316L thermo-elastoplastic mechanical properties from the database of SY SW ELD T M used for thermal-mechanical simulation.

As shown in figure 3.15, longitudinal stress distribution contours are very similar, which indicates the method proposed works very well.

Conclusion

We have first proposed and validated a new thermal boundary condition to obtain a correct temperature prediction even though the steady-state condition is not fulfilled. From the graph above (see figure 3.11), we can see that the new boundary condition gives a great agreement with a reference solution.

For the mechanical aspects, a new steady-state simulation method is presented. In place of calculating steady-state directly by integrating the constitutive equation along the streamline, simulating material motion in a moving reference frame associated with the heat source is proposed and a steady-state will be achieved after some computation steps. Compared with the MRF method, the method proposed is easier to converge but could be more time-consuming because of simulating transient states.

Finally, a 3D welding steady-state simulation has been presented. A comparison of steady-state mechanical simulations has been made. The very similar residual stress contours prove that the method proposed works well as the MRF method for steady-state simulation.

However, one should note that there are several limitations:

1 The projection presented in figure 3.13 should be performed from Gauss points (model A) to Gauss points (model B) directly.

Chapter 4

The nodal integration based finite element method

Introduction

In the last chapter, we discussed solving welding steady-state in a moving reference frame with the heat source. As we present that there is a mechanical field transfer procedures in order to simulate material motion in figure 3.13. In practice, the classical finite element method (FEM) will store physical state variables both at nodes and integration points. The transfer of physical quantities is usually performed with interpolation/extrapolation between nodes and Gauss points, which may lead to lack of precision. In the literature, different field transfer methods have been investigated to deal with different analyses and to satisfy different requirements. Generally, these transfer processes could be classified into three groups [START_REF] Yang | Simulation of continuous damage and fracture in metal-forming processes with 3D mesh adaptive methodology[END_REF]:

1. The mechanical result field at integration point of old mesh (IP old ) will be extrapolated to the nodal points of old mesh (N P old ), and the values at the integration points of the new mesh (IP new ) are computed by interpolation using shape fonction of the old mesh. We can express this process as:

IP old → N P old → IP new .
The value at nodal points of new mesh (N P new ) can be obtained after the discrete extrapolation.

2. Similiar as the first group, instead of computing the values at the integration points of the new mesh (N P old → IP new ), the nodal values of new mesh will be calculated first. The field transfer will be expressed as: IP old → N P old → N P new → IP new .

3. In the third group, mechanical field at integration point in old mesh are directly transferred to integration points of new mesh (IP old → IP new ).

More explications can be found in the works of Yang and Rassineux et al. [START_REF] Yang | Simulation of continuous damage and fracture in metal-forming processes with 3D mesh adaptive methodology[END_REF][START_REF] Yang | A 3d hadaptive local remeshing technique for simulating the initiation and propagation of cracks in ductile materials[END_REF]. These interpolations or extrapolations can be performed continuously by averaging, weighted averaging, or finite element shape function or least square method.

In our work, we focus on filed transfer processes that should be efficient and selfconsistency. In order to simplify the computational complexity of these field transfer processes related to interpolation or extrapolation between the values at nodal points and at integration points, the nodal-integration-based finite element seems an optimal solution as the integration points coincide with nodes. Therefore, the field transfer of the values at nodal points and at integration points can be handled at one time. Moreover, the nodal integration technique does not suffer from locking phenomena due to plastic incompressibility in metal forming simulations.

Finite element simulation of elastoplastic or viscoplastic problems based on the von Mises criterion leads to difficulties associated with the plastic incompressibility imposed by the constitutive equations. This nearly-incompressibility condition is imposed at all integration points where the constitutive equations are solved. If the number of integration points is too large compared to the number of degrees of freedom (all the displacements of the nodes), the number of degrees of freedom becomes insufficient to solve both the incompressibility equations and the equilibrium equations. Then, this leads to a volumetric locking phenomenon. In order to check if the system is able to perform well in the case of nearly-incompressible media, Nagtegaal et al. [START_REF] Rice | On numerically accurate finite element solutions in the fully plastic range[END_REF] suggests to define a constraint ratio as follows:

Constraint ratio = number of dof number of constraints (4.1)

The optimal conditions are fulfilled when constaint ratio in Formula (4.2) is equal to 2 in 2D and 3 in 3D. If we consider a 3D mesh with first-order tetrahedrons, the number of tetrahedrons is approximately equal to five times the number of nodes. As there is one integration point per element and 3 degrees of freedom per node, the constraint ratio is equal to 3/5, which is very bad. Thus, this element cannot be used in practice for such simulations.

Constraint ratio = number of equilibrium equations number of incompressibility equations (

To overcome this difficulty, several finite element formulations have been proposed and are based either on reduced selectively integration schemes or mixed displacementpressure formulations [START_REF] Rice | On numerically accurate finite element solutions in the fully plastic range[END_REF][START_REF] Hughes | The Finite Element Method: Linear Static and Dynamic Finite Element Analysis[END_REF]. In practice, formulations using reduced selectively integration technique are easily applicable only on hexahedral elements (for example, the so-called B-Bar element [START_REF] Hughes | The Finite Element Method: Linear Static and Dynamic Finite Element Analysis[END_REF]). Nevertheless, these elements give rise to heavy meshing operations in order to model complex geometries. The advantage of tetrahedron element is that it can be used for all kinds of geometries and that automatic and robust meshing tools now exist and can be used to mesh very complex geometries.

Mixed formulations [START_REF] Brezzi | On the existence, uniqueness and approximation of saddle-point problems arising from lagrangian multipliers[END_REF][START_REF] Helio | The finite element method with lagrange multipliers on the boundary: circumventing the babuška-brezzi condition[END_REF] can be applied to both hexahedra and tetrahedra by using a Lagrange multiplier and adding other physical unknowns such as pressure (figure 4.2). However, the presence of different kinds of degrees of freedom (displacement and pressure) or even that of additional internal degrees of freedom for certain elements (such as the P1+/P1 element [6,[START_REF] Leblond | Benchmark tests based on the couette viscometer-ii: Thermo-elasto-plastic solid behaviour in small and large strains[END_REF]) can lead to significant extra costs in computation time of non-linear problems. Likewise, high-order elements combined with a reduced integration technique can also be used to avoid the locking problem, but it will lead to a significant increase of the size of the discrete problem and, consequently, of the computation time.

Choosing nodes as integration points is optimal for incompressible or nearly-incompressible materials. Indeed, the constraint ratio is always equal to the optimal values. This idea has been developed in many research works during the last 20 years for 4 node tetrahedral meshes.

The nodal integration technique is one of the many meshless methods, which is one of solution to avoid excessive mesh distortions problem due to large deformations. Numerical integration techniques have therefore been developed for calculating integrals appearing in the weak form of the problem directly from a cloud of points. Among these techniques, the Stabilized Conforming Numerical Integration (SCNI) technique proposed by Chen et al. [START_REF] William Elmer | A stable, meshfree, nodal integration method for nearly incompressible solids[END_REF][START_REF] Chen | A stabilized conforming nodal integration for galerkin mesh-free methods[END_REF][START_REF] Chen | Non-linear version of stabilized conforming nodal integration for galerkin mesh-free methods[END_REF] is extraordinarily efficient. This SCNI scheme satisfies the patch test without using an unfeasible amount of integration points, and there is no volumetric locking problem for incompressible media.

The SCNI method needs to construct a tesselation, in order to apply an integration scheme on its geometrical component. Given a triangulation, each cell is computed by connecting straight lines through the centroid of the triangle and the three midpoints of its edges, as described in figure 4.3. Then different types of shape functions (moving least-squares, local maximum-entropy, linear triangular interpolation) are available [START_REF] Quak | A comparative study on the performance of meshless approximations and their integration[END_REF]. Choosing nodes as integration points in FE methods is not recent, especially in large strain explicit dynamic applications, as coordinate updating with quadrilateral or hexahedral elements may lead to element overlap and tangling. The element choice is a serious matter as the full possibilities offered by adaptive refinement or automatic correction element quality currently cannot be easily exploited with 4-node quadrilateral or 8-node hexahedral elements.

Bonnet et al. [START_REF] Bonet | A simple average nodal pressure tetrahedral element for incompressible and nearly incompressible dynamic explicit applications[END_REF] proposed an average nodal pressure tetrahedral element (ANP) that has been used in explicit dynamics applications involving nearly incompressible materials. The advantage of the tetrahedral element is that the tetrahedron can be easily produced for any geometry by existing powerful generators. However, the actual solutions (linear pressure interpolation and bubble functions) for volumetric locking of tetrahedral elements are unsuitable due to expensive computation time [START_REF] Bell | Suitability of 3-dimensional finite elements for modelling material incompressibility using exact integration[END_REF][START_REF] Thomas | A new finite element formulation for computational fluid dynamics: V. circumventing the babuška-brezzi condition: a stable petrov-galerkin formulation of the stokes problem accommodating equal-order interpolations[END_REF][START_REF] Zienkiewicz | Incompressibility without tears-how to avoid restrictions of mixed formulation[END_REF]. With the average nodal pressure tetrahedron proposed, they avoided volumetric locking but the use of such elements is still computationally costly. Zienkiewicz et al. [START_REF] Zienkiewicz | Triangles and tetrahedra in explicit dynamic codes for solids[END_REF] also avoided the volumetric locking by using a fractional step method. A comparison of these two methods can be found in [START_REF] Bonet | Stability and comparison of different linear tetrahedral formulations for nearly incompressible explicit dynamic applications[END_REF]. Dohrmann et al. [START_REF] Dohrmann | Nodebased uniform strain elements for three-node triangular and four-node tetrahedral meshes[END_REF] introduced a node-based uniform strain element for three-node triangular and four-node tetrahedral meshes. Compared with the element of Bonnet [START_REF] Bonet | A simple average nodal pressure tetrahedral element for incompressible and nearly incompressible dynamic explicit applications[END_REF], the uniform strain elements do not require additional degrees of freedom, what's more, the node-based element is applied not only to the volumetric component of the strain energy (same as Bonnet's element), but also to the deviatoric component, which can improve accuracy and have less sensitivity to shear locking. This idea was firstly applied for the quadrilateral and hexahedron [START_REF] Flanagan | A uniform strain hexahedron and quadrilateral with orthogonal hourglass control[END_REF] and also extensions to other element types [START_REF] Dohrmann | A least-squares approach for uniform strain triangular and tetrahedral finite elements[END_REF][START_REF] Key | A suitable low-order, eight-node tetrahedral finite element for solids[END_REF]. The initial results of node-based uniform strain elements are especially promising. However, the element proposed shows a greatly improved convergence property compared to the standard element for bending dominated problems in small strain elasticity case.

Then, an averaged nodal deformation gradient using linear tetrahedral element has been investigated by Bonnet et al. [START_REF] Bonet | An averaged nodal deformation gradient linear tetrahedral element for large strain explicit dynamic applications[END_REF] for large strain explicit dynamic applications. It is then an extending idea of Dohrmann's works. The gradient calculation is achieved by applying the nodal averaging procedure used for the Jacobian ( [START_REF] Bonet | A simple average nodal pressure tetrahedral element for incompressible and nearly incompressible dynamic explicit applications[END_REF]) to the complete deformation gradient tensor. This nodal integration based formulation can overcome the excessive stiffness of tetrahedrons in bending problems.

Krysl et Zhu et al. [START_REF] Krysl | Locking-free continuum displacement finite elements with nodal integration[END_REF], Puso et al. [START_REF] Puso | Meshfree and finite element nodal integration methods[END_REF][START_REF] Puso | A stabilized nodally integrated tetrahedral[END_REF], Pire et al. [5] have also applied nodal integration technique for finite element method for various tests. The combination of the SCNI technique with the Finite Element method has been the subject of recent work [START_REF] Canales | In-plane/out-of-plane separated representations of updated lagrangian descriptions of viscoplastic flow models in plate domains[END_REF][START_REF] Quak | A comparative study on the performance of meshless approximations and their integration[END_REF]. [START_REF] Chen | A stabilized conforming nodal integration for galerkin mesh-free methods[END_REF] into the finite element method using quadrilateral elements or n-sided polygonal elements to formulate a cell-based smoothed finite element method (CS-FEM or nSFEM) [START_REF] Dai | Free and forced vibration analysis using the smoothed finite element method (sfem)[END_REF][START_REF] Dai | An n-sided polygonal smoothed finite element method (nsfem) for solid mechanics[END_REF][START_REF] Nguyen-Xuan | Smooth finite element methods: Convergence, accuracy and properties[END_REF]. This concept of smoothing domains have been extended to formulate a family of smoothed FEM, such as node-based finite element method (NS-FEM) [START_REF] Liu | A node-based smoothed finite element method (ns-fem) for upper bound solutions to solid mechanics problems[END_REF][START_REF] Nguyen-Thoi | Additional properties of the node-based smoothed finite element method (ns-fem) for solid mechanics problems[END_REF], edge-based finite element method (ES-FEM) [START_REF] Cui | Analysis of plates and shells using an edge-based smoothed finite element method[END_REF][START_REF] Liu | An edge-based smoothed finite element method (es-fem) for static, free and forced vibration analyses of solids[END_REF][START_REF] Nguyen-Xuan | An edge-based smoothed finite element method (es-fem) for analysis of two-dimensional piezoelectric structures[END_REF], face-based finite element method (FS-FEM) [START_REF] Nguyen-Thoi | A face-based smoothed finite element method (fs-fem) for 3d linear and geometrically non-linear solid mechanics problems using 4-node tetrahedral elements[END_REF][START_REF] Nguyen-Thoi | A face-based smoothed finite element method (fs-fem) for visco-elastoplastic analyses of 3d solids using tetrahedral mesh[END_REF]. These smoothed finite elements have been tested for various applications and validated by comparing with classical FEM solutions and analytical solutions.

Liu et al. have combined the strain smoothing technique proposed by Chen

A stable node-based smoothed finite element method has also been implemented(SNS-FEM) [START_REF] Zhang | A fast reanalysis solver for 3d transient thermo-mechanical problems with temperaturedependent materials[END_REF][START_REF] Cui | Steady and transient heat transfer analysis using a stable node-based smoothed finite element method[END_REF] for heat transfer analysis and thermo-elastic analysis. According to the tests, SNS-FEM gives better and smoother solutions than NS-FEM. As we can see in figure 4.4, the stable NS-FEM will create a smoother nodal domain than NS-FEM.

In this chapter, we propose a new 3D finite element method based on arbitrary tetrahedral meshes, nodal integration, and the SCNI technique for solving 3D thermoelastoplastic problems. The aim is to check the method's ability to provide robust and accurate results for the simulation of welding compared to more classical finite element approaches.

We first introduce nodal-integration-based finite element, its advantages and drawbacks. Then different applications and validation tests are also presented.

Presentation of the method

For the sake of simplicity, we place ourselves within the framework of the theory of infinitesimal deformations. As we focus on welding simulations, this assumption is most often sufficient.

Let Ω be the domain of study. The boundary ∂Ω = ∂Ω u ∪ ∂Ω t is decomposed into a part ∂Ω u where the displacement u is prescribed and a part ∂Ω t where surface forces t are imposed.

The finite element method, in its most classical form, allows to calculate an approximation of the displacement field u starting from the virtual work principle which is written, neglecting the dynamic effects:

∀u * Ω u * • bdτ + ∂Ω t u * • tds - Ω ε * : σ dτ = 0 (4.3)
In equation 4.3, u * is a virtual displacement field satisfying u * = 0 on ∂Ω u , b are the volume forces, ε * = 1 2 gradu * + grad T u * , the virtual strain and σ , the stress tensor.

The finite element approximation is a nodal approximation by subdomain. The domain of study is then subdivided in a set of subdomains or finite elements (meshing step), Ω = U elements Ω e with Ω e ∩ Ω f = ∅, connected between them by nodes which support the degrees of freedom of the problem (in general the components of the displacement field). The nodal values of the displacements are then interpolated inside the elements by the shape functions to thus obtain an approximation of the displacement field u on the whole domain Ω.

The integrals appearing in equation (4.3) are then decomposed on each finite element and then calculated by a numerical integration technique (generally using Gauss points); thus, taking the example of the term

Ω ε * : σ dτ : Ω ε * : σ dτ = finite elements Ω e ε * : σ dτ with Ω e ε * : σ dτ = Gauss points ε * g : σ g ω g (4.4)
It is clear from equation 4.4 above that the constitutive equations of the material must be solved at each Gauss point g to calculate the stress tensor σ g . This stress tensor is calculated at each instant from the mechanical state (stresses, internal variables) computed at the previous time step and the strain tensor which is assumed to be known from the displacement field ε * = 1 2 gradv * + grad T v * . Note that this approach obliges to store all the necessary information (internal variables, strains, stresses) at each Gauss point.

The general principle of nodal integration is different in the sense that integrals are not decomposed on each finite element. Thus, for the term Ω ε * : σ dτ, we directly write:

Ω ε * : σ dτ = nodes Ω p ε * : σ dτ = nodes ε * p : σ p ω p (4.5)
In equation (4.5) above, Ω p represents the domain attached to the node p and ω p , its volume. The domains Ω p must constitute a partition of Ω : Ω = nodes Ω p , and be disjoined

(Ω m ∩ Ω n = ∅, m n).
Moreover, in the case of the Gaussian method, the computation of the deformations at each Gauss point results from the interpolation of the displacements directly, the computation of the deformations at a node must be specified.

Definition of nodal domains and nodal strains

The approach chosen is to rely on a mesh of first-order triangles or tetrahedra of any domain of study. The domain associated with a node p is then constituted by the union of the subdomains Ω e p defined as in Figure 4.5 for all the finite elements e containing node p. As suggested by Chen in 2D, the subdomain Ω e p is the domain defined by the node p, the midpoints of the edges connected to node p and the center of gravity of the triangle. In 3D, we propose to define the subdomain Ω e p by the six sub-tetrahedra defined as shown

ε p ij = 1 ω p Ω p ε ij dv = 1 ω p Ω p 1 2 u i,j + u j,i dv = 1 ω p ∂Ω p 1 2 u i n j + u j n i ds (4.6)
Where n i is the component of the outside normal vector, if this last point can be an advantage in the case of meshless methods because it avoids calculating the derivatives of the displacement field, the advantage is less in a finite element context, the computation of the deformations being very usual inside each finite element.

In addition, it can be noted that all the subdomains volumes are equal. The demonstration of this property is immediate for a regular tetraedron and results from the equivalence of the four nodes in this case. In the case of an irregular tetraedron, there exists a unique linear mapping F from a "reference" regular tetraedron onto the actual one.

We therefore have vol

Ω e p = det F vol Ω ref p
. As F is a linear mapping, det F is constant in each element e. The equality of all the subdomains Ω ref p volumes implies that of all the subdomains Ω e p volumes and so vol Ω e p =

1 N e vol (Ω e ) where N e is the number of nodes of the element (3 for a triangle or 4 for a tetraedron). It follows that if the nodal strain components are defined by equation ( 6), then, due to the previous property and the homogeneity of the strain ε e ij within each first-order triangle or tetradron, one has: where S p is the set of elements containing node p. The nodal strains calculated using the SCNI technique on the nodal domains defined above are identical to those coming from the nodal averaging technique used by Dohrmann et al. [START_REF] Dohrmann | Nodebased uniform strain elements for three-node triangular and four-node tetrahedral meshes[END_REF] and in most research works.

ε p ij = 1 ω p Ω p ε ij dv = e∈S p Ω e p ε ij
One could also think to define the nodal strains ε p ij in such a way that the strain field ε nod ij defined by these nodal strains and the shape functions of the elements is equal, in a weak sense, to the strain field calculated using the finite elements. We would then have for each node p :

Ω F p ε nod ij -ε e ij dv = 0 (4.8)
where F p is the shape function associated with node p. As F p vanishes in all the elements not containing node p, it comes:

e∈S p Ω e p F e p         q∈e F e q ε q ij -ε e ij         dv = 0 (4.9)
Concentrating the first member matrix

Ω e
F e p F e q dv on the diagonal and due to the homogeneity of strain ε e ij within each element, it follows:

         e∈S p       Ω e p F e p dv                ε p ij = e∈S p       Ω e p F e p dv       ε e ij = 0 (4.10)
Noticing that

Ω e p F e p dv =
1 N e vol (Ω e ) , we then find that:

ε p ij = e∈S p vol (Ω e ) ε e ij e∈S p vol (Ω e ) (4.11)
The nodal strains thus calculated are again identical to that obtained with the 2 previous techniques.

Definition of nodal thermal strains

The calculation of the thermal strains at a node can simply be obtained from the nodal temperature calculated by the previous thermal analysis. But this option is not consistent with the SCNI approach for which the nodal strains are calculated on average over the nodal volume Ω p (equation (4.6)). Therefore, a second option is to calculate the nodal thermal strain in this way:

ε th p ij = 1 ω p Ω p ε th ij dΩ. ( 4.12)
A last option, that is consistent with the nodal averaging technique, can be to calculate the nodal thermal strains as the average of the thermal strains in the elements containing the node. These three possibilities are compared in the following.

Calculation of stresses and internal variables

Thermal strains have now to be substracted to the nodal strains to get the mechanical part of the nodal strains. Having the nodal mechanical strains, the computation of the stresses and internal variables can now be carried out at each node by solving the constitutive equations of the material.

The method has been implemented in SY SW ELD T M and will be referred as nodal approach in the sequel.

Benefits and drawbacks of the nodal approach

The benefits of the method that has been implemented in SYSWELD TM are the following :

1. The systematic use of tetrahedral meshes.

2. The absence of volumetric locking for incompressible or nearly incompressible media.

3. An easier transfer of the physical quantities: as all the physical quantities are calculated at nodes, interpolations at any point may be performed in a more rigorous way, using the shape functions of the elements. 5. Strong volumetric reduction of the results file: this comes from the fact that the number of nodes in a finite element mesh is generally very much lower than the number of Gauss points. In the case of first-order tetrahedral meshes, the ratio is about 5. The nodal approach nevertheless has some drawbacks. The first one comes from the definition of the nodal strains (formula 4.6). The strains at a node p depend on the displacements of all the nodes belonging to elements containing node p (first neighbor nodes). This choice introduces couplings between each node and its second-neighbor nodes in the tangent matrix (Figure 4.6), while in the classical FEM, couplings only arise between first-neighbor nodes. An increase of the resolution time could, therefore, be expected. In practice, this disadvantage is compensated by the fact that the constitutive equations are now solved at nodes that are less numerous than the Gauss points in the classical FEM.

The second drawback is associated with the possible presence of "hourglass" modes. This zero-energy mode can be remedied by the technique proposed by Puso et al. [START_REF] Puso | Meshfree and finite element nodal integration methods[END_REF][START_REF] Puso | A stabilized nodally integrated tetrahedral[END_REF], which consists in modifying the weak form in the following way:

Ω u * • b dΩ + ∂Ω u * • t dS - Ω ε * : σ dΩ +α                   Ω ε * : σ elast dΩ calculated using nodal integration - Ω ε * : σ elast dΩ calculated using Gaussian integration                   = 0 (4.13)
where α is a small, positive "stabilization parameter", and σ elast denotes the stress tensor calculated assuming a purely elastic behavior. In the expression between parentheses, the first integral is calculated using nodal integration, while Gaussian integration is used for the second.

Applications

Notched tensile specimen

The first classical validation test is to model a notched tensile specimen for compression. The objective is to compare the nodal-integration-based finite element to classical FEM solutions (P1, P1/P1, P1+/P1) on a 3D elastoplastic problem. Here, a 3D numerical model is presented in figure 4.7. As the mesh is irregular, no stabilization parameter is used (α = 0).

In order to take into account the large deformation, an Euler formulation is applied for all elements type. An elastoplastic behavior with isotropic hardening is applied, with module Young E = 200000MP a and Poisson ratio v = 0.3. Yield stress is a function of the cumulated equivalent plastic strain: Firstly, the results obtained from the preliminary analysis of force-axial displacement and CPU time for P1, P1/P1, P1+/P1, and the nodal-integration-based finite element are shown in figure 4.8. The P1 element has poor performance for incompressible media as it is known. The others technique (P1/P1, P1+/P1, nodal-integration-based finite element) provide the same numerical solutions. The figure 4.8 below illustrates that the nodal-integration-based finite element not only gives a correct solution without "locking" phenomena but also is the most efficient among all the chosen elements in terms of computation time.

Secondly, figure 4.9 shows the displacement and axial stress contour results. The nodal-integration-based finite element works as well as the P1+P1 element for stress prediction.

In this study, we have successfully confirmed that the nodal-integration-based finite element can avoid volumetric locking, as mentioned in the literature review. One unan- ticipated finding was that nodal-integration-based finite element is more efficient than actual solutions (P1/P1 and P1+/P1) because the nodal-integration-based finite element can use BFGS (Secant method) solver for non-linear problems while standard Newton method is suggested for P1/P1 and P1+/P1 element for the sake of convergence.

Beam bending

Several reports of Bonnet et al. [START_REF] Bonet | A simple average nodal pressure tetrahedral element for incompressible and nearly incompressible dynamic explicit applications[END_REF][START_REF] Bonet | An averaged nodal deformation gradient linear tetrahedral element for large strain explicit dynamic applications[END_REF][START_REF] Bonet | Stability and comparison of different linear tetrahedral formulations for nearly incompressible explicit dynamic applications[END_REF], and Zienkiewicz et al. [START_REF] Zienkiewicz | Triangles and tetrahedra in explicit dynamic codes for solids[END_REF] have shown that their linear tetrahedral element formulation can overcome the shortcomings in bending. The Objective of this part is to present the performance of the nodal-integrationbased finite element implemented in SY SW ELD T M for static bending simulation. As we know that the linear tetrahedral element suffers the shortcomings in bending dominated problems for large strain. Two meshes have been computed for simulations. Meshes are refined at the edges where there exist high strain and stress gradient. The second-order hexahedral elements can provide the reference result for bending simulation, and nodal-integration-based finite element will be tested. What's more, the P1/P1 element and Q1P0 element are also used as these elements work well for solving the incompressibilities problem. The numerical model (a) in figure 4.11 is a second-order hexahedron elements that contains 1761 nodes and 280 3D elements, and numerical model (b) has 3075 nodes and 2300 3D first-order elements, and numerical model (c) includes 3075 nodes and 13556 tetrahedrons. Figure 4.12 provides displacement contour in Z direction. Table 4.2 presents maximum displacement in Z direction. Compared with the displacement of the Q2 element, Q1P0 element and nodal-integration-based finite element have provided similar displacement distribution, while the P1/P1 element has underestimated the maximum displacement. In this case, the result of the P1/P1 element further supports the conclusion that the standard linear tetrahedral element suffers the shortcomings in bending simulations. The mechanical model of the materials is elastoplastic. The mechanical properties are specified in Table 4.1.

Table 4.1 -Material properties

Module Yuong E = 195122MPa

Poisson ratio v = 0.30

Elastic limit σ Y = 170MPa

Hardening slope constant H = 5000MPa

As shown in Table 4.2, the displacement U z given by Q1P0 element is very approach to that of Q2 element. However, the mesh of the Q1P0 element has about 3 times more integration points. The nodal-integration-based finite element gives 0.04 ′ s difference by using only 3075 integration points.

Numerical simulation of welding

An application to welding simulation, involving very strong thermo-mechanical aspects, will now be presented. Some comparisons of the results with those obtained with more classical FE approaches will be made.

Benchmark model

The aim of this section is to test on a welding simulation the performances of the nodal integration based finite element method with regard to more classical finite element approaches. For nuclear company like FRAMATOME, numerical modeling of welding processes becomes a decision making tool used to speed up and improve the development and qualification of welding and repair techniques. The example we propose comes from the European Network on Neutron Techniques Standardization for Structural Integrity (NeT) Task Group 4. The NeT aims at developing experimental and numerical techniques and standards for the reliable characterization of residual stresses in welded structures [START_REF] Ohms | The european network on neutron techniques standardization for structural integrity net[END_REF]. Task Group 4 deals with the estimation of residual stress fields in a 3-pass slot weld in an 18 mm thick 316L austenitic stainless steel plate. Welding process and type of material are tungsten-inert-gas (TIG) welding and AISI 316L stainless steel [START_REF] Xu | Numerical simulation of a single bead on plate and three pass slot welds in austenitic stainless steel[END_REF]. Dimensions are illustrated in Figure 4.14; the bead length and depth are both somewhat variable, the bead length varying from 74 to 82 mm and the bead depth from 2.1 to 2.7 mm. The welding process parameters are presented in Table-4.3. The simulations are carried out in two steps. A thermal analysis is first performed and provides the temperature distributions during the whole process. Then the calculated temperature distributions are applied as a loading in a thermo-mechanical analysis which gives the resulting stresses and strains. The thermal analysis is carried out using the classical (Gaussian) finite element method on meshes made up of 4-node tetraedra or 8 node hexaedra. Three types of formulation are considered for the mechanical calculation. The first one uses the nodal integration based finite element formulation on tetraedral meshes. The results are compared to those obtained with more classical finite element using mixed displacement-pressure Q1P0 hexaedrons (linear variation of displacements and constant pressure within each element) and P1P1 tetraedrons (linear variations of both displacements and pressure within each element). The comparison of the 3 kinds of simulation is achieved on the residual stresses obtained after the first pass. All the calculations are performed with the finite element code SY SW ELD T M . The welding process parameters are presented in Table 4.3 and all the material properties comes from the SY SW ELD T M database [101].

Meshes

The meshes were created using the VISUAL-MESH TM software. Thanks to the symmetry (see the D-D line in Figure 4.14), only half of the geometry needs to be discretized. The three different meshes considered are presented in Figure 4.15.

The first mesh (A) is intended for element type Q1P0, and the other two (B and C) for the stabilized P1P1 element [START_REF] Feulvarch | A stable p1p1 finite element for finite strain von mises elasto-plasticity[END_REF] and the nodal approach; this is summarized in Table 4.4. Mesh A was generated manually, but both meshes B and C were obtained with an automatic mesh generator. Note that mesh C is more refined than the other two in the weld pass (colored green in Figure 4.15) and in the heat-affected zone.

Thermal simulations

Thermal simulations have been realized for the three different meshes, by using the same parameters given in Table 4. [START_REF] Abdelkhalek | Strip flatness modeling including buckling phenomena during thin strip cold rolling[END_REF] The part is supposed to be initially at 20 C. A two-offsetdouble-ellipsoidal heat source model moving at a constant speed v in the y direction is applied within the elements of the first weld pass, in the same way as Xu et al. [START_REF] Xu | Numerical simulation of a single bead on plate and three pass slot welds in austenitic stainless steel[END_REF]. Equation 4.15 defines the volumetric heat flux (q) inside the front and back regions of the heat source 27 , where these regions are denoted by the subscript 1 and 2 respectively. An element activation/deactivation option is used to simulate the material input in the weld pass.

q 1,2 (x, y, z, t) = 3

√ 3 π √ π f 1,2 Q abc 1,2 exp       -3(y + vt) 2 c 2 1,2       exp -3z 2 b 2 × exp -3(x + x ′ ) 2 a 2 + exp -3(x -x ′ ) 2 a 2 (4.15)
where f 1 = 50.63, f 2 = 25.31, a = 1.799, b = 1.299, c 1 = 1.5, c 2 = 2. These parameters have been adjusted to find by the simulation, the boundary of the experimental melted zone [START_REF] Xu | Numerical simulation of a single bead on plate and three pass slot welds in austenitic stainless steel[END_REF]. The volumic density of power given by equation 4.15 is applied at the Gauss points of the finite elements of the first weld pass. The total power really input in the model may therefore depends on the local mesh where the heat source applies. Therefore, the Q value is calculated at each moment in order to ensure that the total power input is equal to 825 W for the geometry considered.

On all external surfaces of the mesh except the symmetry plane, an exchange coefficient of 10W m -2 ( • C) -1 is applied with an outside temperature of 20 • C. Metallurgical transformations are not included in the simulations since the 316L stainless steel has only one (austenitic) phase.

The temperature distributions obtained with the three meshes during the deposit of the first weld pass are very close, as shown in

Mechanical simulations

The temperature distributions are then put into the mechanical calculation as a loading. The displacement in the x-direction of all the nodes of the symmetry plane are imposed to 0. Three other kinematic conditions are applied to prevent any solid body movement.

A thermo-elastoplastic law with an isotropic hardening is used to describe the AISI316L stainless steel behavior. Young's modulus, yield stress, thermal expansion, and also hardening law depend on temperature. The time of the welding and cooling simulation process is fixed at 1000 seconds in order to obtain a completely cooled state.

The calculation based on the nodal approach has been achieved without any stabilization (α = 0 in equation 4.13) and the nodal thermal strains have been obtained using the nodal temperatures. Table 4.5 gives the CPU times and hardware resources required for all the calculations. All the calculations have been performed on the same computer, with the same precision of convergence on the residual forces and the same number of time steps have been calculated. A BFGS method [11] has been used for solving the non linear problem at each time step for Q1P0 and nodal approach simulations. As the BFGS method did not enable to reach convergence at all the time steps for the P1P1 simulations, a fully Newton-Raphson procedure has been used with this element.

As expected, despite the increase of the bandwidth of the first-member matrix, the CPU time of the nodal approach applied on the meshes B and C is comparable to that of Q1P0 elements and smaller than that of P1P1 elements. The increase of the bandwidth of the nodal approach is thus compensated by the fact that the constitutive equations are solved at nodes and not at Gauss points as for Q1P0 and P1P1 elements. Indeed, because each hexaedron owns 8 Gauss points and each P1P1 tetrahedron one Gauss point, mesh A contains 216000 Gauss points, mesh B, 129277 and mesh C 163798. These values are between 5 and 9 times greater than the number of nodes of mesh B or C.

An other reason which explains that the CPU time of P1P1 elements is greater than that of the nodal approach comes from the fact that with P1P1 elements, each node owns 4 degrees of freedom (3 components of displacement and the pressure) instead of 3 with the nodal approach.

Because all the results of the nodal approach (stresses, strains, internal variables) are stored at nodes instead of Gauss points, the disk space needed by the nodal approach applied on mesh B is 8 to 9 time lower that required by Q1P0 elements on mesh A. But due to the increase of the bandwidth, the RAM needed for the nodal approach is 2 to 3 times bigger than that of the other calculations. Figure 4.17 shows a comparison of the distributions of the longitudinal residual stress obtained with the different approaches after the welding of the first bead. In order to ckeck the results in a more refined way. Figure 4.18 gives the transverse and longitudinal residual stresses according to the depth. The line chosen is located at the middle of the symmetry plane, where high residual stresses are located after the welding process. The residual stresses calculated with Q1P0 mesh are taken as a reference.

The results of the five simulations are very close. Slight oscillations of the residual stresses profiles calculated with the nodal approach on mesh B are observed (Figure 4.18). These oscillations become very small when the mesh is refined in the weld bead and the heat affected zone (mesh C). Calculations have been performed using the stabilization procedure presented in section 2.5 without changing the oscillation level.

Smaller oscillations are also observed on the results obtained with P1/P1 elements while the results given by Q1P0 element are very smooth. But it must be noted that the results shown for the nodal approach come from stresses as calculated in each node while, for the finite element simulations (Q1P0 ans P1/P1), an averaging procedure has been first applied to get the stresses at nodes from the values at Gauss points. Such procedure calculates the stresses at a node by averaging the values obtained at the Gauss point nearest to the node in the elements containing the node.

One way to smooth the results obtained with the nodal approach on a quite coarse mesh is to consider other options of calculation of the thermal strains presented in section of thermal strain defination. So as suggested by Puso et al. [START_REF] Puso | Meshfree and finite element nodal integration methods[END_REF], some averaging procedure can also be applied to the nodal approach results. This procedure consists in first calculating the stresses in each tetrahedron from the stresses calculated at nodes and then, to go back to the nodes by averaging at each node, the stresses so obtained in the elements containing the node.

Stabilization parameter

The stabilization parameter is a solution for "hour-glass" mode, while the portion energy that we have imported could lead to a false solution. Therefore, a study of the effect of the stabilization parameter is still necessary. A zero stabilization (α = 0) has been applied for all the simulations above. The stabilizations tests are carried out with mesh (B).

We have applied α = 0, α = 5 * 10 -4 , α = 5 * 10 -3 as different values. We have extracted the residual stress for the same line in-depth as mentioned in Figure 4.21. For free tetrahedral mesh, the stabilization parameter has a negligible effect or even creates the oscillations. These free tetrahedron meshes generally will not create hour-glass mode; thus stabilization parameter is not recommanded.

Conclusion

A nodal-integration-based finite element method was developed for the numerical solution for incompressibility test, bending problem, and thermo-mechanical problems, and implemented into the SYSWELD TM software. The approach uses linear triangular (in 2D) or tetrahedral (in 3D) meshes, thus permitting the use of standard present-day tools for automatic meshing. As a consequence of the placement of integration points at the nodes instead of the Gauss points of the elements, it avoids volumetric locking phenomena currently encountered with the standard finite element method, for elastoplastic or viscoplastic constitutive equations involving plastic incompressibility.

Nodal subdomains, as required by Chen et al. [START_REF] Chen | A stabilized conforming nodal integration for galerkin mesh-free methods[END_REF][START_REF] Chen | Non-linear version of stabilized conforming nodal integration for galerkin mesh-free methods[END_REF][START_REF] William Elmer | A stable, meshfree, nodal integration method for nearly incompressible solids[END_REF]'s SCNI approach to nodal integration, were defined for both triangular and tetrahedral meshes, as a natural byproduct of the discretization of the domain of study into finite elements. It was shown that with this definition of nodal subdomains, the nodal strains evaluated through the SCNI approach are identical to those calculated using Dohrmann et al. [START_REF] Dohrmann | Nodebased uniform strain elements for three-node triangular and four-node tetrahedral meshes[END_REF]'s nodal averaging technique. A natural weak formulation of the identity of the strain fields in the nodal approach and the standard FEM was also shown to ultimately lead to the same definition of nodal strains. In addition, three different methods were proposed for the calculation of the nodal thermal strains required by thermo-mechanical analyses.

The objective of the notched tensile specimen is a verification of volumetric locking test for nodal integration based finite element. The P1/P1, P1+/P1 are considered as the reference element, as there is no volumetric locking for these two elements. Of course, the P1 element suffers the volumetric locking. The comparison has shown that there is no locking problem for the nodal-based finite elements. What's more, the nodal integration based finite element needs less computation time, which is very encouraging.

The second test is designed for bending simulation. We have used Q2 element (the reference), Q1P0 element, P1/P1 element and nodal integration based finite element. The nodal integration based finit element gives similiar contour like those of Q2 element and Q1P0 element in terme of displacement and stress.

Several calculations of residual stresses induced by some welding process were performed, using both the standard FEM and the nodal approach. The results obtained with the nodal approach are globally consistent with those obtained using the classical FEM with mixed displacement/pressure elements (Q1P0 hexahedra and P1P1 tetrahedra). With the nodal approach, some spatial fluctuations of the residual stresses may be obtained when using a coarse mesh. These oscillations may be reduced and virtually eliminated by refining the mesh, and/or refining the method of calculation of nodal thermal strains, and/or using a "smoothing procedure" of the stress components, similar to that used within the classical FEM when evaluating the average stress components in the elements from their values at Gauss points.

The possible use of linear triangular or tetrahedral meshes, as permitted by the nodal approach without the need for precautions for the prevention of locking phenomena, is a major advantage in the context of numerical simulations of the mechanical consequences induced by welding processes involving complex geometries and/or welding paths. In spite of the increased bandwidth of the left-hand side matrix, the CPU time required by the nodal approach is less than that necessary with the classical FEM with P1P1 elements. Also, convergence of the global "equilibrium iterations" is generally easier, which permits the use of more economical quasi-Newton methods, instead of a full Newton-Raphson method requiring calculation and inversion of a new tangent matrix at each iteration. Finally, with the nodal approach, time-integration of the constitutive equations is performed at the nodes instead of the Gauss points; since the ratio of the total number of nodes over the total number of Gauss points is generally considerably smaller than unity (of the order of 1/5 for a typical large tetrahedral mesh), the disk space necessary to store the results (stresses, strains, internal variables) at each timestep is appreciably reduced, which permits to envisage the simulation of larger welded structures.

In summary, nodal integration based finite element provides promising results and advantages:

1. Without additional new degrees of freedom ( like pressure for P1/P1 element), the nodal approach can avoid the locking phenomena related to plastic incompressibility with standard 4-node tetrahedral elements, and free mesh tool for tetrahedrons is available for complex geometries.

2. The nodal integration based finite element gives promising results in large strain bending simulation, moreover, nodal integration based use less integration points than other elements.

3. According to the welding benchmark test, thanks to the irregularity of the meshes used, adding a stabilization term is not really necessary for the thermal-mechanical simulations.

4. The Nodal integration based finite element works more efficiently than the P1/P1 element as it requires less time computation. However, the nodal integration element needs a reasonable refinement for the high-stress zone.

5. We have also observed that the nodal integration based has a better convergence during the computation. Furthermore, the nodal-integration-based finite element seems less sensitive for mesh tangling than the P1/P1 element.

6. Compared with P1/P1 and Q1P0, the Nodal integration element needs much less disk space for stocking results but more RAM for computation.

Chapter 5

The nodal integration based finite element method with material motion

Simulation of the material motion

In chapter 3, we present a new method to simulate steady-state by simulating material motion in the frame related to the solicitation. With this new method, the mechanical steady-state will be observed after numerous transient analysis steps. However, this new method's mesh must be obtained by translation of a 2D mesh in the welding direction, and the length of each layer must be equal.

In chapter 4, a nodal-integration-based finite element method has been presented and validated by different applications. The nodal-integration-based finite element gives encouraging features. Firstly, the nodal integration based method uses linear tetrahedral mesh, which can be generated easily for all geometry with a meshing tool. Moreover, the nodal-integration-based finite element suffers no volumetric locking problem, better bending performance, more efficient in terms of computing time (compared with P1/P1 element), and less disk space required for storing results. However, the drawback of the nodal-integration-based finite element is that it requires more RAM for storing the stiffness matrix.

In some cases, the FE mesh does not follow the material's motion (case of remeshing, calculation in a moving frame, ...). It is then necessary to « transfer » quantities from the Gauss points of the previous mesh to those of the new mesh. Again this transfer is always based on more or less heuristic and accurate methods. With nodal-integration-based finite element, the transfers are still necessary but may be performed in a more rigorous way, by using the values of the quantities at the nodes and shape functions of the triangle or tetrahedral element. This chapter will present the combination of these two methods for material motion simulation in the moving reference related to solicitation. On the one hand, the linear tetrahedral mesh can be used without locking problem, which can simplify the meshing procedure. On the other hand, the variables transfer for simulating material motion can be performed in a more simple and rigorous manner because all the internal variables (displacement, strain, stress, ...) are computed at nodes. The pre-processor will firstly compute the distance S corresponding to the displacement of material points for a given time step ∆t. Then, a procedure will compute the position of the preceding points for all the nodes. These preceding points can coincide with an existing node, fall into an element, or at the edge of an element.

The preceeding point technique

First of all, we compute each node's displacement S during the time step ∆t in the non-deformed mesh:

S = ∆t * W (5.1)
where W is the material speed in the non-deformed mesh. As we can see in figure 5.1, the point P i is the preceding point of node M i . The internal variables (displacements, strains, stresses, ...) of P i will be computed with the shape function of the element containing P i . If the preceding points coincide with the existing nodes or at edges of the element, the first element found by the pre-processor will be noted. One should note that this transfer procedure will be applied for all the nodes at the beginning of each time step.

The nodes that can not find the element containing his preceding points will keep its internal variables for the next time step, such as the nodes locating at the inlet boundary.

Different motions

Figure 5.2 shows different types of material motions that we will study. We suppose that a material particle's trajectory is the same as his streamline in the non-deformed mesh. Therefore, considering a material particle at a given time t, we find this same particle at the instant t -∆t, by upwinding its streamline.

We consider first the translational motion. For a given time step ∆t, the position ⃗ X P of preceding points (P) for each node (M) can be computed as:

⃗ X P = ⃗ X M -⃗ W ∆t (5.2)
where ⃗ X M is the position of node M, and ⃗ W is the material displacement speed.

Then, we determine the finite element that contains P, and its number and shape functions corresponding to ⃗ X P are stored in memory for a further calculation of internal variables. This idea can also be extended to circular (like pressure vessel and pipe welding simulations) and helical motion (cladding applications, thixoforming simulations). In the case where the material movement is a rotation around an axis passing through the origin with the unit vector ⃗ N = ⃗ Ω ∥ ⃗ Ω∥ ( ⃗ Ω is the rotation speed), the coordinates of the preceding points are computed by applying a rotation of angle ϕ = ∥ ⃗ Ω∥∆t and axis ⃗ N to the vector ⃗ X M .

In general, such a rotation applied to a given vector ⃗ X is obtained in the following way. We start by decomposing ⃗ X in the form:

⃗ X = ( ⃗ X • ⃗ N ) ⃗ N + ( ⃗ X -( ⃗ X • ⃗ N ) ⃗ N ) (5.3) Notice that the vector ( ⃗ X • ⃗ N ) ⃗ N is invariant by the considered rotation. Only the vector ⃗ X ⊥ = ⃗ X -( ⃗ X • ⃗ N ) ⃗ N
, which is orthogonal to the axis of rotation, will be modified. The rotation of angle ϕ and axis ⃗ N applied to ⃗ X ⊥ immediately gives the vector ⃗ Y ⊥ defined by:

⃗ Y ⊥ = cos ϕ ⃗ X ⊥ + sin ϕ ⃗ N ∧ ⃗ X ⊥ (5.4) Finally, noticing that ⃗ N ∧ ⃗ X ⊥ = ⃗ N ∧ ⃗ X the rotation of angle ϕ and axis ⃗ N of any vector ⃗ X leads to ⃗ Y : ⃗ Y = cos ϕ ⃗ X + (1 -cos ϕ)( ⃗ X • ⃗ N ) ⃗ N + sin ϕ( ⃗ N ∧ ⃗ X) (5.5)
These formulas are known as Rodrigues rotation formulas.

By writing the components associated with the above formula, we obtain:

Y i = cos ϕX i + (1 -cos ϕ)N i N j X j + sin ϕε ijk N j X k (5.6)
which gives us the following rotation matrix:

R = cos ϕI + (1 -cos ϕ) ⃗ N ⊗ ⃗ N + sin ϕ            0 -N 3 N 2 N 3 0 -N 1 -N 2 N 1 0            (5.7)
This matrix is thus computed once for all the nodes, and the preceding point's coordinates are determined by:

⃗ X P = R • ⃗ X M (5.8)
In the case of a helical motion, it is necessary to combine translation and rotation. As the translation vector is necessarily collinear with the rotation vector, the combination can be done in any order:

⃗ X P = R • ⃗ X M -⃗ W ∆t (5.9)
The material motion simulation is carried out step by step over time. The only modification to be made in the standard resolution algorithm is that the initial mechanical states of the nodes (displacements, constraints, stresses, ...) must be replaced by the variables of the preceding point at the beginning of each time step. Thus, at the beginning of each time step and for a given instant t, each node's initial mechanical state is exactly the mechanical state of his preceding point at time t -∆t.

The initial mechanical state at time t -∆t of each node is therefore obtained from a simple interpolation of the nodal values of the element containing the preceding point and the shape functions associated with this point. This mechanical state concerns displacements, deformations, stresses, and all internal variables. However, in the presence of material rotation, special attention should be paid to "rotate" the vector and tensor quantities thus calculated and place them in the correct coordinate system. The angle of rotation to take into account is now ∥ ⃗ Ω∥∆t = -ϕ and the corresponding matrix is:

R ′ = cos ϕI + (1 -cos ϕ) ⃗ N ⊗ ⃗ N -sin ϕ            0 -N 3 N 2 N 3 0 -N 1 -N 2 N 1 0           
(5.10) Figure 5.3 provides the algorithm of the material motion simulation procedure. All variables transfer operations take place in the non-deformed mesh. The advantage of working in such a mesh is that the material velocity is constant during the simulation. It is thus not necessary to compute again the position of preceding points (which is time-consuming) at each time step unless ∆t changes.

Mechanical fields transfer

Interpolation techniques associated with the preceeding point technique

First-order interpolation

In order to compute internal variables of preceding points in a simple and efficient way, a first solution is to calculate them using the first-order triangle or tetrahedral shape function and the nodal values. The linear element form functions will be obtained from this approximation. Assuming the values of φ on the nodes of the triangle are known, the coefficients A, B, and C are determined by solving the system:

φ 1 = A + Bx 1 + cy 1 φ 2 = A + Bx 2 + cy 2 φ 3 = A + Bx 3 + cy 3 (5.12)
The Jacobi determinant J is:

J = 1 x 1 y 1 1 x 2 y 2 1 x 3 y 3 (5.13)
where x i and y i are the coordinates of the nodes. Solving 5.12 for A , B , and C leads to:

A = 1 J φ 1 x 1 y 1 φ 2 x 2 y 2 φ 3 x 3 y 3 = (x 2 y 3 -x 3 y 2 ) φ 1 + (x 3 y 1 -x 1 y 3 ) φ 2 + (x 1 y 2 -x 2 y 1 ) φ 3 J B = 1 J 1 φ 1 y 1 1 φ 2 y 2 1 φ 3 y 3 = (y 2 -y 3 ) φ 1 + (y 3 -y 1 ) φ 2 + (y 1 -y 2 ) φ 3 J C = 1 J 1 x 1 φ 1 1 x 2 φ 2 1 x 3 φ 3 = (x 3 -x 2 ) φ 1 + (x 1 -x 3 ) φ 2 + (x 2 -x 1 ) φ 3 J (5.14)
Then the internal variables of P can be computed from equation 5.11, equation 5.14 in the form:

φ p (x, y) = x 2 y 3 -x 3 y 2 + (y 2 -y 3 ) x + (x 3 -x 2 ) y J φ 1 + x 3 y 1 -x 1 y 3 + (y 3 -y 1 ) x + (x 1 -x 3 ) y J φ 2 + x 1 y 2 -x 2 y 1 + (y 1 -y 2 ) x + (x 2 -x 1 ) y J φ 3 (5.15)
We can also express the internal variable P in the reference element as:

φ p (ξ, η) = N 1 φ 1 + N 2 φ 2 + N 3 φ 3 (5.16)
where

N 1 = 1 -ξ -η, N 2 = ξ, N 3 = η.
Another simple and practical way to calculate the contribution of each node for preceding point is to compute the barycentric coordinate (see figure 5.5). The internal values at point P in 2D case can be expressed as:

φ P = w A * φ A + w B * φ B + w C * φ C
w A = S P BC /S ABC w B = S P AC /S ABC w C = S P AB /S ABC (5.17) where φ i represents the physical quantities at node i, and S ABC ,S P BC ,S P AC ,S P AB are the surface of triangles. w A , w B , w C satisfy 0 < w A , w B , w C < 1 and

w A + w B + w C = 1.
Similarly, the barycentric coordinate in 3D can be expressed as: (5.18) where V ABCD ,V P BCD ,V P ACD ,V P ABD , V P ABC are the volume of tetrahedrons. w A , w B , w C , w D satisfy 0 < w A , w B , w C , w D < 1 and w A + w B + w C + w D = 1.

φ P = w A * φ A + w B * φ B + w C * φ C + w D * φ D w A = V P BCD /V ABCD w B = V P ACD /V ABCD w C = V P ABD /V ABCD w D = V P ABC /V ABCD

Second-order interpolation

In order for a better precision and interpolation, we also propose a complete second-order g(⃗ x) interpolation of each quantity g(⃗ x) using the quadratic tetrahedron (or triangle) shape functions:

g(⃗ x) = N N 2 n=1 N n (⃗ x) gn = N N n=1 N n (⃗ x)g n + N N 2 k=N N +1 N k (⃗ x) gk = N s (⃗ x) • G + N m (⃗ x) • G (5.19)
where NN represents the number of element nodes (3 or 4) and NN2 the number of nodes for the corresponding second-order element (6 or 10). g n represents the nodal values of the quantities actually calculated at the nodes. gk are the values estimated at middle edge nodes. The gk values will be computed by reconstructing continuous functions from a set of nodal values1 via a moving least-squares method. We actually have gk = g (⃗ x k ), for

k = N N + 1, ..., N N 2.
Therefore, for all k = N N + 1, ..., N N 2, we should calculate the quantities gk , which can minimize the function: 2 (5.20)

F( G) = N V M i=1 w (⃗ x i ) (g i -g (⃗ x i ))
where NVM is the set of neighboring nodes of the element considered and w(⃗ x i ), a weight function.

Finally, the minimization is carried out by equation 5.21:

N V M i=1 N m (⃗ x i ) T w (⃗ x i ) g i -N s (⃗ x i ) .G -N m (⃗ x i ) . G = 0 (5.21)
equation 5.21 can be expressed as:

       N V M i=1 N m (⃗ x i ) T w (⃗ x i ) N m (⃗ x i )        . G = N V M i=1 N m (⃗ x i ) T w (⃗ x i ) (g i -N s (⃗ x i ) .G) (5.22)
The first member matrix

M =        N V M i=1 N m (⃗ x i ) T w (⃗ x i ) N m (⃗ x i )      
 is a 3x3 or 6x6 matrix, which can be reversed once for all in the case of infinitesimal transformations. In large deformations, it is necessary to take into account the current position of each node and therefore compute the matrix M at each time step.

The evaluation of the quantity at point ⃗ x will, therefore, require computation locally the successive calculations of

Q = N V M i=1 N m (⃗ x i ) T w (⃗ x i ) (g i -N s (⃗ x i ) .G) G = M -1 Q or solve M. G = Q g(⃗ x) = N s (⃗ x).G + N m (⃗ x). G (5.23)
To calculate the shape functions, it is necessary to calculate the barycentric coordinates as in figure 5.6:

The shape functions of the 6-node reference triangle are defined as follows:

N 1 = -α(1 -2α) N 2 = 4ξα N 3 = -ξ(1 -2ξ) N 4 = 4ξη N 5 = -η(1 -2η) N 6 = 4αη (5.24)
and the shape functions of the 10-node tetrahedron are:

N 2 = 4ξα N 1 = -α(1 -2α) N 4 = 4ξη N 3 = -ξ(1 -2ξ) N 6 = 4ηα N 5 = -η(1 -2η) N 7 = 4ζα N 10 = -ζ(1 -2ζ) N 8 = 4ξζ N 9 = 4ηζ (5.25)
Finally, the internal variables of point P can be expressed as:

φ P = N N 2 n=1 N n * φ n (5.26)
where N i is the shape function and φ i are the values at nodes.

The preceeding and subsequent point technique

The current technique consists of seeking his preceding point of each node, which precedes it in the direction of material motion. Thus, for the blue node in figure 5.7, we search for the element (in green) which contains the preceding point and the point itself (indicated by a green star) and the initial mechanical state (at the start of the time step) of the blue point is determined to be that of this point calculated by applying the shape function of the green element.

This technique takes full advantage of the fact that all the mechanical quantities (displacements, deformations, stresses, internal variables) are calculated at the nodes in This technique has the immense advantage of its simplicity. However, it has been reported to create some numerical diffusion related to the variables transfer procedure. This relatively modest numerical diffusion observed with structured meshes in the direction of material motion becomes significantly larger, if the free tetrahedral mesh is applied, especially in the 3D case.

It is clear that this technique has a clear tendency to smooth the transported fields and in particular to reduce the differences between the extreme values.

Each node of the mesh has not only a preceding point but also a subsequent point in the material motion direction. Such as in Figure 5.7, the green node has a subsequent point indicated by a blue star, which is located in the blue element. The initial mechanical state of the point indicated by a blue star is, therefore, that of the green node. We can calculate the initial mechanical states of the subsequent point by using the shape function of the blue triangle, and try to get as close as possible to the mechanical states of the green node.

It is clear that this technique, unlike the previous one, accentuates the differences since we adopt here an approach opposite to the previous one. On the other hand, this technique does not guarantee the existence or the uniqueness of the solution (set of mechanical states at all nodes of the mesh). We will call it "the subsequent point technique".

Therefore, the idea is to combine the preceding point and subsequent point techniques in the least-squares sense to define a technique that will be called "the preceding and subsequent points technique".

To clarify this technique:

v n is the internal variables (strains, stresses, displacements, ...) to be transported to node n at a given time t.

v ′ n is the internal variables (strains, stresses, displacements, ...) transported to the node n after the movement of material between the given time t and t + ∆t.

w p is the internal variables at point p : w p = i N e ip v p = N e p . {v e }, where N e ip is the shape function of material node n at point p.

w ′ q is the internal variables at point q, the material node n will arrive at point q at t + ∆t:

w ′ q = j N f jq v f = N f q • v ′f
The preceding technique gives us internal variables of each node n:

v ′ n = w p = N e p . {v e } (5. 27
)
where p is the preceding point of the node n and e, the element containing this point.

The subsequent technique gives us internal variables of each node n:

v n = w ′ q = N f q • v ′f (5.28)
where q is the subsequent point of the node n and f , the element containing this point.

The problem then consists in finding the values v ′ n which minimize:

F ( v ′ ) = n (1 -α) v ′ n -N e p . {v e } 2 + α v n -N f q . v ′f 2 (5.29)
where α is a weight coefficient for each of the techniques.

We can obtain the solution by solving:

dF = 2 n (1 -α) v ′ n -N e p . {v e } dv ′ n -α N f q • v n -N f q . v ′,f . dv ′f = 0 (5.30)
v ′ is thus solution of the following matrix system:

((1 -α)[I] + α [M s ]) • v ′ = (1 -α) Z p + α{Z} (5.31)
Where [I] is the identity matrix associated with the "preceding point technique", [M s ] is the matrix associated with the "subsequent point technique". m f q = N f q • N f q are the elementary matrices for all the nodes. Z p is the contribution to the second member of " the preceding point technique" and is made up at line number n of the value N e p . {v e }. {Z s } is the contribution to the second member of "the subsequent point technique" and results from the assembly of the vectors z f q = N f q v n .

When the preceding point of node n is outside of the mesh, then we will replaces v ′ n -N e p . {v e } by (v ′ nv n ), at the same time, by adding the term (1α) to the diagonal term of the nth row of the matrix and the term V n on the nth row of Z p . If the node n has a subsequent point, the contribution of this subsequent point is treated in a conventional manner. If the subsequent point of a node n is outside the mesh, this point makes no contribution to the system. One should note that if α = 0, the preceding and subsequent technique is the same as the preceding point technique.

To conclude, various techniques have been implemented. The first-order interpolation and second-order interpolation for computing internal variables of the preceding point have been presented. We have proposed three variable transfer solutions:

1. APPRO 1: the preceding point technique with the first-order interpolation for displacement and the other internal variables.

2. APPRO 2: the preceding point technique with the second-order interpolation for displacement and the other internal variables.

Coupling with thermal analysis

Then, we will try to solve different types of problems. APPRO 1 is always the priority due to his simplicity. If the APPRO 1 can not provide enough precision, APPRO 2 and APPRO 3 will be tested.

Coupling with thermal analysis

For welding steady-state simulation, we proposed the steady-state simulation methods for both thermal and mechanical aspects. However, the temperature simulations for rolling and machining haven't been introduced.

From the physical point of view, it's necessary to take into consideration of the interactions between thermal and mechanical phenomena in the case of rolling and machining simulations. Different from welding simulation, the heat generation in rolling and machining is due to the mechanical dissipation of the workpiece and friction phenomena between the tool and workpiece. Feulvarch et al. [START_REF] Chenegrin | Simulation 3d des transferts de chaleur en perçage par une approche rigid-ale[END_REF][START_REF] Feulvarch | A simple and robust moving mesh technique for the finite element simulation of friction stir welding[END_REF][START_REF] Robe | R-ale simulation of heat transfer during friction stir welding of an aa2xxx/aa7xxx joint on a large process window[END_REF] has developed a method for simulating heat transfer during friction stir welding, and this method allows to take the heat generated by the mechanical dissipation into consideration. This method is used in this study for modeling the temperature during rolling simulation.

Welding steady-state simulation can be found in figure 3.12. Figure 5.9 shows the general procedure for rolling and machining steady-state modeling. According to whether temperature evolution is considered, the procedure can be split into two-part.

Part I can be used to solve the steady-state processes that temperature evolution can be neglected, such as the hot rolling steady-state simulation, wire drawing simulation, and so on.

The temperature plays an important role in certain processes, such as the machining process; therefore, we propose a solution for coupling thermal, mechanical steady-state simulation.

1. We start with Part I simulation, which means the material properties are temperatureindependent. Generally, the material properties at room temperature are chosen if the temperature modeling is in consideration.

2. Once the mechanical steady-state is achieved, the steady-state mesh will be exported.

3. Temperature steady-state simulation will be carried out by using the steady-state mesh obtained by step 1.

4. The steady-state temperature result is applied for the material motion simulation, and the temperature-dependent material properties can be used. 5. Material motion simulation will be continued until a new mechanical steady-state is observed. A new steady-state mesh can be exported.

6. Set up steady-state criteria by the user; for example, temperature evolution simulated in two successive steady-state mesh is inferior 5 degrees or residual stresses are inferior 5 MPa.

7. If the steady-state criteria are satisfied, the simulation will be finished. Otherwise, steps 3-6 will be repeated.

The following sections will introduce some numerical examples to validate the method proposed.

Applications

2D rolling application

Introduction

In this section, we will apply this new method for a 2D cold rolling simulation. The contact condition is defined without friction. At the same time, three transfer techniques will be tested for different mesh types. The aim is to investigate the influence due to the different meshes (structured uniform, structure but non-uniform in operation direction, auto-generate triangle mesh) and transfer solutions. Figure 5.12 presents the mesh for roll (radius r = 25mm) and workpiece (20 mm * 200 mm ). We will first apply this mesh for Lagrange step by step simulation and material motion simulation, which are computed by nodal-integration-based finite element. 3 mm will be laminated, and the moving speed is 1 mm/s. The curves of displacement in the rolling direction and in-depth have also been compared (see figure 5.14). There is no difference in the steady-state zone. Then, in order to reduce computation time, a new mesh is created as we can see in figure 5.15. Compared to the previous mesh (see figure 5.12), this new mesh is only refined nearby solicitation. We consider the first mesh (figure 5.12) as a reference mesh (REF), and the new mesh is a structured non-uniform mesh (STR). All these simulations are carried out with APPRO 1. If we now turn to computation time and CPU time (Table 5.2). Compared with Lagrangian step by step simulation, Material motion simulation with REF mesh requires 180s of computation time to achieve steady-state, while it takes about only 1/4 CPU time of Lagrange formulation. This could be explained as the contact condition is easier to converge in material motion simulation. The simulation with structure mesh needs 190s of computation time, and it only takes 973.02 CPU time. The structured mesh can economize about 15% CPU time than that of REF mesh. What's more, if we want to continue to reduce the CPU time, a shorter numerical model with structured mesh can be also applied for steady-state simulation, which could be more efficient. All the above simulations used APPRO 1 because all the nodes are aligned in the rolling direction. Then an auto-generate tetrahedral mesh (FREE mesh) has also been tested for simulations (see figure 5.17). Similarly, the simulation of Lagrangian step by step simulation and material motion simulations are carried out with the same conditions presented in figure 5.10. Figure 5.19 provides the displacement in rolling direction. APPRO 1 and APPRO 2 are very similar to those of Lagrange step by step simulation. For APPRO 3.2, The displacement decreases in the rolling direction, which is possible due to the subsequent technique. The displacement at curve 1 is not enough representative as all the nodes are aligned in rolling direction. The displacement at curve 2 illustrates that APPRO 2 (second-order interpolation) can reduce numerical diffusion. Table 5.3 provides the computation time and CPU time for steady-state simulations. Among all the simulations, the Lagrange step by step simulation requires the most CPU time. Compared with the simulation with APPRO 1, APPRO 2 requires less CPU time. The explication is perhaps that simulation with APPRO 2 suffers less numerical diffusion. Therefore, it takes less computation time to achieve equilibrium and steady-state. proposed method with APPRO 3.2 190s (0.2s) 1745.5

Conclusion

In this section, we have first tested different mesh (reference mesh, structure mesh, free mesh) and variable transfer techniques (APPRO 1, APPRO 2, APPRO 3.xx).

For the REF mesh and STR mesh, the proposed method gives good agreement compared with Lagrange simulation. In terms of CPU time, the material motion method is more advantages. This advantage comes from the fact that the contact condition is quasi-stable during the simulation. What's more, as all the nodes are aligned in space, the numerical diffusion is unidirectional in the rolling direction and this numerical diffusion in REF and STR meshes is neglectable.

However, when the free mesh is tested, this numerical diffusion becomes two-directional (in-depth and in the rolling direction) and non-negligible. Therefore, APPRO 2 and APPRO 3 are developed for the objective to remedy this numerical diffusion problem. According to figure 5.18 and figure 5.19, APPRO 2 works better than APPRO 1 and APPRO 3 in free triangle mesh.

The choice of variable transfer technique depends on the mesh type. APPRO 1 and is always the priority choice if REF's or STR's type mesh is used. Otherwise, APPRO 2 and APPRO 3.xx can be tested.

Simulation of 3D rolling

Numerical model

In this section, we will extend our method to 3D rolling simulations. A rigid ball will create a groove (0.15 mm) on the sample (plate thickness is 1.39 mm). Figure 5.20 shows the boundary conditions and schema of the process. The same material properties as Table 5.1 are applied. Contact without friction is defined between the ball and the sample. Two numerical models have been prepared (see figure 5.21). The mesh (A) is for the Q1P0 element. The mesh (B) is for the P1/P1 element and the nodal-integration-based finite element. In order to compare the result (displacements, stresses), the mesh (B) has been refined to have more integration points. In the rolling direction, the element length in the hexahedral mesh is 0.207 mm, and that in the tetrahedral mesh is 0.151 mm. Table 5.4 shows the number of integration points for each element type.

Transient simulations comparisons

Among the three simulations, the results of the Q1P0 element are considered as the reference, and we would like to compare the residual stresses simulated by nodal-integrationbased finite element and P1P1 element (solution actual for tetrahedral mesh) to those of the Q1P0 element. In terms of computational efficiency, we would like to make a comparison between the nodal-integration-based finite element and the P1P1 element. A comparison of CPU time and hardware resources has been presented in Table 5.5. The nodal-integration-based finite element requires more RAM due to a larger bandwidth of the stiffness matrix, as each node is coupled not only with his first neighbors but also the second neighbors. The nodal-integration-based finite element needs less disk space for storing the results because there are fewer integrations points than P1P1 elements. In terms of CPU time, the Q1P0 element and nodal-integration-based finite element require almost the same CPU time, while the P1P1 element needs about 3 times more CPU time.

Material motion simulation

The comparisons of transient simulations with different elements have been presented. Then, we will try to simulate this rolling process in a moving reference frame.

The length of model is 45 mm, and 25 mm will be laminated. The rest 20 mm is to ensure the steady-state boundary condition, as we know, if the sphere is too approach to the inlet boundary (inlet boundary for material motion simulation), the boundary effect will become very important. Especially, the plastic deformation should not be observed at the elements of the inlet boundary.

For the material motion simulation procedure, the elements at the inlet boundary represent the initial material state. Therefore, the loading should be placed far away from the inlet boundary. The first attempt is carried out with the same boundary condition as transient simulations, while the nodes of the inlet boundary still suffer a loading that augments progressively. This loading finally leads to deformation plastic non-negligible. A second attempt is to create a zone non-deformable by imposing a new boundary condition, which can ensure zero deformation and zero stress for the inlet boundary (see figure 5.25). For the first steady-state simulation, we start with the mesh (B) as the transient simulation. APPRO 1 is applied for material motion simulations. A free tetrahedral mesh (see figure 5.26) will also be tested for the steady-state simulation. Turning now to stress comparisons, the numerical results with mesh (B) are very similar both in steady-state and the final transient state, which confirms the conclusion that structured mesh suffers less numerical diffusion.

The numerical diffusion can be observed for the numerical result of mesh (C) -APPRO 1, especially for Stress XX. According to figure 5.29, the second-order interpolation APPRO 2 can improve variables transfer quality. Mesh (B) has been used in both transient simulation and material motion simulation, and Table 5.6 shows that material motion simulation requires less CPU time than that of transient simulation. Moreover, APPRO 2 needs more computation for the calculation of second-order interpolation while less CPU time for solving equilibrium equation after the variables transfer, that's why the simulation with APPRO 2 is more efficient.

Conclusion

In this section, we have first presented three transient simulations with different elements: Q1P0 with hexahedral mesh, p1p1 and nodal-integration-based finite element with tetrahedral mesh. The P1/P1 element uses the same mesh as the nodal-integration-based finite element.

The comparison of transient simulation confirms that the Nodal-integration-based finite element gives good simulations results as Q1P0 and it is more efficient than the P1/P1 element. Moreover, the nodal-integration-based finite element requires less disk space for storing transient results but more RAM for storing the stiffness matrix.

The objective is for steady-state simulation. The structured mesh (B) and free mesh (C) are tested. Mesh (B) with APPRO 1 gives similar contours as those computed by transient simulation.

Auto-generate mesh (C) with APPRO 1 shows some numerical diffusion, which can be remedied by APPRO 2.

3D roll forming process simulation

Similarly, we'd like to simulate the 3D roll forming process. Numerical model is shown in figure 5.31. The rotation speed is 0.0333 rad/s, and translation speed is 0.05 mm/s. The same material properties as previous simulations are applied. In this simulation, the structured mesh with APPRO 2 is chosen to reduce numerical diffusion. Figure 5.32 presents the von Mises stress distribution. Simulation is carried out in two steps: firstly, the tool movement is to apply the loading to the workpiece. Secondly, the tool is fixed in space, then, the material motion method is activated to simulate material movement.

In order to study the numerical diffusion due to variable transfer, we decide to show the plastic deformation vs time, because the plastic deformation is a permanent deformation that remains after unloading. According to figure 5.33, the plastic deformation gradually decrease due to numerical diffusion. 3D roll forming simulation is an illustration test of helical material motion. Globally, the results of the simulation are very encouraging, and the numerical diffusion problem can be ameliorated if the mesh is finer. In this part, we will discuss the steady-state simulation of welding process, both MRF method and material motion simulation method are presented. Figure 5.34 shows some boundary conditions, dimension of the numerical model and welding parameters. A double-ellipsoid (Q b1 = 25, Q b2 = 13.5, b1 = 1.3, b2 = 2.6, a = 2.5, c = 3 ) heat source is applied. The total power input is 1007.4 W att. Convection and radiation losses have been taken into account.

Sysweld T M provides thermal and mechanical material properties of 316L stainless steel. The fusion temperature is 1400 degrees. The material properties are temperaturedependent. An elastoplastic with isotropic hardening model is applied in this simulation.

Material motion simulation

With the heat source parameters presented above, the steady-state temperature distribution is shown (see figure 5.35). The temperature distribution is almost the same even if the mesh changes. Therefore, only one temperature figure is presented.

For the MRF method, some specific refinement is required for the sake of convergency. As we can see in figure 5.36, the mesh has been refined at the inlet, outlet boundary, and the domain where we place the heat source. This mesh is obtained by the translation of the 2D section in the welding direction.

For material motion simulation, the first mesh is presented in figure 5.37. All the layers have the same length, and the preceding points of nodes are supposed to coincide The second mesh is a structured mesh. Mesh refinement appears only at the zone where there is the heat source. The objective of structured mesh is to investigate the influence of numerical diffusion existing only in the welding direction. As we can see structured mesh has fewer nodes and 3D elements than those of reference mesh. The same time step and APPRO 1 have been used for two simulations.

A comparison of residual stresses has been presented in figure 5.38. A Lagrange simulation is also performed. Firstly, the material motion simulation with reference mesh (A) gives great agreement with those of MRF method and Lagrange simulation. No obvious differences have been observed. This encouraging report is obtained in the context of without numerical diffusion. Then, the stress contours simulated by mesh (A) and mesh (B) are very similar, which means that numerical diffusion existing only in the welding direction leads to negligible numerical diffusion, according to the stresses contours.

Turning now to comparing the residual stress quantitively, figure 5.39 and figure 5.40 show the curves of residual stresses. The curve yellow is the line where there are high residual stresses. In order to plot the curves of stresses, all the stresses are presented at nodes. An average procedure from gauss points to nodes has been carried out for the results of the MRF method. Figure 5.39 displays the comparisons of the stress distributions of four simulations. Figure 5.39 shows that the stresses are globally very similar. According to the stresses yy in the welding direction, the Nodal approach works very well as the stress curves of step by step simulation almost coincide with those of the MRF method.

From the figures, we can see that the stress of mesh (A) and mesh (B) are almost the same. Thus, the conclusion that the numerical diffusion if only comes from welding direction would lead to negligible influence, can be confirmed again. Therefore, the mesh (B) could be an optimal choice both for the computational efficiency aspect and simulation quality.

Figure 5.40 shows stress distribution in depth. The material motion method based on the nodal integration technique works well. This difference between MRF method and Nodal approach perhaps is perhaps due to the average procedure as there are high residual stress gradients in depth. The average procedure can lead to the modification of real stresses computed at gauss points. We should never forget that the nodal integration based finite element has only 10 integration points in depth while the classical hexahedral element has 18 integration points in this comparison. A better prediction can be produced if there are more integration points in depth. After the reference mesh and structured mesh, a free tetrahedral mesh has also been tested for steady-state simulation (see figure 5.41). As we can see the residual stresses in figure 5.41, the nodal approach and P1/P1 element give almost the same contours with the same mesh (C) in the Lagrangian step by step simulation. While the steadystate simulation with free tetrahedral mesh presents severe numerical diffusion related to first-order interpolation method, this numerical diffusion can be ameliorated by the second-order interpolation.

Circumferential welding

Another application is the circumferential welding, where classical 316L material properties are used. Figure 5.43 presents dimensions of model and boundary condition. Fine mesh is created at the approximation of heat source. The number of nodes and 3D elements is shown.

Thermal steady-state should be solved by diffusion convection equation firstly. Then steady-state temperature distribution can be founded in figure 5.44.

Firstly, we try to use MRF method for solving mechanical steady-state. However, the MRF method can not converge. After that, The material motion method is applied. The stresses at different moments are shown in the figures below.

The material motion simulation method is a step by step simulation. Therefore, the method proposed is more time consuming, but it is easier to converge than MRF method. Compared with Lagrangian simulation, the advantage is that it is unnecessary to use fine mesh along all the heat source trajectory.

Conclusion

In this section, the material motion method is applied for steady-state welding simulation. The simulations of Lagrange formulation based on nodal approach and MRF method can be considered as the reference. Different meshes are tested.

The simulation results of nodal-mesh (A) can avoid the numerical diffusion related to the variable transfer procedure. The comparisons of stresses simulated by MRF method and nodal-mesh (A) show no difference. Therefore, the material motion method is validated for steady-state welding simulations.

The results obtained by nodal-mesh (B) show that the numerical diffusion only existing in the welding direction has negligible influence for steady-state welding simulation. Table 5.7 presents that the mesh (B) can save considerable computing resources (RAM, CPU time, disk space).

The first-order interpolation and second-order interpolation method are applied for the mesh (C). The second-order interpolation method has improved the results' quality but the numerical diffusion cannot be eliminated. Therefore, how to remedy the numerical diffusion due to variable transfer for the free tetrahedral mesh could be the future study. As we can see in Table 4.5, the MRF method is the most efficient method in plate welding simulation. However, the MRF method cannot converge for circumferential welding simulation. The material motion method with selective refinement mesh (B) can be used for steady-state simulation, which can not only greatly shorten computing time but also reduce disk space required for storing results.

3D thermal-mechanical rolling simulation

This numerical example presents the possibility of taking the changing of temperature into consideration for rolling simulation. We have introduced the procedures in figure 5.9; the simulation starts with ambient temperature T = 25 degree (Part I). 316L thermalmechanical properties are used for thermal-mechanical simulation. The basic information is given in figure 5.47. Friction between tool and workpiece is neglected. The symmetrical conditions are used (top-base, left-right. The mesh used is only for a quatre of the complete model. The properties for temperature modeling are shown in Table 5.8 and Table 5.9. A node-to-node formulation for non-matching meshes developed by Feulvarch et al. [START_REF] Geuffrard | Symmetrical node-tonode formulation for thermal contact problems between non-conforming or nonmatching meshes[END_REF] is used.

Once the temperature is obtained, we will continue the material motion simulation with this new temperature load. The final mechanical state (t = 100s) will be considered as the initial state for the following simulation. Figure 5.49 provides the temperature and von Mises distribution. Then, a second steady-state mechanical mesh at t = 200s is exported for the second temperature calculation. The temperature at t = 100s is considered as the initial condition for temperature calculation, and the stress at t = 200s is the initial state for mechanical simulation. Figure 5.50 shows the temperature and von Mises stress. Compared with previous temperature and stress distribution, one can note that the workpiece's temperature has changed while the von Mises stress is almost the same.

The temperature plays an important role in stress and plastic cumulation distribution, as we can see that the von Mises stress decreased about 150 MP a.

To conclude, the Part I presented in figure 5.9 provides solutions for the simulations that neglected the temperature variation due to deformation plastic, such as hot rolling simulation.

Part II is to take temperature into consideration. The thermal-mechanical coupling is carried out in an uncoupled manner, as we can see in figure 5.9. Part II could be continued until the difference of stresses or temperature between two successive cycles satisfies your criteria.

Chapter 6

Conclusion and Perspectives

Conclusion

In this thesis, the research works have two principal parts: the first part concerns the validation of nodal-integration-based finite element method. The second part is the development of material motion simulation method.

Nodal-integration-based finite element method not only allows to use of tetrahedral elements without locking problems related to von Mises criteria but also permits to compute all the internal variables (displacement, strain, stress, ...) at nodes. As we know, the triangle and tetrahedron can be generated automatically by existing meshing tools, especially for complex industrial pieces. Several tests have been studied to validate the nodal-integration-based finite element method.

1. Notched tensile specimen test has shown that the nodal-integration-based finite element method can avoid locking problems (contrary to standard linear tetrahedral elements in von Mise plasticity). Another finding is that nodal-integration-based finite element method is more efficient than the actual solution (P1P1 element, P1+P1 element).

2. The bending test has shown that the P1P1 element is too stiff. However, nodalintegration-based finite element method gives a good performance in bending dominated problems and more efficient than the P1P1 element.

3. TG4 benchmark simulation shows that the nodal-integration-based finite element method works very well for thermo-mechanical problems. The nodal-integrationbased finite element method results are globally consistent with those obtained using Q1P0 and P1P1 elements.

In these tests, the results of the nodal-integration-based finite element method are always compared with a reference solution, such as P1+P1 for locking test, Q2 for bending test, Q1P0 for welding simulation. The nodal-integration-based finite element gives excellent agreement with the reference solution. Moreover, the nodal-integration-based finite element is more efficient and requires less disk space for storing the results than the P1P1 element for tetrahedral meshes, while more RAM is needed for storing the stiffness matrix.

For the following parts, the nodal-integration-based FEM is used for solving the problem in a moving reference frame related to the solicitations. Firstly, the material motion is simulated by applying some preceding points technique at each time step. One should note that we search the position of the preceding points for all the nodes in the non-deformed mesh, in which the velocity of all the nodes is constant. Secondly, at the beginning of each time step, all nodes' internal variables should be replaced by those of preceding points, which could be computed by the first-order shape function and nodal values of the element containing the preceding point. The equilibrium state is calculated at each time step.

With the preceding point technique, the translational, circular, and helical material movement can be solved. A relatively modest numerical diffusion can be observed with structured tetrahedral mesh. However, numerical diffusion becomes significantly important if a free tetrahedral mesh is used. In order to improve the numerical diffusion problem in free meshes, the second-order shape function and "the preceding and subsequent points technique" have been implemented.

Finally, the material motion method is applied to various tests. According to these tests, we can conclude:

1. If a structured mesh is used for material motion simulation, the APPRO 1 is the prior choice. If the important numerical diffusion is observed, APPRO 2 and APPRO 3.XX can be tested. 3. The material motion method can be used for steady-state simulation, while this new method can also provide the information of transient states.

The nodal-integration-based finite element method with material motion has been validated by comparing simulation results with those obtained by the Lagrangian formulation and MRF method. For transient-state simulation, this new method is an incremental method similar to the Lagrangian method, while the new method requires a mesh refinement only near the solicitation. This refinement selective can largely save the computing resource compared to the Lagrangian formulation. For steady-state simulation, this new method is much easier to converge but more time-consuming than MRF method as it needs to simulate transient states, and it permits also to simulate the steady-state of large deformation problem.

Perspectives

Improvements of the proposed methods

The nodal-integration-based finite element method and material motion simulation have been applied for various typical tests, and the results given by the proposed method have a good consistent with those obtained by classical methods. However, there are still some improvements that should be made: As we can see in figure 6.1, there are the symmetrical plan and free surface and the σ 33 is perpendicular to the free surface. The first improvement concerns the stress calculation on the free surface. It should be noted here that with the nodal-integration-based FEM, the computation of the stresses in a node located on the surface takes into account the volume located under this node and thus does not respect precisely the condition on the surface. The finite element method does not respect either the condition on the surface (in fact it respects it but in a weak sense) but it is less annoying because the computation is carried out at the Gauss points which are generally not located on the surface. Therefore, the resolution of the behavioral equations must be by adding the supplementary conditions to be respected on the stresses on the free surface. Bergheau et Leblond has already worked on this point and a solution is under development and validation.

The second point to be improved is to take the phase transformation into account for thermal-mechanical simulation. As it is well-known that phase transformations play an important role in residual stresses and distortions, especially for welding and heat treatment.

The third improvement is vised to remedy numerical diffusion observed in the case of free mesh. The diffuse approximation method proposed by Rassineux et al. [START_REF] Yang | Simulation of continuous damage and fracture in metal-forming processes with 3D mesh adaptive methodology[END_REF][START_REF] Yang | A 3d hadaptive local remeshing technique for simulating the initiation and propagation of cracks in ductile materials[END_REF] works very well for variables transfer. Therefore, applying the diffuse approximation method for material motion simulation could be possible to remedy numerical diffusion.

Machining simulation

The machining process is governed by several phenomena, such as the energy dissipation due to plastic deformation, chips formation, the fracture of chips, energy transferring between chips and tool and the workpiece, and so on.

In the literature, lots of numerical models and finite element or meshless methods have been reported [START_REF] Prieto | Numerical methods for the modelling of chip formation[END_REF], and we will not discuss the advantages and drawbacks of each method in this section. The objective is to present the possibility to simulate the steady-state of machining processes by the material motion method presented in chapter 5. Steady-state simulation is still necessary for specific requirements. For example, in order to predict the remaining cutting tool life for maximizing the utilization of each machining tool, the cutting force that is applied to the cutting tool in steady-state is the main factor for predicting tool wear.

In this section, a simple illustration of steady-state mechanical simulation is presented. Therefore, the phenomena related to material separation criteria are beyond this research consideration. In addition, as the aim of the simulation is to calculate the steady-state of the process, there is no need for a material separation criterion. 1. Advancing the tool until to the end of the pre-crack.

2. Once the tool arrives at the end of the pre-crack, the material motion is simulated until observation of steady-state.

An elastoplastic model with isotropic hardening is chosen. Young's modulus is 144000 MP a, the yield stress is 229 MP a, and the strain hardening is presented in the table 6.2. The friction between the tool and workpiece can also be taken into account (we neglected the friction in our test). Figure 6.3 provides a schema explication for steady-state simulation and also the displacement UY during the simulation. The material motion simulation can be considered a manner of correction to revise the false mechanical states due to unreal material separation. Then, a 3D machining simulation is also presented. Figure 6.4 shows the model dimensions and boundary conditions. SY SW ELD T M has some macro-elements that are designed for the contact between a deformable body (with large relative displacements) and one or several rigid targets possibly mobile. Contact conditions and material properties are the same as 2D simulation. A pre-crack has been created for the zone of separation of material. The presence of cutting tools is for considering heat exchange between the cutting tool and workpiece. Figure 6.5 shows a structured tetrahedral mesh to reduce numerical diffusion related to the variables transfer procedure. Refinement can be found for the chip due to the existence of a high gradient of strain/stress. Figure 6.6 presents von Mises' stress distribution during the simulation. The first step is tool movement until the separation zone. The second step is the material motion simulation. The chip formation can be observed. The steady-state is achieved at 80s. In this section, the nodal integration based material motion simulation method has been successfully applied for machining simulation. The object is to give a simple demonstration and test the possibility of modeling machining processes. The model proposed permits to couple a thermal analysis and a mechanical calculation. Temperature changes during the machining process should be modeled in the future study.

Besides, it is interesting that the nodal-integration-based finite element method could be implemented in conjunction with a remeshing technique or mesh adaptation technique to avoid mesh tangling due to large deformation.

Simulation numérique des états stationnaires associés aux processus thermomécaniques

Résumé

De nombreux procédés de fabrication thermomécanique comme le laminage, le soudage ou encore l'usinage mettent en jeu soit des sollicitations mobiles par rapport à la matière fixe, soit de la matière mobile par rapport à des sollicitations fixes. Dans tous les cas, après un régime transitoire en général assez court, les champs thermiques, métallurgiques et mécaniques associés à ces procédés atteignent un état stationnaire. La recherche de ces états stationnaires à l'aide de la méthode des éléments finis classique nécessite de mettre en oeuvre des modèles complexes et couteux où les sollicitations se déplacent par rapport à la matière (ou l'inverse).

La recherche directe des états stationnaires a fait l'objet de nombreux travaux de recherche ces trente dernières années. Des méthodes sont aujourd'hui disponibles et pour certaines sont proposées dans des codes de calcul du commerce. Ainsi, une option de calcul dite repère mobile proposée par différents auteurs est disponible dans le logiciel SY SW ELD T M . Cette méthode permet de calculer les états thermique, métallurgique et mécanique stationnaires associés à un procédé de soudage, en résolvant un problème de diffusion-convection en thermique et en intégrant, en mécanique, les équations constitutives du comportement du matériau le long des lignes de courant. Si cette méthode a été utilisée avec succès dans de nombreuses applications, elle présente néanmoins quelques limitations. Ainsi le maillage doit être structuré et la convergence des calculs est en général assez lente.

Nous proposons dans cette thèse de résoudre le problème mécanique dans un repère liè aux sollicitations, en nous appuyant sur une méthode de calcul par éléments finis reposant sur l'intégration nodale et la technique SCNI (Stabilized Conforming Numerical Integration). Cette méthode permet l'utilisation de maillages en tétraèdres (ou triangles en 2D) sans rencontrer de problème de verrouillage volumique résultant de l'incompressibilité plastique associée au critère de plasticité de von Mises. Plutôt que de rechercher directement l'état stationnaire, l'idée générale est ici de construire l'état stationnaire à partir d'une analyse transitoire en faisant entrer pas à pas la matière par l'amont et en la faisant sortir par l'aval d'un maillage fixe par rapport aux sollicitations et de taille limitée. L'état (quasi-)stationnaire n'est donc atteint qu'au bout d'un certain temps d'analyse. Les avantages de cette méthode résident dans:

1. L'utilisation d'un maillage libre en tétraèdres ou en triangles (au lieu des maillages structurés), 2. Un maillage raffiné uniquement dans la zone située au voisinage des sollicitations, 3. Une grande robustesse et notamment une bonne convergence des calculs mécaniques non linéaires liée à la résolution des états transitoires.

Après une introduction générale (Chapitre 1) et un état de l'art sur les méthodes existantes (Chapitre 2), nous présentons une approche de simulation du mouvement de matière dans le cadre de la méthode des éléments finis classique sur un problème de soudage (Chapitre 3). Nous y proposons également des conditions aux limites thermiques pertinentes pour calculer directement la distribution de températures en régime stationnaire.

La méthode des éléments finis reposant sur l'intégration nodale est ensuite décrite au Chapitre 4. Les avantages et inconvénients de la méthode sont discutés. La méthode est validée sur une application en grandes déformations élastoplastiques, un problème de flexion et une simulation thermomécanique de soudage.

La méthode des éléments finis reposant sur l'intégration nodale est alors développée pour prendre en compte un mouvement de matière (Chapitre 5). Trois types de mouvement sont considérés : en translation, circulaire et en hélice. Différentes méthodes de transport de champ sont abordées et discutées ainsi que le couplage thermomécanique. Des exemples d'application dans le domaine du laminage et du soudage pour différents mouvements de matière montrent l'efficacité de la méthode développée

3. 1

 1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Thermal steady-state computation . . . . . . . . . . . . . . . . 3.2.1 How to choose boundary conditions . . . . . . . . . . . . . . . . . .

3. 3

 3 Mechanical steady-state computation . . . . . . . . . . . . . .

3. 3 . 1

 31 Simulation of material motion . . . . . . . . . . . . . . . . . . . . .

3. 4

 4 Application to 3D welding simulation . . . . . . . . . . . . . . 3.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Chapter 4 The nodal integration based finite element method Introduction Contents 1.1 INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . 1.1.1 Rolling processes . . . . . . . . . . . . . . . . . . . . . . 1.1.2 Welding processes . . . . . . . . . . . . . . . . . . . . . . 1.1.3 Machining processes . . . . . . . . . . . . . . . . . . . . 1.1.4 Numerical simulation of thermo-mechanical processes 1.1 INTRODUCTION

Figure 1 . 1 -

 11 Figure 1.1 -Example of the rolling process.

Figure 1 . 2 -

 12 Figure 1.2 -Diagram of arc and weld area, in shielded metal arc welding. 1. Welding torch 2. Wire 3. Weld pool 4. Heat-affected zone 5. Base metal.

Figure 1 . 3 -

 13 Figure 1.3 -Fracture problems in Liberty ship related to some defects in the welding technique.

Figure 1 . 4 -

 14 Figure 1.4 -Schematic diagram for chip formation.

Figure 1 . 5 -

 15 Figure 1.5 -Physical phenomena involved in welding simulation -couplings and interactions.

Figure 2 . 1 -

 21 Figure 2.1 -Procedure description for adaptive re-meshing [17].

Figure 2 . 2 -

 22 Figure 2.2 -The eulerian reference frame for rolling and welding simulation

Figure 2 . 3 -

 23 Figure 2.3 -Mixed four-node quadrilateral element [88].

Figure 2 . 4 -

 24 Figure 2.4 -Definition of moving reference [8].

Figure 2 . 5 -

 25 Figure 2.5 -Gauss points sequences.

Figure 2 . 6 -

 26 Figure 2.6 -Schematic representation of Eulerian, Lagrangian, and ArbitraryLagrangian-Eulerian formulation[START_REF] Boman | Finite element simulation of lubricated contact in rolling using the arbitrary lagrangian-eulerian formulation[END_REF].

Figure 2 . 7 -

 27 Figure 2.7 -Configurations involved in the mathematical modeling of a rolled strip [105].

Contents 3 . 1

 31 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Thermal steady-state computation . . . . . . . . . . 3.2.1 How to choose boundary conditions . . . . . . . . . . . .3.3 Mechanical steady-state computation . . . . . . . . 3.3.1 Simulation of material motion . . . . . . . . . . . . . . .

3. 4

 4 Application to 3D welding simulation . . . . . . . . 3.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 3 . 1 -

 31 Figure 3.1 -A moving reference frame related to heat source.

Figure 3 . 2 -

 32 Figure 3.2 -Diagram of material motion simulation in the reference frame.

Figure 3 . 3 -

 33 Figure 3.3 -Numerical model for thermal transient simulation.

Figure 3 . 4 -

 34 Figure 3.4 -A double-ellipsoid heat source.

Figure 3 . 5 -

 35 Figure 3.5 -The steady-state temperature distribution in the step by step simulation.

Figure 3 . 6 -

 36 Figure 3.6 -Upstream model.

Figure 3 . 7 -

 37 Figure 3.7 -Downstream model.

Figure 3 . 8 -

 38 Figure 3.8 -Material flow passes through the computation frame.

Figure 3 . 9 -

 39 Figure 3.9 -The numerical model for steady-state simulation.

Figure 3 .

 3 Figure 3.10 shows the temperature distribution with two different boundary conditions. Compared with the reference solution given by transient step by step (figure 3.5), the simulation with the new boundary condition proposed gives a better prediction. Then, a comparison of the temperature profile is presented in figure 3.11. The plotline chosen locates at the symmetry plan (see figure 3.8). The temperature profile of the new boundary condition gives perfect agreement with the reference.

Figure 3 . 10 -

 310 Figure 3.10 -Steady-state temperature distribution simulated with classical boundary condition (A) and new boundary condition (B).

Figure 3 . 11 -

 311 Figure 3.11 -Comparison of temperature in the welding direction (0.83mm from top surface).

Figure 3 . 12 -

 312 Figure 3.12 -Procedure for welding steady-state simulation.

3. 4 .

 4 Application to 3D welding simulation

Figure 3 . 13 -

 313 Figure 3.13 -The material motion simulation procedure.

Figure 3 . 14 -

 314 Figure 3.14 -Mechanical boudary conditions.

Figure 3 .

 3 Figure 3.15 -Steady-state longitudinal stress distribution, (a) MRF (b) method proposed.

Figure 4 . 1 -

 41 Figure 4.1 -Selective reduced integration scheme in 2D.

Figure 4 . 2 -

 42 Figure 4.2 -P1/P1 and P1+/P1 scheme in 2D.

Figure 4 . 3 -

 43 Figure 4.3 -Geometrical objects used for computation [89].

Figure 4 . 4 -

 44 Figure 4.4 -The schematic of NS-FEM and SNS-FEM in 2D and 3D[34].

4 .

 4 Local remeshing in case of excessive distortion of the elements around a node: one can change the local connectivity of the elements without moving the nodes, thus avoiding any transfer of the mechanical fields.

Figure 4 . 6 -

 46 Figure 4.6 -Couplings between node and its neighbors.

Figure 4 . 7 -

 47 Figure 4.7 -Dimensions of the specimen (unit: mm) and mesh from two different angles (7530 nodes, 33185 elements).

Figure 4 . 8 -

 48 Figure 4.8 -Force-axial displacement and CPU time.

Figure 4 . 9 -

 49 Figure 4.9 -Displacement and stress distribution of P1+/P1 (reference) and nodalintegration-based finite element.

Figure 4 .

 4 10 present the bending model (10*0.5*4 mm) and load and boundary conditions (ex: SYM U y means a condition symmetrical in Y direction).

Figure 4 . 10 -

 410 Figure 4.10 -Load and boundary conditions for bending test.

Figure 4 . 11 -

 411 Figure 4.11 -Numerical model for simulations.

Figure 4 .

 4 Figure 4.13 shows the von Mises stresses. Q2 element and linear tetrahedral proposed and Q1P0 element provide similar stresses' contours, while the P1/P1 element has a poor performance in predicting the von Mises stress.

Figure 4 . 12 -

 412 Figure 4.12 -Displacement contours.

Figure 4 . 13 -

 413 Figure 4.13 -Von Mise stress contours.

Figure 4 . 14 -

 414 Figure 4.14 -Diagram of the 3-pass slot weld sample [107].

Figure 4 . 15 -

 415 Figure 4.15 -Meshes -Top: global meshes -Middle: 2D sections -Bottom: zooms on the weld zone.

  Figure 4.16.

Figure 4 . 16 -

 416 Figure 4.16 -Temperature distributions calculated at t = 34 s with (a) mesh A, (b) mesh B, (c) mesh C.

Figure 4 . 17 -

 417 Figure 4.17 -Longitudinal σ yy residual stress for (a) Q1P0 elements, (b) P1P1 elements with mesh B, (c) P1P1 elements with mesh C, (d) Nodal approach with mesh B, (e) Nodal approach with mesh C (unit: MPa) at t = 1000s.

Figure 4 .

 4 Figure 4.19 gives the residual mean stress (or pressure) distributions obtained on mesh B with the different options proposed. Clearly the results obtained when averaging the temperature on all the elements containing each node (case (b) of figure 4.19) or on the nodal volume (case (c)) are better. The nodal volume averaging option which is consistent with the SCNI approach is therefore recommended.

Figure 4 .

 4 Figure 4.20 give the mean stress (pressure) obtained using this procedure. The

Figure 4 . 18 -

 418 Figure 4.18 -Comparison of residual stresses σ xx and σ yy in-depth.

Figure 4 . 19 -

 419 Figure 4.19 -Mean stress obtained on mesh B with various options for the calculation of nodal thermal strains -(a) Nodal temperature, (b) nodal subvolume averaging, (c) element averaging.

Figure 4 . 20 -

 420 Figure 4.20 -Mean stress distributions -(A) Mesh A with elements Q1P0, (B) Nodal approach with smoothed results.

Figure 4 . 21 -

 421 Figure 4.21 -Comparison of residual stresses σ xx and σ yy in-depth with different stabilization parameter.

Contents 5 . 1

 51 Simulation of the material motion . . . . . . . . . . 5.1.1 The preceeding point technique . . . . . . . . . . . . . . 5.1.2 Different motions . . . . . . . . . . . . . . . . . . . . . . 5.2 Mechanical fields transfer . . . . . . . . . . . . . . . . 5.2.1 Interpolation techniques associated with the preceeding point technique . . . . . . . . . . . . . . . . . . . . . . . 5.2.2 The preceeding and subsequent point technique . . . . 5.3 Coupling with thermal analysis . . . . . . . . . . . . 5.4 Applications . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.4.1 2D rolling application . . . . . . . . . . . . . . . . . . . . 5.4.2 Simulation of 3D rolling . . . . . . . . . . . . . . . . . . . 5.4.3 3D roll forming process simulation . . . . . . . . . . . . 5.4.4 Simulation of welding up to steady-state . . . . . . . . . 5.4.5 3D thermal-mechanical rolling simulation . . . . . . . .

Figure 5 .

 5 Figure 5.1 shows a diagram of the general principle for simulating material motion. The pre-processor will firstly compute the distance S corresponding to the displacement of material points for a given time step ∆t. Then, a procedure will compute the position of the preceding points for all the nodes. These preceding points can coincide with an existing node, fall into an element, or at the edge of an element.

Figure 5 . 1 -

 51 Figure 5.1 -Schema principal for material motion simulation.

Figure 5 . 2 -

 52 Figure 5.2 -Principle of translational, circular and helical motion.

Figure 5 . 3 -

 53 Figure 5.3 -Algorithm for material motion simulation.

Figure 5 . 4 -

 54 Figure 5.4 -Triangular and reference element.

Figure 5 . 5 -

 55 Figure 5.5 -Schema of barycentric coordinate for triangle and tetrahedron.

Figure 5 . 6 -

 56 Figure 5.6 -Barycentric coordinate for second-order element.

Figure 5 . 7 -

 57 Figure 5.7 -Principle of the preceding and the subsequent point.

Figure 5 . 8 -

 58 Figure 5.8 -Principle of the preceding and subsequent technique.

Figure 5 . 9 -

 59 Figure 5.9 -The general procedure for steady-state simulation.

Figure 5 . 10 -

 510 Figure 5.10 -Lagrange formulisme and proposed method, boudary conditions.

Figure 5 .

 5 Figure 5.10 shows two different methods and boundary conditions applied. The Lagrange formulation is considered as the reference. Material properties are presented in Table 5.1. An isotropic hardening model is used.

Figure 5 . 11 -

 511 Figure 5.11 -Stress hardening function of plastic deformation for material properties in Table 5.1.

Figure 5 . 12 -

 512 Figure 5.12 -Mesh for roll and workpiece.

Figure 5 .

 5 Figure 5.13 shows the contours of displacement U x. The first contour is computed by Lagrange step by step. The other contours are simulated by the material motion method. These two methods give similar contours.

Figure 5 . 13 -

 513 Figure 5.13 -The U x displacement simulated by Lagrange formulation and material motion method.

Figure 5 . 14 -

 514 Figure 5.14 -Comparisons of displacements in rolling direction and in-depth at final state.

Figure 5 . 15 -

 515 Figure 5.15 -The non-uniform structured mesh.

Figure 5 .

 5 Figure 5.16 presents the comparisons of displacement obtained by reference mesh and structured non-uniform mesh. The same displacements are observed. The only difference is the computation time required to achieve a steady-state.

Figure 5 . 16 -

 516 Figure 5.16 -Comparisons of displacement in rolling direction.

Figure 5 . 17 -

 517 Figure 5.17 -The auto-generate tetrahedral mesh.

Figure 5 . 18 -

 518 Figure 5.18 -Displacement contours of Lagrange simulation and material motion method with different variable transfer techniques.

Figure 5 . 19 -

 519 Figure 5.19 -Comparisons of displacements in rolling direction at final state.

  The moving speed of ball is 0.25 mm/s. The dimension of model is 45 * 2.34 * 1.39 mm.

Figure 5 . 20 -

 520 Figure 5.20 -Schema of the process and boundary conditions.

Figure 5 . 21 -

 521 Figure 5.21 -Mesh (A) for Q1P0 element, mesh (B) for P1/P1 element, and nodal-integration-based finite element

Figure 5 .

 5 Figure 5.22, figure 5.23, figure 5.24 shows the distribution of stresses. Compared with the results of Q1P0 element, both P1P1 and nodal-integration-based finite element give good agreement. The simulations use the update Lagrange formulation to take large strain into account.

Figure 5 . 22 -

 522 Figure 5.22 -von Mises (unit: MPa) for Q1P0 elements, P1/P1 elements and Nodal-integration-based finite element.

Figure 5 . 23 -

 523 Figure 5.23 -Stress XX (unit: MPa) for Q1P0 elements, P1/P1 elements and Nodal-integration-based finite element.

Figure 5 . 24 -

 524 Figure 5.24 -Stress YY (unit: MPa) for Q1P0 elements, P1/P1 elements, and Nodal-integration-based finite element.

Figure 5 . 25 -

 525 Figure 5.25 -Boundary conditions for steady-state simulations.

Figure 5 . 26 -

 526 Figure 5.26 -Auto-generate tetrahedral mesh.

Figure 5 . 27 -

 527 Figure 5.27 -Steady-state of displacement Uz (unit: mm).

Figure 5 . 28 -

 528 Figure 5.28 -Steady-state of von Mises (unit: MPa).

Figure 5 . 29 -

 529 Figure 5.29 -Steady-state of Stress XX (unit: MPa).

Figure 5 . 30 -

 530 Figure 5.30 -Steady-state of Stress YY (unit: MPa).

Figure 5 . 31 -

 531 Figure 5.31 -Numerical model for roll forming process.

Figure 5 . 32 -

 532 Figure 5.32 -von Mises stress distributions (unit: MPa).

5. 4 . 4 Figure 5 . 34 -

 44534 Figure 5.34 -Welding model and boundary conditions.

Figure 5 . 35 -

 535 Figure 5.35 -Steady-state temperature distribution.

Figure 5 . 36 -

 536 Figure 5.36 -Mesh for MRF method simulation.

Figure 5 . 37 -

 537 Figure 5.37 -Reference mesh [A] and structured mesh [B] for material motion simulation method.

Figure 5 . 38 -

 538 Figure 5.38 -Residual stresses distributions for Lagrange simulation, MRF method simulation, material motion simulation.

Figure 5 . 39 -

 539 Figure 5.39 -Residual stresses in welding direction for MRF-Q1P0, Nodal-mesh(A), Nodal-mesh(B), and Setpbystep-mesh(A).

Figure 5 . 40 -

 540 Figure 5.40 -Residual stresses in depth directopn for MRF-Q1P0, Nodal-mesh(A), Nodal-mesh(B) and Setpbystep-mesh(A).

Figure 5 . 41 -

 541 Figure 5.41 -Free tetrahedral mesh (D), mesh section.

Figure 5 . 42 -

 542 Figure 5.42 -Stress-yy simulated by step by step Lagrangian formulation and material motion simulation by first-order interpolation (E) and by second-order interpolation (F).

Figure 5 .

 5 Figure 5.45, and figure 5.46 provide the stresses distribution at different moments. We can clearly observe that material flows in the mesh. Only one turn has been simulated

Figure 5 . 43 -

 543 Figure 5.43 -Circonferentiel welding model, boundary condition and mesh.

Figure 5 . 44 -

 544 Figure 5.44 -Steady-state temperature distribution.

Figure 5 .

 5 Figure 5.45 -von Mises stress distribution as function of time.

Figure 5 . 46 -

 546 Figure 5.46 -Stress axial σ rr distribution as function of time.

Figure 5 . 47 -

 547 Figure 5.47 -Numerical model for 3D rolling simulation.

Figure 5 .

 5 Figure 5.48 presents the stress of von Mise. We can see that the stress became steadystate after t = 60s. The steady-state mesh at t = 100s is exported for temperature calculation. The geometry is obtained by updating the coordinates of nodes (taking the displacements of nodes into account).

Figure 5 . 48 -

 548 Figure 5.48 -Stress of von Mise (MPa) vs time.

Figure 5 . 49 -

 549 Figure 5.49 -The temperature simulated at t = 100s and von Mises at t = 200s.

Figure 5 . 50 -

 550 Figure 5.50 -The temperature simulated at t = 200s and von Mises at t = 300s.

2 .

 2 APPRO 2 can improve the precision of variables transfer in case of free meshes (as we can see in figure 5.28, figure 5.29 and figure 5.30). But APPRO 2 cannot remedy the numerical diffusion problem if the too strong stress gradient exist such as in the welding simulation example (see figure 5.42).

Figure 6 . 1 -

 61 Figure 6.1 -σ 33 for welding simulation.

Figure 6 . 2 -

 62 Figure 6.2 -Boundary conditions and numerical model for 2D machining simulation.

Figure 6 .

 6 Figure 6.2 shows the numerical model and boundary conditions. A pre-crack is defined and the simulation is carried out in two steps (figure 6.2):

Figure 6 . 3 -

 63 Figure 6.3 -Modelling principle for machining simulation (displacement UY).

Figure 6 . 4 -

 64 Figure 6.4 -Numerical model and boundary conditions for 3D machining simulation.

Figure 6 . 5 -

 65 Figure 6.5 -Tetrahedral mesh for machining simulation.

Figure 6 . 6 -

 66 Figure 6.6 -von Mises' stress vs time for machining simulation.

  .

4. 1

 1 Selective reduced integration scheme in 2D. . . . . . . . . . . . . . . . . . . . . 4.2 P1/P1 and P1+/P1 scheme in 2D. . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Geometrical objects used for computation [89]. . . . . . . . . . . . . . . . . . . 4.4 The schematic of NS-FEM and SNS-FEM in 2D and 3D[34]. . . . . . . . . . . . 4.5 Definition of nodal domain in 2D and 3D. . . . . . . . . . . . . . . . . . . . . . 4.6 Couplings between node and its neighbors. . . . . . . . . . . . . . . . . . . . .

  5. The conclusions and perspectives are presented in chapter 6.

	Chapter 2
	State of the art

Contents 2.1 State of the art . . . . . . . . . . . . . . . . . . . . . . . . . 2.1.1 Lagrangian / updated Lagrangian formulation . . . . . . 2.1.2 Eulerian method . . . . . . . . . . . . . . . . . . . . . . . 2.1.3 Displacement-based reference frame formulation . . . . 2.1.4 Moving Reference Frame method with integration of constitutive equations along the streamlines . . . . . . . . . 2.1.5 Arbitary Lagrangian-Eulerian method . . . . . . . . . . . 2.1.6 Mixed Eulerian-Lagrangian method . . . . . . . . . . . . 2.2 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Table 2 .

 2 1 -Influence of elasticity on computation time[START_REF] Hacquin | Modelisation thermomecanique tridimensionnelle du laminage couplage bande/cylindres[END_REF]

			Newton-Raphson
	material behavior	iterations number	interations num-	time
			ber
	Viscoplasic	20	106	63 min
	elasto-viscoplasic	20	2135	9 h 10 min

Table 4 .

 4

		2 -Maximum displacement in Z direction
	element type	integration point	U z	∆z	∆z |U z |
	Q2 element	5600	-3.95	reference	reference
	Q1P0 element	18400	-3.88	-0.07	-0.0177
	P1/P1 element	13556	-3.38	-0.57	-0.144
	Nodal integra-
	tion based finite	3075	-4.11	0.16	0.04
	element

Table 4 .

 4 3 -Welding process parameters TG4[START_REF] Xu | Numerical simulation of a single bead on plate and three pass slot welds in austenitic stainless steel[END_REF].

		Puissance	Energy	Traveling	Bead	Interpass
		(W )	(J/mm)	speed	length	tempera-
				(mm.s -1 )	(mm)	ture ( • C)
	Pass 1	1650	1300	1.27	74	22
	Pass 2	1463	1150	1.27	76	58
	Pass 3	1388	1100	1.27	82	60

Table 4 .

 4 4 -Meshes and elements used

	Element & mesh	Mesh A	Mesh B	Mesh C
	Q1P0 element	Yes	-	-
	P1P1 element	-	Yes	Yes
	Nodal integration	-	Yes	Yes

Table 4 .

 4 5 -Computing times and hardware resources for pass 1

	Type of simulation &	CPU time (h)	RAM (Go) Disk space per
	resources			time step (Mo)
	Q1P0 mesh (A)	0.96	0.47	46.7
	P1P1 mesh (B)	2.27	0.51	29.9
	P1P1 mesh (C)	4.02	0.76	41.5
	Nodal mesh (B)	2.07	0.97	7.1
	Nodal mesh (C)	3.37	1.5	9.9

Table 5 .

 5 1 -Material properties for rolling simulation

	material properties	Values
	young's modulus	144,000 Mpa
	Poisson's ratio	0.3
	elastic limit	365 Mpa
	hardening parameter	see	Figure
		5.11

Table 5 .

 5 2 -Computation time to achieve steady-state and CPU time

	Simulations	Computating time	CPU time
		(time step)
	Lagrange step by step	170s (≈0.06s)	4246.22
	Proposed method with REF mesh	180s (0.2s)	1148.65
	Proposed method with STR mesh	190s (0.2s)	973.02

Table 5 .

 5 3 -Computation time and CPU time for FREE mesh

	Simulations	computating time	CPU time
		(time step)
	Lagrange step by step	170s (0.1s)	3780.27
	proposed method with APPRO 1	190s (0.2s)	1630.5
	proposed method with APPRO 2	180s (0.2s)	1428.9

Table 5 .

 5

		4 -Number of integration points
		Q1P0	Nodal integration	P1/P1
			based
	In rolling direction 241	162	161
	In depth	12	11	10
	In width	20	19	18

Table 5 .

 5

	5 -Comparisons of CPU time, hardware resources
	Simulations	RAM (Go) disk space (Go) CPU time
	Q1P0 element	0.43	1.5	17258
	P1/P1 element	1.2	2.2	63324
	Nodal-integration-based element	2.0	0.31	19569

Table 5 .

 5

	6 -Comparisons of CPU time, hardware resources
		RAM (Go)	CPU time
	Mesh (B) -step by step	2.0	19569
	Mesh (B) -APPRO 1	2.0	14125
	Mesh (C) -APPRO 1	1.1	11595
	Mesh (C) -APPRO 2	1.1	9554

Table 5 .

 5 7 -Computing time and hardware resources

	Type of simulation &	CPU time	RAM (Go) Disk space per
	resources			time step (Mo)
	MRF-Q1P0	3286	0.657	91.4
	MEL-mesh (A)	41034	3.6	21.04
	MEL-mesh (B)	13679	1.5	9.95
	stepbystep-mesh (A) 43647	3.6	21.24

Table 5 .

 5 8 -Material properties for tool

	Temperature	20	100 300 500 700 900 1100
	(°C)
	conductivity	0.11 0.105 0.098 0.090 0.082 0.075 0.066
	[W/(mm.K)]
	Heat capacity	220 365 290 320 331 337 338
	[J/(kg.K)]

Table 5 .

 5

		9 -Material properties for workpiece
	Temperature (°C)	20	100 300 500 700 900 1100
	conductivity	0.012 0.013 0.016 0.019 0.022 0.026 0.029
	[W/(mm.K)]
	Heat capacity*density	275 293 337 381 424 464 504
	[ * 10 -5 J/(mm 3 .K)]
	Parameter K	1207.9 1025 855 616 408 228.9 85

Table 6 .

 6 1 -Work hardening relationship

	plastc deforma-	0	0.01 0.02 0.03 0.04 0.05 0.1	0.2	0.3
	tion
	Hardening	229 251 264 274 284 292 329 388 436
	stress

Each layer of mesh has the same length, which is similar to the mesh used in the Lagrangian formulation. The mesh refined only near the source position is not possible.

This method can solve the problem of translatory motion. The rotational and helical motion is impossible.

the nodes of the element that contains the preceding points and the nodes of neighbor's element.

APPRO 3.xx: the preceding and subsequent point technique with the first-order interpolation for both displacement and the other internal variables. xx represents the fraction of the subsequent point technique in this mixed approach.

Remerciements

in Figure 4.5 by the node p, the middle of one edge connected to node p, the center of gravity of one face and the center of gravity of the tetrahedron. Regarding the definition of nodal strains, Chen et al. [START_REF] Chen | A stabilized conforming nodal integration for galerkin mesh-free methods[END_REF][START_REF] Chen | Non-linear version of stabilized conforming nodal integration for galerkin mesh-free methods[END_REF] propose to calculate them on average over the volume Ω p . By applying the generalized divergence theorem, the authors show that these strains can be expressed from an integral on the surface (or contour in 2D) encompassing the domain Ω p .

Numerical simulation of steady states associated with thermomechanical processes

Abstract

In the numerous thermomechanical manufacturing processes such as rolling, welding, or even machining involve either moving loads with respect to the fixed material or moving material with respect to fixed loads. In all cases, after a transient regime which is generally quite short, the thermal, metallurgical, and mechanical fields associated with these processes reach a steady state. The search for these stationary states using the classical finite element method requires the implementation of complex and expensive models where the loads move with respect to the material (or vice versa).

The steady-state simulation in one increment has been the subject of much researches over the past thirty years. Methods are now available and some are integrated into calculation codes commercial. Thus, a so-called Moving Reference Frame method proposed by various authors is available in the SY SW ELD T M software. This method makes it possible to calculate the steady-state of thermal, metallurgical, and mechanical states associated with a welding process, by solving a thermal diffusion-convection problem in thermal-metallurgy and by integrating, in mechanics, the constitutive equations of the material along the streamline. Moreover, this method has been used successfully in many applications, it nevertheless has some limitations. Thus the mesh must be structured and the convergence of computations is generally quite slow.

In this thesis, we propose to solve the mechanical problem in a frame linked to the solicitations, by relying on a finite element calculation method based on nodal integration and the SCNI (Stabilized Conforming Numerical Integration) technique. This method allows the use of tetrahedron meshes (or 2D triangles) without encountering a locking problem resulting from the plastic incompressibility associated with the von Mises plasticity criterion. Rather than directly calculating the steady-state, the general idea here is to construct the steady-state from a transient analysis by bringing material step by step upstream and by making it exit downstream of a fixed mesh related to the solicitations and of the limited mesh size. The (pseudo-) steady-state is therefore only achieved after certain steps of analysis. The advantages of this method lie in:

1. The use of a tetrahedral or triangles mesh (instead of structured meshes), 2. A refined mesh is only needed in the area where the solicitation locates, 3. The method proposed shows the robustness and good convergence of nonlinear mechanical calculations because of the resolution of the transient states.

Apart from a general introduction (Chapter 1) and a state of the art on the existing methods (Chapter 2), we present an approach of simulation of the movement of material within the framework of the classical finite element method on a welding problem (Chapter 3). We also provide relevant thermal boundary conditions for directly calculating the steady-state of temperature distribution.

The finite element method based on the nodal integration technique is then described in Chapter 4. The advantages and disadvantages of the method are discussed. The nodalintegration-based finite element is validated by comparing its simulation results with classical finite element methods in large elastoplastic strains, a bending problem, and a thermomechanical simulation of welding.

The nodal-integration-based finite element is then developed and applied to simulate material motion (Chapter 5). Three types of movement are considered: translational, circular, and helical. Different methods of field transport are approached and discussed as well as thermomechanical coupling. Examples of applications in the field of rolling and welding for different movements of material show the efficiency of the developed method.

Perspectives for this work are presented in Chapter 6. The envisaged perspectives aim, on the one hand, to improve the proposed method and on the other hand, to develop the