Saber Amdouni

Éric Dalissier

Mathieu Fabre

Amin Karoui

Khadidja Ben

Imen Mekkaoui

Rafik Ben

Rabii Mlika

Keywords:

For the guardians of the kurdish people, the beauty of our lands.

For the courageous mountains of Kurdistan Pour les gardiens du peuple kurde, la beauté de nos terres.

Pour les courageuses montagnes du Kurdistan

Je dédie cette thèse à la mémoire de mon grand-père dont j'aurais tellement aimé qu'il soit présent mais la nature en a décidé autrement.

Je souhaite dédier en particulier cette thèse à mon père, ma mère et tous les autres membres de ma famille. Malgré la distance qui nous sépare, je n'ai jamais manqué de soutien et d'encouragements de leurs part.

Introduction

Un ou plusieurs corps déformables en mouvement peuvent entrer en contact en un point, le long d'une ligne ou d'une surface. Les systèmes aux dérivées partielles sont particulièrement bien adaptés pour modéliser des mouvements complexes de corps soumis à des conditions aux limites faisant intervenir les conditions de contact. Ces problèmes de contact occupent un rôle prépondérant dans de nombreux domaines comme par exemple dans l'industrie automobile où le contact entre le pneu et la chaussée (cf. Figure 1) ainsi que les crash tests (cf. Figure 2) intéressent particulièrement les ingénieurs ou encore en génie civil où les ouvrages d'arts de plus en plus complexes tels que les ponts ou encore les gratte ciel nécessitent une étude préalable des matériaux en contact qui vont devoir subir d'importants efforts.

Figure 1: La simulation de roulement du pneu par Brinkmeier et al [START_REF] Brinkmeier | A finite element approach for the simulation of tire rolling noise[END_REF] La modélisation des phénomènes de contact pose de sérieuses difficultés du point de vue mathématique. Plus précisément, suite à l'application d'une charge sur deux corps en contact, la surface de contact devient inconnue a priori. Ainsi, la modélisation des problèmes de contact conduit à des équations non-linéaires équivalentes à un système d'inéquations à faible régularité Présentation de la thèse de la solution. En présence de frottement et dans le cas multi-corps, le problème de contact devient plus compliqué. La formulation d'un problème bien posé présente un défi majeur pour le théoricien surtout dans le cas de dimensions supérieures à un et de problèmes complexes (frottement, multi-corps...). D'un point de vue numérique, l'établissement d'un schéma robuste et efficace reste encore un challenge important.

Un peu d'histoire

Certaines découvertes archéologiques faites en Egypte montrent que les bâtisseurs des pyramides s'intéressaient déjà aux problèmes de contact et de frottement. Ainsi des textes expliquant partiellement certaines techniques employées pour transporter des pierres de plusieurs tonnes sur des centaines de métres ont été mis à jour récemment par des archéologues. Par ailleurs, il existe également un certain nombre de bas reliefs comme celui trouvé à Ninive (cf. Figure 3) représentant le transport de divers objets de grande taille et mettant en évidence que ces problèmes intéressent l'homme depuis la nuit des temps.

Les premiers travaux significatifs dans le domaine du contact et du frottement ont été effectués par L. De Vinci au XVIème siècle. Il a établi les deux lois fondamentales de frottement; la première indique que le frottement est proportionnel à la force normale et la seconde met en évidence que les zones de contact n'ont pas d'effet sur le frottement. G. Amontons a redécouvert les deux lois fondamentales de frottement au XVII siècle qui n'ont pas été publiées par Léonard de Vinci. Il faut attendre le XVIIIème siècle et L. Euler qui en dérivant l'angle au repos d'une masse sur un plan incliné fait la distinction entre les frottements statique et cinétique. L. Euler écrivit dans [START_REF] Euler | Sur la diminution de la résistance du frottement[END_REF]: "L'expérience nous ayant fait voir que la force du frottement est toujours égale à une certaine partie de la pression, dont un corps est pressé contre la surface, sur laquelle il se meut, de sorte que le frottement ne dépend ni de la grandeur de la base, dont le corps touche la surface, ni du degré de vitesse, il n'est pas diffcile de déterminer l'effet du frottement dans toutes sortes de machines par le moyen du calcul: vu que le frottement doit être regardé comme une force constante, qui est toujours directement contraire à la direction du mouvement, et qui agit dans une direction, qui passe par le plan de l'attouchement des corps qui se meuvent l'un sur l'autre". Ch.-A. de Coulomb a formulé mathématiquement les lois de Amontons dans la deuxième moitié du XVIIIème siécle. Il a ainsi déterminé l'ampleur du frottement lors de mouvements de glissement et il a étudié également l'influence de quatre facteurs principaux sur le frottement: la nature des matériaux et leurs revêtements de surface, l'étendue de la surface, la pression normale et la durée de temps pendant laquelle les surfaces sont restées en contact (temps de repos). De plus, il a pris en considération dans ses travaux l'influence de la vitesse de glissement, de la température et de l'humidité [START_REF] Coulomb | Théorie des machines simples, en ayant égard au frottement de leurs parties et la raideur des cordages[END_REF]. Au XIXème siècle, H. R. Hertz dans [START_REF] Hertz | Uber die beruhrung fester elastischer korper[END_REF] s'est intéressé à des surfaces sans frottement et des solides élastiques linéaires non-conformes. Cependant, seulement quelques problèmes de contact peuvent être résolus analytiquement, H. R. Hertz s'est essentiellement penché sur des formulations simples et à des modéles idéalistes qui ne sont pas suffisamment réguliers pour aboutir à une solution Présentation de la thèse dans les cas complexes. Néanmoins des contributions dans cette direction ont été obtenues dans [START_REF] Galin | Contact problems in the theory of elasticity[END_REF][START_REF] Timoshenko | Theory of Elasticity[END_REF][START_REF] Gladwell | Contact problems in the classical theory of elasticity[END_REF][START_REF] Johnson | Contact Mechanics[END_REF][START_REF] Doyen | Time-integration schemes for the finite element dynamic Signorini problem[END_REF]. Enfin, A. Signorini a introduit une formulation mathématique pour le problème général de l'équilibre d'un corps élastique en contact sans frottement avec une structure rigide [START_REF] Signorini | Sopra alcune questioni di statica dei sistemi continui[END_REF]. G. Fichera a appelé les conditions de contact introduites par A. Signorini, les conditions de Signorini. Il s'est par ailleurs intéressé à l'existence de l'équilibre statique par minimisation de l'énergie totale sur un convexe dans [START_REF] Fichera | Problemi elastostatici con vincoli unilaterali: il problema di signorini con ambique condizioni al contorno[END_REF]. En particulier, il est le premier à avoir formulé le problème de contact en utilisant les inéquations variationnelles. Cette formulation mathématique du contact est à l'origine de nombreux travaux utilisant différentes théories comme les inéquations variationnelles, l'analyse de Fourier, la théorie du point fixe ou encore la méthode des opérateurs d'hystérésis.

Etat de l'art

Les travaux mathématiques s'articulent autour de la formulation de la condition de contact introduite par A. Signorini. Cette condition peut être relaxée et remplacée par un problème approché faisant intervenir une pénalisation appelée encore compliance normale (cf. [START_REF] Martins | Existence and uniqueness results for dynamic contact problems with nonlinear normal and friction interface laws[END_REF][START_REF] Klarbring | Frictional contact problems with normal compliance[END_REF][START_REF] Klarbring | A global existence result for the quasistatic frictional contact problem with normal compliance[END_REF][START_REF] Andersson | A quasistatic frictional problem with a normal compliance penalization term[END_REF][START_REF] Andersson | A quasistatic frictional problem with normal compliance[END_REF][START_REF] Eck | Existence results for the static contact problem with Coulomb friction[END_REF][START_REF] Eck | Existence results for the semicoercive static contact problem with coulomb friction[END_REF][START_REF] Bock | On hyperbolic contact problems[END_REF][START_REF] Sofonea | Analysis of a viscoelastic contact problem with multivalued normal compliance and unilateral constraint[END_REF][START_REF] Martins | On a quasi static contact problem in viscoelsaticity with normal compliance[END_REF]). Le problème ainsi obtenu est beaucoup plus simple à étudier, il est ainsi aisé de prouver l'existence et l'unicité d'une solution en utilisant différentes théories comme les inéquations variationnelles ou encore la méthode du point fixe. Ensuite, toute la difficulté réside dans le passage à la limite dans le problème approché qui s'obtient généralement en effectuant des estimations a priori. Il existe par ailleurs d'autres approches permettant d'obtenir l'existence ainsi que l'unicité pour des problèmes de contact comme par exemple la reformulation du problème en utilisant les inclusions différentielles ou encore la discrétisation. Nous ferons dans la suite une synthèse des résultats existants pour les problèmes de contact statiques et dynamiques ainsi que des méthodes numériques permettant d'obtenir une solution approchée pour ces problèmes.

Le premier résultat d'existence d'une solution pour un problème statique modélisant un corps élastique en contact avec un obstacle rigide sans frottement est obtenu dans [START_REF] Fichera | Problemi elastostatici con vincoli unilaterali: il problema di signorini con ambique condizioni al contorno[END_REF]. Quant au problème avec frottement de Tresca, il a été étudié quelques années plus tard dans [START_REF] Dautray | Analyse mathématique et calcul numérique pour les sciences et les techniques[END_REF]. Le problème de Signorini avec frottement de Coulomb a été traité dans [START_REF] Necas | On the solution of the variational inequality to the signorini problem with small friction[END_REF]. Plus précisément, sous l'hypothèse que le coefficient de frottement est suffisamment petit, l'existence d'une solution a été obtenue en utilisant l'algorithme itératif de Panagiotopoulos [START_REF] Panagiotopoulos | A nonlinear programming approach to the unilateral contact-, and friction-boundary value problem in the theory of elasticity[END_REF] ainsi que le théorème du point fixe. Par ailleurs, les méthodes de "shift" (cf. [START_REF] Agmon | Estimates near boundary for solutions of elliptic partial differential equations satisfying general boundary conditions[END_REF][START_REF] Fichera | Existence theorems in elasticity. boundary value problems of elasticity with unilateral constraints[END_REF]) ont permis de mettre en évidence un résultat de régularité. Certaines techniques employées dans [START_REF] Necas | On the solution of the variational inequality to the signorini problem with small friction[END_REF] ont été utilisées ensuite pour traiter des domaines coercifs [START_REF] Jarusek | Contact problems with bounded friction. Coercive case[END_REF] et semi-coercifs [START_REF] Jarusek | Contact problems with bounded friction. Semicoercive case[END_REF]. Ces résultats ont été redémontrés par des méthodes de pénalisation dans [START_REF] Eck | Existence results for the static contact problem with Coulomb friction[END_REF][START_REF] Eck | Existence results for the semicoercive static contact problem with coulomb friction[END_REF]. Remarquons que l'existence de solutions pour le problème de Signorini statique avec une régularisation du terme de frottement a été prouvée dans [START_REF] Cocu | Existence of solutions of Signorini problems with friction[END_REF]. En ce qui concerne la non-unicité de la solution de ce problème avec un coefficient de frottement suffisamment grand a été mis en évidence dans [START_REF] Hild | An example of nonuniqueness for the continuous static unilateral contact model with coulomb friction[END_REF][START_REF] Hild | Non-unique slipping in the Coulomb friction model in two-dimensional linear elasticity[END_REF]. D'autres contre-exemples d'unicité de solution ont été présentés récemment pour le cas statique et quasi-statique dans [START_REF] Ballard | Steady sliding frictional contact problems in linear elasticity[END_REF] et dans [START_REF] Ballard | A counter-example to uniqueness in quasi-static elastic contact problems with small friction[END_REF], respectivement. Par ailleurs, dans [START_REF] Martins | On an example of nonexistence of solution to a quasistatic frictional contact problem[END_REF] un exemple de non-existence pour une problème quasistatique avec frottement est étudié. Enfin, sous des hypothèses sur le déplacement tangentiel et pour un coefficient de frottement suffisamment petit, un résultat d'unicité a été démontré dans [START_REF] Renard | A uniqueness criterion for the Signorini problem with Coulomb friction[END_REF]. Le lecteur est également invité à consulter [START_REF] Eck | Unilateral contact problems : variational methods and existence theorems[END_REF] où un certain nombre de résultats sur les problèmes de contact statiques et quasi-statiques sont présentés. L'étude mathématique de problèmes de contact dynamiques est généralement plus délicate que celle des problèmes de contact statiques ou quasi-statiques. Les difficultés mathématiques sont liées à la nature hyperbolique des problèmes élastodynamiques combinée à la contrainte exprimée en déplacement conduisant à des irréguliarités intrinsèques du problème. Notons que dans certains travaux (cf. [START_REF] Dautray | Analyse mathématique et calcul numérique pour les sciences et les techniques[END_REF]), la contrainte en déplacement est remplacée par la contrainte en vitesse levant ainsi une difficulté majeure à l'étude de problèmes de contact dynamiques. Commençons par nous intéresser aux problèmes de contact dynamiques monodimensionnel. Sous l'hypothèse de choc élastique parfait, l'existence et l'unicité pour un problème de corde vibrante avec obstacle rigide et un obstacle ponctuel et avec une loi d'impact arbitraire sont établis respectivement dans [START_REF] Amerio | Study of the motion of a string vibrating against an obstacle[END_REF] et [START_REF] Amerio | On the motion of a string vibrating through a moving ring with a continuously variable diameter[END_REF]. Ensuite, des résultats d'existence et d'unicité ont été obtenus pour une corde vibrante avec un obstacle continu dans [START_REF] Bamberger | New results on the vibrating string with a continuous obstacle[END_REF], avec un obstacle concave dans [START_REF] Schatzman | A hyperbolic problem of second order with unilateral constraints: The vibrating string with a concave obstacle[END_REF], avec un obstacle ponctuel dans [START_REF] Schatzman | Un problème hyperbolique du 2ème ordre avec contrainte unilatérale: la corde vibrante avec obstacle ponctuel[END_REF] et avec obstacle adhésif dans [START_REF] Citrini | On the motion of a string vibrating against a gluing obstacle[END_REF]. Enfin le mouvement de deux cordes vibrantes se rejoignant en un point est traité dans [START_REF] Citrini | Some unilateral problems for the vibrating string equation[END_REF]. Quant aux problèmes de contact dynamiques en dimension supérieure à 1, il existe très peu de résultats significatifs dans la litérature. Ainsi l'existence et l'unicité d'une solution est établie pour l'équation des ondes avec des contraintes unilatérales dans un demi-espace en utilisant des techniques microlocales dans [START_REF] Lebeau | A wave problem in a half-space with a unilateral constraint at the boundary[END_REF], il n'a pas été cependant possible d'étendre ces techniques à domaines bornés. L'existence de la solution pour un domaine borné utilisant une technique de compacité par compensation est étudiée dans [START_REF] Kim | A boundary thin obstacle problem for a wave equation[END_REF]. L'existence de la solution est etablie pour un problème de contact multidimensionnelle, en dynamique et quasi-statique, avec prise en compte d'un terme thermo-élastique dans [START_REF] Muñoz Rivera | Multidimensional contact problems in thermoelasticity[END_REF]. La question de l'unicité et de la conservation de l'énergie restent des problèmes ouverts en dimension supérieure à 1. Certaines difficultés rencontrées pour des problèmes de contact elastodynamiques peuvent être levées en ajoutant de la viscosité dans les équations (cf. [START_REF] Jarusek | Contact problems with given time-dependent friction force in linear viscoelasticity[END_REF][START_REF] Jarusek | Variational inequality for a viscous drum vibrating in the presence of an obstacle[END_REF][START_REF] Jarusek | Solvability of unilateral hyperbolic problems involving viscoelasticity via penalization[END_REF][START_REF] Jarusek | Dynamic contact problems with given friction for viscoelastic bodies[END_REF][START_REF] Jarusek | Dynamic contact problems with small Coulomb friction for viscoelastic bodies. Existence of solutions[END_REF][START_REF] Eck | Existence of solutions for the dynamic frictional contact problem of isotropic viscoelastic bodies[END_REF]). Par ailleurs, les espaces fonctionnels caractérisant les traces ont été identifiés dans [START_REF] Petrov | Mathematical results on existence for viscoelastodynamic problems with unilateral constraints[END_REF] pour l'équation des ondes amorties ainsi que pour l'evolution d'un matériau viscoélastique de Kelvin-Voigt, les deux problèmes ont des conditions unilatérales aux bords. Cependant ces résultats de régularité ne permettent pas de déduire que les pertes d'énergie sont purement visqueuses. Il est toutefois possible sous certaines hypothèses supplémentaires de construire une solution permettant de conclure que les pertes d'énergie sont purement visqueuses (cf. [START_REF] Petrov | Viscoélastodynamique monodimensionnelle avec conditions de Signorini[END_REF][START_REF] Petrov | A pseudifferential linear complementarity problem related to one dimensional viscoelastic model with Signorini conditions[END_REF]).

Comme nous l'avons évoqué précédemment, l'unicité d'un problème de contact reste ouvert en dimension strictement supérieure à un et pour un domaine quelconque. C'est pour cette raison que les approches numériques pour résoudre les problèmes complexes sont utilisés. En effet, ces approches ne cessent de se développer grâce à l'evolution rapide des capacités des ordinateurs. En dépit de ce développement, la modélisation correcte des problèmes de contact reste encore un challenge à cause de la non-linéarité du contact et du frottement .

La méthode des élèments finis est largement utilisée pour simuler numériquement des problèmes de contact et consiste à discrétiser la géométrie du corps en un ensemble d'éléments de taille finie. La première application de la méthode des éléments finis pour les problèmes de contact a été faite dans [START_REF] Wilson | Finite element analysis of elastic contact problems using differential displacements[END_REF][START_REF] Chan | A finite element method for contact problems of solid bodies-Part i. theory and validation[END_REF][START_REF] Conry | A mathematical programming method for design of elastic bodies in contact[END_REF]. Dans ces travaux, les auteurs se sont intéressés aux problèmes d'élasticité sans frottement. Ils ont été suivis par plusieurs autres travaux, citons notamment Présentation de la thèse [START_REF] Hughes | A finite element method for a class of contact-impact problems[END_REF][START_REF] Francavilla | A note on numerical computation of elastic contact problems[END_REF]. A noter que la première analyse mathématique de convergence a été effectuée dans [START_REF] Haslinger | Contact between elastic bodies. II. Finite element analysis[END_REF]. La solution du problème de contact avec frottement par éléments finis est traitée dans [START_REF] Fredriksson | Finite element solution of surface nonlinearities in structural mechanics with special emphasis to contact and fracture mechanics problems[END_REF]. Dans [START_REF] Sachdeva | A finite element solution for the two-dimensional elastic contact problems with friction[END_REF], les auteurs ont développé une procédure utilisée dans [START_REF] Francavilla | A note on numerical computation of elastic contact problems[END_REF] permettant de tenir compte des effets de frottement. D'autres analyses par éléments finis ont été réalisées afin de décrire le problème de contact avec frottement en élastostatique dans [START_REF] Campos | A numerical analysis of a class of contact problems with friction in elastostatics[END_REF] et en élastodynamique dans [START_REF] Martins | A numerical analysis of a class of problems in elastodynamics with friction[END_REF].

Imposer la condition de contact dans le cadre de la méthode des élements finis (autrement dit transformer les inégalités variationnelles en égalités) est un autre sujet important pour les aspects numériques. Les deux méthodes les plus utilisées sont les multiplicateurs de Lagrange et la méthode de pénalisation. La premiére méthode a été introduite dans [START_REF] Hughes | A finite element method for a class of contact-impact problems[END_REF], et elle consiste à écrire la contrainte de contact en utilisant la force de contact. Une méthode pour inclure les effets de frottement en utilisant des facteurs de frottement de Coulomb en termes de multiplicateurs de Lagrange est introduite dans [START_REF] Okamoto | Finite element incremental contact analysis with various frictional conditions[END_REF]. Les premiers travaux a avoir utilisé la méthode de pénalisation pour les problèmes de contact sont [START_REF] Chan | A finite element method for contact problems of solid bodies-Part i. theory and validation[END_REF][START_REF] Babuska | The finite element method with penalty[END_REF][START_REF] Fredriksson | Finite element solution of surface nonlinearities in structural mechanics with special emphasis to contact and fracture mechanics problems[END_REF][START_REF] Kikuchi | Convergence of a penalty-finite element approximation for an obstacle problem[END_REF]. Dans cette méthode, l'interpénétration de deux surfaces de contact est possible et la force de contact normale est liée à la pénétration par un paramètre arbitraire de pénalisation. Pour plus d'informations sur les avantages et les inconvénients des méthodes, de multiplicateurs de Lagrange et de pénalisation, voir [START_REF] Kikuchi | Contact problems in elasticity: a study of variational inequalities and finite element methods[END_REF]. On peut citer aussi les autres techniques qui existent dans la littérature, comme la méthode de la flexibilité [START_REF] Francavilla | A note on numerical computation of elastic contact problems[END_REF], l'approche de programmation mathématique [START_REF] Hung | Frictionless contact of elastic bodies by finite element method and mathematical programming technique[END_REF][START_REF] Klarbring | A mathematical programming approach to contact problems with friction and varying contact surface[END_REF], la méthode du Lagrangien augmenté [START_REF] Simo | An augmented Lagrangian treatment of contact problems involving friction[END_REF]. La méthode des éléments finis classiques combine une discrétisation en espace et un schéma d'intégration en temps pour résoudre le problème semi-discret. Dans le cadre de problèmes linéaires, un schéma converge s'il vérifie la stabilité linéaire (par exemple, zérostabilité et A-stabilité [START_REF] Lambert | Computational methods in ordinary differential equations[END_REF]) et la condition de consistance, lorsque la taille du pas de temps tend vers zéro. Dans les problèmes élastodynamiques linéaires simples, on peut utiliser un schéma explicite conditionnellement stable. Pour les cas complexes (par exemple pour les matériaux incompressibles), les schémas implicites inconditionnellement stables ont été développés dans [START_REF] Newmark | A method of computational for structural dynamics[END_REF][START_REF] Hilber | Improved numerical dissipation for time integration algorithms in structural dynamics[END_REF]. Ces schémas perdent leurs stabilités inconditionnelles dans les cas des problèmes nonlinéaires. Pour ce type de problèmes, d'autres critères de stabilité doivent être utilisés, comme par exemple, l'analyse énergétique du problème. Il est donc très difficile de construire un schéma numérique stable pour les problèmes dynamiques avec contraintes unilatérales. Par conséquent, deux critères importants pour la stabilité d'une méthode sont utilisées pour justifier les simulations numériques pour les problèmes de contact. Ces critères sont le comportement de l'énergie et les oscillations parasites du déplacement et des forces de contact [START_REF] Hauret | Energy-controlling time integration methods for nonlinear elastodynamics and low-velocity impact[END_REF]. Après nombreux efforts, une méthode efficace et robuste a été proposée. Une possibilité consiste à ajouter des termes d'amortissement mais nous avons une perte de la précision de la solution. Une autre possibilité est d'impliciter la force de contact [START_REF] Carpenter | Lagrange constraints for transient finite element surface contact[END_REF] mais cela conduit à une perte d'énergie, même avec un raffinement du pas de temps. Certaines méthodes d'intégration temporelle sont proposées dans les articles [START_REF] Laursen | Design of energy conserving algorithms for frictionless dynamic contact problems[END_REF][START_REF] Chawla | Energy consistent algorithms for frictional contact problems[END_REF][START_REF] Laursen | Improved implicit integrators for transient impact problemsgeometric admissibility within the conserving framework[END_REF][START_REF] Barboteu | An energy-conserving algorithm for nonlinear elastodynamic contact problems -extension to a frictional dissipation phenomenon. Analysis and Simulation of Contact Problems[END_REF], elles conservent l'énergie mais n'évitent pas complètement des oscillations du déplacement et de la force de contact. Une autre façon d'éviter les oscillations parasites est d'imposer la compatibilité de la vitesse et de l'accélération des points de contact grâce aux vecteurs d'erreur correspondants (voir par exemple [START_REF] Lee | A numerical solution for dynamic contact problems satisfying the velocity and acceleration compatibilities on the contact surface[END_REF][START_REF] Armero | Formulation and analysis of conserving algorithms for frictionless dynamic contact/impact problems[END_REF]). Cette méthode n'est pas totalement satisfaisante. Il est bien connu que la semi-discrétisation à l'aide des éléments finis classiques conduit à un problème mal posé qui est provoqué par la multiplicité de solution du problème [START_REF] Khenous | Mass redistribution method for finite element contact problems in elastodynamics[END_REF]. L'ajout d'une loi d'impact supplémentaire, permet d'obtenir un problème semi-discret bien posé [START_REF] Paoli | Time discretization of vibro-impact[END_REF], mais cela conduit à des oscillations parasites. Une autre approche consistant à supprimer la masse au niveau du noeud de contact est proposée dans [START_REF] Khenous | Mass redistribution method for finite element contact problems in elastodynamics[END_REF] appelée la méthode de redistribution de masse, cette méthode évite les oscillations parasites au bord de contact. Le problème semi-discret, en utilisant la méthode de redistribution de masse, est un problème bien posé sans oscillations parasites et la solution conserve l'énergie (voir [START_REF] Khenous | Mass redistribution method for finite element contact problems in elastodynamics[END_REF] et chapitre 1). La convergence du problème semi-discret pour des matériaux visco-élastiques pour la dimension deux et trois est établie dans [START_REF] Doyen | Convergence of a space semi-discrete modified mass method for the dynamic Signorini problem[END_REF]. Nous avons établi la convergence de la méthode des matériaux élastiques en dimension un (voir chapitre 1).

Formulations mathématiques

La condition de contact introduite par A. Signorini s'écrit [START_REF] Brinkmeier | A finite element approach for the simulation of tire rolling noise[END_REF] u ν ≤ 0, r ν ≤ 0, u ν r ν = 0, où u ν et r ν sont respectivement le déplacement normal et la force de contact. La condition [START_REF] Brinkmeier | A finite element approach for the simulation of tire rolling noise[END_REF] indique que si le contact n'est pas établi, u ν < 0 alors r ν = 0, autrement dit la réaction entre le matériau et l'obstacle rigide n'existe pas. Lorsque le contact est établi, (1) permet de déduire que u ν = 0 et r ν < 0 en d'autres termes une réaction existe entre le matériau et l'obstacle rigide. Nous utiliserons dans suite les notations u t

ε kl (u) = 1 2
( ∂u k ∂x l + ∂u l ∂x k
) .

Nous supposerons que le matériau est homogène et isotrope. En utilisant la loi de Hooke, le tenseur d'élasticité s'exprime alors en fonction des constantes de Lamé λ et µ comme suit

A ijkl def = λδ ij δ kl + 2µδ kl δ jk .
Par conséquent, le tenseur des contraintes peut se réécrire de la manière suivante:

σ(u) = λTr(ε)I + 2µε,
u(•, 0) = u 0 et u t (•, 0) = v 0 dans Ω d , et des conditions aux limites u = α sur Γ Dir d × (0, T ) et σ(u)ν = g sur Γ Neu d × (0, T ), (4a) 0 ≥ ν T u ⊥ ν T (σ(u)ν) ≤ 0 et τ i T (σ(u)ν) = 0, i = 1, 2 sur Γ Sig d × (0, T ), (4b)
où ν et τ désignent la normale unitaire extérieure et la tangente unitaire, α est un vecteur donné et g est une force de surfacique dépendant de l'espace, (•) T est la transposée de (•) L'orthogonalité a ici le sens suivant; si nous avons suffisamment de régularité, cela signifie que le produit ν T u(ν T σ(u)ν) est nul presque partout sur la frontière. Si ce n'est pas le cas, l'inégalité (4b) est intégrée sur un ensemble approprié de fonctions de test, ce qui conduit à une formulation faible pour la condition unilatérale. En particulier, Γ Dir d = ∅ dans le cas où nous considérons le problème d'une barre en chute libre heurtant un obstacle rigide.

Nous allons maintenant définir la formulation faible associée au problème (2)-( 4). Nous introduisons au préalable le convexe

K def = {u ∈ U : u(•, t) ∈ V p.p. t et ν T u ≤ 0 p.p. sur Γ Sig d }, où V def = {u ∈ H 1 (Ω d ) : u = 0 p.p. sur Γ Dir d } et U def = {u ∈ L 2 (0, T ; V) : u t ∈ L 2 (0, T ; L 2 (Ω d ))} avec X désignant l'espace X d .
Nous obtenons la formulation faible associée à (2)-(4) en multipliant (2) par v-u, v ∈ K, puis en intégrant formellement cette expression sur Ω d , on obtient:

         trouver u : [0, T ] → K telle que ∫ Ω d ρu tt •(v-u) dx + ∫ Ω d A ijkl ε ij (u)ε kl (v -u) dx ≥ ∫ Ω d •(v -u) dx + ∫ Γ Neu d g•(v -u) dx, pour tout v ∈ K.
Dans les trois premiers chapitres de cette thèse, nous allons nous restreindre au cas monodimensionel qui après adimensionnement s'écrit comme suit:

(5)

u tt -u xx = dans Ω 1 × (0, T ), avec les conditions initiales (6) u(•, 0) = u 0 et u t (•, 0) = v 0 sur Ω 1 , et des conditions aux limites (7) u(L, •) = α, u x (L, •) = g et 0 ≤ u(0, •) ⊥ u x (0, •) ≤ 0 sur [0, T ].
Comme précédemment, nous allons définir la formulation faible associée au problème ( 5)- [START_REF] Gladwell | Contact problems in the classical theory of elasticity[END_REF].

Soit le convexe

K def = {u ∈ U : u(•, t) ∈ V p.p. t et u(0) ≥ 0}, où V def = {u ∈ H 1 (0, L) : u(L) = 0} et U def = {u ∈ L 2 (0, T ; V) : u t ∈ L 2 (0, T ; L 2 (0, L))}.
Par conséquent, la formulation faible associée à (2)-( 4) est obtenue en multipliant (2) par v -u, v ∈ K, puis en intégrant cette expression formellement sur (0, L) × (0, T ), ce qui nous donne (8)

                 trouver u ∈ K tel que - ∫ L 0 v 0 (v(•, 0)-u 0 ) dx - ∫ T 0 ∫ L 0 u t (v t -u t ) dx dt + ∫ T 0 ∫ L 0 u x (v x -u x ) dx dt ≥ ∫ T 0 (v-u), pour tous v ∈ K pour lesquels il existe ζ > 0 avec v = u pour t ≥ T -ζ.

Résumé des travaux

Cette thèse s'insére dans la problèmatique d'une bonne approximation numérique des problèmes de contact en élastodynamique. En particulier, nous nous attachons à caractériser les propriétés de la méthode de redistribution de masse (stabilité et convergence) qui consiste à redistribuer la masse des noeuds de contact sur les autres noeuds. La thèse est organisée en quatre chapitres. Les trois premiers chapitres sont consacrés à l'étude du mouvement d'une barre pouvant heurter un obstacle rigide. La méthode de redistribution de masse est introduite et son efficacité est validée numériquement. Le dernier chapitre de cette thèse est dédié au problème en dimension 2 et 3. Toutes les simulations numériques présentées ont été effectuées en utilisant la bibliothèque d'éléments finis Getfem++ (cf. [START_REF] Renard | Getfem++. An Open Source generic C++ library for finite element methods[END_REF]).

Chapitre 1: Convergence de la méthode de redistribution de masse pour l'équation des ondes monodimensionnelle avec une contrainte unilatérale à la frontière.

Nous nous intéressons dans ce chapitre au mouvement d'une barre élastique de longueur L dont l'une des extrémités est encastrée et l'autre est libre de bouger tant qu'elle ne touche pas un obstacle rigide (cf. Figure 5). On considère le déplacement vertical u(x, t) de cette barre

L x = 0
def = x + t et η def = x -t,
il est possible de décomposer la solution u en la somme de deux fonctions p et q de la manière suivante:

u(ξ, η) = p(ξ) + q(η), ce qui conduit à l'expression u(x, t) = p(x+t) + q(x-t).
Les conditions initiales ainsi que la condition de contact peuvent être réécrites comme suit

p(x) + q(x) = u 0 (x) et p (x) -q (x) = v 0 (x), 0 ≤ p(t) + q(-t) ⊥ p (t) + q (-t) ≤ 0 et p(L+t) + q(L-t) = 0.
La sous-différentielle J N de la fonction indicatrice qui est une application multivoque définie par

J N [x] def =        {0} si x < 0, [0, +∞) si x = 0, ∅ si x > 0,
permet de reformuler le problème ( 5)-( 7) et d'écrire

f (t) ∈ -J N [f ](t) -2p (t) p.p. t ∈ (0, L), (9a) f (0) = -u 0 (0). (9b)
Les propriétés de monotonie et de convexité de l'opérateur conduisent au théorème d'existence et d'unicité pour le problème (9) explicité ci-dessous. Théorème 1. Supposons que p(t) ∈ W 1,1 (0, L). Alors le problème de Cauchy (9) admet une solution unique f (t) ∈ W 1,1 (0, L).

