These last years spent in Nantes as a bachelor, master, and PhD student have been a real pleasure and will leave me full of good memories.

First, I would like to thank the directors of this thesis, Laurent Stainier and Thomas Heuzé, for trusting me for these research. They both guided me very patiently throughout this work, especially during the reviewing of manuscripts. In addition, I must emphasize the very good musical tastes of Laurent Stainier and the gastronomic tastes of Thomas Heuzé, which have been the subject of interesting discussions. Owing to their moral values and scientific skills, I have been very happy to work with them.

Second, I am grateful to several faculty members met when I was a master student at Nantes University, who then became colleagues: Jean-Christophe Thomas, Mathilde Chevreuil, Marc François, Anh Le van, Rabah Bouzidi and Laurent Gornet. Their courses showed me how much computational mechanics can be interesting though challenging.

Third, I also want to thank my "building F colleagues": Baptiste and Erwan (along with the other CBEM organizers), Pierre, Rohit, Jorge, Benoît, Romain, Simon, George, Tauno, Abdullah, Zoltán, Cécile, Daria, Quentin, Steven, Fabien, Hong-Duc, Elie, Rhaphaël, Adrien and Ataollah, as well as the "building T people": Julie, Marie and Nicolas. They all contributed to an atmosphere which helps during difficult times. This reminded me the good old days at the university with Jan, Cocaud, Laporte, Le Mouel, Macé, "Le Commodore" Varin and Claude.

Fourth, I am thankful to Cynthia and Marine who managed the administrative difficulties that I encountered. Moreover, I thank Anthony who saved me by lending me his laptop the day of my defense.

Fifth, I would like to thank Morgane for considerably supporting me and for undergoing my mood swings. I also must thank Dorothée, Hanaë and Gaël for lending me their house so that I could write peacefully.

At last, I am grateful to my friends and family for coming to my defense and to the celebration of my birthday the next day. They all supported me during these three years, especially my sister Emmanuelle by the end of her pregnancy.

Chapter 1 Introduction 1.1 General context

The work presented in this thesis is concerned with hyperbolic systems of conservation laws, used in mechanics to describe the propagation of mechanical waves. A wide variety of engineering problems including among others acoustics, aerodynamics or impacts, which are of major importance for plenty of applications, are modeled with this class of mathematical equations. Although these applications may involve solid and fluid media whose constitutive responses differ, the focus is here on solid mechanics. More specifically, an important class of models is that of dissipative solid behaviors depending on the loading history undergone. Such materials are implicated in high-speed forming techniques, crash-proof design or in the reliability of structures.

In these solids the waves propagate while carrying the information about the loading, interacting with each other and reflecting on the boundary so that complex structures arise. The accurate assessment of irreversible deformations in dissipative solids therefore requires the correct description of those waves as well as the ability to account for their interactions. Moreover, the time scale governing the propagation of waves may be different from that of other phenomena also involved in a deformation, as for viscous effects in solids for instance. At last, the problem may be further complicated by possibly large displacements, rotations or strains, undergone by the solid.

Given the complexity of the equations, the numerical simulation provides a framework for approximating the solutions of a model. It is a way to virtually experiment and highlight phenomena that are implicitly described by a model but not necessarily observed. The simulation can then be seen as a way to make the theory explicit. Furthermore, owing to feasibility reasons or to the limited amount of information that can be measured by instrumental devices, experimental works leave some gaps which can be filled by numerical investigations. Finally, the simulation allows understanding the models and therefore, to better design so that subsequently, a given problem is addressed as well as possible. Nevertheless, numerical methods are based on parameters that require interactions with experimental works, especially for the constitutive models in solid mechanics that become more and more complex. Furthermore, non-physical parameters must sometimes be defined by the user in order to obtain results that are close to a phenomenon one is focusing on.

Naturally, the numerical simulation is not a straightforward undertaking, particularly in solid mechanics for the reasons mentioned above. More specifically, the solution of hyperbolic problems in finite deforming dissipative solids is still an open scientific issue. Indeed, the combination of large deformations with a wave structure can lead to complex loading paths that are not even observed at the model level. Most of the time, these aspects are not directly tackled in constitutive numerical integrators, which however provide rather satisfactory results.

Numerical methods for hyperbolic problems in solid mechanics

The numerical simulation of initial boundary value problems (IBVP), to which hyperbolic ones belong, has been and is still widely performed in solid mechanics with the Finite Element Method (FEM) [START_REF] Belytschko | Nonlinear finite elements for continua and structures[END_REF]. This approach is based on the discretization of a computational domain into a set of elements with simple geometries used to solve the equations. An approximate solution is thus built by means of a combination of polynomial functions whose degree defines the approximation order. FEM became attractive due to its ability to handle low or high-order approximations, and to easily deal with complex geometries and nonlinear constitutive models as those mentioned above. Nevertheless, difficulties may be encountered if large deformations occur when the finite elements deform with the domain, according to a material description of the motion (Lagrangian approach). Indeed, the method is less efficient and accurate when the elements are highly distorted or entangled so that re-meshing techniques and projection steps must be employed. These issues can be avoided by using a spatial description of the motion consisting in viewing the elements vertices as fixed points of the space (Eulerian approach). However interface tracking techniques and diffusing convection steps are required in order to follow the boundaries and transport internal variables which is less convenient for solid mechanics. Alternatively, arbitrary Lagrangian Eulerian (ALE) methods aim at meeting advantages of both approaches by freeing itself of their respective limits by distinguishing the motion of the mesh to these of material points. This type of strategy nonetheless also requires re-meshing or re-zoning procedures that can be costly for fine meshes. They also require diffusive projection steps of internal variables for solid media. In addition to problems caused by finite deformations, classical explicit time integrators used in solid dynamics with FEM introduce high frequency noise in the vicinity of discontinuities. Such regions of high gradient may be caused by the waves arising in the solutions of hyperbolic problems. The removal of these spurious numerical oscillations with additional viscosities is difficult to achieve without loss of accuracy, and can be troublesome for the wave tracking.

On the other hand, the Finite Volume Method (FVM) [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF], initially developed for fluid dynamics, provided until the 90s piece-wise constant or piece-wise linear approximate solutions in cells that discretize a continuum. The extension to very high-order has since been proposed by increasing the stencil of the scheme (see WENO [START_REF] Liu | Weighted essentially non-oscillatory schemes[END_REF]). This class of methods can embed tools based on the Total Variation Diminishing (TVD) stability condition [START_REF] Harten | High resolution schemes for hyperbolic conservation laws[END_REF], thus ensuring that no numerical spurious oscillations arise in the solutions. The formulation moreover lies on a conservative form leading to the same order of accuracy for all components of the unknown vector U. In particular, one shows that both velocity and gradients arise in U in solid mechanics. This point makes a big difference with respect to methods that do not use a formulation written as a differential system of order one, namely FEM, for which the convergence rate of gradients is one order less than that of displacement. To some extent, the writing of solid mechanics equations in the form of conservation laws amounts to a mixed approach, well-known in FEM. FVM formulations moreover rely on numerical fluxes that enable to account for the characteristic structure of hyperbolic equations. Hence, finite volumes enable an accurate tracking of the path of waves although the most widely used approximations are second-order. Recent studies furthermore extended these approaches to solid mechanics for problems involving history-dependent models [START_REF] Favrie | Dynamics of shock waves in elastic-plastic solids[END_REF][START_REF] Heuzé | Lax-Wendroff and TVD finite volume methods for unidimensional thermomechanical numerical simulations of impacts on elastic-plastic solids[END_REF], or finite deformations with a Lagrangian formulation [START_REF] Lee | Development of a cell centred upwind finite volume algorithm for a new conservation law formulation in structural dynamics[END_REF][START_REF] Haider | A first-order hyperbolic framework for large strain computational solid dynamics: An upwind cell centred total lagrangian scheme[END_REF]. The latter methods are nevertheless grid-based techniques for which the numerical difficulties linked to mesh occur.

The discontinuous Galerkin (DG) approximation [START_REF] Reed | Triangular mesh methods for the neutron transport equation[END_REF] makes possible to build numerical schemes that benefit from both FEM and FVM. In DG-methods, the approximate solution is sought as a combination of piece-wise polynomial functions whose supports are dictated by the discretization used. Therefore, the combination of the DG approximation with the finite element formulation (DGFEM) yields a local high approximation order, the same order being achievable for both velocity and gradients if a conservative form is used. Furthermore, numerical fluxes at the interface between elements, which enable to take into account the characteristic structure of hyperbolic systems, arise from the introduction of DG approximation. Several development steps, aiming at removing the numerical noise appearing in DGFEM solution, followed the works on FVM in order to reach Total Variation Diminishing in the Means (TVDM) and Total Variation Bounded (TVB) high order schemes that ensure convergence to entropy-satisfying solutions [START_REF] Cockburn | Discontinuous galerkin methods for convection-dominated problems[END_REF]. However, in spite of the piece-wise continuous approximation, the methods fail to accurately capture discontinuities owing to the time discretization. Adopting a similar approach, the Time-Discontinuous Galerkin (TDG) [START_REF] Hughes | Space-time finite element methods for elastodynamics: formulations and error estimates[END_REF] and later, the Space-Time Discontinuous Galerkin method (SDG) [START_REF] Richter | An explicit finite element method for the wave equation[END_REF], relaxed the continuity of fields in the time domain. By discretizing the entire space-time domain as a possibly unstructured mesh, SDG avoids the difficulties related to the time integrators and hence, enables the following of waves. Nevertheless, although these approaches can easily handle meshadaption strategies due to the relaxation of fields continuity, it does not eliminate the mesh tangling problems for Lagrangian formulations [START_REF] Vilar | A high-order Discontinuous Galerkin discretization for solving two-dimensional Lagrangian hydrodynamics[END_REF].

In order to address the loss of accuracy caused by mesh distortion, another class of numerical methods based on a space discretization using a collection of points has been developed in parallel to the above approaches [START_REF] Belytschko | Meshless methods: An overview and recent developments[END_REF][START_REF] Nguyen | Meshless methods: A review and computer implementation aspects[END_REF]. In contrast to finite volumes or finite elements, mesh-free methods represent a spatial domain by means of a collection of points that are given a support allowing them to interact with each other. A wide variety of mesh-free methods such as the Smoothed-Particle Hydrodynamics [START_REF] Gingold | Smoothed particle hydrodynamics: theory and application to nonspherical stars[END_REF] or the Element-Free Galerkin [START_REF] Belytschko | Element-free galerkin methods[END_REF], have thus been constructed.

Particle-In-Cell methods (PIC) [START_REF] Harlow | The particle-in-cell method for numerical solution of problems in fluid dynamics[END_REF] are, on the other hand, based on particles that move in a computational mesh while carrying the fields of a problem. The underlying grid is used in order to compute an approximate solution that is projected and stored at particles. Hence, the background mesh can be discarded at each time step and re-constructed for computational convenience. The application of PIC to solid mechanics led to the Material Point Method (MPM) in which the constitutive equations are managed at particles. As a result MPM can be seen as a mesh-free extension of FEM with quadrature points moving in elements. Overcoming the storage of the approximate solutions based on elements or cells enables to remove mesh entanglement instabilities. It is nevertheless well-known that PIC exhibits numerical dissipation that can be reduced at the cost of spurious oscillations [START_REF] Burgess | Mass matrix formulation of the flip particle-in-cell method[END_REF].

In light of the above brief overview, it appears that in spite of the plenty of existing numerical methods, none can easily deal with all the difficulties involved by hyperbolic problems in finite deforming dissipative solids. First, the large deformations often met in solid mechanics problems make tricky the employment of mesh-based approaches. Second, the waves propagating in media can be accurately followed providing that the scheme used computes solutions devoid of spurious oscillations and too much numerical diffusion. At last, accounting for the characteristic structure of the solution of hyperbolic problems within a numerical method yields results closer to the expected output of the model. That is, the solution of the wave structure is in general not directly tackled in solid mechanics solvers in such a way that the approaches are not devoted to the accurate description of waves.

Approach followed in this work

Point of view of this thesis

This work is based on the idea that the accurate numerical solution of hyperbolic problems relies on robust and efficient discretization techniques and the ability to embed information about some particular solution of the model into the numerical scheme. The underlying concepts arise from the Godunov method [START_REF] Godunov | Finite difference method for numerical computation of discontinuous solutions of the equations of fluid dynamics[END_REF] that first proposed to account for the characteristic structure of the solution of hyperbolic problems within numerical schemes. This approach, as we shall see later, enables significant improvements over some other methods. Although the amount of information provided is limited by the computational cost, new model reduction tools can now allow to "compress" the amount of information and, in turn, to reconsider the use of more complex Riemann solvers in order to improve the accuracy of numerical solutions.

Two scientific objectives have therefore been pursued in the present work:

• the development of a promising discretization for solid mechanics problems involving finite deformations

• the identification of the response of two-dimensional elastic-plastic solids to dynamic step loading

The strategy adopted

First, the numerical scheme developed here provides a material description of the motion so as to handle historydependent constitutive models while avoiding the shortcomings of ALE and Eulerian methods. In addition, in order not to suffer from mesh entanglement, it seems better to turn to a mesh-free method. This method is wanted to mimic the physics of hyperbolic problems by accounting for their intrinsic structure.