Nous allons maintenant établir que (8) et (9) sont équivalents. Pour cela nous divisons le domaine [0, L] × [0, T ] en quatre domaines I-IV comme représenté sur la Figure 6. Dans chacun des domaines, une expression de la solution u en fonction de p et q est obtenue conduisant au théorème suivant: Théorème 2. Supposons que u 0 and v 0 appartiennnent respectivement à K et L 2 (0, L). Alors f est solution du problème [START_REF] Doyen | Time-integration schemes for the finite element dynamic Signorini problem[END_REF] si et seulement si u est une solution de (8) pour T = L Par ailleurs, la solution u conserve l'énergie, c'est-à-dire nous avons

I II III IV x t (L, L) (L, 0) (0, L)
E(t) = 1 4 ∫ L 0 |u 0 (x)+v 0 (x)| 2 dx + 1 4 ∫ L 0 |u 0 (x)-v 0 (x)| 2 dx est constante.
Nous allons introduire maintenant des problèmes discretisés en temps et en espace approchant la formulation faible [START_REF] Johnson | Contact Mechanics[END_REF]. Dans un premier temps, un problème semi-discretisé en espace est introduit. Ce problème utilise la méthode des éléments finis affines de Lagrange associé à une méthode des multiplicateurs de Lagrange. Plus précisément, la solution u est approchée par

u h (x, t) = n-1 ∑ j=0 u j (t)ϕ j (x), où la fonction de base ϕ i (x) est definie comme suit ϕ i (x) def = 1-|x-a i | h lorsque x ∈ [a max(i-1,0) , a i+1 ] et ϕ i (x)
               trouver u h : [0, T ] → V h et λ : [0, T ] → R tels que v h ∈ V h ∫ L 0 u h,tt v h dx + ∫ L 0 u h,x v h,x = -λv h (0) p.p. t ∈ [0, T ], 0 ≤ u h (0, •) ⊥ λ ≤ 0 p.p. t ∈ [0, T ], u h (•, 0) = u 0 h et u h,t (•, 0) = v 0 h .
Notons que [START_REF] Signorini | Sopra alcune questioni di statica dei sistemi continui[END_REF] peut être réécrit de la manière suivante:

(11)            trouver U : [0, T ] → R n et λ : [0, T ] → R tels que M U tt + SU = -λe 0 p.p. t ∈ [0, T ], 0 ≤ u 0 ⊥ λ ≤ 0 p.p. t ∈ [0, T ], U (0) = U 0 et U t (0) = V 0 ,
où M et S sont appelées respectivement les matrices de masse et de rigidité et e 0 def = (1, 0, . . . , 0) T . Nous pouvons remarquer que [START_REF] Fichera | Problemi elastostatici con vincoli unilaterali: il problema di signorini con ambique condizioni al contorno[END_REF] n'est a priori pas bien posé car l'équation au noeud de contact peut avoir plusieurs solutions. La méthode de redistribution de masse est alors introduite. Cette méthode est basée sur une redistribution de masse de manière à ce qu'il n'y ait pas d'inertie au niveau du noeud de contact, supprimant ainsi l'opérateur différentiel en temps au noeud de contact. En d'autres termes, la matrice de masse M est remplacée par la matrice suivante

M mod def = ( 0 0 0 M ) , où Mij = M i+1,j+1 pour tous i, j = 1, . . . , n -2.
Il est alors commode d'introduire les notations suivantes:

S = ( S 00 C T C S ) et U = ( u 0 Ū ) , Présentation de la thèse où Ū def = (u 1 , . . . , u n-1 ) T et Sij def = S i+1,j+1 avec C def = ∫ Ω ϕ i+1 ϕ 0 dx, i = 0, . . . , n -2.
Avec ces notations et en se servant de la méthode de redistribution de masse, le problème se réécrit comme suit

(12)                  trouver U : [0, T ] → R n et λ : [0, T ] → R tels que M Ūtt + S Ū = -Cu 0 p.p. t ∈ [0, T ], S 00 u 0 + C T Ū = -λ p.p. t ∈ [0, T ], 0 ≤ u 0 ⊥ λ ≤ 0 p.p. t ∈ [0, T ], U (0) = U 0 et U t (0) = V 0 .
Puisque la matrice M n'est pas singulière, on déduit que [START_REF] Martins | Existence and uniqueness results for dynamic contact problems with nonlinear normal and friction interface laws[END_REF] possède une unique solution (U, λ) qui est lipschitzienne. Par ailleurs, cette solution conserve l'énergie. En choisissant de manière judicieuse un projecteur qui conserve la contrainte au point de contact puis en effectuant des estimations a priori, le théorème suivant en découle:

Théorème 3. Supposons que lim h→0 ( u h (•, 0)-u 0 H 1 (0,L) + u h,t (•, 0)-v 0 L 2 (0,L) ) = 0.
Alors la solution du problème approché [START_REF] Martins | Existence and uniqueness results for dynamic contact problems with nonlinear normal and friction interface laws[END_REF] converge vers la solution unique de la formulation faible [START_REF] Johnson | Contact Mechanics[END_REF].

Enfin la dernière partie du premier chapitre est consacrée à quelques tests numériques qui confirment l'efficacité de la méthode de redistribution de masse. En effet, nous observons sur la Figure 7 (gauche) que plus les pas de temps ∆t et d'espace ∆x sont raffinés plus la solution approchée est proche de la solution exacte. Nous pouvons faire la même observation pour l'évolution d'énergie (cf. Figure 7 (droite)).

Figure 7: Convergence numérique de la solution associée au problème avec la méthode de redistribution de masse vers la solution exacte au point de contact x = 0 (à gauche). Convergence numérique de l'évolution de l'énergie (à droite).

12 Chapitre 2: Approximations numériques d'un problème de contact élastodynamique unidimensionnel basé sur la méthode de redistribution de masse.

Nous comparons ici différents schémas numériques que nous avons développés pour le problème introduit dans le chapitre précédent avec cette fois-ci un second membre non nul (x, t), dépendant de x et de t. Nous analysons également le cas de la chute libre d'une barre qui heurte un obstacle rigide. Pour ces deux problèmes, une solution explicite est mise en évidence puis une comparaison est effectuée avec les solutions approcheés obtenues pour différents schémas numériques utilisant notamment la méthode de redistribution de masse.

Nous commençons par établir que le problème ( 5)-( 7) possède une solution unique. Pour cela, nous introduisons ū qui est la solution unique du problème auxiliaire ( 5)-( 7) où la condition contact a été remplacée par une condition Dirichlet homogène. Ensuite, nous définissons v def = u-ū comme la solution du problème

v tt (x, t) -v xx (x, t) = 0, (x, t) ∈ (0, L) × (0, T ), v(x, 0) = 0 et v t (x, 0) = 0, x ∈ (0, L), 0 ≤ v(0, t) ⊥ v x (0, t) + ūx (0, t) ≤ 0 et v(L, t) = 0, t ∈ [0, T ].
En utilisant la méthode des caractéristiques ainsi que la notation w(t) def = -p(t) -q(-t), nous sommes en mesure de montrer que le problème précédent est équivalent à un problème d'inclusion différentielle. Ce problème est ensuite régularisé puis des estimations a priori sont obtenues permettant de conclure que notre problème admet une unique solution u ∈ L ∞ (0, T ; V) ∩ W 1,∞ (0, T ; H) sous des hypothèses de régularité sur les données initiales ainsi que sur la force externe .

La stratégie adoptée pour obtenir le problème semi-discrétisé en espace est la même que celle présentée dans le premier chapitre; la seule différence vient du second membre qui est ici différent de zéro. Plus précisément, la méthode des éléments finis affines de Lagrange classique combinée avec la matrice de masse modifiée ou la matrice de masse non modifiée (appelée matrice de masse standard dans la suite) est étudiée. Des méthodes d'intégration en temps permettant d'éviter des oscillations artificielles pour le déplacement et pour la force de contact conservant l'énergie sont proposées. Nous avons ainsi utilisé trois méthodes d'intégration temporelle.

Les méthodes de Newmark:

                 trouver U n+1 h : [0, T ] → R n et λ n+1 : [0, T ] → R tel que: U n+1 h = U n h + ∆tU n h,t + ( 1 2 -β ) ∆t 2 U n h,tt + β∆t 2 U n+1 h,tt , U n+1 h,t = U n h,t + (1-γ)∆tU n h,tt + γ∆tU n+1 h,tt , M U n+1 h,tt + SU n+1 h = -λ n+1 e 0 + F n+1 , 0 ≤ u n+1 0 ⊥ λ n+1 ≤ 0. Présentation de la thèse La θ-méthode:                  trouver U n+1 h : [0, T ] → R n et λ n+1 : [0, T ] → R tel que: U n+1 h = U n h + ∆t((1-θ)U n h,t +θU n+1 h,t ), U n+1 h,t = U n h,t + ∆t((1-θ)U n h,tt +θU n+1 h,tt ), M U n+1 h,tt + SU n+1 h = -λ n+1 e 0 + F n+1 , 0 ≤ u n+1 0 ⊥ λ n+1 ≤ 0.
La méthode de Paoli-Schatzman:

                 trouver U n+1 h : [0, T ] → R n et λ n : [0, T ] → R tel que: M (U n+1 h -2U n h +U n-1 h ) ∆t 2 + S(βU n+1 h +(1-2β)U n h +βU n-1 h ) = -λ n e 0 pour tous n ≥ 2, 0 ≤ u n,e 0 = u n+1 0 + eu n-1 0 1 + e ⊥ λ n ≤ 0, U 0 et U 1 donné.
Chaque méthode proposée ci-dessus est testée avec une matrice de masse standard ainsi qu'une matrice de masse modifiée puis elle est comparée aux solutions analytiques. Notons que la solution analytique obtenue pour la chute libre d'une barre n'est pas suffisamment régulière, ce qui ne permet pas de mettre en évidence l'efficacité de la méthode de redistribution de masse. Concernant le second problème traité dans ce chapitre, une solution explicite périodique a été mise en évidence. Plus précisément, nous avons déterminé la solution en utilisant la méthode des caractéristiques sur trois phases du phénomène; avant, durant et après le contact. Dans les simulations numériques présentées ci-dessous, nous avons calculé la solution approchée en utilisant la méthode de Crank-Nicolson. Tout d'abord, nous avons comparé trois solutions approchées avec la solution analytique, en utilisant la matrice de masse standard et en choisissant des pas d'espace ∆x et de temps ∆t de telle sorte qu'ils tendent vers zéro (pour les solutions approchées). Puis ensuite, nous avons effectué un travail analogue en remplaçant la matrice de masse standard par la matrice de masse modifée. Enfin nous avons comparé les résultats obtenus pour les deux matrices de masse pour le déplacement de tous les points de la barre.

Le même travail a été effectué pour la force de contact ainsi que pour l'énergie. Nous avons également vérifié la convergence vers la solution analytique (cf. Figure 8). Nous présentons cidessus des courbes de convergence pour le déplacement dans les normes L p (0, T ; L 2 (0, L)) et L p (0, T ; H 1 (0, L)) avec p = 2, ∞, pour la force de contact et pour le multiplicateur de Lagrange dans la norme L 2 (0, T ) et pour l'énergie dans les normes L p (0, T ) avec p = 2, ∞, pour les méthodes de masse standard et modifée (cf. Figure 9). Nous avons mené les mêmes investigations pour des méthodes de Newmark en prenant comme paramètres

(β, γ) = ( 1 2 , 1 2 ) et (β, γ) = ( 1 2 , 1
), pour la méthode d'Euler implicite et pour la méthode Paoli-Schatzman avec des coefficients de restitution e = 0, 1 2 , 1. Les tests numériques ont montré que la méthode de redistribution de masse élimine les oscillations du déplacement, réduit les oscillations des forces de contact et conserve l'énergie à de très petites oscillations près. De plus la méthode de redistribution de masse améliore les taux de convergence. Chapitre 3 : Sur l'efficacité d'une méthode de redistribution de masse pondérée pour un problème élastodynamique avec des contraintes unilatérales.

Nous considérons, comme dans le chapitre 1, le mouvement d'une barre élastique dont l'une de ses extrémités est encastrée et l'autre est libre de tout mouvement tant qu'elle ne touche pas un obstacle rigide. La formulation mathématique de ce problème est donnée par ( 5)- [START_REF] Gladwell | Contact problems in the classical theory of elasticity[END_REF]. Ensuite le problème semi-discretisé en espace est introduit en utilisant encore une fois la méthode des éléments finis affines de Lagrange associée à un multiplicateur de Lagrange. La différence fondamentale avec les travaux précédents réside dans le fait que la masse enlevée au noeud de contact est maintenant redistribuée sur les autres noeuds et non plus juste supprimée (cf. chapitres 1 et 2). Par conséquent, la masse totale du matériau est préservée. Par ailleurs, une fonction pondérée permettant de redistribuer la masse du noeud de contact sur les autres noeuds tout en éliminant la masse du noeud de contact est introduite. La matrice de masse modifiée M mod est maintenant définie par

M mod ij def = ∫ L 0 ϕ i ϕ j w h,ij dx tel que w h,0i = w h,i0 pour tout i = 0, . . . , n -1 et w h,ij = 1 + o(h)
for all i, j = 1, . . . , n -1 où w h est la fonction pondérée mentionnée plus haut. En utilisant les notations du chapitre 1 et en se servant de la méthode de redistribution de masse, la formulation faible (8) est approchée par (13)

(a) u h -u L ∞ (0,T ;H) (b) u h -u L 2 (0,T ;H)
                 trouver U : [0, T ] → R n et λ : [0, T ] → R such that M Ūtt + S Ū = -Cu 0 + F p.p. t ∈ [0, T ], S 00 u 0 + C T Ū = -λ + F 0 a.e. t ∈ [0, T ], 0 ≤ u 0 ⊥ λ ≤ 0 a.e. t ∈ [0, T ], U (0) = U 0 et U t (0) = V 0 , où F def = ( ∫ L 0 f ϕ 1 dx, . . . , ∫ L 0 f ϕ n-1 dx)
                     trouver u h : [0, T ] → V h telle que pour tout v h ∈ K ∩ V h ∫ T 0 ∫ L 0 u h,tt (x, t)(v h (x, t)-u h (x, t))w h (x) dx dt + ∫ L 0 ∫ T 0 u h,x (x, t)(v h,x (x, t)-u h,x (x, t)) dx dt ≥ ∫ L 0 ∫ T 0 f (x, t)(v h (x, t)-u h (x, t)) dx dt.
Le résultat de convergence est résumé dans le théorème suivant:

Théorème 5. Supposons que lim h→0 ( u h (•, 0)-u 0 H 1 (0,L) + u h,t (•, 0)-v 0 L 2 (0,L) ) = 0.
Alors la solution du problème approché [START_REF] Klarbring | A global existence result for the quasistatic frictional contact problem with normal compliance[END_REF] converge vers l'unique solution de la formulation faible (8) lorsque h tend vers 0.

Dans les simulations numériques, nous avons considéré trois répartitions de masses correspondants respectivement à la non répartition de la masse du noeud de contact (Mod 1), à la répartition de la masse du noeud de contact uniformément sur les autres noeuds (Mod 2) ainsi qu'à la répartition de la masse du noeud de contact sur le noeud précédant le noeud de contact (Mod 3) (cf. Figure 10). La solution approchée est d'autant meilleure que la répartition de la masse du noeud de contact est faite sur un autre noeud proche du noeud de contact. En particulier, la meilleure approximation est obtenue dans le cas où toute la répartition de masse est faite sur le noeud qui précéde le contact (cf. Figure 10

Présentation de la thèse

Chapitre 4 : Etude numérique de la convergence de la méthode de redistribution de masse pour les problèmes de contact élastodynamique L'objectif du chapitre 4 est l'étude numérique de la convergence de la méthode de redistribution de masse pour des problèmes de contact élastodynamique en dimension deux et trois. Comme dans les chapitres précédents, nous considérons le contact d'une barre élastique et d'un obstacle rigide dont l'une des extrémités de la barre est encastrée et l'autre peut se déplacer librement tant qu'elle ne touche pas un obstacle rigide. La formulation mathématique de ce problème est donnée par (2)-(4). Ensuite, nous avons utilisé la méthode des éléments finis pour obtenir le problème semi-discrétisé ainsi que la matrice de masse modifiée. Par conséquent, le problème approché de (2)-(4) s'écrit

           trouver U N h : [0, T ] → N , U N ⊥ h : [0, T ] → N ⊥ et λ : [0, T ] → R nc telle que M mod U N ⊥ h,tt + S(U N h + U N ⊥ h ) = F + B T λ pour p.p. t ∈ [0, T ], 0 ≥ ν T i U N h ⊥ λ i ≤ 0 pour tous i ∈ I c et pour p.p. t ∈ [0, T ], U h (0) = U 0 h et U h,t (0) = V 0 h .
Les méthodes d'intégration en temps de Newmark ainsi que de Euler implicite ont été utilisées pour résoudre le problème semi-discret. Comme dans le second chapitre, chaque schéma est testé avec les matrices de masse standard et modifiée. En dimension deux, nous pouvons comparer le déplacement normal au niveau des noeuds de contact, la force de contact et l'évolution de l'énergie pour les deux matrices de masse (cf. Figures 11 et 12). En raison de la croissance rapide du nombre de degrés de liberté en dimension trois, nous avons pu obtenir des courbes de convergence avec seulement trois points. Cependant les courbes de convergence confirment l'efficacité de la matrice de masse modifiée.

Contrairement au problème monodimensionnel, nous n'avons pas pu déterminer une solution analytique associée à (2)-(4). Toutefois une solution de référence peut être obtenue en prenant ∆t et ∆x très petits. Nous présentons des courbes de convergence, en dimension deux, pour le déplacement, la force de contact ainsi que l'énergie.

Présentation de la thèse

(a) U n h -U L ∞ (0,T ;H) (b) U n h -U L 2 (0,T ;H)
a) U n h -U L ∞ (0,T ;H) (b) U n h -U L 2 (0,T ;H)

Conclusion et perspectives

Les principales contributions de ce travail de thèse sont l'étude théorique et numérique de la méthode de redistribution de masse en élastodynamique avec impact. Une attention particulière a été portée à l'étude numérique des taux de convergence. Plus précisément, nous avons mis en évidence que la méthode de redistribution de masse était particulièrement bien adaptée pour approcher des problèmes elastodynamiques avec des contraintes unilatérales. Les perspectives sont nombreuses à court et à long terme. Il serait intéressant de voir si les méthodes numériques développées dans cette thèse donnent des résultats similaires pour d'autres matériaux comme par exemple les matériaux à mémoire de forme ou encore les matériaux ferromagnétiques. Par ailleurs, les estimations d'erreur en temps et en espace représentent un challenge important. En effet, dans le cas monodimensionnel le problème continu est bien posé, une estimation d'erreur de la discrétisation devrait être possible à obtenir. La principale difficulté consiste à déterminer un lien entre la formulation éléments finis et la formulation en terme d'inclusion différentielle. Enfin, il semble que l'introduction d'une méthode à poids devrait permettre une écriture unique du problème approché conduisant ensuite à une estimation d'erreur en espace.

Introduction

We consider an elastic bar of length L vibrating vertically. More precisely, one end of the bar is free to move, as long as it does not hit a material obstacle, while the other end is clamped (see Fig. 1.1). The obstacle constrains the displacement of the extremity to be greater than or equal to 0.

We describe now the mathematical situation. We assume that the material of the bar is homogeneous and satisfies the theory of small deformations. Let x be the spatial coordinate along the bar, with the origin at the material obstacle, let u(x, t) be the displacement at time t ∈ [0, T ], T > 0 , of the material point of spatial coordinate x ∈ [0, L]. Then the mathematical problem can be formulated as follows:

(1.1) ü(x, t) -u (x, t) = 0, (x, t) ∈ (0, L) × (0, T ), Chapter 1. Convergence of mass redistribution method for the one-dimensional wave equation with a unilateral constraint at the boundary with Cauchy initial data (1.2) u(x, 0) = u 0 (x) and u(x, 0) = v 0 (x), x ∈ (0, L), and Signorini and Dirichlet boundary conditions at x = 0 and x = L, respectively,

(1.3) 0 ≤ u(0, t) ⊥ u (0, t) ≤ 0 and u(L, t) = 0, t ∈ [0, T ],
where u def = ∂u ∂t and u def = ∂u ∂x . The orthogonality has the natural meaning; namely if we have enough regularity, it means that at the boundary x = 0 the product u(0, •)u (0, •) vanishes for almost every t. If it is not the case, the above inequality is integrated on an appropriate set of test functions, leading to a weak formulation for the unilateral condition. Observe that from mathematical viewpoint, the Signorini conditions mean that when the bar touches the obstacle in x = 0, its reaction can be only upwards, so that u (0, t) ≤ 0 on the set {t : u(0, t) = 0}. While in the case where the bar does not touch the obstacle, its end is free to move. More precisely, we have u (0, t) = 0 on the set {t : u(0, t) > 0}.

We suppose that the initial displacement u 0 belongs to the Sobolev space H 1 (0, L) and satisfies the compatibility conditions, i.e. u 0 (L) = 0 and u 0 (0) ≥ 0 and the initial velocity v 0 belongs to L 2 (0, L).

We describe now the weak formulation of the problem. To this aim, we denote by K the following convex set:

K def = {u ∈ H 2 : u(•, t) ∈ K for almost every t}, where H 2 def = {u ∈ L 2 (0, T ; V ) : u ∈ L 2 (0, T ; L 2 (0, L))} and K def = {u ∈ V : u(0) ≥ 0} with V def = {u ∈ H 1 (0, L) : u(L) = 0}.
Then the weak formulation associated to (1.1)-(1.3) is obtained by multiplying (1.1) by v -u, v ∈ K and by integrating formally this result over Q T def = (0, L) × (0, T ):

(1.4)          find u ∈ K such that - ∫ L 0 v 0 (v(•, 0)-u 0 ) dx - ∫ Q T u( v-u) dx dt + ∫ Q T u (v -u ) dx dt ≥ 0 for all v ∈ K for which there exists ζ > 0 with v = u for t ≥ T -ζ.

Existence and uniqueness results by using the characteristic method

The weak formulation (1.4) is derived from [START_REF] Schatzman | Numerical approximation of a wave equation with unilateral constraints[END_REF] where the contact conditions are given in a slightly different context. Existence and uniqueness results are obtained for a vibrating string with concave obstacle in one-dimensional space in [START_REF] Schatzman | A hyperbolic problem of second order with unilateral constraints: The vibrating string with a concave obstacle[END_REF] and for a wave equation with unilateral constraint at the boundary in a half space of R N in [START_REF] Lebeau | A wave problem in a half-space with a unilateral constraint at the boundary[END_REF]. An existence result for a wave equation in a C 2 regular bounded domain constrained by an obstacle at the boundary in R N for N ≥ 2 is proven in [START_REF] Kim | A boundary thin obstacle problem for a wave equation[END_REF].

The paper is organized as follows. In Section 1.2, the problem (1.1)-(1.3) is reformulated as a differential inclusion problem by using characteristic method, which is a crucial ingredient to prove the uniqueness result. Then, the rest of this section is devoted to the proof of existence and uniqueness results as well as to the energy balance. In Section 1.3, the equivalence between the weak formulation associated to (1.1)-(1.3) and the differential inclusion obtained in Section 1.2 is established. Then, in Section 1.4, a mass redistribution method is introduced and its convergence is proved. This method is based on a redistribution of the body mass such that there is no inertia at the contact node (see [START_REF] Khenous | Mass redistribution method for finite element contact problems in elastodynamics[END_REF][START_REF] Hager | Analysis of a space-time discretization for dynamic elasticity problems based on mass-free surface elements[END_REF]). Finally, some numerical examples are reported and analyzed in Section 1.5. More precisely, the energy with and without the mass redistribution method are compared as well as the approximated solution associated to the mass redistribution method and an exact solution.

Existence and uniqueness results by using the characteristic method

This section is devoted to the proof of existence and uniqueness results for the problem (1.1)- (1.3). The first step consists in rewriting (1.1)-(1.3) as a differential inclusion problem by using the characteristic method. To this aim, we introduce the following notations:

ξ def = x + t and η def = x -t.
Therefore the chain rule gives

∂ 2 u ∂x 2 = ∂ 2 u ∂ξ 2 + 2 ∂ 2 u ∂ξ∂η + ∂ 2 u ∂η 2 and ∂ 2 u ∂t 2 = ∂ 2 u ∂ξ 2 -2 ∂ 2 u ∂ξ∂η + ∂ 2 u ∂η 2 ,
which by using (1.1) implies that ∂ 2 u ∂ξ∂η vanishes. Thus we may conclude that

u(ξ, η) = p(ξ) + q(η),
where p and q are two differentiable functions such that

(1.5) u(x, t) = p(x+t) + q(x-t).
In particular, taking t = 0 and using the initial data u 0 and v 0 , we get (1.6) p(x) + q(x) = u 0 (x) and p (x) -q (x) = v 0 (x), which by integration gives (1.7)

p(ξ) = u 0 (ξ) 2 + 1 2
∫ ξ 0 v 0 (x) dx and q(η) = u 0 (η) 2 - 1 2
∫ η 0 v 0 (x) dx
Chapter 1. Convergence of mass redistribution method for the one-dimensional wave equation with a unilateral constraint at the boundary for all ξ and η belonging to [0, L]. According to (1.5), the boundary conditions (1.3) can be rewritten as follows:

(1.8) 0 ≤ p(t) + q(-t) ⊥ p (t) + q (-t) ≤ 0 and p(L+t) + q(L-t) = 0 for all t belonging to [0, T ]. Thanks to the above identity, we may extend p(t) for all t ∈ [L, 2L],

i.e. we have p(L+t) = -q(L-t) for all t belonging to [0, L]. If we choose t = L + t, we get p(t ) = -q(2L-t ). We already have the solution for q(t) with 0 ≤ t ≤ L and if L ≤ t ≤ 2L, we can obtain p(t ) by observing that 0 ≤ 2L -t ≤ L and by using q(2L-t), it follows that (1.9)

p(t) = - u 0 (2L-t) 2 + 1 2 ∫ 2L-t 0 v 0 (x) dx for all t belonging to [L, 2L].
Let us introduce the multivalued function J N : R → P( R) defined by

J N (x) def =        {0} if x < 0, [0, +∞) if x = 0, ∅ if x > 0,
where P( R) is the set of all subsets of R. More precisely, J N (x) is the subdifferential of the indicator function ∂I (-∞,0] (x) defined by

I (-∞,0] (x) def = { 0 if x ∈ (-∞, 0], +∞ if x / ∈ (-∞, 0].
Obviously, I (-∞,0] is a lower semi-continuous and convex function, for further details the reader is referred to [START_REF] Brézis | Opérateurs maximaux monotones et semi-groupes de contractions dans les espaces de Hilbert[END_REF]. Then, the inequalities in (1.8) can be rewritten as follows (1.10) p (t) + q (-t) ∈ -J N (-p(t)-q(-t))).

Note that at this stage, q(η), η ∈ [-L, 0], is the unique unknown of (1.10). We define now

(1.11) f (t) def = -p(t) -q(-t).
We insert (1.11) into (1.10) to get

f (t) ∈ -J N (f (t)) -2p (t).
Finally, we find the following Cauchy problem

f (t) ∈ -J N (f (t)) -2p (t) a.e. t ∈ (0, L), (1.12a) f (0) = -u 0 (0). (1.12b)
Observe that the Cauchy problem (1.12) formally is equivalent to (1.1)- (1.3). Note that the existence and uniqueness results in half-space, with some appropriate conditions on u 0 and v 0 , were established in [START_REF] Lebeau | A wave problem in a half-space with a unilateral constraint at the boundary[END_REF]. The proof of Theorem 1 is rather classical. However for the reader convenience, this proof is given in the Appendix A.

Existence and uniqueness results by using the characteristic method

Theorem 1 (Existence and uniqueness results). Assume that p belongs to W 1,1 (0, L). Then the Cauchy problem (1.12) admits a unique absolutely continuous solution.

We introduce now some new notations: let X and ∂X be the interior and boundary of the set X, respectively, and let

I def = {t ∈ [0, L] : f (t) = 0} and J def = {t ∈ [0, L] : f (t) < 0}.
In the sequel, the notations for the constants introduced in the proofs are only valid in the proof. The aim of the next lemma is to prove further regularity results for the solution f of Problem (1.12).

Lemma 1 (Regularity result). Assume that p belongs to H 1 (0, L). Then the solution f to Problem (1.12) belongs to H 1 (0, L).

Proof. Note that Theorem 1 implies that f is bounded in W 1,1 (0, L). It follows that f (t) = 0 on I and f (t) = -2p (t) on J.

Clearly, we have

∫ J |f (t)| 2 dt = ∫ J 4|p (t)| 2 dt < ∞ and ∫ I |f (t)| 2 dt = 0.
Observe that if t is an accumulation point of ∂I, we may deduce that there exists a sequence t n belonging to I such that t n → t so that f (t n ) = 0 and f (t) = 0. Therefore, we deduce that f is bounded in L 2 (0, L), which implies that f belongs to H 1 (0, L). This concludes the proof.

It is convenient to define the following spaces:

H def = L 2 (0, L) and V def = {u ∈ H 1 (0, L) : u(L) = 0}
endowed with the norms • H and • V . Let (•, •) and a(•, •) the scalar products in H and V , respectively. This allows to define

H 2 def = {u ∈ L 2 (0, T ; V ) : u ∈ L 2 (0, T ; H)} endowed with the norm u H 2 def = ( ∫ T 0 ( u(•, t) 2 V + u(•, t) 2 H ) dt ) 1/2 ,
and •, • the duality corresponding product between H 2 and H 2 . We observe that H 2 → C 0 ([0, T ]; H) (see [START_REF] Simon | Compact sets in the space L p (0, T ; B)[END_REF]).

Chapter 1. Convergence of mass redistribution method for the one-dimensional wave equation with a unilateral constraint at the boundary Lemma 2. Assume that u belongs to H 2 and ü -u (defined in the sense of distributions on Q T ) is square integrable. Then we have for all ε > 0,

u ∈ C 0 ([0, T ]; H 1 (ε, L)) ∩ C 0 ([0, L]; H 1 (ε, T -ε)), u ∈ C 0 ([0, T ]; L 2 (ε, L)) and u ∈ C 0 ([0, L]; L 2 (ε, T -ε)).
Proof. The proof is obtained by using the same techniques detailed in [START_REF] Schatzman | Numerical approximation of a wave equation with unilateral constraints[END_REF]. Since it is quite a routine to adapt this proof to our case, we let the verification to the reader.

The aim of the next lemma is to obtain some further regularity results for the solution u to (1.4).

Lemma 3. Let u be the solution to (1.4). Then for all ε > 0,

u ∈ C 0 ([0, L]; H 1 (0, T -ε)) and u ∈ C 0 ([0, L]; L 2 (0, T -ε)).
Proof. The proof of this result exploits the local energy identity inside Q T , the reader can find a detailed proof in the Appendix as well as in [START_REF] Schatzman | Numerical approximation of a wave equation with unilateral constraints[END_REF] where a succinct proof is given.

We deal now with the energy balance. More precisely, we prove below that the energy associated to (1.12) given by

E(t) = 1 2 ∫ L 0 ( |u (t, x)| 2 +| u(t, x)| 2 ) dx is constant with respect to time t.
Lemma 4 (Energy balance). Assume that p is bounded in W 1,1 (0, L). Then the solution u to Problem (1.12) is energy conserving.

Proof. We observe first that (1.5) gives

E(t) = 1 2 ∫ L 0 ( (p (x+t)+q (x-t)) 2 +(p (x+t)-q (x-t)) 2 ) dx = ∫ L 0 |p (x+t)| 2 dx + ∫ L 0 |q (x-t)| 2 dx. (1.13)
We evaluate now the two integrals on the right hand side of (1.13). We note first that (1.7), (1.9) and (1.11) lead to