It is therefore proposed here to mix the above features by extending the material point method to the Discontinuous Galerkin approximation. The Discontinuous Galerkin approximation provides an appealing framework for describing moving discontinuities such as waves propagating in solids, and the potential ability of increasing approximation order. Moreover, the numerical fluxes naturally arising from the use of DG approximation introduces the characteristic structure on the one hand, and takes advantage of the work done in the context of finite volumes on the other hand. Furthermore, the use of an arbitrary grid motivates the choice of the MPM due to the convenience it allows in computing the numerical fluxes at the interfaces between finite elements by means of an approximate Riemann solver [START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics: a practical introduction[END_REF]. A balance between diffusion and oscillations exists in PIC methods in such a way that there is some freedom in regard to the MPM setting used. It is however thought here that the use of DG approximation leads to a reduction of the influence domain of the particles and hence, to a reduction of the numerical diffusion. Therefore, it is preferred here to avoid oscillations at the risk of introducing diffusion, which should be limited by the DG approximation. As a first development step of the Discontinuous Galerkin Material Point Method, we restrict our attention to space-DG. At last, particular attention is paid here to discontinuous solutions, the extension of the method to higher-order approximation for regular ones will be the purpose of future works.

The approach described previously should be able to account for the characteristic structure of hyperbolic problems. Indeed, the point of view adopted here is that a numerical scheme can properly mimic the solution providing that a sufficient amount of information about the model is available. Nonetheless, it is not the case for all the constitutive models considered in this work. More specifically, while the response of one-dimensional elastoplastic solids is well-known, there are some lacks in describing the behavior of such materials in more dimensions. As a consequence, the characteristic structure of the solution of hyperbolic problems in elasticplastic solids under small strains is also investigated in the present work. This preliminary work is expected to pave the way for the solution of similar problems within the large deformations framework.

Organization of the manuscript

This manuscript starts with a review of some basics about general hyperbolic problems in chapter 2. Then, the governing equations of solid mechanics including balance laws and constitutive models are recalled so that hyperbolic systems of solid mechanics are written. Thus, applying the tools of characteristic analysis, it will be seen that wave solutions naturally arise from the equations. In particular, the exact solutions of Riemann problems in linear elastic and elastic-plastic solids are recalled and a plane wave solution in a hyperelastic medium is developed. Given the complexity in computing the exact solution of Riemann problems for nonlinear problems, a well-known approximate Riemann solver is at last presented.

A historical review of PIC methods is made and both Eulerian and Lagrangian formulations of the MPM are derived in 3. After illustrating some shortcomings of the latter scheme on a simple test case, the extension of the MPM to the discontinuous Galerkin approximation is developed. The DGMPM is then derived with a total Lagrangian formulation. As we shall see, this new approach enables the use of fractional-step methods for non-homogeneous systems as well as Riemann solvers for the computation of intercell fluxes. The remainder of the chapter concerns the numerical analysis of the DGMPM in terms of stability and accuracy.

Chapter 4 is devoted to the comparison of the performance of the method with other existing schemes and exact solutions. To begin with, problems falling in the linearized geometrical framework will be considered in one and two space dimensions. More specifically, the constitutive models assumed include elasticity, elasto-viscoplasticity and elastoplasticity. Then, hyperbolic problems in one-dimensional and two-dimensional hyperelastic media are looked at.

The simulations performed on elastic-plastic solids in chapter 4 emphasize that the numerical solutions can be improved by using approximate elastic-plastic Riemann solvers, thus introducing the understanding one has of the physics into the scheme. The purpose of chapter 5 is therefore to give a contribution to the solutions of hyperbolic problems in two-dimensional elastic-plastic solids undergoing dynamic step-loadings. In particular, the understanding of the structure of such problems, along with the study of loading paths in simple waves, are addressed. Indeed, though a lot of numerical results of this class of problems may be found in the literature, the structure of the solutions is rather unknown, as shown by the survey in the chapter. Thus, the proposed characteristic analysis of plane strain and plane stress problems highlights the combined-stress wave nature of the solutions. Then, the mathematical study of loading paths followed across plastic waves, supplemented by numerical results, emphasizes typical responses of solids to given loading conditions. Finally, the idea of benefit from the identified stress paths in numerical schemes by means of a dedicated Riemann solver is discussed.

Chapter 2

Hyperbolic partial differential equations for solid dynamics Introduction

In this chapter, generalities about systems of Partial Differential Equations (PDEs) are first introduced. More specifically, the notions of characteristics and hyperbolicity for PDEs are introduced in section 2.1 for first-order systems, which will be of particular interest in the remainder of the manuscript. Then, introduction of balance laws of solid dynamics and derivation of constitutive equations from the thermodynamics in section 2.2, lead to first-order hyperbolic PDEs. By using tools introduced in section 2.1, the characteristic analysis of these systems is carried out in section 2.3 in order to derive exact solutions of particular problems in section 2.4. These solutions allow the highlighting of different types of waves: (i) discontinuous waves, governed by the Rankine-Hugoniot jump condition, within one-dimensional linear elastic and elastic-plastic media (ii) shock waves, also following the Rankine-Hugoniot condition, and simple waves, within a non-linear problem (one-dimensional strain state in a Saint-Venant-Kirchhoff hyperelastic medium). At last, strategies enabling the computation of approximate solutions of non-linear problems are reviewed in section 2.5.

Generalities -Hyperbolic partial differential equations

A wide variety of physical phenomena are modeled in the space-time domain by partial differential equations. The purpose of this section is to review generalities about PDEs and suited strategies depending on the nature of the equations so that solutions can be derived.

General concepts

A system of partial differential equations can be written by means of a vector operator G of independent and dependent variables (x 1 , ..., x N ) and (u 1 , ..., u I ):

G x 1 , ..., x N , u 1 , ..., u I , ∂u 1 ∂x 1 , ..., ∂ M u I ∂x M N = 0 (2.1)
The dimension of the system is given by the size I of the array U T = [u 1 , ..., u I ] ∈ R I , referred to as the unknown vector. The highest derivative of the unknown vector in the system defines the order of the system M .

In equation (2.1) and in what follows, sans-serif symbols refer to matrices while calligraphic symbols stand for column arrays. Furthermore, the partial derivatives of a quantity u with respect to a variable x may be written u x when there is no ambiguity. Making use of index notation and the convention of implicit summation over repeated indices, a system of partial differential equations reads: Coefficients matrices A p and the vector S may depend on independent variables and the unknown vector (x 1 , ..., x N , U) leading to different types of partial differential systems. Namely, whether those terms are functions of the x k or not leads respectively to linear systems with variable coefficients or to linear systems with constant coefficients. The system remains linear if S depends linearly on U, and is semi-linear if the relation is non-linear. Finally, if A p depends on the vector U and its derivatives up to order M -1, the system is called quasi-linear.

N k=1 M p=1 A p ij ∂ p U j ∂x p k + S i = 0
The Cauchy problem consists in finding a solution U of system (2.2) that satisfies a set of given prescribed values. Geometrically speaking, the solution of such a problem can be seen as the building of a hyper-surface of R I+N , hence the term of integral surface for U. Such a problem can be reduced to that of solving a first-order Cauchy problem by using suitable changes of variables [22, p.54], we will therefore focus on first-order PDEs.

Notion of characteristics -Hyperbolic problems

The theorem of Cauchy-Kowalewski locally ensures the existence of solutions of a Cauchy problem for partial differential systems and is based on the restrictive requirement of analytic coefficient matrices and initial data (see [22, p.46]). The case of first-order systems however, only requires continuity and differentiability conditions and is based on the concept of characteristics, which makes the development of a solution more flexible.

First-order quasi-linear equations

To illustrate the aforementioned notions, we consider the first-order quasi-linear PDE with independent variables x and t: au x + bu t = c (

where coefficients a and b are such that a 2 + b 2 = 0. Given values of u are prescribed along a curve defined by means of a parameter η in the (x, t) plane as C 0 : (x(η), t(η)), so that u(x(η), t(η)) draws a curve C in the space (x, t, u). We assume that C 0 is regular, namely dx dη 2 + dt dη 2 = 0, and that one of the derivatives, say dt dη , does not vanish. u(x, t) that contains the initial curve and satisfies (2.3). Thus, one seeks the partial derivatives u x and u t of u on C in order to extend the given data into a strip in the neighborhood of the initial curve. The total derivative of u along C being: du dη = u x dx dη + u t dt dη one gets, after multiplying the previous equation by dη, the following relation between the partial derivatives of u on the initial curve:

u t = du dt -u x dx dt Then, equation (2.
3) can be rewritten as:

(a -b dx dt )u x = c -b du dt t∈C0 (2.4)
The right-hand side of equation (2.4) is known along C so that the Cauchy problem admits a unique solution

u x if and only if: dx dt = a b (2.5)
An initial curve satisfying the condition (2.5) is a non-characteristic curve and enables to uniquely determine a solution of the Cauchy problem. On the other hand, an initial curve defined such that dx dt = a b is a characteristic curve and yields infinitely many solutions [23, p.65].

Geometrical representation of characteristic curves

Consider a partial differential equation of the form (2.3), and prescribed values of u along a characteristic curve C . Since one cannot find a unique solution of the Cauchy problem in this case, infinitely many integral surfaces

u (i) (x, t) with normal vectors n (i) = [u (i) x , u (i) t , -1] can intersect C . Those integral surfaces satisfy equation (2.
3) and hence, n (i) • w = 0 where w = [a, b, c], so that the set of tangent planes to solutions u (i) (x, t) forms a fan whose axis is w. This situation is depicted in figure 2.2 for which an initial characteristic straight line C is contained by integral surfaces satisfying a PDE of the form (2.3).

x t u C u (1) u (2) u (3) Characteristic line elements, tangent to all integral surfaces u (i) (x, t) are then defined as:

  dx dt du   =   a b c   (2.6)
Introduction of a parameter η and integration of equation (2.6) yield a one-parameter family of characteristic curves of the PDE:

x = x(η) ; t = t(η) ; u = u(η)
Hence, a characteristic curve is tangent at every point to all the integral surfaces, and an infinity of integral surfaces cross one characteristic curve. As a consequence, if the initial curve is a characteristic curve, infinitely many integral surfaces contain it so that the Cauchy problem cannot be solved. However, the following statement holds [23, p.63]:

Theorem 1 (Courant) Every surface u(x, t) generated by a one-parameter family of characteristic curves is an integral surface. Conversely, every integral surface is generated by a one-parameter family of characteristic curves.

This theorem will be used in what follows to solve the Cauchy problem.

First-order quasi-linear systems

The concept of characteristic curves is now extended to first-order quasi-linear systems of dimension I. Consider the following system written in matrix form:

A t (x, t, U) U t + A x (x, t, U) U x + S = 0 (2.7)
Similarly to quasi-linear PDEs, given values of U are prescribed along a regular curve C 0 : (x(η), t(η)) defining an initial curve U(x(η), t(η)) of the (x, t, U) space. The Cauchy problem consists in finding all the derivatives of U(x, t) such that equation (2.7) is satisfied in the vicinity of C . Making use of the total derivative of U along the initial curve:

dU dη = U x dx dη + U t dt dη → U t = dU dt t∈C0 -U x dx dt system (2.7
) can be rewritten:

A x -λA t U x + S + A t dU dt t∈C0 = 0 (2.8)
where:

λ = dx dt (2.9)
With S and A t dU dt t∈C0 known along C 0 , the Cauchy problem admits a unique solution U x along C if the determinant of the system does not vanish, that is:

D = A x -λA t = 0 (2.10)
where D is called the characteristic determinant of system (2.7). If D does not have real roots along C 0 , the problem is said elliptic and the Cauchy problem can be solved. Indeed, in that case the knowledge of U along the initial curve allows the computation of derivatives and hence, the building of an integral strip defined by U, U x , U t . If the characteristic determinant admits I real roots on the other hand, system (2.8) can no longer be solved. Those eigenvalues come along with left and right eigenvectors respectively defined as:

L k i A x ij = λ k L k i A t ij ; A x ij R k j = λ k A t ij R k j k = 1, ..., I (2.11)
Remark 1 Note that eigenvectors can be stored as matrices R and L where R ij = R j i and L ij = L i j .

Definition 1 A first-order system of I partial differential equations is said hyperbolic if it admits real eigenvalues and I independent eigenvectors [START_REF] Courant | Methods of mathematical physics[END_REF]. For those problems, one can draw a set of one-parameter families of characteristic curves ϕ k in the (x, t) plane by integrating the relations

λ k = dx/dt (1 ≤ k ≤ I).
Example 1 Consider the first-order system with variable coefficients

x 0 0 -x ∂ ∂t U 1 U 2 + ∂ ∂x U 1 U 2 = 0 0
whose characteristic determinant (2.10) is:

(1 -λx)(1 + λx) = 0
We thus have two solutions λ 1,2 = ±1/x leading, by integration of (2.9), to two one-parameter families of characteristic curves:

t 1 (x) = 1 2 x 2 + c 1 and t 2 (x) = - 1 2 x 2 + c 2
Those curves are drawn in figure 2.3a for several values of integration constants c 1 and c 2 .

Example 2 Consider now the first-order system with constant coefficients

1 0 0 2 ∂ ∂t U 1 U 2 + ∂ ∂x U 1 U 2 = 0 0
whose eigenvalues, according to equation (2.10) satisfy

(1 -λ)(1 -2λ) = 0
Two real roots exist λ 1 = 1 ; λ 2 = 1/2, leading by integration of (2.9) to two one-parameter families of straight lines:

t 1 (x) = x + c 1 and t 2 (x) = 2x + c 2
Unlike example 1, coefficient matrices do not depend on independent variables, thus yielding characteristic straight lines in the (x, t) plane (see 2.3b).

(a) Example 1:

λ 1,2 = ±1/x x t c 1 = 0 c 1 = -1 c 1 = -2 c 1 = -3 c 2 = 5 c 2 = 6 c 2 = 7 c 2 = 8 (b) Example 2: λ 1 = 1 and λ 2 = 1/2 x t c 1 = 0 c 1 = -1 c 1 = -2 c 1 = -3 c 2 = 0 c 2 = 1 c 2 = 2 c 2 = 3

The method of characteristics

As theorem 1 states, an integral surface is generated by a one-parameter family of characteristic curves. Therefore the knowledge of those curves enables the building of the solution of the Cauchy problem. Indeed, the projection of the quasi-linear system (2.7) onto the left eigenbasis or left characteristic basis leads to:

L k A t U t + A x U x + L k S = 0
where L k satisfies (2.11), and hence:

L k A t (U t + λ k U x ) + L k S = 0
In this equation, the directional derivative of U along the kth characteristic curve ϕ k arises, namely:

dU dt t∈ϕ k = U t + λ k U x
Thus, along a characteristic curve a system of partial differential equations reduces to a system of Ordinary Differential Equations (ODEs) composed of the following characteristic equations:

L k A t dU dt + S = 0 (2.12)
Integration of equations (2.12) yields a set of integral curves from which the Cauchy problem can be solved. Indeed, the solution at a point of the (x, t) plane can be determined by tracing backward the characteristic curves to the initial curve and integrating ODEs (2.12) along those paths according to the method of characteristics. Note that if S is zero, then U is constant along characteristic curves.

To illustrate the method, let us consider again the quasi-linear system of example 1 for which the Cauchy problem is built by prescribing initial conditions along the x-axis. Note that "initial conditions" have now a physical meaning since they are defined at t = 0, the Cauchy problem is then an Initial Value Problem (IVP). Through a point (x * , t * ) pass two characteristic curves, each belonging to a different one-parameter family. The solution at this point can be determined by integrating the ODE corresponding to the first (resp. second) eigenvalue of the system between (x 1 , 0) (resp. (x 2 , 0)) and (x * , t * ). The singularity of hyperbolic problems can hence be circumvented by using the characteristic structure in order to determine a unique solution. We see that

x t λ 1 λ 2 x 1
x 2 Domain of dependence (x * , t * ) Figure 2.4: Domain of dependence of the solution at (x * , t * ) for the system of example 1. only a segment of the initial curve has an influence on the solution at a given point. Namely, the intersections of the initial curve and characteristic curves with the highest and the lowest slopes define the domain of dependence of the solution at this point (see figure 2.4). This property of hyperbolic problems implies the existence of waves that propagate information at finite speeds corresponding to the eigenvalues of a quasi-linear form. The theory presented so far will be applied in what follows to solid mechanics.

Governing equations of solid mechanics

The mathematical laws describing the deformation of a solid body are summarized in this section. First, the kinematic laws governing the motion of each material point belonging to a solid are considered. Then

GOVERNING EQUATIONS OF SOLID MECHANICS

strain measures are associated with internal forces through the thermodynamic framework so that constitutive equations are derived. For a more exhaustive review of governing equations, see for instance [START_REF] Belytschko | Nonlinear finite elements for continua and structures[END_REF][START_REF] Marsden | Mathematical foundations of elasticity[END_REF][START_REF] Truesdell | The non-linear field theories of mechanics[END_REF][START_REF] Simo | Computational inelasticity[END_REF].

Kinematic laws -Strain measures

Consider a three-dimensional solid with volume denoted by Ω ⊂ R 3 bounded by the surface ∂Ω. This body undergoes external forces that can either be localized on a part of the external surface of the body (i.e. surface forces) or act in the whole solid domain (i.e. volume forces). Due to the presence of such loads, the domain may change within the time interval τ = [0, T ] and is hence written as a function of time Ω(t) (t ∈ τ ). The state of the solid at time t = 0, corresponding to a non-deformed state with volume Ω(t = 0) = Ω 0 , is referred to as the initial configuration. Some problems require the use of a reference configuration that can be deformed and to which equations are referred. In what follows, the reference and initial configurations are identical. At a given time t > 0, the volume is Ω(t) = Ω t and the state of the solid corresponds to the current configuration. These configurations are depicted in figure 2.5. The points of R 3 are located with Eulerian or

E 1 , e 1 E 3 , e 3 E 2 , e 2 Ω 0 Ω(t) • X • φ(X, t) Figure 2
.5: Deformation of a solid body between a reference state Ω 0 to a subsequent state Ω t . spatial coordinates x = x i e i while material particles in the reference configuration are located with Lagrangian coordinates X = X α E α . At time t, the particle initially located at X may have moved to a different position given by the smooth mapping φ(X, t) = φ i (X, t)e i , providing the path of every particle of the solid during the deformation. Thus, the current volume of the solid is defined by means of Eulerian coordinates and the deformation function φ(X, t) as: Ω(t) = {x ∈ R 3 : x = φ(X, t), X, t ∈ Ω 0 × τ }. Note that in the above definitions Greek indices are used for quantities evaluated in the reference configuration whereas Latin ones refer to quantities defined in the current configuration.

The displacement and velocity vectors of a particle at time t are respectively:

u(X, t) = φ(X, t) -X ∀ X, t ∈ Ω 0 × τ (2.13) v(X, t) = ∂φ ∂t (X, t) = φ(X, t) ∀ X, t ∈ Ω 0 × τ (2.14)
where the superposed dot denotes the material time derivative. The second-order two-point deformation gradient tensor is defined as:

F = ∇ 0 φ(X, t) (2.15)
where ∇ 0 (•) is the gradient operator in the reference configuration. This tensor can also be written by using equation (2.13):

F = ∇ 0 u(X, t) + I (2.16)
with I, the second-order identity tensor. The deformation gradient tensor characterizes the variations of lengths, angles, areas and volumes. Indeed, the infinitesimal vector, oriented surface and volume elements in the reference configuration, respectively denoted by dX, dS and dV , transform respectively to:

dx i = F iα dX α ds i = JF -1 αi dS α dv = JdV (2.17)
in the current configuration. Transport equations (2.17) involve the determinant of the deformation gradient J = det(F ) > 0, also called the Jacobian of the deformation.

Since it accounts for changes in lengths and angles (i.e. the change of shape of a body), the deformation gradient is one strain measure among others. For instance, one also defines the right Cauchy-Green and the Green-Lagrange tensors as:

C = F T F ; E = 1 2 (C -I)
respectively. Making use of equation (2.16), the Green-Lagrange tensor reads:

E = 1 2 (∇ 0 u + ∇ 0 u T + ∇ 0 u T ∇ 0 u)
In particular, when a deformation involves displacement vectors such that ∇ 0 u 1, the last term of the previous definition can be neglected, leading to:

E ≈ 1 2 (∇ 0 u + ∇ 0 u T ) = ε (2.18)
with ε the symmetric linearized strain tensor. Such deformations fall in the small strains framework and are characterized by small strains but possibly large displacements. Furthermore, when the deformation leads to a displacement vector u L 1, where L is a characteristic length of the domain, reference and current configurations are considered as identical within equations of the Initial Boundary Value Problem (IBVP). The aforementioned situations correspond to the linearized geometrical framework or infinitesimal theory.

Balance equations

The time derivative of equations (2.15) and (2.18) combined with the definition of the velocity field (2.14) yield respectively:

Ḟ -∇ 0 v = 0 ε -∇ s v = 0 (2.19)
where ∇ s (•) denotes the symmetric gradient operator. By rewriting the gradient operators as:

∇ 0 v = ∇ 0 • (v ⊗ I) ∇ s v = 1 2 ∇ • (v ⊗ I + I ⊗ v)
with ∇ 0 • (•) and ∇ • (•), the right divergence operators in reference and current configurations respectively. With these forms of gradient operators, geometrical relations (2.19) can be written as kinematic or geometrical balance laws [START_REF] Plohr | A conservative eulerian formulation of the equations for elastic flow[END_REF][START_REF] Haider | A first-order hyperbolic framework for large strain computational solid dynamics: An upwind cell centred total lagrangian scheme[END_REF][START_REF] Bonet | A first order hyperbolic framework for large strain computational solid dynamics. part i: Total lagrangian isothermal elasticity[END_REF][START_REF] Lee | A new jameson-schmidt-turkel smooth particle hydrodynamics algorithm for large strain explicit fast dynamics[END_REF]:

Ḟ -∇ 0 • (v ⊗ I) = 0 (2.20) ε - 1 2 ∇ • (v ⊗ I + I ⊗ v) = 0 (2.21)
Then, assuming that the mass of some amount of matter remains constant during the deformation, one writes the conservation of mass in integral form:

Ω ρdv = Ω0 ρ 0 dV ∀ t ∈ τ, ∀ Ω 0
which, with the third transport formula reads:

Ω0 (Jρ -ρ 0 ) dV = 0 ∀ Ω 0 (2.22)
Since equation (2.22) holds regardless of the volume Ω 0 , the integrand must vanish so that the local conservation of mass is written:

ρ (φ(X, t), t) = ρ 0 (X) J ∀ X, t ∈ Ω 0 × τ
Furthermore, Newton's second law states the equilibrium between inertia and external forces undergone by a solid Ω. In the current configuration this conservation law consists of the translational and rotational balances, also known as linear momentum and angular momentum balance equations, which are respectively:

d dt Ω ρvdv = ∂Ω tds + Ω ρbdv ∀ t ∈ τ, ∀ Ω (2.23a) d dt Ω ρx × vdv = ∂Ω x × tds + Ω ρx × bdv ∀ t ∈ τ, ∀ Ω (2.23b)
where t and b denote surface and volume forces and the cross operator denotes the vector product. The secondorder Cauchy stress tensor σ is then introduced by using Cauchy's theorem t = σ • n where n is the outward normal vector to the surface element ds.

Theorem 2 (Ostrogradski)

The divergence theorem relates the flow of a quantity through a closed surface ∂Ω to the divergence of this quantity inside the volume Ω delimited by ∂Ω:

∂Ω (•) • n ds = Ω ∇ • (•) dv ∀ Ω (2.24)

Definition 2

The Piola transform T P of a second-order tensor T is defined as:

T P = JT • F -T
and satisfies:

∇ 0 • T P = J∇ • T
The conservation of linear momentum (2.23a), combined with the volume transport theorem (2.17), reads in the reference configuration:

d dt Ω0 ρ 0 v dV = Ω0 J∇ • σ dV + Ω0 ρ 0 b dV ∀ t ∈ τ, ∀ Ω 0
or, by using definition 2:

Ω0 (ρ 0 v -∇ 0 • Π -ρ 0 b) dV = 0 ∀ t ∈ τ, ∀ Ω 0 (2.25)
where the first Piola-Kirchhoff stress tensor (PK1) Π = Jσ • F -T is the Piola transform of Cauchy stress tensor. Thus, the vanishing of the integrand in equation (2.25) yields the balance equation of the Lagrangian linear momentum:

ρ 0 v -∇ 0 • Π = ρ 0 b ∀ X, t ∈ Ω 0 × τ (2.26)
or equivalently for the current configuration:

ρ v -∇ • σ = ρb ∀ x, t ∈ Ω × τ (2.27)
On the other hand, the conservation of angular momentum (2.23b) leads to the symmetry of Cauchy stress tensor σ = σ T , or equivalently from the definition of PK1 tensor, Π • F T = Π T • F [START_REF] Marsden | Mathematical foundations of elasticity[END_REF]. We complete the set of balance laws by considering the first law of thermodynamics. This law is a balance between the rates of change of kinetic and internal energies, the power of external forces, and the amount of heat entering the system as volume or surface heat sources.

d dt Ω 1 2 ρv • v + ρe dv = ∂Ω (σ • n) • v ds + Ω ρb • v dv + Ω ρr dv - ∂Ω q • n ds ∀ t ∈ τ, ∀ Ω
where q is the outward heat flux vector, r is a volume heat source and e is the internal energy density. The divergence theorem (2.24) yields:

d dt Ω 1 2 ρv • v + ρe dv = Ω (∇ • (σ • v) + ρb • v) dv + Ω ρr dv - ∂Ω q • n ds ∀ t ∈ τ, ∀ Ω
The transport of this relation in the reference configuration and introduction of the Lagrangian linear momentum (2.26) and of kinetic conservation laws (2.20) lead to:

Ω0 ρ 0 ėdV = Ω0 Π : Ḟ dV + Ω0 (ρ 0 r -∇ 0 • Q) dV ∀ t ∈ τ
where Q = Jq • F -1 is the Lagrangian heat flux vector. One thus deduces the balance equation of internal energy in the reference configuration:

ρ 0 ė -Π : Ḟ + ∇ 0 • Q = ρ 0 r ∀ X, t ∈ Ω 0 × τ (2.28)
Finally, the small strain version of equation (2.28) is:

ρ ė -σ : ε + ∇ • q = ρr ∀ x, t ∈ Ω × τ (2.29)
Strain and stress are then conjugate fields through an energy function. The former are referred to as state variables describing the evolution of the thermodynamic system while the latter are thermodynamic forces governed by constitutive equations. In what follows, such constitutive equations are derived.

Constitutive equations -Thermodynamics

The closure of the continuum equations is given by constitutive equations for the stress. Once and for all, we consider here constitutive models within the Generalized Standard Materials (GSM) framework [START_REF] Halphen | Sur les matériaux standard généralisés[END_REF].

The general hyperelasticity formulation

First, the Clausius-Duhem inequality resulting from combination of first and second laws of thermodynamics, reads:

Π : Ḟ + ρ 0 (θ η -ė) D int - 1 θ q • ∇ 0 θ D th ≥ 0 ∀ X, t ∈ Ω 0 × τ (2.30)
where θ and η denote the temperature and the entropy, and D int and D th are respectively the mechanical and thermal dissipations. The relation (2.30) becomes an equality for reversible processes and a strict inequality for irreversible ones. Furthermore, a widely used assumption consists in considering that mechanical and thermal dissipations simultaneously satisfy non-negativeness. Note that the Fourier's law of conduction is based on the non-negativeness of the thermal dissipation and leads to the following definition of the heat flux vector in order to ensure the positiveness of the thermal dissipation:

q = -k • ∇ 0 θ
where k is a positive-definite second-order tensor.

We assume that the internal energy density is a function of strain, entropy and additional internal variables

V p (1 ≤ p ≤ N )
, describing irreversible processes. The Helmholtz free energy density on the other hand, defined as the Legendre transform of internal energy, is a function of temperature and not of entropy:

ψ (F , θ, V) = e (F , η, V) -θη.
The free energy density is supposed objective or frame indifferent [26, p.255], concave with respect to temperature and convex with respect to other variables. The mechanical dissipation thus reads:

D int = Π : Ḟ -ρ 0 ψ + η θ
or, by introducing the time derivative of the Helmholtz free energy density ψ = ∂ψ ∂F :

Ḟ + ∂ψ ∂θ θ + ∂ψ ∂V V D int = Π -ρ 0 ∂ψ ∂F : Ḟ -ρ 0 ∂ψ ∂θ + η θ -ρ 0 ∂ψ ∂V V
Since the mechanical dissipation must be non-negative regardless of the nature of the transformation, it must in particular vanish for a reversible isothermal process (i.e. θ = const) for which every additional internal variable is constant (i.e. V = 0). With these considerations, we are left with the relation:

Π -ρ 0 ∂ψ ∂F : Ḟ = 0 2.2. GOVERNING EQUATIONS OF SOLID MECHANICS
holding regardless of the deformation, and hence:

ρ 0 ∂ψ ∂F = Π (2.31)
A material is said hyperelastic if there exists a stored energy density function ρ 0 ψ from which can be derived the first Piola-Kirchhoff stress tensor [24, p.8]. Similar considerations lead to the state laws for entropy and are assumed for additional thermodynamic forces associated with internal variables V:

∂ψ ∂θ = -η ; ρ 0 ∂ψ ∂V = -A (2.32)
Remark 2 Temperature has been introduced as a state variable and requires the first principle of thermodynamics, rewritten as the heat equation, in order to close the system:

ρ 0 C θ = ρ 0 r -∇ 0 • Q -ρ 0 ∂ψ ∂V V + θ ∂Π ∂θ : Ḟ - ∂A ∂θ V
Nevertheless, we will restrict our attention in the following to isothermal deformations so that temperature can be omitted and internal energy balance equations (2.28) or (2.29) are not considered.

For isothermal reversible deformations in hyperelastic solids, the time derivative of equation (2.31) leads to:

Π = ρ 0 ∂ 2 ψ ∂F 2 : Ḟ = H :

Ḟ

where H is the fourth-order tangent modulus tensor (major symmetric). The above discussion is now specified to constitutive models that will be used in the remainder of the manuscript.

Example 3 (Nearly incompressible Neo-Hookean)

The nearly incompressible neo-Hookean hyperelastic model is well-suited to describe rubber-like materials and is based on the polyconvex stored energy function (i.e. convex with respect to all its arguments):

ρ 0 ψ(J, F ) = κ 2 (J -1) 2 + µ 2 J -2/3 (F : F ) -3
where κ is the bulk modulus and µ the Lamé shear modulus. The first Piola-Kirchhoff stress and the acoustic tensor A ij = N α H iαjβ N β are for this model [START_REF] Haider | A first-order hyperbolic framework for large strain computational solid dynamics: An upwind cell centred total lagrangian scheme[END_REF]:

Π = µJ -2/3 F - 1 3 (F : F ) F -T + κ (J -1) H and A = 5 9 µJ -8/3 (F : F ) + κ (H • N ) ⊗ (H • N ) + µJ -2/3 I - 2 3 µJ -5/3 [(H • N ) ⊗ (F • N ) + (F • N ) ⊗ (H • N )] (2.33)
where H = JF -T is the adjoint tensor of the deformation gradient. The polyconvexity of this model ensures the positive definiteness of the acoustic tensor (2.33) [START_REF] Kluth | Analyse mathématique et numérique de systèmes hyperélastiques et introduction de la plasticité[END_REF].

Example 4 (Saint-Venant-Kirchhoff) The Saint-Venant-Kirchhoff hyperelastic model is based on the stored energy function:

ρ 0 ψ = 1 8 F T F -I : C : F T F -I (2.34)
where C is the fourth-order elasticity tensor defined as:

C iαjβ = λδ iα δ jβ + µ (δ ij δ αβ + δ iβ δ jα )
, with Lamé parameters (λ, µ). By differentiating the stored energy function (2.34) with respect to the deformation gradient, the PK1 stress is:

Π = 1 2 λ (F : F -3) F + µF F T F -I
whose derivative with respect to F yields the tangent modulus:

H = λ F ⊗ F + 1 2 (F : F -3)I + µ (B -I) (2.35) with B iαjβ = F kα F kβ δ ij + F jα F iβ + F iµ F jµ δ αβ
and the fourth order identity tensor I iαjβ = δ ij δ αβ . The previous tangent modulus finally leads to the acoustic tensor:

A = λ (F • N ) ⊗ (F • N ) + 1 2 (F : F -3)I +µ (F • N ) • (F • N )I + (F • N ) ⊗ (F • N ) + F • F T -I
Even though this model can lead to non-physical solutions, as we shall see in section 2.4.4, it will be used for a one-dimensional strain problem for it enables the development of an exact solution.

The infinitesimal theory formulation

The linearized geometrical framework leads, by assuming the existence of a Helmholtz free energy density ψ that depends on the temperature θ, the infinitesimal strain tensor ε and additional internal variables V, to the following relation [START_REF] Simo | Computational inelasticity[END_REF]Ch.2]:

ρ ∂ψ ∂ε = σ (2.36)
The infinitesimal strain tensor is further assumed to be additively decomposed into an elastic and a plastic part:

ε = ε e + ε p .
Then, with irreversible deformations due to plastic strains, the mechanical dissipation reads:

D int = σ : εp -ρ ∂ψ ∂V V ≥ 0
A yield condition is defined by means of function f (σ, A) so that the elastic domain E in forces space (σ, A) corresponds to:

E = {(σ, A) | f (σ, A) ≤ 0} (2.37)
According to the GSM framework [START_REF] Halphen | Sur les matériaux standard généralisés[END_REF] we assume the existence of a dissipation pseudo-potential Φ(σ, V), convex with respect to thermodynamic forces and vanishing at the origin of the (σ, V) space. This pseudo-potential enables the derivation of the plastic flow and hardening rules:

εp = ∂Φ ∂f ∂f ∂σ (2.38a) V = - ∂Φ ∂f ∂f ∂A (2.38b)
where ∂Φ ∂f = ṗ is the equivalent plastic strain rate. An example of a model used for metals is the plastic J 2 flow theory where the elastic domain is here described by a yield function that depends on the deviatoric part of the Cauchy stress, s = σ -1 3 trace σI, through its second invariant J 2 (s) = 1 2 s : s. In addition, a set of internal variables and associated forces describing the plastic hardening of the material is used

{V, A} = {[ε p , p] , [s -Y , -R(r)]
} in order to define the von-Mises yield surface:

f (σ, A) = 3 2 s -Y -(R(r) + σ y ) ≡ 0 (2.39)
where σ y is the tensile yield stress and r is to be defined. In deviatoric stress space, the von-Mises yield surface is a circle whose center and radius are Y and R(r). These thermodynamic forces hence respectively describe the displacement of the elastic domain center due to kinematic hardening, and the evolution of its radius due to isotropic hardening. Setting R(r) (resp. Y ) to zero amounts to specializing the yield surface (2.39) to kinematic (resp. isotropic) hardening. Then, flow rules (2.38a) and (2.38b) applied to the yield function (2.39) lead to:

εp = ṗ 3 2 s -Y s -Y = ṗ 3 2 m (2.40a) ṙ = ṗ (2.40b)
where r is identified from (2.40b) as the equivalent plastic strain and m is referred to as the plastic flow direction.

Next, assuming a Prager-Ziegler linear kinematic hardening [26, p.91], the Helmholtz free energy density takes the form:

ρψ = 1 2 ε e : C : ε e + 2 3 Cε p : ε p + H(p) (2.41)
where H(p) is defined so that H (p) = C is the hardening modulus, and C is the fourth-order elastic stiffness tensor (major and minor symmetric) defined for isotropic materials as

C ijkl = λδ ij δ kl + µ (δ ik δ jl + δ il δ jk ).
Thermodynamic forces are finally related to internal variables by means of equation (2.32), that is:

Y = 2 3 Cε p R = H (p)
We consider in what follows isothermal deformations of isotropic solids, that may be irreversible by specifying the above developments to some well-known small strain constitutive models.

Example 5 (Linear elasticity)

The simplest case that is considered hereinafter does not involve irreversible deformations, and hence additional internal variables (i.e. ε p ≡ 0), and is referred to as linear elasticity. The combination of equations (2.36) and (2.41) then leads to Hooke's law:

σ = C : ε or in rate form: σ = C : ε (2.43)
The elastic acoustic tensor is further defined as:

A elast ij = n k C kijl n l = λn i n j + µ (n i n j + δ ij )
Example 6 (Elastoplasticity) Rate-independent plasticity or elastoplasticity is based on the assumption that admissible thermodynamic forces lie within or on the boundary of the elastic domain (2.37). The equivalent plastic strain rate becomes a Lagrange multiplier in order to ensure f (σ, A) ≤ 0 and must obey the Kühn-Tucker compatibility conditions:

ṗ ≥ 0 ; f ≤ 0 ; ṗf = 0
The equivalent plastic strain rate, is determined by the consistency condition ḟ = ∂f ∂(s-Y ) : ( ṡ -Ẏ ) -∂f ∂R Ṙ = 0 that leads to:

3 2 m : Ẏ + Ṙ = 3 2 m : σ
Then, combination of the above equation with the elastic law σ = C : ( εεp ) and equation (2.40a) yields:

ṗ = 3 2 2µ 3µ + (C + R ) m : ε = 3 2 2µ 3µ + (C + R ) m : ε (2.44)
At last, equations (2.40a) and (2.44) can be successively introduced in the elastic law so that one gets [26, eq (2.2.22)]:

σ = C - 6µ 2 3µ + (C + R ) m ⊗ m : ε = C ep : ε (2.45)
with C ep the elastoplastic tangent modulus. The elastoplastic acoustic tensor is defined as:

A ep ij = n k C ep ikjl n l = A elast ij - 6µ 2 3µ + (C + R ) (n k m ik )(m jl n l )
which is positive-definite for positive linear hardening (C > 0 ; R > 0).

Example 7 (Elasto-viscoplasticity)

Viscoplasticity or rate-dependent plasticity can be seen as a regularization of rate-independent plasticity that relaxes the condition f (σ, A) ≤ 0 and thus leads to admissible thermodynamic forces lying outside the elastic domain [26, p.58]. Viscoplasticity provides on the other hand an explicit definition of the equivalent plastic strain, for example the Perzyna or Sokolowskii-Malvern model [START_REF] Perzyna | Fundamental problems in viscoplasticity[END_REF] is governed by:

ṗ = f γ n (2.46)
where

• = •+|•| 2
is the positive part function, and n and γ are parameters. Hence, the plastic fluxes εp , V are completely determined by (2.38a) and (2.38b). It then comes out that rate-dependent plasticity is driven by the elastic law:

σ = C : ( ε -εp )
in which εp is given by the combination of equations (2.38a) and (2.46), namely:

εp = f γ n 3 2 m (2.47)

The general formulation

Balance and constitutive equations obtained previously are now summarized for various classes of materials and regimes of deformation. Recall that the deformations are assumed isothermal and that history effects are considered within the infinitesimal theory only.

Hyperelasticity

The system of conservation laws for problems involving hyperelastic solids is composed of kinematic laws (2.20) and the balance equation of Lagrangian linear momentum (2.26), repeated here for convenience:

Ḟ -∇ 0 • (v ⊗ I) = 0 ρ 0 v -∇ 0 • Π = ρ 0 b
Assuming a Cartesian coordinate system, this system can be written in conservative form:

U t + D α=1 ∂F • E α ∂X α = S (2.48)
where the vector of conserved quantities U, flux vectors F • E α and the source term S are:

U = ρ 0 v F ; F • E α = -Π • E α -v ⊗ E α ; S = ρ 0 b 0 (2.49)
A quasi-linear system may then be built by introducing an auxiliary vector Q = v Π and using the chain rule according to [START_REF] Trangenstein | A higher-order godunov method for modeling finite deformation in elasticplastic solids[END_REF]:

∂Q ∂t + ∂U ∂Q -1 ∂F • E α ∂Q ∂Q ∂X α = ∂U ∂Q -1

S = S

In the quasi-linear form, the derivative of the vector of conserved quantities with respect to the auxiliary vector leads to the diagonal matrices:

∂U ∂Q = I 0 3 0 3 ∂F ∂Π ⇒ ∂U ∂Q -1 = I 0 3 0 3 ∂Π ∂F
where the tangent modulus ∂Π ∂F = H arises and 0 p is a pth-order zero tensor. Moreover, the derivative of the flux vectors with respect to the auxiliary vector reads:

∂F • E α ∂Q = - 0 2 1 ρ0 I ⊗ E α I E α 0 4 21

GOVERNING EQUATIONS OF SOLID MECHANICS

in which the operator I E α is the transpose on second and third indices of the classical tensor product, namely: I E α = δ jk e j ⊗ E α ⊗ e k . Finally, the quasi-linear form associated with hyperelastic problems is:

Q t + A α ∂Q ∂X α = S with: A α = - 0 2 1 ρ0 I ⊗ E α H • E α 0 4 (2.50)
where the dependence on Q of matrices A α (Q) has been omitted for simplicity.

Linear elasticity and elasto-viscoplasticity

The governing equations of elasticity and elasto-viscoplasticity within the linearized geometrical framework consist of the kinematic law (2.21), the balance equation of linear momentum (2.27) and the elastic law (2.43):

ε -∇ • v ⊗ I + I ⊗ v 2 = 0 ρ v -∇ • σ = ρb σ -C : ( ε -εp ) = 0
Combining kinematic and elastic laws and considering again a Cartesian coordinates system yields, for a homogeneous media (i.e. ∇ρ = 0), the following conservative form:

Q t + D i=1 ∂F • e i ∂x i = S (2.51)
with conserved quantities, flux and source term vectors respectively defined as:

Q = v σ ; F • e i = -1 ρ σ • e i -C : v⊗ei+ei⊗v 2 ; S = b -C : εp
Note that here, the direct writing of the conservative form in terms of v and σ is made possible by the linearity of the elasticity tensor, avoiding thus the introduction of an auxiliary vector. The quasi-linear form of equation (2.51) is derived by means of the chain rule:

Q t + A i ∂Q ∂x i = S (2.52)
where:

A i = ∂F • e i ∂Q = - 0 3 1 ρ I ⊗ e i C • e i 0 4
in which symmetries of the elastic stiffness tensor have been used. Since the elastic stiffness tensor is constant, system (2.52) is linear for elasticity and semi or non-linear for elasto-viscoplasticity depending on the flow rule (2.38a). Moreover, the source term arising due to the viscoplastic flow rule (2.47) can be written in terms of a relaxation term S and a relaxation time τ = (γ/σ y ) n as S = S/τ [START_REF] Heuzé | Simulation of impacts on elastic-viscoplastic solids with the flux-difference splitting finite volume method applied to non-uniform quadrilateral meshes[END_REF], so that the system (2.52) can be identified to a relaxation system [START_REF] Liu | Hyperbolic conservation laws with relaxation[END_REF]. In the asymptotic limit τ → 0 or in the vanishing viscosity limit, system (2.52) tends to the equilibrium system corresponding to elastoplasticity [START_REF] Heuzé | Simulation of impacts on elastic-viscoplastic solids with the flux-difference splitting finite volume method applied to non-uniform quadrilateral meshes[END_REF].

Elastoplasticity

The writing of a conservative form for elastoplasticity is similar to what was done for hyperelastic solids. Indeed, the system composed of kinematic laws (2.21) and the balance equation of linear momentum (2.27):

ε -∇ • v ⊗ I + I ⊗ v 2 = 0 ρ v -∇ • σ = ρb 22 CHAPTER 2.
HYPERBOLIC PARTIAL DIFFERENTIAL EQUATIONS FOR SOLID DYNAMICS can be written as:

U t + D i=1 ∂F • e i ∂x i = S (2.53)
where the conserved quantities, flux and source term vectors are:

U = v ε ; F • e i = -1 ρ σ • e i -v⊗ei+ei⊗v 2 ; S = b 0 (2.54)
Analogously to hyperelasticity, a quasi-linear form involving the elastoplastic tangent modulus (5.2c) is derived by means of the auxiliary vector Q = v σ and the chain rule:

Q t + A i ∂Q ∂x i = S with: A i = - 0 2 1 ρ I ⊗ e i C ep • e i 0 4
(2.55)

Characteristic analysis -Structure of solutions

The quasi-linear forms written above enable the particularization of the theory of first order quasi-linear systems developed in section 2.1 to solid mechanics problems.

For the sake of simplicity, studies of finite deformations and linearized geometrical frameworks will be condensed in this part by using the generic notation of stress S and vectors written in the reference configuration. Furthermore, instead of considering multi-dimensional systems of conservation laws, we will focus without loss of generality on the quasi-linear form (2.50) projected on an arbitrary direction

N = [E 1 , E 2 , E 3 ] [2, p.425-426].
In this direction, one has:

Q t + J ∂Q ∂X N = S (2.56)
in which X N = X • N and the Jacobian matrix J = A α N α of dimension m arise. In three dimensions, the non-symmetrical PK1 tensor in quasi-linear form (2.50) yields a Jacobian matrix of dimension m = 3 + 9 while systems (2.52) and (2.55) involving the Cauchy tensor, lead to m = 3 + 6. The characteristic analysis of system (2.56) is therefore equivalent to that of linear combinations of matrices A α . With the previous developments, the Jacobian matrix reads:

J = - 0 2 1 ρ0 I ⊗ N H • N 0 4
in which H is either the hyperelastic or elastoplastic tangent modulus, or the elastic stiffness tensor depending on the case considered. The characteristic structure of the problem is given by the m eigenvalues c K and associated left eigenvectors L K = v K , S K of the Jacobian matrix satisfying:

L K (J -c K I) = 0
where I is the m × m identity matrix. Thus, for non-zero eigenvalues one gets:

-S K : H • N -c K v K = 0 (2.57a) - 1 ρ 0 v K ⊗ N -c K S K = 0 (2.57b)
Substitution of S K obtained from (2.57b) in (2.57a) leads to:

(v K ⊗ N ) : H • N -ρ 0 c 2 K v K = 0 (2.58)
which is the left eigensystem of the acoustic tensor A ij = N α Hiαjβ N β . Due to the symmetry of A, system (2.58) is equivalent to the right eigensystem:

N α Hiαjβ N β -ρ 0 c 2 K δ ij v K j = 0
or alternatively, with the eigenvalues ω p and associated eigenvectors of the acoustic tensor l p (p = 1, 2, 3):

(A -ω p I) • l p = 0
The condition for system (2.56) to be hyperbolic (real eigenvalues and independent eigenvectors) is thus ensured by the positive definiteness of the acoustic tensor, also known as the strong ellipticity condition [START_REF] Marsden | Mathematical foundations of elasticity[END_REF]:

(n ⊗ N ) : H : (n ⊗ N ) > 0 ∀N , n ∈ R 3 ; N , n = 0
If the condition holds, the acoustic tensor admits 3 couples of eigenvalue-eigenvector {ω p , l p } leading to 6 couples {c K , L K } for the Jacobian matrix, the 6 other eigenvalues being null [START_REF] Kluth | Analyse mathématique et numérique de systèmes hyperélastiques et introduction de la plasticité[END_REF]. The couples {c K , L K } are referred to as left characteristic fields. The left eigenvectors associated with non-zero eigenvalues of the Jacobian matrix are obtained by using equation (2.57b) so that the following 6 eigenfields of the quasi-linear form (2.56) can be defined:

± ω p ρ 0 ; ±ρ 0 ω p ρ 0 l p , -l p ⊗ N , p = 1, 2, 3
At last, one has to find six independent left eigenvectors associated with the null eigenvalue of multiplicity 6 by solving equation (2.57a) for the null eigenvalue:

S K : H • N = 0, K = 1, ..., 6
Following the same procedure for right eigenvectors R K = v K S K , the Jacobian matrix right eigensystem reads:

-1

ρ 0 S K • N -c K v K = 0 -H : v K ⊗ N -c K S K = 0
which leads to the right characteristic fields associated with the non-null eigenvalues:

± ω p ρ 0 ; ± ω p ρ 0 l p , -H : (l p ⊗ N ) , p = 1, 2, 3 (2.60)
In equation (2.60), {ω p , l p } still denotes the eigenfields of the acoustic tensor. Moreover, the 6 independent right eigenvectors associated with the zero eigenvalue required to complete the set of right characteristic fields must satisfy:

S K • N = 0, K = 1, ..., 6
Remark 3 Since the right-hand side of equation (2.56) is not involved in the characteristic analysis, linear elasticity and elasto-viscoplasticity in small strains share the same characteristic structure.

Remark 4

In the case of a vanishing source term S, the specialization of characteristic equations (2.12) to system (2.56) leads to:

L K • dQ = 0, K = 1, ..., 6
meaning that the solution is constant along each characteristic straight line with slope ξ = c K . Such solutions Q(ξ) that only depend on the ray ξ are called self-similar solutions.

Some solutions of Riemann problems

The characteristic analysis carried out above is now applied to specific solid mechanics problems. Both linear and non-linear problems, whose solutions involve several types of waves, are considered. As we shall see, the different characteristic structures involved within linear elastic, elastoplastic and hyperelastic solids require the use of different techniques in order to develop exact solutions. First a particular type of IVP, of particular interest in this manuscript, is introduced.

The Riemann problem

A Riemann problem is a Cauchy problem with piecewise constant initial data. In particular, the Riemann problem based on the conservative form (2.48) for hyperelastic solids, in the arbitrary direction N = N α E α , takes the form:

U t + ∂F • N ∂X N = S, U(X N , t = 0) = U L if X N < 0 U(X N , t = 0) = U R if X N > 0 (2.61)
Analogously, for small strains one writes the Riemann problem corresponding to conservative forms (2.51) or (2.53) in the direction n = n i e i :

Q t + ∂F • n ∂x n = S, Q(x n , t = 0) = Q L if x n < 0 Q(x n , t = 0) = Q R if x n > 0 (2.

Linear elastodynamics problems

A homogeneous hyperbolic system of dimension m is considered in a linear elastic solid so that a Riemann problem of the form (2.62) is written.

Characteristic variables -Waves solution

By introducing a set of characteristic variables

P = R -1 Q (R ij = R j i )
, the quasi-linear form of system (2.62) reads:

∂P i ∂t + c i ∂P i ∂x = 0 P i (x, t = 0) = P L i if x < 0 P i (x, t = 0) = P R i if x > 0 (2.63) with C ij = c i δ ij , the matrix of eigenvalues so that J ij R j K = R K i C Kj .
The solution of this problem is straightforward since it corresponds to a superposition of scalar linear advection equations, namely, the initial profile P i (x, t = 0) simply propagates with speed c i as depicted in figure 2.6. Thus, the solution P i (x, t) at a given point is given by tracing backward the characteristic of slope c i passing through this point to the x-axis, that 2.4. SOME SOLUTIONS OF RIEMANN PROBLEMS is: P i (x, t) = P i (xc i t, 0) [21, p.52]. The vector Q is then determined by inverting the relation: [START_REF] Leveque | A study of numerical methods for hyperbolic conservation laws with stiff source terms[END_REF]) is an eigenvector expansion with coefficients P R,L i from which we see that Q is a linear superposition of m waves, each having the shape R i P i (x, 0). Noticing that for given values of x and t, there exists one characteristic I such that xc i t > 0 for all i ≤ I, and xc i t < 0 for all i ≥ I + 1, equation (2.64) can be rewritten [21, p.56]:

x t c i t 1 t 2 t 0 Pi(x, t2) x c i t 2 Pi(x, t1) x c i t 1 Pi(x, t0) x P L i P R i
Q(x, t) = m i=1 R i P i (x -c i t, 0) ⇒            Q(x < 0, 0) = Q L = m i=1 R i P L i Q(x > 0, 0) = Q R = m i=1 R i P R i (2.64) Equation (2.
Q = I i=1 R i P R i + m i=I+1 R i P L i (2.65)
or, by introducing the expansions of initial data (2.65):

Q = m i=1 R i P R i - m i=I+1 R i P R i -P L i = Q R - m i=I+1 R i P R i -P L i Q = m i=1 R i P L i + I i=1 R i P R i -P L i = Q L + I i=1 R i P R i -P L i
These equations are equivalent to jump conditions across multiple discontinuous waves:

Q -Q R = - m i=I+1 R i δ i (2.66) Q -Q L = I i=1 R i δ i (2.67)
where Q(x, t) is the state lying in the region of the (x, t) plane delimited by the characteristics I and I + 1, and R i δ i the jump carried by the ith wave. The coefficients δ i = P R i -P L i are weighting coefficients involved in wave strengths W i = R i δ i , that can be computed from the expansions of initial conditions by solving:

Q R -Q L = m i=1 R i δ i = Rδ (2.68)
We see that the solutions of Riemann problem (2.62) consist of discontinuous waves emanating from the origin of the (x, t) plane. Across such discontinuous waves, the following condition is satisfied [START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics: a practical introduction[END_REF]:

Definition 3
The Rankine-Hugoniot condition is satisfied across a discontinuous wave of speed s i associated with the ith characteristic field, which is a solution of the hyperbolic system Q t + F(Q) x = 0:

F = s i Q (2.69)
where • denotes the jump operator across the discontinuity. More specifically, shock waves that will be met for non-linear problems in section 2.4.4 also satisfy the Rankine-Hugoniot condition.

Solution of the elastic bar problem

The above discussion is now specified to a one-dimensional elastic medium, x ∈ [-l, l], of density ρ undergoing one-dimensional stress and strain states within the infinitesimal framework: ε = ε e 1 ⊗ e 1 ; σ = σ e 1 ⊗ e 1 . As a consequence, the bar hypothesis holds with v = ve 1 . Neglecting body forces and introducing Young's modulus E such that σ = Eε, the Riemann problem takes the form (2.62) with conserved quantities and flux vector:

Q = v σ ; F = -1 ρ σ -Ev
along with Riemann-type data on the horizontal velocity (i.e. v(x < 0, 0) = v L ; v(x > 0, 0) = v R ) as initial conditions. In addition, the solid is assumed to be initially unstressed. The eigenvalues, left and right eigenvectors of the corresponding Jacobian matrix are:

c = ± E ρ ; L p = [ρc p , -1] ; R p = 1 -ρc p
The characteristic structure of the solution consisting of two elastic discontinuities emanating from the origin of the (x, t) plane is depicted in figure 2.7. The solution Q * lying in the region bounded by the two elastic waves is computed by means of equation (2.67) after solving (2.68) for the wave strength coefficients. For a system of dimension 2 by writing Q R -Q L = ∆Q those wave strengths read:

x t c -c Q R Q L Q *
δ = 1 R 1 1 R 2 2 -R 2 1 R 1 2 R 2 2 ∆Q 1 -R 2 1 ∆Q 2 R 1 1 ∆Q 2 -R 1 2 ∆Q 1 (2.70)
and more specifically for a bar:

δ = 1 2ρc
ρc∆v + ∆σ ρc∆v -∆σ Hence, equation (2.67) yields the solution:

Q * = Q L + R 1 δ 1 = σ R -σ L 2ρc + v R +v L 2 ρc v R -v L 2 + σ R +σ L 2

Elastic-plastic media in the geometrical linearized limit

The previous problem is now extended to elastoplastic media by considering the same bar made of a linear hardening material of tensile yield stress σ y . For such a solid domain, a Riemann problem of the form (2.62) can be written by means of the following conserved quantities and flux vectors [START_REF] Heuzé | Lax-Wendroff and TVD finite volume methods for unidimensional thermomechanical numerical simulations of impacts on elastic-plastic solids[END_REF]:

Q = v σ ; F = -1 ρ σ -Hv
where H = E for elastic loadings while H = dσ/dε is the tangent modulus for elastic-plastic evolutions. In addition, Riemann-type data on the horizontal velocity (i.e. v(x < 0, 0) = v L ; v(x > 0, 0) = v R ) are used as initial conditions, so that plastic flow may occur.

The discontinuity of H across the plastic threshold prevents here the direct derivation of the solution by using the approach followed for linear elasticity. Indeed, two sets of characteristic speeds and associated eigenvectors must be considered, that is [START_REF] Heuzé | Lax-Wendroff and TVD finite volume methods for unidimensional thermomechanical numerical simulations of impacts on elastic-plastic solids[END_REF][START_REF] Wang | Foundations of stress waves[END_REF]:

c = ± H ρ ; L p = [ρc p , -1] ; R p = 1
-ρc p so that waves are referred to as elastic or plastic waves. Whether a plastic wave appears or not hence depends on the tangent modulus and subsequently, on the yield function. A predictor-corrector procedure must thus be followed by first solving an elastic Riemann problem whose resulting (trial) solution Q is tested against the yield criterion on both sides x < 0 or x > 0. In general, possibly different yield stresses, plastic strains and hardening parameters in left and right regions lead to a yield criterion that may be violated or not, and hence to one, two or no plastic waves that must be added as a correction to the original problem. If neither f L nor f R leads to a violation of the criterion, the problem is elastic and the trial state is the solution. Otherwise, the plastic correction is performed by computing the stress in regions of the (x, t) plane bounded by elastic and plastic waves ( QL,R in figure 2.8a) so that the yield function satisfies f L,R = 0. Then, the velocity and elastic (a) Characteristic structure wave strengths in yielding regions are given by solving:

x t c c p -c -c p t 1 t 2 Q L Q R QL QR Q * (b) Stress profile in the bar σ(x, t 1 ) x σ y σ(x, t 2 ) x σ y
QR = Q R -δ 1 E R 1 E QL = Q L + δ 2 E R 2 E
At last a plastic Riemann solver is used in order to compute the solution of the problem Q * by solving successively the system (2.68) for plastic wave strengths, and either system (2.66) or (2.67), that is:

δ P = R -1 QR - QL ⇒ Q * =    QR -δ 2 P R 2 P QL + δ 1 P R 1 P Figure 2
.8a shows the characteristic structure of the solution of the Riemann problem in a homogeneous medium considered here, involving two plastic waves, and figure 2.8b the corresponding stress field in the bar.

Remark 5 Note that the above solver also applies to plane wave problems characterized by one-dimensional strain and multi-dimensional stress states by considering different wave speeds.

Hyperelastic media: A Saint-Venant-Kirchhoff solution

The approaches followed above are no longer possible for problems involving a non-linear Jacobian matrix. Indeed, the writing of the Riemann problem (2.63) in terms of characteristic variables is valid for right eigenvectors that are constant. Moreover, as we shall see with an example, the characteristic structure of such problems can be more complex and depend on the initial data. Consider a hyperelastic medium made of a Saint-Venant-Kirchhoff material, infinite in directions E 2 and E 3 , and semi-infinite in direction E 1 (i.e. X 1 ∈ [0, +∞[) in the reference configuration. This medium suddenly undergoes a load at (X 1 = X = 0, t = 0) in direction E 1 so that the deformation gradient and the PK1 tensors are respectively of the form:

F = F e 1 ⊗ E 1 + e 2 ⊗ E 2 + e 3 ⊗ E 3 (2.71) Π = Π 11 e 1 ⊗ E 1 + Π 22 (e 2 ⊗ E 2 + e 3 ⊗ E 3 )
which corresponds to a plane wave solution. We assume that F (0, t) = F is given, leading to a Picard problem [36, p.20] involving both initial and boundary conditions with neglected body forces:

Q t + ∂F • N ∂X N = 0, Q(X N , t = 0) = Q R if X N > 0 F (0, t) = F (2.72)
with N = E 1 and:

Q = v F ; F = -1 ρ0 Π -v
where Π = Π 11 . The quasi-linear form is written by using the chain rule: Q t + ∂F ∂Q ∂Q ∂X = 0 so that the Jacobian matrix reads:

J = ∂F ∂Q = - 0 H1111 ρ0 1 0
The tangent modulus of the SVK model (2.35) yields the following characteristic fields:

           c 1 = - λ + 2µ 2ρ 0 (3F 2 -1) c 2 = λ + 2µ 2ρ 0 (3F 2 -1) ; L p = [1 , -c p ] ; R p = -c p 1 (2.73)

Remark 6

The non-linear flux function of the SVK model yields characteristic fields depending on the strain state and possibly complex celerities leading to a loss of hyperbolicity of the problem for F < 1 3 .

Suppose now that initial data are given so that F > F R . The resulting characteristic speeds then satisfy c 2 ( F ) > c 2 (F R ) and the two families of characteristics collide in the right region of the (x, t) plane (figure 2.9a).

On the other hand, F < F R yields characteristics moving away from each other in the right region according to

c 2 ( F ) < c 2 (F R ) (figure 2.9b
). These two situations respectively correspond to a shock and a simple wave. Note that this characteristic structure is similar to that resulting from the dam-break problem with shallow water equations, the following developments are hence very close to those of [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF]Ch.13].

(a) Right-going shock wave

Shock waves

A shock wave is a discontinuous wave satisfying the Rankine-Hugoniot condition (2.69):

-1

ρ 0 Π -Π i = s (v -v i ) (2.74) -(v -v i ) = s F -F i (2.75)
where the shock speed s is to be defined and i ∈ {L, R}. In what follows, F is considered as an unknown so that a relation connecting Q i to a set of solutions Q through a shock wave can be developed. Substitution of s from equation (2.75) and introduction in equation (2.74) where Π = λ+2µ 2 F 3 -F yield:

s = - v -v i F -F i v -v i = ± λ + 2µ 2ρ 0 ( F -F i ) F 3 -F -(F 3 i -F i ) (2.76)
In addition to the Rankine-Hugoniot condition, one has to consider the Lax entropy conditions, stating that characteristic curves collide in a shock wave [2, p.268]:

c( F ) < s < c(F i ) (2.77)
The Lax entropy condition implies that the square root in equation (2.76) is real, leading to two families of curves in the phase plane (F, v). When considering an infinitesimal jump (i.e. F = F i ± with → 0), each of these curves is expected to identify with one of the jump conditions derived for the linear case (2.66) of (2.67). Thus, equation (2.76) reads:

v -v i = ± λ + 2µ 2ρ 0 [(F i + ) 3 -(F i + ) -(F 3 i -F i )]
where, with → 0:

(F i + ) 3 ≈ F 3 i (1 + 3 F i ) so that: v F = v i F i + ± λ+2µ 2ρ0 [3F 2 i -1] 1
The minus sign yields equation (2.66) associated with the right-going wave and therefore corresponds to a rightgoing shock wave. On the other hand, the plus sign stands for equation (2.67) and left-going shocks. Finally, the Rankine-Hugoniot condition across a left-going shock and a right-going shock respectively lead to:

v -v L = λ + 2µ 2ρ 0 ( F -F L ) F 3 -F -(F 3 L -F L ) (2.78) v -v R = - λ + 2µ 2ρ 0 ( F -F R ) F 3 -F -(F 3 R -F R ) (2.79)

Simple waves

In order to study the evolution of fields within the region bounded by characteristics that move away in figure 2.9b, let's write the left-going characteristic equation through it with

L 1 = [1, -c 1 ]: dv -c 1 (F )dF = 0 (2.80)
The complete set of states Q connected to Q R through a simple wave is obtained by integration of equation (2.80). Note that this integration results in a smooth evolution of fields inside a simple wave even for discontinuous initial conditions, unlike shocks. Moreover, the vanishing right-hand side of the conservative form of the Picard problem (2.72) yields a similarity solution. The particular case of a simple wave constant along each ray ξ = x/t corresponds to a rarefaction wave [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF]. Integration of equation (2.80) is performed by using the change of variables: F → ch(x)/ √ 3 so that one gets:

dv = - λ + 2µ 6ρ 0 sh(x) 2 dx (2.81)
where, the hyperbolic cosine ch(x) and sine sh(x) satisfy: ch(x) 2 sh(x) 2 = 1. Equation (2.81) can be easily integrated with the exponential form of hyperbolic sine:

sh(x) = e x -e -x 2
, thus yielding:

v -v R = - 1 4 λ + 2µ 6ρ 0 [sh(2x) -2x -(sh(2x R ) -2x R )]
At last, the inverse change of variables leads to the relation:

v -v R = - λ + 2µ 24ρ 0 √ 3 F 3F 2 -1 -F R 3F 2 R -1 -ln √ 3F + √ 3F 2 -1 √ 3F R + 3F 2 R -1 (2.82)
In a similar manner, dvc 2 (F )dF = 0 must hold through a left-going rarefaction wave so that:

v -v L = λ + 2µ 24ρ 0 √ 3 F 3F 2 -1 -F L 3F 2 L -1 -ln √ 3F + √ 3F 2 -1 √ 3F L + 3F 2 L -1 (2.83)
Equations (2.82) and (2.83) correspond to integral curves that connect initial conditions to a set of solutions through a right-going or a left-going rarefaction respectively.

Solution of the Riemann problem

The above developments are now generalized by considering the Riemann problem in an infinite medium:

Q t + ∂F • N ∂X N = 0,        Q(X N , t = 0) = Q L = v = 0 F L if X N < 0 Q(X N , t = 0) = Q R = v = 0 F R if X N > 0 (2.84)
such that the plane wave state (2.71) holds. As for the Picard problem, initial conditions influence the characteristic structure of the solution. Indeed, if initial conditions are given such that F L < F R , left-going characteristics will collide while right-going ones will move away from one another (see figure 2.10a). In that case, the first and second characteristic fields are respectively referred to as a 1-shock and a 2-rarefaction. Conversely, if F L > F R the solution corresponds to a 1-rarefaction and a 2-shock (figure 2.10b).

For the 1-shock-2-rarefaction solution one then seeks a state Q that is connected to Q L and Q R through a shock wave and a rarefaction wave respectively. Hence, Q must satisfy equations (2.78) and (2.82), that is: Analogously, the 1-rarefaction-2-shock solution is given by the solution of equations (2.83) and (2.79):

             v -v L = λ + 2µ 2ρ 0 (F -F L ) [F 3 -F -(F 3 L -F L )] v -v R = - λ + 2µ 24ρ 0 √ 3 F 3F 2 -1 -F R 3F 2 R -1 -ln √ 3F + √ 3F 2 -1 √ 3F R + 3F 2 R -1 (2.85) (a) F L < F R X t (b) F L > F R X t
             v -v L = λ + 2µ 24ρ 0 √ 3 F 3F 2 -1 -F L 3F 2 L -1 -ln √ 3F + √ 3F 2 -1 √ 3F L + 3F 2 L -1 v -v R = - λ + 2µ 2ρ 0 (F -F R ) [F 3 -F -(F 3 R -F R )] (2.86)
Once one of these systems is solved, the solution Q is known everywhere except inside the rarefaction fan. Nevertheless, in this region Q only varies with the ray ξ = c i (F ) and hence, the solution inside an i-rarefaction wave satisfies:

1 1.2 1.4 1.6 1.8 2 0 0.5 1 1.5 • • Q L Q R F v(m/s) (a) Initial data: F L = 1 ; F R = 2 Shock Rarefaction 1 1.2 1.4 1.6 1.8 2 -1.5 -1 -0.5 0 • • Q R Q L F v(m/s) (b) Initial data: F L = 2 ; F R = 1 Rarefaction Shock
ξ = ± λ + 2µ 2ρ 0 (3F 2 -1) ⇒ F (ξ) = 2ρ 0 3(λ + 2µ) ξ 2 -1
The curves corresponding to equations (2.85) and (2.86) are depicted in figures 2.11a and 2.11b for parameter values such that λ+2µ ρ0 = 1. In both cases, the solution Q(x, t) in the region bounded by the shock and the rarefaction fan is given by the intersection of curves in the phase plane.

Remark 7 Note that the above developments are based on a constitutive model that leads to a concave flux function F

T = -1 ρ0 Π ; v (i.e. ∂ 2 Π ∂F 2 > 0)
. As a consequence, the characteristic speeds are monotonically increasing functions of the deformation gradient (in absolute value). Since the dependence of characteristic speeds on the deformation gradient governs the wave pattern (i.e.: either a rarefaction or a shock wave), the structure of the solution differs from other constitutive laws with convex flux function such as the neo-Hookean model. Comparisons of (F, Π 11 ) and (F, |c|) are shown in figures 2.12a and 2.12b as an illustration of previous remarks. At last, figure 2.12a shows the non-physical behavior of Saint-Venant-Kirchhoff model for high-compression loads that lead to a stress tensor tending to zero.

0 0.5 1 1.5 2 -2 -1 0 1 •10 12

Approximate-State Riemann solvers

It has been seen in the previous section that the complete solution of a Riemann problem in solid dynamics is possible for simple problems. However, such a solution may become complicated for multi-dimensional problems or for other non-linear problems. Numerical methods such as upwind or Godunov-based methods [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF] require the solution of many Riemann problems within a discretized medium. When dealing with non-linear problems, the exact solution of those problems may increase drastically the computational cost, making the numerical scheme unappealing. Moreover, numerical procedures often require only little information about the solution of Riemann problems and do not need the complete solution. In that context, alternative procedures have been developed in order to take into account the characteristic structure of a hyperbolic system by computing an approximate solution of Riemann problems. Approximate Riemann solvers developed for Computational Fluid Dynamics allow to extract information for either flux functions (HLL, HLLC, Roe and Osher approximate Riemann solvers [START_REF] Trangenstein | Numerical solution of hyperbolic partial differential equations[END_REF], [START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics: a practical introduction[END_REF]) or for vectors of conserved quantities (approximate-state Riemann solver [21, Ch.9], [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF]Ch.22]). Some of these have been applied to specific problems in solid mechanics problems such as the Osher approximate solver (see [START_REF] Lee | Development of a cell centred upwind finite volume algorithm for a new conservation law formulation in structural dynamics[END_REF] and [START_REF] Haider | A first-order hyperbolic framework for large strain computational solid dynamics: An upwind cell centred total lagrangian scheme[END_REF]) or the HLLC approximate solver (see [START_REF] Ortega | Numerical simulation of elastic-plastic solid mechanics using an eulerian stretch tensor approach and hlld riemann solver[END_REF]) for hyperelasticity . We recall here the formulation of the approximate-state Riemann solver for solid mechanics. The approach is then applied to the non-linear problem of section 2.4.4.

General ideas

As in the previous section, we consider the Riemann problem in the space direction N :

Q t + J (Q) ∂Q ∂X N = 0, Q(X N , t = 0) = Q L if X N < 0 Q(X N , t = 0) = Q R if X N > 0 (2.87)

APPROXIMATE-STATE RIEMANN SOLVERS

The approach for developing an approximate-state Riemann solver consists in linearizing the problem (2.87) by approximating J in the vicinity of .15]. Note that this approximation is valid for small jumps in initial data (i.e Q L ≈ Q R ) and that J must ensure hyperbolicity of the system, namely J has real eigenvalues and a complete set of independent eigenvectors. The approximate matrix also satisfies the consistency condition:

Q L and Q R by a constant matrix J = J Q L , Q R [2, Ch
J (Q, Q) = J (Q)
Such a matrix can be defined by using the definition of right eigenvectors and characteristic speeds JR = RC ⇒ J = RCR -1 in which left-going (resp. right-going) characteristics and associated eigenvectors are assumed to depend on Q L (resp. on Q R ) only. Namely, one writes:

R = R 1 (Q L ), • • • , R I (Q L ), R I+1 (Q R ), • • • , R m (Q R ) C =         c 1 (Q L ) • • • c I (Q L ) c I+1 (Q R ) • • • c m (Q R )        
where c I (Q) and m are the highest negative eigenvalue and the dimension of the Jacobian matrix. At last, the linearized Riemann problem thus written enables the determination of every state vector Q(x, t) by following the procedure described in section 2.4.2 for linear problems, recalled here for convenience for a system of dimension m:

Q R -Q L = m i=1 R i δ i Q(x, t) = Q R - m i=I+1 R i δ i Q(x, t) = Q L + I i=1 R i δ i
where the point (x, t) lies in the region bounded by the characteristics I and I + 1.

Remark 8

Note that since one can define a complete set of independent eigenvectors of the Jacobian matrix, the matrix R is non-singular so that J can be uniquely determined.

Application: Hyperelastic plane wave

We finish this section with an illustration of the approximate Riemann solver by considering the plane wave problem in the Saint-Venant-Kirchhoff medium treated in section 2.4.4. Recall that the eigenvalues and right eigenvectors matrices read for that problem:

C = -c 0 0 c ; R = c -c 1 1 , c = λ + 2µ 2ρ 0 (3F 2 -1)
Hence, the linearized problem is written with:

C = -c L 0 0 c R ; R = c L -c R 1 1
In section 2.4.2, the expression of the wave strengths vector δ has been established for general linear systems of dimension 2 (see equation (2.70)):

δ = 1 c R + c L c R ∆F + ∆v c L ∆F -∆v
leading to the solution Q between the two discontinuous waves:

Q = Q L + δ 1 R 1 = v L F L + δ 1 c L 1 or Q = Q R -δ 2 R 2 = v R F R -δ 2 -c R 1
Substitution of δ 1,2 from the second equations into the first provides straight line equations in the phase plane (F, v):

v = v L + c L (F -F L ) ; v = v R + c R (F R -F )
The intersection of those straight lines in the phase plane corresponds to the approximate solution. Figure

1 1.2 1.4 1.6 1.8 2 0 0.5 1 1.5 2 2.5 • • Q L Q R F v(m/s) (a) Initial data: F L = 1 ; F R = 2
Shock Rarefaction

• • Q L Q R F v(m/s) (b) Initial data: F L = 1 ; F R = 1.1 Shock Rarefaction Figure 2
.13: Comparison of approximate (dashed lines) and exact (solid lines) solution for a one-dimensional strain problem in a Saint-Venant-Kirchhoff hyperelastic material 2.13 shows comparisons of approximate and exact solutions for various initial data, all leading to a 1-shock-2-rarefaction exact solution. As expected, approximate and exact solutions are different and get closer for small initial discontinuities, satisfying the linearization assumption Q L ≈ Q R . As a consequence, a big initial discontinuity is considered in figure 2.13a so that the approximation error is larger than that of figure 2.13b for which initial data are based on a weak jump.

Conclusion

It has been seen in this chapter that solid dynamics balance equations can be written as a first order hyperbolic system whose theory has been recalled in section 2.1. Indeed, the thermodynamics framework assuming generalized standard materials combined with conservation laws allowed in section 2.2 the building of conservative or quasi-linear forms. Those systems of partial differential equations admit non-complex eigenvalues and independent eigenvectors provided that some requirements on the stored energy function are satisfied (positive definiteness of the acoustic tensor). Then, the characteristic analysis of the quasi-linear form in section 2.3 enabled the highlighting of specific wave types involved in the solutions of dynamic problems, that is: discontinuous, shock and simple waves. Even though exact solutions of linear and non-linear problems have been developed in section 2.4, it is not possible in general, hence the introduction of approximate-state Riemann solvers in section 2.5. This solution strategy will be used in what follows as an element of the Discontinuous Galerkin Material Point Method, which is the object of the next chapter.

Introduction

It has been highlighted in the previous chapter that even though exact solutions of hyperbolic systems have been derived, it is not possible in general. It is indeed well known that in addition to the mathematical complexity of PDEs, physics often involves multidimensional problems in domains with complex geometries that cannot be solved analytically. Numerical strategies may then be employed in order to compute an approximate solution to such problems. Since the early 50's, plenty of numerical methods consisting mainly of mesh-based schemes, that is methods subdividing a complex domain into elements of simple shapes in which an approximate solution is sought, have been developed (finite element [START_REF] Belytschko | Nonlinear finite elements for continua and structures[END_REF], finite volume [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF] etc.). Problems involving very large deformations may lead, however, to numerical difficulties when this kind of approach is used with a material description (Lagrangian formulation) due to severe mesh distortions. Alternatively, Eulerian methods avoid mesh entanglement by building an approximate solution of a PDE system on a fixed mesh that corresponds to a discretized control volume. Nevertheless, interface tracking techniques and convection steps are required in Eulerian approaches in order to follow the boundaries and transport internal variables, which is less convenient for solid than for fluid mechanics because of history dependent constitutive behaviors.

The Material Point Method (MPM) [START_REF] Sulsky | A particle method for history-dependent materials[END_REF] mixes the advantages of both Lagrangian and Eulerian methods in order to circumvent mesh entanglement. However, the MPM suffers from numerical dissipation and oscillations that make it unable to accurately capture discontinuous waves traveling in solids. These limitations have been the object of research that led to significant improvements of the method. Nevertheless, the MPM is still unable to capture discontinuities. This is the purpose of this chapter.

In what follows, a brief historical review of developments that led to the MPM formulation is made and the original formulation is recalled in section 3.1. Next, after emphasizing some shortcomings of the method, an extension of the MPM to the Discontinuous Galerkin (DG) approximation is proposed in section 3.2. At last, the numerical analysis of the Discontinuous Galerkin Material Point Method (DGMPM) is performed in section 3.3 in terms of convergence and stability. These studies of the numerical scheme show that the Courant number may be set to unity in specific cases, at the cost of first-order accuracy in velocity.

The material point method

The early developments that led to the original MPM started with the Particle-In-Cell method (PIC) for fluid dynamics problems [START_REF] Harlow | The particle-in-cell method for numerical solution of problems in fluid dynamics[END_REF]. The novelty brought by PIC was the representation of a fluid by a collection of moving particles inside a background control volume subdivided into cells. Every single particle is given a constant mass and a position which is updated based on the velocity field resulting from the solution of the discrete linear momentum balance equation on the fixed background mesh. On the other hand, energy and pressure are stored at cell centroids during the whole computation as element-wise constant fields (C 0 approximation). During the movement of Lagrangian points, fields are temporarily projected back and forth between the grid and particles so that the latter can transport information from one cell to another. This procedure is referred to as the convective phase. PIC thus enables the merging of Lagrangian and Eulerian techniques since the moving particles with ascribed masses allow the simulation of problems involving several highly deformed fluids without mesh distortion.

In spite of the good results provided by PIC, the numerical diffusion it suffers from has been addressed by two different ways. First, the method has been extended to a fully Lagrangian version [START_REF] Mccrory | Growth and saturation of instability of spherical implosions driven by laser or charged particle beams[END_REF] by storing not only mass and position but all the fields at particles. Second, a new projection procedure from the grid to particles has been proposed in order to reach second-order accuracy in space of the convective phase [START_REF] Nishiguchi | Second-order fluid particle scheme[END_REF]. Rather than ascribing cell velocity to particles as done in PIC, Lagrangian point velocities are updated based on the change of linear momentum resulting from the solution of the discrete system. Merging those two improvements yielded the so-called FLuid Implicit Particle method (FLIP) [START_REF] Brackbill | Flip: A method for adaptively zoned, particle-in-cell calculations of fluid flows in two dimensions[END_REF]. This approach enables the reduction of numerical diffusion but introduces noise in the vicinity of discontinuities [START_REF] Burgess | Mass matrix formulation of the flip particle-in-cell method[END_REF]. It is worth noticing that in this new PIC formulation, the background mesh is used for solving balance equations only, thus providing an adaptive feature to the numerical scheme.

Even though particles carry the whole history of the problem, FLIP has been essentially used until the 90's to model history-independent constitutive models which were dealt with at the grid level. The first application of the method to history-dependent materials was made in the context of solid mechanics [START_REF] Sulsky | A particle method for history-dependent materials[END_REF]. Lagrangian particles are then seen as material points since they are used to store every field of the problem and to compute constitutive equations, thus leading to the Material Point Method that is derived hereinafter.

Derivation of the MPM

Consider a solid domain with volume Ω t bounded by the surface ∂Ω t , subject to traction forces and prescribed velocity on its boundaries within the time interval τ = [0, T ]:

σ • n = T d on ∂Ω σ t v = ṽ on ∂Ω v t
where n is the outward normal vector to ∂Ω t and the set of boundaries satisfies

∂Ω t = ∂Ω σ t ∪ ∂Ω v t .

Weak formulation of the continuum problem

The weak form of the Eulerian balance equation of linear momentum (2.27) is written based on the Galerkin approach and the following function spaces:

V 1 = {u ∈ H 1 (Ω t )} ; V 1 h = {u ∈ V 1 |u = ṽ on ∂Ω v t } ; V 1 h,0 = {u ∈ V 1 h |u = 0 on ∂Ω v t } where H 1 (Ω t ) denotes the Sobolev space [1, Ch.4], [43, Ch.1]. Multiplication of equation (2.27) by a test function w ∈ V 1
h,0 and integration over Ω t yields, after integration by parts, the following weak form:

Find v ∈ V 1 h such that Ωt ρ v • w dv + Ωt σ : ∇w dv - ∂Ω σ t T d • w ds = Ωt ρb • w dv ∀ w ∈ V 1 h,0 , ∀ t ∈ τ (3.1)

The MPM discretization

The continuum body is discretized into a set of N p material points in an arbitrary Cartesian grid. The grid is here supposed to be made of N n nodes and E non-overlapping cells or elements. The set of faces separating cells that contain particles and empty ones defines the boundary of the mesh (see figure 3.1 for a two-dimensional case).

Analogously to the Finite Element Method (FEM) [START_REF] Belytschko | Nonlinear finite elements for continua and structures[END_REF], the velocity is approximated on the MPM background mesh by: where v i is the velocity of the ith grid node, and S i (x) the (linear) shape function attached to it. One feature of the original MPM is that the mass density is approximated in the grid based on the mass carried by each material point:

v(x, t) = S j (x)v j (t) (3.2) 3.1. THE MATERIAL POINT METHOD Ω t • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • boundary of the domain
ρ (x, t) = Np p=1 m p δ(φ(X p , t) -x) = Np p=1 m p δ(x p -x) (3.3)
where the Dirac delta function δ is often referred to as the characteristic function of material points, and x p is the position of the pth particle. In what follows, the dependence on time will be omitted for simplicity and p and (i, j) stand respectively for material points and nodes. By introducing successively in the weak form (3.1) the specific Cauchy tensor ρσ = σ and the approximation of mass density (3.3), the integration property of the delta function allows writing:

Find v ∈ V 1 h such that Np p=1 m p [ v(x p ) • w(x p ) + σ(x p ) : ∇w(x p )] = Np p=1 m p b(x p ) • w(x p ) + ∂Ω σ t T d • w ds ∀ w ∈ V 1 h,0 , ∀ t ∈ τ
Then, with the finite element approximation (3.2), the weak form reads:

Find v ∈ V 1 h such that w i Np p=1 m p S ip S jp vj + ∇S ip • σp = w i   Np p=1 m p S ip b p + ∂Ω σ t S i (x)T d ds   ∀ w ∈ V 1 h,0 , ∀ t ∈ τ (3.4)
in which S ip = S i (x p ) and σp = σ(x p ). Next, the arbitrariness of the test function w leads to the semi-discrete equation in matrix form:

M ij vj = f i int + f i ext (3.5)
where the definition of the mass matrix M ij and internal and external forces vectors f i int and f i ext are:

M ij = Np p=1 m p S ip S jp (3.6a) f i int = - Np p=1 m p ∇S ip • σp (3.6b) f i ext = Np p=1 m p S ip b p + ∂Ω σ t S i (x)T d ds (3.6c)
In addition, note that particles play the role of integration points in the semi-discrete form so that the mass matrix depends on the positions of material points in the grid and must be computed at each time step. Hence, the MPM can be seen as an extension of FEM with moving Gauss' points. However, this quadrature rule results in a consistent mass matrix M ij that may be singular when only one material point lies in an element due to reduced integration. The use of a linear combination of the (positive-definite) diagonally lumped mass matrix M L i = p S ip m p , and the (positive semi-definite) consistent mass matrix is then recommended [START_REF] Love | An unconditionally stable, energy-momentum consistent implementation of the material-point method[END_REF]. Nevertheless, since no parameter values are prescribed for this combination, the lumped mass matrix is widely used in MPM simulations in spite of additional dissipation of kinetic energy it introduces [START_REF] Burgess | Mass matrix formulation of the flip particle-in-cell method[END_REF].

At last, system (3.5) is solved for nodal accelerations vi , from which time integration allows updating nodal velocities. Hence, the time interval τ is discretized in N T sub-intervals of size ∆t n such that N T n=1 ∆t n = T and time integration is performed with an explicit forward Euler algorithm. The discrete equations are thus written:

v i,n+1 -v i,n ∆t n = vi (3.7)
where the superscripts • k,l denote the time step l at which the kth nodal field is evaluated.

Recall that the background grid inherited from FLIP is arbitrary so that fields must be projected between particles and nodes to (i) solve the linear system (3.5) and (ii) advect the solution by moving material points. The first mapping step, required to build the discrete form, aims at satisfying the conservation of linear momentum:

M L,n i v i,n = Np p=1 S ip m p v p,n (3.8)
that must be solved for every v i . Once the nodal accelerations are calculated from system (3.5), nodal velocity increments can be computed to update the velocity of each material point by following FLIP approach [START_REF] Nishiguchi | Second-order fluid particle scheme[END_REF]:

v p,n+1 = v p,n + ∆t n Nn i=1 S ip vi (3.9)
On the other hand, updated nodal velocities resulting from the solution of (3.7) are used to update the particle positions:

x p,n+1 = x p,n + ∆t n Nn i=1 S ip v i,n+1 (3.10)
Remark 9 Equations (3.9) and (3.10) lead to a double definition of the velocity since the time derivative of the displacement of a given particle does not correspond to its velocity. This can be seen in a finite differences sense by combining equations (3.10) and (3.7):

x p,n+1 -x p,n ∆t n = Nn i=1 S ip v i,n+1 = Nn i=1 S ip v i,n + ∆t vi (3.11)
Next, on gets from equation (3.9):

∆t n Nn i=1 S ip vi = v p,n+1 -v p,n
which, once introduced in (3.11), yields:

x p,n+1 -x p,n ∆t n = Nn i=1 S ip v i,n + v p,n+1 -v p,n
Considering the mapping (3.9), the right-hand side of this equation does not identify to the updated material point velocity, it then comes out that ẋp,n+1 = v p,n+1 . The computation of constitutive equations is performed at material points within the MPM, thanks to the updated nodal velocity field and the gradient of the shape functions. One has therefore some freedom since stress can be updated right after either the resolution of (3.7) or the projection (3.8). The first option yields the Update Stress Last (USL) algorithm while the second implementation is called Update Stress First (USF) [START_REF] Bardenhagen | Energy conservation error in the material point method for solid mechanics[END_REF].

Boundary conditions enforcement

Writing the weak form of linear momentum balance on the grid implies that Neumann boundary conditions are enforced on nodes through the external forces vector (3.6c). Next, since deformation gradient is updated by means of the nodal velocity field, Dirichlet boundary conditions are also dealt with at nodes. For USL, it is done once the solution is updated on the grid so that stresses can be computed and material point kinematics can be updated. On the other hand, the USF formulation first requires applying Dirichlet boundary conditions to compute the stress, and then again after solving discrete equations for particle kinematics [START_REF] Buzzi | Caveats on the implementation of the generalized material point method[END_REF]. Furthermore, the enforcement of boundary conditions is still a challenging aspect of the material point method (as the other mesh-free methods) for problems involving complex geometries, see for instance [START_REF] Michael | Imposition of essential boundary conditions in the material point method[END_REF].

Solution scheme summary

All the ingredients derived above are now summarized to write the procedure followed within the MPM for one time step. Position and velocity vectors x n and v n are supposed known at every material point that discretize the continuum body in a computational grid at time t n . USL formulation further requires the knowledge of the specific stress tensor σn at material points. The MPM solution scheme then consists of the following steps: USF formulation: enforcement of Dirichlet boundary conditions on the mesh and integration of constitutive equations so that σn is computed at material points.

(c) Evaluation of internal forces from equation (3.6b).

(d) Solution of the semi-discrete and discrete forms (3.5) and (3.7) so that nodal accelerations vi and updated velocities v i,n+1 are determined. Enforcement of Dirichlet boundary conditions.

(e) Update material point velocities and positions with equations (3.9) and (3.10) respectively. At this point the mesh has virtually moved, but since fields have been transferred back to particles, the underlying grid can be discarded and rebuilt for computational convenience, thus involving the projection (b) at the subsequent time increment.

USL formulation: computation of constitutive equations so that σn+1 is stored at material points.

Other Lagrangian formulations of the MPM

The above formulation of the MPM, which is widely used, provides a material description of the motion while being based on Eulerian quantities. Alternatively, one can solve a discrete system written on a reference configuration, provided the use of a suitable stress measure [START_REF] Love | An unconditionally stable, energy-momentum consistent implementation of the material-point method[END_REF]. Although both formulations lead to a material description of the deformation, they are distinguished hereinafter by referring to them as the Eulerian formulation and the Lagrangian formulation respectively. The derivation of the latter approach is now briefly recalled.

Discretized equations

An approximate solution of the balance equation of Lagrangian linear momentum is sought within the time interval τ = [0, T ], in a solid with volume Ω 0 and boundary ∂Ω 0 in the reference configuration. The body is subject to traction forces and prescribed velocity on its boundaries:

Π • N = T d,0 on ∂Ω Π 0 v = ṽ on ∂Ω v 0
where N is the outward normal vector to ∂Ω 0 and the set of boundaries satisfies ∂Ω 0 = ∂Ω Π 0 ∩ ∂Ω v 0 . Thus, the weak form of equation (2.26) is written by means of the Galerkin method. After integration by parts, one gets:

Find v ∈ V 1 h such that Ω0 ρ 0 v • w dv + Ω0 Π : ∇ 0 w dv - ∂Ω Π 0 T d,0 • w ds = Ω0 ρ 0 b • w dv ∀ w ∈ V 1 h,0 , ∀ t ∈ τ
The solid is represented in the reference configuration with a collection of N p material points whose masses, combined with the Dirac delta function, define the reference mass density: ρ 0 (X) = Np p=1 m p δ (X p -X). Thus, the specific PK1 stress tensor ρ 0 Π = Π is introduced so that the weak form reads:

Find v ∈ V 1 h such that Np p=1 m p vp • w p + Πp : ∇ 0 w p - ∂Ω Π 0 T d,0 • w ds = Np p=1 m p b p • w p ∀ w ∈ V 1 h,0 , ∀ t ∈ τ
Next, shape functions defined at nodes of a background grid and arbitrariness of the test field w lead to the semi-discrete form:

M ij vj = f i int + f i ext (3.12)
Since the particles are considered as integration points, the diagonally lumped mass matrix is used in order to circumvent an eventual singularity. Hence, the lumped mass matrix and forces vectors involved in the semi-discrete system (3.12) equation read:

M L i = Np p=1 m p S ip (3.13a) f i int = - Np p=1 m p ∇ 0 S ip • Πp (3.13b) f i ext = Np p=1 m p S ip b p + ∂Ω Π 0 S i (x)T d,0 ds (3.13c)
At last, the explicit time discretization of equation (3.12) yields the discrete system that is solved at nodes to advance the velocity from time increment t n to t n+1 :

M L i v j,n+1 -v j,n ∆t n = f i int + f i ext (3.14)

Solution scheme

Suppose that position and velocity vectors, as well as the specific PK1 tensor, are stored at material points at time t n in the volume Ω n , that is: x p,n , v p,n and Πp,n are known ∀p = 1, ..., N p . In the reference configuration, which may be defined as the the last computed configuration according to [START_REF] Love | An unconditionally stable, energy-momentum consistent implementation of the material-point method[END_REF], one first has to compute the mass matrix (3.13a). This updated Lagrangian formulation of the MPM bases the computation of the mass matrix on the shape functions evaluated at particle positions x p,n . Alternatively, the reference configuration may be kept constant during simulations so that the mass matrix is computed once and for all based on X p . It is however worth noticing that such a total Lagrangian formulation of the MPM does not, to our knowledge, appear in the literature.

Next, analogously to the Eulerian formulation of the MPM, the mapping of the velocity from particles to nodes aims at ensuring the conservation of Lagrangian linear momentum by solving the following system for the nodal velocities:

M L i v i,n = Np p=1 S ip m p v p,n
Once the force vectors are built from equations (3.13b) and (3.13c), acceleration and velocity are computed on the grid by successively solving systems (3.12) and (3.14). Then, according to the FLIP mapping (3.9), these fields are used to update the material point kinematics (i.e. compute v p,n+1 and x p,n+1 ). Moreover, the gradient of the shape functions provides the deformation gradient thanks to the following explicit time discretization:

F p,n+1 -F p,n ∆t = ∇ 0 v p,n+1 (3.15)
The above equations yield directly the updated deformation gradient for the total Lagrangian formulation. On the other hand, the use of the velocity gradient tensor is made in the updated Lagrangian formulation, which is approximated as:

∇ n v ≈ Np p=1 ∇S ip (x p,n ) ⊗ v i,n
where ∇ n (•) = ∂v ∂x x=φ(X,t n ) denotes the gradient operator in the configuration Ω n . The chain rule then allows identifying this tensor in the right-hand side of equation (3.15):

∇ 0 v p,n+1 = ∇ n v p,n+1 • ∂φ(X, t n ) ∂X = ∇ n v p,n+1 • F p,n
so that the deformation gradient at time t n+1 in the updated Lagrangian formulation results from:

F p,n+1 -F p,n ∆t = ∇ n v p,n+1 • F p,n
Finally, constitutive equations can be dealt with at the particle level. Note that the choice between USL and USF algorithms for the integration of constitutive equations also arises with Lagrangian formulations developed above.

Although the material point method presented above has been successfully applied to a wide range of complex engineering problems involving finite deformations [START_REF] Wiȩkowski | The material point method in large strain engineering problems[END_REF], this approach suffers from limitations which are discussed below.

Shortcomings of the MPM

The grid-crossing instability

The research on MPM mainly focused so far on the instability arising from the computation of internal forces when material points move from one cell to another, the so-called grid-crossing error [START_REF] Bardenhagen | The generalized interpolation material point method[END_REF]. While particle positions have no influence in total Lagrangian approaches, the discontinuity of the gradient the of shape functions across element interfaces yields tremendous noise in the numerical solution of both Eulerian and updated Lagrangian MPMs. The Generalized Interpolation Material Point Method (GIMP) addresses this issue by modifying the particle characteristic function, thus widening the domain of influence of material points [START_REF] Bardenhagen | The generalized interpolation material point method[END_REF]. By doing so, every particle is given a domain having a constant shape or not, that must be tracked during the computation, leading respectively to the uniform GIMP (uGIMP) or contiguous particle GIMP (cpGIMP) formulations [START_REF] Zhang | The Material Point Method: A Continuum-based Particle Method for Extreme Loading Cases[END_REF]. In the cpGIMP algorithm, diagonal entries of the deformation gradient are used to update the "shapes" of particles. However, these methods involve other difficulties related to the tracking of the deforming domains of material points in the Eulerian grid, and may also suffer from mesh entanglement [START_REF] Zhang | Material point method enhanced by modified gradient of shape function[END_REF]. The Convected Particle Domain Interpolation (CPDI) [START_REF] Sadeghirad | A convected particle domain interpolation technique to extend applicability of the material point method for problems involving massive deformations[END_REF], accounts for shear deformations and rotations of particle domains and also considers alternative grid basis functions, thus providing an improvement of GIMP. Other approaches consisting in the direct modification of the approximation basis have then been followed to tackle the grid-crossing error. First, since the midpoint rule leads to significant quadrature error for discontinuous functions over particle domains, the use of quadratic or cubic B-Spline as nodal shape functions has been proposed in [START_REF] Steffen | Analysis and reduction of quadrature errors in the material point method (mpm)[END_REF]. Then, the B-Spline Material Point Method (BSMPM) [START_REF] Tielen | High-order material point method[END_REF], based on continuous gradients, allows circumventing the grid-crossing instability. Second, the classical approximation basis is used and modified gradients with extended supports are introduced in the Dual Domain Material Point Method (DDMPM) [START_REF] Zhang | Material point method enhanced by modified gradient of shape function[END_REF]. BSMPM and DDMPM both solve the grid-crossing error by widening the domain of influence of material points since B-Spline functions span more cells than Lagrange polynomials, and the support of the modified gradient of DDMPM is larger than that of the shape functions.

Remark 10 Since it enables to get rid of the grid-crossing instability, the total Lagrangian formulation will be used for all applications of the MPM in this work.

The back-mapping instability

Numerical oscillations that are due to the updating of the velocity field at material points [START_REF] Burgess | Mass matrix formulation of the flip particle-in-cell method[END_REF] appear in the MPM solution, even for problems that do not involve the grid-crossing instability. Although a total Lagrangian formulation is considered, such a mapping of the velocity is required in order to, for instance, employ meshadaption techniques. As an illustration, consider a one-dimensional elastic bar of length L = 6 m with Young's modulus E = 200 GP a and mass density ρ = 7800 kg • m -3 within the infinitesimal theory. The bar is initially in a free stress state and Riemann-type initial velocities are prescribed along the bar, that is:

v = v 0 > 0 for x ∈ [0, L/2[ and v = -v 0 for x ∈]L/2, L].
Both ends of the bar are traction free so that this problem is equivalent to an impact of two elastic bars moving toward each other. The bar is discretized by using the MPM in a computational grid composed of 50 regular cells, each containing either one or two material points (1ppc or 2ppc). Single material points are centered in cells while particles sharing an element are placed symmetrically with respect to elements centers and are regularly spaced. The exact solution of this problem (see section 2.4.2) consists of two elastic discontinuities emanating from the middle of the bar and propagating leftward and rightward. In figure 3.2, numerical results provided by the USL and USF implementations of the MPM are compared to the exact solution before reflection of the waves on the boundaries. Note that here and in what follows, the stress and the velocity are extracted at material points. Figures 3.2a and 3.2b show that both USL and USF solutions oscillate once the elastic front passed. Those oscillations are much more significant in the USF solutions regardless of the number of particles used. Indeed, the noise in USL results quickly reduces so that the correct stress level is reached in the middle region of the bar. Moreover, even though the 1ppc and 2ppc discretizations provide different results (slightly different for USL), neither of them enables the removal of oscillations. Figure 3.2c shows, on the other hand, the evolution of numerical total energies resulting from an MPM approximation of kinetic and strain energies:

e kin = 1 2 Ω ρv • v dΩ ≈ 1 2 Np p=1 m p v p • v p e strain = 1 2 Ω ρσ : ∇u dΩ ≈ 1 2 Np p=1 m p σp : ∇u p
One can see that the USL formulation dissipates more energy than the USF for this problem. These results differ from observations made in [START_REF] Bardenhagen | Energy conservation error in the material point method for solid mechanics[END_REF] in which no significant differences for problems involving smooth solutions are noticed. At last, figure 3.3 shows the numerical velocities at various times. The same observations as before are made by looking at figures 3.3a and 3.3b, namely numerical oscillations in both USL and USF solutions with bigger peaks in the USF one. Furthermore, a notable numerical artifact that prevents the proper assessment of the velocity occurs at the middle of the bar when two particles per cell are used. This can be partially explained by the oscillating stress field that introduces noise in internal forces as well. However, the velocity field given by the USL computation "locks" in the central region although the stress field converges in time to the exact solution (figures 3.2a and 3.2b), which implies that this issue is not due to the computation of internal forces, but most likely to the convective phase. Since the back-mapping is known to be responsible for noise introduced in the MPM solution [START_REF] Burgess | Mass matrix formulation of the flip particle-in-cell method[END_REF], we propose to project the velocity onto particles by a simple interpolation, as it is done in PIC:

v p,n+1 = Nn i=1 S ip v i,n+1 (3.16)
Particle positions are next updated as in FLIP, with the nodal velocity (3.10), thus amounting to deforming Lagrangian particles with the grid and exactly satisfying the relation ẋp = v p , which is not the case in MPM. It is then worth noticing that the weak form (3.4) is based on the definition (3.16), this projection seems therefore more consistent. The comparison between MPM using either PIC or FLIP mapping and exact velocities is plotted in figures 3.4a and 3.4b As can be seen, the PIC projection eliminates the noise in numerical solutions as well as the velocity locking in the central region of the bar. In addition, stress results provided by the USL using the PIC mapping are also close to the exact solution (see figures 3.5a and 3.5b). However, the resolution of discontinuities (i.e. the sharpness of the solution) allowed by the "modified" USL-MPM is slightly lower than that of the original formulation owing to numerical dissipation as shown in the energy plot in figure 3.5c. This result was expected since the FLIP projection was preferred to PIC one for it enabled a reduction of numerical dissipation [START_REF] Nishiguchi | Second-order fluid particle scheme[END_REF].

Strategy for reducing oscillations and diffusion

We are concerned with the accurate solution of hyperbolic problems in solids undergoing finite deformations.

Although the MPM enables an efficient management of large strains, the oscillations it suffers from do not allow accurate capture of waves propagating in a medium. The above numerical simulations however suggest that noise can be removed by using an interpolation instead of an update of material point velocity, at the cost of additional diffusion. Thus, the use of a linear combination of PIC and FLIP projections to update particles velocity was proposed in [START_REF] Nairn | Numerical simulation of orthogonal cutting using the material point method[END_REF], which can be seen as an artificial damping. Further, a tunable mapping procedure based on a parameter m that completely eliminates noise in the MPM solution has been proposed in the Extended PIC of order m (XPIC(m)) [START_REF] Hammerquist | A new method for material point method particle updates that reduces noise and enhances stability[END_REF]. A classical interpolation is selected for m = 1 whereas the mapping tends to FLIP as m → ∞ while avoiding oscillations without dissipating energy. Nevertheless, the method is still unable to capture discontinuities due to numerical diffusion. The point of view adopted in this thesis is that numerical diffusion is essentially due to the projection of fields from nodes to particles which spreads the information. Hence, a reduction of the domain of influence of nodes might result in lower dissipation. As a consequence, the reintroduction of the PIC mapping combined with the Discontinuous Galerkin approximation (DG) within the MPM is proposed here. These two features are expected to respectively avoid spurious oscillations and limit numerical diffusion.

Extension of the MPM to discontinuous Galerkin approximation

After a brief historical review of DG methods, the Discontinuous Galerkin Material Point Method is derived within the large strain framework with a total Lagrangian formulation. It will be seen that this new numerical approach makes use of the approximate-state Riemann solver developed in section 2.5 to compute intercell terms whose purpose is to connect elements together. At last, the DGMPM solution scheme will be provided for hyperbolic problems.

The discontinuous Galerkin approximation

The DG approximation was first introduced in the context of the finite element method for the solution of the neutron transport equation [START_REF] Reed | Triangular mesh methods for the neutron transport equation[END_REF]. This hyperbolic equation describes the advection of the angular flux which describes the probable number of neutrons at a given location. Since neutrons can lie in a cell of a finite element mesh while its neighbors are empty, the need for a discontinuous approximation of the primal field across element interfaces within a FEM context arises. Hence, an approximate solution was sought by the Galerkin method, in a domain discretized with triangular elements by means of Lagrange polynomials that can be discontinuous across the cells. This approach amounts to duplicate the nodes of the mesh so that the support of each shape function reduces to one finite element. This early work have launched a series of developments of the Discontinuous Galerkin Finite Element Method (DGFEM) for parabolic [START_REF] Arnold | An interior penalty finite element method with discontinuous elements[END_REF], elliptic [START_REF] Hansbo | Discontinuous galerkin methods for incompressible and nearly incompressible elasticity by nitsche's method[END_REF][START_REF] Noels | A general discontinuous galerkin method for finite hyperelasticity. formulation and numerical applications[END_REF], and hyperbolic problems [START_REF] Cockburn | Discontinuous galerkin methods for convection-dominated problems[END_REF][START_REF] Vilar | A high-order Discontinuous Galerkin discretization for solving two-dimensional Lagrangian hydrodynamics[END_REF]. Indeed, the DGFEM gained more and more popularity since the 80's, even for problems that do not involve discontinuities, on account of its ability to locally handle high-order approximation and its highly parallelizable nature. Research conducted in the context of hyperbolic problems, of particular interest here, enabled the introduction of numerical tools developed for Finite Volume Methods (FVM) within finite element schemes. Namely, the use of suitable slope limiters [START_REF] Van Leer | Towards the ultimate conservative difference scheme. ii. monotonicity and conservation combined in a second-order scheme[END_REF] based on the total variation [START_REF] Harten | High resolution schemes for hyperbolic conservation laws[END_REF] enables the formulation of flexible numerical methods in which good resolution of discontinuities is possible without destroying the accuracy in smooth regions. Furthermore, these approaches can easily handle mesh-adaption strategies thanks to the relaxation of the fields continuity. Nevertheless mesh tangling problems do not vanish. Thus, the introduction of DG approximation in the MPM should lead to a numerical method that benefits from both FEM and FVM features and enables local high-order approximation while avoiding mesh entanglement instabilities.

Derivation of the DGMPM

Consider again a continuum solid body with volume Ω t within the time interval τ . The DGMPM is expected to provide a material description of a deformation so that an approximate solution of a Lagrangian system of conservation laws written in conservative form is sought. Recall that such a conservative form for some vector of conserved quantities U reads, in Cartesian coordinates:

∂U ∂t + D α=1 ∂F • E α ∂X α = S ∀X, t ∈ Ω 0 × τ (3.17)

The DGMPM discretization

As for MPM, a continuum body Ω t is discretized within the time interval τ into a set of N p material points in an arbitrary Cartesian grid made of N n nodes and E non-overlapping cells of volume Ω e . The boundary of the domain is again defined by the set of edges separating empty cells from those containing particles (see figure 3.1 for a two-dimensional example). In addition, the reference mass density is described in the computational grid by means of the Dirac delta function and particle masses:

ρ 0 (X) = Np p=1 m p δ (X p -X) (3.18)
In a similar manner to FEM and MPM, the vector of conserved quantities is approximated on the mesh by:

U(X, t) = Nn i=1 S i (X)U i (t) (3.19)
with U i the vector of conserved quantities at node i, and S i (X) the shape function attached to it. Note that the convention of denoting particle and nodal fields by p and (i, j) still holds in this section.

Weak formulation of the continuum problem

Multiplying equation (3.17) by a test function V yields the weak formulation of the problem:

Find U ∈ V 1 h such that Ωt ∂U ∂t V dΩ + Ωt ∂F α ∂X α V dΩ = Ωt SV dΩ ∀ V, t ∈ V 1 h × τ
The key idea of DG methods is to allow the jump of fields across mesh elements faces by using broken polynomial spaces for the approximate solution [43, Ch.1]:

V k = {V ∈ H k (Ω e )} ; V k h = {V ∈ P k (Ω e )} ⊂ V k with H k (Ω e
), the Sobolev space and P k (Ω e ), the space of polynomials of degree k in Ω e . We restrict our attention here to linear polynomials (k = 1). These broken polynomial spaces allow the rewriting of the weak form element-wise. After integration by parts, one gets:

Find U ∈ V 1 h such that Ω e ∂U ∂t V dΩ - Ω e F α ∂V ∂X α dΩ + ∂Ω e (F • N ) V dΓ = Ω e SV dΩ ∀ V, e, t ∈ V 1 h × [1, E] × τ
where ∂Ω e is the boundary of the eth element with outward normal vector N . The dot operator F • N denotes the inner product between the outward normal vector and every component of the flux, thus yielding the intercell flux, written F N for simplicity. Next, the introduction of specific fields:

U = ρ 0 Ū ; F α = ρ 0 Fα ; S = ρ 0 S
combined with the definition of mass density (3.18), leads to the following Lagrangian formulation:

Np p=1 m p ∂ Ū ∂t V -Fα ∂V ∂X α -SV |X=X p + ∂Ω e F N V dΓ = 0 ∀ V, e, t ∈ V 1 h × [1, E] × τ (3.20)
As a first development step of the method, only the total Lagrangian formulation is considered in order to eliminate the grid crossing instability. Introduction of the DGMPM approximation (3.19) and arbitrariness of the test field in the weak form (3.20) finally provide the semi-discrete system that must be solved on the grid:

Np p=1 S ip m p S jp ∂ Ūj ∂t - ∂S ip ∂X α m p S jp Fj α -S ip m p Sp + ∂Ω e S i (X)F N dΓ = 0 ∀ e, t ∈ [1, E] × τ (3.21)
or, in matrix form:

M ij ∂ Ūj ∂t -K α ij Fj α -S i + Fi = 0
Here again, particles play the role of integration points in volume integrals owing to the Dirac delta function. Hence, the consistent mass matrix M ij may also be singular due to reduced integration so that the diagonally lumped mass matrix M L i is used. [START_REF] Jacobs | High-order nodal discontinuous galerkin particle-in-cell method on unstructured grids[END_REF] and [START_REF] Stindl | Comparison of coupling techniques in a high-order discontinuous galerkin-based particle-in-cell solver[END_REF] in which different projections of fields between the grid and particles are used. Although those methods allow local high-order approximation, particles do not carry every field so that the DGPIC, as the original PIC, cannot be considered as a fully Lagrangian approach. In addition, the use of the Gauss quadrature rule for volume integrals of the weak form makes the approach different from that developed in the following.

Remark 11 An extension of PIC to DG approximation for the solution of Maxwell's equations is proposed in

The discrete system is derived by discretizing the time interval τ into N t subintervals and using the explicit forward Euler method:

M L i Ūi,n+1 - Ūi,n ∆t n = K α ij Fj,n α + S i,n - Fi,n (3.22)
where again, the superscripts (•) k,l denote a field evaluated at node k and time step l. Note that in general the source term S may depend on the vector of conserved quantities, hence the superscript n in equation (3.22). Alternatively, a second-order Runge-Kutta (RK2) explicit time discretization may be employed, leading to the following two-stage discrete form:

M L i Ūi,n+1/2 - Ūi,n ∆t n = 1 2 K α ij Fj,n α + S i,n - Fi,n M L i Ūi,n+1 - Ūi,n ∆t n = K α ij Fj,n+1/2 α + S i,n+1/2 - Fi,n+1/2 (3.23)
Remark [START_REF] Richter | An explicit finite element method for the wave equation[END_REF] We chose here one existing two-stage second order Runge-Kutta method among others. See for instance [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF]Sec. 10.4.2] for a Total Variation Diminishing version of the RK2 time discretization.

Non-homogeneous hyperbolic system

Solid mechanics equations may lead to source terms in the conservative form even for neglected body forces. For instance, the writing of equations in cylindrical and spherical coordinates systems, in which the gradient operators involve terms that are not derivatives, yields a right-hand side in system (3.17) that depends on U and called a geometric source term [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF]Ch.17]. This is also the case for relaxation systems as that identified in section 2.2.4 for elasto-viscoplasticity. In the latter case, a small relaxation time τ compared to the time scale governed by the convective part (i.e. the characteristic speeds of the problem) leads to a stiff system [START_REF] Jin | Runge-kutta methods for hyperbolic conservation laws with stiff relaxation terms[END_REF]. The solution of this class of problems on underresolved grids, that is τ ∆t, requires the use of asymptotically convergent, accurate and stable numerical methods in order to avoid nonphysical results [START_REF] Leveque | A study of numerical methods for hyperbolic conservation laws with stiff source terms[END_REF]. Thus, the numerical scheme must enable the use of underresolved grids whose time step is dictated by the Courant-Friedrichs-Lewy (CFL) stability condition. The method also has to be asymptotic preserving in the limit τ → 0 for fixed ∆t and ∆X [START_REF] Jin | Runge-kutta methods for hyperbolic conservation laws with stiff relaxation terms[END_REF], that is the limiting scheme is a good discretization of the equilibrium system even if the source term is underresolved. At last, the order of accuracy in the stiff limit and stability at the discrete level must be preserved [START_REF] Pareschi | Implicit-explicit runge-kutta schemes and applications to hyperbolic systems with relaxation[END_REF], and the numerical scheme must be able to preserve the steady state numerically [START_REF] Dumbser | Finite volume schemes of very high order of accuracy for stiff hyperbolic balance laws[END_REF].

A commonly used approach to solve non-homogeneous systems consists in solving alternatively a homogeneous PDE system and a system of ODEs, namely: (∆t) and F (∆t) respectively, the discrete solution reads [START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics: a practical introduction[END_REF]:

∂U ∂t + D α=1 ∂F • E α ∂X α = 0 (3.24a) dU dt = S (3.
U n+1 = F (∆t) H (∆t) (U n ) (3.25)
The fractional-step method (3.25), known as Godunov's splitting, is only first-order accurate in time when H and F are at least first-order accurate solution operators. On the other hand, Strang splitting:

U n+1 = F (∆t/2) H (∆t) F (∆t/2) (U n ) (3.26)
is second-order accurate if each solution operator is at least second-order accurate [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF].

Two sub-problems are thus solved separately at each time step, the solution of the first one being used as initial conditions in the second. Thus, fractional-step methods (3.25) and (3.26) enable to take advantage of efficient tools already developed for homogeneous systems of conservation laws and of ODE solvers. However, it is worth noticing that even though the use of an implicit solver for equation (3.24b) eliminates the influence of the relaxation time on the CFL condition, an L-stable solver as the Euler implicit method [START_REF] Ascher | Computer methods for ordinary differential equations and differentialalgebraic equations[END_REF] is required to compute the correct solution.

Intercell fluxes

The DGMPM solution scheme reduces to the discrete solution operator H (∆t) , or H for simplicity, which in turn, aims at finding approximate similarity solutions since it only solves homogeneous systems (recall remark 4 in section 2.3). For a DGMPM space-time discretization made of N p material points and N T time increments, those solutions can be written:

Q p,n+1 = H Q j,n p = 1..., N p ; j = 1..., N n ; n = 0, ..., N T -1 (3.27)
where the set of nodes j having an influence on Q p defines the stencil of the method.

Definition 4 A numerical scheme is said monotone if it satisfies:

∂H ∂Q j ≥ 0 ∀j The following statement then holds: Theorem 3 (Godunov [START_REF] Godunov | Finite difference method for numerical computation of discontinuous solutions of the equations of fluid dynamics[END_REF]) Monotone linear numerical schemes can be at most first-order accurate.

Intercell fluxes of the weak form propagate information across cells by taking into account the different values that fields can take on each side of the interface. DG methods for hyperbolic problems are based on the requirement of ensuring monotonicity of the scheme for piecewise constant approximations [START_REF] Cockburn | Discontinuous galerkin methods for convection-dominated problems[END_REF]. Such a numerical method is monotone for flux functions F N that are Lipschitz continuous, consistent and monotone, namely, they must be E-fluxes [START_REF] Osher | Riemann solvers, the entropy condition, and difference[END_REF]. One possibility, which is widely used and adopted here, is the Godunov flux function.

The Godunov flux

The Godunov method [START_REF] Godunov | Finite difference method for numerical computation of discontinuous solutions of the equations of fluid dynamics[END_REF] has been proposed in the context of finite difference schemes in which the piecewise constant approximation of the solution naturally allows the definition of a local Riemann problem at cell interfaces. That is, two cells i and i + 1 with interface having normal vector N , define a Riemann problem in the direction X N = X • N whose stationary solution is used to compute the intercell numerical flux F • N . The stationary solution U * is the similarity solution along the vertical characteristic in the (X N , t) plane of the Riemann problem:

∂U ∂t + ∂F N ∂X N = 0 U(X N , 0) = U X - N if X N < 0 U X + N if X N > 0 (3.28)
where U X - N and U X + N are the states lying infinitely close to the interface in cells i and i + 1 respectively. Godunov's method hence allows accounting for the complete wave structure of the solution within the numerical scheme. Since this method is known to be based on E-fluxes, intercell fluxes involved in boundary integrals of the weak form of DG-methods can be computed as Godunov's ones. To this end, Riemann problems are defined at cells faces by considering that initial data are piecewise constant even for high-order approximations. Note that this strategy is also followed in FVM when high-order reconstruction techniques of the fields are used.

Let us consider as an illustration of the computation of interface fluxes within the DGMPM, the twodimensional case depicted in figure 3.6. One Riemann problem per interface is considered in order to avoid a dramatic increase in computational times. Thus, by averaging nodal fields on each side of the interface, one obtains mean downwind and upwind states U X - N and U X + N that correspond to the initial conditions of the Riemann problem. Furthermore, as mentioned in section 2.5, the exact (time consuming) solution of the Riemann problem is not necessarily desired since only the stationary solution is kept. As a consequence, the stationary states are approximated here by means of an approximate-state Riemann solver, and the corresponding Godunov fluxes are determined.

1 2 4 3 1 2 4 3 X N W L = W1+W4 2 W R = W2+W3

Remark 13

Hyperbolic systems having zero eigenvalues lead to stationary waves, propagating at zero celerity, across which the solution of the Riemann problem may have discontinuities. In that case, the associated Godunov flux is also discontinuous so that one stationary state must be considered on both sides of the characteristic. This approach leads to the computation of two fluxes, each contributing to one cell only.

Recall that hyperelasticity and elastoplasticity conservative forms (equations (2.49) and (2.54)) involve strains in the vector of conserved quantities, and stresses in the flux vector. Hence, the calculation of Godunov fluxes requires the integration of constitutive laws that can be time-consuming for non-linear problems. Nevertheless, the introduction of an auxiliary vector of conserved quantities Q and the Riemann problem it provides:

∂Q ∂t + J ∂Q ∂X N = 0 Q(X N , 0) = Q X - N if X N < 0 Q X + N if X N ≥ 0 (3.29)
avoids the computation of constitutive equations. Indeed, the auxiliary vectors introduced in section 2.2.4 are rearrangements of the flux components, that is F N (U) = LQ with L some localization matrix. Then, the intercell flux is F N (U * ) = LQ * where Q * is the stationary solution of Riemann problem (3.29). Note that this approach requires that stress as well as velocity and strain are projected back and forth between material points and nodes within the DGMPM.

Remark 14 Moreover, solving Riemann problems for the auxiliary vector Q * does not avoid the integration of constitutive equations on the grid when a multi-stage time integrator as the RK2 is used. Indeed, every subsequent step 0 ≤ k ≤ 1 of the time integrator requires an initial auxiliary vector Q n+k , which is computed from the vector of conserved quantities U n+k .

Transverse corrections

The method derived above for the computation of normal fluxes can be viewed as the Donor-Cell Upwind (DCU) method [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF] in which only contributions from upwind cells sharing an edge (in two dimensions) with the current one are considered. For multidimensional problems waves can travel in several directions such that contributions coming from corner cells must be taken into account in order to improve accuracy and stability of the numerical scheme. The Corner Transport Upwind (CTU) method [START_REF] Colella | Multidimensional upwind methods for hyperbolic conservation laws[END_REF] consists in considering contributions propagating in bias and coming from upwind cells sharing only a node (in two dimensions) with another. This approach allows improving the Courant condition especially for solid mechanics problems for which strain components are coupled through Poisson's effect. At each cell interface, one defines left-going and right-going fluctuations as:

A -(∆U) = F N (U * ) -F N (U X - N ) ; A + (∆U) = F N (U X + N ) -F N (U * )
Let's consider the patch of grid cells shown in figure 3.7, and focus on the edge denoted (i) whose local The spectral analysis of the corresponding Jacobian matrix, carried out in [START_REF] Kluth | Analyse mathématique et numérique de systèmes hyperélastiques et introduction de la plasticité[END_REF], leads to right eigenvectors that are be written R i W hereinafter. The negative normal fluctuation is, for instance, decomposed onto the characteristic basis associated with edge (j) as:

A i,+ (∆W) A i,-(∆W) (i) N i (k) (m) (j) (l) N k N j N m N l B k,+ A i,-(∆W) B j,+ A i,-(∆W) B m,+ A i,+ (∆W) B l,+ A i,+ (
A i,-(∆U) = M m=1 β m R j,m W where R j,m
W is based on the normal vector N j but also on different tangent moduli between grid cell L and its neighbor T. Since only waves with positive characteristic speeds with respect to the orientation defined by the outward normal vector to the considered edge will contribute to the transverse fluctuation, only the positive operator B + is used:

B j,+ A i,-(∆U) = M m=1 cm>0 c m β m R j,m W (3.30)
with c m , the celerity of wave m. An additional numerical flux defined at edges is hence built from these transverse fluctuations:

F j,tran = ∆t 2∆X j B j,+ A i,-(∆U)
which contributes to the flux between cell L and T (∆X j being the length of edge (j)). In-going transverse corrections must be counted positively while outgoing ones must be subtracted from intercell fluxes. At last, numerical fluxes computed at interfaces can be then integrated over cells faces in order to complete the discrete system.

Remark 15 For linear elasticity, there is no need for an auxiliary vector of conserved quantities so that transverse contributions are computed with the same eigenbasis as that of the approximate-state Riemann solver.

Boundary conditions enforcement

As in MPM [START_REF] Love | An unconditionally stable, energy-momentum consistent implementation of the material-point method[END_REF][START_REF] Steffen | Examination and analysis of implementation choices within the material point method (MPM)[END_REF], boundary conditions (BC) are treated at nodes that compose the domain boundary. By introducing ghost nodes on boundary interfaces, one can widen the use of the approximate-state Riemann solver in order to enforce both Dirichlet and Neumann BC. Similarly to finite volumes [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF], a state vector Q G is ascribed to those nodes so that the stationary solutions of Riemann problems are consistent with boundary conditions. Dirichlet BC enforcement for instance consists of (i) setting stress and free velocity components at ghost nodes equal to that of the associated interior node and (ii) solving the approximate Riemann problem for the velocity that must be enforced on ghost nodes knowing the stationary one that is equal to that of the BC. The same procedure holds for Neumann boundary conditions by interchanging velocity and stress.

Notice that the use of the auxiliary vector of conserved quantities so that state vectors contain stress allows the easy enforcement of Neumann boundary conditions. Transverse corrections further require the introduction

G c G B G L N G c G B G L N

Mesh nodes

Material points of corner ghost nodes equivalent to finite volume corner cells [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF]. Consider the two-dimensional case depicted in figure 3.8 in which boundary edges are represented by red lines. First, an inverse Riemann problem is solved between the corner ghost node G c and one regular ghost node, say G B , so that the stationary solution of the Riemann problem between those ghost nodes corresponds to the boundary condition holding on the vertical edge. Second, this procedure is repeated between the corner ghost G c node and G L in order to enforce the boundary condition holding on the bottom boundary between them.

Remark 16

The solution of the Riemann problem involving ghost nodes is made possible by extending material properties of the adjacent cell so that the characteristic structure of the solution can be computed. This implies that the deformation of interior nodes is duplicated at ghost nodes for problems such as hyperelasticity, whose eigenstructure depends on the deformation gradient. Hence, for such problems ghost nodes may carry stress and strain that are not related by constitutive equations. However, the deformation gradient at ghost nodes has no physical sense and is only used to compute the correct wave speeds.

DGMPM solution scheme

Let us assume that the vector of specific conserved quantities Ūn as well as the auxiliary vector Q n are known at every material point that discretizes a continuum body Ω 0 , in a grid made of N n nodes at a given time t n . The computational procedure followed within the DGMPM between two time steps n and n + 1 can now be summarized. We consider cases based on the use of an auxiliary vector, the others being only particular cases. The procedure then reads:

M L i Ūi = Np p=1 S ip m p Ūp and M L i Q i = Np p=1 S ip m p Q p (3.31)
to be solved for each Ūi and Q i respectively. The projection of fields from particles to nodes hence follows the weighted least squares interpolation used in MPM. (g) Back-mapping: as motivated at the end of section 3.1, the nodal updated solution is projected to material points with the classical interpolation as in PIC:

Ūp,n+1 = N i=1 S ip Ūi,n+1 (h)
Material point kinematics and constitutive model: The new solution Ūp,n+1 allows incrementing the deformation ϕ(X, t) and updating stress components, which will be used in the auxiliary vector for the next time step, through hyperelastic constitutive equations:

ϕ p,n+1 = X p + ∆tv p,n+1 Π p,n+1 = ∂Ψ ∂F (F p,n+1 )
The grid may then be discarded and reconstructed and in particular by means of adaptive algorithms applied in the reference configuration, in order to improve wave front tracking in the current one. The reconstruction of the grid is responsible for the presence of step (a).

Let's now recall or highlight significant differences between the DGMPM and the original MPM schemes. First, the use of conservation laws (3.17) instead of the momentum equation in the weak form implies that both velocity and gradients are solved at nodes making this new approach close to finite volume methods, which provides the same order of accuracy for both fields. Next, since the deformation gradient is no longer calculated with shape function gradients, the task of choosing between USF and USL algorithms vanishes. In that sense the DGMPM scheme is simpler. At last, the solution of Riemann problems at every edge of the mesh increases computational time. Fortunately, the use of discontinuous Galerkin approximation makes this numerical method highly parallelizable [START_REF] Cockburn | Discontinuous galerkin methods for convection-dominated problems[END_REF].

The numerical scheme derived above is analyzed in terms of stability and convergence in the next section.

Numerical analysis of the DGMPM

Following [START_REF] Hirsch | Numerical computation of internal and external flows: The fundamentals of computational fluid dynamics[END_REF], the DGMPM discretization of scalar linear advection problems are now written in a finite difference sense. The equations thus obtained are the starting point for von Neumann linear stability analyses. First, the one-dimensional problem is considered and the equations of the DGMPM space discretization combined with both forward Euler and RK2 explicit time integrations are derived. Second, the two-dimensional equations are written using the DGMPM space discretization along with the explicit forward Euler time integration only.

One-dimensional stability analysis

Model equation -Space discretization

We consider the scalar linear advection equation for an arbitrary quantity q = ρq moving at the constant speed a ∈ R + * in a homogeneous one-dimensional medium of length l:

∂ q ∂t + ∂ f ∂X = 0
with specific flux function f = aq, leading to the quasi-linear form:

∂ q ∂t + a ∂ q ∂X = 0 (3.32)
Equation (3.32) is discretized with the discontinuous Galerkin material point method. Thus, the medium is divided with N p material points arbitrarily distributed in E two-node elements of constant length ∆X (figure 3.9). The grid is such that at least one particle lies in every cell during the computation in order to ensure that there is no hole in the bar. Moreover, periodic boundary conditions are considered to simplify the analysis.

1 2 1 2k -1 2k k 2E -1 2E E Figure 3
.9: One-dimensional mesh made of E elements of constant length ∆X = l E . Since fields are carried by particles, we seek the scheme equation that gives the solution at a material point for a given time step, with respect to the solutions at other particles at the prior time step. In this section, Latin and Greek symbols are respectively devoted to nodes and material points. Since we consider here scalar quantities, the information on nodes and particles can be written as subscripts without ambiguity with vector components. In addition, to make the distinction between fields easier, upper and lower case Latin symbols are used for particle and nodal fields respectively. Hence, the solution at material point α and time step n reads Qn α whereas that of node i is qn i . Then, the cell containing the particle α will be denoted by c(α) so that the nodes interacting with this particle are 2c(α) -1 and 2c(α). At last, the linear shape functions defined in element c(α) are:

S 2c(α)-1 (X) = X 2c(α) -X ∆X S 2c(α) (X) = X -X 2c(α)-1 ∆X X ∈ X 2c(α)-1 , X 2c(α)
and S iα or S i,α correspond to the shape function of node i evaluated at the position of the αth material point.

Scheme equation: Euler time discretization

The method followed in order to write the scheme equation consists in tracing backward the numerical procedure described in section 3.2 in order to get an expression of the form (3.27) for the material point α:

Qn+1 α = H Qn β β = 1, .
., N p Quantities at time t n+1 are obtained by interpolating nodal solutions of the discrete equation (3.22) in the cell containing the αth particle:

Qn+1 α = S 2c(α)-1,α qn+1 2c(α)-1 + S 2c(α),α qn+1 2c(α) (3.33)
The linear scalar advection equation leads to the interface flux function F N = (aq * )N , where q * is the stationary solution of Riemann's problem at a cell interface and N = ±1 the outward unit normal. The discrete form (3.22) thus reads for both nodes of cell c(α):

qn+1 2c(α)-1 = qn 2c(α)-1 + ∆t M L 2c(α)-1 K 2c(α)-1,j aq n j -aρq * 2c(α)-1 N 2c(α)-1 qn+1 2c(α) = qn 2c(α) + ∆t M L 2c(α) K 2c(α),j aq n j -aρq * 2c(α) N 2c(α) (3.34)
where in a one-dimensional grid, the outward unit vectors are N 2c(α)-1 = -1 and N 2c(α) = 1. Discrete equations (3.34) are then simplified by first considering that the mass density is defined in a cell as the ratio of total nodal mass and the element volume:

ρ(X) = M L 2c-1 + M L 2c ∆X = N c p µ=1 m µ ∆X , X ∈ [X 2c-1 , X 2c ]
where N c p is the number of particles in cell c and m µ , the mass carried by the µth material point. For homogeneous medium, we further assume that the mass is uniformly distributed between particles so that the previous definition reduces to ρ = N c p m c /∆X, with m c the mass carried by particles lying in c. Second, linear shape functions lead to the following the lumped mass matrix and pseudo-stiffness matrices:

M L i = Np µ=1 S iµ m µ = m c(i) Np µ=1 S iµ K 2c(α)-1,j = Np µ=1 ∂S 2c(α)-1,µ ∂X m µ S jµ = -m c(i) Np µ=1 S jµ ∆X K 2c(α),j = Np µ=1 ∂S 2c(α),µ ∂X m µ S jµ = m c(i) Np µ=1

S jµ ∆X

The discontinuous approximation basis moreover yields a bloc diagonal pseudo-stiffness matrix so that one can write:

K ij qn j = K i,2c(i)-1 qn 2c(i)-1 + K i,2c(i) qn 2c(i)
Third, a right-going wave (i.e. a > 0) leads to a stationary solution of the Riemann problem equal to the state of the upwind node of an interface, that is:

q * 2c(α)-1 = ρq n 2c(α)-2 = N c(α) p m c(α) ∆X qn 2c(α)-2 q * 2c(α) = ρq n 2c(α) = N c(α) p m c(α) ∆X qn 2c(α)
Therefore, gathering all the previous considerations, equations (3.34) read:

qn+1 2c(α)-1 = qn 2c(α)-1 - a∆t ∆X    N c(α) p µ=1 S 2c(α)-1,µ qn 2c(α)-1 + S 2c(α),µ qn 2c(α) -N c(α) p qn 2c(α)-2 N c(α) p µ=1 S 2c(α)-1,µ    qn+1 2c(α) = qn 2c(α) + a∆t ∆X    N c(α) p µ=1 S 2c(α)-1,µ qn 2c(α)-1 + S 2c(α),µ qn 2c(α) -N c(α) p qn 2c(α) N c(α) p µ=1 S 2c(α),µ    (3.35)
where the volume flux contributions

Qn+1 α = S 2c(α)-1,α q n 2c(α)-1 + S 2c(α),α 1 - a∆t ∆X N c(α) p µ S 2c(α),µ q n 2c(α) + N c(α) p a∆t ∆X S 2c(α)-1,α µ S 2c(α)-1,µ q n 2c(α)-2 + a∆t ∆X S 2c(α),α µ S 2c(α),µ - S 2c(α)-1,α µ S 2c(α)-1,µ f n c(α) (3.36)
Then, the solution at nodes and time step n in equation (3.36) result from the projection between particles and the grid (3.31):

qn i = β S iβ m β Qn β γ S iγ m γ = β S iβ qn β β S iβ (3.37)
In particular, volume flux contributions can be written:

f n c = N c p µ=1 S 2c-1,µ β S 2c-1,β qn β γ S 2c-1,γ + S 2c,µ β S 2c,β qn β γ S 2c,γ = Np β=1 (S 2c-1,β + S 2c,β ) qn β (3.38)
Thus, introduction of mappings (3.37) and (3.38) in equation (3.36) and permutation of sums over β and i lead after some simplifications to the scheme equation:

Qn+1 α = Np β=1 Qn β S 2c(α)-1,α S 2c(α)-1,β µ S 2c(α)-1,µ + S 2c(α),α S 2c(α),β µ S 2c(α),µ - a∆t ∆X N c(α) p S 2c(α),α µ S 2c(α),µ S 2c(α),β µ S 2c(α),µ + a∆t ∆X N c(α) p S 2c(α)-1,α µ S 2c(α)-1,µ S 2c(α)-2,β µ S 2c(α)-2,µ + a∆t ∆X S 2c(α),α µ S 2c(α),µ - S 2c(α)-1,α µ S 2c(α)-1,µ S 2c(α)-1,β + S 2c(α),β (3.39)
Note that the last term of formula (3.39) is non-zero if particles β and α share the same cell, and in that case the parenthesis is one. Hence, the scheme equation can be rewritten as:

Qn+1 α = Np β=1 Qn β      2E i=1 S iβ S iα µ S iµ + N c(α) p a∆t ∆X    S 2c(α)-1,α µ S 2c(α)-1,µ S 2c(α)-2,β µ S 2c(α)-2,µ - S 2c(α),α S 2c(α),β µ S 2c(α),µ 2    + a∆t ∆X S 2c(β),α µ S 2c(β),µ - S 2c(β)-1,α µ S 2c(β)-1,µ (3.40)
The first (resp. second) bracket in equation (3.40) involves shape functions that are non zero if material points β and α lie in adjacent cells (resp. the same cell). Hence, the numerical domain of dependence of the DGMPM for the scalar linear advection equation covers two cells regardless of the number of material points. It is worth noticing that in the particular case of one material point lying in every cell, the convective phase (3.37) simplifies to: qn i = Qn α for c(i) = c(α). The mapping therefore amounts to a piece-wise constant reconstruction of the field on the background grid that makes the DGMPM equivalent to the First Order Upwind (FOU) method. However, this no longer holds for other distributions of material points within the computational grid.

Scheme equation: RK2 time discretization

The discrete system on the grid resulting from the second-order Runge-Kutta time integration consists of the two-stage procedure (3.23) which, for the one-dimensional scalar linear advection equation, particularizes at node i as:

qn+1/2 i = qn i + 1 2 ∆t M L i a   2E j=1 K i,j qn j -q * ,n i N i   (no sum on i) (3.41a) qn+1 i = qn i + ∆t M L i a   2E j=1 K i,j qn+1/2 j -q * ,n+1/2 i N i   (no sum on i) (3.41b)
This procedure can be seen as a recursive use of the Euler scheme (3.35) with suitable time step sizes. The first stage (3.41a) yields the intermediate nodal fields in cell c(α):

qn+1/2 2c(α)-1 = qn 2c(α)-1 - a∆t 2∆X   f n c(α) -N c(α) p qn 2c(α)-2 N c(α) p µ=1 S 2c(α)-1,µ   qn+1/2 2c(α) = qn 2c(α) + a∆t 2∆X   f n c(α) -N c(α) p qn 2c(α) N c(α) p µ=1 S 2c(α),µ   (3.42)
and the second stage (3.41b) leads to the expression of nodal quantities at the end of the time step:

qn+1 2c(α)-1 = qn 2c(α)-1 - a∆t ∆X    N c(α) p γ=1 S 2c(α)-1,γ qn+1/2 2c(α)-1 + S 2c(α),γ qn+1/2 2c(α) -N c(α) p qn+1/2 2c(α)-2 N c(α) p µ=1 S 2c(α)-1,µ    qn+1 2c(α) = qn 2c(α) + a∆t ∆X    N c(α) p γ=1 S 2c(α)-1,γ qn+1/2 2c(α)-1 + S 2c(α),γ qn+1/2 2c(α) -N c(α) p qn+1/2 2c(α) N c(α) p µ=1 S 2c(α),µ    (3.43)
Then, introduction of the interpolation from nodes to particles (3.33) in equations (3.43) leads to the solution at material point α and time step n + 1:

Qn+1 α =S 2c(α)-1,α qn 2c(α)-1 - a∆t ∆X S 2c(α)-1,α -S 2c(α),α γ S 2c(α)-1,γ µ S 2c(α),µ qn+1/2 2c(α)-1 + S 2c(α),α qn 2c(α) + a∆t ∆X S 2c(α),α -S 2c(α)-1,α γ S 2c(α),γ µ S 2c(α)-1,µ -N c(α) p S 2c(α),α µ S 2c(α),µ qn+1/2 2c(α) + N c(α) p a∆t ∆X S 2c(α)-1,α µ S 2c(α)-1,µ qn+1/2 2c(α)-2
Nodal values q n+1/2 i are provided by the first stage of RK2 algorithm and can be substituted in the second stage:

Qn+1 α = S 2c(α)-1,α qn 2c(α)-1 + S 2c(α),α qn 2c(α) - a∆t ∆X S 2c(α)-1,α -S 2c(α),α µ S 2c(α)-1,µ µ S 2c(α),µ   qn 2c(α)-1 - a∆t 2∆X   f n c(α) -N c(α) p qn 2c(α)-2 µ S 2c(α)-1,µ     + a∆t ∆X S 2c(α),α 1 - N c(α) p µ S 2c(α),µ -S 2c(α)-1,α µ S 2c(α),µ µ S 2c(α)-1,µ   qn 2c(α) + a∆t 2∆X   f n c(α) -N c(α) p qn 2c(α) µ S 2c(α),µ     +N c(α) p a∆t ∆X S 2c(α)-1,α µ S 2c(α)-1,µ   qn 2c(α)-2 + a∆t 2∆X   f n c(α)-1 -N c(α) p qn 2c(α)-2 µ S 2c(α)-2,µ     (3.44)
Note that the solution of the downstream node of the adjacent cell q n+1/2 2c(α)-2 results from the second equation of the set (3.42). Therefore, by rearranging formula (3.44) as:

Qn+1 α = S 2c(α)-1,α - a∆t ∆X S 2c(α)-1,α -S 2c(α),α µ S 2c(α)-1,µ µ S 2c(α),µ qn 2c(α)-1 + S 2c(α),α + a∆t ∆X S 2c(α),α 1 - N c(α) p µ S 2c(α),µ -S 2c(α)-1,α µ S 2c(α),µ µ S 2c(α)-1,µ qn 2c(α) + 1 2 a∆t ∆X 2   N c(α) p S 2c(α)-1,α µ S 2c(α)-1,µ - S 2c(α),α µ S 2c(α),µ + S 2c(α),α N c(α) p µ S 2c(α),µ 2   qn 2c(α) +N c(α) p a∆t ∆X S 2c(α)-1,α µ S 2c(α)-1,µ 1 - a∆t 2∆X 1 + N c(α) p µ S 2c(α)-2,µ + a∆t 2∆X S 2c(α),α µ S 2c(α),µ qn 2c(α)-2 - 1 2 a∆t ∆X 2 N c(α) p S 2c(α),α µ S 2c(α),µ 2 f n c(α) + 1 2 a∆t ∆X 2 N c(α) p S 2c(α)-1,α µ S 2c(α)-1,µ f n c(α)-1 µ S 2c(α)-2,µ
the use of mapping equations (3.37) and (3.38) allows to write:

Qn+1 α = β Qn β S 2c(α)-1,β µ S 2c(α)-1,µ S 2c(α)-1,α - a∆t ∆X S 2c(α)-1,α -S 2c(α),α γ S 2c(α)-1,γ µ S 2c(α),µ + S 2c(α),β µ S 2c(α),µ S 2c(α),α + a∆t ∆X S 2c(α),α 1 - N c(α) p µ S 2c(α),µ -S 2c(α)-1,α γ S 2c(α),γ µ S 2c(α)-1,µ + S 2c(α),β µ S 2c(α),µ 1 2 a∆t ∆X 2   N c(α) p S 2c(α)-1,α µ S 2c(α)-1,µ - S 2c(α),α µ S 2c(α),µ + S 2c(α),α N c(α) p µ S 2c(α),µ 2   + N c(α) p S 2c(α)-2,β µ S 2c(α)-2,µ a∆t ∆X S 2c(α)-1,α µ S 2c(α)-1,µ 1 - a∆t 2∆X 1 + N c(α) p µ S 2c(α)-2,µ + a∆t 2∆X S 2c(α),α µ S 2c(α),µ + 1 2 a∆t ∆X 2 N c(α) p    S 2c(α)-2,β + S 2c(α)-3,β µ S 2c(α)-1,µ µ S 2c(α)-2,µ S 2c(α)-1,α - S 2c(α)-1,β + S 2c(α),β µ S 2c(α),µ 2 S 2c(α),α        
Once the previous formula is simplified, the one-dimensional scheme equation of the DGMPM with the RK2 time discretization can be written:

Qn+1 α = β Qn β i S iβ S iα µ S iµ + a∆t ∆X S 2c(β),α µ S 2c(β),µ - S 2c(β)-1,α µ S 2c(β)-1,µ + a∆t ∆X N c(α) p    S 2c(α)-1,α µ S 2c(α)-1,µ S 2c(α)-2,β µ S 2c(α)-2,µ - S 2c(α),α S 2c(α),β µ S 2c(α),µ 2    + 1 2 a∆t ∆X 2 N c(α) p S 2c(α),β µ S 2c(α),µ - S 2c(α)-2,β µ S 2c(α)-2,µ S 2c(α)-1,α µ S 2c(α)-1,µ - S 2c(α),α µ S 2c(α),µ + 1 2 a∆t ∆X 2 N c(α) p S 2c(β),α µ S 2c(α),µ 2 N c(α) p S 2c(α),β µ S 2c(α),µ -1 + 1 2 a∆t ∆X 2 S 2c(β)+1,α N c(α) p µ S 2c(α)-1,µ µ S 2c(α)-2,µ 1 -N c(α) p S 2c(α)-2,β µ S 2c(α)-2,µ (3.45)
The three first terms of the latter scheme equation correspond to that obtained for the Euler algorithm (3.40) while the second order terms are provided by the two-stage time integration. The brackets in those higher-order corrections vanish when only one point is in each cell of the grid, so that the scheme is, as that of the previous section, equivalent to the FOU method.

The von Neumann linear stability analysis

The scheme equations obtained above are written for simplicity:

Qn+1 α = Np β=1 D αβ Qn β (3.46)
Moreover, the computational domain is repeated periodically by mapping it to the domain [-l, 0] so that the solution at material point α and time step n is expanded into a discrete Fourier basis of 2E + 1 harmonics over the domain X ∈ [-l, l]:

Qn α = E j=-E A n j e iαkj ∆X
with A n j , the magnitude of the jth harmonic at time step n, i = √ -1, and k j a wave number. Introduction of this expansion in equation (3.46) yields:

A n+1 j e iαkj ∆X = Np β=1 A n j D αβ e iβkj ∆X ∀j = -E, ..., E
The amplification factor between two time steps at a given point is defined as:

A n+1 j A n j = Np β=1 e i(β-α)kj ∆X D αβ ∀j = -E, ..., E (3.47)
A necessary condition to ensure the stability of a numerical scheme is that the amplification factor must be lower than or equal to one in modulus: A n+1 /A n ≤ 1. This upper bound prevents an increasing error during the computation. For expression (3.47), this leads to:

Np β=1 e i(β-α)kj ∆X D αβ ≤ Np β=1 e i(β-α)kj ∆X D αβ = Np β=1 |D αβ | ∀j = -E, ..., E
where the triangle inequality, and the unit modulus of the complex number e i(β-α)kj ∆X have been used. Hence, the Courant number must be set so that the following condition is satisfied for all material points:

Np β=1 |D αβ | ≤ 1 ∀ α = 1, ..., N p (3.48)
Note however that the use of the triangle inequality leads to a more severe constraint than the one really holding. As a consequence, the Courant number can be set in practice to higher values than that resulting from the solution of (3.48).

According to scheme equations (3.40) or (3.45), the stability condition (3.48) can be very hard to find analytically for general discretizations. Nevertheless, though the infinity of possible material point distributions prevents the explicit derivation of a general stability condition, the optimal CFL number satisfying the equality in (3.48) can be found numerically. Some configurations are studied in table 3.1 where the critical Courant number resulting from the two time discretizations studied above are compared. Those results have been obtained by using the same particle distribution in every element of a one-dimensional regular mesh. First, the DGMPM scheme is stable for all CFL numbers lower than or equal to one when cells contain one particle only regardless of their positions. This property characterizes the FOU method that is retrieved with this discretization. Then, when many material points share the same element, several configurations are considered:

(i) particles are positioned symmetrically with respect to cell centers and regularly spaced in the mesh. This space discretization is referred to as the natural configuration in the following.

(ii) particles in the natural configuration are all shifted by u = ∆X/10.

(iii) the same as (ii) with u so that one material point overlaps every left node of the cells.

(iv) the same as (iii) for right nodes.

(v) particles are placed symmetrically with respect to cell centers but not regularly spaced in the mesh. Material points in the left half of cells are shifted by u 1 = -∆X/10 while those in the right half are shifted by u 2 = ∆X/10.

(vi) the same as (v) with the first and last particles overlapping left and right nodes of cells respectively. First, table 3.1 shows that the natural configuration leads to a better stability bound for the RK2 integrator while the CFL number allowed by using Euler time discretization decreases with increasing numbers of particles per cell. Second, moving every point rightward in the mesh (i.e. cases (ii) and (iv)) causes a drop in the critical Courant number for both RK2 and Euler algorithms. In particular, stability conditions for more than two particles per cell provided by the RK2 are, in those cases, more restrictive than that of the Euler. Third, the leftward shift (i.e. case (iii)) leads to an improvement of the stability condition for Euler time integration compared to that in the natural configuration. The CFL number provided by RK2 integration also decreases due to the shift while remaining higher than the Euler. At last, particle distributions conserving the symmetry with respect to cell centers (i.e. cases (v) and (vi)), yield the optimal stability condition for the RK2 integrator while the Euler CFL depends on the spacing between material points. More specifically, the Euler algorithm leads to a vanishing CFL for the case (vi), thus preventing any simulation. The distribution of material points, and therefore the resulting quadrature, has an influence on the stability of the scheme. The optimal Courant number can be reached even with the Euler discretization, whereas the classical DGFEM scheme developed in [START_REF] Chavent | A finite-element method for the 1-d water flooding problem with gravity[END_REF] is restricted to condition ∆t/∆X = O( √ ∆X). This limitation has been addressed by introducing slope limiters in order to remove non-physical oscillations in the vicinity of sharp solutions while providing high-order accuracy in smooth regions [START_REF] Chavent | The local projection P 0 -P 1 -discontinuous-galerkin finite element method for scalar conservation laws[END_REF]. However, the stability of the method is still bounded by CF L ≤ 1/2 and the scheme is first-order accurate. The use of a second-order Runge-Kutta [START_REF] Cockburn | The runge-kutta local projection P 1 -discontinuous-galerkin finite element method for scalar conservation laws[END_REF] allows second-order accuracy of the scheme, but the stability condition then reduced to CF L ≤ 1/3. It is worth noting that space-time DGFEM formulations [START_REF] Richter | An explicit finite element method for the wave equation[END_REF][START_REF] Abedi | A space-time discontinuous galerkin method for linearized elastodynamics with element-wise momentum balance[END_REF] provided more recently the ability to relax constraints of pure space DGFEM and obtain a critical CFL of 1.

Number of particles Position of particles in cell c

Space convergence analysis of the one-dimensional scheme

We now propose to focus numerically on convergence properties of DGMPM schemes for linear elasticity problems in homogeneous solids undergoing one-dimensional stress and strain states so that the bar assumption under small strains holds. The following system of one-dimensional linear scalar advection equations is then considered:

∂ ∂t v σ - ∂ ∂X 1 ρ σ Ev = 0 0 (3.49)
in a bar of length l = 1 m, with Young's modulus E = 2 × 10 11 P a, and mass density ρ = 7800 kg • m -3 . The bar is assumed to be initially motionless in a free-stress state, that is: v(X, 0) = 0 ; σ(X, 0) = 0 ∀X ∈ [0, l], and subject to the following boundary conditions on its left and right ends respectively:

σ(X = 0, t) = σ sin 0.4πt T σ(X = l, t) = 0, ∀t ∈ [0, T ] (3.50)
where T = l/c is the time taken by waves to travel from one end of the bar to the other at the sound speed c = E ρ . The combination of the hyperbolic system (3.49) and boundary conditions (3.50) yields a Picard problem whose exact solution is [36, Ch.2]:

σ(X, t) = σ sin 0.4π(t -X/c) T v(X, t) = - σ ρc sin 0.4π(t -X/c) T , ∀ X, t ∈ [0, l] × [0, T ] (3.51)
The continuum system (3.49) is discretized with both MPM (the USF formulation is selected to limit the numerical dissipation) and DGMPM by using 2 i cells (i = 2, 3, ..., 7) each containing N p = 2, 3, 4, 8 regularly spaced material points so that the regular grid size is ∆X = l/2 i . The refinement of the mesh operates on both nodes and particles in such a way that the number and position of material points per cell initially selected is held constant. Configurations involving different numbers of particles per element are studied separately and the CFL number can therefore be set once and for all (optimal for DGMPM, and 0.5 for the MPM). Numerical solutions are then compared to the exact one (3.51) according to the relative L 2 norm of the error, also referred to as the L 2 error :

σ (t) =   ∆X Np α=1 (σ α -σ exact (X α , t)) 2 σ exact (X α , t) 2   1 2 v (t) =   ∆X Np α=1 (v α -v exact (X α , t)) 2 v exact (X α , t) 2   1 2
The errors are computed at time t = 0.5T so that no reflection of waves on the right boundary occur and the stress amplitude is set at σ = 4 × 10 4 P a. Figure 3.10 shows the evolution of the L 2 errors σ with regard to the grid size for the MPM and the DGMPM using Euler and RK2 time discretizations for several numbers of particles per element. We see that, despite the rates of convergence of stress are close to one (obtained by means of non-linear least squares fitting and reported in table 3.2), the DGMPM-Euler and DGMPM-RK2 exhibit lower rates compared to those of the MPM. Second, an increase in the number of particles per cell provides slightly higher convergence rates, according to table 3.2, and shifts DGMPM curves upward as can be seen in figures 3.10a to 3.10d. On the other hand, the MPM error curves do not exhibit this behavior as observed in [START_REF] Tielen | High-order material point method[END_REF]. At last, for coarse grids (i.e. ∆X = l/2 i for i ≤ 3) the DGMPM-Euler leads to the lowest error but its smaller rate of convergence makes it quickly less accurate than the two other methods. Similar behavior can be seen for the velocity, whose convergence curves are depicted in figure 3.11 for the same discretizations. Again, the DGMPM shows lower rates of convergence than the MPM, and increasing the number of particles per cell leads to more error in DGMPM solutions. The same order of accuracy for velocity and stress resulting from the DGMPM (see table 3.2) was expected since the weak form of system (3.49) leads to the same approximation for both fields. The MPM is on the other hand, as FEM, characterized by a higher order of convergence for The loss of accuracy for velocity in DGMPM can be explained by the mapping procedure it uses from nodes to particles since the FLIP mapping has been introduced in order to reach second-order accuracy [START_REF] Nishiguchi | Second-order fluid particle scheme[END_REF]. Convergence properties of the MPM using the PIC projection shown in table 3.2 and figure 3.11 confirm the previous remark. Indeed, modifying the way fields are projected from nodes to particles does not affect dramatically the rate of convergence of stress but leads to a significant decrease in velocity accuracy. We now propose to set the Courant number to 0.1 for the MPM and the DGMPM-Euler in order to compare the accuracy of those first order in time schemes with the same time steps. A comparison between the convergence curves obtained for a CFL set to 0.1 and the previous results (i.e. figures 3.10 and 3.11) is made in figure 3.12. The results first show that the order of accuracy of the DGMPM-Euler is also lower than that of the MPM when the same CFL is used in the methods. Second, as can be seen in figure 3.12, the reduction of the time step yields additional error which can be explained by the numerical diffusion introduced in the solutions. As a consequence, the curves resulting from computations performed with a Courant number set to 0.1 are above those obtained previously. However, a lower CFL enables an improvement of the rates of convergence of the methods. Indeed, the rates of convergence for both velocity and stress are in that case 0.66 for the DGMPM-Euler, and 1.00 for the MPM. Since the DGMPM aims at capturing non-regular solutions (discontinuities), its first-order of convergence is sufficient. Indeed, the achievable accuracy for such solutions is at most one [2, p.149]. Nevertheless, one might be interested in solving problems that involve smooth solutions so that the need for higher-order accuracy arises. Although linear shape functions have been employed so far, the method may handle higher-order approximation by using Legendre polynomials or B-Spline functions for instance (higher-order Lagrange polynomials being to be avoided for moving particles [START_REF] Tielen | High-order material point method[END_REF]). Furthermore, DG features enable the local increase in the approximation order in regions where the solution is smooth. Note however that the particle-based quadrature rule limits the shape function order one can select since few material points may lead to reduced integration. This issue can however be avoided by means of function reconstruction techniques which build a (linear, quadratic etc.) function on the grid that is based on the values a field takes at particle locations and which can be evaluated everywhere in a cell. In particular, the use of moving least squares [START_REF] Ming | Improving the Material Point Method[END_REF] or spline interpolation [START_REF] Tielen | High-order material point method[END_REF][START_REF] Tielen | A high order material point method[END_REF] within the MPM, enables evaluation of a reconstructed function at Gauss point locations in order to improve the quadrature rule. However, the introduction of such reconstruction techniques in the DGMPM might lead to a restrictive stability condition as it would make the method closer to DGFEM.

Finally, it is worth noticing that only regularly-spaced material points have been considered in the above convergence analysis. An additional means of achieving high-order accuracy can be derived from the following theorem applying to the one-dimensional scalar linear advection equation with a > 0:

Theorem 4 (Roe [21]) A numerical scheme of the form:

Q n+1 α = ∞ β=-∞ D β Q n α+β is pth order accurate (p ≥ 0)

in space and time if and only if:

β β k D β = -a ∆t ∆X k 0 ≤ k ≤ p
The complexity of scheme equations (3.40) and (3.45) developed previously, and the lack of such equations for MPM formulations prevent the straightforward calculation of convergence rates. However, one can imagine the building of adaptive strategies of the set of particles based on theorem 4, so that a given order of accuracy is reached. If such a distribution of material point exists, the corresponding Courant number can also be computed, thus ensuring both stability and accuracy of the numerical scheme.

Two-dimensional stability analysis

Model equation -Space discretization

We now move on to the scalar linear advection equation for an arbitrary quantity q = ρq moving at constant speeds a, b ∈ R + * in the X and Y directions respectively in a homogeneous two-dimensional medium:

∂ q ∂t + a ∂ q ∂X + b ∂ q ∂Y = 0 (3.53)
Note that the flux functions are in this case f X = aq and f Y = bq. The physical domain [0, l] × [0, h] is discretized with N p material points arbitrarily distributed in a Cartesian grid made of E four-node bilinear elements with constant size ∆X × ∆Y . With the nodal positions in cell C denoted by

X C i = X C i Y C i , as depicted ξ η 1 -1 -1 1 X Y X 1 X 2 X 3 X 4 X ξ
ξ = 2 X -X C 1 ∆X -1 ; dξ = 2 dX ∆X η = 2 Y -Y C 1 ∆Y -1 ; dη = 2 dY ∆Y (3.54)
Horizontal and vertical edges lengths are distinguished here in spite of the Cartesian nature of the grid in order to easily extend the following study to rectilinear grids. Again, there are no empty cells inside the physical domain so that no hole is generated, and periodic boundary conditions are considered to simplify the analysis.

Two-dimensional scheme equation

One can imagine the combination of the DGMPM discretization with a multi-stage time integration, as proposed for one-dimensional problems. The analysis of the DGMPM scheme for two-dimensional problems carried out here however only considers the Euler time discretization. Therefore, the updated solution at material point α is obtained by interpolation of nodal solutions satisfying equation (3.53), discretized as:

qn+1 i = qn i + ∆t M L i K X ij aq n j + K Y ij bq n j -f * i (3.55)
so that:

Qn+1 α = 4E i=1 S iα qn+1 i (3.56)
The lumped mass matrix in equation (3.55) has the same expression as in the one-dimensional case that depends on the shape functions of the four-node bilinear element: M L i = µ m µ S iµ . Making use of parent coordinates (3.54), the pseudo-stiffness matrices read:

K X ij = β ∂S iβ ∂X m β S jβ = 2 ∆X β ∂S iβ ∂ξ m β S jβ K Y ij = β ∂S iβ ∂Y m β S jβ = 2 ∆Y β ∂S iβ ∂η m β S jβ
As for the one-dimensional case, the homogeneous medium yields the same mass for every particle so that, by writing ∂(•) ∂ξ = ∂ ξ (•), one gets:

K X ij M L i = 2 ∆X λ ∂ ξ S iλ S jλ γ S iγ K Y ij M L i = 2 ∆Y λ ∂ η S iλ S jλ γ S iγ
The nodal solutions at time n being given by the projection qC,n

i = β S iβ Qn β γ
Siγ , volume fluxes of the discrete form can be rewritten as:

a K X ij M L i qn j = β Qn β 2 ∆X a λ ∂ ξ S iλ j S jλ S jβ γ S iγ µ S jµ b K Y ij M L i qn j = β Qn β 2 ∆Y b λ ∂ η S iλ j S jλ S jβ γ S iγ µ S jµ (3.57)
Then, the nodal interface flux f * i results from the integration of Godunov fluxes along edges connected to the node, according to the weak form (3.21). Referring to a quantity defined at an interface by means of superscripts in parentheses, the Godunov flux corresponding to equation (3.53) at interface (i) is:

f (i) = c n q (i) U = c n q (i) D f N (q (i) D ) -c n (q (i) D -q (i) U ) A + U/D (3.58)
where c n is the speed in the normal direction to the interface (i.e. b for horizontal and a for vertical edges). Equation (3.58) further involves state vectors q (i) U and q (i) D obtained by averaging nodal values connected to interface (i) on upwind and downwind sides respectively, and the right-going fluctuation A + U/D . The CTU is adopted by subtracting from fluxes (3.58) transverse corrections based on those fluctuations according to equation (3.30):

B + A + U/D = c t c n (q (i) D -q (i)
U ) with c t the speed in the tangent direction to the interface. The final expression of intercell fluxes is hence:

f (i) = c n q (i) U -c t c n ∆t 2∆X (i) (q (i) D -q (i) U ) (3.
B + A + BL/B B + A + BL/L f (3) f (2)
f (1) f (4) ∆X∆Y . Thus, introduction of the particle fields projection yields the following expressions for interface fluxes:

f (1) = b q B,n 3 + q B,n 4 2 -ab ∆t 2∆Y q B,n 1 + q B,n 4 2 - q BL,n 2 + q BL,n 3 2 f (2) = a q C,n 2 + q C,n 3 2 -ab ∆t 2∆X q C,n 1 + q C,n 2 2 - q B,n 3 + q B,n 4 2 f (3) = b q C,n 3 + q C,n 4 2 -ab ∆t 2∆Y q C,n 1 + q C,n 4 2 - q L,n 2 + q L,n 3 2 f (4) = a q L,n 3 + q L,n 4 2 -ab ∆t 2∆X q L,n 1 + q L,n 2 2 - q BL,n 3 + q BL,n
f (1) = Np β=1 Qn β bN C p m C 2∆X∆Y S B 3β γ S B 3γ + S B 4β γ S B 4γ -a ∆t 2∆Y S B 1β γ S B 1γ + S B 4β γ S B 4γ - S BL 2β γ S BL 2γ - S BL 3β γ S BL 3γ f (2) = Np β=1 Qn β aN C p m C 2∆X∆Y S C 2β γ S C 2γ + S C 3β γ S C 3γ -b ∆t 2∆X S C 1β γ S C 1γ + S C 2β γ S C 2γ - S B 3β γ S B 3γ - S B 4β γ S B 4γ f (3) = Np β=1 Qn β bN C p m C 2∆X∆Y S C 3β γ S C 3γ + S C 4β γ S C 4γ -a ∆t 2∆Y S C 1β γ S C 1γ + S C 4β γ S C 4γ - S L 2β γ S L 2γ - S L 3β γ S L 3γ f (4) = Np β=1 Qn β aN C p m C 2∆X∆Y S L 3β γ S L 3γ + S L 4β γ S L 4γ -b ∆t 2∆X S L 1β γ S L 1γ + S L 2β γ S L 2γ - S BL 3β γ S BL 3γ - S BL 4β γ S L 4γ
written for simplicity:

f (i) = Np β Qn β c n N C p m C 2∆X∆Y φ (i) β + φ (i),T β (3.60)
In the last expression, φ (i) is devoted to normal contributions while φ (i),T stands for transverse corrections at interface (i). Numerical fluxes considered above are based on normal vectors oriented in the direction of the stream (see figure 3.14). Nodal interface fluxes on the other hand, as defined in the semi-discrete system:

f * i = Γ S i (X)f * N dΓ
are based on the outgoing flux to an element so that f (1) and f (4) must be counted negatively. The integral for cell C is then: (4) dY which can be computed analytically using parent coordinates (3.54):

f * i = - X C 2 X C 1 S i (X, Y C 1 )f (1) dX + Y C 3 Y C 2 S i (X C 2 , Y )f (2) dY + X C 3 X C 2 S i (X, Y C 3 )f (3) dX - Y C 4 Y C 1 S i (X C 1 , Y )f
f * 1 = - 1 2 ∆Xf (1) + ∆Y f (4) ; f * 2 = - 1 2 ∆Xf (1) -∆Y f (2) f * 3 = 1 2 ∆Xf (3) + ∆Y f (2) ; f * 4 = 1 2 ∆Xf (3) -∆Y f (4)
A condensed way of writing those fluxes is adopted by means of the middle point of edge (j) with coordinates

X (j)
1/2 , at which the shape functions are:

S i (X (j) 1/2 ) =    1 2
if node i belongs to edge (j) 0 otherwise.

In addition, components of the outward normal vector to edges

N (i) X and N (i)
Y allow taking into account different signs of intercell fluxes in the Cartesian grid. One thus writes:

f * i = 1 2 edges j 2S i (X (j) 1/2 ) ∆Y N (j) X + ∆XN (j) Y f (j)
which, combined with equation (3.60) leads to:

f * i = Np β Qn β edges j S i (X (j) 1/2 ) a∆Y N (j) X + b∆XN (j) Y N C p m C 2∆X∆Y φ (i) β + φ (i),T β
These terms are divided by the lumped mass matrix in the discrete form:

f * i M L i = β Qn β µ S iµ edges j=1 1 2 S i (X (j) 1/2 )N C p m C aN (j) X ∆X + bN (j) Y ∆Y φ (j) β + φ (j),T β (3.62)
At last, gathering the mapping of updated nodal quantities to the particles (3.56), expressions of volume fluxes (3.57) and intercell fluxes (3.62), the updated value at material point α contained in cell C reads:

Qn+1 α = Np β=1 Qn β 4E i=1 S iα µ S iµ    S iβ + 2 4E j=1 S jβ γ S jγ Np λ=1 S jλ a ∆t ∆X ∂ ξ S iλ + b ∆t ∆Y ∂ η S iλ - 1 2 edges k=1 S i (X (k) 1/2 )N C p a ∆t ∆X N (j) X + b ∆t ∆Y N (j) Y φ (k) β + φ (k),T β (3.63)
Recall that transverse contributions φ (j),T β depend on ∆t, thus providing second-order corrections in the twodimensional scheme equation (3.63), that can also be rewritten as:

Qn+1 α = Np β=1 Qn β D αβ (3.64)

The von Neumann linear stability analysis

Analogously to the one-dimensional case, the solution at a material point can be expanded into a discrete Fourier basis over the domain [-l, l] × [-h, h]. We consider here a structured distribution of particles made of N p = N X p × N Y p material points so that one can denote the solution at particles by Qαβ , where α and β are the row and column of material point indices. For one arbitrary Fourier mode, one has [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF]Ch.20]:

Qn αβ = A n jq e i(αkj +βkq)∆X
where k j and k q are wave numbers. Then, the amplification factor reads:

A n+1 jq A n jq = N X p γ=1 N Y p µ=1 e i([α-γ]kj +[β-µ]kq)∆X D αβ,γµ with i = √ -1.
The requirement that the absolute value of the amplification factor is lower than or equal to one leads to the following stability condition:

A n+1 jq A n jq = N X p γ=1 N Y p µ=1 e i([α-γ]kj +[β-µ]kq)∆X D αβ,γµ ≤ 1 ⇔ N X p γ=1 N Y p µ=1 |D αβ,γµ | ≤ 1
or more simply:

Np β=1 |D αβ | ≤ 1 ∀α = 1, ..., N p (3.65)
Again, the single particle-per-cell discretization leads to a piece-wise constant reconstruction of the field on the computational grid after the projection from material points to nodes, thus providing the first order upwind method. This method is known to be bounded by the stability requirements [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF]Ch.20 Configurations involving more particles in the computational grid cells are then studied numerically by assuming the same material points distribution in every element. Furthermore, we consider only regular cells ∆Y = ∆X and wave speeds satisfying a ≥ b > 0, so that Courant the number is a∆t/∆X. The scheme equation (3.64) can then be written as a function of the CFL number by means of the speed ratio a/b. Hence, the maximal Courant number satisfying the stability condition (3.65) also depends on the speed ratio. Evolutions of the CFL numbers corresponding to several distributions of particles in a two-dimensional grid are gathered in tables 3.3 and 3.4 for the DGMPM scheme using DCU and CTU methods. The first column of these tables shows the positions of material points inside cells for discretizations based on 2 or 4 particles per element.

The space discretization leading to 2 particles lying in every cell of the mesh is such that within an element, the two material points are both either on the horizontal axis or on the vertical axis of the cell, corresponding respectively to the results reported in the first and second rows of table 3.3. Two situations are then to be distinguished:

• Material points are regularly-spaced within the grid and placed symmetrically two-by-two with respect to cell centers. These distributions are drawn in the first column of table 3.3 by using blue circles to represent particles.

• Material points still satisfy symmetry in cells, but are no longer regularly-spaced in the mesh. In that case, particles are drawn with red crosses. First, the results of table 3.3 show that the CFL number exhibits a non-linear dependence on the speed ratio a/b that asymptotically approaches some value which depends on the particle distribution. Second, we see that a reduction of spacing between particles conserving the symmetry between them with respect to cell centers, as for the one-dimensional case, yields an increase in the critical Courant number for both the DCU and CTU approaches. Third, whether particles lie on the horizontal axis or the vertical axis of cells has a great influence on the critical Courant number one can expect. Hence, the configurations of the second row of table 3.3 yield higher CFL numbers for given speed ratios. It then appears that in order to improve the stability of the scheme, one must use a lower number of material points in the direction of the dominating wave speed than in the perpendicular one. For a Cartesian distribution of particles

Particles in cells

N X p × N Y p this corresponds to N Y p > N X p if a > b, and N X p > N Y p if b > a.
At last, it is worth noticing that the improvement brought by the CTU is much less significant than in the case of one single particle-per-cell discretization for which the Courant number can be set to one according to equations (3.66a) and (3.66b).

We now move on to cases for which grid cells each contain 4 material points, by considering a square shaped distribution of particles in every element whose centers coincide with cell centroids. This pattern can be contracted or simply translated without change of shape as depicted in the first column of table 3.4. Two configurations are gathered in each row of the table and are distinguished by using either blue circles or red crosses for material points. Again, we observe that the increase in CFL number enabled by the use of the CTU approach is less important than in the case of one particle-per-cell. Next, as for the one-dimensional cases and configurations studied in table 3.3, we see that the closer particles are from cell centers, the higher the CFL number. Finally, it can be seen from the two last rows of table 3.4 that the translation of the square of particles inside elements does not have great influence on the evolution of the Courant number with respect to the speed ratio.

Particles in cells

Conclusion

In this chapter, the formulation of the MPM has been recalled and some drawbacks of the method when applied to hyperbolic problems have been emphasized in section 3.1. It has been seen that the projection of the nodal velocity field to material points inherited from FLIP, though it reduces the numerical dissipation, introduces noise in the solution. Hence, an alternative method using both the PIC mapping procedure and Discontinuous Galerkin approximation has been proposed in section 3.2 in order to avoid spurious oscillations. The resulting DGMPM is based on the weak form of a system of conservation laws written element by element in an arbitrary computational grid in which particles can move. Interface fluxes involved in boundary integrals of the weak form result from the solution of Riemann problems at cell interfaces computed thanks to an approximate Riemann solver (see section 3.2.4). This method combines thus the strength of Finite Element and Finite Volume methods.

The numerical analysis of the scheme applied to the solution of one and two-dimensional linear scalar advection equations performed in section 3.3 led to the ability to determine the maximal Courant number ensuring stability for a given discretization. This property allows to fully exploiting the ability of the method to, for instance, rebuild the grid by employing adaptive mesh techniques on the reference configuration in order to accurately track waves in the current one. Indeed, after such a reconstruction, the number and positions of material points in grid cells can change and one must properly adapt the CFL number so that the scheme remains stable. An advantage over the original MPM, for which the critical Courant number cannot be computed, is hence highlighted. It is however worth noticing that the MPM seems less dependent on the particle distribution so that the CFL is usually set to an arbitrary value (0.5 or 0.7 for one and two-dimensional problems). The DGMPM is on the other hand, characterized by a CFL number that can be set at one in particular cases.

In addition to the stability of the method, the convergence properties of the DGMPM have been compared to that of the original MPM on a one-dimensional elastic problem. While the MPM, as FEM, shows a second-order accuracy in velocity and first-order in stress, the DGMPM exhibits a first-order accuracy for both fields. The loss of accuracy has been attributed to the back-mapping used in the DGMPM and strategies to handle higherorder approximations have been proposed. However, the purpose of this work being the accurate capturing of waves that can be non-regular, high-order accuracy goes beyond the scope of this thesis.

In the following, attention is paid to the capturing of waves without oscillations unlike what has been highlighted for the MPM in section 3.1. The object of the next chapter is the illustration of the DGMPM performances with one and two-dimensional simulations.

Chapter 4

Numerical Results

Introduction

The Discontinuous Galerkin Material Point Method is now illustrated by comparing some of its solutions to those of other numerical schemes and to exact solutions. First, problems falling in the small deformation framework are considered in section 4.1. More specifically, DGMPM and MPM solutions of the Riemann problem in an elastic bar are compared in section 4.1.1. The methods are then applied to the solution of problems involving multi-dimensional stress and one-dimensional strain states (plane wave problem) in section 4.1.2, and a plane strain state in sections 4.1.3 and 4.1.4. Comparisons with MPM, FVM, FEM and exact (when existing) solutions are shown for elastic, elastic-viscoplastic and elastoplastic solids. These simulations highlight the ability of the DGMPM to track sharp solutions and the possibility of using dedicated numerical tools to accurately deal with history-dependent constitutive models, which are of particular interest for the applications targeted by the method.

Second, in section 4.2 attention is paid to waves propagating in finite deforming solids, for which history effects are not considered. For that purpose, DGMPM simulations performed on plane wave problems in a hyperelastic Saint-Venant-Kirchhoff medium are compared to MPM and exact solutions in section 4.2.1. Then, a comparison between DGMPM, FEM and MPM solutions of a plane strain state problem in a neo-Hookean solid is proposed in section 4.2.2.

Although several constitutive models are assumed in this chapter, elastic, viscous and plastic properties considered are the same for all materials. Table 4.1 summarizes the values of Young's modulus E, Poisson's ratio ν and reference mass density ρ 0 . In addition, linear isotropic or kinematic hardening of modulus C and tensile yield stress σ y are assumed for plastic evolutions, along with the viscosity γ and sensitivity n for viscous ones. At last, no body forces are considered here. E = 2 × 10 11 P a σ y = 4 × 10 8 P a n = 4.37

ν = 0.3 C = 10 9 P a γ = σ y × τ 1/n ρ 0 = 7800 kg.m -3

Linearized geometrical framework

Riemann problem in an isotropic elastic bar

To begin with, let's focus on the problem that illustrated some shortcomings of the MPM in section 3.1 and motivated the development of the DGMPM. We thus consider a bar of length l = 6 m in direction e 1 in which the Cauchy stress and infinitesimal strain tensors are of the form:

σ = σ e 1 ⊗ e 1 ε = ε e 1 ⊗ e 1
The stress is initially zero everywhere and Riemann-type initial conditions on the axial velocity

v = v • e 1 are prescribed: v = v 0 > 0 for x 1 ∈ [0, L/2[ and v = -v 0 for x 1 ∈ ]L/2, L].
In addition, both ends of the domain are traction free. The exact solution of this problem [36, Ch.1], recalled in section 2.4.2, consists of two elastic discontinuities propagating leftward and rightward in the bar at constant speeds c = ± E/ρ. The discretization of the domain lies on a regular background grid made of 50 cells containing material points distributed so that two situations are distinguished: each cell contains one particle that coincides with the element centroid for the 1ppc discretization; each cell contains two particles symmetrically placed with respect to element centers and regularly spaced in the grid for the 2ppc discretization.

The problem is solved on the one hand with the MPM-USL in which the nodal velocity is based on either FLIP or PIC mappings (CFL=0.5), and with the DGMPM coupled with both Euler (CFL=1 with 1 ppc and CFL=0.5 for 2ppc) and RK2 (2ppc only and CFL=1) time integration on the other hand. the numerical solutions at two different times in terms of stress and velocity, compared to exact ones. First, since the Courant number can be set to one for the DGMPM-Euler with 1ppc, the method is able to capture the discontinuities and yields solutions fitting perfectly the analytical ones. The same property holds for the DGMPM-RK2 with 2ppc while for the same discretization the DGMPM-Euler is restricted by a lower CFL number that prevents the accurate resolution of waves. As identified in section 3.1, the use of the PIC mapping within the DGMPM leads to solutions that do not exhibit oscillations. In addition, this projection of updated fields from nodes to particles eliminates the locking of velocity in the central region that can be seen in USL-FLIP solutions. Moreover, the introduction of the DG approximation within the USL-PIC leads to a reduction of numerical diffusion, though less significant than that permitted by using FLIP mapping as originally proposed in the MPM. This can be seen in figure 4.2 in which the evolution of total energies resulting from every numerical scheme is depicted. The situations for which the CFL number is set to unity for DGMPM formulations obviously yields an exact conservation of the total energy during the computation while other results suffer from dissipation.

Since, the number of material points per cell has little influence on USL results, the MPM is from now only used with the 1ppc discretization. Furthermore, the PIC mapping has been used here within the MPM for comparison purposes and is no longer considered in the remainder of the simulations.

Plane wave in a history-dependent material

We now consider a infinite medium in directions e 2 and e 3 , and length l = 6 m in direction e 1 . Riemann-type initial conditions similar to those treated above are assumed to yield the following infinitesimal strain and Cauchy stress tensors:

ε = εe 1 ⊗ e 1 σ = σ L e 1 ⊗ e 1 + σ T (e 2 ⊗ e 2 + e 3 ⊗ e 3 )
which correspond to the plane wave case. In that configuration, a relation exists between longitudinal and transverse stress components σ L and σ T in such a way that a one-dimensional hyperbolic system is solved for σ L = σ, and the transverse component σ T is computed subsequently. In this section, the behavior of the DGMPM on relaxation systems is looked at on a solid made of an elastic-viscoplastic material following the Perzyna model with linear kinematic hardening [START_REF] Perzyna | Fundamental problems in viscoplasticity[END_REF]. In the asymptotic limit τ = (γ/σ y ) n → 0, where τ is the relaxation time, the computed elastic-viscoplastic solution should tend to the elastoplastic one derived in [START_REF] Heuzé | Lax-Wendroff and TVD finite volume methods for unidimensional thermomechanical numerical simulations of impacts on elastic-plastic solids[END_REF]. The writing of the viscosity as a function of the relaxation parameter in table 4.1 enables the tuning of the stiffness of the hyperbolic system by setting different values of τ .

The solid is initially in a free stress state and the initial velocity is set so that plastic flow occurs:

v 0 = 2 Y H ρc L
where Y H = (λ + 2µ)σ y /2µ denotes the Hugoniot elastic limit, c L = (λ + 2µ)/ρ is the elastic pressure wave speed, and (λ, µ) are Lamé's constants. Both ends of the medium are traction free so that rightward and leftward compressive elastic waves reflect as unloading waves that interact with the incident plastic ones [START_REF] Heuzé | Simulation of impacts on elastic-viscoplastic solids with the flux-difference splitting finite volume method applied to non-uniform quadrilateral meshes[END_REF].

Elastoviscoplasticity

The elastic-viscoplastic problem is solved with the MPM using both USL and USF formulations, the DGMPM-Euler with Godunov splitting, and the DGMPM-RK2 coupled to Strang splitting. The latter formulation is however not used for stiff systems since this fractional method is known to fail to assess the correct solution in those cases [START_REF] Heuzé | Simulation of impacts on elastic-viscoplastic solids with the flux-difference splitting finite volume method applied to non-uniform quadrilateral meshes[END_REF][START_REF] Leveque | A study of numerical methods for hyperbolic conservation laws with stiff source terms[END_REF]. The ODE systems resulting from fractional approaches are discretized with an implicit backward Euler scheme for Godunov and a backward differentiation formula of order 3 for Strang splitting. The viscoplastic flow rule is then integrated explicitly at the end of the time step to update viscoplastic strains. On the other hand, constitutive equations are integrated with a radial return algorithm [START_REF] Simo | Computational inelasticity[END_REF] in the MPM. First, the relaxation system is considered in a non-stiff setting characterized by a relaxation time bigger than the time step governed by the convection part, that is τ = 50∆t. Figure 4.3 shows a comparison of numerical stress and plastic strain with the exact solutions of the elastoplastic limit. For this non-stiff configuration, viscous effects lead to much smoother solutions compared to the elastic-plastic one as can be seen in figures 4.3a. First, USL and USF results are quite similar in terms of stress, up to some oscillations on the plastic plateau that lead to different assessments of the viscoplastic strain. Then, Godunov splitting with 1ppc and Strang splitting provide solutions in which one can distinguish elastic waves from the viscoplastic flow while it is not possible in the MPM solution and in DGMPM-Euler results with 2ppc due to numerical diffusion. Furthermore, local overshoots appear in the viscoplastic strain computed with the MPM and the DGMPM-Euler when one ppc is used (see figure 4.3a). The former can be explained by the velocity locking at the middle of the bar and the latter can be eliminated by integrating implicitly the viscoplastic flow rule together with the source term [START_REF] Heuzé | Simulation of impacts on elastic-viscoplastic solids with the flux-difference splitting finite volume method applied to non-uniform quadrilateral meshes[END_REF]. The same remarks as before can be made with a lower relaxation time τ = ∆t × 10 -2 as can be seen in figure 4.4. However, numerical solutions are sharper and get closer to the elastic-plastic stress and strain. Furthermore, the decrease in the relaxation time leads, as expected, to an apparent yield stress tending to that of the elastoplastic solution for Godunov splitting with 1ppc and Strang splitting. Once again, the overshoots arising in the DGMPM-Euler viscoplastic strain can be removed by integrating the flow rule implicitly. This is also the case for the spurious oscillations that can be observed in the stress solution in figure 4.4a [START_REF] Heuzé | Simulation of impacts on elastic-viscoplastic solids with the flux-difference splitting finite volume method applied to non-uniform quadrilateral meshes[END_REF].

As for FVM, the DGMPM can benefit from splitting methods to compute source terms. Though these methods are convenient and easy to implement, some of them can lead to non-physical solutions in the stiff limit in such a way that the solution of the source term pollutes that of the convective part (i.e. as shown for Strang splitting [START_REF] Heuzé | Simulation of impacts on elastic-viscoplastic solids with the flux-difference splitting finite volume method applied to non-uniform quadrilateral meshes[END_REF]). The Godunov splitting is robust and leads to acceptable results but is only first order accurate. Higher-order of accuracy may require the employment of high-order time integrators for stiff systems, such as implicit-explicit IMEX Runge-Kutta schemes [START_REF] Pareschi | Implicit-explicit runge-kutta schemes and applications to hyperbolic systems with relaxation[END_REF]. These questions do not arise for MPM since classical elastic-viscoplastic integrators are used (which may still have convergence troubles in the stiff limit).

Rather than solving the relaxation system in the stiff limit, the solution of the elastic-plastic problem can be computed, for one-dimensional solids, based on dedicated approximate Riemann solvers. Such approaches allow better accuracy while avoiding numerical difficulties related to the computation of stiff source terms. An elastoplastic Riemann solver is used for the simulations presented below.

Elastoplasticity

The conservative and associated quasi-linear forms governing elastoplasticity under small strains derived in section 2.2.4 are used within the DGMPM to write the discrete equations and to compute intercell fluxes respectively. Those two systems are recalled here for convenience:

U t + D i=1 ∂F • e i ∂x i = 0 (4.1a) Q t + A i ∂Q ∂x i = 0 with U = v ε , F • e i = -1 ρ σ • e i -v⊗ei+ei⊗v 2 and Q = v σ
So far, the solution of Riemann problems at cell faces has been carried out by means of a Riemann solver that only takes into account the characteristic structure of the elastic problem. Such a solver may also be used to solve hyperbolic problems in elastic-plastic solids by considering only elastic characteristics for the computation of numerical fluxes. The vector of conserved quantities U resulting from the solution of discrete equations is then used to integrate the plastic flow at the particle level by means of a radial return algorithm. On the other hand, the Riemann solver can be based on exact solutions for linear hardening one-dimensional media [START_REF] Heuzé | Lax-Wendroff and TVD finite volume methods for unidimensional thermomechanical numerical simulations of impacts on elastic-plastic solids[END_REF][START_REF] Wang | Foundations of stress waves[END_REF] by following a prediction-correction procedure: (i) the stationary state is computed based on elastic characteristics only; (ii) the characteristic structure is corrected by adding plastic waves in regions where the yield criterion is violated as described in section 2.4.3. Approximate elastoplastic Riemann solvers, in which non-linear waves are solved as discontinuous ones, can also be employed for non-linear hardening onedimensional media based on the same solution scheme. This type of solver, by accounting for both elastic and plastic waves, should yield more accurate results but does not avoid the integration of plastic flow by radial return algorithms if a conservative form (4.1a) involving ε is used.

The solutions of the DGMPM-Euler scheme, combined with the two approaches discussed above, are compared to MPM and exact solutions on the 1ppc space discretization in figure 4.5. It can be seen in figure 4.5a that the use of an elastic Riemann solver in DGMPM yields slight oscillations on the plastic plateau before the reflection of waves. The same behavior is noticed for the MPM solution. Furthermore, the plastic strain is properly assessed by the DGMPM combined with the elastic-plastic solver while the field is overestimated, in absolute value, by the MPM and the DGMPM using the elastic Riemann solver. Next, additional noise is introduced in numerical solutions by the elastic unloading waves except if the elastic-plastic Riemann solver is employed (see figures 4.5b). In the latter case the solution does not exhibit spurious oscillations. It is however noteworthy that the discontinuity generated in the plastic strain is rather well solved by the DGMPM though all the numerical schemes fail to capture the elastic unloading waves.

We now propose to compare the DGMPM solution based on the elastoplastic Riemann solver to FEM and FVM solutions. For that purpose, finite elements and finite volumes coincide with DGMPM grid cells (a) t = 4.17 × 10 -4 s.

•10 9 (b) t = 9.38 × 10 -4 s.

0 1 2 3 4 5 6 -3 -2.5 -2 -1.5 -1
-0.5

0 •10 -3 x(m) ε p 0 1 2 3 4 5 6
•10 -3

x(m) usl 1ppc dgmpm (elastic-plastic Riemann solver) dgmpm (elastic Riemann solver) exact (b) t = 9.38 × 10 -4 s. for which only the 1ppc discretization is considered. The comparison is therefore made by extracting finite element solutions at integration points, consistently with finite volumes centroids and particles. Figure 4.6 shows the solution of explicit P1-finite elements based on a lumped mass matrix and a radial return algorithm, and a second-order TVD finite volume method using Superbee flux limiters (SB) and the forward Euler time discretization [START_REF] Heuzé | Lax-Wendroff and TVD finite volume methods for unidimensional thermomechanical numerical simulations of impacts on elastic-plastic solids[END_REF]. The introduction of wave limiters in numerical schemes allows to avoids spurious oscillations while steepening the solution near discontinuities by introducing some amount of artificial viscosity based on a Total Variation criterion. The finite volume scheme also makes use of an elastoplastic Riemann solver in such a way that both elastic and plastic waves are limited.

In figure 4.6a, the incident elastic waves are perfectly captured by all methods due to a CFL number set to unity. On the other hand, plastic fronts are steeper for FEM and FVM solutions than for DGMPM, although FEM oscillations yield an overestimated plastic strain. After reflection of waves at the free boundaries (figure 4.6b), additional noise appears in the FEM solution. Those oscillations however occur mainly on elastic waves so that the finite element plastic strain does not exhibit more overshoots. Moreover, the elastic unloading waves are differently solved by FVM and DGMPM. As a result, the plastic strain profiles in figures 4.6a and 4.6b are not identical, the best solution being provided by second-order TVD finite volume. Thus, limiting elastic waves

-1.4 -1.2 -1 -0.8 -0.6 -0.4 -0.2 0 •10 9
σ(P a) (a) t = 4.17 × 10 -4 s.

•10 9 (b) t = 9.38 × 10 -4 s.

0 1 2 3 4 5 6 -2.5 -2 -1.5 -1
-0.5

0 •10 -3 x(m) ε p 0 1 2 3 4 5 6
•10 -3

x(m) dgmpm fem fvm (SB) exact as well as plastic ones, as is allowed by an elastoplastic Riemann solver, enables a significant improvement in the accuracy for plastic waves.

The results presented above show the better accuracy allowed by the introduction of an elastic-plastic Riemann solver. Although the characteristic structure of the solution for one-dimensional media considered here is known, this is not the case in general. Hence, the two-dimensional, three-dimensional or hyperelasticplastic conservation laws are usually solved with finite volume by means of elastic Riemann solvers combined with constitutive update algorithms at cell centroids [START_REF] Haider | A first-order hyperbolic framework for large strain computational solid dynamics: An upwind cell centred total lagrangian scheme[END_REF][START_REF] Sambasivan | A cell-centered lagrangian finite volume approach for computing elasto-plastic response of solids in cylindrical axisymmetric geometries[END_REF][START_REF] Maire | A nominally second-order cell-centered lagrangian scheme for simulating elastic-plastic flows on two-dimensional unstructured grids[END_REF]. Alternatively, the computation of interface fluxes can take into account plastic flow if the Riemann problem is written based on the conservative form (4.1a). Indeed, in that case the flux corresponding to the stationary solution can be calculated through the integration of constitutive equations. Nevertheless, constitutive updates at each face of a grid drastically increase the computational cost. Further research on the characteristic structure of solutions of hyperbolic problems in (hyper)elastic-plastic solids would enable the development of approximate Riemann solvers for more complex constitutive models.

Plane strain problem -Elasticity

We now consider an elastic solid with square section of dimension l = 3m in the (e 1 , e 2 ) plane, infinite in the direction e 3 so that the plane strain assumption holds (i.e. ε 33 = ε 13 = ε 23 = 0). The solid suddenly undergoes a tensile load on a part of its left boundary (see figure 4.7a) leading to shear and pressure waves propagating in the medium until they reflect at the right end. The MPM and the DGMPM are compared to a Q1-finite (a) Geometry and boundary conditions element (bilinear approximation) solution coupled with a central difference explicit time integrator, computed with the code Cast3M [START_REF]User's manual[END_REF]. Since an upper bound of the CFL number cannot be determined for the DGMPM scheme using RK2 time discretization for two-dimensional problems, we consider the DGMPM-Euler scheme in the remainder of the manuscript. More specifically, the following simulations have been carried out by means of the CTU method. The domain is discretized such that material points are equivalent to finite element nodes: l × l ≡ 28 × 28 particles and nodes. Moreover, two arbitrary grids are used for the DGMPM so that either one or four material points lie in every cell according to the situations depicted in figure 4.7b. Analogously to one-dimensional problems, the figures show that DGMPM solutions do not suffer from spurious oscillations while FEM and MPM ones exhibit numerical noise. The 4ppc discretization leads to a decrease in the CFL number which, in turn, yields a less accurate resolution of the jump discontinuity carried by the longitudinal pressure wave than for 1ppc (see figure 4.9a). Next, the interaction of shear and pressure waves traveling in the medium leads to a curved stress profile upstream of the discontinuity that is captured differently by all methods. On the other hand, the cylindrical profile of the longitudinal stress is quite well

l = 3 m a = 1 m σ • e 1 =   σ d 0 0   e 1
0 0.5 1 1.5 2 2.5 3 0 1 2 3 4 5
•10 8

x(m)

σ 11 (P a) (a) t = 2.5 × 10 -4 s 0 0.5 1 1.5 2 2.5 3
•10 8

x(m) described by FEM and DGMPM with 1ppc, even after reflection at the fixed boundary, as can be seen in figure 4.8. The smoothness of DGMPM solutions using 4ppc and the oscillations in MPM results however prevent from distinguishing this structure.

(b) t = 1.0 × 10 -3 s fem ctu 1ppc ctu 4ppc mpm
Figure 4.9b also shows that the stress level on the right boundary of the domain differs from one method to another. Furthermore, MPM and FEM solutions after the passage of the reflected waves oscillate a lot compared to those of DGMPM schemes. In particular, the stress levels greatly differ after the reflection of the pressure wave.

Plane strain problem -Elastoplasticity

The above solid is now assumed to be made of an elastic-plastic material with linear isotropic hardening. A tensile impact of amplitude σ d = 2σ y leading to plastic flow is considered.

Comparison between FEM (Cast3M), MPM and DGMPM using either one or four particles per element is made. The computation of intercell fluxes within the DGMPM is based on an elastic Riemann solver, and plastic flow is integrated by means of a radial return algorithm. The evolution of longitudinal stress σ 11 and plastic strain ε p 11 are depicted in figure 4.10. It can first be seen that similar remarks as above can be made about the longitudinal stress though the integration of plastic flow leads to fewer oscillations in the finite element stress. Next, DGMPM longitudinal plastic strains before the reflection of the pressure wave (figure 4.10a) are quite close to each other. On the other hand, the FEM plastic strain curve is above the others and the MPM one is in advance compared to the other solutions. After the reflection (figure 4.10b), the stress profiles are rather close to each other. However, the MPM solution is far higher than the others on the right end of the domain. The same observations can be made on the longitudinal plastic strain which is quite similar for both DGMPM and FEM solutions, though the DGMPM maximum value computed with one ppc is higher than the others. 4.11 furthermore shows FEM and DGMPM solutions in terms of equivalent plastic strain p, before and after reflection of the longitudinal pressure wave on the right end of the domain. MPM results for the same time steps are plotted separately in figure 4.12 due to the large amplitude of equivalent plastic strain. A concentration of plastic strain occurs at the interface between loaded and traction free parts of the left boundary of the domain in DGMPM solutions. Such concentrations in the high gradient region can also be seen in MPM solutions depicted in figure 4.12 so that this phenomenon seems to be a singularity owed to boundary conditions. Although the final DGMPM and FEM profiles of equivalent plastic strain have the same shape, the results highlight the fact that different numerical methods yield different assessments of irreversible deformations. This point is crucial for the accurate simulation of solid mechanics applications such as forming techniques or crash problems.

Large strain framework 4.2.1 Plane wave in a one-dimensional hyperelastic medium

A plane wave state is considered in a semi-infinite medium of length l = 6m in direction e 1 , made of a hyperelastic Saint-Venant-Kirchhoff material:

F = F e 1 ⊗ e 1 + e 2 ⊗ e 2 + e 3 ⊗ e 3 Π = Πe 1 ⊗ e 1 + Π r (e 2 ⊗ e 2 + e 3 ⊗ e 3 )
where Π = 2µ+λ 2 F (F 2 -1), Π r = λ 2 (F 2 -1) and (µ, λ) are the Lamé's coefficients. A traction force is enforced on the left boundary of the solid initially at rest: Π • (-e 1 ) = Π d e 1 . The exact solutions of the Picard problem thus formulated have been developed in section 2.4.4. With characteristic speeds depending on the deformation gradient, a compressive (resp. tensile) load leads to a rarefaction (resp. shock) wave traveling in the medium. Recall that this non-intuitive behavior is due to the concave nature of the SVK fluxes (see remark 7 in section 2.4.4). Both cases are considered hereinafter before reflection on the right end. Moreover, it has been established that the problem is no longer hyperbolic if the deformation gradient is such that F < 1 3 (see remark 6 in section 2.4.4). Hence, we consider here loading conditions that do not yield a loss of hyperbolicity.

The one-dimensional medium is discretized by using either 100 or 200 material points lying in 100 regular grid cells. The 1ppc and 2ppc discretizations used here are the same as before as well as material parameters of table 4.1.

Compressive impact on a SVK medium

To begin with, the body is submitted to a compressive load on its left end so that a rarefaction wave propagates rightward. The total and updated Lagrangian formulations of the MPM are used along with DGMPM schemes. First of all, the compressive load is set to Π d = 4 × 10 8 P a and numerical solutions are compared in figure 4.13 at two different times. Though the low load amplitude applied results in a solution which looks like a discontinuity, it avoids grid crossing. First, both solutions of updated and total Lagrangian MPM formulations oscillate. It is noteworthy that these two stress profiles are superimposed except at the left end of the domain where an instability occurs in the updated Lagrangian solution due to the displacement of particles. Second, DGMPM solutions still behave well and show good agreement with the exact solution. However, the enforcement of a higher load amplitude is more relevant in order to see the error made by the DGMPM. The applied load is therefore raised to Π d = 2×10 10 P a.

Significant differences between numerical methods and the exact solution are then visible in figure 4.14. First, grid crossing now occurs so that the MPM results based on the updated Lagrangian formulation are no longer correct. Hence, this formulation is omitted from now on. On the other hand, the total Lagrangian MPM scheme yields a solution which does not oscillate more than for linear cases. Next, as expected from the use of an approximate-state Riemann solver for non-linear problems, DGMPM stresses, though close to the analytical solution, no longer fit it exactly (see figure 4.14b). Furthermore, the use of DGMPM-Euler with 2ppc leads to a smoother solution.

One can imagine the employment of an exact Riemann solver based on the exact solution in order to reduce the error made in the DGMPM. However, such an implementation requires the solution of a non-linear problem at each cell interface and is very costly. Moreover, an exact solution is available for a SVK material but this is not the case for other constitutive models, which prevents the generic use of an exact solver.

Tensile impact on a SVK medium

A tensile load of magnitude Π d = -2 × 10 10 P a that gives rise to a shock wave in the medium is now studied.

Recall that for this problem, the characteristic speeds are: c = ± λ+2µ 2ρ0 (3F 2 -1) (see equation (2.73) in section 2.4.4). It then follows that the higher the deformation gradient, the faster waves propagate. Hence, the celerity of the fastest wave is given by the deformation gradient on the upwind side and must be used to adapt the time step and hence, satisfy the CFL condition for fixed Courant number. This situation did not arise in the above case of a rarefaction wave since the highest deformation gradient value was (constant) on the downwind side.

Total Lagrangian MPM and DGMPM are used to solve this problem, and compared to the exact solution in figure 4.15. The error made in the MPM solution grows during the computation, implying that the method is no longer stable in spite of the CFL number set to 0.5, which ensured stability for the above problems. Moreover, slight oscillations also appear after the shock in DGMPM solutions based on a unit Courant number (see figure 4.15a) but not in the DGMPM-Euler using 2ppc due to the lower Courant number used. Nevertheless, the oscillations do not increase with time as can be seen in figure 4.15b, meaning that this noise is not an instability. Next, despite a CFL number set to one, the DGMPM-Euler and DGMPM-RK2 do not perfectly capture the discontinuity, though a good behavior is shown. Indeed, according to Lax entropy condition (2.77), the wave celerity used for the time step evaluation is an upper bound of the shock speed so that such a wave cannot be captured. Nevertheless, the loading conditions are deliberately extreme for the purpose of visualization.

A lower tensile load leads to a numerical solution closer to analytical results as can be seen in figure 4.16, which shows the comparison for Π d = -2 × 10 9 P a. In addition, the oscillations appearing in DGMPM-Euler with 1ppc and DGMPM-RK2 with 2ppc solutions are almost eliminated in figure 4.16a. The MPM solution on the other hand, still exhibits an increasing error although a Courant number ensuring the stability in the scheme for linear problems is used.

Problems in two space dimensions

Plane wave in a two-dimensional hyperelastic medium

We move to multi-dimensional simulations by considering an infinite medium in directions e 2 and e 3 , and of width l in direction e 1 (see figure 4.17), with Riemann-type data on the initial velocity (v 1 (x) = v 0 for x < l/2 and v 1 (x) = -v 0 for x > l/2 with v 0 = 100 m/s). Both ends of the domain are traction-free. The solid is now made of a compressible hyperelastic neo-Hookean material which, unlike the SVK model, is based on a polyconvex stored energy function so that hyperbolicity is ensured whatever F (see section 2.2.3). While it has been seen that a compressive load on a SVK material leads to a rarefaction wave, this is not the case for the neo-Hookean model for which a compressive load involves a shock wave (see remark 7).

h = 2.28 × 10 -1 m l = 6 m e 1 e 2
Figure 4.17: Geometry and loading conditions for the plane wave problem in a twodimensional solid.

The domain is modeled by a finite medium, with zero shear stress and transverse velocity components prescribed on the top, bottom and right boundaries (figure 4.17). In addition zero out-of-plane velocity and strain components are assumed so that plane waves can be simulated. The solid is discretized with 80 × 4 material points in a regular grid made of either 80 × 4 or 40 × 2 cells so that the 1ppc and 4ppc discretizations are used. The simulation has been performed with the DGMPM-DCU, since BCs reproduce a plane wave that does not involve Poisson's effect, and compared to MPM solutions. A similar problem is considered in [START_REF] Renaud | A Discontinuous Galerkin Material Point Method for the solution of impact problems in solid dynamics[END_REF] with a traction force on the left boundary applied on a SVK material, for which a comparison with the exact solution shows a good agreement.

The isovalues of the longitudinal PK1 stress Π 11 can be seen in figure 4.18 at two different times before the wave reaches the left and right ends of the domain. As for one-dimensional problems, oscillations appear in MPM solutions but not in DGMPM ones so that the maximum amplitude of stress in the former is mush higher that in the latter. In contrast to the previous simulations, the lower CFL numbers used for this two-dimensional problem lead to a less accurate resolution of the shock wave with additional numerical diffusion. Nevertheless, the one-dimensional dependency of the solution is not disturbed within the numerical schemes. The evolutions of total energy during DGMPM and MPM computations are compared in figure 4. [START_REF] Burgess | Mass matrix formulation of the flip particle-in-cell method[END_REF]. As for the elastic case, it can be seen that the DGMPM is more dissipative than the MPM. It is however worth noticing that both methods dissipate more energy than in the linear elastic case.

Two-dimensional plane strain problem

The plane strain problem studied in sections 4.1.3 and 4.1.4 is now considered in a compressible hyperelastic neo-Hookean material submitted to an imposed velocity v 1 = -1000 m/s on the bottom part of its left end. The solid is discretized such that material points are equivalent to Q1-finite element nodes. Thus, the plate is represented with l × h ≡ 28 × 28 material points, only with the 1ppc configuration. The finite element computation is performed with the software Abaqus [START_REF] Smith | ABAQUS/Standard User's Manual[END_REF] using an explicit time discretization with no artificial viscosity added. These numerical results are compared to those obtained from MPM and DGMPM using CTU computations. The Courant number is set to unity in DGMPM and to 0.5 in MPM leading to average time steps ∆t CT U = 1.41 × 10 -5 s and ∆t M P M = 6.13 × 10 -6 s, whereas the constant time step used in the FEM simulation is ∆t F EM = 1.27 × 10 -5 s. Figure 4.20 shows numerical results in terms of the Cauchy stress tensor isovalues exported from Abaqus to the software Paraview [START_REF] Ayachit | The paraview guide: a parallel visualization application[END_REF] with the code developed in [START_REF] Liu | Paraview visualization of abaqus output on the mechanical deformation of complex microstructures[END_REF], particularized to the present two-dimensional plane strain case. Cauchy stress is plotted on the current configuration in such a way that figure 4.20 also enables the comparison of the deformed shape of the body. At the beginning of the computation (first row in figure 4.20), stress profiles are quite similar despite slight oscillations visible in FEM and MPM solutions. This can also be seen in figure 4.21, in which stress is plotted along the bottom boundary of the domain. However, the MPM solution exhibits, as for small strain problems, a concentration of stress in the high gradients region on the left boundary. It is worth noticing that the DGMPM shows the same behavior that cannot be seen here due to the MPM stress values which are much higher. The deformed shapes of the plate resulting from the three numerical approaches hence remain close, except at the junction of the loaded and free zones of the left edge. When the pressure wave reflects on the fixed boundary at time t = 5.0 × 10 -4 s (second row in figures 4. 20 and 4.22), the stress profiles are still similar, though FEM and MPM solutions oscillate even more. These spurious oscillations are more significant in the velocity fields depicted in figure 4.22 arbitrary grid (i.e. containing one material point that belongs to the right end of the domain) that is about to be reached by the wave through the upwind interface. The intercell flux on the upwind interface resulting from the discontinuity, and subsequently the conserved quantities vector resulting from the solution of the discrete system on the grid, are non-zero. In particular, the horizontal velocity at upwind nodes of the boundary cell does not vanish while that of the downwind edge satisfies the homogeneous Dirichlet condition. Hence, the interpolation of the velocity from nodes to the particle yields a non-zero field at the material point level. Note that this holds for the MPM as well in which the enforcement of boundary conditions is still a challenging question [START_REF] Steffen | Examination and analysis of implementation choices within the material point method (MPM)[END_REF].

Nevertheless, no significant displacements of particles can be seen on the right end in MPM and DGMPM solutions in figures 4.20 and 4.22. At last, oscillations remain in FEM and MPM solutions until the end of the simulation. Since the velocity field depicted in figures 4.22 and 4.23 is used to update the shape of the solid in FEM, the numerical noise yields final configurations that are slightly different. On the other hand, updating particle positions with the grid velocity within the MPM allows better results than if the oscillating material point velocity is used.

Conclusion

The Discontinuous Galerkin Material Point Method has been applied to hyperbolic problems of solid mechanics. It has first been shown on a Riemann problem in a linear elastic bar under small strains in section 4.1.1, that the method is able to capture the exact solution that consists of two elastic discontinuities propagating in the medium. Indeed, the stability properties of the one-dimensional schemes derived in section 3.3.1 enable, for particular space discretizations, the use of a CFL number set to unity. Nevertheless, once the optimal stability condition is lost, that is CF L < 1, the method is slightly more diffusive than the MPM. Next, the solution of problems in history-dependent solids (sections 4.1.2) have shown that efficient tools can be embedded into the method in order to deal with (visco)plastic flows. In particular, approximate elastic-plastic Riemann solvers can be employed, provided that the characteristic structure of the problem is known. In addition, the results of section 4.2.1 highlight that the total Lagrangian formulation of the DGMPM allows circumventing the eventual grid-crossing occurring in updated Lagrangian MPM for problems involving waves in finite deforming solids. As a consequence, the numerical scheme also provides solutions that are close to exact ones for non-linear problems. Moreover, the arbitrariness of the grid can be fully exploited by employing adaptive mesh techniques on the reference configuration so as to track accurately waves in the current configuration for problems involving complex geometries. At last, the two-dimensional simulations performed in sections 4.1.3, 4.1.4 and 4.2.2 show that DGMPM results are in good agreement with FEM while eliminating spurious oscillations.

Recall that the approach followed in this work consisted of: (i) removing the spurious oscillations that appear in MPM solutions by reintroducing the PIC mapping; (ii) reducing the numerical diffusion thus introduced by means of the Discontinuous Galerkin approximation. The results presented in this chapter showed that these purposes have been fulfilled (see section 4.1.1 for point (ii)). It should be noted, however, that the adaptation of the method to high-order space approximations, as well as the employment of the RK2 time discretization for two-dimensional problems, are expected to further reduce the numerical diffusion exhibited by the scheme.

Moreover, numerical tools used for one-dimensional problems can be generalized for two-dimensional ones. This is the case of splitting procedures that should enable the DGMPM to accurately follow waves traveling in elastic-viscoplastic solids in 2D. As discussed in section 4.1.2, elastoplastic problems can be solved as the stiff limit of elastoviscoplastic ones, provided that suitable ODE integrators are employed. On the other hand, the use of approximate elastic-plastic Riemann solvers avoids the difficulties related to the integration of stiff ODEs, and improves the resolution of plastic waves by means of limiters on both elastic and plastic waves [START_REF] Heuzé | Lax-Wendroff and TVD finite volume methods for unidimensional thermomechanical numerical simulations of impacts on elastic-plastic solids[END_REF]. However, this kind of solver does not exist for general two-dimensional problems owing to the lack of knowledge of the characteristic structure of the solutions for multi-dimensional elastic-plastic problems. Such solvers for the computation of intercell fluxes, which would improve the computation of plastic flows and hence a better assessment of residual stresses and strains, are discussed in the next chapter.

Introduction

It has been shown throughout this manuscript that hyperbolic problems in solid mechanics are solved in a different manner depending on the numerical method employed. In particular, irreversible deformations which are usually numerically computed based on well-known constitutive integrators, may greatly differ from one scheme to another even for one-dimensional problems. However, the accurate assessment of residual stresses and strains are of major importance for many industrial applications such as high-speed metal forming, crashproof design or the study of the impact of earthquakes on structures. The simulations performed in chapter 4 emphasized the improvements enabled by the knowledge of the characteristic structure of the solutions of conservation laws, especially for elastoplastic solids. Nevertheless, the use of an elastic-plastic approximate Riemann solver is so far only possible for problems in one space dimension.

The purpose of this chapter is to identify typical behavior of the solutions of two-dimensional elastoplasticity problems under small strains. It is believed that the knowledge of these solutions will allow, through the better understanding of their mathematical features, the building of approximate numerical solutions embedding a sufficient amount of information in order to mimic the analytical behavior. This will be possible at a low computational cost provided that some key-properties of the exact solutions are clearly identified.

This chapter is organized as follows. A brief historical review of the solution of dynamic problems in two-dimensional elastic-plastic solids is made in section 5.1. Then, the equations of plasticity are recalled in section 5.2 so that the characteristic analysis, followed by the application of the method of characteristics, can be carried out. In section 5.3, attention is paid to the evolution of stress components inside simple waves that might propagate by means of a mathematical study of the ODEs satisfied within these waves. Since the developments rapidly become cumbersome, the analysis is supplemented with numerical results in section 5.4. At last, some identified trends are discussed at the end of the chapter in order to use them for building a dedicated Riemann solver.

Historical review

Until the 50s, research on dynamic problems in elastic-plastic solids were focused on uni-axial stress or strain, pure bending or pure torsion loading conditions [START_REF] Taylor | The plastic wave in a wire extended by an impact load[END_REF][START_REF] Von Kármán | On the propagation of plastic deformation in solids[END_REF], and were carried out for materials characterization purposes. The first references that brought some understanding about the response of linearly hardening solids to combined shear and pressure loads are those of Rakhmatulin [START_REF] Rakhmatulin | On the propagation of elastic-plastic waves owing to combined loading[END_REF] and Cristescu [START_REF] Cristescu | On the propagation of elasto-plastic waves for combined stresses[END_REF]. These early analytical investigations on plane stress impacts in the plastic regime led to the conclusion that elastic waves, as well as plastic combined-stress simple waves, can propagate in two-dimensional solids. While the former were well-known, the latter were shown to fall into two families: the fast waves and the slow waves.

Later, Bleich and Nelson [START_REF] Bleich | Plane waves in an elastic-plastic half-space due to combined surface pressure and shear[END_REF] considered superimposed plane and shear waves in an ideally elastic-plastic material submitted to step loads. It has then been highlighted that different loading cases yield different characteristic structures of the solution of a Picard problem, thus revealing the complexity of plastic flows in more than one dimension. The same conclusions have been drawn by Clifton [START_REF] Clifton | An analysis of combined longitudinal and torsional plastic waves in a thin-walled tube[END_REF] for hardening materials under tension-torsion, who furthermore studied the influence of plastic pre-loading on the solution. This contribution established the existence of loading paths through the simple waves arising from the characteristic analysis of the hyperbolic system. Indeed, the combined-stress wave nature lies in ODEs which govern the evolution of stress components within the simple waves. The integration of these equations of the form dσ 11 = ψdσ 12 allows the building of curves that connect the applied stress state of the Picard problem (σ d 11 , σ d 12 ) to the initial state of the medium. It has been for instance shown that if a solid is acted upon by a traction force such that σ d 11 = 0 and σ d 12 lies outside the elastic domain, only an elastic shear discontinuity, followed by a slow simple wave, propagates.

Conversely, other loading conditions may lead to the combination of an elastic pressure discontinuity and a fast wave, possibly followed by a slow wave. Another notable conclusion is that the combined loading paths followed inside simple waves may lead to plastic unloading, whereas only elastic unloading occurs in the one-dimensional theory.

Experimental data collected on a thin-walled tube submitted to a dynamic tensile load [START_REF] Lipkin | An experimental study of combined longitudinal and torsional plastic waves in a thin-walled tube[END_REF][START_REF] Lipkin | Plastic waves of combined stresses due to longitudinal impact of a pretorqued tube-part 1: Experimental results[END_REF] confirmed the existence of two distinct families of simple waves, both involving combined stress paths. These works nevertheless exhibited some discrepancies with the theory which have been attributed to the assumption made on the von-Mises yield surface. As a matter of fact, a constant strain region lying between the fast and slow waves that is predicted by the theory [START_REF] Clifton | An analysis of combined longitudinal and torsional plastic waves in a thin-walled tube[END_REF] could not be seen in experimental results. However, by following the endochronic theory of plasticity [START_REF] Valanis | A theory of viscoplasticity without a yield surface part i -general theory[END_REF] which does not require the introduction of a yield surface, Wu and Lin [START_REF] Han-Chin | Combined plastic waves in a thin-walled tube[END_REF] obtained numerical results that better fit the experimental data provided by Lipkin and Clifton [START_REF] Lipkin | Plastic waves of combined stresses due to longitudinal impact of a pretorqued tube-part 1: Experimental results[END_REF]. The good agreement showed between numerical and experimental results [START_REF] Han-Chin | Combined plastic waves in a thin-walled tube[END_REF] thus confirmed the theory.

Ting and Nan [START_REF] Ting | Plane waves due to combined compressive and shear stresses in a half space[END_REF] then generalized the work of Bleich and Nelson to hardening materials and Ting [START_REF] Ting | Elastic-plastic boundaries in plane and cylindrical wave propagation of combined stresses[END_REF] widened that of Clifton to more complex loadings, that is a superimposition of one plane wave and two shear waves states. Once again, the mathematical study of the ODE system governing the stress evolution inside fast and slow simple waves led to the construction of loading paths in stress space that depend on the external loads. A review of governing equations for all the cases depending on one space dimension considered above can be found in [START_REF] Nowacki | Stress waves in non-elastic solids[END_REF].

The information on characteristic structures thus provided has then been used by Lin and Ballman [START_REF] Lin | A riemann solver and a second-order godunov method for elastic-plastic wave propagation in solids[END_REF] for the development of an iterative Riemann solver. This procedure is based on successive guesses of the stress state lying in the stationary region so that the loading paths predicted by the theory of Clifton [START_REF] Clifton | An analysis of combined longitudinal and torsional plastic waves in a thin-walled tube[END_REF] can be integrated numerically until convergence. The implementation of this solver within a second-order Godunov scheme provided results that were in good agreement with the exact solutions. Nevertheless, the theoretical investigations mentioned above restrict the development of such numerical tools to problems that depend on one space dimension. Clifton tackled the solution of plane strain problems in elastic-plastic solids by looking for bi-characteristics [START_REF] Clifton | Analysis of dynamic deformation of elastic-plastic solids under condition of plane strain[END_REF] in order to build finite difference schemes that account for plastic waves. The point of view adopted here is that one can benefit from the simplifications introduced by the writing of Riemann problems in an arbitrary direction. Indeed, the method of characteristics rather than the more complex method of bi-characteristics can be employed with the quasi-linear forms presented in chapter 2.

On the other hand, the existence of plastic shocks in solids under plane wave assumptions has been investigated by several authors. First, Mandel [START_REF] Mandel | Ondes plastiques dans un milieu indéfini à trois dimensions[END_REF] showed the existence of stable plastic shocks in three-dimensional elastoplastic media. In this work, Hugoniot curves are built by assuming that the internal variables followed a radial loading path through a plastic shock. Lee and Germain [START_REF] Germain | On shock waves in elastic-plastic solids[END_REF] considered that Hugoniot curves in elastic-plastic solids cannot be constructed without studying the internal structure of the shock. Thus, an elastic-viscoplastic continuum problem is solved by magnifying the narrow region in the vicinity of the shock in which the fields vary sharply. The shock solution was then taken as the limit when viscosity tends to zero. A study of the internal structure of the shock has also been made by Stolz [START_REF] Dequiedt | Propagation of a shock discontinuity in an elasto-plastic material: constitutive relations[END_REF]. In the latter approach, the Hugoniot conditions across a shock moving at constant speed are derived by doing an asymptotic analysis. The author thus provided existence and uniqueness conditions for a shock in compression provided that elastic stiffening dominates the (concave) hardening saturation. Nevertheless, according to Mandel [START_REF] Mandel | Ondes de choc longitudinal dans un milieu élastoplastique[END_REF], such an analysis of the internal structure of the shock is not required to build Hugoniot curves, provided one chooses ε p 1 as internal variable and not the specific work w. However, the propagation of plastic shocks is still an open scientific issue and subject to debate.

In what follows, simple waves are considered in elastic-plastic solids with natural initial conditions by assuming a concave hardening law with no stiffening so that plastic shocks do not arise.

Elastic-plastic wave structure in two space dimensions

Key points of the chapter: The simple wave structure in two-dimensional solids is studied by writing the governing equations of plane strain and plane stress problems in a unique generic framework. Thus, the works carried out on simple waves mentioned above, formulated in terms of elastoplastic softnesses, appear as particular cases of the present investigation. As a result, we shall refer to those existing solutions in order to confirm our development as much as possible.

Governing equations

We are concerned with linear isotropic hardening materials whose elastic domain is given by the von-Mises yield surface, under isothermal deformations in the linearized geometrical framework. The balance equation of linear momentum with neglected body forces, and the geometrical balance equations are:

     ρ v -∇ • σ = 0 ε -∇ • v ⊗ I + I ⊗ v 2 = 0 (5.1)
In addition, the elastic-plastic constitutive equations derived from thermodynamics in section 2.2.3 are recalled here:

                         f (σ, A) = 3 2 s -(R(p) + σ y ) ≡ 0, with s = σ - 1 3 trace σI R(p) = C p σ = C elast -β s ⊗ s : ε = C ep : ε β = 6µ 2 3µ + C × 1 s : s A ep ij = A elast ij -β(n k s ki )(s jl n l )
(5.2a)

(5.2b)

(5.2c)

(5.2d)

(5.2e)

In the expression of the von-Mises yield function (5.2a), the (positive) linear isotropic hardening law (5.2b) is considered. Moreover, the elastoplastic acoustic tensor (5.2e) is decomposed as an elastic part A elast ij and a plastic part depending on the direction of the plastic flow through the coefficient β (5.2d). By inverting the (isotropic) elasticity tensor C involved in equation (5.2c), the following elastic law is written in the isotropic case:

ε e = 1 + ν E σ - ν E trace σI (5.3)
with Young's modulus E and Poisson's ratio ν.

The quasi-linear form of the sets of equations (5.1) and (5.2) in a Cartesian coordinate system and an arbitrary direction n is:

Q t + J ∂Q ∂x n = 0 (5.4)
where

x n = x•n, Q = v σ
, and J is the Jacobian matrix. It has been shown in section 2.3 that the 3 eigenvalues ω p and eigenvectors l p of the acoustic tensor lead to 6 left characteristic fields of the Jacobian matrix {c K ; L K } according to:

± ω p ρ ; ±ρ ω p ρ l p , -l p ⊗ n , p = 1, 2, 3 (5.5)
In addition, three independent left eigenvectors associated with the zero eigenvalue of system (5.4), which is of multiplicity 3, are found by solving:

σ K : (C ep • n) = 0, K = 1, 2, 3 (5.6)
The present formulation differs from those of Bleich [START_REF] Bleich | Plane waves in an elastic-plastic half-space due to combined surface pressure and shear[END_REF], Clifton [START_REF] Clifton | An analysis of combined longitudinal and torsional plastic waves in a thin-walled tube[END_REF], and hence these of Ting and Nan [START_REF] Ting | Plane waves due to combined compressive and shear stresses in a half space[END_REF] and Ting [START_REF] Ting | Elastic-plastic boundaries in plane and cylindrical wave propagation of combined stresses[END_REF], in that equation (5.4) is based on the elastoplastic stiffnesses rather than softnesses. As a consequence, it will be seen in what follows that the equations can be easily specialized to plane strain and plane stress cases.

Problems in two space dimensions

We now focus on the solid domain

x 1 × x 2 × x 3 ∈ [0, ∞[×[-h, h] × [-e,
e] in a Cartesian coordinate system, where e and h are arbitrary lengths. It is assumed that all quantities depend solely on x 1 and x 2 except the velocity component v 3 that may depend on x 3 . In particular, this is the case for e h. The solid is under plane strain conditions, that is ε • e 3 = 0, if the velocity v does not depend on x 3 and if v 3 vanishes. Thus, combining the additive partition of the infinitesimal strain tensor: ε = ε e + ε p , with the elastic law (5.3) and the kinematic condition ε 33 = 0, one gets:

σ 33 = ν (σ 11 + σ 22 ) -Eε p 33 (5.7)
Hence, the quasi-linear form (5.4) reduces for plane strain problems to a system of dimension 5 with unknowns v 1 , v 2 , σ 11 , σ 12 , and σ 22 .

Alternatively, a plane stress state (σ • e 3 = 0) is assumed if the planes x 3 = ±h are traction free and e h. As a result, the stress component σ 33 can be removed from the system (5.4). Nevertheless, the tangent modulus must account for the vanishing out-of-plane stress component. Specialization of equation ( 5 The characteristic structure of the problem is then given by the associated acoustic tensor A ep = n • C ep • n. The removal of σ 33 from system (5.4) for both plane strains and plane stresses allows solving the problem in a two-dimensional setting. Then, generically denoting the acoustic tensor by A, the characteristic structures are given by the eigenvalues:

ω 1 = 1 2 A 11 + A 22 + (A 11 -A 22 ) 2 + 4A 12 2
(5.9a)

ω 2 = 1 2 A 11 + A 22 -(A 11 -A 22 ) 2 + 4A 12 2
(5.9b) and the associated eigenvectors:

l 1 = [A 22 -ω 1 , -A 12 ] ; l 2 = [-A 12 , A 11 -ω 2 ]
(5.10)

From equation (5.5), we see that two families of waves with celerities c f = ± ω 1 /ρ and c s = ± ω 2 /ρ may travel in the domain. These waves are respectively referred to as fast and slow waves. Note that subtracting equations (5.9a) and (5.9b) leads to:

ρc 2 f -ρc 2 s = (A 11 -A 22 ) 2 + 4A 12 2 ≥ 0 (5.11)
Hence, the characteristic speed associated with fast waves is always greater than or equal to that of slow waves.

The four left eigenfields of the Jacobian matrix thus read:

±c f ; L c ± f = ±ρc f l 1 , -l 1 ⊗ n ±c s ; L c ± s = ±ρc s l 2 , -l 2 ⊗ n
where L c + f and L c - f are associated with the right-going and left-going fast waves respectively. The same goes for L c + s and L c - s . Furthermore, one stationary wave associated with the zero eigenvalue of the Jacobian matrix, and whose left eigenvector satisfies equation (5.6), has to be added:

L 0 T =          v 0 1 v 0 2 σ 0 11 σ 0 22 σ 0 12          =          0 0 (C 121i C 222j -C 221i C 122j ) n i n j (C 111i C 122j -C 112i C 121j ) n i n j (C 112i C 221j -C 111i C 222j ) ninj 2          =       0 0 α 11 α 22 α 12       (5.13)
with C = C ep for plain strain and C = C ep for plane stress. It has been seen in section 2.4.4 that the solution of non-linear problems may contain shock and/or simple waves. Nevertheless, we restrict here to simple waves by assuming that: (i) the characteristic speeds satisfy

c 1 ≥ c f ≥ c 2 ≥ c s ,
where c 1 and c 2 are the speeds of elastic pressure and shear discontinuities respectively; (ii) c f and c s monotonically decrease with the hardening of the material; (iii) the computational domain is in an initial natural, plastic strain free state.

The characteristic equations L K • dQ = 0 are then written: Integration of equations (5.14a) to (5.14e) leads to integral curves through simple waves in which several stress components vary, hence the name combined-stress simple waves [START_REF] Cristescu | On the propagation of elasto-plastic waves for combined stresses[END_REF]. Following [START_REF] Clifton | An analysis of combined longitudinal and torsional plastic waves in a thin-walled tube[END_REF], the method of characteristics is applied by combining equations (5.14a) to (5.14e). The approach consists in tracing every characteristic from some downstream point of a wave where the state vector Q is known, to an upstream point where the solution is sought. Figures 5.1a and 5.1b schematically illustrate the method for slow and fast simple waves in which the state is known along the head wave and is looked for at point P lying on the tail wave. The integral curves through slow and fast simple waves are derived in the next section.

ρc f l 1 • dv -l 1 i n j dσ ij = 0 along dx/dt = c f (5.14a) -ρc f l 1 • dv -l 1 i n j dσ ij = 0 along dx/dt = -c f (5.14b) ρc s l 2 • dv -l 2 i n j dσ ij =

Integral curves through simple waves

The right-going slow waves are first looked at by adding equations (5.14a) and (5.14b):

l 1
i n j dσ ij = 0 Given the geometry of the problem, the vector n may be reduced to e 1 or e 2 . It therefore comes out:

dσ 11 = - l 1 2 l 1 1 dσ 12 = ψ s 1 dσ 12
for n = e 1 (5.15a)

dσ 22 = - l 1 1 l 1 2 dσ 12 = ψ s 2 dσ 12
for n = e 2 (5.15b)

where ψ s 1 and ψ s 2 are functions of all components of σ. Moreover, the s and f superscripts stand for slow and fast waves respectively in the remainder of the manuscript. With the above equations, the characteristic equation related to the contact wave (5.14e) reads:

dσ 22 = - ψ s 1 α 11 + α 12 α 22 dσ 12
for n = e 1 (5.16a)

dσ 11 = - ψ s 2 α 22 + α 12 α 11 dσ 12 for n = e 2 (5.16b)
The sets of equations (5.15a)-(5.16a) and (5.15b)-(5.16b) show the combined-stress nature of slow simple waves.

On the other hand, the subtraction of equations (5.14a) and (5.14b) leads to:

dv 1 = ψ s 1 dv 2 = 1 ψ s 2 dv 2
which, once combined with equations (5.15a)-(5.15b) and introduced in (5.14d), yields after simplifications: for n = e 2 (5.17b)

dv 1 = -

Remark 17

The integral curves through a left-going slow wave result from the combination of equations (5.15a)-(5.15b) introduced in (5.14c) rather than (5.14d). Therefore, the only difference lies in the signs in equations (5.17a) and (5.17b).

Similar results are obtained for right-going fast simple waves by using l 2 instead of l 1 and c f rather than c s . Hence, the evolution in slow and fast waves is governed by the loading functions:

ψ s 1 = - l 1 2 l 1 1 n=e1 , ψ s 2 = - l 1 1 l 1 2 n=e2 , ψ f 1 = - l 2 2 l 2 1 n=e1 , ψ f 2 = - l 2 1 l 2 2 n=e2
(5.18)

The equations satisfied across right-going slow and fast simple waves are summarized in table 5.1. The integration of equations gathered in table 5.1 should provide the complete solution of a given problem by means of integral curves, or loading paths. For instance, the velocity resulting from the passage of right-going waves in the direction e 1 obeys:

v 1 = v 0 1 - σ σ 0 dσ 11 ρc 2 ; v 2 = v 0 2 - σ σ 0 dσ 12 ρc 2 (5.19)
where the zero superscript denotes the downstream state. Nevertheless, Clifton [START_REF] Clifton | An analysis of combined longitudinal and torsional plastic waves in a thin-walled tube[END_REF] emphasized that depending on the loading conditions, only one simple wave or both may arise in the solution. Therefore, it is crucial to identify the stress path followed to properly compute integrals (5.19). This is the purpose of the next section.

dv 1 = -dσ12 ρc 2 s dv 1 = -dσ11 ρc 2 f dv 1 = -dσ12 ρc 2 f dv 2 = -dσ12 ρc 2 s dv 2 = -dσ22 ρc 2 s dv 2 = -dσ12 ρc 2 f dv 2 = -dσ22 ρc 2 f dσ 11 =

Loading paths through simple waves

Properties of the loading paths

The stress paths followed within slow and fast simple waves are governed by the mathematical properties of the loading functions (5.18). Before specializing the discussion to plane stress and plane strain cases, some general properties holding regardless of the loading conditions are highlighted. First, the functions satisfy the orthogonality properties:

ψ s 1 ψ f 1 = -1 and ψ s 2 ψ f 2 = -1.
Indeed, considering the left eigenvectors of the acoustic tensor given in equation (5.10), the product ψ s 1 ψ f 1 reads:

ψ s 1 ψ f 1 = l 1 2 l 1 1 l 2 2 l 2 1
Since the eigenvectors of symmetric second-order tensors all satisfy l 1 •l 2 = 0, it comes out that the above product is equal to -1. Whereas this orthogonality has already been noticed for particular plane strain and plane stress cases [START_REF] Clifton | An analysis of combined longitudinal and torsional plastic waves in a thin-walled tube[END_REF][START_REF] Ting | Plane waves due to combined compressive and shear stresses in a half space[END_REF], the generic formulation proposed here shows that this is valid for all problems in two space dimensions. As a result, the study can be restricted to one function in each direction, say ψ s 1 and ψ s 2 .

Second, if the function ψ s 1 vanishes at some point of stress space, the projection in the (σ 11 , σ 12 ) plane of the loading path followed within a slow wave is vertical according to the ODE (5.15a) (i.e dσ 11 = 0). Conversely, if ψ s 1 → ∞, the loading path is horizontal in the (σ 11 , σ 12 ) plane (i.e dσ 12 = 0). These situations respectively correspond to:

ψ s 1 = 0 ⇔ A 12 = 0 (5.20a) ψ s 1 → ∞ ⇔ A 22 -ω 1 = 0 (5.20b)
In particular, if A 12 = 0 equation (5.20b) reads:

A 22 -ω 1 = 1 2 A 22 -A 11 -(A 11 -A 22 ) 2 + 4A 2 12 = -A 11 -A 22
where • denotes the positive part operator. Hence, if A 12 = 0 and A 11 = A 22 , one has ψ s 1 = 0 and hence ψ f 1 → -∞. If moreover A 11 = A 22 , both components of the eigenvectors vanish and the functions ψ s 1 and ψ f 1 are undetermined. At last, it follows from equation (5.11) that the simultaneous satisfaction of conditions (5.20a) and (5.20b) leads to characteristic speeds of simple waves that are identical. Hence, the situation c f = c s corresponds to a loss of hyperbolicity of the system.

Analogously, the function ψ s 2 is such that:

ψ s 2 → ∞ ⇔ A 12 = 0 (5.21a) ψ s 2 = 0 ⇔ A 22 -ω 1 = 0 (5.21b)
Therefore, if both conditions (5.21a) and (5.21b) are satisfied on A 12 and A 22ω 1 , the system is no longer hyperbolic with characteristic speeds of fast and slow waves that are identical.

ρc 2 f = 1 2 λ + 3µ -βs 2 11 + λ + µ -βs 2 11
Two cases are to be considered:

(i) if βs 2 11 < λ + µ, the expression further reduces to:

ρc 2 s = µ ρc 2 f = λ + 2µ -βs 2 11
The characteristic speed of slow waves is therefore equivalent to that of elastic shear waves for plane strain

c s = c 2 = µ/ρ.
(ii) if λ + µβs 2 11 < 0, the characteristic speeds read:

ρc 2 s = λ + 2µ -βs 2 11
ρc 2 f = µ Therefore, the celerity of fast waves reduces to that of elastic shear waves. Note, however, that the characteristic speed of slow waves remains real if and only if λ + 2µ > βs 2 11 . One then gets the following bounds: λ + 2µ > βs 2 11 > λ + µ. At last, the equality βs 2 11 = λ + µ leads to A ep 11 -A ep 22 = 0 and hence, to undetermined loading functions.

The plane stress case

The elastoplastic tangent modulus under consideration is now that given in equation (5.8). Let's first consider ψ s 1 related to the vector n = e 1 . Thus: (5.27c)

A ep 11 = C ep 1111 - (C ep 1133 ) 2 C ep 3333 = λ + 2µ -
In order to ensure the hyperbolicity of the system, the components of the acoustic tensor also have to be defined, that is C ep 3333 > 0. This condition leads to:

λ + 2µ -βs 2 33 > 0 ⇔ s 2 33 < λ + 2µ β
Second, from equation (5.27c), A ep 12 admits two roots in terms of the components of the deviatoric stress tensor, namely:

s 12 = 0 ; s 11 = λ λ + 2µ s 33
In terms of the components of the Cauchy stress tensor, these conditions read:

σ 12 = 0 ; σ 11 = 2µ 3λ + 4µ σ 22 (5.28)
If on the other hand the vector n = e 2 is considered, the acoustic tensor components read: The complexity introduced by the plane stress tangent modulus prevents finding other singular configurations for the hyperbolic system. In particular, it is difficult to deal with the equation A ep 11 = A ep 22 due to the expressions given in equations (5.27a) and (5.27b). Nevertheless, since the stress state s 12 = 0 also constitutes a singular point for plane stress, the same remarks as those made for the plane strain loading path hold. Namely, σ 12 becomes constant if it falls to zero along the loading path followed inside a fast (resp. slow) wave propagating in direction e 1 (resp. e 2 ) as summarized in table 5.2.

A ep 11 = C ep 1212 - (C ep 1233 ) 2 C ep 3333 = µ -

Numerical integration of stress paths

Although some properties of the simple waves have been emphasized in section 5.3, the complexity of the equations prevents the complete characterization of the loading paths followed. In order to get additional information on the evolution of the stress states, the systems of ODEs gathered in table 5.1 are numerically integrated for plane stress and plane strain loadings. In particular, the thin-walled tube problem considered by Clifton [90] is first looked at so that the above developments can be validated. Next, the plane stress and plane strain cases are treated. The values of the elastic properties considered here are those used in the previous chapter (see table 4.1). On the other hand, the tensile yield stress σ y = 1 × 10 8 P a and the hardening modulus C = 1 × 10 8 P a are set arbitrarily. Finally, we restrict to positive shear stress σ 12 ≥ 0.

Thin-walled tube problem

Consider the semi-infinite domain in the Cartesian coordinate system:

x 1 × x 2 × x 3 ∈ [0, ∞[×[-h, h] × [-e, e],
being acted upon by a traction vector T d at x 1 = 0 and free surfaces x 2 = ±h and x 3 = ±e. Only the first two components of T d are non-null so that the stress and strain tensors within the medium are of the form:

σ =   σ 11 σ 12 0 0 0 sym 0   ; ε =   ε 11 ε 12 0 ε 22 0 sym ε 33  
By using the following mapping of coordinates: (1, 2, 3) → (z, θ, r), such a state corresponds also to that holding in a hollow cylinder with radius and length much bigger that its thickness, submitted to combined longitudinal and torsional loads. Hence the name thin-walled tube problem. As a particular plane stress case, the set of ODEs along characteristics derived in section 5.3 applies by taking into account the vanishing stress component

σ 22 : σ22 = C ep 22ij εij = 0 i, j = {1, 2} ⇒ C ep 2222 ε22 = -C ep 22ij εij ij = {11, 12, 21}
where C ep is the plane stress tangent modulus (5.8). Thus, inverting the above equation and introducing it in the constitutive equation, we are left with the following law:

σij = C ep ijkl εkl - C ep ij22 C ep 22kl C ep 2222 εkl = C ep ijkl εkl ij, kl = {11, 12, 21} (5.30)
The characteristic analysis of the hyperbolic system based on this tangent modulus also leads to loading paths followed across slow and fast waves, involving however two components of stress rather than three. For the sake of simplicity, the stress components are denoted by σ 11 = σ and σ 12 = τ whereas the velocity components reads

v 1 = u and v 2 = v.
Thus, the ODEs governing the evolution of stress components inside the waves of combined-stress read:

dσ = ψ s,f dτ (5.31)
where the loading functions ψ s,f depend on the component of the acoustic tensor that corresponds to the tangent modulus (5.30). Equations (5.31) as well as those of Clifton [START_REF] Clifton | An analysis of combined longitudinal and torsional plastic waves in a thin-walled tube[END_REF] have been numerically integrated, starting from several arbitrary points lying on the initial yield surface. Since the loading functions are odd functions of σ and τ [START_REF] Clifton | An analysis of combined longitudinal and torsional plastic waves in a thin-walled tube[END_REF], τ (σ) and σ(τ ) are even functions and hence, the loading paths exhibit symmetries with respect to τ and σ axes. Therefore, the study is restricted to the quarter-plane (σ > 0, τ > 0). Figure 5.2 shows one stress path resulting from the integration of the ODE related to right-going fast waves with σ used as a driving parameter. The initial stress state lies on the yield surface at σ = 0 and the ODE is discretized by means of the backward Euler method, the integration being performed until the stress reaches the value σ = σ y . The path is respectively depicted in the stress space and in the deviatoric plane in figures 5.2a and 5.2b. The deviatoric plane projection is obtained by drawing the paths in the eigenstress space and projecting them onto the plane perpendicular to the hydrostatic axis σ 1 = σ 2 = σ 3 . In this plane, the von-Mises yield surface is a circle drawn with dashed lines. As observed by Clifton, the path inside fast waves first follows the initial yield surface up to the intersection with the σ-axis. Then, the loading path is such that dτ = 0 while σ increases as far as hyperbolicity holds, that is for c f < c 2 = µ/ρ [START_REF] Clifton | An analysis of combined longitudinal and torsional plastic waves in a thin-walled tube[END_REF]. Notice that these conclusions are similar to those made in the previous section. The ODEs derived in section 5. Comparison between the results obtained from equations (5.31) and these of [START_REF] Clifton | An analysis of combined longitudinal and torsional plastic waves in a thin-walled tube[END_REF].

Adopting the same approach with τ as driving parameter, some stress paths through slow waves have been reported in figure 5.3. Since fast waves lead to loading paths following the initial yield surface, the orthogonality property of the loading functions implies that those of slow waves move away from it. This is seen in figure 5.3a. This property holds in the (σ, τ ) plane but not in the deviatoric plane, as can be seen in figure 5.3b, since the quasi-linear form (5.4) and hence, the ODEs, are not written in terms of s 1 , s 2 , s 3 . Comparison between the results obtained from equations (5.31) (cross markers) and those of [START_REF] Clifton | An analysis of combined longitudinal and torsional plastic waves in a thin-walled tube[END_REF] (solid lines).

The behaviors highlighted above allow the solution of the Picard problem in a thin-walled cylinder, that is:

• initial conditions σ(x, t = 0) = 0, v(x, t = 0) = 0

• step-loading boundary conditions σ(x 1 = 0, t) = σ d and τ (x 1 = 0, t) = τ d Indeed, with given (σ d , τ d ) outside of the initial yield surface, one can first integrate backward the loading path through the slowest wave. Two situations are then distinguished: (i) if the integration leads to some point of the initial yield surface, which can be reached by elastic discontinuities, the solution is complete.

(ii) if the slow wave connects (σ d , τ d ) to the σ-axis at some point lying outside of the initial yield surface, then a fast wave must be integrated backward to the initial elastic domain. Indeed, analogously to the results of table 5.2, it is shown in [START_REF] Clifton | An analysis of combined longitudinal and torsional plastic waves in a thin-walled tube[END_REF] that the paths followed through slow waves (resp. fast waves) are perpendicular (resp. parallel) to the σ-axis. As a result and by virtue of the symmetries with respect to the σ and τ axes of the loading paths, the initial yield surface can be reached by considering a fast wave.

At last, the cases τ d = 0 and σ d = 0 respectively lead to one single fast wave and one single slow wave. Once the characteristic structure of the problem has been determined (i.e. one fast wave, one slow wave, or both), the complete set of ODEs can be integrated in order to compute the solution. It is worth emphasizing the complexity introduced by waves of combined-stress since the characteristic structure of the solution of a Picard problem now depends on the boundary conditions. Hence, for developing a Riemann solver that would provide the stationary solution, additional computational effort must be made.

The iterative procedure proposed by Lin and Ballman [START_REF] Lin | A riemann solver and a second-order godunov method for elastic-plastic wave propagation in solids[END_REF] to solve Riemann problems is based on the above considerations. The left and right initial conditions of that problem satisfy equations similar to (5.19):

u * = u L + σ * σ L dσ ρc ; v * = v L + σ * σ L dτ ρc (5.32a) u * = u R - σ * σ R dσ ρc ; v * = v R - σ * σ R dτ ρc (5.32b)
where the asterisk denotes the stationary state of the Riemann problem. First, a stress state (σ, τ ) is assumed to be connected to σ L and σ R (see figure 5.4 for the illustration of the method). By taking into account the

u σ u L σ L u 1 u R σ R u 2 σ σ v τ v L τ L v 1 v R τ R v 2 τ τ Figure 5
.4: Schematic representation of the iterative Riemann solver proposed in [START_REF] Lin | A riemann solver and a second-order godunov method for elastic-plastic wave propagation in solids[END_REF]. conclusions of Clifton the loading paths followed can be identified so that equations (5.32a) and (5.32b) can be integrated. Thus, virtual integral curves are built in (u, σ) and (v, τ ) planes as depicted with dashed lines in figure 5.4. Second, the intersection of the curves joining respectively v L to v 1 and v R to v 2 gives a stress state ( σ, τ ) that is used to apply the procedure again until some criterion v 1v 2 ≤ is achieved. At last, the state obtained ( v, σ) corresponds to the stationary state of the Riemann problem and can be used to compute numerical fluxes at cell interfaces. Notice that in this procedure, the intersection of integral curves is found by means of the tangent lines approximation at σ so that this solver does not fully account for the exact solution. Therefore, it appears that the loading paths identified in [START_REF] Clifton | An analysis of combined longitudinal and torsional plastic waves in a thin-walled tube[END_REF] are crucial in order to determine the wave pattern corresponding to a guessed stationary state. This is exactly what is lacking for the more general plane strain and plane stress problems.

Plane stress

We now move on to a more general plane stress case for which the stress σ 22 is not zero. Although the equations of section 5.3 have been derived for two directions of propagation, attention is paid here to n = e 1 only. Thus, the system of ODEs considered reads (see table 5.1):

dσ 11 = ψ s,f
1 dσ 12 (5.33)

dσ 22 = - ψ s,f 1 α 11 + α 12 α 22 dσ 12
One path through a fast simple wave is first looked at by assuming an initially free-stress state, brought to the yield surface at the point σ 11 = σ 22 = 0. Equations (5.33) planes, while the projection in the deviatoric plane can be seen in figure 5.5b. In addition, figure 5.5a shows the evolution of the characteristic speed associated with the fast wave along the path by means of a color gradient. Thus, it can be seen that for the loadings under consideration, the wave celerity is a decreasing function of the stress so that the simple wave solutions are valid. Next, it appears that σ 12 is an even function of σ 11 and σ 22 .

At last, analogously to the thin-walled cylinder solution, the stress components follow the initial yield surface, which is obvious in the deviatoric plane (figure 5.5b). Furthermore, according to the property (5.28) the stress path must be horizontal in the (σ 11 , σ 12 ) and (σ 22 , σ 12 ) planes, once the σ 11 -axis is reached. As depicted in figure 5.5b, this point corresponds to a direction of pure shear in the deviatoric plane. Nevertheless, the numerical integration of ODEs once the shear stress σ 12 vanishes is not possible owing to an indeterminacy of the loading function ψ f 1 that has not been identified so far. We now focus on the stress evolution inside slow waves. The same procedure is followed for several starting points on the initial yield surface. In addition, various initial values are considered for σ 22 since, even for a solid in a free stress state at t = 0, a fast wave may lead to σ 22 = 0. The loading paths thus obtained for the (a) Projections of loading paths in (σ 11 , σ 12 ) and (σ 22 , σ 12 ) planes (arbitrary) initial values σ 22 = -5.8 × 10 7 P a, σ 22 = 0 and σ 22 = 5.8 × 10 7 P a, are respectively depicted in figures 5.6, 5.7 and 5.8. The projections in the stress space and the deviatoric plane are shown. The evolution of the characteristic speed associated with the slow wave can also be seen by means of a color gradient. Once again, the simple wave solution appears to be valid with the considered loading conditions. Furthermore, one can see that the stress paths are now more complex since, for instance, no symmetry appears in the (σ 11 , σ 12 ) plane. Moreover, the behavior is even more complex in the (σ 22 , σ 12 ) plane in which the variations first mainly concern σ 22 and next, the slopes of curves roughly change so that the paths are almost vertical. This sharp change in slopes is also notable in the deviatoric plane in figure 5.6b. On the other hand, similar observations can be made for the other initial values of σ 22 as can be seen in figures 5.7 and 5.8. However, the paths depicted in the (σ 22 , σ 12 ) plane in figure 5.8a follow a direction opposite (a) Projections of loading paths in (σ 11 , σ 12 ) and (σ 22 , σ 12 ) planes to those corresponding to the initial condition σ 22 = -5.8 × 10 7 P a. The previous remark is also valid in the deviatoric plane in figure 5.8b. Indeed, the integral curves first describe clockwise curved lines until the break in slopes occurs, after which a behavior close to straight lines is seen. It is worth noticing that the driving parameter used for slow waves has not been chosen arbitrarily. As a matter of fact, the numerical integration does not go well when driven with σ 11 . Furthermore, numerical issues occur if the starting point is such that σ 12 = 0. More generally, the loading paths resulting from the integration of ODEs governing the behavior inside simple waves in plane stress can be summarized as follows. Whereas the integral curves inside a fast wave first exhibit a phase in which the stress is restricted to the initial yield surface, the passage of a slow wave makes the stress leave the elastic domain almost instantaneously. It is moreover noteworthy that the shear stress component σ 12 undergoes the biggest variation through a slow wave, in spite of the visible combined-stress nature of the corresponding paths. In addition, roughs change in the slopes of the integral curves associated with slow waves occur. However, such phenomena might also be observed for fast waves once the shear stress vanishes but have not been highlighted in figure 5.5 due to numerical issues in integrating an undetermined function. Furthermore, after the slopes of slow wave integral curves broke, the paths are straights the deviatoric plane.

Plane strain

Assuming that a solid initially at rest undergoes external loads leading to a plane strain case, the previous approach is now repeated. However, the derivation of the hyperbolic system in a two-dimensional setting relies on the writing of the out-of-plane stress component as a function of plastic strain. Hence, the integral curves associated with simple waves are integrated implicitly, along with the plastic flow. To do so, the consistency condition (see section 2.2.3) of the von-Mises yield surface (5.2a) is written:

ḟ = 0 ⇔ 3 2 s : σ s = C ṗ
In addition, combination of the plastic flow rule (2.40a): εp = ṗ along with the ODE related to the out-of-plane component:

dσ 33 = ν (dσ 11 + dσ 22 ) -Edε p 33
Since a fast wave propagates faster than a slow one, a material particle is first acted upon by the effects of the former. Thus, figure 5.9 shows the evolution of stresses resulting from numerical integration using σ 11 as a driving parameter by rearranging equations (5.34). Several starting points on the initial yield surface are considered. The evolution of the celerity of fast waves along the integral curve confirms the validity of the simple wave solution. Furthermore, the starting points are chosen in such a way that a symmetry of the loading path with respect to σ 11 = 0 and σ 22 = 0 planes is notable. Although the stress paths depicted in figure 5.9a are rather different to those resulting from a fast wave in plane stress, the behavior in the deviatoric plane is similar as can be seen in figure 5.9b. Indeed, the fact that the computed loading paths through a fast wave is parallel to the initial yield surface is obvious when looking at the deviatoric plane. More specifically, the von-Mises circle is traced by the integral curve even when the shear component σ 12 is zero (see paths 1 and 6 in figure 5.9). Notice that the integration of the loading path in the plane σ 12 = 0 is here possible contrary to the plane stress. The integral curves of figure 5.9 however exhibit a cusp which is not explained (see the two external arrows pointing toward axes of pure tensile/compression).

On the other hand, some loading paths resulting from the integration of slow wave ODEs are depicted in figures 5.10, 5.11 and 5.12. Analogously to the plane stress case, three initial values σ 22 = -1.3 × 10 8 P a, σ 22 = 0 and σ 22 = 1.3 × 10 8 P a are considered since a fast wave may modify the initial free-stress state. The integration of loading paths through slow waves in plane strain is performed by using σ 12 as a driving parameter. However, numerical difficulties arise owing to the characteristic speed associated with slow waves which start increasing rather than decreasing at some point along the path. In order to circumvent this issue, the last stress state leading to a decreasing celerity is used as an initial condition for a second integration driven by means of σ 11 . The final value is set so that the variation of σ 11 (i.e increasing or decreasing) undergone up to that singularity is continued. This strategy allows carrying on the integration further. Nevertheless, the same problem of increasing characteristic speed again occurs and the computation must be aborted. The integral (a) Projections of loading paths in (σ 11 , σ 12 ) and (σ 22 , σ 12 ) planes Whereas σ 11 varies little as shown by the projection of the path in the (σ 11 , σ 12 ) plane, this is not the case for σ 22 (see figure 5.10a). Indeed, the projections of the integral curves in the (σ 22 , σ 12 ) plane exhibit complex paths similar to those obtained for plane stress. In contrast, it can be seen in figure 5.10b that the paths first follow the initial yield surface and next a direction of pure shear in the deviatoric plane, which differs from plane stress solutions. Even though the paths for plane stress and plane strain have similar shapes in the stress space, one cannot expect the same observations in the deviatoric plane due to the out-of-plane stress component. Namely, the loading paths are restricted to the (σ 1 , σ 2 ) plane under plane stress, while they can take values in the whole space (σ 1 , σ 2 , σ 3 ) under plane strain.

The same behavior is observed in figures 5.11 and 5.12 which respectively show the loading paths resulting from the zero and the positive initial values of σ 22 . Nevertheless, the integral curves in figure 5.11 reveal that numerical issues occur "faster" than in the previous case. Indeed, the characteristic speeds quickly start increasing so that the stress paths depicted are short. Furthermore, the projection in the (σ 22 , σ 12 ) plane of the integral curves 2 and 3 show that the slow wave mainly influences σ 22 . The projections in the deviatoric plane in figure 5.11 however show that the stress paths first follow the initial yield surface until the direction of pure shear is reached, and next follow the radial direction. In addition, considering figures 5.10, 5.11 and 5.12, it can be seen that if the initial condition on σ 11 is greater than the value corresponding to the maximum shear stress on the initial yield surface, σ 22 increases along the loading path. Conversely, σ 22 decreases along the loading path if σ 11 is initially lower than the value that corresponds to the maximum shear stress on the initial yield surface. Note that the same goes for σ 11 . Thus, it seems that the property sign(dσ 22 ) = sign(dσ 11 ) holds along the loading paths followed inside a slow simple wave in plane strain, though this has not been proved mathematically. Generally speaking, it appears that for the ranges of stress considered here, the hardening of the material is mainly due to slow simple waves for plane strain cases. Indeed, the latter may lead to radial loading paths that greatly increase the radius of the von-Mises cylinder in principal stress space, whereas the integral curves corresponding to fast waves are restricted to the initial yield surface. Notice, however, that the above results have been obtained by using a rather low hardening modulus.

The curves resulting from the use of the hardening parameter C = 1 × 10 10 P a exhibit slight differences. First, the integral curves of fast waves are depicted in figure 5.13. As before, the integral curves follow the initial yield surface but then branch off to reach a direction of pure tensile/compressive loading. Moreover, the visible cusps in the loading paths 2, 3, 4 and 5, which already arise in the solutions based on a lower hardening modulus, indicate that the path followed inside fast waves converges to a direction of pure tension/compression under plane strain. The increase in hardening parameter also allows the elimination of the integration issues that occur for a lower one. As a result, the stress paths followed inside slow simple waves depicted in figure 5.14 lead to stress states lying further outside of the initial elastic domain. Moreover, for every initial value of σ 22 considered in σ y -σ y the figures 5.14a, 5.14b and 5.14c, the slopes of the integral curves no longer break but smoothly vary to reach the direction of pure shear in the deviatoric plane.

Remark 18

The use of a higher hardening parameter for slow waves under plane stress also leads to smoother paths in the deviatoric plane. On the other hand, the integral curves associated with the fast waves under plane stress (slightly) move away from the yield surface up to a direction of pure shear. Nevertheless, at that point numerical difficulties occur due to the indeterminacy of the loading functions already mentioned.

Conclusion

Summary of the chapter

In this chapter, the characteristic structure of the solution of hyperbolic problems in elastic-plastic solids in two space dimensions has been highlighted. It is known since the 50s that plastic flow in two-dimensional solids yields two families of waves whose speeds depend on the stress state, the slow and fast waves. In addition, these plastic waves may have an impact on all stress components in contrast to elastic discontinuities, hence the name of combined-stress waves. During the 60s, attention has been paid to simple waves in particular two-dimensional problems thus providing, among others, solutions of Picard problems in an elastic-plastic medium undergoing step loadings [START_REF] Clifton | An analysis of combined longitudinal and torsional plastic waves in a thin-walled tube[END_REF][START_REF] Ting | Plane waves due to combined compressive and shear stresses in a half space[END_REF][START_REF] Ting | Plastic wave propagation in linearly work-hardening materials[END_REF]. The singular nature of such problems lies in the fact that the characteristic structure of the solution depends on the external loading undergone. Indeed, it has been shown [START_REF] Clifton | An analysis of combined longitudinal and torsional plastic waves in a thin-walled tube[END_REF] that boundary conditions can lead to plastic flow involving one fast, one slow, or both simple waves. Therefore, it is crucial to be able to identify typical stress paths followed in each simple wave in order to link the initial data to a given stress state, and subsequently to determine the occurring wave pattern.

Based on these works, an iterative Riemann solver [START_REF] Lin | A riemann solver and a second-order godunov method for elastic-plastic wave propagation in solids[END_REF], whose procedure has been recalled in section 5.4, has been developed for the numerical solution of the thin-walled tube problem. This solver relies on the ability to connect a stationary state to initial data by a characteristic wave pattern. Following this approach, identifying characteristic wave patterns for general elastoplastic problems in two space dimensions should allow to enrich the numerical solution with the knowledge of physics. For that purpose, the characteristic analysis of twodimensional problems in elastic-plastic materials with linear isotropic hardening under plane strain and plane stress, in projection in an arbitrary direction of space, has been carried out in section 5.2. Fast and slow waves are also involved in the solution so that applying the method of characteristic through the simple waves provides a system of ODEs. Integration of this system leads to integral curves in terms of velocity and stress components that are followed inside the combined-stress waves. Specializing the ODEs to one direction of a Cartesian grid, it has been shown in section 5.3 that the loading paths satisfied through slow and fast waves are perpendicular in the stress space for both plane strain and plane stress. Moreover, it has been established that the stress paths exhibit particular behavior in the space (σ 11 , σ 22 , σ 12 ), that is dσ 11 = 0, dσ 12 = 0 or dσ 22 = 0, for special values of the components of the acoustic tensor. These situations are achieved for different stress states depending on whether the problem involves plane stresses of plane strains as shown in section 5.3.

The complexity of the ODEs derived in section 5.2 prevents identifying all the singularities which may occur along the loading paths. Hence, the mathematical analysis has been supplemented with numerical results consisting of the integration of stress paths from arbitrary initial stress values lying on the initial yield surface, for the particular direction e 1 . First, in section 5.4.1 the loading paths resulting from the integration of the ODEs derived in section 5.2 have been compared to those of Clifton [START_REF] Clifton | An analysis of combined longitudinal and torsional plastic waves in a thin-walled tube[END_REF]. The two different formulations, respectively based on elastoplastic stiffnesses and softnesses, show good agreement.

Second, the evolution of stress components across fast and slow waves under plane stress has been looked at in section 5.4.2. It appears that though the loading paths are rather complex in stress space through a fast wave, the stress evolution in the deviatoric plane is restricted to the initial yield surface until one direction of pure shear is reached. A singularity then occurs so that the numerical integration cannot be pursued. On the other hand, the loading paths resulting from the integration of ODEs satisfied inside a slow wave exhibit complex shapes along which σ 11 varies much less than the other stress components.

Third, the plane strain case has been considered in section 5.4.3. Once again, the integral curves inside a fast wave show complex shapes in stress space, and an evolution restricted to the initial yield surface in the deviatoric plane. In that case, however, the paths may follow a direction of pure tension/compression in the latter plane so that the plastic flow is radial for high values of the hardening modulus. In contrast, the paths inside slow waves first rotate on the yield surface and then lead to a stress state of pure shear in the deviatoric plane.

Towards a two-dimensional elastoplastic Riemann solver

The physical structures emphasized in this chapter enable a better understanding of the propagation of waves in two-dimensional elastoplastic media, although further investigations are required. On the other hand, the loading paths followed in fast and slow simple waves can be used in order to improve the numerical simulation of these problems.

One possibility is to generalize the approach proposed by Lin and Ballman [START_REF] Lin | A riemann solver and a second-order godunov method for elastic-plastic wave propagation in solids[END_REF] based on the clues provided above. The idea would be to successively assume stationary states of the Riemann problem in terms of stress σ 11 , σ 12 and σ 22 in order to build stress paths starting from the initial data. Namely, considering the direction e 1 , the loading paths followed through a slow wave can be integrated backward starting from the guessed state. Then, different situations may occur:

(1-a) the curve thus obtained crosses the initial yield surface at a point where σ 22 satisfies the initial data. In that case, the elastic discontinuities led to that stress state so that the characteristic structure corresponds to that depicted in figure 5.15a.

(1-b) if on the other hand the point reached on the initial yield does not satisfy the initial stress σ 22 , a fast wave is added in order to browse the initial yield surface until the initial data is recovered. This situation is depicted in figure 5.15b.

(2-a) the curve resulting from the reverse integration across a slow wave intersects the plane σ 12 = 0. Then, assuming that the paths of slow waves are symmetric with respect to that plane, a fast wave is added in order to reach the initial yield surface at the initial value of σ 22 . Indeed, the fast waves have been shown to yield horizontal paths in the (σ 11 , σ 12 ) plane, in such a way that only that type of wave enables the achievement of the initial elastic domain. This also corresponds to figure 5.15b.

(2-b) if at last, the guessed state is such that σ 12 = 0, a fast wave allows reaching the initial yield surface as depicted in figure 5.15c. Notice however that the above elementary loading paths are based on strong assumptions about the symmetry of the loading paths that have not been shown so far. As a result, additional work must be performed in order to develop this approach and to introduce it in numerical methods. Moreover, the hardening of the material may modify the behavior of the loading paths and have not been considered yet. At last, the generalization of the approach followed in this chapter to more complex hardening models (kinematic, nonlinear etc.) and other yield surfaces would be very interesting for the understanding of the physics.

Chapter 6

Conclusion and future works

General conclusion

The purpose of this work was the development of a numerical method allowing the accurate solution of hyperbolic problems in solid mechanics. This class of mathematical problems governs the propagation of waves in finite deforming media. Given the ability to deal with large deformations enabled by mesh-free methods and the possibility of describing discontinuous solutions provided by the DG approximation, it has been proposed to merge those two approaches. Furthermore, particular attention has been paid throughout this manuscript to the characteristic structure of the solution of hyperbolic problems in order to: (i) provide qualitative and quantitative comparisons between numerical solutions and exact ones; (ii) allow numerical schemes to mimic the behavior of physical systems.

In chapter 2 the balance laws and constitutive equations of solid mechanics have been derived and written as a system of conservation laws in conservative form. The characteristic analysis of the associated quasi-linear form then leads to the formulation of a hyperbolicity condition of the system in terms of the eigenstructure of the acoustic tensor.

Exact solutions of one-dimensional problems have then been presented. First, the existing solutions of Riemann problems in elastic and elastoplastic solids with linear hardening under small strains emphasized that discontinuous waves may propagate. Second, the derivation of exact solutions of Picard problems in a one-dimensional Saint-Venant-Kirchhoff hyperelastic medium highlighted the existence of shocks and rarefaction waves in solid dynamics. Given the complexity hyperbolic systems may present, especially for more than one space dimension or other constitutive models, a well-known approximate-state Riemann solver has been presented. Such solvers represent an ideal way to take into account the characteristic structure of the solutions within numerical schemes and are therefore widely used in finite volume.

Eulerian and Lagrangian formulations of the Material Point Method have been recalled in chapter 3. This method constitutes the starting point of the present work due to its ability to avoid the mesh instabilities caused by large deformations often met in solid mechanics. Nevertheless, the oscillating solutions provided by this numerical scheme prevent accurately following the waves. The aforementioned numerical noise can be eliminated by going back to the early version of the approach: the Particle-In-Cell method, at the cost of additional numerical diffusion. Thus, the approach proposed consisted of: (i) removing the oscillations by re-introducing within the MPM the diffusive mapping between material points and grid nodes used in PIC; (ii) reducing the diffusion thus introduced by means of the discontinuous Galerkin approximation.

The Discontinuous Galerkin Material Point Method has therefore been derived with a total Lagrangian formulation. The method inherits appealing features of the finite element, finite volume and material point methods. First, it is based on the weak form of a system of conservation laws so that the same order of accuracy is achieved for both velocity and gradients. Second, numerical fluxes defined at element interfaces arising in the weak form are computed with the solution of Riemann problems by means of the approximate-state Riemann solver presented in chapter 2. This approach not only allows the introduction of the characteristic structure in the numerical solution, but it also enables taking into account transverse corrections by adapting the Corner Transport Upwind method. Third, the DGMPM solution scheme can be combined with ODE integrators for source terms if fractional-step methods are considered. The last three points are common with FVM. Then, the FEM approximation using polynomial shape functions is used, thus providing a convenient framework for high-order approximation. In particular, as for other DG-methods, the order of approximation may be modified element-wise. At last, the arbitrary grid of the MPM, in which the particles move while carrying fields and internal variables, enables the employment of mesh-adaption strategies without the need for additional diffusive projection steps. Notice however that the last aspect is not considered in this thesis.

Then, the numerical analysis of DGMPM schemes showed that the stability properties of the method may be better than those of the original MPM depending on the space discretization used. Namely, if the DGMPM semi-discrete system is discretized in time with the forward Euler algorithm, the optimal Courant number can be achieved when one particle lies in every cell of a one-dimensional grid. The same result holds if the second-order Runge-Kutta time discretization is used, and even for more than one material point per cell depending on their positions. Even though the condition CFL=1 is not valid in general, a formula for evaluating the critical Courant number allows the ensuring of the stability of the scheme, whereas it is not possible for the MPM. Analogously, the two-dimensional DGMPM scheme exhibits the optimal stability condition providing that one particle per cell is used along with the CTU algorithm for transverse corrections. It has been also shown that the stability features of the two-dimensional scheme depend on the ratio of the characteristic speeds involved in the linear advection equation. On the other hand, the convergence analysis of the method on a one-dimensional linear elasticity problem under small strains showed that by using first-order shape functions, only first-order accuracy can be achieved for velocity and gradients. Although one could expect second-order convergence since the DGMPM approximation is similar to that of FEM, the PIC mapping of the updated nodal velocity from the grid to particles used yields first-order accuracy.

Chapter 4 was devoted to the illustration of the DGMPM on simulations of solid dynamics problems. Comparisons with other numerical methods and exact solutions showed very good behavior. More specifically, since the optimal CFL number can be achieved, the DGMPM is able to capture the discontinuities arising in the Riemann problem in an elastic bar within the linearized geometrical framework. On the other hand, it has been shown that when this condition no longer holds, that is CFL < 1, the method suffers less from diffusion than the original MPM using the PIC mapping. Furthermore, the same behavior has been observed for a twodimensional problem in an elastic medium. The DGMPM moreover provides good results on problems involving large deformations in one and two-dimensional solids. Hence, the method fulfills its objectives.

Then, history-dependent constitutive models have been considered within the infinitesimal theory. First, the solutions of a plane wave in an elasto-viscoplastic solid resulting from DGMPM, MPM, FVM and FEM have been compared. In the stiff limit, that is when the viscosity is close to zero, the above solutions should tend to the elastic-plastic one. The simulations based on fractional-step methods combined with the DGMPM highlighted that the use of suitable ODE integrators is required in order to compute the correct solution and the ability to achieve high-order accuracy. Alternatively, dedicated approximate Riemann solvers can be employed within the DGMPM for elastoplasticity so as to take into account elastic as well as plastic characteristics and hence, to get more accurate solutions. Nevertheless, the structure of the solutions of hyperbolic problems in elastic-plastic materials is only known for one-dimensional solids in such a way that approximate Riemann solvers do not exist for multi-dimensional cases.

The simulations in history-dependent materials performed in chapter 4, as well as the point of view adopted in this thesis, motivated the developments presented in chapter 5. Indeed, the need for a better understanding of the characteristic structure of elastic-plastic hyperbolic problems arises from the lack of approximate Riemann solvers in more than one-dimension. Although some references have proposed the spectral analysis of hyperbolic systems in elastoplastic solids under small strains since the 60s, the problems treated were limited to particular cases of plane stress and plane strain conditions. The present approach consisted in rewriting the existing formulations by means of the elastoplastic stiffnesses to provide a unified framework for all twodimensional problems in the Cartesian coordinates system. It is well-known that the particularity of elastoplastic problems in more than one space dimension lies in the fact that the wave pattern of the solution time discretization could be embedded in the method in order to reduce the numerical diffusion when more than one particle lies in a cell, providing that a von-Neumann analysis of two-dimensional DGMPM schemes is performed.

At last, some properties of the method which have not be considered so far must be investigated, that is, the conservation of angular momentum and the performances for bending dominated problems in the nearlyincompressible limit. The DGMPM formulation might then require the use of an alternative projection of the nodal solution that conserves the angular momentum, and the introduction of the Jacobian J as a conserved variable to solve bending dominated nearly-incompressible deformations [START_REF] Bonet | A first order hyperbolic framework for large strain computational solid dynamics. part i: Total lagrangian isothermal elasticity[END_REF].

(ii) On the solution of hyperbolic systems in elastoplastic solids

The study of the characteristic structure of hyperbolic systems in elastic-plastic solids allows perspectives.

First, the iterative procedure inspired by that of Lin and Ballman [START_REF] Lin | A riemann solver and a second-order godunov method for elastic-plastic wave propagation in solids[END_REF] and proposed in chapter 5 has to be implemented. Moreover, the analytical developments proposed here can be prosecuted in order to get additional information about the solution so that an approximate Riemann solver can be developed.

Second, the framework provided allows to take into account more general hardening laws that is, kinematic, combined isotropic and kinematic, and nonlinear hardenings. Nevertheless, the resulting equations should be much more complex and prevent the building of iterative or approximate Riemann solvers.

At last, plastic shocks, which raised scientific questions that are still unanswered, have to be studied for they are involved in problems of dynamic impacts in elastic-plastic solids.

1. 3 .

 3 APPROACH FOLLOWED IN THIS WORK

CHAPTER 2 .

 2 HYPERBOLIC PARTIAL DIFFERENTIAL EQUATIONS FOR SOLID DYNAMICSor equivalently, in matrix form:

Figure 2 .Figure 2 . 1 :

 221 Figure 2.1: Example of initial curve C in the (x, t, u) space and its projection C 0 in the (x, t) plane.

Figure 2 . 2 :

 22 Figure 2.2: Examples of integral surfaces passing through the same curve C defined such that t = constant and u = constant along C .

Figure 2 . 3 :

 23 Figure 2.3: Family of characteristic curves corresponding to the eigenvalues of the first-order systems given in examples 1 and 2.

Figure 2 . 6 :

 26 Figure 2.6: Solution to linear advection equation of the quantity P i with characteristic speed c i .

Figure 2 . 7 :

 27 Figure 2.7: Solution to Riemann problem (2.62) for an elastic bar.

Figure 2 . 8 :

 28 Figure 2.8: Example of a solution of a Riemann problem in a homogeneous elastoplastic bar with linear hardening and initial plastic strain ε p (x, 0) = 0.

Figure 2 . 9 :

 29 Figure 2.9: Solutions of the Picard problem (2.72) depending on initial and boundary data.

Figure 2 . 10 :

 210 Figure 2.10: General wave patterns arising in the solution of the Riemann problem (2.84) depending on initial data. (a): 1-shock-2-rarefaction. (b): 1-rarefaction-2-shock.

Figure 2 . 11 :

 211 Figure 2.11: Set of connected states Q to initial data through shock and rarefaction waves with v L = v R = 0 in both cases: (a) 1-shock, 2-rarefaction solution ; (b) 1-rarefaction, 2-shock.

Figure 2 . 12 :

 212 Figure 2.12: Comparison of neo-Hookean and Saint-Venant-Kirchhoff hyperelastic models.

Figure 3 . 1 :

 31 Figure 3.1: Representation of a continuum body by a set of material points in R 2 .

  (a) Computation of the consistent and lumped mass matrices as well as external forces (Neumann boundary conditions) from equations (3.6a) and (3.6c). (b) Projection of fields to the grid (3.8) and enforce Dirichlet boundary conditions on the nodal velocity.

  time t = 1.69 × 10 -4 s. time t = 4.23 × 10 -4 s.

Figure 3 . 2 :

 32 Figure 3.2: MPM solutions of the bars impact problem for various discretizations. (a)-(b) comparison of stress computed with USF and USL formulations and exact solutions (c) numerical total energy evolution. Courant number: CFL= 0.7 ; v 0 = 1 200 E ρ .

  time t = 1.69 × 10 -4 s. time t = 4.23 × 10 -4 s.

Figure 3 . 3 :

 33 Figure 3.3: Comparison between exact, USF and USL velocities of the bars impact problem for various discretizations.

  v(m/s) (a) time t = 1.69 × 10 -4 s. v(m/s) (b) time t = 4.23 × 10 -4 s.

Figure 3 . 4 :

 34 Figure 3.4: Comparison between USL algorithm with classical or PIC mapping and exact velocities of the bars impact problem for various discretizations.

Figure 3 . 5 :

 35 Figure 3.5: Comparison of USL solutions using either PIC or FLIP projection and exact solution of the bars impact problem for various discretizations. (a)-(b) comparison of stress profiles at various times (c) evolution of total energy. CFL=0.7 ; v 0 = 1 200 E ρ .

2 Figure 3 . 6 :

 236 Figure 3.6: Duplication of nodes at an interface and building of initial conditions of the Riemann problem (2D).

Figure 3 . 7 :

 37 Figure 3.7: Normal and transverse fluctuations defined from edge i.normal vector is N i . The Riemann problem defined at this edge gives rise to normal fluctuations A i,-(∆U) and A i,+ (∆U) contributing to cells L and R respectively. These terms lead to the computation of transverse fluctuations giving contribution to neighboring cells across edges (j) and (k) for cell L, and across edges (m) and (l) for cell R. Transverse fluctuations are computed by projecting normal fluctuations onto the characteristic basis associated with the Riemann problem (3.28) defined on the adjacent edge, hence the name transverse Riemann solver. The spectral analysis of the corresponding Jacobian matrix, carried out in[START_REF] Kluth | Analyse mathématique et numérique de systèmes hyperélastiques et introduction de la plasticité[END_REF], leads to right eigenvectors that are be written R i W hereinafter. The negative normal fluctuation is, for instance, decomposed onto the characteristic basis associated with edge (j) as:

Figure 3 . 8 :

 38 Figure 3.8: Corner ghost nodes in a two-dimensional DGMPM mesh.

  (a) The lumped mass and pseudo-stiffness matrices M L ij and K α ij of the semi-discrete form (3.21) are computed. (b) The discrete equation (3.22) and Riemann problems at cell interfaces (3.29) require a projection of fields onto the grid:

  (c) The specific flux vectors Fi α involved in the equations are computed from Q i knowing ρ 0 , thus avoiding the computation of constitutive equations. (d) Enforce the boundary conditions on ghost nodes. (e) Computation of interface fluxes: 1-Build the state vectors Q X ± N based on Q i,n where i denotes the nodes belonging to the face on both sides of an interface. 2-Compute the stationary solution Q * by means of the approximate-state Riemann solver. 3-Calculate the corresponding Godunov flux F N (Q * ) by either using the DCU or the CTU approach. (f) Advance the solution in time by solving the discrete equation (3.22) at each node.

1 +

 1 S 2c,µ qn 2c are written for simplicity f n c , and Courant number a∆t/∆X arises. Introduction of these equations in the updated material point solution (3.33) leads after some simplifications to:

Figure 3 . 10 :

 310 Figure 3.10: Evolution of the L 2 errors in stress with regard to the number of particles lying in every cell.

Figure 3 . 11 :

 311 Figure 3.11: Evolution of the L 2 error in velocity with regard to the number of particles lying in every cell.

50 Figure 3 . 12 :

 50312 Figure 3.12: Evolution of the L 2 errors in stress and velocity for two particles per cell. CFL set to 0.1.

Figure 3 . 13 :

 313 Figure 3.13: Parent and current configurations of a rectangular four-node bilinear element in figure 3.13, the current location X of an arbitrary point in cell C maps to the parent coordinates (ξ, η) in the domain [-1, 1] × [-1, 1] according to:

59 )Figure 3 .

 593 Figure 3.14 shows transverse corrections in the cell C based on fluctuations coming from Bottom (B), Left (L), Bottom Left (BL) neighbor elements. The use of the numbering of interfaces and nodes adopted in figure 3.14

Figure 3 . 14 :

 314 Figure 3.14: Two-dimensional patch of cells of constant size ∆X × ∆Y . allows the specialization of equation (3.59) to intercell fluxes of cell C:

  at time step n and node i of cell C. Denoting the number of particles in cell C and the mass they carry by N C p and m C respectively, the mass density reads ρ = N C p m C

Table 3 . 3 :

 33 Values of critical Courant number a ∆t ∆X for two-dimensional DGMPM scheme using either DCU or CTU with respect to the locations of the 2 material points lying in every cell as a function of the speed ratio a/b.

Table 3 . 4 :

 34 Values of critical Courant number a ∆t ∆X for two-dimensional DGMPM scheme using either DCU or CTU with respect to the material points distribution as a function of the speed ratio a/b.

Figure 4

 4 time t = 1.21 × 10 -4 s. •10 8 (b) time t = 4.84 × 10 -4 s.

Figure 4 . 1 :

 41 Figure 4.1: Stress and velocity solutions of the Riemann problem in an isotropic elastic bar at two different times (columns a and b). Comparison between DGMPM coupled with Euler or RK2 time integration, MPM-USL formulation using either PIC or FLIP mapping, and the exact solution for an initial velocity set to v 0 = c 2000 .

Figure 4 . 2 :

 42 Figure 4.2: Evolution of total energy e for DGMPM and MPM-USL solutions on the Riemann problem in an elastic bar.

  t = 4.17 × 10 -4 s.•10 9(b) t = 9.38 × 10 -4 s.

Figure 4 . 3 :

 43 Figure 4.3: Plastic strain and longitudinal stress resulting from MPM and DGMPM simulations before (column a) and after (column b) reflection of incident plane waves at the free boundaries. Non-stiff problem: τ = 50∆t.

Figure 4 . 4 :

 44 Figure 4.4: Plastic strain and longitudinal stress resulting from MPM and DGMPM simulations before (column a) and after (column b) reflection of incident elastic waves at the free boundaries. Stiff problem: τ = ∆t × 10 -2 .

Figure 4 . 5 :

 45 Figure 4.5: Plane wave solution of the Riemann problem in an elastoplastic material with linear hardening. Comparison between DGMPM-Euler using either an elastic or an elastoplastic Riemann solver, MPM, and exact solutions in terms of stress and plastic strain.

Figure 4 . 6 :

 46 Figure 4.6: Stress and plastic strain solutions of the Riemann problem in a one-dimensional elastoplastic medium: comparison between FEM, Superbee FVM and DGMPM-Euler using an elastoplastic Riemann solver. All CF L set to 1.

2 Figure 4 . 7 :

 247 Figure 4.7: Geometry, loading and boundary conditions for the tensile impact problem on a two-dimensional elastic medium.

Figure 4 Figure 4 . 8 :

 448 Figure 4.8: Isovalues of longitudinal stress σ 11 solution of the tensile impact problem in a two-dimensional elastic medium. Comparison between FEM (CFL=0.9), DGMPM-CTU using 1ppc (CFL=1) or 4ppc (CFL=0.23), and MPM using 1ppc (CFL=0.7). the longitudinal stress σ 11 in the two-dimensional medium at two different times with the traction force set to σ d = 200 M P a. The two instants at which the solutions are depicted correspond to incident and reflected

Figure 4 . 9 :

 49 Figure 4.9: Evolution of longitudinal stress σ 11 along the bottom boundary of the elastic square plate. Comparison between FEM (CFL=0.9), DGMPM-CTU using 1ppc (CFL=1) or 4ppc (CFL=0.23), and MPM using 1ppc (CFL=0.7).

Figure 4 . 10 :

 410 Figure 4.10: Evolution of longitudinal stress σ 11 and plastic strain ε p11 along the bottom boundary of the elastic-plastic square. Comparison between FEM (CFL=0.9), DGMPM using 1ppc (CFL=1) and DGMPM using 4ppc (CFL=0.23) solutions.

Figure

  Figure 4.11 furthermore shows FEM and DGMPM solutions in terms of equivalent plastic strain p, before and after reflection of the longitudinal pressure wave on the right end of the domain. MPM results for the same time steps are plotted separately in figure 4.12 due to the large amplitude of equivalent plastic strain.A concentration of plastic strain occurs at the interface between loaded and traction free parts of the left boundary of the domain in DGMPM solutions. Such concentrations in the high gradient region can also be seen in MPM solutions depicted in figure4.12 so that this phenomenon seems to be a singularity owed to boundary conditions. Although the final DGMPM and FEM profiles of equivalent plastic strain have the same shape, the results highlight the fact that different numerical methods yield different assessments of irreversible

84

 84

Figure 4 . 11 :Figure 4 . 12 :

 411412 Figure 4.11: Isovalues of equivalent plastic strain p in an elastic-plastic plate linear isotropic hardening material at two different times. Comparison between FEM (CFL=0.9), DGMPM using 1ppc (CFL=1) and DGMPM using 4ppc (CFL=0.23) solutions.(a) time t = 3.5 × 10 -4 s

Figure 4 . 13 :

 413 Figure 4.13: First Piola-Kirchhoff stress along a horizontal line of the one-dimensional hyperelastic medium at different times: solution of the compressive impact problem Π d = 4× 10 8 P a. Comparison between updated and total Lagrangian MPMs (CFL=0.5), DGMPM-Euler with 1ppc (CFL=1) or 2ppc (CFL=0.5), DGMPM-RK2 with 2ppc (CFL=1), and the exact solution.

Figure 4 . 14 :

 414 Figure 4.14: First Piola-Kirchhoff stress along a horizontal line of the one-dimensional hyperelastic medium at different times: solution of the compressive impact problem Π d = 2 × 10 10 P a. Comparison between MPM (CFL=0.5), DGMPM-Euler with 1ppc (CFL=1) or 2ppc (CFL=0.5), DGMPM-RK2 with 2ppc (CFL=1), and the exact solution.

Figure 4 . 15 :

 415 Figure 4.15: First Piola-Kirchhoff stress solution of the tensile impact problem at different times: Π d = -2 × 10 10 P a. Comparison between MPM (CFL=0.5), DGMPM-Euler with 1ppc (CFL=1) or 2ppc (CFL=0.5), DGMPM-RK2 with 2ppc (CFL=1), and the exact solution.

  t = 4.18 × 10 -4 s.

Figure 4 . 16 :

 416 Figure 4.16: First Piola-Kirchhoff stress solution of the tensile impact problem at different times: Π d = -2×10 9 P a. Comparison between MPM (CFL=0.5), DGMPM-Euler with 1ppc (CFL=1) or 2ppc (CFL=0.5), DGMPM-RK2 with 2ppc (CFL=1), and the exact solution.

  dgmpm (1ppc) (a) t = 9.94 × 10 -4 s dgmpm (1pc) (b) t = 2.54 × 10 -4 s

Figure 4 . 18 :

 418 Figure 4.18: Isovalues of the longitudinal PK1 stress Π 11 solution of the plane wave problem in a two-dimensional compressible neo-Hookean material with Riemann-type data on the initial velocity. Comparison between DGMPM-DCU with 1ppc (CFL=0.57) or 4ppc (CFL=0.23) and MPM solutions also using 1ppc and 4ppc (CFL=0.7).

Figure 4 . 19 :

 419 Figure 4.19: Evolution of total energy for the plane wave problem in a two-dimensional compressible neo-Hookean material.

90

 90

Figure 4 . 20 :Figure 4 . 21 :

 420421 Figure 4.20: Isovalues of Cauchy stress tensor component σ 11 in a two-dimensional plate made of a neo-Hookean material, submitted to a velocity v • e 1 = -1000 m/s on a part of its left end.

Figure 4 . 22 :Figure 4 . 23 :

 422423 Figure 4.22: Isovalues of velocity component v 1 in a two-dimensional plate made of a neo-Hookean material, submitted to a velocity v • e 1 = -1000 m/s on a part of its left end.

  .2c) to σ 33 yields: σ33 = C ep 33ij εij = 0 from which one writes: C ep 3333 ε33 = -C ep 33ij εij i, j = {1, 2} Hence, the constitutive equations are rewritten by means of a two-dimensional tangent modulus C ep : σij = C ep ijkl εkl -C ep ij33 C ep 33kl C ep 3333 εkl = C ep ijkl εkl i, j, k, l = {1, 2} (5.8)

  0 along dx/dt = c s (5.14c) -ρc s l 2 • dvl 2 i n j dσ ij = 0 along dx/dt = -c s (5.14d) α 11 dσ 11 + α 12 dσ 12 + α 22 dσ 22 = 0 along dx/dt = 0 (5.14e)

Figure 5 . 1 :

 51 Figure 5.1: The method of characteristics through slow and fast simple waves in the (x n , t) plane.

5. 3 .= e 1 n = e 2 n = e 1 n = e 2 dv 1 = -dσ11 ρc 2 s

 3121212 LOADING PATHS THROUGH SIMPLE WAVESRight-going slow waveRight-going fast wave n

3 Figure 5 . 2 :

 352 Figure 5.2: Stress path followed in a fast simple wave for the thin-walled tube problem.Comparison between the results obtained from equations (5.31) and these of[START_REF] Clifton | An analysis of combined longitudinal and torsional plastic waves in a thin-walled tube[END_REF].

3 Figure 5 . 3 :

 353 Figure 5.3: Stress paths followed in a slow simple wave for the thin-walled tube problem.Comparison between the results obtained from equations (5.31) (cross markers) and those of[START_REF] Clifton | An analysis of combined longitudinal and torsional plastic waves in a thin-walled tube[END_REF] (solid lines).

Figure 5 . 5 :

 55 Figure 5.5: Loading paths through a fast simple wave starting from the initial yield surface with initial condition σ 11 = σ 22 = 0 in directions of decreasing and increasing σ 11 .

Figure 5 . 6 :

 56 Figure 5.6: Stress paths in a slow simple wave for various starting points lying on the initial yield surface for σ 22 = -5.8×10 7 P a. Projections in the stress space (figure a) and deviatoric plane (figure b).

Figure 5 . 7 :

 57 Figure 5.7: Stress paths in a slow simple wave for various starting points lying on the initial yield surface for σ 22 = 0. Projections in the stress space (figure a) and deviatoric plane (figure b).

109 5 . 4 .Figure 5 . 8 :

 5458 Figure 5.8: Stress paths in a slow simple wave for various starting points lying on the initial yield surface for σ 22 = 5.8 × 10 7 P a. Projections in the stress space (figure a) and deviatoric plane (figure b).

3 eFigure 5 . 9 :

 359 Figure 5.9: Loading paths through a fast simple wave with initial condition σ 22 = 0 for different starting points on the initial yield surface.

Figure 5 . 10 :

 510 Figure 5.10: Loading paths through slow simple waves for different starting points on the initial yield surface for the initial condition σ 22 = -1.3 × 10 8 P a. curves depicted in figure 5.10 results from the negative initial value σ 22 = -1.3 × 10 8 P a for several starting

Figure 5 . 11 :

 511 Figure 5.11: Loading paths through slow simple waves for different starting points on the initial yield surface for the initial condition σ 22 = 0.

Figure 5 . 12 :

 512 Figure 5.12: Loading paths through slow simple waves for different starting points on the initial yield surface for the initial condition σ 22 = 1.3 × 10 8 P a.

3 Figure 5 . 13 :

 3513 Figure 5.13: Fast simple wave solutions of the plane strain problem C = 1 × 10 10 P a.

s 1 s 2 s 3 (a) σ 22 = - 1 . 3 × 10 8 3 (c) σ 22 = 1 . 3 ×

 32213832213 10 8 P a.

Figure 5 . 14 :

 514 Figure 5.14: Slow simple wave solutions of the previous problem with C = 1 × 10 10 P a.

c 2 Figure 5 . 15 :

 2515 Figure 5.15: Characteristic structures possibly occurring in two-dimensional elastic-plastic solids.

Table 3 .

 3 2: Order of accuracy of MPM and DGMPM with regard to the number of particles per cell velocity than stress owing to the use of shape function derivatives for the computation of gradients. Hence, the MPM rate of convergence in terms of velocity is getting closer to two.

	PPC DGMPM-Euler DGMPM-RK2	MPM	MPM-PIC
		σ	v	σ	v	σ	v	σ	v
	2 6 10 20	0.63 0.66 0.67 0.68	0.63 0.66 0.67 0.67	0.80 0.91 0.93 0.95	0.80 0.91 0.93 0.95	0.88 1.46 0.94 0.85 0.91 1.62 1.07 0.92 0.92 1.62 1.10 0.92 0.92 1.61 1.12 0.93

Table 4 .

 4 1: Material parameters. The viscosity is expressed as a function of the relaxation time τ characterizing relaxation systems (see section 2.2.4).

Table 5 .

 5 1: Summary of the ODEs satisfied inside right-going slow and fast simple waves.

Acknowledgments

Chapter 3

The Discontinuous Galerkin Material Point Method Chapter 5

Contribution to the solution of elastic-plastic hyperbolic problems in two space dimensions

According to the ODEs of table 5.1, the particular values of the loading functions ψ s,f i through the simple waves propagating in direction e i for i = {1, 2}, provide information about the loading paths in stress space. First, ψ s,f i = 0 leads to dσ ii = 0 (no sum on i) so that the longitudinal stress is constant within the simple wave. Conversely, with loading functions tending to infinity, the stress σ 12 does not vary. Notice that the coefficients α ij of the left eigenvector of the Jacobian matrix associated with the zero eigenvalue (5.13) also have to be regarded. Nevertheless, those terms resulting from products of the components of the elastoplastic tangent modulus have complex expressions and are assumed to have non-zero values in the remainder of the manuscript.

The above discussions are now specified to plane strain and plane stress, for which loading conditions leading to A 12 = 0 and A 11 -A 22 = 0 are identified.

The plane strain case

The case of plane strain is first considered by using the elastoplastic tangent modulus so that the components of the acoustic tensor for n = e 1 read:

(5.22a)

(5.22b)

The associated eigenvalues are then: Condition s 12 = 0: According to equations (5.23a) and (5.23b), the eigenvalues of the acoustic tensor become:

Condition s 11 = 0: Considering the relation (5.7) between stress components for plane strain, one writes:

so that s 11 = 0 is equivalent to:

In contrast to what has been seen previously, the functions ψ cannot be undetermined in the case s 11 = 0 since the equation A ep 11 -A ep 22 = λ + µ + βs 2 12 = 0 does not admit real solutions. However, the stress state (5.24) yields the following characteristic speeds:

ρc 2 f = λ + 2µ so that the celerity of fast waves is equivalent to that of elastic pressure waves under plane strain c f = (λ + 2µ)/ρ = c 1 .

The same analysis can be carried out in the direction n = e 2 by considering the following acoustic tensor components:

The characteristic speeds are then: With these expressions, the same remarks as for n = e 1 can obviously be made by replacing s 11 with s 22 . Among the above results, the most significant arises from the condition s 12 = 0. Indeed, it has been seen that A ep 12 = 0 leads to ψ s 1 = 0 and ψ s 2 → ∞ in such a way that the corresponding loading paths in the (σ 11 , σ 12 ) plane are respectively vertical and horizontal. Under the orthogonality property of the loading functions, the stress path followed in a fast wave propagating in the direction e 1 is horizontal in the same plane. Hence, if the path through a fast wave intersects the plane σ 12 = 0, the shear stress component remains constant afterwards. The same result holds for the slow wave propagating in the direction e 2 . The above conclusion are summarized in table 5.2. depends on the external loads. More precisely, one must be able to connect the stationary stress state of a Picard problem to the initial data through loading paths followed inside fast and slow simple waves of combined-stress.

Stress path in (σ

Hence, the mathematical study of the loading paths resulting from the characteristic analysis of the hyperbolic system, has been performed. This analysis resulted in the identification of typical behaviors of these stress paths. Nevertheless, the complexity of the equations prevents the complete characterization of the stress evolution inside the simple waves. As a consequence, a numerical study of the loading paths based on the implicit integration of the ODEs governing the evolution of the system within the waves has been carried out for plane strain and plane stress problems. First, the particularization of the equations to those governing the thin-walled tube problem treated in the literature, allowed to validate the present approach. Second, the loading paths followed inside fast and slow waves for problems under plane strain and plane stress have been plotted so that additional features have been highlighted. Finally, some ideas in order to develop an iterative Riemann solver have been drawn based on the identified behaviors of the loading paths. This aspect of the present work aimed at building numerical schemes that are based on a robust discretization, and embed a sufficient amount of information in order to mimic the physical response. Though interesting features of the loading paths have been emphasized, additional effort is required so as to introduce the characteristic structure by means of an approximate Riemann solver.

Future work

This thesis has consisted of: (i) the foundation of a computational framework for the simulation of hyperbolic problems in finite deforming solids which enables several extensions; (ii) the study of the characteristic structure of the solution of elastoplastic hyperbolic problems in two space dimensions. Each of these aspects comes with medium/long term perspectives.

(i) Improving the DGMPM Although the DGMPM has been developed in order to capture discontinuities so that its first-order accuracy is sufficient, the method should be adapted so as to achieve higher-order convergence in regions where the solution is smooth. Higher-order accuracy should be achieved by using shape functions based on high-order polynomials (padaption). Notice furthermore that the polynomial order can be set arbitrarily within the background mesh due to the discontinuous Galerkin approximation. Nevertheless, the particle-based quadrature of the DGMPM weak form may limit the polynomial order for a given space discretization in order to avoid problems due to reduced integration. Indeed, it has been seen that the diagonal lumping of the mass matrix can balance the reduced integration with one particle and linear shape functions, but nothing can be stated for higher-order polynomials. Hence, analogously to recent developments of the MPM [START_REF] Tielen | High-order material point method[END_REF][START_REF] Tielen | A high order material point method[END_REF][START_REF] Ming | Improving the Material Point Method[END_REF], functions reconstruction techniques such as moving least squares or spline interpolation may be employed in order to use Gauss quadrature even for few particles in a cell. However, such an implementation would be equivalent to the DGFEM which suffers from a very restrictive CFL condition that prevents the capturing of discontinuities, which is prohibitive. One solution to circumvent this issue would be to combine the p-adaption with h-adaption by refining the arbitrary grid in the vicinity of discontinuities. The underlying idea is to use a thin layer of elements containing one particle along the discontinuity, for which the CFL number can be set to unity. The thickness h of that layer would be set so that the time step ∆t s corresponding to the most restrictive Courant number in the smooth regions CFL s , allows the capture of the discontinuity, that is: ∆t s = CFL s h c , where c is the speed of the fastest wave. Alternatively, one can imagine a particle-adaption, while ensuring the conservation of mass, so that the particle-based integration can be adapted to the regularity of the solution. That is: remove material points in such a way that one particle per cell is used near a discontinuity; add particles in elements in which the polynomial order requires so. Furthermore, one could benefit from Roe's theorem 4 in order to adapt the distribution of particles so that a given order of accuracy is achieved for a given polynomial order, as suggested in section 3.3. Nonetheless, it is not guaranteed that such a configuration can be found.

Slope or wave limiters have not been considered so far within the DGMPM due to the non-oscillatory solutions provided by the scheme. Nevertheless, the extension of the method to higher-order approximation may require the introduction of such numerical tools that are well-known in DG-methods [START_REF] Cockburn | Discontinuous galerkin methods for convection-dominated problems[END_REF].

The second-order Runge-Kutta time integrator has not been considered for two-dimensional problems, due to the lack of a formula to determine the critical CFL number for a given space discretization. Such a two-stages Titre : Extension de la Méthode des Points Matériels à l'approximation de Galerkin Discontinue : Application aux problèmes hyperboliques en mécanique des solides. Abstract : In this thesis, the material point method (MPM) is extended to the discontinuous Galerkin approximation (DG) and applied to hyperbolic problems in solid mechanics. The resulting method (DGMPM) aims at accurately following waves in finite-deforming solids whose constitutive models may depend on the loading history. Merging finite volumes and finite elements methods, the DGMPM takes advantage of an arbitrary computational grid in which fluxes are evaluated at element faces by means of approximate Riemann solvers. This class of solvers enables the introduction of the characteristic structure of the solutions of hyperbolic partial differential equations within the numerical scheme.

Convergence and stability analyses, along with one and two-dimensional numerical simulations, demonstrate that this approach enhances the MPM ability to track waves.

On the other hand, a second purpose has been followed: it consists in identifying the response of two-dimensional elastoplastic solids to dynamic step-loadings in order to improve numerical results on these problems. Although some studies investigated similar questions, only particular cases have been treated. Thus, a generic framework for the study of the propagation of simple waves in elastic-plastic solids under plane stress and plane strain problems is proposed in this thesis. The loading paths followed inside those simple waves are further analyzed.