p (x+t) = { 1 2 ( u 0 (x+t)+v 0 (x+t) ) , 0 ≤ x + t ≤ L, 1 2
( u 0 (2L-(x+t))-v 0 (2L-(x+t)) ) , L ≤ x + t ≤ 2L, and
q (x-t) = { 1 2 ( u 0 (x-t)-v 0 (x-t) ) , 0 ≤ x -t ≤ L, ( f (-(x-t))+p (-(x-t)) ) , -L ≤ x -t ≤ 0.
On the one hand, we may deduce that

(1.14) ∫ L 0 |p (x+t)| 2 dx = 1 4 ∫ L-t 0 |u 0 (x+t)+v 0 (x+t)| 2 dx + 1 4 ∫ L L-t |u 0 (2L-x-t)-v 0 (2L-x-t)| 2 dx = 1 4
∫ L t |u 0 (x)+v 0 (x)| 2 dx + 1 4 ∫ L L-t |u 0 (x)-v 0 (x)| 2 dx.

The equivalence between the variational formulation and the differential inclusion

On the other hand, by using the same kind of arguments as in the proof of Lemma 1, we may obtain

|q (η)| 2 = |f (-η)+p (-η)| 2 = { |0+p (-η)| 2 a.e. on I, |-2p (-η)+p (-η)| 2 on J.
∫ L 0 |q (x-t)| 2 dx = ∫ t 0 |p (t-x)| 2 dx + 1 4 ∫ L t |u 0 (x-t)+v 0 (x-t)| 2 dx = 1 4
∫ t 0 |u 0 (t-x)+v 0 (t-x)| 2 dx + 1 4 ∫ L t |u 0 (x-t)+v 0 (x-t)| 2 dx = 1 4 ∫ t 0 |u 0 (x)+v 0 (x)| 2 dx + 1 4 ∫ L-t 0 |u 0 (x)-v 0 (x)| 2 dx.
Inserting (1.14) and (1.15) into (1.13), we get

E(t) = 1 4 ∫ L 0 |u 0 (x)+v 0 (x)| 2 dx + 1 4 ∫ L 0 |u 0 (x)-v 0 (x)| 2 dx,
and the conclusion is clear.

The equivalence between the variational formulation and the differential inclusion

The present section is dedicated to prove the equivalence between the weak formulation (1.4) and the differential inclusion (1.12). Consequently together with the results obtained in the previous section, it is possible to deduce that (1.4) possesses a unique solution.

Lemma 5. Assume that u 0 and v 0 belong to K def = {u ∈ V : u(0, •) ≥ 0} and H, respectively. Then the unique solution to Problem (1.12) defines a weak solution to (1.4) for T = L.

Proof. The idea of the proof consists to split the domain Q T into four regions according to Figure 1.2 and to use the expression of the solution on each region to show that u and u belong to L 2 (0, T ; V ) and to L 2 (0, T ; H), respectively (see Fig. 1

.2).

Let us go into the details. On the one hand, we observe that

x + t ∈ [0, L] and x -t ∈ [-L, 0] in the region I while x + t ∈ [L, 2L] and x -t ∈ [-L, 0] in the region II and we have u(x, t) = p(x+t) -p(t-x) -f (t-x), (1.16a) u (x, t) = p (x+t) + p (t-x) + f (t-x), (1.16b) u(x, t) = p (x+t) -p (t-x) -f (t-x). (1.16c) Since u 0 ∈ K, v 0 ∈ H and (1.7)-(1.9
) hold, we may infer that p belongs to H 1 (0, 2L) in the regions I and II. Besides, Lemma 1 implies that f belongs to H 1 (0, L). According to (1.16), we conclude that (1.17)

∫ L 0 ∫ L x ( |u(x, t)| 2 +|u (x, t)| 2 +| u(x, t)| 2 ) dt dx < +∞.
Chapter 1. Convergence of mass redistribution method for the one-dimensional wave equation with a unilateral constraint at the boundary On the other hand, we note that

I II III IV x t (L, L) (L, 0) (0, L)
x + t ∈ [L, 2L] and x -t ∈ [0, L] in the region III while x + t ∈ [0, L] and x -t ∈ [0, L] in the region IV and we have u(x, t) = p(x+t) + q(x-t), (1.18a) u (x, t) = p (x+t) + q (x-t), (1.18b) u(x, t) = p (x+t) -q (x-t). (1.18c)
Still using the fact that u 0 ∈ K, v 0 ∈ H and (1.7)-(1.9) hold, we may infer that p and q belong to H 1 (0, L) in the regions III and IV. Thanks to (1.18), we may infer that

(1.19) ∫ L 0 ∫ x 0 ( |u(x, t)| 2 +|u (x, t)| 2 +| u(x, t)| 2 ) dt dx < +∞.
Therefore, it follows from (1.17) and (1.19) that u and u are bounded in L 2 (0, T ; V ) and L 2 (0, T ; H), respectively. We deduce from Lemma 2 that u(0, •) and u (0, •) belong to L 2 (0, T ) and to L 2 (ε, T -ε) for all ε > 0, respectively, and u(•, 0) and u(•, 0) belong to L 2 (0, L) and to L 2 (ε, L) for all ε > 0, respectively. It remains to verify that (1.4) holds. To this aim, we observe by using the notations introduced above that

(1.20) - ∫ L 0 v 0 (v(x, 0)-u(x, 0)) dx - ∫ Q T u(x, t)( v-u)(x, t) dx dt + ∫ Q T u (x, t)(v -u )(x, t) dx dt = - ∫ L 0 v 0 (v(x, 0)-u(x, 0)) dx - ∫ Q T (p (x+t)-q (x-t))( v(x, t)-(p (x+t)-q (x-t))) dx dt + ∫ Q T (p (x+t)+q (x-t))(v (x, t)-(p (x+t)+q (x-t))) dx dt
for all v belonging to K such that there exists ζ > 0 with v = u for t ≥ L -ζ. We evaluate each integral on the right hand side of (1.20). Thanks to (1.6), we have

(1.21) ∫ L 0 v 0 (v(x, 0)-u(x, 0)) dx = ∫ L 0 (p (x)-q (x))(v(x, 0)-(p(x)+q(x))) dx.

The equivalence between the variational formulation and the differential inclusion

The second integral on the right hand side of (1.20) is integrated by parts with respect to t, we get

(1.22) ∫ Q T (p (x+t)-q (x-t))( v(x, t)-(p (x+t)-q (x-t))) dx dt = - ∫ L 0 (p (x)-q (x))(v(x, 0)-(p(x)+q(x))) dx - ∫ Q T (p (x+t)+q (x-t))(v(x, t)-(p(x+t)+q(x-t))) dx dt,
while the third one is integrated by parts with respect to x, we find

(1.23) ∫ Q T (p (x+t)+q (x-t))(v (x, t)-(p (x+t)+q (x-t))) dx dt = - ∫ T 0 (p (t)+q (-t))(v(0, t)-(p(t)+q(-t))) dt - ∫ Q T (p (x+t)+q (x-t))(v(x, t)- ( p(x+t)+q(x-t))) dx dt.
We substitute (1.21)-(1.23) into (1.20) and according to (1.8), we find

- ∫ T 0 (p (t)+q (-t))(v(0, t)-(p(t)+q(-t))) dt ≥ 0,
which implies that (1.4) holds.

Lemma 6. The weak solution (1.4) for T = L defines a solution to Cauchy problem (1.12).

Proof. Let u be a solution to (1.4), it follows from Lemma 2 that u has traces in L 2 (ε, L) on {t = 0} × (ε, L) and it comes that v 0 = u(•, 0) makes sense. We choose ψ = v -u such that ψ belongs to the space of infinitely differentiable functions on Q T with a compact support D(Q T ).

Then it follows that

- ∫ T 0 ( u, ψ) dt + ∫ T 0 a(u, ψ) dt = 0, for all ψ belonging to D(Q T ) where a(u, ψ) def = ∫ L 0 u ψ dx. This gives ü, ψ -u , ψ = 0 for all ψ belonging to D(Q T ).
Here the duality product between D (Q T ) and D(Q T ) is denoted by •, • . Therefore, we may deduce that ü -u vanishes in the sense of distributions in Q T . Thus we have

(1.24) ∫ Q T (-u ψ+u ψ ) dx dt = 0
for all ψ belonging to D(Q T ). We introduce the following notations: µ

+ def = x + t, µ - def = x -t and ∆ denotes the region bounded by the lines µ + = -µ -, µ + = 2L -µ -, µ -= µ + and µ -= µ + -2T in the plane (µ -, µ + ). Hence we have (1.25) u = ∂u ∂µ + - ∂u ∂µ - , u = ∂u ∂µ + + ∂u ∂µ - , ψ = ∂ψ ∂µ + - ∂ψ ∂µ - and ψ = ∂ψ ∂µ + + ∂ψ ∂µ - .
Chapter 1. Convergence of mass redistribution method for the one-dimensional wave equation with a unilateral constraint at the boundary Carrying (1.25) into (1.24), we find

- ∫ ∆ ( ∂u ∂µ + - ∂u ∂µ - )( ∂ψ ∂µ + - ∂ψ ∂µ - ) dµ -dµ + + ∫ ∆ ( ∂u ∂µ + + ∂u ∂µ - )( ∂ψ ∂µ + + ∂ψ ∂µ - ) dµ -dµ + = 0
for all ψ belonging to D(∆). We observe that (1.26)

∫ ∆ ∂u ∂µ + ∂ψ ∂µ - dµ -dµ + + ∫ ∆ ∂u ∂µ - ∂ψ ∂µ + dµ -dµ + = 0
for all ψ belonging to D(∆). The first term in (1.26) is integrated by parts with respect to µ - while the second one is integrated by parts with respect to µ + to get

- ∂ 2 u ∂µ + ∂µ - , ψ - ∂ 2 u ∂µ -∂µ + , ψ = 0,
for all ψ belonging to D(∆). Since we have

∂ 2 u ∂µ + ∂µ - , ψ = u, ∂ 2 ψ ∂µ + ∂µ - = ∂ 2 u ∂µ -∂µ + , ψ , it follows that ∂ 2 u ∂µ + ∂µ - , ψ = 0
for all ψ belonging to D(∆). Then we conclude that ∂ 2 u ∂µ + ∂µ -vanishes which holds if and only if u = p(µ + ) + q(µ -).

We observe that Lemma 2 and Lemma 3 imply that u(•, 0) = v 0 and u (0, •) belong to H and to L 2 (0, T -ε), respectively. According to Theorem 3 (given in Appendix B), the following Green's formulas make sense (1.27)

∫ T 0 ∫ L 0 (ü-u )ψ dx dt = - ∫ T 0 ( ( u, ψ)-a(u, ψ) ) dt + ∫ T 0 u (L, t)ψ(L, t) dt - ∫ T 0 u (0, t)ψ(0, t) dt + ( u(x, T ), ψ(x, T ) ) - ( v 0 , ψ(x, 0) )
for all ψ belonging to H 2 . We insert (1.27) into (1.4) and we choose v = u + ψ, we obtain

(1.28) - ∫ T 0 u (0, t)(v(0, t)-u(0, t)) dt ≥ 0
for all v belonging to K. Thanks to (1.5), we may deduce that (1.28) is equivalent to

(1.29) - ∫ T 0 (p (t)+q (-t))(v(0, t)-(p(t)+q(-t))) dt ≥ 0 for all v belonging to K. Since p(t) + q(-t) ≥ 0, it follows that v(0, t) = p(t) + q(-t) + α(t)
1.4. Convergence of mass redistribution method for all α(t) ≥ 0. Therefore we may infer from (1.29) that -

∫ T 0 (p (t)+q (-t))α(t) dt ≥ 0
for all α(t) ≥ 0, which implies that p (t) + q (-t) ≤ 0 for almost every t ∈ (0, T ). Finally we choose v(0, t) = 0 and v(0, t) = 2(p(t)+q(-t)) in (1.29), we get

∫ T 0 (p (t)+q (-t))(p(t)+q(-t)) dt = 0,
which allows us to infer that (p (t)+q (-t))(p(t) + q(-t)) vanishes for almost every t ∈ (0, T ). This concludes the proof.

Convergence of mass redistribution method

The semi-discretized problem by using finite elements is not well-posed which emphasizes some instabilities of time integration schemes (see [START_REF] Moreau | Liaisons unilatérales sans frottement et chocs inélastiques[END_REF][START_REF] Moreau | Nonsmooth mechanics and applications[END_REF][START_REF] Paoli | Schéma numérique pour un modèle de vibrations avec contraintes unilatérales et perte d'énergie aux impacts, en dimension finie[END_REF][START_REF] Paoli | Approximation et existence en vibro-impact[END_REF][START_REF] Khenous | Mass redistribution method for finite element contact problems in elastodynamics[END_REF]). In the literature many different approaches were elaborated to overcome this difficulty. For instance the uniqueness for an impact law of rigid bodies can be recovered by introducing a restitution coefficient (see [START_REF] Paoli | Time discretization of vibro-impact[END_REF]). However, this approach is not totally satisfactory for deformable bodies. Indeed the presence of oscillations due to displacement and to normal stress on the contact boundary induces some difficulties in the construction of energy conserving schemes (see [START_REF] Laursen | Design of energy conserving algorithms for frictionless dynamic contact problems[END_REF][START_REF] Paoli | Time discretization of vibro-impact[END_REF][START_REF] Laursen | Improved implicit integrators for transient impact problemsgeometric admissibility within the conserving framework[END_REF]). Another approach consists in using the penalty method which introduces some oscillations as well but which can be reduced with a damping technique (see [START_REF] Schweizerhof | Efficiency refinements of contact strategies and algorithms in explicit finite element programming[END_REF]). One of the key points to avoid oscillations is to use the mass redistribution method, the reader is referred to [START_REF] Khenous | Mass redistribution method for finite element contact problems in elastodynamics[END_REF] and the references therein.

We first approximate (1.1)-(1.3) by using the Lagrange affine finite element method. To this aim, we define h = L n where n is an integer and

V h def = {v h ∈ C 0 ([0, L]) : v h | [a i ,a i+1 ] ∈ P 1 , i = 0, . . . , n -1, v h (L) = 0}.
Here, a i def = ih, i = 0, 1, . . . , n, and P 1 is the space of polynomials of degree less than or equal to 1. A classical basis of V h is given by the sequence of shape functions ϕ i ∈ V h for i = 0, 1, . . . , n -1, defined by

ϕ i (x) def = { 1 -|x-a i | h if x ∈ [a max(i-1,0) , a i+1 ], 0 otherwise. Note that ϕ i (a j ) = δ ij , j = 0, 1, . . . , n -1, i.e. δ ij = 1 if i = j
and δ ij vanishes otherwise (δ is the Kronecker symbol). We approximate the solution u belonging to V to the weak formulation (1.4) by

u h (x, t) = n-1 ∑ j=0 u j (t)ϕ j (x).
Chapter 1. Convergence of mass redistribution method for the one-dimensional wave equation with a unilateral constraint at the boundary Consequently, we have u i = u h (a i ), i = 0, 1, . . . , n -1. The weak formulation (1.4) is approximated as follows

(P u h )                find u h : [0, T ] → V h and λ : [0, T ] → R such that for all v h ∈ V h ∫ L 0 üh v h dx + a(u h , v h ) = -λv h (0) a.e. t ∈ [0, T ], 0 ≤ u h (0, •) ⊥ λ ≤ 0 a.e. t ∈ [0, T ], u h (•, 0) = u 0h and uh (•, 0) = v 0h ,
where u 0h and v 0h belong to V h such that

(1.30) lim h→0 ( u 0h -u 0 V + v 0h -v 0 H ) = 0,
where λ is the Lagrange multiplier. We introduce now the following notations:

U def = (u 0 , . . . , u n-1 )
T and e 0 def = (1, 0, . . . , 0) T . The corresponding algebraic formulation is given by

(P U λ )            find U : [0, T ] → R n and λ : [0, T ] → R such that M Ü + SU = -λe 0 a.e. t ∈ [0, T ], 0 ≤ u 0 ⊥ λ ≤ 0 a.e. t ∈ [0, T ], U (0) = U 0 and U (0) = V 0 ,
where M and S denote the mass and stiffness matrices, respectively;

M ij def = ∫ L 0 ϕ i ϕ j dx and S ij def = a(ϕ i , ϕ j ) = 1 h             1 -1 0 • • • • • • 0 -1 2 -1 . . . . . . 0 -1 2 . . . . . . . . . . . . . . . . . . . . . . . . 0 . . . . . . -1 2 -1 0 • • • • • • 0 -1 1             , for all i, j ∈ [0, n -1].
We define now the modified mass matrix as follows:

M mod ij def = ∫ L h ϕ i ϕ j dx.
Clearly, we may observe that M 0i = M i0 = 0 for all i = 0, . . . , n -1. Therefore, the modified mass matrix reads

M mod def = ( 0 0 0 M ) .
Note that Mij = M i+1,j+1 for all i, j = 1, . . . , n -2. We introduce now the following notations:

Ū def = (u 1 , . . . , u n-1 ) T and Sij def = S i+1,j+1 with C def = ∫ Ω ϕ i+1 ϕ 0 dx, i = 0, . . . , n -2.
Observe that C = (S 10 , 0, . . . , 0) T . Thus by using the above notations, we have

S = ( S 00 C T C S ) and U = ( u 0 Ū ) .

Convergence of mass redistribution method

This leads to an algebraic formulation of the semi-discrete approximation with mass redistribution method given by

(P mod U λ )                  find U : [0, T ] → R n and λ : [0, T ] → R such that M Ü + S Ū = -Cu 0 a.e. t ∈ [0, T ], S 00 u 0 + C T Ū = -λ a.e. t ∈ [0, T ], 0 ≤ u 0 ⊥ λ ≤ 0 a.e. t ∈ [0, T ], U (0) = U 0 and U (0) = V 0 .
It follows that

u 0 = ( -λ-C T Ū S 00 ) = ( -λ-S 10 u 1 S 00
) .

If S 10 u 1 ≥ 0 then the compatibility condition gives u 0 = 0, so λ = (C T Ū ) -otherwise we have

λ = 0. This implies that u 0 = ( -S 10 u 1 S 00
) + , and then we may conclude that (P mod U λ ) is equivalent to the following second order Lipschitz continuous ordinary differential equation:

(P mod Ū )          find Ū : [0, T ] → R n-1 such that M Ü + S Ū = -C ( -S 10 u 1 S 00 ) + a.e. t ∈ [0, T ], U (0) = U 0 and U (0) = V 0 .
Lemma 7 (Existence and uniqueness results for (P mod U λ )). Problem (P mod U λ ) admits a unique solution (U, λ) which is Lipschitz continuous.

Proof. We use the fact that M is not a singular matrix as well as the same techniques detailed in [START_REF] Crouzeix | Analyse numérique des équations différentielles[END_REF] to establish that (P mod Ū ) possesses a unique Lipschitz continuous solution. On the other hand, we may deduce from (P mod U λ ) that u 0 = ( -S 10 u 1 S 00 ) + and λ = (C T Ū ) -. This allows us to conclude that u 0 and λ are also Lipschitz continuous and then the conclusion is clear.

We deal with the energy balance and we establish the energy conservation of the solution to problem (P mod U λ ) . More precisely, the discrete energy associated to problem (P mod U λ ) is given by (1.31)

E(t) = 1 2 ( U T M mod U +U T SU )(t). Lemma 8. The solution (U, λ) to problem (P mod U λ ) is energy conserving. Proof. We observe first that U T M mod Ü + U T SU = -U T λe 0 .
Therefore, we integrate this expression over (0, t) to get

∀t ∈ [0, T ] : E(t) -E(0) = - ∫ t 0 u0 (s)λ(s) ds. Let us define ω def = {t ∈ [0, T ] : u 0 (t) > 0}.
On the one hand, the contact conditions imply that λ = 0 on ω. On the other hand, the continuity of λ on [0, T ] gives that λ = 0 on ω where ω is the closure to ω. Furthermore, u0 vanishes in the interior of the set [0, T ]\ω. Hence u0 λ = 0 on [0, T ] and we conclude that

E(t) = E(0) for all t ∈ [0, T ].
Chapter 1. Convergence of mass redistribution method for the one-dimensional wave equation with a unilateral constraint at the boundary

We observe that (P mod U λ ) is equivalent to

(P mod u h )          find u h : [0, T ] → V h such that for all v h ∈ K h ∫ T 0 ( ∫ L h üh (v h -u h ) dx + a(u h , v h -u h ) ) dt ≥ 0, u h (•, 0) = u 0h and uh (•, 0) = v 0h .
We establish below the convergence of the solution u h of (P mod u h ) to the solution of (1.4) by using some ideas developed in [START_REF] Schatzman | Numerical approximation of a wave equation with unilateral constraints[END_REF].

Theorem 2. Assume that (1.30) holds. Then, the solution u h of (P mod u h ) converges strongly in H 2 to the unique solution of (1.4) as h tends to 0.

Proof. We observe that (1.32)

∫ L 0 | uh (x, t)| 2 dx = ∫ h 0 | uh (x, t)| 2 dx + ∫ L h | uh (x, t)| 2 dx.
We evaluate now the right hand side of (1.32). To this aim, we note that u 0 (t) = -S 10 S 00

( u 1 (t) ) + implies that (1.33) | u0 (t)| ≤ | u1 (t)|,
since -S 10 S 00 = 1. Therefore, by using (1.33) and Cauchy Schwarz's inequality, we may deduce that there exists

C 0 > 0 such that ∫ h 0 | uh (x, t)| 2 dx = ∫ h 0 | u1 (t)ϕ 1 (x)+ u0 (t)ϕ 0 (x)| 2 dx ≤ C 0 h| u1 (t)| 2 .
Furthermore, the energy conservation of Lemma 8 implies that there exists

C 1 > 0 such that ∫ L h | uh (x, t)| 2 dx ≤ C 1 .
Consequently, we deduce that

∫ (i+1)h ih | uh (x, t)| 2 dx = ∫ (i+1)h ih | ui+1 (t)ϕ i+1 (x)+ ui (t)ϕ i (x)| 2 dx ≤ C 1 for i = 1, . . . , n -1.
We conclude that for sufficiently small h, we get

| ui+1 (t)| 2 + | ui (t)| 2 ≤ 6C 1 h for i = 1, . . . , n -1. Therefore | u1 (t)| 2 ≤ 6C 1 h which implies that uh (•, t) L 2 (0,h
) is bounded independently of h and then uh (•, t) H is also bounded independently of h. By using Lemma 8, we can prove that u h V is bounded. It follows that there exists C > 0 independent of t such that sup

t∈[0,T ] ( u h (•, t) V + uh (•, t) H ) ≤ C.

Convergence of mass redistribution method

Then, it is possible to extract a subsequence, still denoted by u h , such that

u h u in L ∞ (0, T ; V ) weak * , (1.34a) uh u in L ∞ (0, T ; H) weak * . (1.34b) Let us define H ∞ def = {u ∈ L ∞ (0, T ; V ) : u ∈ L ∞ (0, T ; H)}
endowed with the following norm

u H∞ def = ess sup t∈[0,T ] ( u(•, t) V + u(•, t) H ) .
We may infer from (1.34) that

u h u in H ∞ weak * .
Since for all α < 1 2 the following injections [START_REF] Schatzman | Numerical approximation of a wave equation with unilateral constraints[END_REF]), where → is continuous embedding and → → is compact embedding, we get

H ∞ → C 0, 1 2 (Q T ) → → C 0,α (Q T ) hold (see
u h → u in C 0,α (Q T )
for all α < 1 2 . We observe that both u h (t) and u(t) belong to K. In order to prove that the limit u satisfies (1.4), it is necessary to choose convenient test functions. We approximate the elements of K before projecting them onto V h . Indeed, the L 2 projection does not conserve the constraint at x = 0, and therefore, the elements of K need another approximation in order to satisfy the constraint strictly. To this aim, let v be an element of K which is equal to u for t ≥ T -ε and let

v β (x, t) def =      u(x, t) + 1 β ∫ t+β t (v(x, s)-u(x, s)) ds + k(β)(L-x)g(t) if t ≤ T -β, u(x, t) if t ≥ T -β, for all β ≤ ε 4 .
Here g(t) is a smooth, positive function (see Figure 1.3) and satisfying

g(t) 1 t 0 T -ǫ 2 T -ǫ 4 T Figure 1.3: Smooth and positive function g(t). g(t) def = { 1 if t ∈ [0, T -ε 2 ], 0 if t ∈ [T -ε 4 , T ].
Chapter 1. Convergence of mass redistribution method for the one-dimensional wave equation with a unilateral constraint at the boundary The next step consists to choose adequately k(β). Since u belongs to C 0, 1 2 (Q T ), we may deduce that there exists C > 0, such that

u(0, t) -β -1 ∫ t+β t u(0, s) ds ≤ β -1 ∫ t+β t |u(0, t)-u(0, s)| ds ≤ C u H∞ β -1 ∫ β 0 √ s ds = 2 3 C u H∞ √ β.
Furthermore, we have the following inequality

v β (0, t) ≥ β -1 ∫ t+β t v(0, s) ds - 2 3 C u H∞ √ β + k(β)Lg(t)
for all t ≤ T -ε 2 . Choosing k(β) = 5 3L C u H∞ √ β, we get v β (0, t) ≥ C u H∞ √ β for all t ≤ T -ε 2 .
On the other hand, we have

v β (x, t) = u(x, t) + k(β)(L-x)g(t) for all t belonging to [T -ε 2 , T -β]. Hence v β belongs to K and in the other hand v β belongs to L ∞ (0, T ; V ) because v β (•, t)-u(•, t) V ≤ k(β)Lg(t) + 1 √ β v-u H 2 .
We denote by Q h the orthogonal projection onto V h with respect to the scalar product in H such that Q h z-z V → 0 when h → 0 for all z ∈ V (see [START_REF] Ern | Theory and practice of finite elements[END_REF]). The Sobolev injections imply that there exists a sequence α h converging to zero as h tends to zero such that

Q h z-z C 0 ≤ α h z V ,
for all z belonging to V with lim h→0 α h = 0. The test function is defined as follows:

(1.35) v h (•, t) def = u h (•, t) + Q h (v β -u)(•, t)
for all t belonging to [0, T ]. By using a continuity argument, v h (0, •) ≥ 0 for h small enough. Carrying (1.35) into (P mod u h ) and using the integration by parts, we find

- ∫ L 0 uh (•, 0)Q h (v β -u)(•, 0) dx + ∫ h 0 uh (•, 0)Q h (v β -u)(•, 0) dx - ∫ T 0 ∫ L 0 uh (•, t)Q h ( vβ -u)(•, t) dx dt + ∫ T 0 ∫ h 0 uh (•, t)Q h ( vβ -u)(•, t) dx dt + ∫ T 0 a(u h (•, t), Q h (v β -u)(•, t) dt ≥ 0. Since (v β -u)(•, t) is bounded in H ∞ ,
the above integration makes sense. Thus we may pass to the limit when h tends to zero. Since we have

Q h ( vβ -u) → ( vβ -u) in L 2 (0, T ; H) and Q h (v β -u) → (v β -u) in L 2 (0, T ; V ).

Numerical examples

Then, we conclude that

- ∫ L 0 v0 (v β -u)(•, 0) dx - ∫ L 0 ∫ T 0 u(•, t)( vβ -u)(•, t) dt dx + ∫ T 0 a(u(•, t), (v β -u)(•, t)) dt ≥ 0.
We pass now to the limit with respect to β so we obtain variational formulation (1.4).

On the one hand, we observe that Lemma 4 leads to (1.36)

∫ L 0 (| u(•, t)| 2 +|u (•, t)| 2 ) dx = ∫ L 0 (|v 0 | 2 +|u 0 | 2 ) dx.
On the other hand, Lemma 8 implies that

∫ L h | uh (•, t)| 2 dx + a(u h (•, t), u h (•, t)) = ∫ L h | uh (•, 0)| 2 dx + a(u h (•, 0), u h (•, 0)),
which by using (1.30) and (1.34) gives

(1.37) lim h→0 ∫ L h (| uh (•, t)| 2 +|u h (•, t)| 2 ) dx = ∫ L h (|v 0 | 2 +|u 0 | 2 ) dx.
Therefore from (1.36) and (1.37)

lim h→0 ∫ L h (| uh (•, t)| 2 +|u h (•, t)| 2 ) dx = ∫ L 0 (| u(•, t)| 2 +|u (•, t)| 2 ) dx.
Since u h converges weakly to u in H ∞ and u h H∞ converges to u H∞ and since H ∞ → H 2 then we conclude that u h converges strongly to u in H 2 .

Numerical examples

We perform a finite element discretization in space and we use a classical Newmark time stepping method. This leads to consider the following problem:

(1.38)                  U n+1 = U n + ∆t U n + ( 1 2 -β ) ∆t 2 Ü n + β∆t 2 Ü n+1 , U n+1 = U n + (1-γ)∆t Ü n + γ∆t Ü n+1 , M Ü n+1 + SU n+1 = -λ n+1 e 0 , 0 ≤ u n+1 0 ⊥ λ n+1 ≤ 0, U (0) = U 0 and U (0) = V 0 ,
where β ∈]0, 1/2[ and γ ∈]1/2, 1[ are the classical parameters of the Newmark scheme. Note that if β = 0.25 and γ = 0.5, the scheme (1.38) is the so-called Crank-Nicholson scheme which is an implicit, unconditionally stable and second-order accurate scheme for elastodynamic problems without contact conditions and moreover it is energy conserving (see [START_REF] Khenous | Mass redistribution method for finite element contact problems in elastodynamics[END_REF]). On the other hand, it is well known that the space-semi discretization of contact problems in elastodynamics presents some numerical instabilities (see [START_REF] Hughes | A finite element method for a class of contact-impact problems[END_REF]) which can be avoided by using a modified mass method (see [START_REF] Khenous | Mass redistribution method for finite element contact problems in elastodynamics[END_REF] and the references therein). We make below some comparisons between two different Chapter 1. Convergence of mass redistribution method for the one-dimensional wave equation with a unilateral constraint at the boundary approaches; the one using a standard mass matrix and the one using a modified mass matrix. The parameters used in the numerical simulations are the space step ∆x = 0.1, the time step ∆t = 0.01, the initial displacement u 0 (x) = 0.5x -0.5, the initial velocity v 0 (x) = 0 and the Dirichlet value u(L, t) = 0.45 with L = 1 and T = 5. The numerical experiments are performed by employing the finite element library Getfem++ (see [START_REF] Renard | Getfem++. An Open Source generic C++ library for finite element methods[END_REF]). In particular, the generalized Newton algorithm has been used to compute the unique solution of (1.38) (see [START_REF] Renard | Generalized Newton's methods for the approximation and resolution of frictional contact problems in elasticity[END_REF][START_REF] Curnier | A generalized Newton method for contact problems with friction[END_REF]). The numerical results show that when the constraint is active, small oscillations occur in the case where M is a standard mass matrix (see Figures 1.4 We present now some numerical results obtained for an undeformed elastic bar which is dropped with some initial velocity on a rigid obstacle and we compare the exact solution to the approximated one obtained by using the mass redistribution method. More precisely, we assume that this bar falls from a height u 0 , with an initial velocity -v 0 and under the gravity g ≥ 0. Furthermore both ends of the bar are free to move, as long as the bar does not hit a rigid obstacle. The length and the Young modulus of the bar are denoted by L and E, respectively. Let u(x, t) be the displacement at time t of the material point of spatial coordinate x ∈ [0, L] and the contact pressure equal to the normal stress -Eu (0, •). Then the mathematical problem can be formulated as follows:

(1.39) ü(x, t) -Eu (x, t) = -(g+Eu (0, t)), (x, t) ∈ (0, L) × (0, T ),
with Cauchy initial data

(1.40) u(•, 0) = u 0 and u(•, 0) = -v 0 ,
and Signorini and Neumann boundary conditions at x = 0 and x = L, respectively, for t > 0

(1.41) 0 ≤ u(0, •) ⊥ Eu (0, •) ≤ 0 and u (L, •) = 0.
The existence and uniqueness results for (1.39)-(1.41) is obtained by rewriting this problem as a differential inclusion problem and then by using the same techniques detailed in the proof of Theorem 1. Since it is quite a routine to adapt this proof to the case considered here, we leave the verification to the reader. In order to calculate the analytical solution to problem (1.39)-(1.41), we distinguish three phases, namely before the contact, during the contact and after the contact. To this aim, we choose v 0 = 0 and g > 0 so that the bar can make several impacts.

The bar reaches the rigid obstacle at time

t 1 = √ 2u 0
g with the velocity equal to √ 2u 0 g. After the impact, the bar stays in contact with the rigid body; as soon as a shock wave travels from bottom of the bar to the top and vice versa then the bar takes off. The wave takes a time L √ E Chapter 1. Convergence of mass redistribution method for the one-dimensional wave equation with a unilateral constraint at the boundary to travel along the bar. The impacts occur at time

t 4k+1 = 3 L √ E + 16k L √ E , t 4k+2 = t 4k+1 + 2 L √ E , t 4k+3 = t 4k+1 + 8 L √ E , t 4k+4 = t 4k+1 + 10 L √ E .
We introduce also the following notations:

h 1 (x, t) = - √ 2u 0 g min ( x √ E , L √ E -t- L √ E ) + ∞ ∑ n=1 2g ELν 3 n (cos( √ Eν n t)-1) sin(ν n x), h 2 (x, t) = u 0 - 1 2 g ( t- √ 2u 0 g ) 2 - 2gL 2 3E + ∞ ∑ n=1 4g Eλ 2 n cos( √ Eλ n t) cos(λ n x), with ν n = (n-1 2 ) π L and λ n = n π L .
Then, the explicit solution reads as

(1.42) u(x, t) =                  u 0 -1 2 gt 2 if t ≤ t 1 , h 1 (x, t-t 4k+1 ) if t 4k+1 < t ≤ t 4k+2 , h 2 (x, t-t 4k+2 ) if t 4k+2 < t ≤ t 4k+3 , h 1 (x, t 4k+4 -t) if t 4k+3 < t ≤ t 4k+4 , u 0 -1 2 g ( t-t 4k+4 - √ 2u 0 g ) 2 if t 4k+4 < t ≤ t 4(k+1)+1 .
Here, some details are omitted, the reader is referred to [START_REF] Wriggers | Computational contact mechanics[END_REF][START_REF] Doyen | Time-integration schemes for the finite element dynamic Signorini problem[END_REF] for a detailed explanation. We choose L = 10, T = 6.5, E = 900, g = 10, the initial data u 0 (x) = 5, v 0 (x) = 0 and Neumann value u (L, t) = 0. The Newmark time stepping method with β = 0.25 and γ = 0.5 is used to evaluate the approximated solution. Let us emphasize that if the space step ∆x and time step ∆t tend to 0, the approximated solution obtained by using the mass redistribution method (P mod U λ ) converges to the solution of (1.39)-(1.41) explicitly given by (1.42) (see Figure 1.7 (left)). On the other hand, we can write at least formally an energy relation for (1.39): we multiply this equation by u, we integrate by parts over Q τ , τ ∈ [0, T ], we get

1 2 ∫ L 0 | u(•, τ )| 2 dx + 1 2 ∫ L 0 | √ Eu (•, τ )| 2 dx = - ∫ Qτ g u dx dt
for all τ belonging to [0, T ]. Observe that the energy tends to be conserved when the space step ∆x and the time step ∆t tend to 0 and the energy decreases otherwise (see Figure 1.7 (right)).

Numerical examples

Appendix A

The aim of this section is to give the proofs of Theorem

F (t, y) = -J N (y) -2p (t).
The multivalued map F : [0, L] × (-∞, 0] → P( R) has closed convex values and F is measurable with respect to its second variable. We prove now that F is upper semi-continuous with respect to its second variable which is equivalent to establish that

J N is upper semi-continuous. Note that if A ⊂ R, J -1 N (A) = [0, +∞) or J -1 N (A) = {0} or J -1 N (A) = ∅ which are closed sets. According to Definition [104, p.4, Def. 1.1], J N is upper semi-continuous.
We verify that there is a function

r(t, y) = c(t)(1+|y|) with c ∈ L 1 (0, L) such that F (t, y) ∩ r(t, y) B1 (0) ∩ T (-∞,0] (y) = ∅ on [0, L] × (-∞, 0],
where B 1 (0) is the ball of radius 1 at the origin and T (-∞,0] (y) is the tangent cone on y, the reader is referred to [START_REF] Deimling | of de Gruyter Series in Nonlinear Analysis and Applications[END_REF] for further details. Indeed, we distinguish two cases, on the one hand, if y belongs to the interior of (-∞, 0], T (-∞,0] (y) = R and F (t, y) = -2p (t), we choose c(t) = 1 + |2p (t)| and on the other hand, if y = 0, T (-∞,0] (0) = (-∞, 0] and F (t, y) = (-∞, -2p (t)], we choose c(t) = 1 + |2p (t)|. Therefore the existence of solution to (1.12) follows.

The uniqueness result comes from the monotonicity of J N , namely J N is the subdifferential of a convex, lower semi-continuous and proper function, the reader is referred to [START_REF] Brézis | Opérateurs maximaux monotones et semi-groupes de contractions dans les espaces de Hilbert[END_REF] for further details.

Chapter 1. Convergence of mass redistribution method for the one-dimensional wave equation with a unilateral constraint at the boundary Proof of Lemma 3. We note that

(1.43) (u 2 + u2 ) -2 ∂ ∂t (u u) = 0 in the sense of distributions.
Hence we integrate (1.43) over

[x 0 , x 1 ] × [t 0 , t 1 ], with 0 < x 0 < x 1 < L and 0 < t 0 < t 1 < T , to get (1.44) ∫ t 1 t 0 ((u 2 + u2 )(x 1 , t)-(u 2 + u2 )(x 0 , t)) dt = ∫ x 1 x 0 ((2u u)(x, t 1 )-(2u u)(x, t 0 )) dx.
According to Lemma 2, the right hand side of (1.44) is bounded independently of x 0 , x 1 , t 0 , t 1 as long as 0 < x0 x 0 < x 1 < L. We integrate now (1.44) with respect to x 0 over [x 0 , L], we may deduce that there exists C > 0 independent of x 1 , x0 such that

(L-x 0 ) ∫ t 1 t 0 (u 2 + u2 )(x 1 , t) dt ≤ ∫ L x0 ∫ t 1 t 0 (u 2 + u2 )(x 0 , t) dt dx 0 + C(L-x 0 ), which implies that x → ∫ t 1 t 0 (u 2 + u2 )(x 1 , t) dt is bounded on [x 0 , L] independently of t 0 and t 1 , it follows that x → ∫ T 0 (u 2 + u2 )(x, t) dt is bounded on [x 0 , L]. Let v be the solution of the following problem (1.45)        v -v = 0 on (0, x 0 ) × (0, T ), v(x 0 , t) = u(x 0 , t) and v (x 0 , t) = u (x 0 , t) for all t ∈ [0, T ], v(x, 0) = u 0 (x) and v(x, T ) = u(x, T ) for all x ∈ [0, x 0 ].
These conditions are illustrated in Figure 1.8. Since u (x 0 , •) and u(x 0 , •) belong to L 2 (0, T ),

0 T v(x, T ) = u(x, T ) v(x 0 , t) = u(x 0 , t) v ′ (x 0 , t) = u ′ (x 0 , t) v(x, 0) = u 0 (x) v(x, 0) = v 0 (x)
x 0 we may infer that there exists a unique solution to (1.45). More precisely, w = u -v satisfies (1.45) with homogeneous boundary initial conditions and the existence and uniqueness theorem in [START_REF] Ball | Strongly continuous semigroups, weak solutions, and the variation of constants formula[END_REF] holds. Furthermore v = u on (0, x 0 ) × (0, T ) and in particular we have v(x, 0) = v 0 (x).

We solve (1.45) by employing a classical characteristic method. To this aim, it is convenient to introduce the following notations:

β 1 def = -x + t and β 2 def = -x -t.
We may deduce that ∂ 2 v ∂β 1 ∂β 2 vanishes which implies that v(x, t) = f (β 1 ) + g(β 2 ). Notice that the general solution for all of points in the rectangle (0, x 0 ) × (0, T ) does not exist. Then we split the 1.5. Numerical examples rectangle into three regions by using characteristic lines as it is shown on Figure 1.9. We look for the solution in each region. More precisely, in region I and according to the initial condition of problem (1.45) in x = x 0 for x ≤ x 0 , we get

(1.46) v(x, t) = 1 2 (u(x 0 , t+(x 0 -x))+u(x 0 , t-(x 0 -x))) - 1 2 ∫ t+(x 0 -x) t-(x 0 -x) u (x 0 , ζ) dζ.
Observe that Figure 1.9 gives a better interpretation of this phenomenon. Indeed, the interval used will be the intersection of the line x = x 0 with the forward wave cone at (x, t) which is the region between the two straight lines having a slope of ±1 but directed upwards from an origin (x, t). The forward wave cone at the point (x, t) will enclose all those points (x 0 , ζ) which motion will be influenced by what occurred at the point x at the time t. Concerning the region II, we

0 T v(x, T ) = u(x, T ) v(x 0 , t) = u(x 0 , t) v ′ (x 0 , t) = u ′ (x 0 , t) v(x, 0) = u 0 (x) v(x, 0) = v 0 (x)
T v(x, T ) = u(x, T ) v(x 0 , t) = u(x 0 , t) v ′ (x 0 , t) = u ′ (x 0 , t) v(x, 0) = u 0 (x) v(x, 0) = v 0 (x) x 0 I II III A Figure 1.10: Characteristics in the region II. It follows that v (x A , t A ) -v(x A , t A ) = u 0 (x A +t A ) -v 0 (x A +t A ), (1.47a) v (x A , t A ) + v(x A , t A ) = u (x 0 , t A +(x 0 -x A )) + u(x 0 , t A +(x 0 -x A )). (1.47b) Therefore (1.47) leads to (1.48) v (x A , t A ) = 1 2 (u 0 (x A +t A )-v 0 (x A +t A )+u (x 0 , t A +(x 0 -x A ))+ u( x 0 , t A +(x 0 -x A ))), and
v(x A , t A ) = 1 2 (-u 0 (x A +t A )+v 0 (x A +t A ) + u (x 0 , t A +(x 0 -x A ))+ u(x 0 , t A +(x 0 -x A ))).
Chapter 1. Convergence of mass redistribution method for the one-dimensional wave equation with a unilateral constraint at the boundary

We have the solution for all the points located in the regions I and II. Concerning the solution in region III, we need some further regularity result to conclude. We obtained some regularity results for u(x 0 , t), u (x 0 , t) and u(x 0 , t) in Lemma 2. Besides by using (1.48) and (1.46), it is possible to deduce that v belongs to C 0 ([0, x 0 ]; L 2 (0, T -x 0 )) which implies that v belongs to C 0 ([0, x 0 ]; H 1 (0, T -x 0 )). As x 0 is arbitrary small, the conclusion is clear.

Appendix B

We establish below a Green's formula that is crucial in the proof of Lemma 6. To this aim, let us introduce a linear topological space D and Hilbert spaces V, H, Z and S with their topological duals denoted by D , V , H , Z and S .

Suppose that D is contained in V and it is dense in H. Here H is identified with its dual, namely we have H = H . Furthermore, Suppose that V is contained in H with finer topology and we denote by V 0 , the closure of D in V such that

D → V 0 → H = H → V 0
with dense embedding. More precisely, D is an abstraction of the usual space D(Q T ) of test function. In our case the spaces V, S and Z denote the admissible displacement, stresses and boundary values containing traces of element of V, respectively. We also introduce a linear operator A ∈ L(V, S) and its restriction to V 0 denoted by

A 0 ∈ L(V 0 , S) such that ∀v ∈ V 0 , Av = A 0 v.
Let A * ∈ L(S , V ) be the adjoint of the operator A, defined by

τ, Av S ×S = A * τ, v V ×V . Then ∀τ ∈ S , ∀v ∈ V 0 : τ, Av S ×S = A * 0 τ, v V ×V . Let us define S (A * 0 ) def = {τ ∈ S : A * 0 τ
∈ H}, then we have the following trace theorem.

Theorem 3. Suppose that V 0 is the kernel of a surjective map γ ∈ L(V, Z) from V onto Z.
Then there exists a unique linear operator π ∈ L(S (A * 0 ), Z ) such that the following Green's formula holds:

∀τ ∈ S (A * 0 ), ∀v ∈ V : τ, Av S ×S -A * 0 τ, v V ×V = πτ, γv Z ×Z .
The detailed proof of (3) is given in [START_REF] Aubin | Approximation of elliptic boundary-value problems[END_REF]. In particular, we are interested in this work to the case where

S = S def = {(u 1 , u 2 ) ∈ L 2 (Q T )}, D def = D(Q T ), V def = H 1 (Q T ), V 0 = H 1 0 (Q T ), V 0 = H -1 (Q T ), H def = L 2 (Q T ), S (A * 0 ) = {(u 1 , u 2 ) ∈ S : ∂ ∂t u 1 -∂ ∂x u 2 ∈ L 2 (Q T ) in the sense of distributions}, A : u → ( ∂ ∂t u, -∂ ∂x u), Z def = H 1/2 (∂Q T ) and Z = H -1/2 (∂Q T )
, where (∂Q T ) is the boundary of Q T and the trace operator γ : V → Z. Then there exists a unique π ∈ L(S

(A * 0 ), Z ) such that ∫ Q T ( τ 1 ∂ ∂t v + τ 2 ∂ ∂x v ) dt dx - ∫ Q T ( ∂ ∂t τ 1 - ∂ ∂x τ 2 ) v dt dx = πτ, γv Z ×Z , 1.5. Numerical examples for all v ∈ V and (τ 1 , τ 2 ) ∈ S (A * 0 ).
Then by density argument, we have

πτ, γv Z ×Z = ∫ ∂Q T (τ.n)v ds,
where n is the outward unit normal to Q T , when τ and v are regular enough.

Appendix C

Multivalued differential equations

In this Section, we recall some basic results for multivalued differential equations, the reader refer to [START_REF] Barbu | Nonlinear semigroups and differential equations in Banach spaces[END_REF][START_REF] Deimling | of de Gruyter Series in Nonlinear Analysis and Applications[END_REF] for further details. A multivalued map is a map which values are subsets of a vector space X. Often these are particular subsets, for example closed convex subsets. These are seen as extensions of ordinary functions on X, and elementary operations on multivalued map are deduced from the elementary operations on X.

Definition 1. Let Xand Y be two Banach spaces. Then we denote P(X) the set of all subsets of X. Given another set D = ∅ a subset of Y , a map F : D → P(X)\{∅} is called a multivalued map.

A single-valued map means

F (ω) = {f (ω)} on D. A function f : D → X with f (ω) ∈ F (ω) is called a selection of F .
We describe now the main elementary operations on the multivalued maps. Firstly, for two subsets A, B of X, we define

A + λB = {a + λb : a ∈ A, b ∈ B} and x + B = {x} + B.
So, for example, B r (x 0 ) = x 0 + rB 1 (0) and

A + B r (0) = {x ∈ X : ρ(x, A) ≤ r} with ρ(x, A) = inf a∈A |x -a|.
The addition between the multivalued maps is defined as follows:

(F + λG)(x) = F (x) + λG(x) = {f + λg : f ∈ F (x), g ∈ G(x)}.
Hausdorff distance between two sets A and B is defined by

d H (A, B) = max { sup x∈A ρ(x, B), sup x∈B ρ(x, A) } .
For A ⊂ X we denote by

A = sup y∈A y X . The convex hull of A ⊂ X is the smallest convex set containing A, i.e. convA = { m ∑ i=1 λ i x i : m ∈ N, λ i ∈ [0, 1], m ∑ i=1 λ i = 1, x i ∈ A for i = 1, . . . , m } .
The property of upper semicontinuity property of multivalued maps plays an important role to obtain some existence results for differential inclusions.

Chapter 1. Convergence of mass redistribution method for the one-dimensional wave equation with a unilateral constraint at the boundary Definition 2. Let X, Y be Banach spaces and ∅ = D ⊂ Y , a map F : D → P(X)\{∅} is said to be upper semi-continuous (usc), if

F -1 (A) is closed in D whenever A ⊂ X is closed.
The condition of being usc may become more transparent in terms of sequences since it means:

if (ω n ) ⊂ D, A ⊂ X is closed, ω n → ω 0 ∈ D and F (ω n ) ∩ A = ∅ for all n ≥ 1 then also F (ω 0 ) ∩ A = ∅.
So, intuitively, F (ω 0 ) must be at least as large as the limit of the F (ω n ).

A prototype of an usc multivalued map, defined by

F (t) =      {0} for t < 0 {0, 1} for t = 0 {1} for t > 0
Later, it will often be essential to have convex values, here, this can simply be achieved by letting F (0) = [0, 1], i.e. by filling in the gap at the point of discontinuity. This situation can easily be generalized, given f : D → X, by

F (x) = ∩ δ>0 convf (B δ (x) ∩ D) for x ∈ D.
Obviously, f is a selection of F , F has closed convex values, and

F (ω 0 ) = {f (ω 0 )} if f is continuous at ω 0 . However, if f is locally compact (i.e. for all ω 0 ∈ D there exists r = r(ω 0 ) > 0 such that f (B r (ω 0 ) ∩ D) is compact) then F (x)
is usc (see [START_REF] Deimling | of de Gruyter Series in Nonlinear Analysis and Applications[END_REF]).

Definition 3. [104, p.35 Sec 4.3] Let X be a real Hilbert space with inner product (•, •), a map F : D ⊂ X → P(X)\{∅} is said to be monotone if

(F x -F y , x -y) ≥ 0 on D × D,
which is understood as (u -v, x -y) ≥ 0 for all u ∈ F x and all v ∈ F y . In case -F is monotone we say that F is dissipative. Note that G δ -set is the name for intersections of countably many open sets.

Definition 4. Let X be a Hilbert space with inner product (•, •), Such a solution exists if F does not grow too fast with respect to x, namely

F : D ⊂ X → P(X)\{∅} is said semi-Lipschitz if there is k ∈ R such that (F (x) -F (y), x -y) ≤ k x -y 2 X on D × D. Given X = R n , J = [0, T ] ⊂ R, be closed D ⊂ X,
F (t, x) = sup{ y X : y ∈ F (t, x)} ≤ c(t)(1 + x X ) on J × D, with c ∈ L 1 (J) and F (t, x) ∩ T D (x) = ∅ on J × D,
where

T D (x) is a tangent cone on X T D (x) = { y ∈ X : lim λ→0 + λ -1 ρ(x + λy, D) = 0 } for x ∈ D.
We introduce the concept of measurability in multivalued maps Definition 5. Let (D, A) be a measurable space, a multivalued map F :

D → P(R n )\{∅} is called measurable if F -1 (B) ∈ A for all the open sets B ⊂ R n .
There are some existence results of Cauchy problem (1.49). We present some of them. The basic existence result for autonomous problems i.e. f independent of t is the following Lemma:

(a) F (t, x) ∩ T D (x) = ∅ on J × D, (b) F (t, x) ≤ c(t)(1 + x X ) on J × D, for c ∈ L 1 (J), (c) F (t, •) is usc for all t ∈ J, (d ) F (•, x) is measurable for all x ∈ D.
Then for any initial u 0 ∈ D, the Cauchy problem (1.49) admits at least one absolutely continuous solution on J.

Another important result

(a) F (t, x) ∩ r(t, x)B 1 (0) ∩ T D (x) = ∅ on J × D with r(t, x) = c(t)(1 + x X ). and c ∈ L 1 (J), (b) F (t, •) is usc for all t ∈ J, (c) F (•, x) is measurable for all x ∈ D.
Then for any initial u 0 ∈ D, the Cauchy problem (1.49) admits at least one absolutely continuous solution on J.

After proving the existence, uniqueness was investigated. Let X be a real Hilbert space with inner product (•, •) and

F : J × D ⊂ X → P(X)\{∅} be monotone. If there exist k ∈ L 1 (J) such that (1.50) ( F (t, x) -F (t, y), x -y ) ≤ k(t) x -y 2 X , on J × D × D,
i.e. F be a semi-Lipschitz function. This ensures the uniqueness of the solution.

Theorem 5. [104, p. 138 Thm. 10.5] Let X be a real Hilbert space, J = [0, T ] ⊂ R, and F : J × D → P(X)\{∅} have closed convex values and be such that

(a) F (t, •) is usc, (b) F (•, x) has a strongly measurable selection, (c) condition ( F (t, x) -F (t, y), x -y ) ≤ ω(t, x -y X ) x -y X holds with ω(t, ρ) = k(t)ρ and k ∈ L 1 (J), (d ) F (t, x) ≤ c(t)(1 + x X ) on J × D, for c ∈ L 1 (J),
are satisfied. Then Cauchy problem (1.49) has a unique solution on J.

Some functional results

We recall some functional results. First, we start with some definitions and results on semigroups. Let X be a Banach space with norm • . We denote by L(X) the Banach algebra of all bounded linear operators on X endowed with the operator norm, still denoted by • . Definition 6. A family (S(t)) t 0 of bounded linear operators on a Banach space X is called a (one-parameter) semi-group on X if it satisfies in the functional equation

(i) S(t + s) = S(t)S(s)
for all t, s ≥ 0,

(ii) S(0) = I.
If (i) and (ii) hold even for all t, s ∈ R, we call (S(t)) t∈R a (one-parameter) group on X.

Definition 7. A one-parameter semi-group (S(t)) t 0 on a Banach space X is called uniformly continuous (or norm continuous ) if

R + t → S(t) ∈ L(X),
is continuous with respect to the uniform operator topology on L(X).

Numerical examples

The linear operator A defined by

D(A) = { x ∈ X : lim t→0 S(t)x -x t exists } ,
and

Ax = lim t→0 S(t)x -x t = d + S(t)x dt | t=0 for x ∈ D(A)
is the infinitesimal generator of the semi-group S(t). D(A) denotes the domain of A.

(A + λ) -k M (λ -β) -k , λ > β, k 1.
The set of all operators A satisfying the conditions (i) and (ii) will be denoted by g(M, β). Let {S(t)} t 0 be a unitary semi-group and S(t) 1, S(t) is called a contraction semi-group. If A ∈ g(1, 0) then the operator -A is the infinitesimal generator of a semi-group which can be written as S(t) = e -tA (see [START_REF] Kato | Perturbation theory for linear operators[END_REF]).

For u 0 ∈ D(A), there exists a unique function

u in C 1 ([0, T ]; X) ∩ C 0 ([0, T ]; D(A)) such that (1.51) du dt (t) + Au(t) = 0 for all t ∈ (0, T ),
with initial condition

u(0) = u 0 , then u(t) = S(t)u 0 .
This unique function u is called a strong solution. For more details we refer the reader to [START_REF] Dautray | Analyse mathématique et calcul numérique pour les sciences et les techniques[END_REF].

Let A ∈ g(M, β). The semi-group S(t) = e -tA gives the solution of the differential equation (1.51) in the form u(t) = S(t)u(0). Let us consider the inhomogeneous differential equation for f ∈ C 0 (0, T ; X)

(1.52)    du dt (t) + Au(t) = f (t) for all t ∈ (0, T ), u(0) = u 0 .
The solution can be expressed with S(t) and depends linearly on u 0 and f . We have

d ds (S(t -s)u(s) = S(t -s)f (s))
Integrating this expression over (0, t), we obtain

(1.53)      u(t) = S(t)u 0 + ∫ t 0 S(t -s)f (s) ds, u(0) = u 0 .
In particular this implies that the solution of (1.52) is uniquely determined by u(0). Note that if u is the unique solution to (1.52) and if f ∈ C 0 ([0, T ]; X), then Au(t) ∈ C 0 ([0, T ]; X) and the equation (1.52) is verified in C 0 ([0, T ]; X). For more details we refer the reader to [START_REF] Kato | Perturbation theory for linear operators[END_REF]).

If S * (t) is defined by

S(t)u, u * = u,
d dt u(t), u * = u(t), A * u * + f (t), u * , for almost all t ∈ [0, T ].
For the existence and uniqueness of this weak solution, we have the following Theorem:

Theorem 8. [105]
For each u 0 ∈ X there exists a unique weak solution u(t) of (1.52) satisfying initial condition u(0) = u 0 if and only if A is the generator of a strongly continuous semi-groups S(t) of bounded linear operators on X, and in this case u(t) is given by (1.53).

Theorem 9. [Stone's theorem [START_REF] Dautray | Analyse mathématique et calcul numérique pour les sciences et les techniques[END_REF]] Let A be a closed operator of domain D(A) dense in X, such that A = iH with H self-adjoint, then A is the infinitesimal generator of a unitary group in X.

In the following Lemma, if u belongs to H 2 and 2u is square integrable, then we define the traces of u on {t = 0} ∪ {t = T } and similarly traces u on {x = 0} ∪ {x = L}.

Lemma 10. Let u belongs to H 2 , and assume that 2u = ü -u (defined in the sense of distributions on Q T ) is square integrable. Then, for ε > 0, we have

u ∈ C 0 ( [0, T ]; H 1 (ε, L) ) and u ∈ C 0 ( [0, T ]; L 2 (ε, L) ) , (1.54a) u ∈ C 0 ( [0, L]; H 1 (ε, T -ε) ) and u ∈ C 0 ( [0, L]; L 2 (ε, T -ε) ) . (1.54b) Proof. Let ψ belong to C ∞ ([0, L]) be such that ψ = { 1 on [ε, L], 0 on [0, ε/2].
Let us define ũ = uψ, then we may deduce

2ũ = f ∈ L 2 (Q T ), ũ(0, t) = ũ(L, t) = 0.
Chapter 1. Convergence of mass redistribution method for the one-dimensional wave equation with a unilateral constraint at the boundary If we let

X = H 1 0 (0, L) × L 2 (0, L) and D(A) = [H 1 0 (0, L) ∩ H 2 (0, L)] × H 1 0 (0, L),
and with

A [ u v ] = - [ v u ] for all û = [ u v ] ∈ D(A),
then A is skew-adjoint and therefore, with Stone's theorem, A generates a unitary group S(t) on X [START_REF] Schatzman | Numerical approximation of a wave equation with unilateral constraints[END_REF].

If we let û = [ ũ u] and F = [ 0 f] ,
then by Theorem 8, there is a unique weak solution û of dû dt + Aû = F , and according to Definition 9, and using the group property of S, û belongs to C 0 ([0, T ]; X) which proves (1.54a). To establish

(1.54b), first let χ ∈ C ∞ ([0, L]) be such that χ = { 1 on [ε, T -ε], 0 on [0, ε/2] ∪ [T -ε/2, T ].
Let us define ũ = uχ.

Then, we can argue exactly as above.

Chapter 2

Numerical approximations of a one dimensional elastodynamic contact problem based on mass redistribution method

This chapter is submitted for publication [START_REF] Dabaghi | Numerical approximations of elastodynamic contact problems based on mass redistribution method[END_REF] Contents

Introduction

This paper aims to give some new numerical results for the dynamical evolution with impact of a linearly elastic body and a rigid obstacle by using the mass redistribution method introduced in [START_REF] Khenous | Mass redistribution method for finite element contact problems in elastodynamics[END_REF]. The situations involving contact abound in industrial or biomedical problems. For this reason a considerable engineering and mathematical literature has been devoted to contact problems. Usually, with the hypotheses of small deformations, the contact phenomena is modelled by using the so-called Signorini's boundary conditions in displacement, which are based on a linearization of the physically meaningful non penetrability of masses. The elastodynamic contact problems are typically stated as hyperbolic type variational inequalities in Sobolev spaces, and a few existence of solutions results has been established. The reader is referred to [START_REF] Schatzman | A hyperbolic problem of second order with unilateral constraints: The vibrating string with a concave obstacle[END_REF][START_REF] Lebeau | A wave problem in a half-space with a unilateral constraint at the boundary[END_REF][START_REF] Kim | A boundary thin obstacle problem for a wave equation[END_REF] and also the comprehensive monograph [START_REF] Eck | Unilateral contact problems : variational methods and existence theorems[END_REF]. The uniqueness result was successfully investigated for a wave equation in a half-space with a unilateral contact at the boundary in [START_REF] Lebeau | A wave problem in a half-space with a unilateral constraint at the boundary[END_REF]. In particular, the uniqueness still remains an open problem in other frameworks. The fundamental mathematical difficulties are related to the intrinsic non-smoothness of the problem coming from the Signorini boundary conditions. Finally, notice that the existence of a strong solution for viscoelastodynamic problems with unilateral constraints is proved in [START_REF] Petrov | Mathematical results on existence for viscoelastodynamic problems with unilateral constraints[END_REF].

On the other hand, another challenging task consists to develop efficient numerical methods for solving dynamical contact problems. Newmark's method is one of the most popular numerical solvers which is usually combined with finite elements methods. Unfortunately, for almost all values of the parameters of this method, unphysical energy blow-ups occur during the time integration as well as numerical instabilities at the dynamical contact boundaries. To overcome these difficulties, some numerical methods based on the Newmark scheme for solving impact problems and which are energy dissipative have been introduced in [START_REF] Carpenter | Lagrange constraints for transient finite element surface contact[END_REF]. Although stable, those methods lead to an important energy loss at impact which do not vanish when the time step decreases. Laursen et al. in [START_REF] Laursen | Improved implicit integrators for transient impact problemsgeometric admissibility within the conserving framework[END_REF][START_REF] Laursen | Design of energy conserving algorithms for frictionless dynamic contact problems[END_REF][START_REF] Chawla | Energy consistent algorithms for frictional contact problems[END_REF][START_REF] Laursen | Computational contact and impact mechanics. Fundamentals of modeling interfacial phenomena in nonlinear finite element analysis[END_REF] have designed time integration schemes of Newmark type which are energy conserving. However, these schemes are unable to prevent some spurious oscillations of the displacement and of the contact stress on the contact boundary. These unphysical oscillations are avoided by the numerical methods developed in [START_REF] Deuflhard | A contact-stabilized Newmark method for dynamical contact problems[END_REF] (see also the reference therein), but these methods are still energy dissipative. Another approach consisting to remove the mass at the contact nodes is investigated in [START_REF] Khenous | Mass redistribution method for finite element contact problems in elastodynamics[END_REF]. This mass redistribution method prevents the oscillations at the contact boundaries mentioned above, and it can be proven that the solution is unique and energy conserving, see [START_REF] Dabaghi | Convergence of mass redistribution method for the one-dimensional wave equation with a unilateral constraint at the boundary[END_REF]. Furthermore, some different variants has been proposed in [START_REF] Hager | A stable energy-conserving approach for frictional contact problems based on quadrature formulas[END_REF] with an estimate of the space approximation for linear elasticity. It has also been extended for thin structures in [START_REF] Renard | The singular dynamic method for constrained second order hyperbolic equations: application to dynamic contact problems[END_REF][START_REF] Pozzolini | Some energy conservative schemes for vibro-impacts of a beam on rigid obstacles[END_REF][START_REF] Tkachuk | Hybrid-mixed discretization of elastodynamic contact problems using consistent singular mass matrices[END_REF] considering different discretizations for the displacement and the velocity.

The paper is organised as follows. In Section 2.2, the mathematical models are presented and a new proof of existence and uniqueness results is proposed. A space semi-discretization based on a redistribution of mass is introduced in Section 2.3 and then some time integration method are introduced. In Sections 2.4 and 2.5, two benchmark problems with their analytical solutions are exhibited, one being new. More precisely, Section 2.4 is concerned by an elastodynamic problem with Signorini and Neumann boundary while Section 2.5 deals with an elastodynamic problem with Signorini and Dirichlet boundary conditions. Finally, numerical experiments for 2.2. Description of the models these benchmark problems are reported and analysed. In particular, convergence rates for the displacement, Lagrange multiplier and energy evolution show the efficiency of mass redistribution method.

Description of the models

We consider a model for an elastic bar of length L vibrating longitudinally. One end of this bar is free to move, as long as it does not hit a material obstacle, while the other end can be clamped or free to move. The obstacle constrains the displacement of the extremity to be non negative. The material of the bar is supposed to be homogeneous and the theory of small deformations is considered. Let us denote by x the spatial coordinate along the bar, with the origin at the material obstacle, let u(x, t) be the displacement at time t ∈ [0, T ], T > 0 of the material point of spatial coordinate x ∈ [0, L]. Let f (x, t) denotes a density of external forces, depending on time and space. In the case where the bar is clamped at one end (see Figure 2.5), the mathematical problem is formulated as follows:

(2.1)

u tt (x, t) -u xx (x, t) = f (x, t), (x, t) ∈ (0, L) × (0, T ),
with Cauchy initial data

(2.2) u(x, 0) = u 0 (x) and u t (x, 0) = v 0 (x), x ∈ (0, L),
and Signorini and Dirichlet boundary conditions at x = 0 and x = L, respectively,

(2.3) 0 ≤ u(0, t) ⊥ u x (0, t) ≤ 0 and u(L, t) = 0, t ∈ [0, T ].
While in the case where the bar is free to move at one end (see Figure 2.1), the motion is governed by the momentum equilibrium equation (2.1) together with Cauchy initial data (2.2) and Signorini and Neumann boundary conditions at x = 0 and x = L, respectively,

(2.4) 0 ≤ u(0, t) ⊥ u x (0, t) ≤ 0 and u x (L, t) = 0, t ∈ [0, T ].
Here u t def = ∂u ∂t and u x def = ∂u ∂x . The orthogonality has the natural meaning; namely if we have enough regularity, it means that at the boundary x = 0 the product u(0, •)u x (0, •) vanishes almost everywhere in time. If it is not the case, the above inequality is integrated on an appropriate set of test functions, leading to a weak formulation for the unilateral condition. Observe that from mathematical viewpoint, Signorini's condition means that when the bar touches the obstacle at x = 0, its reaction can be only upwards, so that u x (0, •) ≤ 0 on the set {t : u(0, •) = 0}. While in the case where the bar does not touch the obstacle, its end is free to move. More precisely, we have u x (0, •) = 0 on the set {t : u(0, •) > 0}. We suppose that the initial displacement u 0 belongs to the H 1 (0, L) and satisfies the compatibility conditions, i.e. u 0 (L) = 0 and u 0 (0) ≥ 0, the initial velocity v 0 belongs to L 2 (0, L) and the density of external forces f belongs to L 2 (0, T ; L 2 (Ω)).

We describe now the weak formulations of both the two problems. For that purpose, it is convenient to introduce the following notations:

V Dir def = {u ∈ H 1 (0, L) : u(L) = 0}, V Neu def = H 1 (0, L) and H def = L 2 (0, L).
We denote by K the following convex set:

K def = {u ∈ H 2 : u(•, t) ∈ K for almost every t},
         find u ∈ K such that - ∫ L 0 v 0 (v(•, 0)-u 0 ) dx- ∫ QT u t (v t -u t ) dx dt+ ∫ QT u x (v x -u x ) dx dt≥ ∫ QT f (v-u) dx dt for all v ∈ K for which there exists ζ > 0 with v = u for t ≥ T -ζ.
Notice that the sole difference of the weak formulation associated to (2.1)-(2.3) and to (2.1)-(2.2), (2.4) comes from the definition of the convex set K.

It should also be noted that existence and uniqueness results are obtained for a similar situation of a vibrating string with concave obstacle in one dimensional space in [START_REF] Schatzman | A hyperbolic problem of second order with unilateral constraints: The vibrating string with a concave obstacle[END_REF] and also for a wave equations with unilateral constraint at the boundary in a half space of R N in [START_REF] Lebeau | A wave problem in a half-space with a unilateral constraint at the boundary[END_REF]. An existence result for a wave equation in a C 2 -regular bounded domain constrained by an obstacle at the boundary in R N for N ≥ 2 is proven in [START_REF] Kim | A boundary thin obstacle problem for a wave equation[END_REF].

Theorem 10 (Existence and uniqueness results). Assume that u

0 ∈ V, v 0 ∈ H and f ∈ L 2 (0, T ; H). Then Problem (2.1)-(2.3) admits a unique solution u ∈ L ∞ (0, T ; V)∩W 1,∞ (0, T ; H).
Proof. The first step consists to introduce an auxiliary problem. Namely, let ū be the solution to Problem (2.1) with initial data (2.2) and boundary conditions (2.3) where Signorini's boundary condition is replaced by Dirichlet's one, i.e. ū(0, t) = 0, t ∈ [0, T ]. Note that the auxiliary problem possesses a unique solution ū belonging to C 0 ([0, T ]; H 2 (0, [START_REF] Lions | Problèmes aux limites non homogènes et applications[END_REF][START_REF] Brezis | Functional analysis, Sobolev spaces and partial differential equations[END_REF]). Then we denote by

L) ∩ H 1 0 (0, L)) ∩ C 1 ([0, T ]; H 1 0 (0, L)) ∩ C 2 ([0, T ]; H) (see
v def = u -ū a solution of (2.6) v tt (x, t) -v xx (x, t) = 0, (x, t) ∈ (0, L) × (0, T ),
with Cauchy initial data (2.7) v(x, 0) = 0 and v t (x, 0) = 0, x ∈ (0, L), and Signorini and Dirichlet boundary conditions at x = 0 and x = L, respectively,

(2.8) 0 ≤ v(0, t) ⊥ v x (0, t) + ūx (0, t) ≤ 0 and v(L, t) = 0, t ∈ [0, T ].

Define g(t)

def = ūx (0, t). The second step consists to rewrite (2.6)-(2.8) as a differential inclusion problem by using the characteristic method. To this aim, it is convenient to introduce the characteristic coordinates

ξ def = x + t and η def = x -t.
Therefore the chain rule gives

∂ 2 v ∂x 2 = ∂ 2 v ∂ξ 2 + 2 ∂ 2 v ∂ξ∂η + ∂ 2 v ∂η 2 and ∂ 2 v ∂t 2 = ∂ 2 v ∂ξ 2 -2 ∂ 2 v ∂ξ∂η + ∂ 2 v ∂η 2 .
According to (2.1), ∂ 2 v ∂ξ∂η vanishes and it follows that v(ξ, η) = φ(ξ) + ψ(η), where φ and ψ are two differentiable functions such that (2.9) v(x, t) = φ(x+t) + ψ(x-t).

Description of the models

In particular, taking t = 0 and using the Cauchy initial data, we get (2.10) φ(x) + ψ(x) = 0 and φ (x) -ψ (x) = 0, which by integration implies that there exist two constants C φ and C ψ such that φ(ξ) = C φ and ψ(η) = C ψ for all ξ and η belonging to [0, L]. According to (2.9), the boundary conditions (2.3) can be rewritten as follows:

( We introduce the multivalued function J N : R → P( R) defined by

J N [x] def =        {0} if x < 0, [0, +∞) if x = 0, ∅ if x > 0,
where P( R) is the set of all subsets of R. More precisely, J N (x) is the subdifferential of the indicator function defined by

I (-∞,0] (x) def = { 0 if x ∈ (-∞, 0], +∞ if x / ∈ (-∞, 0].
Obviously, I (-∞,0] is a lower semi-continuous and convex function, for further details the reader is referred to [START_REF] Brézis | Opérateurs maximaux monotones et semi-groupes de contractions dans les espaces de Hilbert[END_REF]. Then, the inequalities in (2.11) can be rewritten as follows (2.12)

φ (t) + ψ (-t) + g(t) ∈ -J N [-φ(t)-ψ(-t)].
Note that at this stage, ψ(η), η ∈ [-L, L], is the unique unknown of (2.12). Let us define now

(2.13) w(t) def = -φ(t) -ψ(-t).
We insert (2.13) into (2.12) to get

w (t) ∈ -J N [w(t)] -g(t).
Finally, we find the following Cauchy problem

w (t) ∈ -J N [w(t)] -g(t) a.e. t ∈ (0, L), (2.14a) w(0) = 0. (2.14b) Chapter 2.

Numerical approximations of a one dimensional elastodynamic contact problem based on mass redistribution method

Then we look for a solution of the following problem:

w (t) + h(t) = -g(t) a.e. t ∈ (0, L), (2.15a) h(t) ∈ J N [w(t)] a.e. t ∈ (0, L), (2.15b) w(0) = 0. (2.15c)
We observe that (2.15) is equivalent to (2.11), the verification is left to the reader. We regularize the problem (2.15), we construct a solution of the regularized problem, and we show the existence of a solution by passing to the limit with respect to regularity parameter.

We regularize the problem (2.15) by using Yosida's regularization. More precisely, we introduce first a sequence {g ε } ε>0 such that g ε belongs to C 0 ([0, L]) for all ε > 0 and (2.16)

g ε → g in H as ε 0.
We may approximate (2.15) by the following problem:

w ε (t) + J N,ε [w ε (t)] = -g ε (t) a.e. t ∈ (0, L), (2.17a) w ε (0) = 0. (2.17b) where J N,ε def = 1 ε (1-J ε )
denotes the Yosida approximation with the resolvant defined by J ε def = (1+εJ N ) -1 . Note that J N,ε is a monotone and Lipschitz continuous mapping defined on all H. Clearly the contraction principle step by step in time allows us to prove that (2.15) possesses a unique solution w ε belonging to C 1 ([0, L]).

We establish now some a priori estimates, which later will enable to infer the existence of a weak solution. To this aim, we multiply (2.17a) by w ε (t), and we integrate this expression over (0, L) to get (2.18)

∫ L 0 |w ε (t)| 2 dt + ∫ L 0 J N,ε [w ε (t)]w ε (t) dt = - ∫ L 0 g ε (t)w ε (t) dt.

Let us define now

I ε (-∞,0] (w) def = min z∈H { 1 2ε |w-z| 2 + I (-∞,0] (z)
} for all ε > 0 and w ∈ H. Since J N is the subdifferential of a proper, convex and lower semicontinuous function I (-∞,0] , it follows that I ε (-∞,0] is convex, Fréchet differentiable in H and its subdifferential ∂I ε (-∞,0] coincides with J N,ε . Furthermore we have

∀z ∈ H : I ε (-∞,0] (z) I (-∞,0] (z) as ε 0.
Then we may deduce from (2.18) that

(2.19) 1 2 w ε 2 H + I ε (-∞,0] (w ε (L)) ≤ 1 2 g ε 2 H .
We observe that by definition I ε (-∞,0] (w ε (L)) ≥ 0 which follows from (2.19) that there exists

C 1 > 0, independent of ε > 0, such that (2.20) w ε H ≤ C 1 .

Description of the models

Furthermore we have

(2.21) |w ε (τ )| 2 ≤ 2L w ε 2 H
for all τ ∈ (0, L). Then we deduce from (2.19) and (2.21) that there exists C 2 > 0, independent of ε > 0, such that

(2.22) |w ε (τ )| 2 + w ε 2 H ≤ C 2
for all τ ∈ (0, L). Since J ε is a contraction on all H (see [START_REF] Brézis | Opérateurs maximaux monotones et semi-groupes de contractions dans les espaces de Hilbert[END_REF]), it follows from (2.22) that there exists C 3 > 0, independent of ε > 0, such that (2.23)

w ε V + J ε w ε V ≤ C 3 .
We observe that (2.17a) implies that (2.24)

J N,ε [w ε ] 2 H ≤ 2 w ε 2 H + 2 g ε 2 H .
The above a priori estimates allow us to infer that there exists w and h belonging to V∩L ∞ (0, L), respectively, and passing to subsequences, if necessary, we find

J ε w ε w in L ∞ (0, L) weak * , (2.25a) J N,ε [w ε ]
h in H weak, (2.25b) as ε tends to 0. Since V → H with compact embedding, it follows that

J ε w ε → w in H, (2.26a) w ε → w in H and w ε w in V weak, (2.26b)
as ε tends to 0. It remains to prove that w and h satisfy (2.15). We integrate (2.17a) over (0, τ ), we get (2.27)

∫ τ 0 ( w ε (t)+J N,ε [w ε (t)]+g ε (t) ) dt = 0
for all τ ∈ (0, L). Thanks to (2.25) and (2.26), we may pass to the limit in all the terms of (2.27), namely we find

∫ τ 0 ( w (t)+h(t)+g(t) ) dt = 0 for all τ ∈ (0, L), which is equivalent to (2.15a). Let us note that J N,ε [w ε ] ∈ J N [J ε w ε ],
for some technical details, the reader is referred to [START_REF] Brézis | Opérateurs maximaux monotones et semi-groupes de contractions dans les espaces de Hilbert[END_REF]. Then (2.25b) and (2.26a) enable us to deduce that (2.28) lim sup

ε 0 ∫ L 0 J N,ε [w ε (t)]J ε w ε (t) dt ≤ ∫ L 0 h(t)w(t) dt,
which implies that (2.15b) holds. Then we may conclude that there exists a solution v to Problem (2.6)-(2.8).

Chapter 2. Numerical approximations of a one dimensional elastodynamic contact problem based on mass redistribution method

The uniqueness result comes from the monotonicity of J N , namely J N is the subdifferential of a convex, lower semi-continuous and proper function (see [START_REF] Brézis | Opérateurs maximaux monotones et semi-groupes de contractions dans les espaces de Hilbert[END_REF] for further details).

The existence and uniqueness results to Problem (2.1)-(2.2) together with (2.4) is established by using analogous approach developed above for Problem (2.1)-(2.3). The main difference comes from the choice of the auxiliary problem. Let ū be the unique solution to Problem (2.1) with initial data (2.2) and boundary conditions (2.4) where the Signorini boundary conditions are replaced by Dirichlet boundary conditions, i.e. ū(0, t) = 0, t ∈ [0, T ]. Hence the rest of the proof is quite a routine, the verification is left to the reader.

Finally, we consider the following energy which is associated to (2.5):

E(t) = 1 2 ∫ L 0 ( |u t (x, t)| 2 +|u x (x, t)| 2 ) dx.
This energy is constant with respect to time t when the density forces f vanishes.

Note that formulations (2.5) and (2.15) are equivalent (in the sense that they have the same weak solutions), since the proof given in [START_REF] Dabaghi | Convergence of mass redistribution method for the one-dimensional wave equation with a unilateral constraint at the boundary[END_REF] can be straightforwardly adapted.

Discretization

Classical finite element discretization

We introduce now a semi-discrete problem in space associated to (2.5) . To this aim, we choose a set of parameters h def = L m (mesh size) having in mind the limit m → +∞, where m is an integer, and let V Neu

h def = {v h ∈ C 0 ([0, L]) : v h | [a i ,a i+1 ] ∈ P 1 , i = 0, . . . , m -1} and V Dir h def = {v h ∈ V Neu h : v h (L) = 0}.
Here, a i def = ih, i = 0, 1, . . . , m and P 1 is the space of polynomials of degree less than or equal to 1. A classical basis of V h is given by the sequence of shape functions ϕ i ∈ V h for i = 0, 1, . . . , J, defined by

ϕ i (x) def = { 1 -|x-a i | h if x ∈ [a max(i-1,0) , a min(i+1,m) ], 0 otherwise,
where J = m -1 in the case where V h = V Dir h and J = m in the case where V h = V Neu h . Observe that ϕ i (a j ) = δ ij , j = 0, 1, . . . , m (δ is Kronecker's symbol). We approximate the solution u belonging to V to the weak formulation (2.5) by

u h (x, t) = J ∑ j=0 u j (t)ϕ j (x).
Consequently, we have u i = u h (a i ), i = 0, 1, . . . , J. We shall consider two strategies to approximate the dynamic contact problem. The first one is a classical finite element semi-discretization

Discretization

which uses a multiplier. It reads as follows:

(2.29)

               find u h : [0, T ] → V h and λ : [0, T ] → R such that for all v h ∈ V h ∫ L 0 u h,tt v h dx + ∫ L 0 u h,x v h,x dx = -λv h (0) + ∫ L 0 f v h dx a.e. t ∈ [0, T ], 0 ≤ u h (0, •) ⊥ λ ≤ 0 a.e. t ∈ [0, T ], u h (•, 0) = u 0 h and u h,t (•, 0) = v 0 h ,
where u 0 h and v 0 h belong to V h and λ is the Lagrange multiplier being here the contact force. Furthermore, the approximated problem (2.29) can be rewritten in an algebraic formulation form as follows:

(2.30)

           find U h : [0, T ] → R n and λ : [0, T ] → R such that M U h,tt + SU h = -λe 0 + F a.e. t ∈ [0, T ], 0 ≤ u 0 ⊥ λ ≤ 0 a.e. t ∈ [0, T ], U h (0) = U 0 h and U h,t (0) = V 0 h ,
where M and S denote the mass and stiffness matrices, respectively, and

U def = (u 0 , . . . , u n-1 ) T , e 0 def = (1, 0, . . . , 0) T and F def = ( ∫ L 0 f ϕ 0 dx, . . . , ∫ L 0 f ϕ n-1 dx).

Finite element discretization with the mass redistribution method

An alternative to the standard discretization is to consider a trivial mass redistribution method which consists to replace the mass matrix M in (2.30) by a modified mass matrix M mod defined by

M mod ij def = ∫ L h ϕ i ϕ j dx for all i, j ∈ [0, n -1].
Let us observe that there are more sophisticated mass redistribution methods, see [START_REF] Hager | Analysis of a space-time discretization for dynamic elasticity problems based on mass-free surface elements[END_REF]. The modified mass matrix has the form

M mod def = ( 0 0 0 M ) ,
where M admits an unmodified part Mij = M i+1,j+1 , for all i, j = 1, . . . , n-2. As a consequence and in a certain sense, the node at the contact boundary evolves in a quasi-static way. It has been established in [START_REF] Dabaghi | Convergence of mass redistribution method for the one-dimensional wave equation with a unilateral constraint at the boundary[END_REF] in a slightly different context that the approximated solution using the mass redistribution method converges to the unique solution to (2.5).

Finally, note that the discrete energy associated to problem (2.30) is given by

(2.31) E h (t) = 1 2 (U T h,t M U h,t +U T h SU h -U T h,t F )(t).
We introduce now the time discretization. In order to fix the notations, let the time interval [0, T ] be divided by n + 1 discrete time-points such that 0 = t 0 < t 1 < . . . < t n = T . Furthermore the discrete quantities U n h , U n h,t , U n h,tt and λ n are assumed to be given by algorithmic approximations Chapter 2. Numerical approximations of a one dimensional elastodynamic contact problem based on mass redistribution method of the displacement U h (t n ), the velocity U h,t (t n ) and the acceleration U h,tt (t n ) and the Lagrange multiplier λ(t n ), respectively. Some time-stepping schemes allowing to obtain an approximated solution to Problem (2.1)-(2.4) are introduced below and their efficiency is discussed and analyzed in the next sections.

The most wide-spread time-stepping scheme for solving contact problems is the family of classical Newmark methods proposed by Newmark in [START_REF] Newmark | A method of computational for structural dynamics[END_REF]. The algorithms form a subset of the Hilber-Hughes-Taylor (HHT) family of temporal integrators, the reader is referred to [START_REF] Hilber | Improved numerical dissipation for time integration algorithms in structural dynamics[END_REF] for further details. The underlying concept of the discretizations are Taylor expansions of displacements and velocities neglecting terms of higher order. For the contact problem introduced above, the discrete evolution is described by the finite difference equations:

U n+1 h = U n h + ∆tU n h,t + ( 1 2 -β ) ∆t 2 U n h,tt + β∆t 2 U n+1 h,tt , (2.32a) U n+1 h,t = U n h,t + (1-γ)∆tU n h,tt + γ∆tU n+1 h,tt , (2.32b) M U n+1 h,tt + SU n+1 h = -λ n+1 e 0 + F n+1 , (2.32c) 0 ≤ u n+1 0 ⊥ λ n+1 ≤ 0, (2.32d)
where ∆t is a given time step and (β, γ) are the algorithmic parameters. Furthermore U 0 h , U 0 h,t and λ 0 are given and U 0 h,tt is evaluated by using (2.32c). Notice that the Newmark family contains many well-known and widely-used algorithms which correspond to different choices of the algorithmic parameters (β, γ) (see [START_REF] Hughes | Linear static and dynamic finite element analysis[END_REF][START_REF] Laursen | Computational contact and impact mechanics. Fundamentals of modeling interfacial phenomena in nonlinear finite element analysis[END_REF]). This method is unconditionally stable for linear elastodynamic problem for γ ≥ 1 2 and β ≥ 1 4 ( 1 2 +γ) 2 (see [START_REF] Hughes | Linear static and dynamic finite element analysis[END_REF][START_REF] Krenk | Energy conservation in Newmark based time integration algorithms[END_REF]). One of the most used method to study dynamic problems is the Crank-Nicolson method, also called trapezoidal or average acceleration method, which is obtained by setting the parameters (β, γ) = ( 14 , 1 2 ). This method is second-order consistent and unconditionally stable in the unconstrained case. Furthermore, the total energy of the discrete evolution is preserved, for the purely elastic case, see [START_REF] Laursen | Computational contact and impact mechanics. Fundamentals of modeling interfacial phenomena in nonlinear finite element analysis[END_REF]. That is the reason why the Crank-Nicolson method is commonly used in the community of computational mechanics. However, the situation should be examined in the case of contact constraints. Indeed the order of accuracy is degraded, for further details, the reader is referred to [START_REF] Hughes | Linear static and dynamic finite element analysis[END_REF][START_REF] Krenk | Energy conservation in Newmark based time integration algorithms[END_REF][START_REF] Grosu | Stability of semidiscrete formulations for elastodynamics at small time steps[END_REF]. Let us define now the energy evolution by

∆E n h def = E n+1 h -E n h
, where E n h is assumed to be given by an algorithmic approximation of the energy E h (t n ). In particular, the energy evolution associated to (2.32) for (β, γ) = ( 14 , 1 2 ) is given by (2.33)

∆E n h = -(V n h ) T Λ n γ e 0 ,
where

V n h def = U n+1 h -U n h and Λ n γ def = (1-γ)λ n + γλ n+1 .
Let us observe that in the presence of permanent contact, the energy is strictly conserved, while the release of an existing contact decreases the energy and the detection of a new contact increases the energy. However, regain and loss of energy do not balance and there is no guarantee that the energy will stay bounded during the time integration.

Another approach to study elastodynamic contact problems consists to use the θ-method.

The wave equation with Signorini and Neumann boundary conditions

More precisely, the discrete evolution is described by the finite difference equations:

U n+1 h = U n h + ∆t((1-θ)U n h,t +θU n+1 h,t ), (2.34a) U n+1 h,t = U n h,t + ∆t((1-θ)U n h,tt +θU n+1 h,tt ), (2.34b) M U n+1 h,tt + SU n+1 h = -λ n+1 e 0 + F n+1 , (2.34c) 0 ≤ u n+1 0 ⊥ λ n+1 ≤ 0, (2.34d) with θ ∈ [0, 1]. Furthermore U 0
h , U 0 h,t and λ 0 are given and U 0 h,tt is evaluated by using (2.34c). The energy evolution associated to (2.34) is given by

(2.35) ∆E n h = ( 1 2 -θ ) (V n h,t ) T M V n h,t + ( 1 2 -θ ) (V n h ) T SV n h -(V n h ) T Λ n θ e 0 .
Notice that (2.34) with θ = 0 and θ = 1 are called the forward and backward Euler method, respectively, while with θ = 1 2 , we obtain Crank-Nicolson's method. Finally, we focus on the so-called Paoli-Schatzman method that consists to fix the contact constraint at an intermediate time step. Indeed the method proposed below is a slight modification of Paoli-Schatzman method (see [START_REF] Paoli | Time discretization of vibro-impact[END_REF][START_REF] Paoli | A numerical scheme for impact problems. I ,II. The multidimensional case[END_REF]) which takes into account the kernel of the modified mass matrix. A simple application of Paoli-Schatzman method based on Newmark scheme to Problem (2.30) with γ = 1 2 leads to (2.36)

                 find U n+1 h : [0, T ] → R n and λ n : [0, T ] → R such that: M (U n+1 h -2U n h +U n-1 h ) ∆t 2 + S(βU n+1 h +(1-2β)U n h +βU n-1 h ) = -λ n e 0 for all n ≥ 2, 0 ≤ u n,e 0 = u n+1 0 + eu n-1 0 1 + e ⊥ λ n ≤ 0, U 0 and U 1 given.
Here e belongs to [0, 1] and is aimed to be interpreted as a restitution coefficient. Note that U 0 h and U 1 h are given data and U 1 h can be evaluated by a one step scheme. We may observe that taking M = M mod in (2.36), we are not able to resolve the problem on the kernel of M mod . That is the reason why, SU n-1 h as well as SU n h are projected on the orthogonal of the kernel of M . The energy evolution associated to (2.36) is given for β = 1 4 by

(2.37)

∆E n h = 1+e 2 λ n u n-1 0 .
Notice that the numerical simulations were performed by employing the finite element library Getfem++ (see [START_REF] Renard | Getfem++. An Open Source generic C++ library for finite element methods[END_REF]). For the reader convenience, the energy evolutions (2.33), (2.35) and (2.37) are justified in the Appendix.

The wave equation with Signorini and Neumann boundary conditions

Analytical solution

We consider in this section the analytical solution proposed in [START_REF] Wriggers | Computational contact mechanics[END_REF]. It is the solution to the equation of motion (2.1) together with Cauchy initial data (2.2) and boundary conditions (2.4) in the case where the density of external forces f (x, t) = -(g+u x (0, t)) with g > 0 being the gravity acceleration. Let us describe the solution of this problem (see Figure 2.1). Before the impact, the bar is undeformed and its initial velocity v 0 vanishes. The bar reaches the rigid obstacle at

time t 1 = √ 2u 0
g with the velocity √ 2u 0 g. After the impact, the bar stays in contact as long as the shock wave travels from bottom to the top of the bar and vice versa and then it takes off. Since the velocity of the shock wave c is assumed to be equal to 1, the bar stays in contact with the rigid obstacle a time L. Finally, notice that impacts occur at t 4k+1 = 3L + 16kL, t 4k+2 = t 4k+1 + 2L, t 4k+3 = t 4k+1 + 8L and t 4k+4 = t 4k+1 + 10L. We may deduce that the

(2.38) u(x, t) =                  u 0 -1 2 gt 2 if t ≤ t 1 , h 1 (x, t-t 4k+1 ) if t 4k+1 < t ≤ t 4k+2 , h 2 (x, t-t 4k+2 ) if t 4k+2 < t ≤ t 4k+3 , h 1 (x, t 4k+4 -t) if t 4k+3 < t ≤ t 4k+4 , u 0 -1 2 g ( t-t 4k+4 - √ 2u 0 g ) 2 if t 4k+4 < t ≤ t 4(k+1)+1 .
where we have used the following notations

h 1 (x, t) def = - √ 2u 0 g min(x, L-|t-L|) + +∞ ∑ n=1 2g Lν 3 n (cos(ν n t)-1) sin(ν n x), (2.39) h 2 (x, t) def = u 0 - 1 2 g ( t- √ 2u 0 g ) 2 - 2gL 2 3 + +∞ ∑ n=1 4g λ 2 n cos(λ n t) cos(λ n x), (2.40) with ν n def = (n-1 2 ) π L and λ n def = n π L .

Numerical simulations

The parameters used in the numerical simulations are L = 1, g = 0.11, u 0 (x) = 0.5, v 0 (x) = 0 and u x (L, •) = 0. We evaluate the approximated solution (U n h , λ n ) to Problem (2.1)-(2.3) by using Crank-Nicolson method that consists to choose (β, γ) = ( 14 , 1 2 ) in (2.32) for the standard and modified mass matrices. Then we compare the solutions obtained by employing Crank-Nicolson method to the analytical one (u, λ) in different points of the bar (see Figure 2.2). In the case where the constraint is active, the numerical experiments highlighted some spurious oscillations occurring for the approximated solution (U n h , λ n ) obtained by solving (2.32) with the standard mass matrix M (see Figure 2.2 (left)) while these oscillations disappeared when the standard mass matrix is replaced by the modified mass matrix, i.e. M = M mod (see Figure 2.2 (right)). However the approximated solutions obtained by both methods converge to the analytical solution when both the space ∆x and time ∆t steps tend to 0 in L p (0, T ; H) norms with p = 2, +∞, see Figure 2.3 (a constant ratio ∆x ∆t = 1 has been used).

The wave equation with Signorini and Neumann boundary conditions

Surprisingly, the error curves in the norms u h -u L p (0,T ;V) , with p = 2, +∞ shown in the same figure are both diverging for the standard and modified mass matrices. The explanation of this

(a) u h -u L ∞ (0,T ;H) (b) u h -u L 2 (0,T ;H) (c) u h -u L ∞ (0,T ;V) (d) u h -u L 2 (0,T ;V)

The wave equation with Signorini and Dirichlet boundary conditions

divergence is illustrated in Figure 2.4 where the space derivative of the displacement is shown for both the analytical and the approximated solution. Indeed, Figure 2.4 (left) illustrates the fact that the series of the derivative corresponding to (2.40) do not converge. This is why it seems not to be possible to use this analytical solution to estimate the error in L p (0, T ; V) norm. For that reason, another analytical solution is exhibited in the next section.

The wave equation with Signorini and Dirichlet boundary conditions

Analytical solution

In this section, we propose a new analytical piecewise affine and periodic solution to the problem composed by the momentum equilibrium equation (2.1) together with Cauchy initial data (2.2) and boundary conditions (2.3) in the case where the density of external forces f (x, t) vanishes. Suppose that the initial displacement and velocity are given by u 0 (x) = 1 2 -x 2 and v 0 (x) = 0, respectively. The length of bar is L = 1, namely, the bar is compressed at t = 0. Let us describe the solution of this problem (see Figure 2.5). Before the impact, the bar, which is clamped at one end, elongates until it reaches the rigid obstacle at t 1 = 1. After the impact, the bar stays in contact for the time t ∈ (t 1 , t 2 ) and it takes off at time t 2 = 2. We use characteristics method to determine the solution of the above problem. To this aim, we split each domain defined by (0, L) × (t i , t i+1 ), i = 0, 1, 2, corresponding to the phases before the impact, during the impact and after the impact, respectively, into four regions according to Figure 2.6. We choose below L = 1 and t i = i, i = 0, . . . , 3. Once again by using the characteristics lines, we divide the domain (1, 2)×(0, 1) into four regions. Then the solution (2.41b) calculated in the region IV allows us to evaluate the initial conditions at time t 1 = 1, namely we found u(•, 1) = 0 and u t (•, 1) = -1 2 . Therefore we obtain

x = 0 L t = 0 t 2 = 2 t 3 = 3 t 1 = 1
v = 0 u = 1-x 2 t 0 = 0 v = -1/2 u = 0 t 1 = 1 v = 1/2 u = 0 t 2 = 2 v = 0 u = 1-x 2 t 3 = 3
u(x, t) = -1 2 (t 1 -t) in the region I, (2.42a) u(x, t) = -x 2 in the region II, (2.42b) u(x, t) = 1 2 (x-1) in the region III, (2.42c) u(x, t) = 1 2 (t-t 1 -1) in the region IV. (2.42d)
After the impact, we divide the domain (2, 3) × (0, 1) into four regions. It follows from (2.42d) that u(•, 2) = 0 and u t (•, 2) = 1 2 which implies that u(x, t) = 1 2 (t-t 2 ) in the regions I, II, u(x, t) = 1 2 (1-x) in the regions III, IV.

Furthermore we observe that u(•, 3) = u 0 and u t (•, 3) = v 0 . Then we conclude that the solution is periodic of period 3. Note that the Lagrange multiplier λ is equal to the contact force u x (0, •).

Numerical simulations

We evaluate below (U n h , λ n ) by using different time-stepping methods, namely the Newmark method (2.32), the backward Euler method (2.34) with θ = 1 and the Paoli-Schatzman method (2.36). More precisely, for each time-stepping method, the approximated solutions (U n h , λ n ) are obtained by solving (2.32) for different values of the couple (β, γ), (2.36) for β = 1 4 and different values for the parameter e and (2.34) in the particular case where θ = 1, respectively. The computation are also performed for both the standard mass matrix M and the modified one M = M mod and they are compared to the analytical solution (u, λ) exhibited in Section 2. [START_REF] Galin | Contact problems in the theory of elasticity[END_REF] presented for standard matrix as well as for modified mass matrix. Notice that the slope of error curves for nonlinear problems reflect both the regularity of solution and the ability of numerical schemes to reproduce some technical inequalities. Finally, the convergence rates for different time-stepping methods, analyzed in details above, are summarized in Tables 2.1 and 2.2 for the reader convenience. These convergence rates were obtained for both standard (STD) and modified (MOD) mass matrices.

Scheme Method u h -u L ∞ (0,T ;H) u h -u L 2 (0,T ;H) u h -u L ∞ (0,T ;V) u h -u L 2 (0,T ;V) Crank-
Scheme Method λ h -λ L 2 (0,T ) Oscillations E h -E L ∞ (0,T ) E h -E L 2 (0,

Chapter 2. Numerical approximations of a one dimensional elastodynamic contact problem based on mass redistribution method

One can see in Table 2.1, that the rate of convergence of order one methods (Newmark I and backward Euler) in L p (0, T ; H) norms of the displacement (with p = 2, +∞) are limited to approximatively 1 2 for both the standard and modified mass matrices. Concerning the other methods, which are order two methods at least for the linear part of the problem, the rate of convergence is still approximatively 1 2 for the standard mass matrix but is improved to approximatively 1 for the modified mass matrix. Note that the Paoli-Schatzman scheme is not a fully order two method, except for e = 1, due to the discretization of the contact condition.

The advantage is less pronounced in L p (0, T ; V) norms of the displacement, but it is even more important in L 2 (0, T ) norm of the contact stress (see Table 2.2) since convergence does not occur for all the methods with the standard mass matrix.

An important remark is that, despite the low regularity of the exact solution, order two methods perform better than order one methods and should be preferred to approximate elastodynamic problem with impact.

Concerning the efficiency of the mass redistribution method, we can conclude that, even though it represents an additional approximation compared to the standard approximation, the convergence rates on the displacement and on the contact force are improved and the method reduces the potential spurious oscillations. The mass redistribution method appears then to be a robust method to approximate elastodynamic problem with impact.

The Crank-Nicolson method

We focus here on Newmark's method (2.32) in the particular case where (β, γ) = ( 14 , 1 2 ). The analytical solution (u, λ) exhibited in Section 2.5.1 and the approximated solutions (U n h , λ n ) obtained for different time and space steps are represented on Figure 2.7. We observe some spurious oscillations on Figure 2.7 (left) for the solution (U n h , λ n ) to Problem (2.32) after the contact takes place while these oscillations disappear when the standard mass matrix is replaced by the modified one, see Figure 2.7 (right).

It can easily be seen on Figure 2.8 that the scheme with the standard mass matrix is unstable with a rapidly growing energy while the one with the modified mass matrix is almost conservative as the space ∆x and time ∆t steps tend to 0. The error curves for (U n h , λ n ) highlighted that the norms u h -u L p (0,T ;H) , with p = 2, +∞, and λ n -λ H converge to 0 when M = M mod as n tends to +∞, which is unfortunately not the case for the standard mass matrix (see Figure 2.9). Finally, Figure 2.10 shows that the norm E h -E L p (0,T ) , p = 2, +∞, converges when M = M mod while this norm diverges when the standard mass matrix is considered.

The wave equation with Signorini and Dirichlet boundary conditions

(a) u h -u L ∞ (0,T ;H) (b) u h -u L 2 (0,T ;H) (c) u h -u L ∞ (0,T ;V) (d) u h -u L 2 (0,T ;V) (e) λ n -λ L 2 (0,T )
Figure 2.9: Comparison of the error curves obtained by using the standard and modified mass matrices.

The wave equation with Signorini and Dirichlet boundary conditions

(a) E h -E L ∞ (0,T ) (b) E h -E L 2 (0,T )
Figure 2.10: Comparison of the error curves of the energy obtained by using the standard and modified mass matrices.

The Newmark method I

We deal now with Newmark's method (2.32) in the particular case where (β, γ) = ( 12 , 1). We observe some oscillations on Figure 2.11 (left) for the Lagrange multiplier λ n associated to Problem (2.32) after the contact takes place while these oscillations are much less preeminent when the standard mass matrix is replaced by the modified one, see Figure 2.11 (right). Furthermore, the present scheme is dissipative and stable.

The Newmark method II

We consider Newmark's method (2.32) in the case where (β, γ) = ( 12 , 1 2 ). Some oscillations on Figure 2.13 (left) for the Lagrange multiplier λ n associated to Problem (2.32) can be observed after the contact takes place while these oscillations are much less important when the modified mass matrix is considered in (2.32), see Figure 2.13 (right). Note that here the scales for the Lagrange multiplier associated to standard and modified mass matrices are different, see We deal here with backward Euler's method (i.e. (2.34) with θ = 1). Since the present scheme is dissipative and stable (∆E n h < 0, see (2.35)), the approximated solutions (U n h , λ n ) obtained by using the standard mass as well as modified mass matrices do not oscillate as it can be observed in the previous numerical simulations. We are interested here in Paoli-Schatzman's method with (β, e) = ( 14 , 0). Note that the stability result straightforwardly follows from [START_REF] Dumont | Vibrations of a beam between obstacles. Convergence of a fully discretized approximation[END_REF]. We focus here on Paoli-Schatzman's method with (β, e) = ( 1 4 , 1 2 ). Observe that the stability result is still an open problem for e = 0. Some spurious oscillations can be observed for the approximated Lagrange multiplier as well as for the approximated energy. We consider now Paoli-Schatzman's method with (β, e) = ( 14 , 1). Note that this is a second order scheme.

Appendix D

The aim of the present section consists to establish the energy evolutions associated to schemes (2.32), (2.34) and (2.36). Let us observe that the energy evolution, when the external force f = 0, is given by (2.43)

∆E n h = 1 2 (U n+1 h,t -U n h,t ) T M (U n+1 h,t +U n h,t ) + 1 2 (U n+1 h -U n h ) T S(U n+1 h +U n h ).
We prove now that the energy evolution associated to the Crank-Nicolson method is given by(2.33). To this aim, we observe that (2.32a) and (2.32b) lead to (2.44)

M (U n+1 h -U n h -∆tU n h,t ) = -∆t 2 (( 1
2 -β ) SU n h +βSU n+1 h ) -∆t 2 (( 1
2 -β ) λ n e 0 +βλ n+1 e 0 ) ,

and

(2.45)

M (U n+1 h,t -U n h,t ) = -∆t((1-γ)SU n h +γSU n+1 h ) -∆t((1-γ)λ n e 0 +γλ n+1 e 0 ).
Then we multiply (2.45) by -∆t 2 and we subtract this expression to (

M (U n+1 h,t +U n h,t ) = 2 ∆t M (U n+1 h -U n h ) + 2∆t (( β- γ 2 ) S(U n+1 h -U n h ) ) + 2∆t (( β- γ 2 ) (λ n+1 e 0 -λ n e 0 ) 2.44), we get (2.46)
) .

By using (2.46) in (2.43), it follows that

∆E n h = 1 ∆t (U n+1 h,t -U n h,t ) T M (U n+1 h -U n h ) + ∆t ( β- γ 2 ) (U n+1 h,t -U n h,t ) T S(U n+1 h -U n h ) + ∆t ( β- γ 2 ) (U n+1 h,t -U n h,t ) T λ n+1 e 0 -λ n e 0 ) + 1 2 (U n+1 h -U n h ) T S(U n+1 h +U n h ).
Since M is a symmetric matrix, it follows from (2.45) that

∆E = ( 1 2 -γ)(U n+1 h -U n h ) T S(U n+1 h -U n h ) -(U n+1 h -U n h ) T ( (1 -γ)λ n e 0 + γλ n+1 e 0 ) ∆t(β - γ 2 )(U n+1 h,t -U n h,t ) T S(U n+1 h -U n h ) + ∆t(β - γ 2 )(U n+1 h,t -U n h,t ) T (λ n+1 e 0 -λ n e 0 ).
By setting the parameters (β, γ) = ( 1 4 , 1 2 ) (2.33) holds. Let us establish the energy evolution associated to the θ-method. We infer from (2.34a) and (2.34c) that (2.47)

M (U n+1 h,t -U n h,t ) = (1-θ)∆t(-SU n h -λ n e 0 ) + θ∆t(-SU n+1 h -λ n+1 e 0 ).
Then using (2.47) in (2.43), we find

∆E n h = - ∆t 2 (U n+1 h,t +U n h,t ) T ((1-θ)SU n h +θSU n+1 h +((1-θ)λ n +θλ n+1 )e 0 ) + 1 2 (U n+1 h -U n h ) T S(U n+1 h +U n h ).

The wave equation with Signorini and Dirichlet boundary conditions

Note that from (2.34a) ∆t 2 (U n+1 h,t +U n h,t ) = U n+1 h -U n h + ∆t ( 1 2 -θ ) (U n+1 h,t -U n h,t ), it follows that ∆E n h = ( 1 2 -θ ) (U n+1 h -U n h ) T S(U n+1 h -U n h ) -(U n+1 h -U n h ) T (((1-θ)λ n +θλ n+1 )e 0 ) + ∆t ( 1 2 -θ ) (U n+1 h,t -U n h,t ) T ((1-θ)(-SU n h -λ n e 0 ) + θ(-SU n+1 h -λ n+1 e 0 )).
Clearly we have ∆t

((1-θ)(-SU n h -λ n e 0 ) + θ(-SU n+1 h -λ n+1 e 0 )) = M (U n+1 h,t -U n h,t
), which allows us to conclude that (2.35) holds.

Finally we establish the energy evolution associated to (2.36) in the particular case where β = 1 4 . To this aim, let us introduce the energy evolution as follows:

∆E n h = [ 1 2 U T h,t M U h,t + 1 2 U T h SU h ] n+ 1 2 n-1 2
which leads to

∆E n h = 1 2 (U n+ 1 2 h,t -U n-1 2 h,t ) T M (U n+ 1 2 h,t +U n-1 2 h,t ) + 1 2 (U n+ 1 2 h -U n-1 2 h ) T S(U n+ 1 2 h +U n-1 2 h ),
where U

n+ 1 2 h def = 1 2 (U n h +U n+1 h ) and U n-1 2 h def = 1 2 (U n h +U n-1 h ).
Then it comes that

U n+ 1 2 h + U n-1 2 h = 1 2 (U n+1 h +2U n h +U n-1 h ) and U n+ 1 2 h -U n-1 2 h = 1 2 (U n+1 h -U n-1 h ),
and also

U n+1 h = U n h + ∆tU n+ 1 2 h,t and U n h = U n-1 h + ∆tU n-1 2 h,t , so we have U n+ 1 2 h,t -U n-1 2 h,t = U n+1 h -2U n h +U n-1 h ∆t and U n+ 1 2 h,t + U n-1 2 h,t = U n+1 h -U n-1 h ∆t .
Obviously, we get

∆E n h = 1 2 (U n+1 h -U n-1 h ) T M ( U n+1 h -2U n h +U n-1 h ∆t 2 ) + S ( U n+1 h +2U n h +U n-1 h 4
) .

According to (2.36) together with β = 1 4 , we deduce that

∆E n h = 1 2 (U n+1 h -U n-1 h ) T (-λ n+1 e 0 ) = 1 2 (u n+1 0 -u n-1 0 )(-λ n+1 ). Since u n+1 0 -u n-1 0 = (1+e)u n,e 0 -u n-1 0
, we find Convergence curves for the Newmark method II

∆E n h = 1 2 (1+e) ( (-λ n+1 u n,e 0 ) -(-λ n+1 u n-1 0 ) ) = 1 2 (1+e)λ n+1 u n-1 0 .

The wave equation with Signorini and Dirichlet boundary conditions

(a) In E h -E L ∞ (0,T ) (b) In E h -E L 2 (0,T )
(a) U n h -u L ∞ (0,T ;L 2 (0,L)) (b) U n h -u L 2 (0,T ;L 2 (0,L)) (c) U n h -u L ∞ (0,T ;H 1 (0,L))
Chapter 2. Numerical approximations of a one dimensional elastodynamic contact problem based on mass redistribution method

(d) U n h -u L 2 (0,T ;H 1 (0,L)) (e) λ n -λ L 2 (0,T )
Figure 2.25: Comparison of the convergence curves obtained by using standard and modified mass matrices. Convergence curves for the backward Euler method

(a) E h -E L ∞ (0,T ) (b) E h -E L 2 (0,T )
(a) U n h -u L ∞ (0,T ;L 2 (0,L)) (b) U n h -u L 2 (0,T ;L 2 (0,L))
2.5. The wave equation with Signorini and Dirichlet boundary conditions

(c) U n h -u L ∞ (0,T ;H 1 (0,L)) (d) U n h -u L 2 (0,T ;H 1 (0,L)) (e) λ n -λ L 2 (0,T )
a) E h -E L ∞ (0,T ) (b) E h -E L 2 (0,T )
Figure 2.28: Comparison of the convergence curves of the energy obtained by using standard and modified mass matrices.

Chapter 2. Numerical approximations of a one dimensional elastodynamic contact problem based on mass redistribution method

Convergence curves for the Paoli-Schatzman method I

(a) U n h -u L ∞ (0,T ;L 2 (0,L)) (b) U n h -u L 2 (0,T ;L 2 (0,L)) (c) U n h -u L ∞ (0,T ;H 1 (0,L)) (d) U n h -u L 2 (0,T ;H 1 (0,L)) (e) λ n -λ L 2 (0,T )
Figure 2.29: Comparison of the convergence curves obtained by using standard and modified mass matrices. Convergence curves for the Paoli-Schatzman method II

The wave equation with Signorini and Dirichlet boundary conditions

(a) E h -E L ∞ (0,T ) (b) E h -E L 2 (0,T )
(a) U n h -u L ∞ (0,T ;L 2 (0,L)) (b) U n h -u L 2 (0,T ;L 2 (0,L)) (c) U n h -u L ∞ (0,T ;H 1 (0,L)) (d) U n h -u L 2 (0,T ;H 1 (0,L))
Chapter 2. Numerical approximations of a one dimensional elastodynamic contact problem based on mass redistribution method

(e) λ n -λ L 2 (0,T )
Figure 2.31: Comparison of the convergence curves obtained by using standard and modified mass matrices. Convergence curves for the Paoli-Schatzman method III

(a) E h -E L ∞ (0,T ) (b) E h -E L 2 (0,T )
(a) U n h -u L ∞ (0,T ;L 2 (0,L)) (b) U n h -u L 2 (0,T ;L 2 (0,L))
2.5. The wave equation with Signorini and Dirichlet boundary conditions

(c) U n h -u L ∞ (0,T ;H 1 (0,L)) (d) U n h -u L 2 (0,T ;H 1 (0,L)) (e) λ n -λ L 2 (0,T )
a) E h -E L ∞ (0,T ) (b) E h -E L 2 (0,T )
Figure 2.34: Comparison of the convergence curves of the energy obtained by using standard and modified mass matrices.

Chapter 3

Efficiency of a weighted mass redistribution method for an elastodynamic problem with unilateral constraints

Introduction

The motion of an elastic bar of length L which is free to move as long as it it does not hit a material obstacle is studied, see Figure 3.1. The assumptions of small deformations are assumed and the material of the bar is supposed to be homogeneous. Let u(x, t) be the displacement

at time t ∈ [0, T ], T > 0 of the material point of spatial coordinate x ∈ [0, L]. Let f (x, t)
denotes a density of external forces, depending on time and space. The mathematical problem is formulated as follows:

(3.1)

u tt (x, t) -u xx (x, t) = f (x, t), (x, t) ∈ (0, L) × (0, T ),
with Cauchy initial data

(3.2) u(x, 0) = u 0 (x) and u t (x, 0) = v 0 (x), x ∈ (0, L),
) 0 ≤ u(0, t) ⊥ u x (0, t) ≤ 0 and u(L, t) = 0, t ∈ [0, T ]. (3.3
The orthogonality has the natural meaning: an appropriate duality product between two terms of relation vanishes.

Let us describe the weak formulation associated to (

. For that purpose, it is convenient to introduce the following notations:

V def = {u ∈ H 1 (0, L) : u(L) = 0}, H def = L 2 (0, L), V def = {u ∈ L 2 (0, T ; V) : u t ∈ L 2 (0, T ; H)} and the convex set K def = {u ∈ V : u(0) ≥ 0}.
                 find u ∈ K such that - ∫ L 0 v 0 (x)(v(x, 0)-u 0 (x)) dx - ∫ Q T u t (x, t)(v t (x, t)-u t (x, t)) dx dt+ ∫ Q T u x (x, t)(v x (x, t)-u x (x, t)) dx dt ≥ ∫ Q T f (x, t)(v(x, t)-u(x, t)) dx dt for all v ∈ K for which there exists ζ > 0 with v = u for t ≥ T -ζ.
Existence and uniqueness results are obtained for a similar situation of a vibrating string with concave obstacle in one dimensional space in [START_REF] Schatzman | A hyperbolic problem of second order with unilateral constraints: The vibrating string with a concave obstacle[END_REF] and also for a wave equation with unilateral constraint at the boundary in a half-space of R N in [START_REF] Lebeau | A wave problem in a half-space with a unilateral constraint at the boundary[END_REF]. An existence result for a wave equation in a C 2 -regular bounded domain constrained by an obstacle at the boundary in R N for N ≥ 2 is proven in [START_REF] Kim | A boundary thin obstacle problem for a wave equation[END_REF]. The reader is also referred to [START_REF] Dabaghi | Convergence of mass redistribution method for the one-dimensional wave equation with a unilateral constraint at the boundary[END_REF].

The paper is organized as follows. In Section 3.2, a space semi-discretization based on mass redistribution method is presented. This mass redistribution method consists to transfer the mass of the contact node on the other nodes while the inertia vanishes at the contact node. The error estimate in time as well as the convergence result are established. A benchmark problem is introduced in Section 3.3 and its analytical solution is exhibited. Then numerical experiments for some space-time discretizations like the Crank-Nicolson or the backward Euler methods are reported. These numerical experiments highlight that the choice of the nodes where the mass is transferred plays a crucial role to get a better approximate solution.

Finite element discretization and convergence of the mass redistributions method

Finite element discretization and convergence of the mass redistributions method

This section is devoted to semi-discrete problems in space associated to (3.4) by using the mass redistributions method, see [START_REF] Khenous | Mass redistribution method for finite element contact problems in elastodynamics[END_REF][START_REF] Dabaghi | Convergence of mass redistribution method for the one-dimensional wave equation with a unilateral constraint at the boundary[END_REF]. More precisely, the mass redistributions method consists in transferring the mass of the contact node on the other nodes implying that the node at the contact boundary evolves in a quasi-static way. To this aim, we choose a set of parameters h def = L m (mesh size) having in mind the limit m tends to +∞, where m is an integer, and let

V h def = {v h ∈ C 0 ([0, L]) : v h | [ih,(i+1)h] ∈ P 1 , i = 0, . . . , m -1, v h (L) = 0}
where P 1 is the space of polynomials of degree less than or equal to 1. A classical basis of V h is given by the sequence of shape functions ϕ i ∈ V h for i = 0, 1, . . . , J, defined by

ϕ i (x) def = { 1 -|x-ih| h if x ∈ [a max(i-1,0) , a min(i+1,m) ], 0 otherwise,
where J = m -1. Clearly, ϕ i (jh) = δ ij , j = 0, 1, . . . , m, where δ denotes Kronecker's symbol.

Thus the solution u belonging to V to the weak formulation (3.4) is approximated by

u h (x, t) = J ∑ j=0 u j (t)ϕ j (x).
Consequently, u i = u h (a i ), i = 0, 1, . . . , n -1 and the weak formulation (3.4) is approximated by

(P u h )                find u h : [0, T ] → V h and λ h : [0, T ] → R such that for all v h ∈ V h ∫ L 0 u h,tt v h dx + ∫ L 0 u h,x v h,x dx = -λ h v h (0) + ∫ L 0 f v h dx a.e. t ∈ [0, T ], 0 ≤ u h (0, •) ⊥ λ h ≤ 0 a.e. t ∈ [0, T ], u h (•, 0) = u 0 h and u h,t (•, 0) = v 0 h ,
where u 0 h and v 0 h belong to V h and λ is the Lagrange multiplier being here the contact force. The corresponding algebraic formulation is given by

(P U h λ h )            find U h : [0, T ] → R n and λ h : [0, T ] → R such that M U h,tt + SU h = -λe 0 + F a.e. t ∈ [0, T ], 0 ≤ u 0 ⊥ λ h ≤ 0 a.e. t ∈ [0, T ], U h (0) = U 0 and U h,t (0) = V 0 ,
where M and S denote the mass and stiffness matrices, respectively, and

U h def = (u 0 , . . . , u n-1 ) T , e 0 def = (1, 0, . . . , 0) T and F def = ( ∫ L 0 f ϕ 0 dx, . . . , ∫ L 0 f ϕ n-1 dx) T with (•)
T being the transpose of a tensor.

An alternative approach to the standard discretization presented above is to consider the mass redistributions method which consists to replace the mass matrix M in (P U h λ h ) by a modified mass matrix M mod defined by

M mod ij def = ∫ L 0 ϕ i ϕ j w h,ij dx such that w h,0i = w h,i0 for all
Chapter 3. Efficiency of a weighted mass redistribution method for an elastodynamic problem with unilateral constraints i = 0, . . . , n -1 and w h,ij = 1 + o(h) for all i, j = 1, . . . , n -1. Clearly this modified mass matrix reads

M mod def = ( 0 0 0 M ) .
This leads to an algebraic formulation of the semi-discrete approximation with mass redistribution method given by

(P mod U h λ h )                  find U h : [0, T ] → R n and λ : [0, T ] → R such that M Ūh,tt + S Ūh = -Φ 0 u 0 + F a.e. t ∈ [0, T ], S 00 u 0 + Φ T 0 Ūh = -λ h + F 0 a.e. t ∈ [0, T ], 0 ≤ u 0 ⊥ λ h ≤ 0 a.e. t ∈ [0, T ], U h (0) = U 0 and U h,t (0) = V 0 .
Here

Ū def = (u 1 , . . . , u n-1 ) T , Sij def = S i+1,j+1 with Φ 0 def = ∫ Ω ϕ i+1 ϕ 0 dx, i = 0, . . . , n -2, F def = ( ∫ L 0 f ϕ 1 dx, . . . , ∫ L 0 f ϕ n-1 dx) T . Since Φ 0 = (S 10 , 0, . . . , 0) T , it comes that u 0 = ( -λ h -Φ T 0 Ūh S 00 ) = ( -λ h -S 10 u 1 S 00
)

. If we assume that S 10 u 1 ≥ 0 then the compatibility condition implies that u 0 = 0 and λ h = (Φ T 0 Ūh ) -otherwise λ = 0. This gives that u 0 = ( -S 10 u 1 S 00

) + , and then we may conclude that (P mod U h λ h ) is equivalent to the following second order Lipschitz continuous ordinary differential equation:

(P mod Ūh )          find Ūh : [0, T ] → R n-1 such that M Ūh,tt + S Ūh = -Φ 0 ( -S 10 u 1 S 00 ) + + F a.e. t ∈ [0, T ], U h (0) = U 0 and U h,t (0) = V 0 .
It is convenient to introduce the following notations:

Ūh def = ( Ūh , Ūh,t ) T and G( Ūh ) def = ( Ūh,t , -M -1 S Ūh + M -1 D Ūh + M -1 F ) T with D Ūh def = ((u 1 ) + /h, 0, . . . , 0) T . Therefore (P mod

Ūh

) can be rewritten as follows:

(P mod Ūh )      find Ūh : [0, T ] → R n-1 × R n-1 such that Ūh,t = G( Ūh ) a.e. t ∈ [0, T ], Ūh (0) = (U 0 , V 0 ) T .
Observe that G( Ūh ) is Lipschitz continuous. Indeed, let Ū1 h and Ū2 h be two solutions to Problem (P mod

Ūh

), then we have

G( Ū1 h ) -G( Ū2 h ) 2 H ≤ max(1, 1/h 2 , ||| S||| 2 ) Ū1 h -Ū2 h 2 H ,
where ||| • ||| denotes the induced matrix norm. Existence and uniqueness results for the problem (P mod

Ūh

) follow from the Lipschitz continuity of G( Ūh ), for further details the reader is referred to [START_REF] Crouzeix | Analyse numérique des équations différentielles[END_REF]. In particular, we have Ūh ∈ C 1 ([0, T ]).

Lemma 11. Let ∆t be the time step and N = T ∆t ∈ N for given T > 0. Then the truncation error for the Crank-Nicolson method to solve the semi-discrete problem (P mod

Ūh

) is of order ∆t.

Finite element discretization and convergence of the mass redistributions method

Proof. Let the time discretization [0, T ] be divided by n + 1 discrete time-points such that 0 = t 0 < t 1 < . . . < t n = T . Furthermore the discrete quantities Ūn h and Ūn h,t are assumed to be given by algorithmic approximations of the displacement Ūh (t n ) and the velocity Ūh,t (t n ), respectively. Notice that the Crank-Nicolson method is based on the trapezoidal rule in time. Thus the straightforward application of the Crank-Nicolson method to Problem (P mod

Ūh

) is a combination of the forward Euler method at n and the backward Euler method at n + 1, we obtain

(3.5) Ūn+1 h -Ūn h ∆t = 1 2 (G( Ūn+1 h ) -G( Ū n h )).
Since Ū ∈ C 1 ([0, T ]), we may deduce from (P mod Ū ) that

(3.6) Ūh (t n+1 ) -Ūh (t n ) ∆t = 1 2 (G( Ūh (t n+1 )) -G( Ūh (t n ))) + ε n (∆t),
where ε n is the local truncation error such that

(3.7) ε n (∆t) H ≤ C F ∆t 2 with C F def = sup [t 0 ,t 0 +T ] F ( Ū) H .
Let us define the global truncation error ēn def = Ūn h -Ūh (t n ). Then we substract (3.6) from (3.5) to get ēn+1 = ēn + ∆t 2

( (G( Ūn+1 h ) -G( Ūh (t n+1 ))) -(G( Ūn h ) -G( Ūh (t n ))) ) -ε n (∆t 2 ).
By using the Euclidean inner product (•, •), we get

(3.8) (ē n+1 , ēn+1 ) =(ē n+1 , ēn ) + ∆t 2 ( ēn+1 , G( Ūn+1 h ) -G( Ūh (t n+1 )) ) - ∆t 2 ( ēn+1 , G( Ūn h ) -G( Ūh (t n )) ) - ( ēn+1 , ε n (∆t 2 )
) .

Let us define

) ēn+1 H ≤ ( 1 + C h S ∆t 2
) ēn H + ε n (∆t 2 ) H .
According to (3.7), for ∆t sufficiently small, we find

(3.9) ēn+1 H ≤ ( 1 + C h S ∆t C h,∆t S ) ēn H + C F ∆t 2 2C h,∆t S where C h,∆t S def = 1 - C h S ∆t
2 . We may infer from (3.9) that

ēn+1 H ≤ ( 1 + C h S ∆t C h,∆t S ) n+1 ē0 H + C F 2 n ∑ i=0 ( 1 + C h S ∆t C h,∆t S ) i ( ∆t 2 C h,∆t S ) ,
Chapter 3. Efficiency of a weighted mass redistribution method for an elastodynamic problem with unilateral constraints which implies that

ēn+1 H ≤ ( 1 + C h S ∆t C h,∆t S ) n+1 ē0 H + C F 2 ( 1 - ( 1 + C h S ∆t C h,∆t S ) n+1 1 - ( 1 + C h S ∆t C h,∆t S ) ) ( ∆t 2 C h,∆t

S

) .

Then the Taylor expansion leads to

ēn+1 H ≤ exp ( C h S T C h,∆t S ) ē0 H + C F ∆t 2C h S ( 1 -exp ( C h S T C h,∆t S ))
.

It follows that

ēn+1 H ≤ exp ( C h S T C h,∆t S )( ē0 H + C F ∆t 2C h S )
.

Let us denote the spectral radius of Id+ S by ρ(Id+ S). Clearly, we have C h S ≥ max(1, 1/h)ρ(Id+ S) ≥ max(1, 1/h) and it comes that

ēn+1 H ≤ exp ( C h S T C h,∆t S )( ē0 H + C F ∆t max(1, 1/h)
)

.

which concludes the proof.

On the other, we may observe that (P mod U λ ) is equivalent to

(P mod u h )                  find u h : [0, T ] → V h such that for all v h ∈ K ∩ V h ∫ T 0 ∫ L h u h,tt (x, t)(v h (x, t)-u h (x, t))w h (x) dx dt + ∫ Q T u h,x (x, t)(v h,x (x, t)-u h,x (x, t)) dx dt ≥ ∫ Q T f (x, t)(v h (x, t)-u h (x, t)) dx dt, u h (x, 0) = u 0 h (x) and u h,t (x, 0) = v 0 h (x).
We assume that the initial data u 0 h and u h,t satisfy

(3.10) lim h→0 ( u 0 h -u 0 V + v 0 h -v 0 H ) = 0.
The convergence of the solution u h of (P mod u h ) to the solution of (3.4) is proved below. To this aim, the same techniques detailed in the proof of theorem 4.3 in [START_REF] Dabaghi | Convergence of mass redistribution method for the one-dimensional wave equation with a unilateral constraint at the boundary[END_REF] are used. The main difference comes from the weight functions that are introduced in the present work. Indeed the mass of the contact node is redistributed on the other nodes by using these weight functions while this mass vanishes in [START_REF] Dabaghi | Convergence of mass redistribution method for the one-dimensional wave equation with a unilateral constraint at the boundary[END_REF]. The reader is also referred to [START_REF] Schatzman | Numerical approximation of a wave equation with unilateral constraints[END_REF]. Furthermore, we may observe that for all τ ∈ [0, T ], the following energy relation (3.11)

∫ L 0 (|u t (x, τ )| 2 +|u x (x, τ )| 2 ) dx = ∫ L 0 (|v 0 (x)| 2 +|u 0 x (x)| 2 ) dx + ∫ Qτ f (x, t)u t (x, t) dx dt holds.

Finite element discretization and convergence of the mass redistributions method

Theorem 11. Assume that (3.10) holds. Then, the solution u h of (P mod u h ) converges in the strong topology of V to the unique solution of (3.4) as h tends to 0.

Proof. We observe that the energy relation (3.11) and the Gronwall lemma yields that there exists C > 0 independent of τ such that

sup t∈[0,T ] ( u h (•, τ ) V + u h,t (•, τ ) H ) ≤ C.
Then we may extract a subsequence, still denoted by u h , such that

u h u in L ∞ (0, T ; V) weak * , (3.12a) u h,t u t in L ∞ (0, T ; H) weak * . (3.12b) Let us define H ∞ def = {u ∈ L ∞ (0, T ; V) : u t ∈ L ∞ (0, T ; H)} endowed with the norm u H∞ def = ess sup t∈[0,T ] ( u(•, t) V + u t (•, t) H )
. Then, we may deduce from (3.12) that

(3.13) u h u in H ∞ weak * .
Notice that for all α < 1 2 , we have [START_REF] Schatzman | Numerical approximation of a wave equation with unilateral constraints[END_REF]), where → and → → denote the continuous and compact embeddings, respectively. Finally we find

H ∞ → C 0, 1 2 (Q T ) → → C 0,α (Q T ) hold (see
(3.14) u h → u in C 0,α (Q T )
for all α < 1 2 . Furthermore, u h and u belong to K. Our aim consists to establish that the limit u satisfies (3.4). However the elements of K are not smooth enough in time, then they should be approximated before being projected onto V h . Indeed this projection violates the constraint at x = 0, and therefore, the elements of K need another approximation in order to satisfy the constraint strictly. Assume that v ∈ K such that v = u for t ≥ T -ε and for η ≤ ε 4 , let us define

v η (x, t) def =    u(x, t) + 1 η ∫ t+η t (v(x, s)-u(x, s)) ds + k(η)(L-x)ψ(t) if t ≤ T -η, u(x, t) if t ≥ T -η,
where ψ is a smooth and positive function defined by

ψ = 1 on [0, T -ε 2 ] and ψ = 0 on [T -ε 4 , T ]. We precise now how the parameter k(η) to ensure that v η ∈ K ∩ L ∞ (0, T ; V) is chosen. Since u ∈ C 0, 1 2 (Q T ), it follows that there exists C k > 0, such that u(0, t) - 1 η ∫ t+η t u(0, s) ds ≤ 1 η ∫ t+η t |u(0, t)-u(0, s)| ds ≤ C k u H∞ η ∫ η 0 √ s ds = 2 3 C k u H∞ √ η.
Then for all t ≤ T -ε 2 , we obtain

v η (0, t) ≥ 1 η ∫ t+η t v(0, s) ds - 2 3 C k u H∞ √ η + k(η)Lψ(t).
Chapter 3. Efficiency of a weighted mass redistribution method for an elastodynamic problem with unilateral constraints

Observe that the choice k(η) = 5 3L C k u H∞ √ η ensures that for all t ≤ T -ε 2 , we get

v η (0, t) ≥ C k u H∞ √ η.
Let Proj V h be the orthogonal projection onto V h with respect to the scalar product in H. Since the sequence Proj V h satisfies Proj V h z-z V → 0 when h → 0 for all z ∈ V, the Sobolev injections imply that there exists a sequence α h converging to zero as h tends to zero such that for all z belonging to V, we have

Proj V h z-z C 0 ≤ α h z V .
The next step consists in choosing an adequate test function. Let us define

(3.15) v h (•, t) def = u h (•, t) + Proj V h (v η -u)(•, t)
for all t belonging to [0, T ]. Clearly, v h ∈ K∩V h , for all t, and for all h small enough. Introducing (3.15) into (P mod u h ) and according to the integration by parts, it comes that

- ∫ L 0 u h,t (x, 0)Proj V h (v η -u)(x, 0)w h (x) dx - ∫ Q T u h,t (x, t)Proj V h (v η,t -u t )(x, t)w h (x) dx dt + ∫ Q T u h,x (x, t)(Proj V h (v η -u)) x (x, t) dx dt ≥ ∫ Q T f (x, t)Proj V h (v η -u)(x, t) dx dt.
We already know that

Proj V h (v η,t -u t ) → v η,t -u t in L 2 (0, T ; H), (3.16a) Proj V h (v η -u) → v η -u in L 2 (0, T ; V). (3.16b)
According to (3.12)- (3.16), we are able to pass to the limit as h tends to 0, we find

- ∫ L 0 v 0 (v η -u)(x, 0) dx - ∫ Q T u t (x, t)(v η,t -u t )(x, t) dx dt + ∫ Q T u x (x, t)(v η -u)) x (x, t) dx dt ≥ ∫ Q T f (x, t)(v η -u)(x, t) dx dt.
Since it is quite a routine to pass to the limit with respect to η to get the variational formulation (3.4), the verification is left to the reader. It remains to prove that u h converges strongly in V. Notice that (P mod u h ) satisfies the following energy relation

∫ L 0 (|u h,t (x, τ )| 2 w h (x)+|u h,x (x, τ )| 2 ) dx = ∫ L 0 (|v 0 (x)| 2 w h (x)+|u 0 h,x (x)| 2 ) dx + ∫ Qτ f (x, t)u h,t (x, t) dx dt.
According to (3.10) and to the convergence results established above, we get (3.17)

lim h→0 ∫ L 0 (|u h,t (x, τ )| 2 w h (x)+|u h,x (x, τ )| 2 ) dx = ∫ L 0 (|v 0 (x)| 2 +|u 0 h,x (x)| 2 ) dx + ∫ Qτ f (x,

The wave equation with Signorini and Dirichlet boundary conditions

We consider a bar of length L = 1 clamped at one end and compressed at t = 0. The bar elongates under the gravity effect; as soon as it reaches a rigid obstacle at time t 1 then it stays in contact during the time t 2 -t 1 and it takes off at time t 2 , see Figure 3.2. This problem can be described mathematically by (3.1)-(3.3) with the density of external forces f (x, t) = 0.

We first describe how an analytical piecewise affine and periodic solution to our problem can be obtained by using the characteristics method, the reader is referred to [START_REF] Dabaghi | Numerical approximations of elastodynamic contact problems based on mass redistribution method[END_REF] for a detailed explanation. Then approximate solutions for some time-space discretizations and for several mass redistributions are exhibited and their efficiency is discussed.

Analytical solution

The domains considered here are defined by (0, L) × (t i , t i+1 ), i = 0, 1, 2, corresponding to the phases before, during and after the impact, respectively. Each of them are divided into four regions as it is represented on Figure 3.2. We choose below t i = i, i = 0, . . . , 3. The domain (0, 1) × (0, 1) corresponding to the phase before the impact is split into four regions according to the characteristics lines x + t and x -t. Therefore Chapter 3. Efficiency of a weighted mass redistribution method for an elastodynamic problem with unilateral constraints

v = 0 u = 1-x 2 t 0 = 0 v = -1/2 u = 0 t 1 = 1 v = 1/2 u = 0 t 2 = 2 v = 0 u = 1-x 2 t 3 = 3
The domain (1, 2) × (0, 1) corresponding to the phase during the impact is also divided into four regions. Then the solution (3.18) evaluated in the region IV allows us to infer that u t (•, 1) = -1 2 and we conclude that

(3.19) u(x, t) =        t-t 1
2 in the region I, -x 2 in the region II,

x-1 2 in the region III,

t-t 1 -1 2
in the region IV.

The domain (2, 3) × (0, 1) corresponding to the phase after the impact is split into four regions. By using (3.19), we get u(•, 2) = 0 and u t (•, 2) = 1 2 which leads to

u(x, t) = { t-t 2
2 in the regions I, II, 1-x 2 in the regions III, IV.

Since u(•, 3) = u 0 and u t (•, 3) = v 0 , the solution u(x, t) is periodic of period 3. Finally, note that λ = u x (0, •).

Comparisons between different mass redistributions for some timespace discretisations

The time discretization is introduced in this section. To this aim, we divide the time interval [0, T ] by n + 1 discrete time-points such that 0 = t 0 < t 1 < . . . < t n = T and let U n h , U n h,t , U n h,tt and λ n the approximations of the displacement U h (t n ), the velocity U h,t (t n ), the acceleration U h,tt (t n ) and the Lagrange multiplier λ(t n ), respectively. We deal with some approximate solutions to Problem (3.1)-(3.3) obtained by using several time-stepping methods like the Newmark, backward Euler and Paoli-Schatzman methods. For each of these time-stepping methods, approximate solutions (U n h , λ n ) are exhibited for several mass redistributions and they are compared to the analytical solution (u, λ) introduced in Section 3.3.1. In the numerical experiments presented below, Mod1, Mod2 and Mod3 correspond to the mass of the contact node that is not redistributed, that is redistributed on all the other nodes and that is redistributed on the node preceding the contact node, respectively. The efficiency of the mass redistribution method depends on the position of the nodes where the mass is redistributed. Indeed, the numerical experiments highlight that the closer from the contact node the mass is transferred the better the approximate solutions. Then, it is not surprising that the best approximate solution can be expected and indeed it is obtained in the case where all the mass of the contact node is transferred on the node preceding the contact node, see Figures 3.3, 3.5 and 3.7. Note that the numerical simulations presented below were performed by employing the finite element library Getfem++ (see [START_REF] Renard | Getfem++. An Open Source generic C++ library for finite element methods[END_REF]).

The Newmark methods

The Taylor expansions of displacements and velocities neglecting terms of higher order are the underlying concept of the family of Newmark methods, see [START_REF] Newmark | A method of computational for structural dynamics[END_REF]. These methods are unconditionally stable for linear elastodynamic problem for γ ≥ 1 2 and β ≥ 1 4 ( 1 2 +γ) 2 , see [START_REF] Hughes | Linear static and dynamic finite element analysis[END_REF][START_REF] Krenk | Energy conservation in Newmark based time integration algorithms[END_REF], but they are also the most popular time-stepping schemes used to solve contact problems. (3.20)

The wave equation with Signorini and Dirichlet boundary conditions

                 find U n+1 h : [0, T ] → R n and λ n : [0, T ] → R such that: U n+1 h = U n h + ∆tU n h,t + ( 1 2 -β ) ∆t 2 U n h,tt + β∆t 2 U n+1 h,tt , U n+1 h,t = U n h,t + (1-γ)∆tU n h,tt + γ∆tU n+1 h,tt , M U n+1 h,tt + SU n+1 h = -λ n+1 e 0 + F n+1 , 0 ≤ u n+1 0 ⊥ λ n+1 ≤ 0,
where ∆t is a given times step and (β, γ) are the algorithmic parameters, see [START_REF] Hughes | Linear static and dynamic finite element analysis[END_REF][START_REF] Laursen | Computational contact and impact mechanics. Fundamentals of modeling interfacial phenomena in nonlinear finite element analysis[END_REF]. Note that U 0 h , U 0 h,t and λ 0 are given and U 0 h,tt is evaluated by using the third identity in (3.20). We are particularly interested in the case where (β, γ) = ( 14 , 1 2 ). This method is called the Crank-Nicolson method, it is second-order consistent and unconditionally stable in the unconstrained case. However the situation is quite different in the case of contact constraints, indeed the order of accuracy is degraded, for further details, the reader is referred to [START_REF] Hughes | Linear static and dynamic finite element analysis[END_REF][START_REF] Krenk | Energy conservation in Newmark based time integration algorithms[END_REF][START_REF] Grosu | Stability of semidiscrete formulations for elastodynamics at small time steps[END_REF]. The analytical solution (u, λ) exhibited in Section 3.3.1 and the approximate solutions (U n h , λ n ) obtained for different mass redistributions are represented on Figure 3.3. The approximate solution obtained for the redistribution of the contact node mass on the node preceding the contact (see Mod 3 on Figure 3.3) gives much better accuracy than the mass redistribution on all the nodes preceding the contact node (see Mod 2 on Figure 3.3) and no mass redistribution. This highlighted that the choice for the mass redistribution plays a crucial role.

The backward Euler method

We are concerned here with backward Euler's method which for the contact problem (3.1)-(3.3) is described by the following finite difference equations:

(3.21)                  find U n+1 h : [0, T ] → R n and λ n : [0, T ] → R such that: U n+1 h = U n h + ∆tU n+1 h,t , U n+1 h,t = U n h,t + ∆tU n+1 h,tt , M U n+1 h,tt + SU n+1 h = -λ n+1 e 0 + F n+1 , 0 ≤ u n+1 0 ⊥ λ n+1 ≤ 0.
Note that U 0 h , U 0 h,t and λ 0 are given and U 0 h,tt is evaluated by using the third equality in (3.21).

The wave equation with Signorini and Dirichlet boundary conditions

The Paoli-Schatzman methods

We focus on the so-called Paoli-Schatzman method that consists to fix the contact constraint at an intermediate time step. Indeed the method proposed below is a slight modification of Chapter 3. Efficiency of a weighted mass redistribution method for an elastodynamic problem with unilateral constraints Paoli-Schatzman method (see [START_REF] Paoli | Time discretization of vibro-impact[END_REF][START_REF] Paoli | A numerical scheme for impact problems. I ,II. The multidimensional case[END_REF]) which takes into account the kernel of the modified mass matrix. A simple application of Paoli-Schatzman method based on Newmark scheme to Problem (P Ūh ,λ h with γ = 1 2 leads to

(3.22)                  find U n+1 h : [0, T ] → R n and λ n : [0, T ] → R such that: M (U n+1 h -2U n h +U n-1 h ) ∆t 2 + S(βU n+1 h +(1-2β)U n h +βU n-1 h ) = -λ n e 0 for all n ≥ 2, 0 ≤ u n,e 0 = u n+1 0 + eu n-1 0 1 + e ⊥ λ n ≤ 0, U 0 and U 1 given.
Here e belongs to [0, 1] and is aimed to be interpreted as a restitution coefficient. Note that U 0 h and U 1 h are given data and U 1 h can be evaluated by a one step scheme. We may observe that taking M = M mod in (3.22), we are not able to resolve the problem on the kernel of M mod . That is the reason why, SU n-1 h as well as SU n h are projected on the orthogonal of the kernel of M . We are interested here in Paoli-Schatzman's method with (β, e) = ( 14 , 0). Note that the stability result straightforwardly follows from [START_REF] Dumont | Vibrations of a beam between obstacles. Convergence of a fully discretized approximation[END_REF].

Introduction

This note aims to give some numerical results for the evolution of an elastic material being subjected to unilateral boundary conditions by employing the mass redistribution method introduced in [START_REF] Khenous | Mass redistribution method for finite element contact problems in elastodynamics[END_REF]. The contact phenomena is modelled by using the so-called Signorini's boundary conditions in displacement, which are based on a linearization of the physically meaningful non penetrability of masses. A considerable engineering and mathematical literature is devoted to contact problems, however a few existence of solutions results have been established, the reader is referred to [START_REF] Lebeau | A wave problem in a half-space with a unilateral constraint at the boundary[END_REF][START_REF] Kim | A boundary thin obstacle problem for a wave equation[END_REF]. Therefore a challenging task consists to elaborate efficient numerical methods able to approximate this problem. The present work adopts the approach consisting to remove the mass at the contact nodes which prevents the oscillations at the contact boundaries as already observed for the one dimensional elastodynamic contact problems in [START_REF] Dabaghi | Numerical approximations of elastodynamic contact problems based on mass redistribution method[END_REF].

We consider an elastic bar vibrating vertically such that one end of this bar is free to move, as long as it does not hit a material obstacle, while the other end is clamped . The obstacle constrains the displacement of the extremity to be greater than or equal to 0. We have supposed

def = [0, L 1 ] × [0, L 2 ] and Ω 3 def = [0, L 1 ] × [0, L 2 ] × [0, L 3 ] together with L i > 0, i = 1, 2, 3.
We denote by u(x, t) the displacement field at time t ∈ [0, T ], T > 0, of the material point of spatial coordinate x ∈ Ω ⊂ R d . Furthermore, let f (x, t) be the density of external forces depending on time and space. The mathematical problem is formulated as follows:

(4.1) ρu tt -div σ(u) = f in Ω d × (0, T ),
where ρ > 0 is the mass density, (•)

t def = ∂(
0 ≥ ν T u ⊥ ν T (σ(u)ν) ≤ 0 and τ T (σ(u)ν) = 0 on Γ Sig d × (0, T ), (4.3b)
where ν and τ denote the outward unit normal and tangential vectors and (•) T is the transpose of a tensor. Here the orthogonality has the natural meaning; namely if we have enough regularity, it means that the product ν T u(ν T σ(u)ν) vanishes almost everywhere at the contact boundary. If it is not the case, the above inequality is integrated on an appropriate set of test functions, leading to a weak formulation for the unilateral condition. Let us describe now the functional hypotheses on the data; if X is a space of scalar functions, the bold-face notation X denotes the space X d . For the final result, we require the initial displacement u 0 belongs to the H 1 (Ω d ) and satisfies the compatibility conditions, i.e. u 0 = 0 on Γ Dir d and the initial velocity v 0 belongs to L 2 (Ω d ).

It is convenient to introduce the following notations:

V def = {u ∈ H 1 (Ω d ) : u = 0 a.e. on Γ Dir d }, H def = L 2 (Ω d ) and the convex set K def = {u ∈ L 2 (0, T ; V) : u t ∈ L 2 (0, T ; H), u(•, t) ∈ V, ν T u ≤ 0 a.e
         find u : [0, T ] → K such that ∫ Ω ρu tt •(v-u) dx + a(u, v-u) ≥ ∫ Ω f •(v-u) dx, for all v ∈ K, where a(u, v) def = ∫ Ω d A ijkl ε ij (u)ε kl (v) dx.
Observe that existence and uniqueness of the solution to (4.1)-(4.3) is still an open question. However, we have exhibited an approximated solution associated to (4.4) combining the finite element and mass redistribution methods. Furthermore note that the energy associated to (4.4) is given by

E(t) = 1 2 ∫ Ω ( |u t (x, t)| 2 +|∇u(x, t)| 2 ) dx - ∫ Ω f •u t dx.

The finite element approximation of problem

Introducing the Lagrange multiplier λ, then the inequality (4.4) is also formally equivalent to:

             Find u : [0, T ] → V and λ : [0, T ] → Λ ν such that for all v ∈ V ρu tt (t), v V ,V + a(u, v) = ∫ Ω f .v dx + λ ν (t), v ν X ν ,Xν a.e. t ∈ [0, T ], µ ν -λ ν (t), u ν X ν ,Xν ≥ 0 and λ ν (t) ≤ 0, for all µ ν ∈ Λ ν a.e. t ∈ [0, T ], u(•, 0) = u 0 and u t (•, 0) = v 0 , where X ν def = {ω ∈ L 2 (Γ C ) : ∃v ∈ V and ω = v ν a.e. on Γ C }, Λ ν def = {µ ν ∈ X ν : µ ν , v ν X ν ,Xν ≥ 0, for v ν ∈ V, v ν ≤ 0 a.e. on Γ C } and (u ν = u•ν. Here •, • V ,V , •, • X ν ,Xν
are the duality pairing between V and V, X ν and X ν , respectively.

The present note is organized as follows. A space semi-discretization based on a redistribution of mass is introduced in Section 4.2.1. Then two time-integration methods, namely the Crank-Nicolson and backward Euler methods are presented in Section 4.2.2. In Section 4.3, the comparison between the convergence rates for the solution and energy evolution obtained with and without the mass redistribution highlighted the efficiency of mass redistribution method.

The finite element approximation of problem

The semi-discretization in space

We introduce the semi-discrete problem in space associated to (4.4) by using a Lagrange finite element method defined on Ω d . To this aim, we introduce the following notations: Let T h be a regular mesh of

Ω d , V h def = {v h ∈ C 0 (Ω d ) : v h is piecewise linear over each K l ∈ T h \ Γ Dir d , v h = 0 on Γ Dir d }
, n d and n c denote the number of degree of freedom and the number of nodes on Γ Sig d , respectively, and a i , i = 1, . . . , n, denotes the finite element nodes. The basis of V h is defined using the set of shape functions ϕ i ∈ V h , i = 1, . . . , n d . Then the vector of degree of freedom of the finite element field u h (x, t) denoted by

U h = (u 1 (t), . . . , u n d (t)) T such that u h (x, t) = ∑ n d j=1 u j (t)ϕ j (x). Let M ij = ρ ∫ Ω ϕ i •ϕ j dx, S ij = a(ϕ i , ϕ j ) and F i = ∫ Ω f •ϕ i
dx denote the components of the mass and stiffness matrices and the external forces, i, j = 1, . . . , n d , respectively. Furthermore let B def = (ν 1 , . . . , ν nc ) T and λ def = (λ 1 , . . . , λ nc ) T . According to these notations, a finite element semi-discretization of (4.4) with nodal approximation of the contact condition is given by (4.5)

           find U h : [0, T ] → R n d and λ : [0, T ] → R nc such that MU h,tt + SU h = F + B T λ for a.e. t ∈ [0, T ], 0 ≥ ν T i U h ⊥ λ i ≤ 0 for all i ∈ I c and for a.e. t ∈ [0, T ], U h (0) = U 0 h and U h,t (0) = V 0 h , where I c def = {i : a i ∈ Γ Sig d } and such that (ν T u h )(a i , t) = (ν T i U h )(t)
for all i ∈ I c with ν T i ν j = δ i,j and ν i = 1, • being the Euclidean norm. Note that the Lagrange multipliers are indeed the nodal contact equivalent forces. The discrete energy associated to problem (4.5) is given by

(4.6) E h (t) = 1 2 (U T h,t MU h,t +U T h SU h -U T h,t F)(t).
Chapter 4. Numerical study of convergence of the mass redistribution method for elastodynamic contact problems

The multiplicity of the solution to Problem (4.5) allow us to conclude that the considered problem is ill-posed, the reader is referred to [START_REF] Moreau | Liaisons unilatérales sans frottement et chocs inélastiques[END_REF][START_REF] Paoli | Time discretization of vibro-impact[END_REF][START_REF] Khenous | Mass redistribution method for finite element contact problems in elastodynamics[END_REF]) for further details. An alternative approach to the standard discretization presented above is to consider the mass redistributions method which consists to replace the mass matrix M by a modified mass matrix M mod defined by

M mod ij def = ρ ∫ Ω mod d,h ϕ i •ϕ j dx with Ω mod d,h = {K l : K l ∈ T h , K l ∩ Γ Sig d =
(t) = U N h (t) + U N ⊥ h (t)
where N ⊥ is the orthogonal complement of N , and replacing the mass matrix M by the modified one M mod in (4.5), we get (4.7)

           find U N h : [0, T ] → N , U N ⊥ h : [0, T ] → N ⊥ and λ : [0, T ] → R nc such that M mod U N ⊥ h,tt + S(U N h + U N ⊥ h ) = F + B T λ for a.e. t ∈ [0, T ], 0 ≥ ν T i U N h ⊥ λ i ≤ 0 for all i ∈ I c and for a.e. t ∈ [0, T ], U h (0) = U 0 h and U h,t (0) = V 0 h .
Existence and uniqueness of solution associated to (4.7) as well as the energy balance have been established in [START_REF] Dabaghi | Convergence of mass redistribution method for the one-dimensional wave equation with a unilateral constraint at the boundary[END_REF][START_REF] Khenous | Mass redistribution method for finite element contact problems in elastodynamics[END_REF]. Furthermore, the convergence of the mass redistribution method is established in the one-dimensional case in [START_REF] Dabaghi | Convergence of mass redistribution method for the one-dimensional wave equation with a unilateral constraint at the boundary[END_REF] confirming the observations already made in [START_REF] Dabaghi | Numerical approximations of elastodynamic contact problems based on mass redistribution method[END_REF]. Note that the convergence of the mass redistribution is still an open problem in the higher dimension space.

The time integration methods

We introduce now the time discretization. In order to fix the notations, let the time interval [0, T ] be divided by n + 1 discrete time-points such that 0 = t 0 < t 1 < . . . < t n = T . Furthermore the discrete quantities U n h , U n h,t , U n h,tt and λ n are assumed to be given by algorithmic approximations of the displacement U h (t n ), the velocity U h,t (t n ), the acceleration U h,tt (t n ) and the Lagrange multiplier λ(t n ), respectively. Some time-stepping schemes allowing to obtain an approximated solution to Problem (4.5) are introduced below and their efficiency is discussed and analyzed in the next section.

One of the most popular method in the community of computational mechanics is the Crank-Nicolson method. This method is second-order consistent and unconditionally stable in the unconstrained case. Furthermore, the total energy of the discrete evolution is preserved, for the purely elastic case, see [START_REF] Laursen | Computational contact and impact mechanics. Fundamentals of modeling interfacial phenomena in nonlinear finite element analysis[END_REF]. However, the situation is quite different in the contact constraints case, the order of accuracy is degraded, for further details, the reader is referred to [START_REF] Hughes | Linear static and dynamic finite element analysis[END_REF][START_REF] Krenk | Energy conservation in Newmark based time integration algorithms[END_REF][START_REF] Grosu | Stability of semidiscrete formulations for elastodynamics at small time steps[END_REF]. Notice that the discrete evolution associated to (4.5) can be described by the finite difference equations:

(4.8)                      find U n+1 h : [0, T ] → R n d and λ n+1 : [0, T ] → R nc such that U n+1 h = U n h + ∆tU n h,t + ∆t 2 4 (U n h,tt +U n+1 h,tt ), U n+1 h,t = U n h,t + ∆t 2 (U n h,tt +U n+1 h,tt ), MU n+1 h,tt + SU n+1 h = F n+1 + B T λ n+1 , 0 ≥ ν T i U n+1 h ⊥ λ n+1 i ≤ 0 for all i ∈ I c ,

Numerical experiments

where ∆t is a given times step and U 0 h , U 0 h,t and λ 0 are given. Therefore we observe that (4.8) leads to the following algorithm: (4.9)

       find U n+1 h : [0, T ] → R n d and λ n+1 : [0, T ] → R nc such that ( 4M ∆t 2 +S ) U n+1 h = 4M ∆t 2 (U n h +∆tU n+1 h,t ) + MU n h,tt + F n+1 + B T λ n+1 , 0 ≥ ν T i U n+1 h ⊥ λ n+1 i ≤ 0 for all i ∈ I c ,
The energy evolution is defined by

∆E n h def = E n+1 h -E n h
, where E n h is supposed to be given by an algorithmic approximation of the energy E h (t n ). Thus the energy evolution associated to (4.8) by using the Crank-Nicolson method is

∆E n h = -1 2 (U n+1 h -U n h ) T B T (λ n +λ n+1
). Another time integration method to approach the semi-discrete problem (4.5) consists to use the backward Euler method which gives us the following discrete evolution associated to (4.5):

(4.10)                  find U n+1 h : [0, T ] → R n d and λ n+1 : [0, T ] → R nc such that U n+1 h = U n h + ∆tU n+1 h,t , U n+1 h,t = U n h,t + ∆tU n+1 h,tt , MU n+1 h,tt + SU n+1 h = F n+1 + B T λ n+1 , 0 ≥ ν T i U n+1 h ⊥ λ n+1 i ≤ 0 for all i ∈ I c ,
where U 0 h , U 0 h,t and λ 0 are given. Therefore notice that (4.10) enables us to infer the following algorithm:

(4.11)        find U n+1 h : [0, T ] → R n d and λ n+1 : [0, T ] → R nc such that ( 4M ∆t 2 +S ) U n+1 h = 4M ∆t 2 (U n h +∆tU n+1 h,t ) + F n+1 + B T λ n+1 , 0 ≥ ν T i U n+1 h ⊥ λ n+1 i ≤ 0 for all i ∈ I c ,
In this case, the energy evolution associated to (4.10) reads as

∆E n h = -1 2 (U n+1 h,t -U n h,t ) T M (U n+1 h,t -U n h,t )-1 2 (U n+1 h -U n h ) T S(U n+1 h -U n h )-∆t(U n+1 h -U n h ) T B T λ n+1 .
The reader is referred to [START_REF] Dabaghi | Numerical approximations of elastodynamic contact problems based on mass redistribution method[END_REF]Appendix] where the energy evolutions associated to the Crank-Nicolson and backward Euler methods are justified. Note that U h,tt does not appear on the right hand side of (4.11) which is not the case in (4.9) implying that the backward Euler method is dissipative and stable.

Numerical experiments

The parameters used in the numerical simulations are

Ω 2 = [0, 0.1] × [0, 1], Ω 3 = [0, 0.1] × [0, 0.1] × [0, 1], ρ = 1, the Lamé parameters λ = 0.25 and µ = 0.5, u(x 1 , 1, t) = 0.2, x 1 ∈ [0, 0.1]
for dimension two and u(x 1 , x 2 , 1, t) = 0.2, (x 1 , x 2 ) ∈ [0, 0.1] × [0, 0.1] for dimension three. The bar is undeformed at t = 0 ant it is located at the distance of 0.2 from a rigid obstacle. The bar starts to move toward the rigid obstacle with an initial velocity equal to -0.5 (vertical direction) and without any external forces. Note that we used the square and cubic meshes with square and cube elements for two and three dimensional case, respectively.

Chapter 4. Numerical study of convergence of the mass redistribution method for elastodynamic contact problems

The Crank-Nicolson and backward Euler methods are used below to compute (U n h , λ n ). These computations are performed for the standard mass matrix M as well as for the modified mass matrix M = M mod . Contrary to the one-dimensional elastodynamic contact problem treated in [START_REF] Dabaghi | Numerical approximations of elastodynamic contact problems based on mass redistribution method[END_REF], we were not able here to exhibit a true solution associated to (4.1)-(4.3). However the reference solution denoted by (U, λ) can be approximated by (U n h , λ n ) by taking ∆t and ∆x very small. Therefore (U, λ) plays the role of an explicit solution in the numerical simulations presented below. Similarly the energy E is assumed to be equal to E n h for ∆t and ∆x chosen sufficiently small. For the 2D convergence curves, the following parameters are used: ∆x 1 = ∆x 2 = 1 640 and ∆t = 78125 × 10 -4 , to get the reference solution (U, λ) and the energy E. Therefore the approximated solution (U n h , λ n ) and energy E n h are obtained starting by the parameters ∆x 1 = ∆x 2 = 1 20 and ∆t = 25 × 10 -2 which are successively divided by 2. Due to the rapidly growing number of degree of freedom in dimension three, we were only able to compute a reference solution (U, λ) and the energy E with the parameters ∆x 1 = ∆x 2 = 1 160 , ∆x 3 = 1 32 and ∆t = 625 × 10 -3 . The convergence curves have only three points where the first point is obtained with the parameters ∆x 1 = ∆x 2 = 0.1, ∆x 3 = 0.5 and ∆t = 0.1 and the others are obtained by the parameters divided by 2. We may remark that for the Lagrange multipliers, the convergence does not take place because the parameters are not small enough. The numerical simulations were performed by employing the finite element library Getfem++ (see [START_REF] Renard | Getfem++. An Open Source generic C++ library for finite element methods[END_REF]).

The Crank-Nicolson scheme

We focus first on the two-dimensional elastodynamic contact problem for which an approximated solution at the contact nodes as well as some convergence curves were obtained. More precisely, some spurious oscillations can be observed for the normal displacement at the contact nodes and for the Lagrange multiplier, when using the standard mass matrix after the contact takes place, see Figure 4.1. These oscillations do not exist in the case where the standard mass matrix is replaced by the modified one, see Figure 4.1. Note that the scaling in space is different for the Lagrange multiplier evaluated by using standard and modified mass matrices. The scheme with the standard mass matrix is unstable with a rapidly growing energy while the one with the modified mass matrix is almost conservative as the space ∆x and time ∆t steps become small. The error curves for (U n h , λ n ) in Figure 4.3 highlight the convergence of both the displacement and the Lagrange multipliers in L 2 (0, T ; V) and L 2 (0, T ; Γ C ) respectively, which is not the case for the standard mass matrix. while the displacement in the norm L p (0, T ; H), p = 2, ∞, converges in the both cases. The same condition can be drawn for the energy E n h in L p (0, T ) norm for p = 2, +∞. The convergence curves in 3D case are represented on the Figure 4.4. the numerically convergence are very similar to the two-dimensional case except the fact that there is no clear convergence of the displacement in the L 2 (0, T ; V) norm with the modified mass matrix. The evolution of the von Mises stress in a three-dimensional bar by using the standard and modified mass matrices is represented in Figures 4.5 and 4.6. As in the two-dimensional case, some spurious oscillations can be observed in a neighborhood of the contact nodes in the case where the standard matrix is considered. Then these oscillations propagated along the bar (see Figure 4.5). In the case where the standard mass matrix is replaced by the modified one, Chapter 4. Numerical study of convergence of the mass redistribution method for elastodynamic contact problems

The backward Euler method

The normal displacement at the contact nodes and the Lagrange multiplier obtained by using the standard mass as well as modified mass matrices do not oscillate (see Figure 4.7) as it can be observed for the Crank-Nicolson method. This is due to the unconditional stability and dissipativity of Euler method. Furthermore, the error curves converge for both methods (see Figures 4.9

Some complements

The straightforward application of Newmark method to contact Problem (4.5) can be written as follows:

(4.12)

                   Find U n+1 h : [0, T ] → R
+ F n+1 + B T λ n+1 , 0 ≥ ν T i U n+1 h ⊥ λ n+1 i ≤ 0 for all i ∈ I c .
This method is unconditionally stable for linear elastodynamic problem for γ ≥ 1 2 and β ≥ 1 4 ( 1 2 +γ) 2 (see [START_REF] Hughes | Linear static and dynamic finite element analysis[END_REF][START_REF] Krenk | Energy conservation in Newmark based time integration algorithms[END_REF]). The energy evolution associated to (4.12) for F independent of time is given by (4.13)

∆E n h = (1-2γ)(V n h ) T SV n h + ∆t ( β- γ 2 ) (V n h,t ) T SV n h + ∆t ( β- γ 2 ) (V n h,t ) T (B T Λ n ) -∆t(V n h ) T (B T Λ n γ ),
where

V n h def = U n+1 h -U n h , Λ n def = λ n+1 -λ n and Λ n γ def
= (1-γ)λ n + γλ n+1 (see [START_REF] Dabaghi | Numerical approximations of elastodynamic contact problems based on mass redistribution method[END_REF]).

The Newmark method I

We consider now the Newmark method with the couple (β, γ) = ( 1 2 , 1). This method is dissipative and consequently is stable, since from (4.13), we have

∆E n h = -(V n h ) T SV n h -∆t(V n h ) T (B T Λ n 1 ) < 0.
We observe that the normal displacement at the contact node does not oscillate by using the both standard and modified mass matrices while there are some oscillations for the Lagrange multiplier, when the standard matrix is used. These oscillations are reduced in the case of the modified matrix (see Figure 4.11). Note that the energies for both matrices are decreasing (see Figure 4.12). The convergence curves in 2D are represented on the Figure 4.13. We observe that the Lagrange multiplier in the L 2 (0, T ; H) norm does not converge when the standard matrix is used. In the other cases, we have the convergence with almost the same rates. In 3D, there is no clear convergence in L 2 (0, T ; V) norm (see

The Newmark method II

We deal now with the Newmark method with the couple (β, γ) = ( 12 , 1 2 ). We observe some oscillations of the Lagrange multiplier in Figure 4.15 and also a small instability of the energy evolution in Figure 4.16. When the standard matrix is replaced by the modified matrix, these oscillations of Lagrange multiplier are reduced. Moreover the energy is almost conserved (see Figure 4.16). The convergence curves in dimension two illustrate that the Lagrange multiplier does not converge in L 2 (0, T ; H) norm when the standard mass matrix is used (see Figure 4.17). In dimension three, the displacement in the L 2 (0, T ; V) norm does not converge when the standard matrix is used (see Figure 4.18). Chapter 4. Numerical study of convergence of the mass redistribution method for elastodynamic contact problems

The convergence rates

We summarize here the convergence rates of different time integration methods in the two and three-dimensional cases. In two dimensional case, we observe that by using the second order methods, namely the Crank-Nicolson and the Newmark method II, the convergence rates of displacement in L p (0, T ; H) for p = 2, ∞ are still 1 2 for the standard mass matrix. These convergence rates are improved for the modified mass matrix. Moreover the convergence of the Lagrange multipliers by using the second order methods for standard mass matrix do not occur. While for the modified mass matrix, they converge with the rate approximately equal to 1 2 . The use of the stable one order methods, namely the backward Euler and the Newmark method I, gives lower rates of convergence. Furthermore for these methods, there are no substantial differences between the convergence rates obtained by using the standard as well as modified mass matrices. Concerning the three-dimensional case, since the reference solution is computed numerically, therefore the convergence curves are obtained by limited space and time steps. However theses convergence curves confirm the efficiency of the mass modified matrix. We may conclude that the second order time integration methods are less stable than the one order methods presented. Indeed, some spurious oscillations and unfavorable behavior for the energy of system are highlighted. The mass redistribution method reduces these oscillations for the displacement and the contact force. Furthermore, better behaviors for the energy are obtained and the convergence rates are improved for the the displacement and the contact force. We may conclude that 4.4. Some complements the second order time integration methods together with the mass redistribution method is a numerically robust method for approximating the elastodynamic contact problems.

Appendix F

As a complement, let us consider here the case when the bar is dropped toward the rigid obstacle (here Γ D = ∅). The numerical three-dimensional simulations are performed by using the standard and the modified mass matrices and with the Crank-Nicolson scheme. The parameters used here are Ω 3 = [0, 0.1] × [0, 0.1] × [0, 1], ρ = 1, the Lamé parameters λ = 0.25 and µ = 0.5. We supposed that the Neumann conditions vanishes. The bar is undeformed at t = 0 ant it is located at the distance of 0.2 from a rigid obstacle. The bar starts to fall down to the rigid obstacle with an initial velocity equal to -0.5 (vertical direction). Note that any external forces are applied. Here the cubic elements with ∆x = 0.02 are used in numerical experiments. With the standard mass matrix we may observe clearly some oscillations on the body of the bar after the impact takes place. Whereas with the modified mass matrix, the oscillations vanish (see Figures 4.19 and 4.19). Note that the convergence rates for this case validate the convergence of the mass redistribution method. Etude de la convergence des méthodes de redistribution de masse pour les problèmes de contact en élastodynamique

Résumé: Le chapitre 1 porte sur une équation des ondes monodimensionnelle soumise à une condition aux limites unilatérale. Sous des hypothèses de régularité appropriées sur les données initiales, une nouvelle preuve d'existence et d'unicité est proposée. La méthode de redistribution de masse qui repose sur une redistribution de la masse d'un corps de telle sorte qu'il n'y ait pas d'inertie au niveau du noeud de contact est introduite et sa convergence est prouvée. Une approximation de ce problème d'évolution combinant la méthode des éléments finis ainsi que la méthode de redistribution de masse est analysée dans le chapitre 2. Puis deux problèmes ainsi que leurs solutions analytiques respectives (l'une étant nouvelle) sont présentés et des discrétisations possibles en utilisant différentes méthodes d'intégration en temps sont décrites. Enfin, des simulations numériques de ces problèmes sont reportées. Dans le chapitre 3, la masse des noeuds de contact est redistribuée sur les autres noeuds et sa convergence ainsi qu'une estimation de l'erreur en temps sont établies. Ensuite, une solution analytique déjà introduite dans le chapitre 3 est comparée aux approximations obtenues pour plusieurs redistributions de masse possibles mettant ainsi en évidence que plus une redistribution de masse d'un corps se fait à proximité des noeuds de contact meilleures sont les solutions approchées obtenues. Les problèmes de contact élastodynamique en dimension d'espace deux et trois sont étudiés dans le chapitre 4. Comme pour les problèmes de contact monodimensionnels, une solution approchée combinant les éléments finis et la redistribution de masse est exposée. Quelques simulations numériques utilisant des méthodes d'intégration en temps mettent en évidence les propriétés de convergence de la méthode de redistribution de masse.

Mots clés: Elastodynamique, existence, unicité, méthode de redistribution de masse, inégalité variationnelle, contact unilatéral, conservation de l'énergie, solution numérique.

Study of the convergence of the mass redistribution method for the elastodynamic contact problems Abstract: The chapter 1 focuses on a one-dimensional wave equation being subjected to a unilateral boundary condition. Under appropriate regularity assumptions on the initial data, a new proof of existence and uniqueness results is proposed. The mass redistribution method based on a redistribution of the body mass such that there is no inertia at the contact node is introduced and its convergence is proved. An approximation of this evolutionary problem combining the finite element method as well as the mass redistribution method is analyzed in chapter 2. Then two benchmark problems (one being new) with their analytical solutions are presented and some possible discretizations using different time-integration schemes are described. Finally, numerical experiments for these benchmark problems are reported. In chapter 3, the mass of the contact nodes is redistributed on the other nodes and its convergence as well as an error estimate in time are established. Then an analytical solution already introduced in chapter 3 is compared to approximate ones obtained for different choices of mass redistribution highlighting that more a mass redistribution of the body is done near the contact nodes better the approximate solutions are obtained. The two and three-dimensional elastodynamic contact problems are studied in chapter 4.

As for the one-dimensional contact problems, an approximated solution combining the finite element and mass redistribution methods is exhibited. Some numerical experiments using time-integration methods highlighted the convergence properties of the mass redistribution method.

Keywords: Elastodynamics, existence, uniqueness, mass redistribution method, variational inequality, unilateral contact, energy conservation, numerical solution.

Figure 2 :

 2 Figure 2: simulation d'accident de voiture par le groupe de "FEA-Opt Technology" (http://www.fea-optimization.com)

Figure 3 :

 3 Figure 3: Une statue étant transporté sur un traîneau, Ninive, environ 700 avant J.-C., British Museum

Figure 4 :

 4 Figure 4: Les expériences de Coulomb

  considère un domaine borné Ω d ⊂ R d avec d ≤ 3 ainsi que sa frontière ∂Ω d Sig d représentent respectivement les frontières correspondantes aux conditions de Dirichlet, Neumann et au contact unilatéral. On définit le déplacement u(x, t) pour (x, t) ∈ Ω d × [0, T ]. Le problème mathématique qui nous intéresse s'écrit de la manière suivante: (2) ρu tt -div σ(u) = dans Ω d × (0, T ), où ρ est la densité du matériau considéré et σ(u) est le tenseur de contraintes vérifiant σ ij (u) def = A ijkl ε kl (u) dans Ω d , où A ijkl est le tenseur d'élasticité et le tenseur des déformations linéarisé défini par

Figure 5 :

 5 Figure 5: Mouvement d'une barre élastique

Figure 6 :

 6 Figure 6: Les quatre domaines permettant de déterminer u

def = 0

 0 ailleurs. Il en découle que la formulation faible (8) est approchée par[START_REF] Signorini | Sopra alcune questioni di statica dei sistemi continui[END_REF]

Figure 8 :

 8 Figure 8: Comparaison de la solution analytique (u, λ) et des solutions approchées (U n h , λ n ) obtenues en utilisant les matrices de masse standard (à gauche) et modifiée (à droite) au noeud de contact.

Figure 9 :

 9 Figure 9: Comparaison des courbes d'erreurs obtenues avec les matrices de masse standard et modifiées.

4 ,Figure 10 :

 410 Figure 10: Comparaison de la solution analytique (u, λ) et des solutions approchées (U n h , λ n ) par des matrices de masse modifiées dans le noeud de contact.

Figure 11 :

 11 Figure 11: Solutions approchées obtenues en utilisant la méthode de Crank-Nicolson avec des matrices de masse standard (gauche) et modifiées (droite) aux noeuds de contact.

Figure 12 :

 12 Figure 12: Energie associée à la solution approchée pour les matrices de masse standard (gauche) et modifiées (droite) avec le schéma de Crank-Nicolson.

Figure 13 :

 13 Figure 13: Comparaison des courbes de convergence obtenues en utilisant les matrices de masse standard et modifiée avec des éléments carrés et ∆x ∆t = 2 pour le cas 2D.

(

Figure 14 :

 14 Figure 14: Comparaison des courbes de convergence obtenues en utilisant les matrices de masse standard et modifiée avec des éléments cubiques et ∆x ∆t = 2 pour le cas 3D.

L x = 0 Figure 1 . 1 :

 011 Figure 1.1: An elastic bar vibrating on impacting obstacle.

for all η

  belonging to [-L, 0]. It follows that |q (η)| 2 = |p (-η)| 2 almost everywhere on [-L, 0]. Hence we have (1.15)

Figure 1 . 2 :

 12 Figure 1.2: Four regions allowing to determine the value of u.

  and 1.5 (left)) while these oscillations do not exist in the case where M is a modified matrix (see Figures 1.4 and 1.5 (right)). Furthermore, we can observe in Figure1.6 (left), the energy is increasing with the standard mass matrix while with modified mass matrix (right), it is almost conservative.

Figure 1 . 4 :

 14 Figure 1.4: Numerical experiments with standard mass matrix (left) and with modified mass matrix (right).

Figure 1 . 5 :

 15 Figure 1.5: Numerical experiments with standard mass matrix (left) and with modified mass matrix (right) in the contact point x = 0.

Figure 1 . 7 :

 17 Figure 1.7: Numerical convergence of the solution associated to problem with mass redistribution method to the exact solution in the contact point x = 0 (left). Numerical convergence of the energy evolution associated to mass redistribution method to the exact energy (right).

Figure 1 . 8 :

 18 Figure 1.8: Initial and boundary conditions for v on the rectangle (0, x 0 ) × (0, T ).

Figure 1 . 9 :

 19 Figure 1.9: On characteristics in the region I.

  0

Proposition 1 .

 1 [START_REF] Deimling | of de Gruyter Series in Nonlinear Analysis and Applications[END_REF] p. 36 Prop. 4.2] Let X be a real Hilbert space, D ⊂ X with D = ∅, and F : D → P(X)\{∅} be monotone. Then (a) F is locally bounded on D. (b) There exists a G δ -set D 0 ⊂ D, dense in D, such that F | D is single-valued.

0 u

 0 a multivalued map F : J × D → P(X)\{∅} and x 0 ∈ D, we are looking for absolutely continuous solution of (1.49){ u(t) ∈ F (t, u(t)) a.e. on J, u(0) = u 0 ∈ D.1.5. Numerical examplesi.e. u : J → D such that u ∈ L 1 (J, R n ) and u(t) = u 0 + ∫ t (s) ds on J and u ∈ F (t, u(t)) a.e. on J.

Lemma 9 .

 9 [START_REF] Deimling | of de Gruyter Series in Nonlinear Analysis and Applications[END_REF] p. 53 Lem. 5.1] Let X = R n , D ⊂ X be closed and F :D → P(R n )\{∅} satisfy (a) F is usc, F (x) is closed convex for all x ∈ D, (b) F (x) ≤ c(t)(1 + x X ) on D, for some c > 0, (c) F (x) ∩ T D (x) = ∅ on D.Then (1.49) has a solution on R + for every x 0 ∈ D.And the result for non-autonomous problems is:Theorem 4 ([104, p. 58 Cor. 5.2]). Let X = R n , J = [0, T ] ⊂ R, D ⊂ X be closed, and F : J × D → P(X)\{∅} have closed convex values and be such that

Corollary 1 .

 1 [START_REF] Deimling | of de Gruyter Series in Nonlinear Analysis and Applications[END_REF] p. 59 Cor. 5.4] Let X = R n , J = [0, T ] ⊂ R, D ⊂ X be closed, and F : J × D → P(X)\{∅} have closed convex values and be such that Chapter 1. Convergence of mass redistribution method for the one-dimensional wave equation with a unilateral constraint at the boundary

  S * (t)u * for all u ∈ X and u * ∈ X * , 1.5. Numerical examples then S * (t) is a strongly continuous semi-group, whose generator A * is given byD(A * ) = {u * ∈ X * : u → Au, u * is continuous on D(A)in the norm topology of X}, and u, A * u * = Au, u * for all u ∈ D(A) and u * ∈ D(A * ). Definition 9. [105] A function u ∈ C 0 ([0, T ]; X) is a weak solution of (1.52) if and only if for every u * ∈ D(A * ) the function u(t), u * is absolutely continuous on [0, T ] and

Chapter 2 . 2 def=

 22 Numerical approximations of a one dimensional elastodynamic contact problem based on mass redistribution methodwhere H {u ∈ L 2 (0, T ; V) :u t ∈ L 2 (0, T ; H)}, K def = {u ∈ V : u(0) ≥ 0} and V = V Dir or V = V Neudepending on the conditions taken in x = L. Therefore the weak formulations associated to (2.1)-(2.3) and to (2.1)-(2.2) together with (2.4) are obtained by multiplying (2.1) by v -u and by integrating formally this result over Q T def = (0, L) × (0, T ) to get(2.5)

LFigure 2 . 1 :

 21 Figure 2.1: An elastic bar with both ends free to move

Figure 2 . 2 :

 22 Figure 2.2: First analytical solution: comparison of the analytical solution (u, λ) and approximated solutions (u n h , λ n ) by using standard mass matrix (left) and modified mass matrix (right) in the contact.

Figure 2 . 3 :Figure 2 . 4 :

 2324 Figure 2.3: First analytical solution: comparison of the error curves obtained by using standard and modified mass matrices. A constant ratio ∆x ∆t = 1 has been used

Figure 2 . 5 :

 25 Figure 2.5: An elastic bar with one end fixed and the other end free to move

  Before the impact, according to characteristics lines ξ def = x+t and η def = x-t, we divide the domain (0, 1) × (0, 1) into four regions as it is drawn on Figure 2.6. Then we observe that ξ ∈ [0, 1] and η ∈ [0, 1] in the region I, ξ ∈ [0, 1] and η ∈ [-1, 0] in the region II, ξ ∈ [1, 2] and η ∈ [0, 1] in the region III and ξ ∈ [1, 2] and η ∈ [-1, 0] in the region IV. Therefore we may conclude that the solution is given by u(x, t) = 1 2 (1-x) in the regions I, III, (2.41a) u(x, t) = 1 2 (1-t) in the regions II, IV. (2.41b) Chapter 2. Numerical approximations of a one dimensional elastodynamic contact problem based on mass redistribution method

Figure 2 . 6 :

 26 Figure 2.6: The regions allowing to determine the value of u.

Figure 2 . 7 :

 27 Figure 2.7: Comparison of the analytical solution (u, λ) and the approximated solutions (U n h , λ n ) by using the standard (left) and modified (right) mass matrices in the contact node.

Figure 2 . 8 :

 28 Figure 2.8: Comparison of the energy associated to the analytical solution and the energy associated to the approximated ones for the standard (left) and modified (right) mass matrices.

Chapter 2 .

 2 Numerical approximations of a one dimensional elastodynamic contact problem based on mass redistribution method

Chapter 2 .

 2 Figure 2.11: Comparison of the analytical solution (u, λ) and the approximated solutions (U n h , λ n ) by using the standard (left) and modified (right) mass matrices in the contact.

Figure 2 . 12 :

 212 Figure 2.12: Comparison of the energy associated to the analytical solution and the energy associated to the approximated ones for standard (left) and modified (right) mass matrices.

Figure 2

 2

Figure 2 . 13 :

 213 Figure 2.13: Comparison of the analytical solution (u, λ) and the approximated solutions (U n h , λ n ) by using the standard (left) and modified mass matrices (right) in the contact.

Figure 2 . 14 :

 214 Figure 2.14: Comparison of the energy associated to the analytical solution and the energy associated to the approximated ones for the standard (left) and modified (right) mass matrices.

Chapter 2 .

 2 Numerical approximations of a one dimensional elastodynamic contact problem based on mass redistribution method 2.5.2.4 The backward Euler method

Figure 2 . 15 :

 215 Figure 2.15: Comparison of the analytical solution (u, λ) and the approximated solutions (U n h , λ n ) by using the standard (left) and modified (right) mass matrices in the contact.

Figure 2 . 16 :

 216 Figure 2.16: Comparison of the energy associated to the analytical solution and the energy associated to the approximated ones for the standard (left) and modified (right) mass matrices.

2. 5 .

 5 The wave equation with Signorini and Dirichlet boundary conditions 2.5.2.5 The Paoli-Schatzman method I

Figure 2 . 17 :

 217 Figure 2.17: Comparison of the analytical solution (u, λ) and the approximated solutions (U n h , λ n ) by using the standard (left) and modified (right) mass matrices in the contact.

Figure 2 . 18 :

 218 Figure 2.18: Comparison of the energy associated to the analytical solution and the energy associated to the approximated ones for the standard (left) and modified (right) mass matrices.

Chapter 2 .

 2 Numerical approximations of a one dimensional elastodynamic contact problem based on mass redistribution method 2.5.2.6 The Paoli-Schatzman method II

Figure 2 . 19 :

 219 Figure 2.19: Comparison of the analytical solution (u, λ) and the approximated solutions (U n h , λ n ) by using the standard (left) and modified (right) mass matrices in the contact.

Figure 2 . 20 :

 220 Figure 2.20: Comparison of the energy associated to the analytical solution and the energy associated to the approximated ones for the standard (left) and modified (right) mass matrices.

2. 5 .

 5 The wave equation with Signorini and Dirichlet boundary conditions 2.5.2.7 The Paoli-Schatzman method III

Figure 2 . 21 :

 221 Figure 2.21: Comparison of the analytical solution (u, λ) and the approximated solutions (U n h , λ n ) by using standard (left) and modified (right) mass matrices in the contact.

Figure 2 . 22 :

 222 Figure 2.22: Comparison of the energy associated to the analytical solution and the energy associated to the approximated ones for the standard (left) and modified (right) mass matrices.

Figure 2 . 24 :

 224 Figure 2.24: Comparison of the convergence curves of the energy obtained by using standard and modified mass matrices.

Figure 2 . 26 :

 226 Figure 2.26: Comparison of the convergence curves of the energy obtained by using standard and modified mass matrices.

Figure 2 . 27 :

 227 Figure 2.27: Comparison of the convergence curves obtained by using standard and modified mass matrices.

(

Figure 2 . 30 :

 230 Figure 2.30: Comparison of the convergence curves of the energy obtained by using standard and modified mass matrices.

Figure 2 . 32 :

 232 Figure 2.32: Comparison of the convergence curves of the energy obtained by using standard and modified mass matrices.

Figure 2 . 33 :

 233 Figure 2.33: Comparison of the convergence curves obtained by using standard and modified mass matrices.

(

Contents 3 . 1 93 3. 2 95 3. 3 101 3. 3 . 1 101 3. 3 . 2 102 3. 3 . 2 . 2

 319329531013110132102322 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Finite element discretization and convergence of the mass redistributions method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . The wave equation with Signorini and Dirichlet boundary conditions . Analytical solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Comparisons between different mass redistributions for some time-space discretisations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102 3.3.2.1 The Newmark methods . . . . . . . . . . . . . . . . . . . . . . . . The backward Euler method . . . . . . . . . . . . . . . . . . . . . 104 3.3.2.3 The Paoli-Schatzman methods . . . . . . . . . . . . . . . . . . . . 105

Chapter 3 .= 0 Figure 3 . 1 :

 3031 Figure 3.1: An elastic bar vibrating on impacting obstacle.

  Thus the weak formulations associated to (3.1)-(3.3) obtained by multiplying (3.1) by v -u and by integrating formally this result over Q T def = (0, L) × (0, T ) to get (3.4)

Figure 3 . 2 :

 32 Figure 3.2: The regions allowing to determine the value of u.

( 3 .

 3 18) u(x, t) = { 1-x 2 in the regions I, III, 1-t 2 in the regions II, IV.

Figure 3 . 3 :

 33 Figure 3.3: Comparison of the analytical solution (u, λ) and the approximate solutions (U n h , λ n ) by modified mass matrices in the contact node.

Figure 3 . 4 :

 34 Figure 3.4: Comparison of the energy associated to the analytical solution and the energy associated to approximate solutions for modified mass matrices.

Figure 3 . 5 :

 35 Figure 3.5: Comparison of the analytical solution (u, λ) and the approximate solutions (U n h , λ n ) by modified mass matrices in the contact node.

Figure 3 . 6 :

 36 Figure 3.6: Comparison of the energy associated to the analytical solution and the energy associated to approximate solutions for modified mass matrices.

Figure 3 . 7 :

 37 Figure 3.7: Comparison of the analytical solution (u, λ) and the approximate solutions (U n h , λ n ) by modified mass matrices in the contact node.

Chapter 4

 4 Numerical study of convergence of the mass redistribution method for elastodynamic contact problems

Chapter 4 .

 4 Numerical study of convergence of the mass redistribution method for elastodynamic contact problems that the material of the bar is a homogeneous and isotropic. Let us define Ω d ⊂ R d , d = 2, 3, the reference configuration of the bar and its boundary by ∂Ω d Sig d denote the Dirichlet, Neumann and unilateral contact boundaries, respectively. In our case, we assume that Ω 2

  (4.2) u(•, 0) = u 0 and u t (•, 0) = v 0 on Ω d , and the boundary conditions are u = 0 on Γ Dir d × (0, T ) and σ(u)ν = 0 on Γ Neu d × (0, T ), (4.3a)

  ∅} for all i, j = 1, . . . , n d . Note that ker(M mod ) = N where N def = span{ν 1 , . . . , ν nc } denotes the space spanned by ν i for i ∈ I c . Thus employing the identity U h

4. 3 .

 3 Numerical experiments these oscillations do not exist anymore (see Figure4.6).

Figure 4 . 1 :

 41 Figure 4.1: Approximated solutions obtained by using the standard (left) and modified (right) mass matrices in the contact node with Crank-Nicolson scheme.

Figure 4 . 2 :

 42 Figure 4.2: Energy associated to the approximated solution for the standard (left) and modified (right) mass matrices with Crank-Nicolson scheme.

Figure 4 . 3 :

 43 Figure 4.3: Comparison of the convergence curves obtained by the modified and standard mass matrices with square elements and ∆x ∆t = 2 for 2D case by using Crank-Nicolson scheme.

4. 3 .Figure 4 . 4 :

 344 Figure 4.4: Comparison of the convergence curves obtained by the modified and standard mass matrices with square elements and ∆x ∆t = 2 for 3D case by using Crank-Nicolson scheme .

Figure 4 . 5 :

 45 Figure 4.5: Von Mises stress evolution of the deformed bar after the contact by using Crank-Nicolson scheme and with the standard mass matrix (∆x = 0.02 and ∆t = 0.01).

Figure 4 . 6 :

 46 Figure 4.6: Von Mises stress evolution of the deformed bar after the contact by using Crank-Nicolson scheme and with the modified mass matrix (∆x = 0.02 and ∆t = 0.01).

  and 4.10) and the energy are almost the same (see Figure 4.8).

Figure 4 . 7 :

 47 Figure 4.7: Approximated solutions obtained by using the standard (left) and modified (right) mass matrices in the contact node with the backward Euler method.

Figure 4 . 8 :

 48 Figure 4.8: Energy associated to the approximated solution for the standard (left) and modified (right) mass matrices with the backward Euler method.

4. 3 .Figure 4 . 9 :

 349 Figure 4.9: Comparison of the convergence curves obtained by the modified and standard mass matrices with square elements ∆x ∆t = 2 for 2D case by using the backward Euler method .

Figure 4 . 10 :

 410 Figure 4.10: Comparison of the convergence curves obtained by the modified and standard mass matrices with square elements and ∆x ∆t = 2 for 3D case by using the backward Euler method .

4. 4 .

 4 Some complements

Figure 4 . 14

 414

Figure 4 . 11 :

 411 Figure 4.11: Approximated solutions obtained by using the standard (left) and modified (right) mass matrices in the contact node with the Newmark method I.

Figure 4 . 12 :

 412 Figure 4.12: Energy associated to the approximated solution for the standard (left) and modified (right) mass matrices with the Newmark method I.

4. 4 .Figure 4 . 13 :

 4413 Figure 4.13: Comparison of the convergence curves obtained by the modified and standard mass matrices with square elements ∆x ∆t = 2 for 2D case by using the Newmark method I.

Figure 4 . 14 :

 414 Figure 4.14: Comparison of the convergence curves obtained by the modified and standard mass matrices with square elements and ∆x ∆t = 2 for 3D case by using the Newmark method I .

4. 4 .

 4 Some complements

Figure 4 . 15 :

 415 Figure 4.15: Approximated solutions obtained by using the standard (left) and modified (right) mass matrices in the contact node with the Newmark method II.

Figure 4 . 16 :

 416 Figure 4.16: Energy associated to the approximated solution for the standard (left) and modified (right) mass matrices with the Newmark method II.

Figure 4 . 17 :

 417 Figure 4.17: Comparison of the convergence curves obtained by the modified and standard mass matrices with square elements ∆x ∆t = 2 for 2D case by using the Newmark method II.

4. 4 .Figure 4 . 18 :

 4418 Figure 4.18: Comparison of the convergence curves obtained by the modified and standard mass matrices with square elements and ∆x ∆t = 2 for 3D case by using the Newmark method II .

Figure 4 . 19 :

 419 Figure 4.19: Von Mises stress evolution of the deformed bar after the contact by using Crank-Nicolson scheme and with the standard mass matrix (∆x = 0.02 and ∆t = 0.01).

Figure 4 . 20 :

 420 Figure 4.20: Von Mises stress evolution of the deformed bar after the contact by using Crank-Nicolson scheme and with the modified mass matrix (∆x = 0.02 and ∆t = 0.01).

  Présentation de la thèse où δ ij , Tr et I sont respectivement le symbole de Kronecker, la trace et le tenseur identité. Le problème (2) est complété par des conditions initiales données par

	(3)

  Theorem 6.[START_REF] Engel | One-parameter semigroups for linear evolution equations[END_REF] Every uniformly continuous semi-group (S(t)) t 0 on a Banach space X is of the form Chapter 1. Convergence of mass redistribution method for the one-dimensional wave equation with a unilateral constraint at the boundary Let X be a reflexive Banach space with norm • and dual X * . Let A be an unbounded operator in X with domain D(A) and let A be the generator of a strongly continuous semi-group S(t), such that there exist constants M and β with

	(i) S(t)	M e βt ,
	(ii)
	Proposition 3. For every strongly continuous semi-group S(t) t≥0 , there exist constant β ∈ R
	and M ≥ 1 such that
		S(t)	M e βt for all t ≥ 0.

S(t) = e tA for all t 0, for some bounded operator A ∈ L(X).

Theorem 7. [109] A linear operator A is the infinitesimal generator of a uniformly continuous semi-group if and only if A is a bounded linear operator.

Let H be a Hilbert space and for S ∈ L(H) denote by S * its adjoint, i.e. the unique operator satisfying (Sx | y) = (x | S * y) for all x, y ∈ H. Now we take a uniformly continuous group (S(t)) t∈R on H and denote its generator by A, i.e. S(t) = e tA for all t ∈ R.

Since S → S * is continuous on L(H), it follows that the adjoint group (S(t) * ) t∈R is again uniformly continuous and that S(t) * = e tA * for all t ∈ R.

The groups for which all operators S(t) are unitary, i.e. satisfy (S(t)) -1 = (S(t)) * for all t ∈ R, can be characterized as follows:

Proposition 2. The group (e tA ) t∈R is unitary if and only if A is skew-adjoint, i.e. A * = -A. Definition 8. (S(t)) t≥0 is a strongly continuous semi-group if the functional equation, defined in Definition 6, holds and the orbit maps

ξ x : t → ξ x (t) = S(t)x, are continuous from R + into X for every x ∈ X.

Proposition 4. The adjoint semi-group of a strongly continuous semi-group on a reflexive Banach space is again strongly continuous.

  The Crank-Nicolson method . . . . . . . . . . . . . . . . . . . . . 72 2.5.2.2 The Newmark method I . . . . . . . . . . . . . . . . . . . . . . . . 75 2.5.2.3 The Newmark method II . . . . . . . . . . . . . . . . . . . . . . . 76 2.5.2.4 The backward Euler method . . . . . . . . . . . . . . . . . . . . . 78 2.5.2.5 The Paoli-Schatzman method I . . . . . . . . . . . . . . . . . . . 79 2.5.2.6 The Paoli-Schatzman method II . . . . . . . . . . . . . . . . . . . 80 2.5.2.7 The Paoli-Schatzman method III . . . . . . . . . . . . . . . . . . 81 Appendix D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 Chapter 2. Numerical approximations of a one dimensional elastodynamic contact problem based on mass redistribution method

2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 2.2 Description of the models . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 2.3 Discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62 2.3.1 Classical finite element discretization . . . . . . . . . . . . . . . . . . . . . . 62 2.3.2 Finite element discretization with the mass redistribution method . . . . . 63 2.4 The wave equation with Signorini and Neumann boundary conditions 65 2.4.1 Analytical solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 2.4.2 Numerical simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 2.5 The wave equation with Signorini and Dirichlet boundary conditions . 69 2.5.1 Analytical solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 2.5.2 Numerical simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 2.5.2.1 Appendix E . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84

  .1 (see Figures 2.7, 2.11, 2.13, 2.15, 2.17, 2.19 and 2.21). Furthermore, the energy evolution associated to different methods mentioned above are presented in Figures 2.8, 2.12, 2.14, 2.16, 2.18, 2.20 and 2.22. Note that the numerical experiments were still done for different space and time steps with a constant ratio ∆x ∆t = 10. Some error curves for Crank-Nichoslon's method are 2.5. The wave equation with Signorini and Dirichlet boundary conditions

Table 2 .

 2

	Nicolson	STD MOD	0.48045 0.88075	0.47113 0.97113	-0.30273 0.38624	-0.31299 0.36192
	Newmark I	STD	0.51444	0.51768	0.26104	0.24000
	(β, γ) = ( 1 2 , 1)	MOD	0.49236	0.48034	0.28219	0.25943
	Newmark II	STD	0.57503	0.61856	0.34411	0.30100
	(β, γ) = ( 1 2 , 1 2 )	MOD	0.95055	1.00090	0.39734	0.33789
	Backward Euler	STD MOD	0.51804 0.50425	0.51771 0.49024	0.26784 0.27753	0.24847 0.25465
	Paoli-Schatzman I	STD	0.66064	0.66015	0.31204	0.28406
	e = 0	MOD	0.98009	1.00690	0.39570	0.41096
	Paoli-Schatzman II	STD	0.66064	0.66014	0.33914	0.30116
	e = 1 2	MOD	0.97659	1.00170	0.41159	0.42823
	Paoli-Schatzman III	STD	0.58876	0.61110	0.39286	0.35512
	e = 1	MOD	0.97658	1.00200	0.46380	0.46993

1: Convergence rates for the displacement

Table 2 .

 2 T ) 2: Convergence rates for the Lagrange multiplier and energy

	Crank-Nicolson	STD MOD	-0.13546 0.48812	Big Small	-0.61885 0.99486	-0.60737 0.99313
	Newmark I	STD	0.01366	Yes	0.50094	0.49943
	(β, γ) = ( 1 2 , 1)	MOD	0.31778	No	0.51230	0.49687
	Newmark II	STD	-0.01996	Big	1.05720	1.08500
	(β, γ) = ( 1 2 , 1 2 )	MOD	0.49128	Small	0.99887	1.00790
	Backward Euler	STD MOD	0.25559 0.30221	No No	0.50055 0.50918	0.49963 0.49803
	Paoli-Schatzman I	STD	0.28976	Yes	1.01620	1.02650
	e = 0	MOD	0.36928	Yes	0.99331	1.00380
	Paoli-Schatzman II	STD	0.30798	Yes	1.33470	1.00160
	e = 1 2	MOD	0.35044	Yes	1.00550	1.01450
	Paoli-Schatzman III	STD	0.04596	Yes	1.02640	0.99142
	e = 1	MOD	0.34364	Yes	1.01270	1.03150

  t)u t (x, t) dx dt.3.3. The wave equation with Signorini and Dirichlet boundary conditions(x, τ )| 2 +|u x (x, τ )| 2 ) dxSince u h converges weakly to u in H ∞ , u h H∞ converges to u H∞ and H ∞ → H 2 , we may conclude that u h converges strongly to u in H 2 .

	Therefore (3.11) and (3.17) enable us to infer that
		∫ L		∫ L
	lim h→0	0	(|u h,t (x, τ )| 2 w h (x)+|u h,x (x, τ )| 2 ) dx =	0	(|u t

  This chapter is published in the proceedings of WCCM-ECCM-ECFD 2014[START_REF] Dabaghi | Numerical study of convergence of the mass redistribution method for elastodynamic contact problems[END_REF]

Contents 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109 4.2 The finite element approximation of problem . . . . . . . . . . . . . . . 111 4.2.1 The semi-discretization in space . . . . . . . . . . . . . . . . . . . . . . . . 111 4.2.2 The time integration methods . . . . . . . . . . . . . . . . . . . . . . . . . . 112 4.3 Numerical experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113 4.3.1 The Crank-Nicolson scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . 114 4.3.2 The backward Euler method . . . . . . . . . . . . . . . . . . . . . . . . . . 120 4.4 Some complements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123 4.4.1 The Newmark method I . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123 4.4.2 The Newmark method II . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127 4.4.3 The convergence rates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130

  ijkl ε kl (u) where A ijkl and ε kl (u) ( ∂u k ∂x l + ∂u l ∂x k ) are the stress and strain tensors, respectively. Here the summation convention on repeated indices is used. Besides assuming isotropy of the material, the Hooke tensor is defined via Lamé constants λ and µ as A ijkl def = λδ ij δ kl + 2µδ kl δ jk where δ ij denotes the Kronecker symbol. The Cauchy initial data are

	•) ∂t and σ(u) is the stress tensor with σ ij (u)
	def = 1 2

def = A

  . on Γ Sig d for a.e. t}.

	Then the weak formulation associated to (4.1)-(4.3) is given by
	(4.4)

  n d and λ n+1 : [0, T ] → R nc such that ≤ 0 for all i ∈ I c ,where ∆t is a given times step and U 0 h , U 0 h,t and λ 0 are given. The Problem (4.12) leads to the following algorithm:

	U n+1 h U n+1 h,t = U n = U n h + ∆tU n h,t + h,t + (1-γ)∆tU n ( 1 -β 2 h,tt + γ∆tU n+1 ) ∆t 2 U n h,tt + β∆t 2 U n+1 h,tt , h,tt ,
	MU n+1 h,tt + SU n+1 h	= F n+1 + B T λ n+1
	0 ≥ ν T i U n+1 h	⊥ λ n+1 i
	         Find U n+1 β∆t 2 +S )	U n+1 h	=	M β∆t 2 (U n h +∆tU n h,t ) +	1-2β 2β	MU n h,tt

h : [0, T ] → R n d and λ n+1 : [0, T ] → R nc such that ( M

Table 4 .

 4 1: Convergence rates of displacement for the two dimensional case

	Scheme	Method	U n h -u L ∞ (0,T ;H)	U n h -u L 2 (0,T ;H)	U n h -u L 2 (0,T ;V)
	Crank-Nicolson		STD MOD	0.5732 0.9009	0.5756 0.9536	-0.4453 0.3419
	Newmark I (β, γ) = ( 1 2 , 1)		STD MOD	0.6007 0.4711	0.5672 0.5058	0.2563 0.2529
	Newmark II (β, γ) = ( 1 2 , 1 2 )		STD MOD	0.5757 0.8754	0.7893 0.8496	0.3133 0.3160
	Backward Euler		STD MOD	0.6129 0.5650	0.6007 0.5259	0.1953 0.1919
	Scheme		Method	λ h -λ L 2 (0,T ;H)	E h -E L ∞ (0,T )	E h -E L 2 (0,T )
	Crank-Nicolson		STD MOD	-0.2310 0.4416	-0.4855 1.1086	-0.4623 1.1386
	Newmark I (β, γ) = ( 1 2 , 1)		STD MOD	0.0081 0.4782	0.4710 0.4621	0.4621 0.4888
	Newmark II (β, γ) = ( 1 2 , 1 2 )	STD MOD	-0.0186 0.5766	1.1853 1.0864	1.1660 1.1700
	Backward Euler		STD MOD	0.4377 0.4379	0.4744 0.4631	0.4896 0.5005

Table 4 .

 4 2: Convergence rates of Lagrange multiplier and energy for the two dimensional case

Table 4 .

 4 3: Convergence rates of displacement for the three dimensional caseScheme Method E h -E L ∞ (0,T ) E h -E L 2 (0,T )

	Scheme	Method	U n h -u L ∞ (0,T ;H)	U n h -u L 2 (0,T ;H)	U n h -u L 2 (0,T ;V)
	Crank-Nicolson	STD MOD	0.2958 1.2600		0.3731 1.2482	-0.2951 0.4087
	Newmark I (β, γ) = ( 1 2 , 1)	STD MOD	0.2563 0.8113		0.3417 0.8563	-0.2625 0.4073
	Newmark II (β, γ) = ( 1 2 , 1 2 )	STD MOD	0.5205 1.2521		0.4212 1.2203	-0.1859 0.3851
	Backward Euler	STD MOD	0.1535 0.7356		0.1701 0.7507	-0.3097 0.3884
	Crank-Nicolson		STD MOD	-0.8684 1.1567	-1.2925 1.1843
	Newmark I (β, γ) = ( 1 2 , 1)	STD MOD	0.5369 0.7720	0.7181 0.8941
	Newmark II (β, γ) = ( 1 2 , 1 2 )	STD MOD	1.0310 1.0910	1.0670 1.1460
	Backward Euler		STD MOD	0.5079 0.5687	0.6647 0.8684

Table 4 .

 4 4: Convergence rates of energy for the three dimensional case

Cette thèse est accessible à l'adresse : http://theses.insa-lyon.fr/publication/2014ISAL0070/these.pdf © [F. Dabaghi], [2014], INSA de Lyon, tous droits réservés

Cette thèse est accessible à l'adresse : http://theses.insa-lyon.fr/publication/2014ISAL0070/these.pdf © [F. Dabaghi], [2014], INSA de Lyon, tous droits réservés Cette thèse est accessible à l'adresse : http://theses.insa-lyon.fr/publication/2014ISAL0070/these.pdf © [F. Dabaghi], [2014], INSA de Lyon, tous droits réservés

Efficiency of a weighted mass redistribution method for an elastodynamic problem with unilateral constraints

Appendix E

In this section, we are concerned with the convergence curves which are obtained for different time integration methods. These results have been summarized in the Tables 2.1 and 2.2.

Convergence curves for the Newmark method I