Renée El

Melhem Bât

Pascal Blaise

M Stéphane

ÉLECTROTECHNIQUE E E A Électronique

M Philippe Delachartre
email: philippe.delachartre@insa-lyon.fr

M Philippe
email: philippe.normand@univ-lyon1.fr

Mme Sylvie

M Hamamache Kheddouci

M Stéphane Benayoun
email: stephane.benayoun@ec-lyon.fr

Stéphanie Cauvin

M Jocelyn Bonjour

M Christian Montes

Keywords: . . . . . . . . . . . . . . . . . . . . . 1.3 Time-integrators for transient dynamics of continua . . . . . . 1.3.1 Time-integration in smooth dynamics . . . . . . . . . 1.3.2 Time-integration in non-smooth dynamics . . . . . . . 1.3.3 Time-integrators for NSCD with unified notations . explicit time-integrators, non-smooth contact dynamics, symplectic time-integrators, singular mass method, mortar methods

Cette thèse a été réalisée au sein du LaMCoS de l'INSA Lyon; et de l'équipe SIM/ER chez Michelin. Les travaux ont été financés par Michelin, à travers un contrat CIFRE (n • 2017/1555) de l'ANRT.

Je tiens en premier lieu à remercier mes encadrants académiques, David Dureisseix et Anthony Gravouil; et industriels, Gabriel Georges et Thomas Homolle. David et Anthony ont été présents au quotidien, tout en m'accordant une grande confiance dans la manière de mener les travaux. J'apprécie aussi le suivi réalisé par Gabriel et Thomas, à la fois sur les travaux académiques et sur l'implémentation réalisée au sein de MEF++. Ils ont su être présents aux moments nécessaires, tout en me laissant une liberté de choix. Je tiens tout particulièrement à remercier Gabriel, qui a su m'orienter et m'accompagner au quotidien pour l'implémentation dans MEF++. Chez Michelin, je tiens aussi à remercier Philippe. Nos échanges constructifs ont permis d'améliorer la formulation du contact déformable-déformable. Je remercie aussi

: the parts of mortar operators relating the mortar surface to the non-mortar surface 7 vectors with the stiffness of skin elements δ u : a vector with the differences of displacement for skin elements r : a vector with the normal contact impulses

The second part of this work aims then to achieve the conservation of energy by adapting the singular mass matrix to the CD-Lagrange. The formulation is firstly built in 1D, and shows a major improvement for the energy balance. The energy conservation is even reached under a condition on the release time. Then two possible extensions are explored for the 3D case. The first one is a penalty like formulation, fully compatible with large deformations and non-linear materials. It improves the energy balance but decreases the contact stability. The second one ensures stability at contact, but degrades the energy balance and is limited to small deformations.

The third part presents the CD-Lagrange scheme with a mortar formulation for deformable-deformable contact. The explicit feature is lost for the contact problem, but it stays nevertheless only linear. The formulation handles with stability and accuracy a contact with large sliding and friction. An acceleration technique is proposed for solving the contact problem, without any loss of accuracy.

Résumé étendu

Introduction

Concevoir un pneumatique est un processus complexe où interviennent de nombreuses contraintes. Les performances pneumatiques doivent être en accord avec les attentes client, tout en réalisant un compromis entre la technologie, les processus de fabrication disponibles, des coûts de production et de distribution maîtrisés, et les normes en vigueurs. Dans ce contexte, la simulation numérique est un outil puissant. Il permet d'explorer les performances du pneumatique vis à vis de ces contraintes, de manière rapide au sein du processus de conception.

L'usage le plus important de la simulation numérique est la prédiction des performances pour des scénarios d'usage courants. Par exemple, le pneu est simulé sous une charge statique, ou avec une vitesse de roulage constante. La simulation permet de prédire les contraintes et les déformations, qui sont traduites en indicateurs tels que la résistance au roulement, l'empreinte au sol, ou l'angle de dérive. Ces dernières années, la simulation numérique s'est largement imposée pour ces scénarios d'usage. Cela a permis une réduction du temps de développement et des coûts associés. Aujourd'hui les besoins évoluent vers des scénarios plus rares voir extrêmes dans le cycle de vie du pneu. Ils incluent notamment la simulation de chocs sur pneumatiques. L'objectif est alors double. Premièrement le but est de vérifier l'intégrité du pneu dans ces situations et donc de s'assurer de la sécurité. Ceci peut être le cas lors d'un franchissement à haute vitesse d'un nid de poule par exemple, ou d'un trottoir. La simulation permet alors de calculer les contraintes atteintes dans le pneumatique, et de les comparer à celle de rupture. Le deuxième objectif est de mieux comprendre les mécanismes provoquant l'endommagement. En effet, il est difficile d'accéder expérimentalement aux déformations et aux contraintes lors d'un choc. L'instrumentation est difficile, et se limite souvent à quelques endroits en surface ou à une observation par des caméras rapides.

Pour simuler un choc sur pneumatique, la première étape est de concevoir le modèle pneumatique lui-même. Le pneu est en effet composé de plusieurs parties, chacune remplissant un but précis. Pour assurer des performances maximales, chaque partie est composée d'un matériau différent et complexe.

Par exemple, les nappes placées à l'intérieur du pneu sont composées de câbles en acier ou textiles enrobées dans du caoutchouc. Hétérogènes elles permettent de régler finement les propriétés mécaniques dans chaque direction. Leur modèle numérique est complexe de part leur géométrie quasisurfacique et de leur hétérogénéité. Le caoutchouc placé sur la bande de roulement doit lui assurer des performances différentes sur sa plage d'utilisation. Il doit assurer un coefficient de friction suffisant lors du freinage ou des accélérations, tout en limitant les frottements lors des phases de roulement et en garantissant une durabilité élevée. Le modèle numérique est donc non-linéaire, et il fait appel à des lois incompressibles. En plus de ces lois matériaux complexes le modèle numérique doit aussi prendre en compte la finesse géométrique du pneumatique, par exemple pour les sculptures de la bande de roulement.

Les éléments finis sont a méthode la plus utilisée pour les simulations pneumatiques. Sur les modèles 3D, elle requiert des maillages très fins pour représenter la géométrie, et des lois matériaux complexes. Les systèmes obtenus contiennent donc de nombreux degrés de libertés, et sont non-linéaires. Le modèle est encore complexifié par la modélisation du contact, notamment lors de chocs. Plusieurs zones de contact apparaissent : des contacts rigidedéformables entre le sol et le pneu, l'obstacle et le pneu, ou encore la jante et le pneu; et parfois un contact déformable-déformable apparaît à l'intérieur du pneu lorsqu'un pli se forme. Les simulations avec de tels modèles nécessitent une parallélisation des calculs pour être effectuées dans des temps raisonnables. La scalabilité des algorithmes est donc un point d'attention majeur.

Un schéma dédié à la dynamique non-régulière

Pour simuler un choc sur pneumatique, le schéma d'intégration temporel doit être développé dans le cadre de la dynamique non-régulière. Introduit par Moreau et Jean dans [START_REF] Jean | The non-smooth contact dynamics method[END_REF][START_REF] Moreau | Numerical aspects of the sweeping process[END_REF], ce formalisme permet de traiter les événements non-réguliers que sont l'impact et le décollement. En effet en l'absence de contact, le déplacement, solution de la dynamique, est dérivable deux fois en temps partout. Cette régularité facilite l'intégration temporelle. Mais lors d'un impact, la vitesse présente un saut et l'accélération est donc nondéfinie. La dynamique non-régulière permet de formellement prendre en compte ces discontinuités grâce à une mesure définie sur la vitesse. Moreau et Jean introduisent aussi dans ce formalisme une condition de contact particulière. Celle-ci impose le contact non pas sur le déplacement, en interdisant toute pénétration au contact, mais sur la vitesse. La vitesse des noeuds rentrant en contact est annulée pour interdire toute pénétration ultérieure. Cette condition sous sa forme continue en temps est équivalente à celle sur le déplacement. Sous sa forme discrète, elle induit une légère pénétration constante durant le contact qui disparaît à convergence espace-temps.

Un schéma intéressant, développé dans ce formalisme, est le CD-Lagrange introduit dans [START_REF] Fekak | A new heterogeneous asynchronous explicit-implicit time integrator for nonsmooth dynamics[END_REF]. Ce schéma rassemble plusieurs propriétés adaptées à la simulation de chocs sur pneumatiques.

• Intégration temporelle explicite. Le CD-Lagrange utilise le schéma explicite de la différence centrée pour l'intégration temporelle. Les schémas explicites sont particulièrement adaptés aux systèmes nonlinéaires. En effet, l'intégration temporelle ne requiert aucune inversion de système. La différence centrée est de plus symplectique: pour le système discret, elle conserve le moment angulaire à sa valeur exacte, et une valeur approchée de l'énergie. Ceci assure une précision élevée en cas de grandes rotations et de simulations sur un temps long.

• Contact par multiplicateurs de Lagrange. Le schéma CD-Lagrange utilise des multiplicateurs de Lagrange pour imposer le contact. Cette méthode ajoute les forces de contact comme degrés de libertés. Pour les schémas explicites cette méthode est la plus précise.

• Condition de contact sur la vitesse. Le contact est imposé par un couple vitesse-impulsion à la façon de Moreau-Jean. Cette condition rend le contact très stable.

Une limitation majeure de la formulation actuelle du schéma CD-Lagrange est qu'il est limité à un contact rigide-déformable pour garder son aspect explicite.

Cette thèse poursuit deux objectifs principaux:

1. identifier les points d'amélioration des schémas explicites pour la dynamique avec impact, et y apporter de nouvelles solutions;

2. démontrer la faisabilité de la simulation de chocs sur pneumatiques de manière robuste, précise et efficace grâce au formalisme de la dynamique non-régulière.

Le premier objectif est académique. La maturité visée pour les nouveaux résultats proposés ne permettra pas une application dans un cadre industriel. Le second, par contre, est plus appliqué. Il vise par des techniques déjà existantes à concevoir un schéma adapté à la simulation de chocs sur pneumatiques. Ce schéma sera implémenté au sein de l'environnement de développement utilisé chez Michelin. Celui-ci est basé sur MEF++, un code éléments finis développé conjointement par le GIREF (Groupe Interdisciplinaire de Recherche en Éléments Finis) de l'université Laval de Québec, et des industriels partenaires comme Michelin ou Hydro-Québec. MEF++ est un outil de recherche en méthodes numériques, utilisable à la fois dans un contexte académique et industriel. Il se base en effet sur des librairies parallèles qui lui permettent de résoudre les problèmes de taille rencontrée dans l'industrie.

Un benchmark pour la dynamique non-régulière

Le premier apport de cette thèse est de proposer une banque de cas tests.

Ces cas permettent de tester facilement un intégrateur temporel dédié à la dynamique non-régulière. Chacun cible une propriété requise pour garantir un intégrateur adapté. Les cas se veulent simples pour une implémentation facile et ne nécessitant pas de code éléments finis.

La balle rebondissante

Le premier test est celui de la balle rebondissante, focalisé sur un contact rigide-rigide. Une balle, modélisée par un point, rebondit sur un sol rigide. Au contact est introduit un coefficient de restitution. Compris entre 0 et 1, il règle la quantité d'énergie absorbée par l'impact. Pour un coefficient égal à 1, l'impact est élastique : toute l'énergie est restituée. Ce test permet de tester très facilement la conservation de l'énergie. La solution discrète conserve l'énergie si la hauteur atteinte par la balle à chaque rebond est celle initiale. Ce test montre que le CD-Lagrange conserve l'énergie pour un impact élastique. Pour les schémas testés (issus de la dynamique non-régulière), le taux de convergence est d'ordre 1. Le même ordre de convergence est trouvé sur les cas tests suivants. La balle rebondissante permet également de tester le paradoxe de Zénon : une infinité d'impacts en un nombre fini de pas de temps, pour un coefficient compris entre 0 et 1. Ce paradoxe ne peut-être simulé que par des schémas traitant l'équilibre dynamique au sens faible en temps, comme ceux construits dans le formalisme de Moreau-Jean.

Le ressort de Van der Pol

Ce cas test est constitué d'un système masse-ressort, augmenté d'un terme d'amortissement non-linéaire sur le déplacement et proportionnel à la vitesse. Pour rendre le système non-régulier, une contrainte de contact est ajoutée sur la masse. Ce système tend vers un cycle limite dans l'espace des phases. Il permet donc de tester la conservation de l'énergie sur un temps long. Seul les intégrateurs conservant l'énergie se maintiennent sur ce cycle limite, ce qui sera le cas des intégrateurs symplectiques. Ce système démontre aussi l'avantage explicite de la différence centrée: la non-linéarité due au terme d'amortissement est résolue sans aucune complexité algorithmique.

Le ressort tournant

Le ressort tournant est un système masse-ressort linéaire soumis à de grandes rotations. Une contrainte de contact est là encore ajoutée sur la masse pour rendre le système non-régulier. Ce cas test vérifie la conservation du moment angulaire pour un système non-régulier. Le CD-Lagrange conserve exactement le moment angulaire à sa valeur initiale, là où il décroît à chaque impact pour un autre intégrateur tel que le schéma de Moreau-Jean. Ajouté à la conservation de l'énergie, la conservation du moment angulaire est caractéristique d'un intégrateur symplectique.

La barre en impact

Ce cas test introduit une discrétisation éléments finis. Il consiste en une barre 1D discrète, dont une extrémité impacte une frontière rigide. Le contact modélisé est ici de type rigide-déformable. Il comprend donc temporellement trois phases : l'impact, quand la barre rencontre la frontière rigide; le contact; puis le décollement, où la barre quitte le contact. Contrairement à un contact rigide-rigide, impact et décollement sont deux événements nonréguliers distincts. Ce cas test met en avant une limitation de la condition de contact en vitesse. À l'impact, l'énergie cinétique du noeud de contact est annulée. Cela garantie la stabilité au contact, mais empêche la conservation de l'énergie.

Le dôme en impact

Ce cas test est relativement semblable à celui de la barre en impact mais pour un maillage 3D. Un solide 3D, un rectangle surmonté par un demicylindre, impacte une frontière rigide. Le contact, cette fois-ci en deux dimensions, permet d'ajouter de la friction et du glissement. Ce cas met en avant l'efficacité du CD-Lagrange pour un contact rigide-déformable, le contact étant résolu de manière explicite. L'algorithme est alors naturellement parallèle.

Le schéma CD-Lagrange est donc prometteur pour la simulation de chocs sur pneumatiques. Néanmoins il possède deux axes d'amélioration :

• pour un contact rigide-déformable, la conservation de l'énergie lors de l'impact rendrait le schéma symplectique (conservation de l'énergie et du moment angulaire);

• il doit être étendu au contact déformable-déformable sous une formulation efficace.

La masse singulière : vers la conservation de l'énergie

Le troisième chapitre est consacré à adapter la technique de la masse singulière au schéma CD-Lagrange. Cette technique a été introduite dans [START_REF] Khenous | Mass redistribution method for finite element contact problems in elastodynamics[END_REF] et [START_REF] Hager | A stable energy-conserving approach for frictional contact problems based on quadrature formulas[END_REF]. Elle intervient sur le système semi-discret en espace, et consiste à annuler les entrées de la matrice de masse correspondant aux degrés de libertés où s'appliquent les contraintes de contact. Le système semi-discret avec contraintes de contact retrouve alors des propriétés facilitant l'intégration temporelle. En effet, le système avec masse singulière admet une unique solution conservative en énergie, là où le système avec masse consistante peut admettre une infinité de solutions. La masse singulière permet d'améliorer les performances de schémas implicites [START_REF] Dabaghi | A robust finite element redistribution approach for elastodynamic contact problems[END_REF][START_REF] Krause | Presentation and comparison of selected algorithms for dynamic contact based on the Newmark scheme[END_REF]. Instables sur le système consistant, la stabilité est retrouvée pour le système singulier. Le bilan d'énergie, et la stabilité du contact sont aussi améliorés. La technique de la masse singulière est un moyen potentiel d'atteindre la conservation de l'énergie à l'impact pour le CD-Lagrange. En effet, seul l'impact est dissipatif en énergie et l'énergie dissipée alors correspond à l'énergie cinétique des noeuds rentrant en contact. Si ces noeuds n'ont pas de masse, l'énergie dissipée devrait être nulle.

La première étape est d'intégrer la masse singulière dans une formulation 1D du CD-Lagrange. Ceci est réalisé en s'appuyant sur le cas test de la barre en impact. La difficulté principale est d'établir la vitesse des noeuds de contact, celle-ci n'étant plus déterminée par la dynamique. Une loi de contact est proposée. Elle permet d'appliquer une condition de contact de Moreau-Jean, et de calculer la vitesse des noeuds de contact. Cependant elle reste une approximation. Les résultats numériques démontrent qu'avec cette formulation le CD-Lagrange conserve un contact stable. La condition de contact est de plus respectée partout sauf au décollement où les forces de contact sont négatives. Le bilan d'énergie est amélioré : il est conservatif à l'impact, mais devient non-conservatif au décollement. Cependant si le décollement coïncide avec un temps discret, le décollement est conservatif et la condition de contact y est respectée (forces de contact nulles). Ce résultat est un résultat majeur pour la dynamique non-régulière : la conservation de l'énergie est possible sur le système discret.

Cette formulation 1D est ensuite étendue aux cas 3D. La première extension est réalisée par l'ajout d'élément normaux sur la peau du solide où s'applique le contact. Ces éléments correspondent exactement à l'élément de la barre 1D avec le DDL singulier. Cette fois l'élément normal est seulement numérique. La rigidité de ces éléments normaux est un paramètre supplémentaire permettant de régler la réponse de la peau, et d'optimiser le bilan d'énergie. Les résultats numériques montrent un bilan d'énergie amélioré, mais une dégradation de la stabilité du contact. Se rapprochant d'une formulation pénalisée, l'implémentation est facile et totalement compatible avec des grandes déformations ou des matériaux non-linéaires.

Une seconde extension 3D est ensuite proposée. Celle-ci se rapproche plus des formulations existant sur les schémas implicites. Les degrés de libertés à masse singulière sont ceux du maillage initial, et non pas ajoutés sur la peau. La formulation proposée nécessite des propriétés particulières sur le maillage pour pouvoir appliquer la loi de contact, et elle se limite aux petites déformations. Les résultats numériques montrent un contact stable, mais un bilan d'énergie dégradé. Après le décollement les forces de contact ne tendent pas vers une valeur nulle, et apportent de l'énergie au système.

Ceci démontre clairement les limites de l'approximation réalisée par la loi de contact.

Le contact déformable-déformable par les méthodes mortier

Le quatrième chapitre intègre les méthodes mortiers dans le schéma CD-Lagrange. Ces méthodes ont été introduites notamment dans [START_REF] Belgacem | The Mortar finite element method with Lagrange multipliers[END_REF] puis développées dans [START_REF] Popp | Mortar Methods for Computational Contact Mechanics and General Interface Problems[END_REF][START_REF] Puso | A mortar segment-to-segment contact method for large deformation solid mechanics[END_REF][START_REF] Wohlmuth | A Mortar Finite Element Method Using Dual Spaces for the Lagrange Multiplier[END_REF][START_REF] Yang | Two dimensional mortar contact methods for large deformation frictional sliding[END_REF]. Elles permettent de relier deux maillages nonconformes par une projection est effectuée au sens faible. Cette projection rend les méthodes mortiers très robuste pour de grands glissements au niveau du contact. Ceci les rend totalement adaptées pour du contact déformabledéformable dans des problèmes en grandes déformations.

Les méthodes mortier couplées au schéma CD-Lagrange conduisent à un problème de contact qui n'est plus explicite car nécessitant un solveur. Néanmoins grâce à la connaissance préalable de la position et la direction de glissement, le problème de contact reste seulement linéaire. Un solveur est proposé pour résoudre ce problème. Il est basé sur un algorithme de Gauss-Seidel, et donc séquentiel. Ce choix est motivé par des raisons de robustesse. Chaque itération de l'algorithme mène en effet à une approximation de plus en plus proche du résultat final. De plus si le problème de contact déformable-déformable reste petit par rapport à la taille du problème global, les performances du schéma ne devraient pas être impactées.

Pour accélérer la résolution du problème de contact, une technique de diagonalisation (lumping) est appliquée sur les opérateurs mortier comme proposé dans [START_REF] Bussetta | Efficient 3D data transfer operators based on numerical integration[END_REF]. Cette technique permet une diminution de 90% du nombre d'itérations dans la résolution du contact déformable-déformable, et ce à précision égale. Il reste néanmoins à effectuer une vraie mesure de l'efficacité du schéma sur un cas industriel parallèle.

Ces résultats pour les méthodes mortiers et le CD-Lagrange sont démontrés sur un cas test numérique, avec un contact glissant sans grands glissements. Le contact est stable et précis. Le bilan d'énergie présente une perte d'énergie lors de l'impact mais elle reste limitée. La solution obtenue avec les opérateurs mortar diagonalisés converge vers celle provenant des opérateurs standards.

Suivent trois cas tests numériques se rapprochant d'un cas pneu. Le premier modélise un contact glissant en présence de grands glissements. Le bilan d'énergie est alors dégradé par les pertes d'énergie lors des impacts. En effet ceux-ci sont nombreux sur les bords de la zone de contact à cause du glissement. Le second et troisième cas test intègrent de la friction sans puis avec grands glissements. Les résultats numériques démontrent la capacité du CD-Lagrange à modéliser ce type de contact.

Conclusion

Grâce au benchmark proposé, le CD-Lagrange démontre sa capacité à simuler des chocs sur pneumatiques. Grâce à son intégration temporelle explicite, les non-linéarités du système sont facilement traitées, et son aspect symplectique assure une grande précision. La condition de contact sur la vitesse apporte elle stabilité et précision au niveau du contact. Elle permet d'intégrer facilement de la friction. Et surtout elle conduit à un problème de contact explicite dans le cas rigide-déformable. Ceci assure une parallélisation facile du CD-Lagrange, et donc une grande efficacité. Néanmoins deux améliorations sont nécessaires et traitées par la suite. La première serait d'atteindre la conservation de l'énergie à l'impact, ce qui rendrait le schéma symplectique pour un contact déformable. La seconde est d'étendre la formulation du contact à un cas déformable-déformable.

L'adaptation de la masse singulière au schéma CD-Lagrange apporte de nouveaux résultats. La conservation de l'énergie sur le système discret est un résultat majeur, même si celui-ci ne concerne qu'une formulation 1D et est conditionné à un décollement sur un temps discret. La première formulation 3D proposée est intéressante. Elle permet une application facile de la masse singulière, tout en améliorant le bilan d'énergie. Se rapprochant d'une formulation pénalisée, la stabilité du contact est légèrement dégradée. Cette formulation est nouvelle car elle propose un schéma explicite à masse singulière, et une formulation proche de la pénalisation sur un schéma où le contact est écrit en vitesse. La seconde formulation démontre clairement les limites de la loi de contact proposée, et souligne la difficulté d'adapter la masse singulière à un schéma explicite sur la vitesse.

L'extension du CD-Lagrange à cas déformable-déformable se fait de manière adaptée aux simulations de chocs sur pneumatiques. Si l'aspect explicite du schéma est perdu, le problème de contact déformable-déformable ne requiert qu'un solveur itératif linéaire. Un algorithme séquentiel est proposé : plus robuste, il ne devrait pas impacter l'efficacité globale du schéma pour un problème restant petit. Une technique de diagonalisation des opérateurs est appliquée, démontrant une franche accélération de la résolution pour une précision égale. La simulation de chocs sur pneumatique est donc possible avec le schéma CD-Lagrange : il permet de traiter du contact déformabledéformable avec maillages incompatibles et grands glissements, et ce de manière efficace sur des modèles à vocation industrielle.

Les suites directes de ce travail se divisent en deux catégories.

Introduction

Industrial context

Numerical simulation in tire design

Designing a tire is a complex process with several constraints. It aims to line the performances with the customer requirements, while taking into account the available manufacturing processes, the cost of production and the norms of the countries. In this design process, the numerical simulation is a powerful tool and operates at many steps. It helps the tire designer to choose a design by exploring its performances. The first application is to simulate the tire in standard usage scenarios. For example, it predicts the tire deformation under a static load and for several inflate pressures, or for a constant rolling speed. The numerical results give an access to loads and deformations, which are converted into measures of tire performances: rolling resistance, contact patch area, slip angle, or even longevity prediction. In the last decades, the numerical simulation has replaced the experimental measures for these common usages scenarios. As a consequence, the timeliness and the cost of the design process were considerably reduced. Moreover the progress in numerical simulations and the cost reduction of computational power made accessible new scenarios. The integration of fluid-structure interactions allows today to numerically predict the tire performances in presence of water on the road. Or more recently, a need in noise simulation grew with the development of electrical mobility. With the electrical engines, the tires become indeed a non-negligible source of noise.

Today the tire designers want to simulate rare and extreme life scenarios like shocks. Indeed they provoke damages which decrease the tire performances or cause safety issues. A clear comprehension of these scenarios would bring then the optimization further. An example of an accident scenario could be a pothole crossed at high speed. Simulating the shock measures the loads in the tire and verifies if the breaking strength is reached. Moreover the motion accessible at every time and for every part allows a deeper comprehension of the damages. In these extreme scenarios, the role of numerical simulation is then to check the safety of the design. In more common scen-arios, the damages come from the repetition of shocks at lower speed on weaker zones. For an example in urban area delivery, the light vans frequently touch or climb the pavement for parking. The simulation quantifies precisely the loads and the motion involved, and allow an optimization in order to prevent the early damages due to the fatigue. For shocks, the numerical simulation presents a major advantage on the experimental tests besides its lower cost. The experimental tests are indeed particularly difficult. The instrumentation is tough and limits the strain measures to particular places, or through the observation of global deformation with fast cameras.

The numerical simulation plays a role too in the prediction of manufacturing feasibility. The final manufacturing quality depends on the precise execution of each manufacturing steps. The manufacturing includes mechanical steps as the conformation. During this step, the tire is formed into its ring-shape. It includes too thermic transformations especially during the curing, where the tire is heated to its final aspect. The numerical simulation of these steps allow to increase their accuracy, and then the final quality of the tire.

A last application of numerical simulation is for the normative tests. They verify the conformity of pneumatics with standards of safety defined by the authorities. A common one is the "plunger energy" test. It applies on the sidewall or on the tread a load thanks to a plunger to check the tire resistance. As this test is particularly damaging for the tire, the numerical simulation saves cost by verifying conformity of a design without a complete build.

The numerical simulation is then crucial for tire design. By saving cost and time, it accelerates and makes more flexible the design process. It is also a powerful lever of innovation by bringing a fine comprehension of physical phenomenons.

Modelling a shock on a tire

Each part of a tire has a specific purpose which requires distinct material, frame and geometry. The structure is then a complex and heterogeneous assembly which makes the numerical simulation difficult. For example, the tread is the external layer of the tire which ensures the contact with the ground. Its material must be soft enough to ensure friction during braking or acceleration, while being resistant to maximise the tire life-time and limiting the rolling resistance. As the material response changes in the operating range, it is simulated by a non-linear law which integrates incompressibility. Besides the complex material, the tread is sculpted with slots and grooves in an highly precise geometry. The tire contains also surface-shaped parts made in heterogeneous material as the belts. The carcass shape and rigidity is indeed ensured by belts made with steel or textile wires embedded in a rubber matrix. Depending on the direction of traction, the rigidity changes in these belts. Moreover there are superposed inside the tire to finely set the carcass rigidity in a direction. In addition to these difficult to simulate materials, a tire presents a large diversity of stiffness. The bead, the part which contacts the rim, presents the higher stiffness. It is made by a stiff rubber and reinforced with a steel wire. Its role is indeed to seal the tire and the rim even for the highest loads. The softest material is found in the inner liner. It does not have any mechanical purpose as it only ensures the air tightness. This stiffness diversity is a numerical challenge as it involves several orders of magnitude between the fastest and the slowest waves. The accuracy of the simulation is yet conditioned to a fine capturing of each time-scale.

The main simulation technique for tires is the finite element method. It requires a geometrical mesh, driven by the dimensions of the tire parts. For 3D models, the spatial discretization results then in fine meshes with a large number of degrees of freedom.

The numerical complexity is further increased in case of shocks. The shock area requires an high computation accuracy, as it presents the highest loads. The contact geometry can be complex, depending on the obstacle. And in the most extreme situations, a fold appears inside the tire and causes self-contact on the inner layer. These zones requires then a precise simulations of deformations and strains under geometrical and material nonlinearities. Outside the shock zone, contact happens between the tire and the ground, and between the tire and the rim. If they are easier to model, these areas include a large number of degrees of freedom. In the most of cases, a fine model of the ground is not necessary. A simple geometric description together with a macro friction coefficient are accurate enough. But new needs are emerging with finer models for the ground. Some situations requires the deformations of the ground due to the tire, in order to compute the resulting traction. These new scenarios are met on off-road tires like in mining or agriculture.

The model of tires, especially in case of shocks, contains then a large number of degrees of freedom with non-linear equations. Together with the fine time-discretization required for accuracy, the time-cost of a simulation is high. The efficiency of numerical simulations is then crucial for taking part in the design process. Two main levers influence the computation time: the efficiency of algorithms, and the available computational power. In practice both are linked by the parallelism technique. In order to increase the computational power, the simulations are run in parallel on multiple processors. The efficiency of an algorithm and more generally of a finite element software is then conditioned to its capacity to be parallelised.

Goals

This work aims to develop a time-integrator for finite elements method efficient in shocks simulation on tire. It was done in collaboration with Michelin, where the numerical simulation is currently developed toward more precise and efficient simulations of shocks.

The time-integrator requirements are:

• the robustness for dynamics models, with non-linear materials, large deformations, and high-load contact;

• the accuracy especially for the contact quantities;

• and efficiency for being compatible with parallel runs.

These requirements lead to the following starting choices.

1. An explicit time-integrator. The unknowns of the problem must be computed without any system solving. Thanks to this explicit feature, the non-linearities do not bring extra complexity. A main drawback is the conditional stability of explicit time-integrators. If the timediscretization is not fine enough, the time-integrator is unstable. But in shock simulation, the required accuracy leads already to a fine timediscretization for catching all the time-scale.

2.

A contact enforcement by Lagrange multipliers. Several methods coexist for enforcing a contact constraint on a discrete system. The Lagrange multipliers are chosen, being the most accurate especially for explicit time-integrators. Indeed in a shock simulation, the focus is on the contact zone where the higher accuracy is required. The Lagrange multipliers can impact the efficiency because they add degrees of freedom corresponding to the contact forces. This increase must stay small in order to not change the global efficiency of the algorithm.

3.

A formulation for rigid-deformable and deformable-deformable contact. In a shock simulation, the contact involves mainly a rigid and motionless body (the ground or the obstacle) and a deformable body (the tire). But a deformable-deformable contact at the inner layer is possible. If it happens, it becomes the main zone of interest as the most probable place for a damage. The scheme must then deal with both type of contact, rigid-deformable and deformable-deformable, under large deformations.

This work has two main objectives:

1. to explore the improvements of current explicit time-integrators for dynamics with impacts;

Outline

This thesis manuscript is organized in four chapters.

1. The first chapter provides an introduction into the theoretical framework. The first part introduces the solid mechanics problem, and its spatial discretization by the finite elements method. A description of contact constraints on a deformable body follows. They bring non-smooth events in the dynamical problem, which requires a new framework. The non-smooth contact dynamics framework is then introduced. It describes the non-smooth dynamics problem for rigid or deformable bodies discrete in space. This first chapter continues with a review of time-integrators for dynamics. It starts with time-integration in smooth dynamics and then focuses on non-smooth dynamics. Some non-smooth time-integrators are introduced in unified notations, in order to ease the comparison in the next chapter.

2. The second chapter presents numerical test cases which form a benchmark. Each case highlights a required property for a time-integrator designed for non-smooth dynamics. The benchmark identifies the CD-Lagrange scheme as a promising time-integrator for shock simulations. Indeed it is explicit, uses Lagrange multipliers for enforcing the contact, presents an high stability and accuracy especially at contact. Nevertheless two improvements are identified by the benchmark. Its energy balance is conservative except at impact where energy is dissipate. The first improvement in regard of accuracy is then to reach energy conservation. Moreover the CD-Lagrange scheme does not address the deformable-deformable contact with large sliding necessary in a shock simulation. The two next chapters are devoted to improve it in these ways.

3. The third chapter tries to improve the energy balance of CD-Lagrange thanks to the singular mass method. A first formulation is proposed on a 1D problem. It identifies the main obstacles to the singular mass method in explicit time-integrators. Solutions are proposed to overcome these difficulties. The formulation demonstrates a major improvement for the energy balance. Two 3D formulations follow as two different extensions of the 1D formulation.

4. The fourth chapter introduces the mortar method in the CD-Lagrange scheme. It demonstrates its ability to simulate shocks between two deformable bodies, as it handles friction and large motion at contact. The efficiency is explored, and an acceleration technique is proposed for solving the contact problem. Numerical results are presented, not on tires, but with a focus on contact.

Chapter 1

Space and time discretization of impact dynamics problem

This chapter is a global introduction to impact dynamics. The first part presents the dynamical mechanics problem. After a presentation of the general governing equations, the principle of virtual work is introduced. It forms a starting point for the spatial discretization, done by a finite element method. The spatial discretization leads to the discrete-in-space problem.

In addition, the first part briefly exposes the Lagrangian and Hamiltonian formalisms.

The second part adds contact constraints in the problem, and explores the consequences on the equations. Non-smooth events are indeed added in the dynamical problem. This new problem of impacts mechanics is described firstly in the non-smooth Lagrangian formalism. And then follows a description of the non-smooth contact dynamics formalism, more suited to build robust time-integrators.

The third part is devoted to time-integration in order to get space-andtime discrete problems. After a description of time-integration in smooth dynamics, a brief review of non-smooth time-integrators is made. Some of them are presented under unified notations to ease the comparison.

Spatial discretization for solid mechanics

Governing equations

The considered problem is described in fig. 1.1. A deformable body occupies the domain Ω 0 ⊂ R 3 in its initial configuration at time t 0 , and the domain

Ω(t) ⊂ R 3 in a deformed configuration at time t ∈ [t 0 , t f ]. [t 0 , t f ] ⊂ R is the time-interval.
The boundary of the domain in initial configuration ∂Ω 0 is split into two complementary sets Γ D and Γ N : All material points X ∈ R 3 of Ω 0 and all material points x ∈ R 3 in Ω(t) are linked by the bijective total Lagrangian deformation map:

∂Ω 0 = Γ D ∪ Γ N and Γ D ∩ Γ N = ∅ (1.1) Γ D Γ N Ω 0 Γ D (t) Γ N (t) Ω(t) X x Φ t
Φ t : Ω 0 → Ω(t) X → x (1.2) x (X , t) = Φ t (X , t) and X = Φ -1 t (x, t) (1.3)
The displacement is defined on the initial configuration as:

u(X , t) = x (X , t) -X (1.4)
The following tensors describe the deformations (on Ω 0 ):

• F = ∂x ∂X , the deformation gradient, maps an infinitesimal line element ∂X ∈ Ω 0 onto the corresponding one ∂x ∈ Ω(t);

• F -1 , the inverse deformation gradient;

• J = det(F ) > 0, the Jacobian determinant of the deformation;

• C = F t F , the right Cauchy-Green tensor, provides an objective measure of deformation being invariant under any rigid motions;

• E = 1 2 (C -I), the Green-Lagrange strain tensor, measures the strains inside the body and gives a null strain for a null deformation.

The stresses are described on Ω(t) by σ, the Cauchy stress tensor. It links the traction t on a surface to the internal stresses by t = σ • n, n being here the outer normal to the surface. On the initial configuration Ω 0 , the internal stresses are described by two main tensors:

• P = Jσ t F -t , the first Piola-Kirchhoff tensor, not symmetric;

• S = F -1 P , the second Piola-Kirchhoff tensor, symmetric.

Stresses and deformations are linked through constitutive equations or material laws. In a general way these laws are a relation between the deformation, depending on the displacement, and the stresses:

P = f (E) = f (u).
The details of such laws are not addressed here. For more informations see [START_REF] Belytschko | Nonlinear Finite Elements For Continua And Structures 2nd[END_REF]. Nevertheless an interesting case is the hyperelastic material model, in which the stress-strain law derives from a strain energy function Ψ(E). It gathers the Saint Venant-Kirchhoff model, used in its linear form for linear elasticity; but also more complex material like Neo-Hookean or incompressible hyperelastic models used for representing rubber and other materials in a tire.

The boundary Γ D is subjected to Dirichlet's boundary conditions: the displacement û is imposed for all t ∈ [0, t f ]. On Γ N , Neumann's conditions are enforced: the traction t is imposed for all t ∈ [0, t f ].

In order to write the equations describing the body, the following notations denote:

• ρ 0 and ρ the mass density for the initial or the deformed configuration;

• f v,0 and f v the density of external body forces expressed in appropriate coordinates (X or x ); On current configuration Ω, the governing equations are:

• d• dt ,
∇ x • σ + f v = ρ ü(x , t) on Ω × [t 0 , t f ] (1.5) σ = f (E) on Ω × [t 0 , t f ] (1.6) u(x , t) = û(x , t) on Γ D × [t 0 , t f ] (1.7) σ • n = t(x , t) on γ N × [t 0 , t f ] (1.8)
On initial configuration Ω 0 , the governing equations of the problem are:

∇ X • P + f v,0 = ρ 0 ü(X , t) on Ω 0 × [t 0 , t f ] (1.9) P = f (E) on Ω 0 × [t 0 , t f ] (1.10) u(X , t) = û(X , t) on Γ D × [t 0 , t f ] (1.11) σ • n = t(X , t) on γ N × [t 0 , t f ] (1.12)
(1.5) and (1.9) are the dynamics of bodies respectively expressed on deformed and initial configurations; (1.6) and (1.10) are the material laws. They are completed by the Dirichlet's conditions (1.7) and (1.11), and the Neumann's conditions (1.8) and (1.12). These sets of equations are the strong form of the problem.

Γ D Γ N Ω 0 Γ D Γ N Figure 1.2: Spatial discretization of a deformable body 1.1.

The principle of virtual work

The solutions of eqs. (1.9) to (1.12) are in the space:

U = {u ∈ H 1 (Ω) | u(X , t) = û(X , t) on Γ D } (1.13)
completed by the space:

V = {v ∈ H 1 (Ω) | v (X , t) = 0 on Γ D } (1.14) H 1 (Ω)
is the Sobolev space:

H 1 (Ω) = w on Ω, Ω w 2 dx < ∞, Ω ∂w ∂x 2 dx < ∞ (1.15)
The weak form is obtained from the dynamic (1.9) by multiplying by δu ∈ V and integrating over Ω 0 . Thanks to the Gauss divergence theorem and δu = 0 on Γ D , the result is:

∀δu ∈ V, Ω 0 ρ ü • δu dΩ 0 + Ω 0 Tr P (∇δu) t dΩ 0 - Ω 0 f v,0 • δu dΩ 0 - γ N t • δu dγ N = 0 (1.16)
The equation (1.16) is the principle of virtual work (PVW), the starting point of the spatial discretization.

Spatial discretization by the finite elements method

From the weak form to the discrete-in-space system

The domain Ω 0 is approximated by a partitioning of joined geometrical elements: Ω 0 ≈ Ω h = ∪ e Ω e as described on fig. 1.2. The apexes of elements Ω e form the nodes of the mesh. The common linear geometrical elements are:

• for one-dimensional meshes, the line segment (the only one);

• for two-dimensional meshes, triangle or quadrangle;

• for three-dimensional meshes, the 8-node hexahedron or the 4-node tetrahedron.

Each family of elements has particular property which can influence the discrete solution. This aspect is not addressed here, for more informations on meshing see [START_REF] Belytschko | Nonlinear Finite Elements For Continua And Structures 2nd[END_REF]. The space U, defined by eq. (1.13), is approximated by a finite dimensional space U h . Its base is (ψ

I ) × (ψ I ) × (ψ I ) I∈[1.
.N ] with ψ the interpolation functions or shape functions. Each continuous field is then approximated by:

u(X , t) ≈ u h (X , t) = 1 i N [u ix (t) u iy (t) u iz (t)] t ψ i (X )
(1.17)

[u ix (t) u iy (t) u iz (t)]
is the interpolation coefficient associated to the interpolation point i, and ψ i the interpolation function associated to i. The interpolation functions are the same for each space dimension (x, y, z). Low-order polynomials are commonly used. If each [u ix (t) u iy (t) u iz (t)] is known, the discrete field u h is completely determined. The [u ix (t) u iy (t) u iz (t)] forms the degrees of freedom (DOF) of the system, here three by interpolation point.

On Ω h , interpolation points are places on the geometrical elements. A geometrical element together with shape functions and interpolation points define a finite element. For example, the linear isoparametric finite element uses: linear geometrical elements (hexaedron or tetrahedron), interpolation points at the nodes, and linear polynomials as interpolation functions. The nodal values of a discrete field are then directly the interpolation coefficients, and form the degrees of freedom. The linear isoparametric type is used in the following. For finite elements of higher order, the interpolation functions are polynomials with a higher degree and extra interpolation points are added inside the geometrical elements.

The space V, defined by eq. (1.14), is also approximated by a finite dimensional space V h . The same shape functions (ψ I ) × (ψ I ) × (ψ I ) are used to form a basis of V h , but those associated to the DOFs of Γ D are not considered. The discrete-in-space form of the PVW (1.16) is obtained by δu which ranges all basis functions of V h . This gives an equation for each degree of freedom of u h (three per node). So for the node (or DOF) i, δu takes the values:

δu =      δu ix = [1 0 0] t ψ i
for selecting the DOF for x δu iy = [0 1 0] t ψ i for selecting the DOF for y δu iz = [0 0 1] t ψ i for selecting the DOF for z

(1.18)
With these approximations, the PVW (1.16) turns into a matrix system:

M Ü = F ext -F int (U) (1.19)
With :

M 3i..3i+2 3j..3j+2 = Ω h ρψ i ψ j dΩ h    1 0 0 0 1 0 0 0 1    (1.20) F int 3i..3i+2 =    F int ix F int iy F int iz    F ext 3i..3i+2 =    F ext ix F ext iy F ext iz    (1.21) F int ik = Ω h T r P (u)(∇δu ik ) t dΩ h , k ∈ {x, y, z} (1.22)
F ext ik = Ω h f v • δu ik dΩ h + Γ N t • δu ik dΩ h , k ∈ {x, y, z} (1.23)
U k=3i..3i+2 = u ix u iy u iz t (1.24)
Ük=3i..3i+2 = üix üiy üiz

t (1.25)
Note that the place of DOFs in the global vector U depends on the implementation. Here the DOFs are concurrent just for sake of simplicity. With these approximations, the estimation of integrals of PVW are facilitated. The integral over Ω h is computed by adding the contribution of the integrals over each Ω e :

Ω h f (X)dΩ = e Ωe f (X)dΩ (1.26)
In order to compute the integral over elements, a map Ξ e : Ω ref → Ω e is defined between each Ω e (in initial configuration) and a parent or a reference element Ω ref .

For a type of finite elements, this parent element is unique. Each type of integral contribution is calculated once on Ω ref . Thanks to the Jacobian of Ξ e or its gradient, the integral contribution on Ω ref is transferred on each Ω e . The element integral are then added in the global vectors. This process is called the assembly, for more details see [START_REF] Belytschko | Nonlinear Finite Elements For Continua And Structures 2nd[END_REF]. Finally, the discrete problem in space is:

M Ü = F ext -F int (U) t ∈ [t 0 , t f ] (1.27) U(t) Γ D = Û(t) t ∈ [t 0 , t f ] (1.28) U(0) = U 0 U(0) = U0 t ∈ [t 0 , t f ] (1.29)
The eq. (1.27) is the discrete dynamic, which directly integrates the Neumann's conditions in F ext . The eq. (1.28) enforces the Dirichlet's conditions. In order to start the solving, the initial conditions of eq. (1.29) are required.

Here the internal forces depend only on displacement, but for certain material laws a dependence in velocity is possible.

The Lagrangian mechanics formalism

The discrete-in-space eqs. (1.27) to (1.29) can be obtained by an other formalism, the variational discrete mechanics. For details on the formalism see [START_REF] Marsden | Discrete mechanics and variational integrators[END_REF], only the main steps are briefly exposed here. The considered system is already discrete in space, and conservative (no friction or plasticity are considered). The degrees of freedom, gathered in the global vector U, are determined by the discrete Lagrangian:

L(U, U) = T ( U) -V(U) (1.30)
With T ( U) the kinetic energy, and V(U) the potential associated to internal and external stresses:

T ( U) = 1 2 Ut M U (1.31) - ∂V(U) ∂U = F ext -F int (U) (1.32)
S is an action integral defined on time interval [t 0 , t f ] by:

S(U, U) = t f 0 L(U, U)dt (1.33)
A stationary S corresponds to the solution of the discrete-in-space problem of eqs. (1.27) to (1.29). This solution is then defined by:

U ∈ U h | δS(U, U) = 0 (1.34)
By replacing S in eq. (1.34) with its expression (1.33), the system to solve is:

∀δU ∈ V h , t f 0 δU t ∂L ∂U + δ Ut ∂L ∂ U dt = 0 (1.35)
By an integration by parts on the second term, and using δU(t) = 0 for t ∈ {0, t f }, the eq. (1.35) gives:

∀δU ∈ V h , t f 0 δU t ∂L ∂U - d dt ∂L ∂ U dt = 0 ⇔ ∂L ∂U - d dt ∂L ∂ U = 0 (1.36)
The eq. (1.36) gives directly the dynamics of eq. (1.27) with the expressions in eqs. (1.31) and (1.32) of T and V.

The Hamiltonian mechanics

The Hamiltonian mechanics is an equivalent formalism to the Lagrangian mechanics. But it allows to introduce in an attractive way some fundamentals properties of mechanics. Here, only the main ideas are briefly exposed. For more details, see [START_REF] Marsden | Discrete mechanics and variational integrators[END_REF].

The Legendre transform relates the Lagrangian to the Hamiltonian. It introduces the momentum by the expression:

FL : (U, U) → (U, P) = U, ∂L ∂ U (1.37)
The expression of the Hamiltonian is then:

H(U, P) = Ut P -L(U, U) (1.38)
U and P being related by FL. The partial derivatives of H are:

∂H ∂U = -Ṗ (1.39) ∂H ∂P = U (1.40)
The eqs. (1.39)

Ẋ = J ∂H ∂X with J = 0 1 -1 0 (1.42)
The Hamiltonian's equations of eq. (1.42) defines a symplectic form with J.

For the demonstration, the reader will refer to [START_REF] Simo | Exact energy-momentum conserving algorithms and symplectic schemes for nonlinear dynamics[END_REF] or [START_REF] Marsden | Discrete mechanics and variational integrators[END_REF]. Simo et al. in [START_REF] Simo | Exact energy-momentum conserving algorithms and symplectic schemes for nonlinear dynamics[END_REF], and then Kane et al. in [START_REF] Kane | Symplectic-energy-momentum preserving variational integrators[END_REF] identify the symplectic feature as a key one for time-integration. They build discrete time-integrator which keep this property and demonstrate very good properties.

Remark 1. The Hamiltonian and Lagrangian formalism are used in a fundamental theorem of mechanics: the Noether's theorem. The space of solutions is invariant under specific transformations, which express the invariances of the physical model. A translation or a rotation on the coordinates, does not change the physical system, as well as a translation in time. In the Hamiltonian or in the Lagrangian, these invariances lead to the fundamental principles of dynamics: the translational invariance to the dynamics equilibrium, the rotational invariance to the conservation of angular momentum, and the time invariance to the conservation of energy. The Noether's theorem demonstrates then that the dynamics is ruled only by these principles.

A symplectic time-integrator will preserve the corresponding discrete quantities. The discrete system behaves then as the continuous one, which ensures a precise time-integration.

1.2 The problem of impact mechanics

1.2.

Description of contact

Contact kinematics

In order to enforce contact conditions on a deformable body Ω 0 , its boundary ∂Ω 0 is now split into three parts: Γ D , for Dirichlet's conditions; Γ N , for Neumann's ones; and Γ C the part potentially in contact. As before, they form a complementary set:

∂Ω 0 = Γ D ∪ Γ N ∪ Γ C and Γ D ∩ Γ N = Γ D ∩ Γ C = Γ N ∩ Γ C = ∅ (1.43)
When two bodies Ω

0 and Ω

0 are involved, the contact relations link the contact boundaries Γ C . As the contact relations involve the deformed configurations, the counterparts of 1. a rigid-deformable contact between a rigid and motionless body, and a deformable one;

Γ (i) C is γ (i) C in deformed configur- ation. Γ R n γ C (t) Ω(t) Γ D (t) Γ N (t) n γ C,1 (t) Ω 1 (t) Γ D,1 (t) Γ N,1 (t) γ C,2 (t) Ω 2 (t) Γ D,2 (t) Γ N,2 (t)
2. a deformable-deformable contact between two deformable bodies.

They are described on fig. 1.3. Here the name "rigid-deformable" contact is restricted to a motionless rigid body for clarity. The deformable-deformable contact can anyway describe a rigid and a deformable moving body, or even two rigid bodies. Whatever the type of contact, they are described by the same kinematic quantities. The gap measures the distance between contact boundaries γ

C : ∀x (1) ∈ γ (1) C , g(x (1) ) = x (1) -x (2) (x (1) ) • n (1.44)
n is represented on fig. 1.3 as the normal at contact. The normal n must be chosen as either the inner normal of γ

C , or the outward normal of γ

C . In the continuous setting the choice has a small influence because the both definitions merge when contact happens. But in the discrete setting, this choice is much more important [see 80]. x (2) is the closest point on γ

(2)

C from x (1) . Multiple methods exist to get this closest points [see 80] but they are not investigated in this work. In the following, a map is stated between γ

C → γ (2) C x (1) → x (2) = π 21 (x (1) ) (1)
The relative normal velocity v N is the approaching velocity between γ

(1) C and γ

(2)

C : v N (x (1) ) = u(x (1) ) -u(x (2) ) • n (1.46)
For a rigid-deformable contact, as the second body is motionless, v N is then only u(1) projected along n. Both g(x (1) ) and v (x (1) ) are defined relatively to one body, here Ω (1) , and its contact boundary γ

C . In the following, γ

C is denoted γ C as the reference surface and the associated subscripts are omitted. Sometimes the contact boundary γ

C is referred as the slave surface and γ

(2) C as the master surface. The choice of the slave and master bodies depend on spatial discretization, or method used to enforce contact. In [START_REF] Wriggers | Computational Contact Mechanics 2nd[END_REF], some indications are given to lead this choice.

In order to deal with friction, a relative tangential velocity is defined as the part of relative velocity in the tangent plane to n:

v T = [ u(1) -u(2) ] -v N n (1.47)

Contact constraints

The contact is described by the following conditions [see 23] on Ω(t):

• the impenetrability condition, which states that bodies can not penetrate each others,

∀(x , t) ∈ γ C × [0, t f ], g(x ) 0 (1.48)
• the intensility condition, which states that the normal contact stress λ N is only compressive (i.e. no adhesion occurs),

∀(x , t) ∈ γ C × [0, t f ], λ N = (σn) • n 0 (1.49)
• the complementary condition, which reduces the contact states either to an active contact (λ N 0 and g = 0), or a non-active contact (λ N = 0 and g > 0),

∀(x , t) ∈ γ C × [0, t f ], λ N g(x ) = 0 (1.50)
The conditions of eqs. (1.48) to (1.50) are equivalently referred to Karush-Kuhn-Tucker (KKT), Hertz-Signorini-Moreau (HSM), or Signorini's conditions.

The conditions of eqs. (1.48) to (1.50) can be gathered in one complementary form:

∀(x , t) ∈ γ C × [0, t f ], 0 g(x ) ⊥ λ N 0 (1.51)
As showed by Moreau in [START_REF] Moreau | Numerical aspects of the sweeping process[END_REF], the condition of eq. (1.51) is equivalent to the following one:

∀(x , t) ∈ γ C × [0, t f ], if g(x ) > 0, λ N = 0 if g(x ) = 0, 0 v N ⊥ λ N 0 (1.52)
The condition of eq. (1.51) relies on the gap, which directly depends on displacement. The eq. (1.52) gives a condition directly on velocity. One of these formulations is chosen depending on the main unknown of the problem. Sometimes these two forms of Signorini's conditions are referred as impenetrability condition for the form of eq. (1.51), and persistency condition for eq. (1.52). In fact, the eq. (1.51) directly imposes the non-penetration between bodies. And the eq. (1.52) prevents the penetration by stopping the contact boundaries. The conditions eqs. (1.48) to (1.50) govern only the normal part λ N of contact force. No friction occurs in tangential direction to Γ C . This type of contact is called unilateral contact, or frictionless contact.

In order to deal with friction, the tangential direction t is set as the opposite direction of v T defined by eq. (1.47):

t = - v T v T if v T = 0 (1.53)
t is the direction of the tangential contact force λ T :

λ T = λ T t, if λ T 0 (1.54)
v T and λ T are related thanks to the Coulomb law:

∀(x , t) ∈ γ C × [0, t f ], (1.55)
If

λ N = 0, λ T = 0 If λ N > 0, 0 (µλ N -λ T ) ⊥ v T 0 (1.56) With λ T = -αv T , α ∈ R + (1.57)
µ is the friction coefficient depending on materials and surfaces involved in contact. The Coulomb's law of eq. (1.57) is a simple model to compute the stress due to to friction, but certainly the most used. It contains two states:

• sticking, when λ T is lower than the threshold µλ N , no slip happens meaning that the relative tangential velocity is null (v T = 0);

• sliding, when v T > 0 and λ T is equal to its maximal value µλ N whatever the tangential velocity.

1.2.2 Governing equations of impact problem -The nonsmooth Lagrangian mechanics

Γ N Γ R Γ D n γ C Figure 1.4: A discrete in space contact problem
For the developments in this section, the considered problem is described by fig. 1.4. A discrete-in-space deformable body is subjected to an unilateral contact on its boundary γ C . The boundary of the body Ω 0 is divided into three parts: Γ D with Dirichlet's conditions, Γ N with Neumann's ones, and Γ C where a unilateral contact is enforced. Γ D , Γ N and Γ C form a partition of ∂Ω, and denote here both the boundary and the sets of nodes in it. The finite elements are P1-isoparametric in order to localize the degrees of freedom only on nodes. As the contact is only normal, the subscript N is dropped.

On γ C , the deformed configuration of Γ C , each node i must meet the Signorini's conditions. They can be enforced on gap:

∀i ∈ γ C , 0 g(U) i ⊥ λ i 0 (1.58)
or on normal velocity:

∀i ∈ γ C , if g(U) i > 0, λ i = 0 if g(U) i = 0, 0 v i ) ⊥ λ i 0 (1.59)
As the elements are P1-isoparemetric, this nodal enforcement ensures the contact conditions everywhere. But for finite elements of higher order, this is not necessary true.

The vector gathering the normal velocities for each nodes in Γ c can be computed as:

v = L U with L = ∇g(U) (1.60)
The operator L projects the vector over the global domain as U or U on the inward normal n to γ C . The space for solutions C is the restriction of the previous one U h to solutions satisfying contact conditions:

C = {u ∈ U h |g N (u) γ C 0} (1.61)
The Signorini's conditions are enforced for each nodes of γ C by a new term in the Lagrangian:

L (U, U, λ) = T ( U) -V(U) + λ t g N (U) (1.62)
λ is the Lagrange multiplier which enforces the condition on γ C , with λ ∈ (R + ) d with d the number of nodes in γ C . For details on constraint enforcement on discrete Lagrangians, see [START_REF] Marsden | Discrete mechanics and variational integrators[END_REF].

In order to obtain the dynamics from the Lagrangian of eq. (1.62), U needs to be twice-derivable in time. If not, the acceleration is not defined. But a contact problem has at least two times where the velocity is discontinuous:

• firstly at impact time, the nodal velocity along contact normal goes from a negative value to a null one;

• secondly at release time, it passes from a null value to a positive one.

These two non-smooth times must be considered for each node, and each contact phase.

In order to overcome this difficulty, Fetecau et al. extends the variational framework to non-smooth cases in [START_REF] Fetecau | Nonsmooth Lagrangian Mechanics and Variational Collision Integrators[END_REF]. A new action integral is introduced depending on the non-smooth times:

S(U, U, λ, t i ) = t i 0 L (U, U, λ) dt + t f t i L (U, U, λ) dt, U ∈ C (1.63)
Only one non-smooth time is considered here for shortness and clarity, the extension to multiple non-smooth times is straightforward. In order to enforce the contact constraint at t i , Cirak and West [START_REF] Cirak | Decomposition contact response (DCR) for explicit finite element dynamics[END_REF], (see [START_REF] Fekak | A new heterogeneous asynchronous explicit-implicit time integrator for nonsmooth dynamics[END_REF] for an other formalisation) adds a Lagrange multiplier in the action integral of eq. (1.63):

S(U, U, λ, t i ) = t i 0 L (U, U, λ) dt + t f t i L (U, U, λ) dt + λ t t i g(U(t i )) (1.64
) Note that at a non-smooth time, U(t i ) satisfies g(U(t i )) = 0. Indeed, for both impact and release, the body is in contact. For the persistent contact between impact and release, the enforcement of contact constraint is done by the Lagrange multiplier inside the Lagrangian (1.62).

The solution to the problem is given by a minimization of S:

λ, U | ∀ δλ, δU, δt i , δ S(U, U, λ, t i ) = 0 (1.65)
The eq. (1.65) gives:

δ S = [0,t i [ δL (U, U, λ) dt + ]t i ,t f ] δL (U, U, λ) dt + (L t + i -L t - i )δt i + δλ t t i g(U(t i )) + λ t t i δg(U, t i ) With: • δg(U) = ∇g δU • δU = 0, t ∈ {0, t f } • t b ta δL (U, U, λ) dt = t b ta δU t ∂L ∂U - d dt ∂L ∂ U + δλ t g(U) + λ t δg(U) dt + ∂L ∂ U δU ta t b = t b ta δU t ∂L ∂U - d dt ∂L ∂ U + (∇g) t λ + δλ t g(U) dt + ∂L ∂ U δU ta t b
As proven in [START_REF] Fetecau | Nonsmooth Lagrangian Mechanics and Variational Collision Integrators[END_REF] and [START_REF] Cirak | Decomposition contact response (DCR) for explicit finite element dynamics[END_REF], the differential for the constraint g(U, t i ) = 0 is:

δg(U, t i ) = (∇g) (δU + dU t i ) (1.66)
With dU t i , the differential of U at t i :

dU t i = U(t + i )δt i ∀δt i > 0 U(t - i )δt i ∀δt i < 0 (1.67)
Finally, the differential of S is:

δ S = [0,t i [∪]t i ,t f ] δU t ∂L ∂U - d dt ∂L ∂ U + (∇g) t λ + δλ t g(U) dt -δU t ∂L ∂ U t + i t - i -∇g t λ t i + δt i L t + i t - i + dU t t i (∇g t λ t i ) + δλ t t i g(U, t i ) (1.68)
Replacing δ S by its expression of eq. (1.68) in the eq. (1.65) gives the following equations:

∀t ∈ [0, t i [∪]t i , t f ], ∂L ∂U - d dt ∂L ∂ U + (∇g) t λ = 0 (δU) (1.69) ∀t ∈ [0, t i [∪]t i , t f ], g(U) 0 (δλ) (1.70) t = t i , ∂L ∂ U t + i t - i = ∇g t λ t i (δU) (1.71) t = t i , L + Ut ∂L ∂ U t + i t - i = 0 (δt i ) (1.72) t = t i , g(U, t i ) = 0 (δλ) (1.73)
The eq. (1.69) corresponds to the dynamics for smooth time intervals, and the eq.

∀t ∈ [0, t i [∪]t i , t f ], M Ü = F ext -F int (U) + L t λ (δU) (1.74) t = t i , M U(t + i ) -U(t - i ) = L t λ t i (δU) (1.75) t = t i , 1 2 Ut M U t + i t - i = 0 (δt i ) (1.76)
The Lagrange's multipliers λ and λ t i differ: λ is a nodal force, and λ t i an impulsion. In fact the eq. (1.75) defines λ t i as a velocity jump multiplied by a mass. And the eq. (1.74) gives λ as an acceleration multiplied by a mass.

In order to get closer to a physical case, the following problem is considered. A deformable body with only one contact node is impacting a rigid and motionless boundary. There are only one impact time t i , and one release time t r . Between t i and t r , the contact node stays on the rigid boundary. But outside [t i , t r ] no contact happens. This leads to 3 smooth phases: before t i with no contact, between t i and t r with persistent contact, after t r with no contact. The equations describing the dynamics are then:

∀t ∈ [0, t i [, M Ü = F ext -F int (U) (1.77) t = t i , M U(t + i ) -U(t - i ) = L t λ t i (1.78) ∀t ∈]t i , t r [, M Ü = F ext -F int (U) + L t λ (1.79) t = t r , M U(t + r ) -U(t - r ) = L t λ tr (1.80) ∀t ∈]t r , t f ], M Ü = F ext -F int (U) (1.81)
With the two following energy relations for t i and t r , meaning that the kinetic energy is conserved through the jumps:

1 2 Ut M U t + i t - i = 0 1 2 Ut M U t + r t - r = 0 (1.82)

The nonsmooth contact dynamics

The non-smooth contact dynamics (NSCD) is a suitable framework for impact dynamics. It was introduced by Moreau and Jean in [START_REF] Jean | The non-smooth contact dynamics method[END_REF] and [START_REF] Moreau | Numerical aspects of the sweeping process[END_REF]. The dynamics is written with a measure on velocity, defined for both smooth and non-smooth events. And they introduce the Signorini's conditions on velocity of eq. (1.52), in order to enforce the contact. If a time interval [0, t f ] is considered with only one non-smooth event at time t i , the measure of velocity is defined as:

d U = Üdt, t ∈ [0, t i [∪]t i , t f ] U(t + i ) -U(t - i ), t = t i (1.83)
The eq. (1.74) describing the dynamics for smooth phases and eq. (1.75) for non-smooth events can be gathered in one thanks to the measure of eq. (1.83):

Md U = F ext -F int (U) dt + L t dr (1.84)
The measure dr is a measure on contact actions:

dr = λdt, t ∈ [0, t i [∪]t i , t f ] λ t i , t = t i (1.85)
This definition underlines that λ and λ t i are different quantities: λ is a nodal force (in N); and λ t i is an impulse (in N s), homogeneous to a force over a time interval.

Finally the dynamics equation together with the Signorini's conditions on gap are:

Md U = F ext -F int (U) dt + L t dr ∀i ∈ γ C , 0 g(U) i ⊥ λ i 0 (1.86)
and if the Signorini's conditions are enforced on velocity:

Md U = F ext -F int (U) dt + L t dr ∀i ∈ γ C , If g(U) i > 0, (λ) i = 0 Else g(U) i = 0, 0 (v) i ⊥ (λ) i 0 (1.87)
1.3 Time-integrators for transient dynamics of continua

Time-integration in smooth dynamics

Properties of time-integration schemes

The eqs. (1.86) and (1.87) describe a discrete-in-space problem. In order to solve it numerically, a discretization in time is required. The time-integrator realizes this discretization: it links the displacement U, to its derivatives ( U and Ü) for discrete times. The discrete times {t n+1 } form a time grid which covers the interval [t 0 , t f ]. Each discrete time is separated here by a constant time-step h = t n+1 -t n . Thanks to the time-integrator, the dynamics of eq. (1.84) can be expressed with only one variable (displacement or velocity) making the system solvable.

Designing a time-integrator is complex, and still an open subject of research. Due to the diversity of problems and applications, the time-integration schemes are numerous. But they are all designed to meet some crucial properties like stability, accuracy or efficiency. These properties are reached at diverse degrees depending on the required performances of the time-integrator.

Stability. The stability names here the numerical stability. A time-integrator is stable if a small perturbation in the initial state leads to a bounded change in the whole numerical solution. In the linear regime, the properties of the amplification matrix determines generally the stability of the system. In brief, the higher natural frequency of the discrete system sets the stability condition. And the natural frequencies are directly related to the eigenvalues of the amplification matrix. [START_REF] Belytschko | Nonlinear Finite Elements For Continua And Structures 2nd[END_REF] gives more details and some other methods to analyse the stability. An other interesting approach is developed by Krenk in [START_REF] Krenk | Energy conservation in Newmark based time integration algorithms[END_REF]. He analyses the stability thanks to the discrete energy balance, which extends the stability notion beyond the linear regime to the non-linear one. The stability condition is transformed into a stable time-step. If the time-step of the simulation is higher than the stable one, the stability is no more ensured. The mesh is the main parameter for the user to influence the stable time-step. Generally the more accurate is the mesh, the shorter the stable time-step is. Indeed for a given material a shorter characteristic length leads to a higher natural frequency. The time and space discretizations are then linked: finer is the mesh, smaller is the stable time-step. The material characteristics are the second main parameter which determines the stable time-step. But they are generally imposed by the model. Sometimes the stability is included in a more practical notion: the robustness. It refers more to the range of stability for the scheme. For example, a scheme is said robust if it is stable for a wide range of boundary conditions, of materials laws, or for large perturbations in the system... Accuracy. A time-integration scheme approximates the continuous time integration. The accuracy is the error of this approximation. It depends on the time-step h: if the error is in O(h), the scheme is first-order accurate; if the error is in O(h 2 ), the scheme is second-order accurate. The accuracy and stability are related. Even if the stable time-step is large, the simulation time-step is constrained by the accuracy. The simulation computes only the phenomenons whom the characteristic times are greater than the time-step.

Symplecticity. According to the Noether's theorem, a mechanical system described by a Lagrangian conserves three momentums: the linear momentum, the angular momentum and the energy. The three conservation relations correspond to three equations: the dynamics for the conservation of the linear momentum; the equation of angular momentum conservation; and the energy balance. For the discrete mechanics (both in space and time), the symplectic schemes [see 47, 59, for a definition] keep notably this property: discrete angular momentum and discrete energy stay constant over the simulation. But, as mentioned in [START_REF] Kane | Symplectic-energy-momentum preserving variational integrators[END_REF], they are not necessarily equal to the continuous angular momentum and energy. An exact energy-momentum conserving scheme, on the other hand, conserves exactly the energy and angular momentum according to the continuous definition [see 73, for example]. If some exact energy-momentum schemes are not symplectic according to the definition of [START_REF] Kane | Symplectic-energy-momentum preserving variational integrators[END_REF], both families have schemes in common and are then often confused in literature. For a comprehensive introduction about discrete vari-ational mechanics, and a review on some symplectic schemes, the reader will refer to [START_REF] Marsden | Discrete mechanics and variational integrators[END_REF]. As mentionned in section 1.1.3, the symplectic feature for the discrete schemes is a translation of the invariants of continuous Hamiltonian. The Hamiltonian for solid mechanics conserves only the linear momentum, the angular momentum and the energy. But if the thermics plays a role in the physical model, new invariants of the Hamiltonian lead to the laws of thermodynamics. Some recent time-integrators extend the symplectic feature to these invariants as [START_REF] Romero | Algorithms for coupled problems that preserve symmetries and the laws of thermodynamics. Part I: Monolithic integrators and their application to finite strain thermoelasticity[END_REF].

Efficiency. A scheme is said efficient if the computational cost associated to a time-step is low. If the major levers for increasing efficiency depend on computer science,the scheme still plays a role. Indeed the time-integrator sets the shape of the matrix system solved at each time-step. Two types of time-integrators are then distinguished:

• the explicit time-integrators, where no system solving is required because the matrix could be diagonalized;

• the implicit time-integrators, where the matrix is non-diagonal leading to a system solving at each time-step.

For an implicit scheme, the computational cost is then much higher for one time-step. But besides the computational cost of one time-step, the stable time-step influences the computational cost of the whole simulation.

If the explicit schemes are more efficient regarding to one time-step, they are conditionally stable. They require then a minimum number of time-steps for a given time-interval. This limit is overcome by some implicit schemes which are unconditionally stable. Finally, if the whole time-interval is considered, the higher number of time-steps for an explicit scheme can lead to a larger computational cost. The efficiency is influenced also by the linearity of the system solved at each time-step. Some time-integrators lead to non-linear systems, which require time-consuming algorithms as Newton-Raphson's method; while other ones involve only linear system. Even for non-linear models, the explicit time-integrators only involve a diagonal and linear system. This major advantage was early identified to solve non-linear transient systems [START_REF] Belytschko | Explicit algorithms for the nonlinear dynamics of shells[END_REF]. The time-integrator influences also the scalability of the algorithm by its capacity to be parallelized [see 16, 29, 53, for example].

The choice of the time-integrator is then crucial for programming a finite element software, as a first but major lever on efficiency.

The Newmark's time-integrators

The Newmark's schemes are crucial in time-integration for structural dynamics. Introduced by Newmark in [START_REF] Newmark | A Method of Computation for Structural Dynamics[END_REF], they are the most used onese because of their robustness, the simplicity of the algorithm, and the numerous schemes in this family. Newmark defines the time-integrator thanks to two formulas:

U n+1 = U n + h Un + h 2 1 2 -β Ün + β Ün+1 (1.88) Un+1 = Un + h (1 -γ) Ün + γ Ün+1 (1.89)
The first eq. ( 1.88) gives the updated displacement at t n+1 depending on quantities at t n , and the acceleration at t n+1 . The second one gives the updated velocity in the same way. Each equation depends on a weighting term, γ for velocity formula, and β for displacement one. A scheme is determined by the choice of a couple (γ, β). Together with the time-integration relations of the eqs. (1.88) and (1.89), the following discrete dynamics forms the scheme equations:

M Ün+1 = F ext (t n+1 ) -F int (U n+1 , Un+1 ) (1.90)
In the following, the dynamics system is assumed to only have linear elastic materials and no contact. The internal forces are then directly linked to the displacement by the stiffness matrix K, and the discrete dynamics is:

M Ün+1 = F ext (t n+1 ) -KU n+1 (1.91)
Together with the time-integration relations (1.88) and (1.89), it leads to the following linear system:

M Ün+1 = R( Ün , Un , U n ) (1.92) With: R( Ün , Un , U n ) = F ext (t n+1 ) -K U n + h Un + h 2 1 2 -β Ün (1.93) M = M -h 2 βK (1.94)
The matrix M combines the mass matrix M and the stiffness one K. Its shape is generally not diagonal and couples the rows of the system (1.92).

In a time-step, the discrete solution is computed in three main steps:

1. computing Ün+1 by solving the linear system of the eq. (1.92);

2. computing U n+1 with the eq. (1.88)

3. computing Un+1 with the eq. (1.89)

Solving the system (1.92) requires only a linear solver. If the internal forces F int (U n+1 , Un+1 ) are not linear, a non-linear solver is required like the Newton-Raphson method. For details on non-linear models, and more generally on implicit solving see [START_REF] Belytschko | Nonlinear Finite Elements For Continua And Structures 2nd[END_REF].

In the linear regime without damping, the stability analysis leads to a family of unconditionally stable schemes:

β 1 2 γ 1 4 .
For the other ones the stability condition is:

h c = 1 ω max 1 γ 2 -β γ - 1 2 -β - 1 4 , and γ - 1 2 0 (1.95)
with ω max the higher natural frequency of the discrete system (1.92). For a proof see for example [START_REF] Belytschko | Nonlinear Finite Elements For Continua And Structures 2nd[END_REF]. The stability conditions are more complex in case of damping, coupled systems, or non-linearities. The order of convergence is in

O(h 2 ) if γ = 1 2 , and in O(h) if γ = 1 2 [see 41].
The average acceleration method. Corresponding to the couple (γ = 1/2, β = 1/4), this implicit scheme is certainly the most common in computational mechanics. Sometimes it is denoted by an other name: the Crank-Nicholson method. Its numerical properties make it very interesting: unconditionally stable and yet precise with a convergence order of O(h 2 ); and with an energy balance matching the continuous one. As shown by Krenk in [START_REF] Krenk | Energy conservation in Newmark based time integration algorithms[END_REF], the energy balance for the average acceleration method is:

1 2 Ut n+1 M Un+1 + 1 2 U t n+1 KU n+1 t n+1 tn = (U n+1 -U n ) t F ext (t n+1 ) -F int (U n+1 , Un+1 ) (1.96) with: X(t n+1 ) = 1 2 (X n+1 + X n ) (1.97)
The terms in the eq. (1.96) are easily identified:

• E k = 1 2 Ut n+1 M Un+1 , the discrete kinetic energy at t n+1 ; • E int = 1 2 U t n+1 KU n+1
, the discrete internal energy of linear internal forces;

• ∆W ext = (U n+1 -U n ) t F ext (t n+1 ) , the work of external forces on [t n , t n+1 ]; • ∆W int = (U n+1 -U n ) t -F int (U n+1 , Un+1
) , the work of non-linear internal forces.

These terms correspond directly to the physical ones. No extra numerical term appears in this energy balance.

The Central Difference method. An other interesting scheme is the Central Difference method obtained with the coefficients γ = 1/2 and β = 0.

The order of convergence is then in O(h 2 ), and the stability condition is:

h c = 2 ω max (1.98)
If the eq. ( 1.98) gives the theoretical stable time-step for a linear system, it is never used in practice. The maximum frequency ω max is the square root of the maximum eigenvalue of the linear system. The computation of these eigenvalues is too time-consuming for large systems. The stability timestep is then approximate by an other method. Following [START_REF] Belytschko | Nonlinear Finite Elements For Continua And Structures 2nd[END_REF], the element eigenvalue inequality gives an easy access to an upper bound: the maximal eigenvalue is inferior to the maximal eigenvalue of elementary matrices. The estimation is then done element by element with a simple formula based on the Courant-Friedrichs-Lewy criteria [START_REF] Courant | Uber die partiellen Differenzengleichungen der mathematischen Physik[END_REF]:

h c min e l e c e (1.99)
This formula is exact only for 1D meshes: l e is the initial length of element, and c e is the wave velocity depending only on the deformation of the element and material properties [see 9]. The eq. ( 1.99) is extended to 3D meshes by taking l e as the minimal distance between any two nodes of the elements.

According to [START_REF] Belytschko | Nonlinear Finite Elements For Continua And Structures 2nd[END_REF], a correction between 5% and 20% is applied for non-linear systems depending on dynamics. This method is highly suited to finite elements: the stability time-step is approximated by a simple formula looped over all the elements of the mesh. The Newmark's relations of eqs. (1.88) and (1.89) together with the space discrete dynamic of eq. (1.90) leads to the following equations:

M Ün+1 = F ext (t n+1 ) -F int (U n+1 , Un+1 ) (1.100) Un+1 = 1 2h (U n+2 -U n ) (1.101) Ün+1 = 1 h 2 (U n+2 -2U n+1 + U n ) (1.102)
The eqs. (1.100) to (1.102) form the displacement formulation of the CD method. Replacing Ün+1 by its expression of eq. (1.102) in the dynamics (1.100) leads to the following system:

MU n+2 = h 2 F ext (t n+1 ) -F int (U n+1 , Un+1 ) +2MU n+1 -MU n (1.
103)

The linear system (1.103) is written on displacement, but involves three time-steps: U n+2 , U n+1 and U n . For this reason, the following equivalent formulation is preferred.

An other formulation appears indeed by introducing the velocities at midpoint of time-interval [see 9]:

Un+ 3 2 = 1 h (U n+2 -U n+1 ) (1.104)
which gives the following equations:

U n+1 = U n + h Un+ 1 2 (1.105) M Un+ 3 2 -Un+ 1 2 = h F ext (t n+1 ) -F int (U n+1 , Un+ 1 2 ) (1.106) Un+1 = 1 2 Un+ 3 2 + Un+ 1 2 (1.107) Ün+1 = 1 h Un+ 3 2 -Un+ 1 2 (1.108)
The discrete solution can be computed only with the time-integration relation of the eq. ( 1.105), and the discrete dynamics (1.106). The velocity and the acceleration of the eqs. (1.107) and (1.108) are useful only for processing the results.

The linear system comes from the eq. ( 1.106):

M Un+ 3 2 = M Un+ 1 2 + h F ext (t n+1 ) -F int (U n+1 , Un+ 1 2 ) (1.109)
The matrix associated to this linear system is simply the mass matrix M.

A common technique is to diagonalize this matrix by a lumping technique [START_REF] Belytschko | Nonlinear Finite Elements For Continua And Structures 2nd[END_REF][START_REF] Hughes | The finite element method: linear static and dynamic finite element analysis[END_REF]. Indeed with a diagonal M, solving the system of eq. (1.109) is then explicit. The lumping technique are numerous, but they all conserve the total mass. The simplest ones involve only summations of columns in the diagonal term. A diagonal mass matrix can also be directly get from the spatial integration, with a specific location for the Gauss integration points. This approximation does not impact the convergence rate [START_REF] Wu | Lumped mass matrix in explicit finite element method for transient dynamics of elasticity[END_REF], and increases slightly the stable time-step [START_REF] Belytschko | Nonlinear Finite Elements For Continua And Structures 2nd[END_REF]. An other key approximation is made in the integration of the internal forces F int in eq. (1.106). If in general they depend on U n+1 and Un+1 , they are here expressed thanks to Un+ 1

2

. This avoids a dependence on Un+ 3 2 in the right side of the eq. ( 1.106) and preserves the explicit property of the scheme. This approximation was proposed by Belytschko et al., see for example [START_REF] Belytschko | Nonlinear Finite Elements For Continua And Structures 2nd[END_REF]. The CD method stays then explicit and linear even for a non-linear material model, making this scheme well suited to non-linear simulations.

Following [START_REF] Krenk | Energy conservation in Newmark based time integration algorithms[END_REF] or [START_REF] Fekak | Etude de la réponse dynamique des ponts roulants soumis à des chocs multiples pendant un séisme : Co-simulation implicite / explicite multi-échelle en temps pour la dynamique du contact[END_REF], the energy balance between [t n , t n+1 ] for central difference method is:

1 2 Ut n+1 M Un+1 - 1 8 h 2 Üt n+1 M Ün+1 t n+1 tn = (U n+1 -U n ) t F ext (t n+1 ) -F int (U n+1 , Un+ 1 2 ) (1.110)
A new term appears, the complementary energy:

E comp = - 1 8 h 2 Üt n+1 M Ün+1 (1.111)
This numerical term comes from the time-integration. It does not dissipate energy but stores a part of it. Because of it, the energy balance does not match exactly the physical energy. The conservative terms are indeed not exactly equal to the system energy: ] the complementary term disappears. The discrete dynamics of eq. (1.106) is multiplied at left by 1 2 ( Un+

E sys = E k + E int (1.
3 2 + Un+ 1 2 ) t : 1 2 Ut n+ 3 2 M Un+ 3 2 t n+ 3 2 t n+ 1 2 = 1 2 h( Un+ 3 2 + Un+ 1 2 ) t F ext (t n+1 ) -F int (U n+1 , Un+ 1 2 ) (1.113)
In this expression, all the terms are again identifiable with physical ones.

Remark 2. Following [START_REF] Simo | Exact energy-momentum conserving algorithms and symplectic schemes for nonlinear dynamics[END_REF], a time-integrator could be resumed in:

G(Z n+1 , Z n ) = 0 with Z = U M U (1.114)
Its amplification matrix is defined by a linearization of eq. (1.114):

δZ n+1 = AδZ n , A = - ∂G(Z n+1 , Z n ) ∂Z n+1 -1 ∂G(Z n+1 , Z n ) ∂Z n (1.115)
As defined by Simo et al. in [START_REF] Simo | Exact energy-momentum conserving algorithms and symplectic schemes for nonlinear dynamics[END_REF], a time-integrator is symplectic if its amplification matrix verifies:

A t JA = J ⇒ det(A) = 1 (1.116)
with J defined in eq. The central difference method is then the only symplectic time-integrator of Newmark's schemes according to this definition. Indeed [START_REF] West | Variational integrators and the Newmark algorithm for conservative and dissipative mechanical systems[END_REF] demonstrates the symplecticity of the Newmark's schemes according to an other definition.

A brief review over time-integrators for smooth dynamics The Newmark's family. As mentioned above, the Newmark's schemes are the most used in structural dynamics. For the implicit ones, the average acceleration method is often preferred. Indeed in the linear regime it is unconditionally stable, and its discrete energy balance is conservative. But structural dynamics requires sometimes a dissipative time-integrator. In fact the spatial discretization introduces spurious high frequencies which pollute the discrete solution. A suitable feature for the time-integrator is then to dissipate these numerical high frequencies. The Newmark's schemes are dissipative if γ > 1/2 [see 39 , 52]. But γ > 1 2 leads to an accuracy in O(h) [START_REF] Hughes | The finite element method: linear static and dynamic finite element analysis[END_REF]. And the scheme dissipates also the low frequencies [START_REF] Chung | A Time Integration Algorithm for Structural Dynamics With Improved Numerical Dissipation: The Generalizedα Method[END_REF][START_REF] Hilber | Improved numerical dissipation for time integration algorithms in structural dynamics[END_REF], which decreases the accuracy of the discrete solution. As a consequence, Hilber et al. and then Chung & Hulbert propose two modifications for the Newmark's scheme. They keep the second-order accuracy and adjust the dissipation depending on frequencies. In [START_REF] Hilber | Improved numerical dissipation for time integration algorithms in structural dynamics[END_REF], Hilber et al. modifies the discrete dynamics by introducing a weighting coefficient on internal forces. The obtained schemes are called HHT-Newmark. In [START_REF] Chung | A Time Integration Algorithm for Structural Dynamics With Improved Numerical Dissipation: The Generalizedα Method[END_REF], Chung & Hulbert adds a similar coefficient but on acceleration resulting in the α-generalized schemes. These generalisation of Newmark's algorithm present a good stability even for non-linear or coupled systems [START_REF] Farhat | Implicit time integration of a class of constrained hybrid formulations-Part I: Spectral stability theory[END_REF]. If both algorithms are implicit, Hulbert & Chung in [START_REF] Hulbert | Explicit time integration algorithms for structural dynamics with optimal numerical dissipation[END_REF] develop an explicit version of the α-generalized scheme for both damped and undamped systems. More recent explicit schemes with the same properties are presented in [START_REF] Kim | An improved explicit time integration method for linear and nonlinear structural dynamics[END_REF][START_REF] Noh | An explicit time integration scheme for the analysis of wave propagations[END_REF][START_REF] Soares | A novel family of explicit time marching techniques for structural dynamics and wave propagation models[END_REF]. The implicit and explicit time-integrators of the Newmark's family present each advantages. In some cases, it is interesting to have several time-integrators working together on a same discrete problem: an explicit and an implicit ones, or two explicit ones with different time-steps. Firstly introduced in [START_REF] Belytschko | Mixed methods for time integration[END_REF][START_REF] Neal | Explicit-explicit subcycling with noninteger time step ratios for structural dynamic systems[END_REF], these method are developed in [START_REF] Gravouil | Multi-time-step explicit-implicit method for non-linear structural dynamics[END_REF][START_REF] Gravouil | Heterogeneous asynchronous time integrators for computational structural dynamics[END_REF][START_REF] Mahjoubi | A monolithic energy conserving method to couple heterogeneous time integrators with incompatible time steps in structural dynamics[END_REF], and applied to contact problems in [START_REF] Fekak | A new heterogeneous asynchronous explicit-implicit time integrator for nonsmooth dynamics[END_REF].

Variational integrators. Simo, Tarnow and Wong introduce in [START_REF] Simo | The discrete energy-momentum method. Conserving algorithms for nonlinear elastodynamics[END_REF][START_REF] Simo | Exact energy-momentum conserving algorithms and symplectic schemes for nonlinear dynamics[END_REF] a family of exact energy-momentum preserving algorithms. They conserve both the energy and the angular momentum at their continuous definitions. The general design process is: to define firstly a family of schemes based on θ-method which preserves the angular momemtum, and secondly to enforce the energy conservation in a way which depends on the model. This design process is applied on problems from rigid body dynamics in [START_REF] Simo | Exact energy-momentum conserving algorithms and symplectic schemes for nonlinear dynamics[END_REF], and on deformable body (FE method) with non-linear material in [START_REF] Simo | The discrete energy-momentum method. Conserving algorithms for nonlinear elastodynamics[END_REF]. An interesting remark in [START_REF] Simo | Exact energy-momentum conserving algorithms and symplectic schemes for nonlinear dynamics[END_REF] is that the only momentum conserving scheme in Newmark's family is the CD method. It conserves the angular momentum even for non-linear problems.

Following the works on exact energy-momentum algorithms, the variational time-integrators are introduced. In [START_REF] Kane | Symplectic-energy-momentum preserving variational integrators[END_REF], Kane et al. present a formalism to get a time-integrator from a discrete Lagrangian both in space and time. The obtained time-integrator is ensured to be symplectic according to the definition of the authors. [START_REF] Marsden | Discrete mechanics and variational integrators[END_REF][START_REF] West | Variational integrators and the Newmark algorithm for conservative and dissipative mechanical systems[END_REF] give more details on continuous and discrete Lagrangian mechanics. A symplectic time-integrator preserves a discrete energy and a discrete momentum, but not necessarily the continuous ones. These quantities comes directly from the discrete Lagrangian. In [START_REF] West | Variational integrators and the Newmark algorithm for conservative and dissipative mechanical systems[END_REF], the Newmark's time-integrators are deduced from a discrete Lagrangian demonstrating they are symplectic. But the conserved quantities are not the continuous ones. This explains why the Newmark's time-integrator are surprisingly good for long-time simulations and non-linear dynamics [START_REF] West | Variational integrators and the Newmark algorithm for conservative and dissipative mechanical systems[END_REF]. In [START_REF] Marsden | Discrete mechanics and variational integrators[END_REF][START_REF] West | Variational integrators and the Newmark algorithm for conservative and dissipative mechanical systems[END_REF], the variational framework is extended to constrained (but smooth) systems, and even to systems with friction. According to the authors, the symplectic feature is a key property for accurate long time simulations. It ensures that the discrete system behaves close to a continuous physical one.

Time-integration in non-smooth dynamics

Enforcement of contact constraints in dynamics

For enforcing the contact constraints on a discrete system, three methods are mainly used: the penalty, the Lagrange multipliers and the augmented Lagrangian. Initially used in problem minimization, they have been adapted to structural mechanics problems with contact constraints. For example, [START_REF] Popp | Mortar Methods for Computational Contact Mechanics and General Interface Problems[END_REF] and [START_REF] Bussetta | Modelling and resolution of the mechanical contact problem and application in multi-physical context[END_REF] present these techniques in the framework of contact mechanic.

In brief, the penalty consists to apply a contact action proportionally to the penetration. The Lagrange multipliers method adds degrees of freedom to the system. They correspond to the contact actions required to enforce the contact constraints. The augmented Lagrangian joins these two methods: introduced for contact mechanics in [START_REF] Simo | An augmented lagrangian treatment of contact problems involving friction[END_REF], it allows a faster convergence for solving non-linear contact problems. A more recent approach is based on Nitsche's method, for a recent review see [START_REF] Chouly | An Overview of Recent Results on Nitsche's Method for Contact Problems in Geometrically Unfitted Finite Element Methods and Applications[END_REF].

Penalty and Lagrange multipliers have been early mixed with Newmark's time-integrators for addressing contact dynamics. The penalty method decreases the stable time-step, and they are then more suitable to unconditionally stable time-integrators as the AA method. For explicit time-integrators, it is more problematic to decrease the stable time-step. If some explicit schemes use nevertheless the penalty method, the Lagrange multipliers are preferred. Indeed they keep constant the stable time-step. For example, the pinball algorithm [START_REF]Contact-impact by the pinball algorithm with penalty and Lagrangian methods[END_REF] is based on the central difference method written in velocity for the time-integration. As the velocity is the main unknown, the contact is enforced thanks to the persistency condition. An other early explicit time-integrator with Lagrange multipliers is the Carpenter scheme [START_REF] Carpenter | Lagrange constraints for transient finite element surface contact[END_REF]. It is also based on the central difference method, but in the displacement formulation. The contact enforcement is done thanks to the Signorini's conditions on gap. Both algorithms present good results with a stable contact and a good energy behaviour.

The Lagrange multipliers within the implicit Newmark's time-integrators [START_REF] Carpenter | Lagrange constraints for transient finite element surface contact[END_REF] present in contrast a poor stability. For a rigid-deformable contact, the contact stresses largely oscillates [START_REF] Doyen | Time-Integration Schemes for the Finite Element Dynamic Signorini Problem[END_REF]. And for a rigid-rigid contact, the scheme is even unstable [START_REF] Di Stasio | Benchmark cases for robust explicit time integrators in non-smooth transient dynamics[END_REF].

The non-smooth dynamics formalism

A contact enforcement by Lagrange multipliers directly added to the Newmark's scheme forms a not well-defined framework. The Newmark's schemes approximate the integration of a second-order hyperbolic equation. The continuous solution is then a C 2 function regarding to time, with the acceleration defined everywhere. Two types of contact dynamical problems have to be distinguished:

1. a smooth contact problem, where the contact is active at the initial time and until the final one;

2. a non-smooth contact problem, where the contact becomes active and inactive during the simulation.

If the solutions of smooth contact problems are C 2 , the non-smooth contact problems present two type of non-smooth events where the velocity is discontinuous:

• impact, when a degree of freedom comes in contact;

• release, when it leaves the contact.

A clear example of discontinuity in velocity is the impacting bar problem described in [START_REF] Doyen | Time-Integration Schemes for the Finite Element Dynamic Signorini Problem[END_REF]: a 1D bar impacts a rigid boundary at one of its extremities.

The velocity of the contacting end presents to jumps: one at impact, where the velocity goes from a positive value to zero; one at release, from zero to a positive value. In order to address these discontinuities, Moreau in [START_REF] Moreau | Nonsmooth Mechanics and Applications[END_REF] introduces a formalism for impacts based on differential inclusions. This model describes the time-continuous dynamics by measures on velocity and contact actions. Thanks to it, robust time-integrators are developed for non-smooth dynamics.

The Moreau-Jean's family

In [START_REF] Moreau | Numerical aspects of the sweeping process[END_REF] and [START_REF] Jean | The non-smooth contact dynamics method[END_REF], Moreau and Jean apply the previous work [START_REF] Moreau | Nonsmooth Mechanics and Applications[END_REF] for designing a scheme, the non-smooth contact dynamics method. In the following, we refer to it by the Moreau-Jean scheme.

As the authors aim to design a scheme for rigid body dynamics, they add a contact law, the Newton law, to completely describe the problem around an impact:

v N (t + i ) = -e c v N (t - i ) (1.117)
with e c ∈ [0, 1] the restitution coefficient, and v N the relative velocity along the contact normal defined by the eq. (1.46). Indeed the rigid body dynamics does not take into account the deformable response of bodies: impact and release happen at the same time. The contact law of the eq. (1.117) is then a substitute to control the energy dissipated at impact:

• for e c = 1, the impact is elastic meaning the kinetic energy is conserved along the normal direction of impact;

• for e c = 0, the the impact is plastic meaning all the normal kinetic energy is absorbed.

The formulation of dynamics (1.84), used in the non-smooth contact dynamics framework, has the velocity for main unknown. The Signorini's conditions are then written on velocity with the eq. ( 1.52) and augmented by the contact law of the eq. ( 1.117):

∀(x , t) ∈ γ C × [t 0 , t f ],    If g N > 0, λ N = 0 Else g N = 0, 0 v N (t + ) + e c v N (t -) ⊥ λ N (t) 0 (1.118)
The notations are introduced in section 1.2. If the point is not in contact (g N > 0), no contact action is applied on it (λ N = 0). But if the contact is active (g N = 0), the normal contact action λ N enforces the condition v(t + ) = -e c v(t -). If e c = 0 the Moreau-Jean's contact conditions (1.118) matches the persistency condition (1.52). In a way, the persistency condition (strictly equivalent to Signorini's conditions) is a specific case of the Moreau-Jean's conditions. In [START_REF] Jean | The non-smooth contact dynamics method[END_REF][START_REF] Moreau | Numerical aspects of the sweeping process[END_REF], the authors propose a time-discrete form of the eq. (1.118). Under the discrete form, these velocity contact conditions are no more exactly equivalent to Signorini's conditions: bodies slightly penetrate each others. In a way, the contact enforcement focuses on the persistency condition relaxing the impenetrability one. According to [START_REF] Jean | The non-smooth contact dynamics method[END_REF], this is crucial for contact stability. Together with the θ-method for time-integration, it forms the Moreau-Jean scheme. As it is written on velocity, the Coulomb law is easily introduced for tangential contact constraints. The numerical properties of the Moreau-Jean scheme are well suited to non-smooth dynamics: it presents a stable contact both on displacement, velocity or contact actions; and a good energy behaviour [START_REF] Chen | A nonsmooth generalizedα scheme for flexible multibody systems with unilateral constraints[END_REF][START_REF] Di Stasio | Benchmark cases for robust explicit time integrators in non-smooth transient dynamics[END_REF].

Following the work of Moreau and Jean, a major contribution in nonsmooth dynamics is the Paoli-Schatzman scheme introduced in [START_REF] Paoli | A Numerical Scheme for Impact Problems I: The One-dimensional Case[END_REF] for its 1D formulation, and [START_REF] Paoli | A Numerical Scheme for Impact Problems II: The Multidimensional Case[END_REF] for a general formulation. This explicit scheme includes a contact law, like the Moreau-Jean scheme, but written on position. Indeed the time-integrator is the central difference method under its displacement formulation. The major contribution of [START_REF] Paoli | A Numerical Scheme for Impact Problems I: The One-dimensional Case[END_REF][START_REF] Paoli | A Numerical Scheme for Impact Problems II: The Multidimensional Case[END_REF] is an existence result: the scheme converges to a solution. The convergence is proven for a conserved tangential velocity around the impact, meaning that the scheme is not compatible with friction. And it seems touchy to adapt the Coulomb's law: the updated velocity is not accessible during the time-step, but a timestep after due to the CD method on displacement. This is a main drawback for industrial simulations where friction plays a role.

For implicit schemes, Chen et al. in [START_REF] Chen | A nonsmooth generalizedα scheme for flexible multibody systems with unilateral constraints[END_REF] and Schindler & Acary in [71] integrate the Moreau-Jean's contact law in implicit time-integrators. [START_REF] Chen | A nonsmooth generalizedα scheme for flexible multibody systems with unilateral constraints[END_REF] presents a scheme based on α-generalized time-integrator [START_REF] Chung | A Time Integration Algorithm for Structural Dynamics With Improved Numerical Dissipation: The Generalizedα Method[END_REF] and contact enforcement at velocity level with the Moreau-Jean's contact conditions of eq. (1.118). It gathers both the adjustable dissipation of α-generalized timeintegrator and the contact stability of the contact conditions. Note that the Newmark's schemes are contained in this formalism as a specific case of α-generalized scheme. [START_REF] Schindler | Timestepping schemes for nonsmooth dynamics based on discontinuous Galerkin methods: Definition and outlook[END_REF] follows the same philosophy but with discontinuous Galerkin methods [START_REF] Hulbert | Time finite element methods for structural dynamics[END_REF][START_REF] Li | Structural dynamic analysis by a time-discontinuous Galerkin finite element method[END_REF]. It adds then to the Moreau-Jean framework the high order of convergence of Galerkin methods. The main idea for both articles are to split the discrete system between an "impulsive part" which contains the contact terms, and a "non-impulsive" one. The impulsive part concentrates the non-smoothness of the system linked to impact. And then the non-impulsive part can be treated with the time-integrators designed for smooth dynamic. For all schemes based on the contact conditions of Moreau-Jean, the impenetrability condition is only satisfied at convergence. Otherwise the contact presents a residual penetration. A first attempt is made in [START_REF] Acary | Projected event-capturing time-stepping schemes for nonsmooth mechanical systems with unilateral contact and Coulomb's friction[END_REF] to exactly solve the impenetrability and the persistency conditions in the Moreau-Jean scheme. This work is extended in [START_REF] Brüls | Simultaneous enforcement of constraints at position and velocity levels in the nonsmooth generalized-α scheme[END_REF] for the non-smooth generalized-α scheme.

An other explicit time-integrator developed in NSCD framework is the CD-Lagrange scheme. Introduced in [START_REF] Fekak | A new heterogeneous asynchronous explicit-implicit time integrator for nonsmooth dynamics[END_REF], it gathers the central difference method written in its velocity form for the time-integration, and the Moreau-Jean's contact conditions. Contrary to the Paoli-Schatzman scheme, it is fully compatible with a tangential law as Coulomb friction. Its numerical results make it well suited to non-linear problems with non-smooth events as shown in [START_REF] Di Stasio | Benchmark cases for robust explicit time integrators in non-smooth transient dynamics[END_REF] for both rigid or deformable body dynamics. The contact is stable thanks to the use of Moreau-Jean's law. And it keeps the symplectic properties of CD method: the angular momentum is conserved to its exact initial value during the simulation, and the energy balance is conservative for elastic impacts. Moreover for conforming meshes at contact, it is fully explicit even for the contact problem. A key contribution of [START_REF] Fekak | A new heterogeneous asynchronous explicit-implicit time integrator for nonsmooth dynamics[END_REF] is to build the CD-Lagrange scheme both in the Moreau-Jean's formalism and in the variational one for non-smooth dynamics from [START_REF] Fetecau | Nonsmooth Lagrangian Mechanics and Variational Collision Integrators[END_REF].

Remark 3. The surprising numerical performances of the pinball algorithm [START_REF]Contact-impact by the pinball algorithm with penalty and Lagrangian methods[END_REF] and the Carpenter's scheme [START_REF] Carpenter | Lagrange constraints for transient finite element surface contact[END_REF] are easily explained in regard of non-smooth contact dynamics. Indeed the Carpenter's scheme is a specific case of the Paoli-Schatzman's scheme for e c = 0, and the Pinball algorithm is the CD-Lagrange scheme with e c = 0.

Variational time-integrators for non-smooth dynamics

Following the work of Simo & Tarnow [START_REF] Simo | The discrete energy-momentum method. Conserving algorithms for nonlinear elastodynamics[END_REF] on energy-momentum conserving algorithms, Laursen & Chawla in [START_REF] Laursen | Design of energy conserving algorithms for frictionless dynamic contact problems[END_REF] extends the scheme to unilateral contact problems on deformable bodies. It states a discrete persistency condition which ensures the energy conservation. The time-integration is done by a θ-method with θ = 1/2. For computing the contact action, three formulation are possible based on the penalty, the Lagrange multipliers or the augmented Lagrangian methods. As the Moreau-Jean scheme is based too on a θ-method, the two schemes differ only at the contact conditions: both written on velocity, but not at the same discrete times. If the energy is conserved, the contact actions exhibit bounded oscillations [see 28]. The contact condition on velocity leads to a slight interpenetration between contacting bodies similarly to the Moreau-Jean scheme. This drawback is overcomed in [START_REF] Laursen | Improved implicit integrators for transient impact problems?geometric admissibility within the conserving framework[END_REF] by adding a velocity correction which ensures the impenetrability condition. But this velocity correction involves an extra problem to solve under the constraints of energy conservation, and impenetrability condition. Hauret & Le Tallec in [START_REF] Hauret | Energy-controlling time integration methods for nonlinear elastodynamics and low-velocity impact[END_REF] propose a further improvement. A penalty formulation naturally ensures the impenetrability condition, and an energy-controlling integration of the inertial term provides a way to cut the unresolved high frequencies for keeping stability.

The work of Fetecau et al. [START_REF] Fetecau | Nonsmooth Lagrangian Mechanics and Variational Collision Integrators[END_REF] has also to be mentioned. Following [START_REF] Kane | Symplectic-energy-momentum preserving variational integrators[END_REF][START_REF] West | Variational integrators and the Newmark algorithm for conservative and dissipative mechanical systems[END_REF], it proposes a variational formalism for non-smooth dynamics both for a time-continuous and a time-discrete approach. The resulting schemes are then symplectic. In a time-continuous variational approach, Cirak & West proposes in [START_REF] Cirak | Decomposition contact response (DCR) for explicit finite element dynamics[END_REF] an explicit algorithm: the decomposition contact response. It integrates explicitly the dynamics by a predictor-corrector approach. Conservation of momenta and impenetrability condition are separately enforced. A main result is that the conservation of normal kinetic energy around impact is conditioned to an "elastic" impact with a restitution coefficient equal to one.

Leine et al. in [START_REF] Leine | Hamilton's Principle as Variational Inequality for Mechanical Systems with Impact[END_REF] contributes to formalize the Lagrangian of a nonsmooth problem. They propose two ways of writing the Hamiltonian of a system submitted to an inequality constraint. The first form called the "strong Hamiltonian principle" is the one previously used by Fetecau et al. and Cirak & West in [START_REF] Fetecau | Nonsmooth Lagrangian Mechanics and Variational Collision Integrators[END_REF] and [START_REF] Cirak | Decomposition contact response (DCR) for explicit finite element dynamics[END_REF]. The impact time is explicitly written in the time-step. In a one particle example, it naturally leads to a discretein-time form with an incorporated impact law. The restitution coefficient must be equal to one for observing energy conservation. The second form, named "weak Hamiltonian principle", does not required to write explicitly the impact time. The related impulsive part is included by using measures. On a one-particle example, the discrete-in-time form does not include an impact law. An other interesting contribution on discrete Lagrangian with non-smooth inequality constraints is made by Kaufman & Pai in [START_REF] Kaufman | Geometric Numerical Integration of Inequality Constrained, Nonsmooth Hamiltonian Systems[END_REF]. After a review on existing "geometrical" methods for integrating non-smooth Hamiltonian systems, they propose a discrete Lagrangian including nonsmooth inequality constraints and an algorithm in order to solve it. For a variational approach discrete both in space and time, a recent contribution is made by Demoures et al. in [START_REF] Demoures | A multisymplectic integrator for elastodynamic frictionless impact problems[END_REF] with an explicit time-integrator for elastodynamic problems with unilateral contact in 1D. The authors underline that two non-smooth events happen for a deformable contact: impact and release. Their scheme integrates a modified Lagrangian with a null mass for the contact node.

Time-integrators for NSCD with unified notations

In this section, time-integrators for non-smooth dynamics are given under unified notations in order to compare them easily.

The considered problem is described by the fig. 1.4. A meshed deformable body is impacting along Γ C , a motionless rigid boundary Γ R . The discretein-space dynamics (1.84) of this system is also valid for a rigid body problem with several bodies, or even for a set of moving particles. The crucial difference between these problems relies on the integration of internal forces F int . But this aspect is not addressed here. The internal forces obey to a linear elastic model:

F int (U) = KU (1.119)
Note that the integration of internal forces is crucial for time-integration, see for example [START_REF] Hauret | Energy-controlling time integration methods for nonlinear elastodynamics and low-velocity impact[END_REF][START_REF] Simo | The discrete energy-momentum method. Conserving algorithms for nonlinear elastodynamics[END_REF].

In the discrete in space dynamics (1.84), the contact enforcement is done by the Lagrange multipliers r. For the Signorini's conditions, either the displacement or velocity forms are used depending on the time-integrator. The discrete in space dynamics together with Signorini's conditions written in displacement is described by the eq. (1.86); and the discrete-in-space dynamics with Signorini's conditions on velocity by the eq. (1.87). Only a unilateral contact is considered, meaning the contact action involves no friction. If some schemes below are fully compatible with a tangential contact, other ones are not. Considering only a unilateral contact allows to unify the notations.

At contact boundary γ C , the relative normal velocity (1.46) is denoted v N,n+1 and got by multiplying the global velocity Un+1 by the operator L defined in the eq. (1.60). This operator for a contact on a rigid and motionless boundary Γ R is simply a projection on n the contact normal on γ C (see fig. • a time-grid for integer discrete times t n+1 ,

{t n } ≡ t n = t 0 + nh | n ∈ N 0 , n N with N | t f = N h (1.120)
• a time-grid for midpoints or half-steps discrete times

t n+ 1 2 , {t n+ 1 2 } ≡ t n+ 1 2 = t n + 1 2 h | t n ∈ {t n } (1.121)
The time-step h is considered here constant over the simulation. Note that the formulations given below for time-integrators are not valid for a variable time-step [see 9, 59, 83].

The Moreau-Jean scheme

The Moreau-Jean scheme introduced in [START_REF] Jean | The non-smooth contact dynamics method[END_REF][START_REF] Moreau | Numerical aspects of the sweeping process[END_REF] is based on a theta method for time-integration and the Moreau-Jean's contact conditions of eq. (1.118).

The discrete dynamics together with the discrete contact conditions and the time-integration relation are:

U n+1 = U n + h (1 -θ) Un + θ Un+1 (1.122) M( Un+1 -Un ) = h(F ext n+θ -F int n+θ ) + L t n+α r N,n+1 (1.123) ∀i ∈ γ C , (1.124) If (g n+α ) i > 0, (r N,n+1 ) i = 0 Else (g n+α ) i 0, 0 (v N,n+1 ) i + e c (v N,n ) i ⊥ (r N,n+1 ) i 0
With:

• U * n+α = U n + hα Un , a predictor on displacement; • L n+α = L(U * n+α ), g n+α = g(U * n+α )
the projection against the contact normal and the gap, build on the predicted displacement;

• hF ext n+θ = hθF ext (t n+1 ) + h(1 -θ)F ext (t n ), time-integration of external forces; • hF int n+θ = hθF int (U n+1 )+h(1-θ)F int (U n )
, time-integration of internal forces;

• θ ∈ [0, 1], which set the time-integration interpolation.

Each value of θ leads to a time-integrator. The most common are: θ = 1 for the implicit Euler method; or θ = 1/2 for the mid-point rule method.

A major feature of the Moreau-Jean scheme is that the contact is build on the deformed configuration defined by U * n+α , a prediction of the displacement. This deformed configuration is not updated during the time-step, and then both g n+α and L n+α do not depend on U n+1 .

The eq. (1.123) is the discrete dynamics. It comes from the integration of the eq. ( 1.84) by a θ-method between t n and t n+1 . The impulse r n+1 then gathers:

• a non-smooth part, which ensures the jump in velocity at impact

M( Un+1 -Un ) = L t n+α r ns N,n+1
• a smooth part, which is the integration of contact force during a persistent contact over the time-step

0 = h(F ext n+θ -F int n+θ ) + L t n+α r s N,n+1
More details about the splitting between smooth and non-smooth part in such a scheme can be found in [START_REF] Chen | A nonsmooth generalizedα scheme for flexible multibody systems with unilateral constraints[END_REF] or [START_REF] Fekak | A new heterogeneous asynchronous explicit-implicit time integrator for nonsmooth dynamics[END_REF]. The discrete dynamic (1.123) and the time-integration relation of eq. (1.122) form a solvable system on velocity Un+1 . This system is implicit as the internal forces depends on U n+1 which depends itself on Un+1 . With a linear elastic model, the system on Un+1 is:

M Un+1 = R n + L t n+α r N,n+1 (1.125)
With :

• M = M + h 2 θ 2 K , the problem matrix; • R n = M -h 2 θ(1 -θ)K Un -hKU n + hF ext n+θ ,

which depends only of known quantities.

If the contact actions are not considered, the eq. (1.125) defines the free velocity:

Ufree n+1 = M-1 R n (1.126)
In order to get r n+1 , the dynamics (1.125) is condensed on γ C by a leftmultiplication by L n+1 M-1 :

Hr N,n+1 = v N,n+1 -v free N,n+1 (1.127)
With:

• H = L n+α M-1 L t n+α
, the Delassus operator associated to the contact problem;

• v free N,n+1 , v N,n+1 the relative normal free-velocity and velocity. The eq. (1.127) together with the contact conditions of eq. (1.118) form the contact problem: The discrete contact conditions in the eq. ( 1.118) relaxes the impenetrability condition. Indeed the contact is detected by a negative value for the predicted gap g n+α . But if the contact impulse ensures the persistency condition, the update of displacement does not depend on the penetration and then does not ensure the impenetrability condition. As a consequence, the bodies can be interpenetrated when an impact occurs. This penetration stays constant during the contact because of the persistency condition which bans further penetration. And it tends to disappear with a space-time convergence.

Hr N,n+1 = v N,n+1 -v free N,n+1 (1.128) ∀i ∈ γ C , if (g n+α ) i > 0, (r n+1 ) i = 0 if (g n+α ) i 0, 0 (r n+1 ) i ⊥ (v N,n+1 + e c v N,n ) i 0 v N,
Remark 4. Between the continuous form of eq. (1.87) and the discrete one of eq. (1.124), the admissible gap for active contact phases is relaxed from a strictly null value g(U) = 0 to a null or negative value g n+α 0. This discrete approximation of the persistency conditions makes feasible the timeintegration of non-smooth events in a weak sense (see remark 5).

Following [START_REF] Jean | The non-smooth contact dynamics method[END_REF], the operator L n+α = L(U * n+α ) is build on the predicted deformed configuration U * n+α . According to Jean the prediction does not influence a lot the discrete solution. He proposes three possible values for α: α = 0, α = 1/2, α = 1; and recommends α = 0 for θ = 1. But the choice of α can maybe be done by following the equality:

g n+α (U * n+α ) = g n+α (U n ) + hαv N,n (1.129) and g n+α (U n+1 ) = g n+α (U n ) + h(1 -θ)v N,n + hθv N,n+1 (1.130)
Here g n+α names the gap operator on the time-step [t n , t n+1 ], computed on the deformed configuration defined by U * n+α . This operator is then linear on U.

If contact happens at t n+1 , v N,n+1 is null and then

g n+α (U * n+α ) = g n+α (U n+1 ) if α = 1 -θ (1.131)
With this choice, the predicted gap g(U * n+α ) is then the correct one for the node in contact at t n+1 . Moreover it matches the "approximate gap" used by Jean in [START_REF] Jean | The non-smooth contact dynamics method[END_REF].

Algorithm 1 Algorithm for Moreau-Jean scheme

Inputs : U n , Un 1: Ufree n+1 ← M-1 R n Free velocity 2: U * n+α ← U n + hα Un Predicted displacement 3: g n+α ← g(U * n+α ) and L n+α ← L(U * n+α ) Contact treatment 4: r N,n+1 by solving (1.128) 5: Un+1 ← Ufree n+1 + M-1 L t n+α r N,n+1
Velocity correction

6: U n+1 = U n + h (1 -θ) Un + θ Un+1 Displacement update Outputs : U n+1 , Un+1
The algorithm 1 concerns only the case with internal stresses linear on displacement. If the internal stresses depends on velocity or are non-linear on displacement, the steps one to five must be enclosed in a non-linear algorithm as a Newton-Raphson one. It behaves as a predictor-corrector on velocity: the total velocity is the sum of the free-velocity and a velocity correction which enforces the contact constraints. A crucial point in this algorithm is that the contact problem is linear. With the predicted displacement, the contact problem is independent of updated displacement and then it is linear on velocity.

The energy balance for this scheme is established in [START_REF] Acary | Energy conservation and dissipation properties of time-integration methods for nonsmooth elastodynamics with contact[END_REF]. For a system with linear elastic internal forces and submitted to external ones:

1 2 Ut n+1 M Un+1 + 1 2 U t n+1 KU n+1 t n+1 tn = h Ut n+θ F ext n+θ + Ut n+θ L t n+α r N,n+1 + 1 2 -θ Un+1 -Un 2 M + U n+1 -U n 2 K (1.132)
With:

• Un+1 -Un 2 M = ( Un+1 -Un ) t M( Un+1 -Un ) • U n+1 -U n 2 K = (U n+1 -U n ) t K(U n+1 -U n ) • Un+θ = θ Un+1 + (1 -θ) Un
The energy balance (1.132) is conservative under two numerical conditions (with null external forces):

• θ = 1 2 for cancelling the last term on the right side, which corresponds to a numerical dissipation term;

• a contact work Ut n+θ L t n+α r N,n+1 equal to zero. The contact term with θ = 1/2 is equal to:

Ut n+ 1 2 L t n+α r N,n+1 = v t N,n+ 1 2 r N,n+1 = 1 2 (1 -e c )v t N,n r N,n+1 (1.133)
The Moreau-Jean scheme is then conservative for θ = 1 2 and e c = 1. If the role of numerical parameter θ is expected, the role of the restitution coefficient e c is more surprising. Indeed e c is a model parameter, set by the physic of the problem. For example in a problem involving deformable bodies, e c is set to zero which implies an energy loss at impact.

The Moreau-Jean scheme can deal with tangential contact constraints. In [START_REF] Jean | The non-smooth contact dynamics method[END_REF], the Coulomb's law is integrated in the scheme for both deformable and rigid bodies.

In the following, the scheme is set with θ = 1 2 , α = 1 2 : the scheme is conservative, and the predicted gap is the exact one in case of contact.

Remark 5. The Moreau-Jean scheme comes from an integration of the dynamics (1.84) between t n and t n+1 :

t n+1 tn Md U = t n+1 tn F ext -F int (U) dt + t n+1 tn L t dr (1.134) ⇔ M( Un+1 -Un ) = h(F ext n+θ -F int n+θ ) + L t n+α r N,n+1 (1.135)
The measures on velocity and contact actions allow a time-integration in a weak sense of non-smooth events. Whatever their locations in the timeinterval or their number, they are gathered with the smooth phases in integral quantities:

t n+1 tn Md U = M t i,k t - i,k+1 t + i,k Üdt + U(t + i,k+1 ) -U(t - i,k+1 ) ⇔ t n+1 tn Md U = M(U n+1 -U n ) (1.136) t n+1 tn dr = t i,k t - i,k+1 t + i,k λdt + λ t i δt i ⇔ t n+1 tn dr = r N,n+1 (1.137)
with {t i,k } which gathers the ends of time-interval [t n , t n+1 ] and the nonsmooth times inside. The quantities M(U n+1 -U n ) and r n+1 handle both the smooth and non-smooth dynamics. This makes the Moreau-Jean scheme very robust in case of a large number of non-smooth events: the determination of the non-smooth times is not required. The Paoli-Schatzman scheme [START_REF] Paoli | A Numerical Scheme for Impact Problems I: The One-dimensional Case[END_REF][START_REF] Paoli | A Numerical Scheme for Impact Problems II: The Multidimensional Case[END_REF] and the CD-Lagrange scheme [START_REF] Fekak | A new heterogeneous asynchronous explicit-implicit time integrator for nonsmooth dynamics[END_REF] are similarly built, and present then the same property.

The Paoli-Schatzman scheme

Introduced in [START_REF] Paoli | A Numerical Scheme for Impact Problems I: The One-dimensional Case[END_REF] and [START_REF] Paoli | A Numerical Scheme for Impact Problems II: The Multidimensional Case[END_REF] by Paoli & Schatzman this scheme integrates the dynamics with the central difference method written on displacement. It enforces the Signorini's conditions on displacement while integrating the Newton's law. The equations of the scheme are:

1 h 2 M(X n+1 -2X n + X n-1 ) = F ext n -F int (X n ) + L t n λ n+1 (1.138) Un = 1 2h (X n+1 -X n-1 ) (1.139) ∀i ∈ γ C , 0 (g n ) i ⊥ (λ n+1 ) i 0 (1.140)
With:

g n = g X n+1 + e c X n-1 1 + e c (1.141) L n = L(X * n+1 ) (1.142) X * n+1 = h 2 M -1 F ext n -F int (X n ) + 2X n -X n-1 (1.143)
The discrete dynamics (1.138) comes from the dynamics (1.84) integrated by the central difference method. The mass matrix is lumped making it diagonal. The eq. (1.139) is the time-integration relation. The velocity in this scheme is a post-treated quantity. Indeed Un depends on X n+1 , and can then be computed only after the time-step t n to t n+1 . This property makes difficult the integration of the Coulomb's law in the scheme. Indeed the tangential contact force depends on velocity. The contact conditions of eq. (1.140) are quite unusual. They involve:

• X * n+1 , a predictor on position which is the updated position without any contact, it is used to build the normal projection operator of eq. (1.142);

• g n , an estimation of the gap with an approximate Newton's law by a mean on positions involving e c .

For a non-null restitution coefficient, these contact conditions of eq. (1.140) ensure approximately the Newton's law over three time-steps:

If g n+1 -g n-1 = 0, ⇔ L n+1 (X n+2 + e c X n ) -L n-1 (X n + e c X n-2 ) = 0 ⇔ L n+1 (X n+2 -X n ) + e c L n+1 (X n -X n-2 ) = 0 (L n-1 = L n+1 ) ⇔ v N,n+1 + e c v N,n-1 = 0
The Newton's law is ensured exactly between v N,n+1 and v N,n-1 if L n-1 = L n+1 , which means that the normal contact directions are equal. For an impact on a rigid motionless boundary, this is true.

As shown in [START_REF] Doyen | Time-Integration Schemes for the Finite Element Dynamic Signorini Problem[END_REF], both the persistency and the impenetrability conditions are met for e c = 0. Indeed the gap estimation is exactly the geometrical gap as e c = 0, which ensures the impenetrability. Moreover during a contact phase, the impenetrability condition over three time-steps ensures a null acceleration along the contact normal. And then the persistency condition is ensured as well. In [START_REF] Doyen | Time-Integration Schemes for the Finite Element Dynamic Signorini Problem[END_REF], the given energy balance is similar to the one of the central difference method of eq. (1.110):

1 2 Ut n+1 M Un+1 + 1 2 U t n+1 KU n+1 - 1 8 h 2 Üt n+1 M Ün+1 t n+1 tn = (U n+1 -U n ) t 1 2 (L t n+1 λ n+2 + L t n λ n+1 ) (1.144)
The energy balance (1.144) is valid for a conservative system with unilateral contact. But there is no general result about the sign of contact work (the term on the right side). For e c = 0, energy can be injected at impact.

Nevertheless if e c = 0, the contact work is dissipative [START_REF] Doyen | Time-Integration Schemes for the Finite Element Dynamic Signorini Problem[END_REF].

The dynamics (1.138) can be rewritten thanks to

X * n+1 : M(X n+1 -X * n+1 ) = h 2 L t n λ n+1 (1.145)
This equation can be projected on Γ C by a multiplication at left by

L n M -1 .
Together with the contact conditions of eq. (1.140), they form the contact problem:

L n M -1 L n λ n+1 = L n (X n+1 -X * n+1 ) (1.146) ∀i ∈ γ C , 0 (g n ) i ⊥ (λ n+1 ) i 0 Algorithm 2 Algorithm for Paoli-Schatzman scheme Inputs : X n , X n-1 1: X * n+1 ← h 2 M -1 (F ext n -F int (U n )) + 2X n -X n-1 2: g n ← g X n+1 +ecX n-1 1+ec and L n = L(X * n+1 ) 3: λ n+1 by solving (1.146) 4: X n+1 ← X * n+1 + h 2 M -1 L t n λ n+1 Outputs : X n+1
The dynamics (1.145) expressed with the predictor on position leads to the predictor-corrector algorithm 2. This algorithm is explicit if the mass matrix is lumped and the contact problem is explicit. The velocity does not appear as it is a post-treated quantity. The eq. (1.139) gives the formula to get Un according to X n+1 and X n-1 . Here L and the gap operator are build on the predicted displacement X * n+1 . They are then linear on U on the time-step.

Remark For e c = 0, the Paoli-Schatzman scheme is equivalent to the Carpenter scheme [START_REF] Carpenter | Lagrange constraints for transient finite element surface contact[END_REF].

The Laursen-Chawla algorithm

Contrary to the algorithms above, this one concerns only the deformable bodies. It does not integrate the Newton's impact law. Introduced in [START_REF] Laursen | Design of energy conserving algorithms for frictionless dynamic contact problems[END_REF], its major advantage is the energy conservation for both integration of internal forces and contact terms. Laursen & Chawla follows indeed [START_REF] Simo | The discrete energy-momentum method. Conserving algorithms for nonlinear elastodynamics[END_REF] which deals with the integration of internal forces at constant energy. Here the focus is only on the contact conditions which ensure the energy conservation at impact. The internal forces are considered only linear.

The scheme equations are:

M( Un+1 -Un ) = hF n+ 1 2 + L t n+ 1 2 r n+ 1 2
(1.147)

U n+1 = U n + h Un+ 1 2 (1.148) ∀i ∈ γ C , If (g n ) i = 0, (r n+ 1
2

) i = 0 If (g n ) i 0, 0 (r n+ 1 2 ) i ⊥ (v N,n+ 1
2

) i 0 (1.149)
With:

Un+ 1 2 = 1 2 ( Un+1 + Un ) v N,n+ 1 2 = L n+ 1 2 Un+ 1 2 (1.150) U n+ 1 2 = 1 2 (U n+1 + U n ) F n+ 1 2 = - 1 2 K(U n+1 + U n ) (1.151) g n = g(U n ) L n+ 1 2 = L(U n+ 1 2 ) (1.152)
The discrete dynamics (1.147) in term of Un+ 1 2 is:

M + 1 4 h 2 K M Un+ 1 2 = M Un - 1 2 hKU n Rn + 1 2 L t n+ 1 2 r n+ 1 2 (1.153)
Condensed on γ C and together with the contact conditions of eq. (1.149), the eq. (1.153) forms the contact problem:

1 2 L n+ 1 2 M-1 L t n+ 1 2 r n+ 1 2 = v N,n+ 1 2 -v free N,n+ 1 2 (1.154) ∀i ∈ γ C , If (g n ) i = 0, (r n+ 1
2

) i = 0 If (g n ) i 0, 0 (r n+ 1 2 ) i ⊥ (v N,n+ 1
2

) i 0 With: Ufree N,n+ 1 2 = M-1 R n v free N,n+ 1 2 = L n+ 1 2 Ufree N,n+ 1 2 (1.155)
This scheme with linear internal forces is quite similar to the Moreau-Jean scheme with θ = gives:

1 2 Ut n+1 M Un+1 + 1 2 U t n+1 KU n+1 n+1 n = v t N,n+ 1 2 r n+ 1 2 (1.156)
The nullity of contact work (the term on right side) comes directly from the complementary condition of eq. (1.149).

For the algorithm, the main steps are similar to the Moreau-Jean scheme algorithm 1.

The CD-Lagrange scheme

Introduced by Fekak et al. in [START_REF] Fekak | A new heterogeneous asynchronous explicit-implicit time integrator for nonsmooth dynamics[END_REF], this explicit scheme is based on the central difference method written in velocity for time-integration and the contact conditions of Moreau-Jean:

U n+1 = U n + h Un+ 1 2 (1.157) M Un+ 3 2 -Un+ 1 2 = h F ext n+1 -F int (U n+1 , Un+ 1
2

) +L t n+1 r n+ 3 2 (1.158) ∀i ∈ γ C ,    If (g n+1 ) i > 0, (r n+ 3
2

) i = 0, If (g n+1 ) i 0, 0 (r n+ 3
2

) i ⊥ v N,n+ 3 2 + e c v N,n+ 1 2 i 0 (1.159)
With:

g n+1 = g(U n+1 ) L n+1 = L(U n+1 ) (1.160)
The operators of eq. ( 1.160) are build on the deformed configuration defined by the explcit displacement U n+1 . The contact normal n is computed on this deformed configuration, and it defines then the gap operator g and the projection operator L. On a time-step, these operator are not modified as the displacement U n+1 . Contrary to the preceding time-integrators, where g and L are computed on a deformed configuration defined by a predictor, here there are build on the exact discrete displacement U n+1 . The discrete dynamics (1.158) comes from the integration of the dynamics (1.84) between two midpoint times: t n+ 1 2 and t n+ 3

2

. As usual in central difference method, the mass matrix is lumped and the internal forces are evaluated with Un+ 1 2 . For clarity, the internal forces are refereed after as F int (U n+1 ). The contact conditions of eq. (1.159) are expressed on velocity. As for the Moreau-Jean scheme, a slight interpenetration appears at impact. It stays constant during the contact, and decreases when the space-time mesh is refined. The scheme can both deal with deformable and rigid bodies as the contact conditions integrate the Newton's law.

By introducing the free-velocity, Ufree

n+ 3 2 = Un+ 1 2 + hM -1 F ext -F int (U n+1 ) (1.161) the dynamics (1.158) is condensed on Γ C by a multiplication with L n+1 M -1 .
Together with the contact conditions of eq. ( 1.159), it forms the contact problem:

H n+1 r n+ 3 2 = v N,n+ 3 2 -v free N,n+ 3 2 (1.162) ∀i ∈ γ C ,    If (g n+1 ) i > 0, (r n+ 3
2

) i = 0, If (g n+1 ) i 0, 0 (r n+ 3
2

) i ⊥ v N,n+ 3 2 + e c v N,n+ 1 2 i 0 (1.163)
The Delassus operator H n+1 = L n+1 M -1 L t n+1 is diagonal if the contacting body is described by a rigid and motionless boundary. If two deformable bodies are involved, it is diagonal only for conforming meshes [START_REF] Fekak | Etude de la réponse dynamique des ponts roulants soumis à des chocs multiples pendant un séisme : Co-simulation implicite / explicite multi-échelle en temps pour la dynamique du contact[END_REF]. For a contact between a meshed deformable body and a rigid motionless boundary as described by fig. 1.4, the operators L and H are:

∀i ∈ γ C , L i 3i..3i+2
= n ix n iy n iz (1.164)

H ij =    n t i M -1 3i..3i+2 3i..3i+2 n i = M -1 i i = j 0 i = j (1.165) With, M 3i..3i+2 3i..3i+2 =    M -1 i 0 0 0 M -1 i 0 0 0 M -1 i    n i =    n ix n iy n iz   
As the mass matrix is lumped, the operator H is the inverse of mass matrix but for the nodes on Γ C . Its diagonal shape makes the solving of the LCP (1.162) straightforward (here for e c = 0):

∀i ∈ γ C , (r n+ 3
2

) i =    0 if (g n+1 ) i > 0 max 0, -(H ii ) -1 (v free N,n+ 3
2

) i if (g n+1 ) i 0 (1.166)
The energy balance comes directly from the energy balance of eq. (1.113) without impact:

1 2 Ut n+ 3 2 M Un+ 3 2 t n+ 3 2 t n+ 1 2 = 1 2 h( Un+ 3 2 + Un+ 1 2 ) t F ext (t n+1 ) -F int (U n+1 , Un+ 1 2 ) + 1 2 ( Un+ 3 2 + Un+ 1 2 ) t L t N,n+1 r N,n+ 3 2 (1.167)
The contact work is: If the contact involves a deformable solid, e c is set to zero to match the persistency condition. ∆W C,N is then dissipative but only at impact. During contact and at release, the contact does not dissipate any energy: As for the Moreau-Jean scheme [START_REF] Jean | The non-smooth contact dynamics method[END_REF] or the Cirak and West algorithm [START_REF] Cirak | Decomposition contact response (DCR) for explicit finite element dynamics[END_REF], the energy conservation is conditioned to e c = 1 and then impossible for a contact involving a deformable body. Remark 6. The contact work at impact for e c = 0 is directly linked to the loss of kinetic energy when the normal velocity v n+ 3 2 is set to zero. Projecting the dynamics (1.158) on γ C with a left multiplication by LM -1 gives:

∆W C,N = 1 2 (v n+ 3 2 + v n+ 1 2 )r N,n+ 3 2 ⇔∆W C,N = 1 2 (1 -e c )v n+ 1 2 r N,n+ 3
1. at impact, v n+ 1 2 < 0 and r n+ 3 2 > 0, then ∆W C < 0;
(v n+ 3 2 -v n+ 1 2 ) = -hLM -1 F int (U n+1 ) + Hr n+ 3 2 (1.169) As M -1 is diagonal, LM -1 F int (U n+1
) gathers the internal forces f int i projected along n i multiplied by terms of diagonal mass M ii for each node i of γ C . H being diagonal, each line of the eq. (1.169) is independent and corresponds to one node i of γ C . By left multiplying each line by 1 2 (v i,n+ 3 2

+ v i,n+ 1
2

) t H -1
ii , the eq. (1.169) gives an energy balance for each node of γ C :

∀i ∈ γ C , 1 2 v i,n+ 3 2 H -1 ii v i,n+ 3 2 t n+ 3 2 t n+ 1 2 ∆E i,k = - 1 2 h(v i,n+ 3 2 + v i,n+ 1
2

) t (f int i • n i ) ∆W i,int + 1 2 h(v i,n+ 3 2 + v i,n+ 1
2

) t r i,n+ 3 2 ∆W i,C,N (1.170)
The energy balance (1.170) gives for each node of γ C a contact work ∆W i,C,N split in two components:

• ∆E i,k the difference of kinetic energy between t n+ 1 2 and t n+ 3 2 ;

• ∆W i,int the work of internal forces along the contact normal.

If the internal forces along n i do not produce any energy, ∆W i,C,N is directly the loss of kinetic energy due to the persistency condition for node i. The energy loss at impact is then strongly linked to the masses of contact nodes.

Algorithm 3 Algorithm for CD-Lagrange scheme with unilateral contact

Inputs : U n , Un+ 1 2 1: U n+1 ← U n + h Un+ 1 2 Displacement update 2: Ufree n+ 3 2 ← Un+ 1 2 + hM -1 F ext n+1 -F int (U n+1 )
Free-velocity update

3: g n+1 ← g(U n+1 ) and L n+1 ← L(U n+1 ) Contact treatment 4: r n+ 3 2
by solving (1.166)

5: Un+ 3 2 ← Un+ 3 2 + M -1 L t n+1 r n+ 3 2

Velocity correction

Outputs : U n+1 , Un+ 3 2

The algorithm 3 acts like a predictor-corrector on velocity. The displacement is computed once at the beginning of the time-step, and not corrected after. The solving of contact problem (1.162) is then linear, and even explicit if the Delassus operator H is diagonal. In this case, this algorithm is naturally parallel as all steps involve only local quantities.

The discrete Coulomb's law. The CD-Lagrange scheme is fully compatible with a Coulomb's law. We remind that the problem considered involves only a contact between a rigid and motionless boundary and a deformable body (see fig. 1.4). The contact impulse r n+ 3 2 is no more along the inner normal n of γ C . For each node of γ C ,the contact impulse is decomposed into a normal and a tangential part: . A new projection operator L T is defined along t, similarly to the definition of L (now named L N for clarity) along n. L N and L T are orthogonal: L N L t T = L T L t N = 0. The dynamics (1.172) contains now normal and tangential contact impulses:

r i,n+ 3 2 = r i,N,n+ 3 2 n i + r i,T,
M Un+ 3 2 -Un+ 1 2 = h F ext n+1 -F int (U n+1 , Un+ 1
2

) +L t N,n+1 r N,n+ 3 2 + L t T,n+1 r T,n+ 3 2
(1.172)

The dynamics (1.172) is condensed on γ C along t, and together with the discrete Coulomb's law it forms the tangential contact problem:

H T,n+1 r T,n+ 3 2 = v T,n+ 3 2 -v free T,n+ 3 2 (1.173) ∀i ∈ γ C , (1.174)
If (r N,n+ 3
2

) i = 0, (r T,n+ 3
2

) i = 0 Else (r N,n+ 3
2

) i > 0, 0 (µr N,n+ 3 2 -r T,n+ 3
2

) i ⊥ (v T,n+ 3
2

) i 0 (r T,n+ 3
2

) i t i = -α ufree i,n+ 3 2 -( ufree i,n+ 3 2 • n i )n i α, (r T,n+ 3
2

) i ∈ R + (1.175) With: v T,n+ 3 2 = -L T,n+1 Un+ 3 2 (1.176) L T,n+1 = L( Ufree n+ 3
2

) (1.177) H T,n+1 = L T,n+1 M -1 L t T,n+1 (1.178)
The crucial point in this discrete Coulomb's law is that the tangential direction is known thanks to Ufree

n+ 3 2 .
Remark 7. This comes directly from the dynamics (1.172): with the expression of eq. (1.179) (which corresponds to extract the tangential part on γ C ), it follows:

Un+ 3 2 = Ufree n+ 3 2 + M -1 L t N r N,n+ 3 2 + M -1 I T,n+ 3
( Un+ 3 2 -L t N v N,n+ 3 2
) = ( Ufree

n+ 3 2 -L t N v free N,n+ 3
2

) + M -1 L t N r N,n+ 3 2 -L t N H N r N,n+ 3 2 =0 + M -1 I T,n+ 3 2 ⇔ UT,n+ 3 2 = Ufree T,n+ 3 2 + M -1 I T,n+ 3 2 (1.180) UT,n+ 3 2
contains only the tangential part on γ C . As M is diagonal, the eq. (1.180) can be split between nodes of γ C :

∀i ∈ γ C , 1 M i r i,T,n+ 3 2 = ui,T,n+ 3 2 -ufree i,T,n+ 3 2 (1.181)
The continuous Coulomb's law (1.57) gives:

r i,
∀i ∈ gamma C , 1 M i + 1 α r i,T,n+ 3 2 = -ufree i,T,n+ 3 2
The tangential impulses r i,T,n+ 3 2 are then co-linear and opposite to the freevelocities ufree i,T,n+ ) i do not change during the solving. This makes the tangential problem linear. For a diagonal H, it is even explicit and solved similarly to the normal one of eq. (1.166):

∀i ∈ γ C , (1.182) (r T,n+ 3
2

) i =    0 if (r N,n+ 3
2

) i = 0 min µ(r N,n+ 3
2

) i , -(H T,ii ) -1 (v free T,n+ 3
2

) i if (r N,n+ 3
2

) i 0
The energy balance (1.167) is augmented with the tangential contact work to give:

1 2 Ut n+ 3 2 M Un+ 3 2 t n+ 3 2 t n+ 1 2 = 1 2 h( Un+ 3 2 + Un+ 1 2 ) t F ext (t n+1 ) -F int (U n+1 , Un+1 ) 1 2 ( Un+ 3 2 + Un+ 1 2 ) t L t N,n+1 r N,n+ 3 2 + L t T,n+1 r T,n+ 3 2 (1.183)
Algorithm 4 Algorithm for the CD-Lagrange scheme with a rigiddeformable contact and friction

1: U 1 2 ← U0 + 1 2 hM -1 F(U 0 ) Initialisation 2: for (t n+1 , t n+ 3 2 ) ∈ {t n+1 } × {t n+ 3 2 } do 3: U n+1 ← U n + h Un+ 1 2 Displacement update 4:
Ufree

n+ 3 2 ← Un+ 1 2 + hM -1 F ext n+1 -F int (U n+1 )
Free-vel. update 5:

g n+1 ← g(U n+1 ) and L N,n+1 ← L(U n+1 ) Normal contact 6:
r N,n+ 3 2 by solving (1.166)

7:
L T,n+1 ← L( Ufree

Un+ 3 2 ← Un+ 3 2 + M -1 L t N r N,n+ 3 2 + L t T r T,n+ 3 2
Contact corr.

10:

Un+ 3 2 |Γ D ← UD (t n+ 3 2 )

Dirichlet correction 11: end for

The algorithm 4 is slightly modified: only one step is added to the frictionless algorithm 3 for solving the tangential contact problem. For diagonal H N and H T , the algorithm 4 keeps the explicit and parallel features. Its efficiency for a rigid-deformable contact problem is then high compared to the usual implicit algorithms. Indeed their contact problems are non-linear which requires a non-linear solver [START_REF] Jourdan | A Gauss-Seidel like algorithm to solve frictional contact problems[END_REF]. In the algorithm 4, the Dirichlet's boundary conditions are enforced directly on velocity.

Conclusion: the time-integration for the impact dynamics

The description of the discrete-in-space smooth dynamics leads to a symplectic form, which keeps invariant the linear momentum, the angular momentum, and the energy. The fundamental equations of dynamics are an expression of the conservation of these quantities. A robust time-integrator preserves this symplectic feature for the discrete-in-time problem. For the standard Newmark's time-integrators, the central difference method is symplectic in the sense that it conserves the exact angular momentum and an approximate form of energy.

The contact constraints bring non-smoothness in the system which makes the time-integration more complex. A suitable framework is then the nonsmooth contact dynamics formalism from Moreau-Jean. By dealing in a weak sense with the non-smooth events, it leads to robust time-integrators. The Moreau-Jean formalism introduces also a form of Signorini's conditions on velocity, the persistency condition. On its discrete form, it relaxes the impenetrability condition which ensures an high contact stability. The CD-Lagrange is a time-integrator build in the Moreau-Jean formalism. It uses the central difference method for time-integration, which makes it explicit and symplectic. And the contact is enforced on velocity by a Lagrange multiplier, which ensures an high precision and stability.

The next chapter is devoted to numerically demonstrate the properties and benefits of such a time-integrator.

Chapter 2

A benchmark for non-smooth contact dynamics In this chapter, numerical test cases are presented in order to form a benchmark for time-integrators in non-smooth dynamics. Each case aims to highlight a property required for simulating dynamics problems with impacts. The complexity is elementary: most of the benchmark cases contain one or two DOFs. The implementation is then easily done and does not require a FE software. The contact relates only a rigid motionless boundary, and either a rigid or a deformable body. The contact between two deformable bodies is not addressed. Most of the cases do not integrate friction for keeping the systems conservative. The discrete solutions are computed with the time-integrators presented in section 1.3.3 in order to compare them and highlight their properties.

This section is based on a previous published article [START_REF] Di Stasio | Benchmark cases for robust explicit time integrators in non-smooth transient dynamics[END_REF]. But additional results are provided for other time-integrators, and with several DOFs at contact is changed.

The bouncing ball

This first test case focuses on impact. The problem setting is simple: a point rigid ball only submitted to gravity falls against the ground. It is quite common in literature about non-smooth dynamics [START_REF] Chen | A nonsmooth generalizedα scheme for flexible multibody systems with unilateral constraints[END_REF][START_REF] Fetecau | Nonsmooth Lagrangian Mechanics and Variational Collision Integrators[END_REF][START_REF] Schindler | Timestepping schemes for nonsmooth dynamics based on discontinuous Galerkin methods: Definition and outlook[END_REF]. Indeed the continuous solution is analytically known. And its simplicity facilitates the analysis of the time-integrator in regard of a rigid impact.

The problem is described by fig. 2.1. The contact obeys to the Newton's impact law with a restitution coefficient e c ∈ [0..1]. For e c = 1, the impact is elastic: all kinetic energy is given back. For e c = 0, the impact is plastic: all kinetic energy is absorbed through the impact. An interesting case is when 0 < e c < 1. As energy is dissipate at impact, the high reached by the ball at each bounces decreases. And the time between two impacts becomes

m 1 kg g 9.81 m s -2 ż0 0 m s -1 z 0 1 m h 1 × 10 -2 s Table 2.1: Bouncing ball -Numerical values
shorter. The end of the movement illustrates the "Zeno's paradox": an infinite accumulation of closer and closer impacts in a finite time. This timecontinuous phenomenon is not always trivially transposed in the discrete solution.

The analytical solution (for a null initial velocity) gives the following impact times t i :

∀k > 0          t i,k = 2z 0 g 2 × 1 -e k c 1 -e c -1 if e c = 1 t i,k = 2z 0 g (2k -1) if e c = 1 (2.1)
The stopping time for e c < 1 is:

t final = 2z 0 g 1 + e c 1 -e c (2.2)
The problem is set with the values in table 2.1. The time-step h is set to an high value to spotlight the defaults of the schemes. Indeed no stability issue is present here.

The tested schemes are: the Moreau-Jean scheme defined in eqs. (1.122) to (1.124); the Paoli-Schatzman scheme defined in eqs. (1.138) to (1.140); and the CD-Lagrange scheme defined in eqs. (1.157) to (1.159). They all integrate the Newton's impact law, necessary to model the rigid impact. The fig. 2.2 represents the position of the ball over time for e c = 1. For the Paoli-Schatzman scheme, the ball reaches a position higher than the initial one. This illustrates the remark made in section 1.3.3 on the non-conservative contact work: energy can be injected at impact. If the contact work tends to be conservative with time refinement, this should cause stability issue especially in non-linear problems. On the contrary, both CD-Lagrange and Moreau-Jean schemes raise at each bounces the initial height. Indeed with e c = 1 the discrete energy is conserved at impact (see section 1.3.3). For the three schemes, the ball penetrates the boundary at impact. This results from the discrete contact conditions. For the Moreau-Jean and CD-Lagrange schemes, the discrete contact conditions approximate the persistency condition which does not ensure a null penetration. And for the Paoli-Schatzman scheme, if the discrete contact conditions are based on the impenetrability ones, the approximation of the Newton's impact law lets penetration happens.

On fig. 2.3 the restitution coefficient is set to e c = 0.75. Energy is then dissipates at impact, and the end of the movement illustrates the Zeno's paradox. The three schemes simulates correctly the stopping movement.

The fig. 2.4 depicts the position of the ball for e c = 0. In this case, the contact conditions of Paoli-Schatzman (see eq. (1.159)) match exactly the impenetrability conditions. The discrete solution does indeed not penetrate at impact. Otherwise, both CD-Lagrange and Moreau-Jean schemes present a residual penetration after the impact. In order to evaluate the convergence rate, the following Hausdorff measure is used [see 1]:

e = 1 2 h i∈[1,n] |f i -f (t i )| (2.3)
where:

• n is the number of discrete times in interval [t 0 , t f ];

• f i is the discrete solution at time t i ;

• f (t i ) is the analytical or reference solution at time t i .

The norm (2.3) evaluates the distance between the point graph from the discrete solution and the continuous graph from the analytic one. It provides then an error indicator. The fig. 2.5 shows the error between the discrete solutions and the analytical one. The case is set with the values of table 2.1 and e c = 1. For the three schemes the convergence rate is in O(h), which is common for non-smooth dynamics scheme with impacts [see [START_REF] Acary | Higher order event capturing time-stepping schemes for nonsmooth multibody systems with unilateral constraints and impacts[END_REF][START_REF] Chen | A nonsmooth generalizedα scheme for flexible multibody systems with unilateral constraints[END_REF][START_REF] Schindler | Timestepping schemes for nonsmooth dynamics based on discontinuous Galerkin methods: Definition and outlook[END_REF]. The Zeno's paradox and event-driven schemes The non-smooth timeintegrators are sometimes classified between event-driven schemes and timesteppping schemes. An event-driven time-integrator detects the non-smooth events, and integrates the smooth dynamics between them. Some examples can be found in [START_REF] Acary | Numerical Methods for Nonsmooth Dynamical Systems 1st ed[END_REF]. If they seem interesting as the time-integration concerns only smooth dynamics, they are never used in practice. Indeed the time-step is necessarily equal at least at the time between two impacts. In the Zeno's paradox for example, an event-driven scheme will never reach the final time. If this situation is quite specific and extreme, a similar one is met in FE simulations. On a meshed body, non-smooth events correspond to each impact of nodes. If a group of nodes comes in contact at close but different discrete times, the time-integrator generates a succession of small time-steps. This causes an increase in computation time, which can makes the computation too costly. On the contrary, the time-stepping schemes keep a constant time-step whatever the number of non-smooth events. They can simulate the Zeno's paradox; and then close nodal impacts in FE simulations do not increase the computation time. All time-integrators presented here are timestepping, being build in the NSCD framework. The equilibrium equation are integrated in a weak sense in time. The precise times of non-smooth events are not necessary to solve the system.

Conclusion

The bouncing ball underlines the importance of using the nonsmooth contact dynamics framework to design the time-integrators. Firstly it brings stability. Indeed the standard time-integrators diverge at impact on this test case. For example, in [START_REF] Chen | A nonsmooth generalizedα scheme for flexible multibody systems with unilateral constraints[END_REF] discrete solutions with Newmark schemes and Lagrange multipliers are presented. They show a brutal divergence at the first impact. On the contrary, the NSCD time-integrators as the Moreau-Jean, the Paoli-Schatzman or the CD-Lagrange scheme are stable as demonstrated by the preceding results. And secondly, the NSCD framework leads to time-stepping schemes.

The schemes based on velocity contact conditions, the Moreau-Jean and CD-Lagrange schemes, have an advantage on Paoli-Schatzman scheme. They are conservative at impact. On non-linear problems, this is crucial for stability.

The Van der Pol oscillator

The Van der Pol oscillator does not model a mechanical situation. It consists in an oscillator with a non-linear damping term, and a "contact boundary" to introduce non-smooth events. The system equation is:

d ẋ -ξω 0 1 - x 2 x 2 0 ẋdt + ω 2 0 xdt + dr = 0 (2.4)
The "contact impulse" r enforces a position x greater than a lower bound x R which represents a rigid boundary. The contact involves a Newton's impact law with e c = 1 in order to be conservative. The behaviour of the system depends on the two parameters ξ and ω 0 . ξ sets the preponderance

ξ 5 ω 0 1 s -1 ẋ0 1 m s -1 x 0 1 m h 1 × 10 -3 s [t 0 , t f ] [0, 300] s x R -1.2 m
Table 2.2: Van der Pol oscillator -Numerical values of the non-linear terms: for ξ ω 0 , the system tends to the standard linear oscillator.

This system is particularly interesting for long-time simulations. Indeed it tends to a limit cycle in phase space, reached by the discrete solution only if the time-integrator is conservative. Its non-linearity on displacement requires generally a non-linear solver even for this one DOF system. It highlights then two key features of the central difference method. Firstly no non-linear solver is required, the scheme being explicit for the displacement update. And secondly the estimation of internal forces (here the non-linear damping term) thanks to the velocity at t n+ 1 2 makes the dynamics (2.4) fully explicit.

The numerical values for setting the test case are gathered in table 2.2. They set a preponderant non-linear term with high velocities at impact, and a fast convergence to the limit cycle for the solution. The simulation covers a long time interval in order to verify the energy conservation. On this test case, only the CD-Lagrange and the Moreau-Jean schemes are tested. The Moreau-Jean dynamics is solved thanks to a Newton algorithm to handle the non-linearity in displacement. No such algorithm is required by the CD-Lagrange.

The fig. 2.6 represents the discrete positions x over the first 30 s for the Moreau-Jean and the CD-Lagrange schemes. The black line is the contact boundary, where elastic impacts happen. The velocity at impact is at its maximum, which causes large velocity jumps. The non-linearities appear in the discrete position through the sharp changes in gradient. Nevertheless the both schemes are robust. The CD-Lagrange is close from the Moreau-Jean scheme for a lower algorithmic complexity and an extra approximation for evaluating the damping term.

The discrete solution on the full simulation time is depicted in fig. 2.7, with the position on X-axis and velocity on Y-axis. The velocity jump is clearly visible, as the limit cycle. Both scheme stays on the limit cycle thanks to their energy conservative formulations. Remark 8. In fact the energy conservation for the Moreau-Jean scheme is not exact. It depends on the precision of the Newton-Rahpson's algorithm used for solving the dynamics. For the CD-Lagrange, the dynamic is exactly solved as no solver is involved. Here this approximate solving is overcome by setting the Newton-Raphson to the software precision.

For studying convergence, a discrete solution is computed with a short time-step of h = 1 × 10 -6 s. The time-integration is done with the CD-Lagrange, but with the damping term evaluated with the velocity at t n+1 . The formulation becomes implicit in velocity, but the solving is resumed here at a scalar division. This discrete solution is considered as a reference. The error is evaluated thanks to the norm eq. (2.3) as the distance between the discrete solution and the reference one.

A first analysis is done with the settings of table 2.1, but with t f = 100s. The fig. 2.8 shows a convergence rate in O(h) for both schemes as in the bouncing ball case. The error is higher for the CD-Lagrange, maybe because of the extra approximation in the damping term.

A second analysis does no longer involve impacts by removing the rigid boundary. This time three discrete solutions are computed: one with the Moreau-Jean scheme, one with the CD-Lagrange and the damping term velocity at t n+ 1 2 , and one with the CD-Lagrange and the damping term velocity at t n+1 . The convergence rates on fig. 2.9 are in O(h 2 ) for the Moreau-Jean scheme and the CD-Lagrange with the exact damping term. But the CD-Lagrange with the damping term at t n+ 1 2 converges only in O(h). This extra approximation has then a cost on precision. A convergence rate of O(h 2 ) is indeed expected for the Moreau-Jean and the CD-Lagrange

Conclusion

The Van der Pol oscillator highlights the advantages of the CD-Lagrange time-integrator. The explicit displacement and the velocity in internal forces at t n+ 1 2 make the time-integrator explicit even for nonlinear internal forces. On the contrary, the Moreau-Jean scheme requires a non-linear solver even on this one DOF case. For a non-smooth system, they have the same convergence rate. But for a smooth system (no impact), making explicit the damping term decreases the convergence rate of CD-Lagrange scheme. An other highlight is the importance of conservative timeintegrators for a long-time simulation.

The rotating spring

The rotating spring focuses on conservation of angular momentum in large rotations. It consists simply in a mass-spring system in rotation with a circular contact boundary. For this test-case, the time-integrator needs to conserve the energy and the angular momentum under non-smooth events. This feature characterizes a symplectic time-integrator.

The fig. 2.10 describes the system. As it evolves in two dimensions, the contact could be frictional. But here the focus is done only on unilateral contact, as the friction is addressed by the "impacting dome" test case of

m 1 kg k 10 kg s -2 l 0 1 m ẋ0 1e x + 2e y m s -1 x 0 0.8e x m h 1 × 10 -1 s e c 1 x R 1.4 m
Table 2.3: Rotating spring -Numerical values section 2.5. The restitution coefficient is set to e c = 1 for conserving energy through impacts. In [START_REF] Di Stasio | Benchmark cases for robust explicit time integrators in non-smooth transient dynamics[END_REF], discrete solutions with friction are presented for both CD-Lagrange and Moreau-Jean scheme.

The system equations are:

md ẋ + k 1 - l 0 x x = dr n (2.5)
r is the contact impulse, only along the normal n as no friction is considered. The position vector x has two DOF x = (x, y). The eq. (2.5) shows a nonlinearity on x . The Moreau-Jean scheme requires then again a non-linear solver, whereas the CD-Lagrange schemes stays explicit. The table 2.3 gathers the numerical settings. The time-step is the same for the two schemes, and close from the stability limit of the CD-Lagrange. The fig. 2.11 shows the position x over the first 10 s of the simulation. The motion is composed of a general rotation counter clockwise, and of elastic bounces against the rigid boundary. The penetration due to the contact conditions on velocity is visible for the both schemes. On the first 10 s, the discrete solutions stay close. But a closer look on the angular momentum on fig. 2.12 shows a crucial difference. The plotted angular momentum is L = x × v (with × the crossproduct), which is the continuous one. The fig. 2.12 represents it over 100 s. For the CD-Lagrange, it is constant over the simulation staying exactly at its initial value. On the contrary, the angular momentum for Moreau-Jean scheme decreases during the simulation. The θ-method with θ = 1/2 is yet symplectic in the sense given in [START_REF] Kane | Symplectic-energy-momentum preserving variational integrators[END_REF][START_REF] West | Variational integrators and the Newmark algorithm for conservative and dissipative mechanical systems[END_REF]: the preserved angular momentum is not the physical one. Indeed the angular momentum oscillates in smooth phases, around a constant value. But each non-smooth events change this average value, which leads here to a drop of angular momentum. Note that this behaviour is deeply influenced by the precision used for the convergence criteria set in the Newton-Raphson algorithm.

This decrease has no influence on the energy balance. The fig. 2.13 shows the energy balances between [t n+ 1 2 , t n+ 3 2 ] for the CD-Lagrange and [t n , t n+1 ] for the Moreau-Jean scheme. Both system energies are constant during the simulation, as the impacts do not work. The remark 8 is still valid: the energy conservation of Moreau-Jean scheme is conditioned to the precision of the non-linear solver.

Conclusion

The rotating spring demonstrates that the conservation of energy is not enough for large rotations. The exact conservation of angular

The impacting bar

This test case involves a one-dimensional deformable solid, discrete in space thanks to a FE method. It consists in a bar, which impacts a rigid boundary at one of its ends. In literature about deformable dynamics with impact, the impacting bar problem is quite used [for example in 10, 14, 20, 25, 32]. Doyen et al. in [START_REF] Doyen | Time-Integration Schemes for the Finite Element Dynamic Signorini Problem[END_REF] provide a complete review of time-integrators for impact on this case. Indeed the continuous problem has an analytic solution [see 28], which can be compared to the discrete solutions. The impacting bar is described by fig. 2.14. This time the contact requires no restitution coefficient as it involves a deformable solid. The bar is discrete in space by a FE method using P1 elements. All of them have the same length. The material is linear and elastic. The numerical values characterising the case are gathered in table 2.4. The time-step is chosen at 90% of the critical time-step of CD-Lagrange scheme. As the Newton's law is no more required, the Laursen-Chawla scheme of section 1.3.3 can be introduced in the benchmark. The Paoli-Schatzman scheme presents here an interesting property. As e c = 0, its contact law matches the impenetrability conditions. It is also added in the benchmark. The Paoli-Schatzman is in these conditions equivalent to the Carpenter's scheme [START_REF] Carpenter | Lagrange constraints for transient finite element surface contact[END_REF].

The fig. 2.15 gathers the discrete quantities at contact node for the Moreau-Jean and Laursen-Chawla schemes. As mentioned in section 1.3.3, these time-integrators differ only in the contact condition. They use the same form on velocity but at a different discrete time: at t n+1 for the Moreau-Jean scheme, and at t n+ 1 2 for Laursen-Chawla. As demonstrated in [START_REF] Laursen | Design of energy conserving algorithms for frictionless dynamic contact problems[END_REF], this ensures the energy conservation at impact. The fig. 2.15d confirms this result. At impact, the system energy decreases for the Moreau-Jean's solution; but stays constant for Laursen-Chawla. Nevertheless it causes symmetric oscillations on the velocity and contact impulse during contact, visible on figs. 2.15b and 2.15c. Here the discrete velocities for laursen-Chawla are represented at t n+1 , but at t n+ 1 2 they do not present oscillations. Despite these oscillations, the discrete positions on fig. 2.15a are stable during contact. 2.16d shows the sum of the system and complementary energies (all conservative terms). An energy loss is visible at impact. Here it stays small at ∼1% of the initial energy, but for a system with a large number of contact nodes or with multiple impacts this could lead to an important energy loss.

For estimating the error, the norm (2.3) is used to measure the distance between the graphs of the discrete solution and the analytic one at contact node. The figure fig. 2.17 shows the convergence rates, both on position and contact impulse. On position, the four schemes converge in O(h). But on force, the oscillations makes the Laursen-Chawla scheme non-converging.

Conclusion

The impacting bar highlights a major drawback: energy is dissipated at impact. If the energy conservation at impact is achieved by the Laursen-Chawla scheme, this leads to symmetric oscillations on velocity during contact. These oscillations do not interfere with the displacement, which stays constant during contact. But they cause oscillations on impulse, which inhibit the convergence toward the analytical solution.

For the other schemes, the CD-Lagrange and Moreau-Jean have similar properties: the contact is stable but presents a residual penetration. On the contrary, the Paoli

The impacting dome

The impacting dome is a multi-DOF problem in three dimensions. The space discretization is done by a FE method. A "dome", a rectangular cuboid topped by an half cylinder, impacts a rigid boundary. The contact concerns several nodes, and involves friction described by the Coulomb's law (1.57). It highlights the numerical performance of the CD-Lagrange scheme in such a case. 2.5. At t 0 the dome and the cube are distant from a distance g 0 along y . The following results are computed with the mesh 3 of table 2.6. As the contact involves a rigid boundary analytically described, the Delassus operators of contact LCP are diagonal. On this test case, only the CD-Lagrange is tested. It is implemented in MEF++. This software is developed jointly by the GIREF (Laval University) and industrial partners as Michelin. It implements the finite elements method for several types of problems in deformable mechanics. MEF++ aims to be a research software while being enough efficient to address industrial simulations. Indeed, based on PETSc library, it offers a large choice of fast parallel solvers. With a large collection of material laws and physical equations, it can simulate the complex problems met in tire simulations. The fig. 2.19 shows the deformation of the dome at four times: impact, during contact, at maximum normal impulse, and release. The dome is colored according to the magnitude of normal impulse. The arrows represent the tangential impulse: scaled proportionally to the magnitude and in the same direction. Both on fig. 2.19a at impact and on fig. 2.19d at release, the penetration due to contact conditions in velocity is visible. The global deformation shows no instability. The tangential part of contact impulse is as expected along x , the tangential direction of initial velocity.

On fig. 2.20, an average is done for displacements and nodal impulses on the nodes at the top of the dome. The curves represent the normal impulse along -y , and the tangential one along -x (x being the main tangential direction). The displacement is projected along y , the opposite to contact normal. The right vertical axis gives the displacement scale, and the left one the impulse scale. The complementary condition is clearly visible: positive impulses correspond to null velocities. Indeed the top nodes have the same discrete impact time. This is visible on fig. 2.19a. The normal and tangential impulses are stable during contact, with smooth variations. The tangential impulse is inferior to the sliding limit: the normal impulse multiplied by the friction coefficient. No sliding happens.

The energy balance on fig. 2.21 is quite accurate. The energy loss associated to normal contact work is small in front of the tangential work due to friction. Here, an energy loss happens for each impact on a node. With a closer look on normal contact work, two discrete impact times appear. The In order to numerically check the space-time convergence, the dome of fig. 2.18 is meshed with different element sizes. They are gathered in table 2.6, with the critical stability time-steps (exactly computed thanks to eq. (1.98)) and the time-steps used in simulations. The reference solution is computed on the finest mesh (mesh 5). The error is computed as the distance between the graph of a discrete solution and the reference one with the norm (2.3). The looked quantity is the average displacement for the top nodes of the dome. The fig. 2

Conclusion

The impacting dome is a more complex test case as it requires a finite element code. It demonstrates the ability of the CD-Lagrange in case of contact between a deformable and a rigid motionless body: stable contact quantities, and a small energy loss at impact. From an algorithmic point of view, adding friction does not bring any extra complexity. Both normal and tangential LCP are diagonal, which makes their solving local to the nodes. With the lumped mass matrix, the CD-Lagrange algorithm is then naturally parallel. The explicit and parallel features ensure an high computational efficiency.

Conclusion: the CD-Lagrange as a promising scheme

The preceding benchmark demonstrates that the CD-Lagrange scheme is well suited to nonsmooth dynamics. It enforces the contact on velocity with the discrete contact conditions of Moreau-Jean. Thanks to this formulation, the scheme can address rigid-rigid and deformable-rigid contact with an high accuracy. For a rigid-rigid deformable contact, the discrete solution conserves energy through elastic impacts as shown by the bouncing ball in section 2.1. For a deformable-rigid contact as in impacting bar of section 2.4, the discrete persistency contact condition brings an high stability for contact 94 quantities despite a residual penetration. If the accuracy is ensured for the contact quantities, the time-integrator completes it at a global level. The central difference method is indeed a symplectic time-integrator which ensure an accurate long-time integration. It conserves a discrete form of energy and the exact value of angular momentum as shown by the Van Der Pol oscillator in section 2.2, or the rotating spring in section 2.3. Besides its accuracy, the CD-Lagrange is explicit even for the non-linear damped system of the Van Der Pol oscillator of section 2.2. The solving of the discrete dynamics does not require a solver which leads to an high efficiency. If the contact involves an analytically described body or conforming contact meshes, the contact LCPs are explicit too. All operations of the CD-Lagrange are local, and the algorithm is naturally parallel. The impacting dome in section 2.5 shows this efficiency with a case close from a real one: multiple DOFs are in contact with friction, on a deformable meshed body.

But the CD-Lagrange is not yet operable for industrial simulations where the deformable-deformable contact is a key point. The deformable-deformable contact requires to transfer the contact quantities between two meshes. Generally they are not conforming and the quantities must be projected between them. In deformable mechanics, the mortar methods [START_REF] Belgacem | The Mortar finite element method with Lagrange multipliers[END_REF][START_REF] Wohlmuth | A Mortar Finite Element Method Using Dual Spaces for the Lagrange Multiplier[END_REF] handle this with accuracy. But integrating the mortar methods in CD-Lagrange scheme must be made by keeping the accuracy and efficiency.

The other drawback of the CD-Lagrange scheme is the energy loss at impact for a contact involving a deformable body. This energy loss is common to all the schemes based on the velocity contact conditions of Moreau-Jean. The impacting bar of section 2.4 demonstrates that the energy conservation at impact is linked to the persistency condition. If the persistency condition is ensured, the scheme is stable but energy is lost at impact. If the energy conservation is ensured, the velocity is highly oscillating at contact. The remark made in section 1.3.3 can be extended: this is the loss of kinetic energy which ensures the contact stability. An interesting way for energy conservation at impact is the singular mass method [START_REF] Khenous | Mass redistribution method for finite element contact problems in elastodynamics[END_REF] which leads to massless contact nodes.

The two following sections are devoted to explore these two improvements: adapting the singular mass method to the CD-Lagrange scheme in order to ensure energy conservation at impact; and addressing the deformabledeformable contact thanks to the mortar methods.

Chapter 3

The singular mass method and the CD-Lagrange scheme This section is based on the pre-print article: Di Stasio J., Dureisseix D., Georges G., Gravouil A., Homolle T., An explicit time-integrator with singular mass for non-smooth dynamics, , submitted in 'Computational Mechanics' on 11 January 2021. The first three sections follows the pre-print, but the analyses are developed. The fourth section contains new results.

This chapter is devoted to adapt the singular mass method for the CD-Lagrange scheme. The goal is to improve the energy balance by making the impact conservative. After a review of existing schemes with a singular mass, a first singular CD-Lagrange is presented on a 1D formulation. It highlights the difficulties to adapt the singular mass in an explicit scheme, and gives solutions to overcome them. The 1D formulation is then extended to 3D in two different ways.

The singular mass method in existing schemes

The contact constraints change deeply the structure of the dynamical problem. For elastodynamics, the semi-discrete in space problem is an ordinary differential equation (ODE) with a regular solution twice-derivable in time. The time-integration of such a smooth system is then quite facilitated. By adding unilateral contact constraints, the equation becomes a differential algebraic equation (DAE). As a consequence the time-integrators developed for ODEs can not be straightforwardly extended to DAEs, because of stability issues [START_REF] Cardona | Time integration of the equations of motion in mechanism analysis[END_REF][START_REF] Farhat | Implicit time integration of a class of constrained hybrid formulations-Part I: Spectral stability theory[END_REF]. As mentioned by Acary in [START_REF] Acary | Projected event-capturing time-stepping schemes for nonsmooth mechanical systems with unilateral contact and Coulomb's friction[END_REF], reducing the index of the DAEs enhances the stability. For example, the NSCD [START_REF] Jean | Unilaterality and dry friction in the dynamics of rigid body collections in 1st Contact[END_REF][START_REF] Moreau | Nonsmooth Mechanics and Applications[END_REF] performs a reduction of the DAEs index by enforcing the constraint at velocity level. The singular mass method [START_REF] Khenous | Mass redistribution method for finite element contact problems in elastodynamics[END_REF] achieves also an index reduction, and is explored in the following on the CD-Lagrange.

96

Indeed for a discrete-in-space dynamics system under contact constraints, the solution is no longer unique and presents non-smooth events. For example in [START_REF] Khenous | Mass redistribution method for finite element contact problems in elastodynamics[END_REF], an infinite number of solutions is exhibited for a mass-spring system with an unilateral contact. If quite simple, this system is nevertheless critical for structural dynamics because it represents the normal problem for a contacting node. Nevertheless, some uniqueness results can be recovered on the space-and-time discrete system. The work of Paoli & Schatzman in [START_REF] Paoli | A Numerical Scheme for Impact Problems I: The One-dimensional Case[END_REF][START_REF] Paoli | A Numerical Scheme for Impact Problems II: The Multidimensional Case[END_REF] introduces a scheme in the NSCD framework. Their scheme is proven to converge to an unique solution thanks to its restitution coefficient.

An other way to recover the uniqueness for elastodynamics is the singular mass method. Introduced in [START_REF] Khenous | Mass redistribution method for finite element contact problems in elastodynamics[END_REF], it modifies directly the semi-discrete in space system: the DOFs constrained by contact get a zero entries in mass matrix. This makes elliptic the equations for these contact DOFs. Indeed the inertia term associated to acceleration disappears with the cancellation of the mass entry. The proofs of uniqueness can be found in [START_REF] Dabaghi | A weighted finite element mass redistribution method for dynamic contact problems[END_REF][START_REF] Khenous | Mass redistribution method for finite element contact problems in elastodynamics[END_REF][START_REF] Renard | The singular dynamic method for constrained second order hyperbolic equations: Application to dynamic contact problems[END_REF]. Besides the uniqueness result, the persistency condition is automatically fullfilled if the impenetrability one is enforced. This leads to an energy conservative formulation for the semi-discrete system in space [see [START_REF] Khenous | Mass redistribution method for finite element contact problems in elastodynamics[END_REF][START_REF] Renard | The singular dynamic method for constrained second order hyperbolic equations: Application to dynamic contact problems[END_REF].

These new semi-discrete properties bring enough stability to use the timeintegrators from HHT-Newmark family. If they are unstable with a standard mass matrix, the stability is recovered together with an accurate energy balance with the singular mass. Such numerical results are observed in [START_REF] Doyen | Time-Integration Schemes for the Finite Element Dynamic Signorini Problem[END_REF]. For the HHT-Newmark schemes, the energy blows up and the contact stresses highly oscillate with a standard mass matrix. But with a singular mass matrix, the energy balance is only slightly dissipative and the contact stresses become smooth. [START_REF] Dabaghi | A robust finite element redistribution approach for elastodynamic contact problems[END_REF] obtains other interesting numerical results, but this time on the convergence rate. The Crank-Nicholson scheme diverges with a consistent mass matrix but converges with a singular one. For the other schemes, the singular formulation improves at less the convergence rate. Finally in [START_REF] Dabaghi | A weighted finite element mass redistribution method for dynamic contact problems[END_REF][START_REF] Dabaghi | A robust finite element redistribution approach for elastodynamic contact problems[END_REF][START_REF] Khenous | Mass redistribution method for finite element contact problems in elastodynamics[END_REF], a stabilization of contact stresses is observed for Crank-Nicholson, Newmark (β = γ = 1/2) and backward Euler schemes. The energy balance is also improved: the energy blow-up disappears for Crank-Nicholson method, and the energy dissipation is smaller for the Newmark's (β = γ = 1/2) scheme.

This better stability for contact stresses is explained in [START_REF] Krause | Presentation and comparison of selected algorithms for dynamic contact based on the Newmark scheme[END_REF]. A crucial difference exists between the continuous and the space and time discrete model with a consistent mass matrix. On the discrete model, the contact stresses equilibrate both the inertia of contacting nodes and the internal stresses. On the contrary in the continuous model, the contact stresses equilibrate only the internal ones. The continuous contact boundary is indeed a surface without a mass. The singular mass method allows to retrieve a massless contact boundary, but on the discrete model.

Using a singular mass is not only attractive for schemes of Newmark's family. In [START_REF] Hager | A stable energy-conserving approach for frictional contact problems based on quadrature formulas[END_REF], Hager et al. apply the singular mass method in the Laursen-Chawla scheme [START_REF] Laursen | Design of energy conserving algorithms for frictionless dynamic contact problems[END_REF]. Already energy conservative thanks to a discrete persistency contact condition, its main drawback is the numerical bounded oscillations for contact stresses (see section 2.4). These symmetric oscillations keep constant the kinetic energy of contact nodes, where the CD-Lagrange or Moreau-Jean schemes cancel it to stabilize the contact. The assumption in [START_REF] Hager | A stable energy-conserving approach for frictional contact problems based on quadrature formulas[END_REF] is that the oscillations of the Laursen-Chawla scheme should be reduced for massless contact nodes. And that is confirmed by the numerical results.

If the final result is similar: a zero mass entry for contact DOFs, several methods lead to a singular mass matrix. In [START_REF] Khenous | Mass redistribution method for finite element contact problems in elastodynamics[END_REF], the singular mass is obtained through a minimization problem. The first constraint is a null mass for the DOFs along the normal of contact boundary. The second constraint is to keep constant the total mass, the inertial momentum and the center of gravity of the system between the consistent and the singular mass matrices. With this method, the singular mass matrix depends on the deformed configuration as depending of the contact normal. This could represent an extra-computational cost especially with large displacements. In [START_REF] Hager | A stable energy-conserving approach for frictional contact problems based on quadrature formulas[END_REF] the singular mass matrix comes from modified quadrature formulas in the initial assembly. The computational cost is lower as the matrix is computed only once, and directly obtained by an integration. The previous mass quantities of the global system are still preserved. With this method, all DOFs of contact nodes have a null mass entry, not only those along the contact normal. In [START_REF] Renard | The singular dynamic method for constrained second order hyperbolic equations: Application to dynamic contact problems[END_REF], the singular mass matrix results from two different Galerkin discretization spaces for velocity and displacement. Hauret in [START_REF] Hauret | Mixed interpretation and extensions of the equivalent mass matrix approach for elastodynamics with contact[END_REF] and Tkachuk et al. in [START_REF] Tkachuk | Hybrid-mixed discretization of elasto-dynamic contact problems using consistent singular mass matrices[END_REF] continue this method in a variational framework. The displacement, the velocity and the linear momentum are three independent variables, each one with its own discrete functional space. And lastly in [START_REF] Dabaghi | A weighted finite element mass redistribution method for dynamic contact problems[END_REF], weighting coefficients are introduced for the shape functions used to compute the mass matrix. They allow to get several shapes and patterns for the matrix which influence the discrete solution.

The shape and the computation method of the singular mass matrix are not addressed here. Indeed the improvements in stability and energy balance are observed whatever the shape of the matrix. This section is devoted to adapt the singular mass matrix to the CD-Lagrange scheme, and analyses the consequences on the discrete solutions and the energy balance. The first step is the design of a singular mass formulation in 1D on the impacting bar. And then, two extensions are proposed for 3D meshes.

The Impacting bar with a singular mass

In this section, the impacting bar case of section 2.4 is used to design a 1D formulation of the CD-Lagrange with a singular mass. If the case is simple, it highlights the main difficulties and improvements raised by the singular formulation. As mentioned before, this test case is quite used in literature about non-smooth dynamics. Specific results with a singular mass matrix are presented in [START_REF] Dabaghi | A weighted finite element mass redistribution method for dynamic contact problems[END_REF][START_REF] Dabaghi | A robust finite element redistribution approach for elastodynamic contact problems[END_REF]. And the general review of Doyen et al. in [START_REF] Doyen | Time-Integration Schemes for the Finite Element Dynamic Signorini Problem[END_REF] gives also results with a singular mass. The impacting bar is described by fig. 3.1. The problem setting is the same as in section 2.4: the material is elastic and linear and the finite elements are P1. A node corresponds then to one DOF. The initial spatial discretization is modified to get a zero mass entry for the contact node. The mesh is split in two parts: the bulk, which gathers the elements with mass nodes and the skin constituted by the last element with the massless contact node and the connected node on the bulk. The contact condition is applied on γ C , the contact boundary. The interface between the bulk and the skin is γ B .

The singular mass matrix comes here simply from the cancellation of the entry associated to the contact node. Once the mass matrix is assembled and lumped, the contact entry is set to zero. This does not conserve the total mass, but no comparisons will be made between the consistent and singular solutions.

The bulk is described as a standard system: all nodes have a mass. The following notations denote the quantities for the bulk:

• U, U, Ü: displacement, velocity and acceleration vectors gathering DOFs of discrete fields for the bulk;

• X 0 : initial positions of nodes;

• M: mass matrix (consistent as involving only bulk nodes);

• F int (U): internal forces in the bulk;

• F ext : external forces;

• ū: displacement for the DOF on γ B ;

• L: selects the DOF on γ B , ū = LU (in a sense, L is a projection along n on γ B ).

The assemblies of mass matrix and internal forces differ. The global mass matrix M * is assembled on all elements of bulk and skin. Then M * is made singular by setting the entry of the contact node to zero. M denotes the sub-matrix of M * for the entries corresponding to the nodes of bulk. M is then consistent. On the contrary, the internal forces of bulk come from an assembly only for the bulk elements. The force applied by the skin on the bulk is denoted f . The semi-discrete in space system for the bulk is:

M Ü = F ext -F int (U) + Lt f (3.1)
This system is not submitted to any unilateral contact constraint: Ü is well defined as f is the elastodynamics stress in the skin.

The skin system is particular as it has a massless node. It gathers two DOFs: one for the node on γ B and one for the node on γ C . The skin is described by:

• ū: displacement of the node on γ B , ū = LU;

• u: displacement of the contact massless node on γ C , u = LU;

• L: select the massless node of γ C , u = LU (in a sense, L is a projection along n on γ C );

• x 0 : initial position for contact node;

• δ u = uū: difference of displacements for skin nodes, positive if the skin is compressed, negative if stretched;

• l = (x 0 + ū) -( x 0 + u) = l 0 -δ u: thickness of skin, l 0 is the initial thickness, with x0 = LX 0 ;

• k: rigidity of the skin.

The choice of initial skin length and skin stiffness is here straightforward, as taken from the consistent problem. The initial skin length is the element length, and the skin rigidity comes from the rigidity matrix. The skin is described as a static system:

0 0 = -f f - k -k -k k ū u (3.2)
The eq. ( 3.2) added to eq. ( 3.1) form a system on all nodes. But the acceleration and velocity of the contact node are not determined as no inertial term is associated. The eq. ( 3.2) condenses to give the expression of contact force, which is also the force of the skin on the bulk:

f = k δ u (3.3)
The CD-Lagrange, described by eqs. (1.157) and (1.158), performs the time-discretization. The space-and-time discrete system for the bulk is:

U n+1 = U n + h Un+ 1 2 (3.4) M( Un+ 3 2 -Un+ 1 ) = h F ext n+1 -F int (U n+1 ) + Lt r n+ 3 2 (3.5)
In eq. (3.4), the nodal velocities are crucial to determine the updated displacement. They are determined by eq. (3.5) thanks to the inertial term.

No such a relation exists for the skin. The skin equation 3.3 gives only the contact impulse, the skin force integrated over a time-step:

r n+ 3 2 = h f n+1 = h k δ u n+1 (3.6)
For a consistent mass, the CD-Lagrange is completed by the contact conditions on velocity of eq. (1.159). They act on the velocity of γ C , ˙ u n+ 3 2 , which does not appear in preceding equations.

A time-integration relation is then added for the massless contact node:

u n+1 = u n + h ˙ u n+ 1 2 (3.7)
And a contact law is proposed to determine ˙ u n+ 3 2 and enforce the contact conditions. The enforcement is closer as possible of the contact conditions of the CD-Lagrange with consistent mass matrix. But some adjustments are required as the contact node has no inertia.

In order to describe the contact law, the following notations are introduced:

• g n+1 ,the gap at the contact node, equal here to g n+1 = x 0 + u n+1 ;

• Ufree

n+ 3 2
, the free velocity of the bulk, i.e. the bulk velocity without skin action.

The expression of free-velocity comes from the eq. (3.5) with a null skin force:

Ufree n+ 3 2 = Un+ 1 2 + hM -1 F(U n+1 ) (3.8) ufree n+ 3 2 = L Ufree n+ 3 2 (3.9)
Finally the proposed contact law is: are equal only if the skin is at rest (δ u n+1 = 0 and then r n+ 3 2 = 0). And for active contact state, three cases emerge. They are described in table 3.1.

If g n+1 > 0, ˙ u n+ 3 2 = ufree n+ 3 2 (3.10) Else g n+1 0,    If r n+ 3 2 > 0, ˙ u n+ 3 2 = 0 Else r n+ 3 2 0, ˙ u n+ 3 2 = ( ufree n+ 3

Impact g n+1 0

The skin is at initial length, the contact node is stopped.

r n+ 3 2 = 0 ˙ u n+ 3 2 = 0 Contact g n+1 0
The skin is compressed, the contact node does not move.

r n+ 3 2 > 0 ˙ u n+ 3 2 = 0 Release g n+1 0
The skin is at rest or stretched, the contact node moves only to release contact.

r n+ 3 2 0 ˙ u n+ 3 2 0 Table 3.1: Contact law -Quantities on γ C in contact phases
In table 3.1, the skin is considered in its initial rest state at impact. This assumption is verified in the impacting bar problem considered here. But for other initial conditions, the skin tends to relax into rest state during a free-of-contact phase.

The system is finally governed by the following set of equations:

U n+1 = U n + h Un+ 1 2 (3.12) M( Un+ 3 2 -Un+ 1 2 ) = hF(U n+1 ) + Lt r n+ 3 2 (3.13) u n+1 = u n + h ˙ u n+ 1 2 (3.14)
r n+ 3 2 = h k δ u(ū n+1 , u n+1 )
+ Velocity contact law: eq. (3.10) and eq. (3.11)

(3.15)
The unknowns of eqs. (3.12) to (3.15) are the displacement U n+1 and u n+1 , and the velocities Un+ . The displacement are explicitly computed from the velocities thanks to eqs. (3.12) and (3.14). The velocities of bulk come from the dynamics (3.13), and the velocity of contact node from the contact law which depends on the skin state through eqs. (3.15).

The algorithm for one time-step is given in algorithm 5. Its steps are similar to those of the algorithm 3 without singular mass.

Algorithm 5 CD-Lagrange with singular contact mass

1:    U n+1 = U n + h Un+ 1 2 u n+1 = u n + h ˙ u n+ 1 2

Update of position

Velocity correction

The singular mass formulation has interesting energetic properties. The energy balance is done between t n+ 1 2 and t n+ 3 2 . As mentioned in section 1.3.1, this cancels the complementary term which simplifies the analysis of contact work. The first step is to establish the energy balance for the bulk:

∆E k = ∆ Wint + ∆W skin (3.16)
With:

∆E k = 1 2 Ut M U n+ 3 2 n+ 1 2
(kinetic energy) (3.17)

∆ Wint = 1 2 h Un+ 3 2 + Un+ 1 2 t F(U n+1 ) (work of internal forces) (3.18) ∆W skin = 1 2 un+ 3 2 + un+ 1 2 r n+ 3 2 (work of skin force) (3.19)
∆W skin is the work applied by the skin on the bulk, and it can be decomposed into:

∆W skin = ∆ W int + ∆W c (3.20)
With:

∆ W int = - 1 2 δ ˙ u n+ 3 2 + δ ˙ u n+ 1 2 r n+ 3 2 (3.21) ∆W c = 1 2 ˙ u n+ 3 2 + ˙ u n+ 1 2 r n+ 3 2 (3.22)
Here ∆ W int is the 'internal' work of skin, and ∆W c is the work of contact force. This distinction is made for a more comprehensive global energy balance. Gathering eqs. (3.16) and (3.20) gives the energy balance for the whole system:

∆E k = ∆ Wint + ∆ W int + ∆W c (3.23)
The energy conservation depends then on ∆W c . Two conditions make it equal to zero:

• r n+ 3 2 = 0 ⇔ δ u n+1 = 0: the skin is in rest state; • ( ˙ u n+ 3 2 + ˙ u n+ 1 2
) = 0: the contact is active for previous and current time-step.

For the impacting bar problem with an uniform velocity and a null displacement as initial conditions, the skin works only at release and the time-steps after. Indeed, the release is detected thanks to r n+ 3 2 meaning that δ u 0: the skin is stretched. During the time-steps after the release, the skin goes back to its rest state and then δ u = 0. The corresponding work is not null and not necessary dissipative. If this energetic behaviour might seem problematic, a particular case is particularly interesting. Indeed if the skin does not reach a stretched state, the contact work is conservative.

The analysis of energy balance confirms the idea above: the massless contact node makes the impact conservative. And if the release becomes in general not conservative, energy conservation can still be achieved.

Comparison between singular and consistent contact conditions

The contact law described by eqs. (3.10) and (3.11) is close from the Moreau-Jean's velocity conditions used in the consistent CD-Lagrange scheme. If the contact have a mass, they are written with the notations of this section:

If g n+1 > 0, r n+ 3 2 = 0 If g n+1 0, 0 r n+ 3 2 ⊥ ˙ u n+ 3 2 0 r n+ 3 2 is a Lagrange multiplier, which imposes ˙ u n+ 3 2
according to a discrete persistency condition. Two cases are distinguished. is set to zero because the system already meets the contact conditions. The skin is no more constrained.

The velocity contact law for singular mass matrix follows the same principles: a compressed skin maintains the contact by stopping the contact node, and a non-compressed skin allows the node to leave contact. But the consistent conditions have access to ˙ u n+ 3 2 through the dynamic before to set r n+ 3 2 . At the release time-step for the consistent scheme, the dynamics gives a positive velocity (without contact actions) which sets the impulse:

˙ u free n+ 3 2 > 0 ⇒ r n+ 3 2 = 0
r is then positive for all discrete times. But for the singular scheme, the impulse is calculated first and sets the velocity through the contact law:

r n+ 3 2 0 ⇒ ˙ u n+ 3 2 > 0
This violates the positivity of the contact stresses, but it is necessary as no dynamical link exists between r n+ 3 2 and ˙ u n+ 3 2 in the singular case.

A continuous form of contact law ? The contact law of eqs. (3.10) and (3.11) is build for matching the motion at contact for the singular and the consistent discrete systems. But an interesting question is the form of this law on the continuous, or on the semi-discrete systems. The discussion is made here for the 1D impacting bar.

The key feature of the contact law is to approximate the skin velocity ˙ u by the closest node velocity u. For the continuous system, the velocity of the "closest node" is the velocity of the skin. Indeed the space is continuous. In this case, the contact law has then no sense. The velocity of the skin is determined by the dynamics, in the space of solutions which respect the Signorini's conditions.

For the discrete-in-space system, the skin is the one described on fig. 3.1. But as the time is continuous, the contact law would be transformed into:

If g > 0, ˙ u = u (3.24) Else g = 0, If r > 0, ˙ u = 0 Else r = 0, ˙ u = ( u) + (3.25)
The eq. (3.24) expresses the no contact state, where the skin does not deform. As a consequence, r is null only if there is no residual deformations in the skin.

The eq. (3.25) governs the contact phase, where three cases appears successively:

1. the impact where g = 0, r = 0 as the skin is at rest, and the contact node is stopped by ˙ u = 0 as u < 0;

2. the contact where g = 0, r > 0 as the skin is compressed, and ˙ u = 0 for keeping the massless node in contact;

3. the release where g = 0, r = 0 as the skin is no longer in compression but in its rest state, and ˙ u > 0 as u > 0 (the bar is leaving the contact).

This time-continuous contact law does not accept a skin deformation outside the contact phase. This directs toward a skin which would be numerically added to the initial mesh to be active only for the contact phases. Indeed if the skin is a part of the initial mesh, it should accept residual deformations outside the contact. In the discrete system, the time is not continuous and then the release can happen between two discrete times. This is why a negative r is allowed at release and the time-steps after (where the skin returns to its rest state). This temporary violation of the positivity of contact constraints can be seen as the symmetric of the violation of g > 0 in the discrete persistency condition.

Numerical results on impacting bar with a singular mass

General results

The impacting bar described by fig. 3.1 is set with the values of table 3.2. The time-step h used in simulations is set to 0.9 × h CFL . The discrete solution is computed with the singular formulation of the CD-Lagrange time-integrator, and compared to the analytic solution [see 28] at the contact node. On fig. 3.2, the position during contact is stable and the discrete solution stays close from the analytical one. The residual penetration is still visible on the discrete solution. Indeed as for the CD-Lagrange with a consistent mass, the contact node is stopped once the gap is negative.

The contact node velocity is shown on fig. 3.3. It is similar to the discrete velocity got with a consistent mass matrix. At impact, it takes a null value in one time-step when the gap becomes negative. During contact, it stays at zero meaning the contact node is stopped. And at release, it grows rapidly to positive oscillating values. As mentioned in section 2.4, these spurious oscillations after release are due to the spatial discretization. They happen with all schemes on the impacting bar. These numerical results demonstrates The impulse on fig. 3.4 stabilizes rapidly at the analytical value after impact. At release, it reaches a negative value before coming back to zero. This negative value is where the release is detected, and the only difference with the consistent contact conditions.

The fig. 3.5 depicts the energy balance. In order to show the global transfers, the internal energy is the sum of these of bulk and skin. The system energy is the sum of all conservative terms:

E sys = E k + E int + E int .
This graph confirms the energetic behaviour. The impact is conservative, but after the release energy is dissipated when the skin comes back at rest.

The impacting bar is enough simple to set the time-step for matching a discrete time and the release. Doing this, the skin does not reach a stretched state and the scheme should be conservative. Figures 3.6 and 3.7 compare the previous discrete solution (time-step h) and a solution with the release on a discrete time (time-step h release ). The time-step h release is constant but chosen for matching a discrete time and the release. The impulse, on fig. 3.7, with h release is never negative. At release it is exactly equal to zero, where the impulse with h reaches a negative minimum. The consequence on the system energy is visible on fig. 3.6: with h release , the energy is conserved. The singular mass formulation allows then to get the energy conservation on the discrete system. Here setting the time-step to h release is only made in an illustration purpose. Indeed for a more complex problem, looking for h release would be too time-consuming. And an adaptation of the time-step during the simulation would transform the scheme in an event-driven one.

Nevertheless, the impacting bar demonstrates that the singular mass formulation can make the CD-Lagrange energy conservative. Moreover the singular CD-Lagrange keeps the suitable performances for non-smooth dynamics with an high contact stability. Even if the energy conservation is achieved under a strong condition on the release time, it is a major advance for explicit schemes in non-smooth dynamics.

Influence of skin rigidity

The skin is determined by two numerical parameters: l 0 , the initial length; and k the rigidity. On the preceding impacting bar problem, they are chosen thanks to the initial consistent problem. But the impacting bar problem could be build differently. Instead of cancelling the mass of contact node, a skin with a massless node could be added at the end of the bar. This construction requires a choice in l 0 and k. The following is devoted to analyse the influence of these parameters on the discrete solution.

The initial length l 0 has no influence on the discrete solution. The role of l 0 is only to give a physical meaning to the skin by representing it in space. Indeed l 0 does not appear in equations, which depend only on δ u the difference of displacements between γ C and γ B . The only role of l 0 is to ensure that the skin keeps a positive length l. But the expression of f stays valid even for a l < 0.

This discussion about l 0 raises nevertheless an issue for the gap. In the preceding impact bar case, it was calculated between γ C and the rigid boundary γ R . But the position of γ C depends on l 0 . This current gap definition is then not correct, and discussed later on 3D models.

The relevant skin parameter is the skin rigidity k. Firstly, it impacts the stable time-step. On the impacting bar, as the elements are all equal, the stiffness is constant over them and denoted k nominal . If k > k nominal , the stable time-step is decreased. k nominal will be then taken as a superior limit. In the following, several discrete solutions are computed with different k.

The figs. 3.8b and 3.8c shows that the skin response is smoother for a smaller k: the impulse reaches later the stable values after impact and release on fig. 3.8c; and less oscillations happen for velocity on fig. 3.8b. On fig. 3.8c, the maximal impulse during contact is the same for all k. Indeed the force applied by the bulk on the skin does not change. As a consequence on fig. 3.8d, δ u is larger for a smaller k.

On fig. 3.8a, the positions show an influence of k on the release time. k is then a parameter which influences the energy balance: it can bring closer a discrete time and the release. Instead of adjusting the time-step for the whole simulation as done before, k is far more interesting. It is a local parameter which influence directly the energy balance.

In order to estimate the convergence rate, the error is computed as the distance between graphs of the discrete and analytical solutions thanks to the norm [see 1]:

e = 1 2 h i∈[1,n] |f i -f (t i )| (3.26)
With:

• n is the number of discrete times in interval [t 0 , t f ];

• f i is the discrete solution at time t i , here at contact node;

0 1 2 3 10 -4 -1 0 1 2 3 4 5 6 7 8 9 10 -4
(a) Contact node -Position • f (t i ) is the analytical solution at time t i , here at contact node.

0 1 2 3 10 -4 -6 -4 -
The space-time convergence is verified by a mesh refinement: the number of nodes are increased in the bar. The time-step is set for each mesh at 90% of the stable one. The convergence rate on fig. 3.9 is in O(h) for all k, both on position or impulse. This convergence rate is the same as in chapter 2 for cases involving impacts.

Conclusions on 1D formulation

The singular mass formulation can make the CD-Lagrange fully energy conservative for non-smooth dynamics without any degradation of its performances. The massless contact node cancels the inertia of contact boundary, and then the energy loss at impact. The drawback is that the contact velocity is no more determined by the dynamics. This is overcome by a contact law based on skin stress, which enforces a persistency like contact condition. The skin rigidity becomes a numerical parameter which influences the release time without changing the convergence-rate or the contact stability. This extra-parameter makes possible a conservative energy balance, when the release time matches a discrete time. Such a result is an achievement in explicit time-integrators for nonsmooth dynamics. The numerical results of chapter 2 demonstrates that the CD-Lagrange or the Paoli-Schatzman schemes are really accurate. But the energy balance can not be conservative except at space-time convergence. Indeed the energy loss is linked to the inertia of contact node.

In 1D it is quite simple to adjust k to match the release with a discrete time. For a multiple DOF contact the release time is no more unique, which makes difficult to have each releases on a discrete time. As a consequence the goal of the following section is not to achieve energy conservation for multiple DOFs contact. The 1D singular mass CD-Lagrange is extended to a 3D case only for improving the energy balance.

Extension to 3D by normal massless elements

The normal massless elements formulation

This section presents an extension of the 1D formulation designed on the impacting bar problem in section 3.2 to 3D meshes. In impacting bar problem, the skin could be a 1D element added along the normal on the contact boundary. Indeed in this 1D problem, the outer normal to contact boundary is the bar direction. This interpretation leads to a straightforward extension to 3D meshes. The fig. 3.10 describes the contact boundary of a 3D meshed body. The spatial discretization is done by P1 finite elements: a node corresponds then to a DOF. The contact boundary γ B aims to impact a rigid motionless frontier Γ R . In order to ensure the contact, a "skin" is build on γ B : 1D elements are added on each node, along -n the outer normal to γ B . They link the node on γ B and a massless node. Each element is characterized only by a numerical rigidity. The element direction is chosen as -n in order to be collinear to the outer normal to Γ R when contact happens. The persistency condition is applied on massless nodes, which form a discrete set named γ C . The "skin" gathered now all these normal elements and corresponds directly to the skin on the impacting bar, but for each node on γ B . Remark 9. The skin elements have here no length. If a direction is associated, it gives only the direction for applying the contact force. In a sense, the skin is only a numerical artefact.

The "bulk" refers to the initial meshed body. It is composed by the mass nodes which have three DOFs for each components of the displacement field. The material law is not necessary elastic; indeed the following formulation does not depend on it. The same notations as in section 3.2.1 describe the bulk:

• U, U, Ü: vectors gathering the DOFs for displacements, velocities and accelerations of bulk nodes, X 0 is the initial position;

• M, F ext , F int (U): lumped mass matrix, external and internal forces;

• L: projection operator on each node of γ B along n;

• f : a vector of size the number of nodes in γ B with the values of contact forces (along n) for each node.

The semi-discrete in space bulk dynamics is then:

M Ü = F ext -F int (U) + Lt f (3.27)
As f , the scalar normal quantities on γ B are gathered in vectors. With i a subscript identifying a node, they are:

• ( k) i : numerical rigidity of a normal element;

• ( u) i : displacement of a massless node;

• (ū) i : normal displacement for a bulk node, ū = LU;

• (δ u) i = ( uū) i : difference in displacements between the skin nodes and bulk ones;

• ( f ) i = ( k) i (δ u) i : normal force applied on the bulk.
Similarly to 1D case, the equations of skin are gathered in:

f = k • δ u (3.28)
• denotes the Hadamard product (element-wise multiplication).

The time-integration is still performed with the CD-Lagrange, which adds two time-integration relations and gives the discrete in time forms of eqs. (3.27) and (3.28):

U n+1 = U n + h Un+ 1 2 (3.29) M( Un+ 3 2 -Un+ 1 2 ) = hF(U n+1 ) + Lt r n+ 3 2 (3.30)
u n+1 = u n + h ˙ u n+ 1 2 (3.31)
r n+ 3 2 = h k • δ u n+1 (3.32)
Equations (3.29) to (3.32) describe the bulk and skin system, but do not determine the velocities of massless nodes. This is done by a similar contact law to 1D case, but for each node of γ B :

∀i ∈ γ B , If (g n+1 ) i > 0, ( ˙ u n+ 3
2

) i = ( ufree n+ 3
2

) i (3.33) Else (g n+1 ) i 0,    If ( r n+ 3 2 ) i 0, ( ˙ u n+ 3
2

) i = 0 Else ( r n+ 3 2 ) i < 0, ( ˙ u n+ 3
2

) i = ( ufree n+ 3
2

) + i (3.34) 115 Ufree n+ 3 2
is still the bulk velocity without skin actions given by eq. (3.30):

Ufree

n+ 3 2 = Un+ 1 2 + hM -1 F(U n+1 ) (3.35)
The major difference with the 1D case is the definition of gap. In this 3D formulation, the gap is computed between the nodes of γ B and Γ R by a projection along -n. The definition of the contact normal is clear with this choice: n matches the outer normal of Γ r when contact happens, i.e. g N = 0. But the persistency condition is still applied on γ C which is a set of numerical nodes. On γ B , the persistency condition is then not ensured: the penetration changes during contact, and the normal velocities are not null.

Remark 10. This formulation with a totally numerical massless skin is close to a penalty one. Indeed a force is applied on γ B depending on a kind of penetration δ u and a numerical multiplication coefficient k. But they are crucial differences with a standard penalization method: (i) the skin elements stay active after release; (ii) the nodes of γ C are only numerical; (iii) and the contact is ensured by a persistency like condition.

The algorithm 5 is still valid, as the energy balance (3.23). As mentioned before, the release time is no longer unique but specific to each node of γ B . As these nodes are linked through the dynamics of the bulk, finding a release time for each is no longer a local problem. This makes toughly feasible the exact energy conservation. Nevertheless k forms a set of numerical parameters which can improve the energy balance, and set the skin response without any major changes in the solution.

Numerical results

The impacting dome of section 2.5 is used for a numerical validation of the 3D singular mass formulation with normal massless elements. The initial configuration is described by fig. 3.11a, and the table 3.3 gathers the values which set the case. They are similar to those used in section 2.5: the material is not physical but close from a "compressible" rubber.

γ B is formed by the upper nodes of dome, on the cylindrical part. They get normal massless elements along y . Indeed the assumption is made that the contact normal is -y , which seems valid for the top nodes of dome. The k given in table 3.3 is the value for the nodes inside γ B . For those on edge of γ B , the skin rigidity is k/2. k is here the maximal diagonal value of LK Lt , with K the rigidity matrix.

On γ B , a central zone is defined and shown on fig. 3.11b. The impulse is averaged on it with coefficients depending on the position of the node: 1/4 for a node on a corner, 1/2 on a edge, 1 inside. For averaging the displacement and velocity, a standard mean is computed for the highest nodes on the dome in the central zone. The fig. 3.12 represents the deformed mesh for different times. The body is coloured according to the magnitude or r. The numerical skin is not depicted on the graphs. The nodes at the top of the dome are those of γ B . And even if the persistency condition is not satisfied on γ B , the penetration on fig. 3.12c at maximum impulse is quite similar to the one at impact on fig. 3.12a or at release on fig. 3.12d.

In order to visualize the position of γ C , an arbitrary initial length is set to the normal massless elements. The fig. 3.13a depicts the averaged x and x. The maximal penetration stays low: x is not bigger than three times the penetration at impact. But oscillations are visible for x. They are much more visible on u in fig. 3.13b with a large amplitude during contact. Despite these oscillations on velocity, the weighted impulse on central contact zone on fig. 3.13c is quite stable. The normal massless elements decrease then the contact stability of γ B . Indeed the persistency condition is applied on γ C , 117 Here the skin rigidity is arbitrary chosen without any regards to an optimization of the skin response. k being the maximal diagonal value of LK Lt , it can explain the stiff behaviour of skin. The fig. 3.13d compare the conservative terms of energy balance for two discrete solutions: one computed with the normal massless elements formulation, and one with the standard CD-Lagrange. The loss of is 0.3% of the initial energy for the CD-Lagrange with normal massless elements, and 0.6% for the standard CD-Lagrange. Despite the energy losses are small, the singular mass formulation decreases it by half without any optimization on k. This demonstrates its ability to improve the energy balance.

The space-time convergence is checked with two numerical solutions: one as a reference by the consistent CD-Lagrange, and one with normal massless elements. The error is estimated thanks to the norm (3.26) for the nodes of γ B in the central zone. The mesh refinement is similar to section 2.5. The error on fig. 3.14 shows a convergence in O(h) for both the displacement and the impulse. This demonstrates that the normal massless elements solution converges to the consistent one.

Conclusion on normal massless elements

The extension of the 1D formulation to 3D meshes by normal massless elements is quite interesting. Even if it decreases the contact stability compared to the standard CD-Lagrange, it improves the energy balance. Moreover the skin response can be optimized thanks to the k numerical parameter. It can be easily implemented on an existing FE software as it consists only in adding numerical parameters for the contact nodes. Moreover it implies no modifications on the mesh, or in the assembly of the matrices. It is too fully compatible with non-linear materials or large deformations, as it is localized on the contact boundary.

Extension to 3D by the singular mass skin

The singular mass skin formulation

The existing singular mass formulations are quite different of the preceding one with normal massless elements. Indeed the massless nodes are present in the initial mesh, and not added on the contact boundary. Two main singular formulations can be distinguished:

• with a zero mass entry only for the DOFs in the normal direction of contact as proposed in [START_REF] Khenous | Mass redistribution method for finite element contact problems in elastodynamics[END_REF];

• with a zero mass entry for all the DOFs associated to the contact nodes as proposed in [START_REF] Wohlmuth | Variationally consistent discretization schemes and numerical algorithms for contact problems[END_REF].

The first approach is in practice more suitable to small deformations, where the contact normal is constant. The mass matrix modifications are then done once in the initial configuration. In this section, such a formulation is introduced for the CD-Lagrange with first numerical results. The problem is described by fig. 3.15. This time a layer of 3D elements is added on γ B , the contact boundary of the initial mesh. This skin must follow some specific properties:

1. the skin elements are build by a translation of the 2D elements of γ B , along -n with a distance l 0 ;

2. small deformations must be assumed for keeping the deformed configuration equal to the initial one.

A skin with these properties allows to associate directly a node of γ C to a node of γ B by the mean of a projection along n. As the configuration is near constant over the simulation, this association is the same for all discrete times. The mesh is now constituted by the initial one and the added skin. A standard assembly is done over the all mesh in order to obtain the rigidity matrix K, and a lumped mass matrix. Then the mass entries are set to zero for the DOFs along n of the contacting nodes of γ C . The resulting matrix is denoted M. Not all nodes of γ C have necessarily a singular mass entry. Only those in contact during the simulation get a null mass entry along n. The computation of the singular mass matrix is facilitated if n is one of the coordinates directions.

The singular formulation is directly introduced on the space-and-time discrete system. The material is assumed to be elastic and linear. And the DOFs are gathered in global vectors: U for the displacement, U for the velocities. This time, they concern both the DOFs of the bulk and the skin. New notations are introduced to describe the problem:

• L, the operator which projects the global quantities along n on γ C for the nodes with a singular mass entry (the contacting ones);

• L, the operator which projects the global quantities along n on γ B for the nodes associated to those of L;

• B, the operator which select the complementary DOFs of L.

L and L are related by their numbering. A same row in the vectors LU and LU corresponds to the two nodes of γ B and γ C related by the projection along n. B selects all the mass DOFs: all DOFs of the initial mesh and the tangential ones on γ C . The displacement and velocities are decomposed between the mass and massless DOFs:

Ū = BU u = LU U = B t Ū + L t u (3.36)
M, the singular mass matrix, can be decomposed in two sub-matrices:

M = BMB t + LM L t (3.37)
with:

• BMB t the consistent mass matrix for the mass DOFs, diagonal as M is lumped;

• LM L t a matrix with only the null entries of the massless DOFs on γ C .

As small deformations are assumed, this decomposition is valid during the whole simulation. The dynamics of the system is:

M Un+ 3 2 -Un+ 1 2 = hF(U n+1 ) + L t r n+ 3 2 (3.38)
The eq. (3.38) is decomposed thanks to B and L between the mass and the massless DOFs:

[BMB] t Un+ 3 2 -Un+ 1 2 = hBF(U n+1 ) (3.39) 0 = h LF(U n+1 ) + r n+ 3 2 (3.40)
The eq. (3.39) is the dynamic for the mass DOFs, and the eq. (3.40) for the massless ones. In eqs. (3.39) and (3.40) the velocities ˙ u n+ 3 2 are not determined. A contact law is then introduced. This time the free-velocity is obtained in a slight different way by: Ufree

n+ 3 2 = Un+ 1 2 + [BMB t ] -1 hBF(U n+1 ) - Lt r n+ 3 2 (3.41) ūfree n+ 3 2 = L Ufree n+ 3 2 (3.42)
The eq. (3.41) is the velocity of mass DOFs, with a correction which retires the contact impulse on γ B (initially applied on γ C ). This formula is based on the static feature of the skin and its particular building. Indeed the normal stresses on γ C are opposite to those on γ B . The corresponding numbering between L and L allows to cancel the action of the skin on γ B by using r.

With the definition of eq. (3.42), the contact law is similar to the preceding ones. But it concerns only the nodes with a singular entry on γ C , denoted γ Ac.

C :

∀i ∈ γ Ac. C , If (g n+1 ) i > 0, ( ˙ u n+ 3
2

) i = ( ufree n+ 3 2 ) i (3.43)
Else (g n+1 ) i 0,

   If ( r n+ 3 2 ) i 0, ( ˙ u n+ 3
2

) i = 0 Else ( r n+ 3
2

) i < 0, ( ˙ u n+ 3
2

) i = ( ufree n+ 3
2

) + i (3.44)
This time the gap is computed on γ C as the geometrical one.

Remark 11. It is necessary to restrain the massless DOFs to the nodes in contact during the simulation. Indeed the eqs. (3.43) and (3.44) does not ensure a null impulse for a node not in contact. r is the normal part of internal stresses on γ C which are not necessary null for a non-contact phase. This is a major difference with the normal massless elements formulation where r tends to zero for a non-contact phase.

Both the energy balance of eq. (3.23) and the algorithm 5 are still valid. Following the remark above, the contact work is not necessarily null for a non-contact phase.

Numerical results

This formulation is tested on the case described by fig. 3.16a and set with the values of table 3.4. The body measures 1 m along x and z , and 0.5 m along y . The upper layer of element is the skin. It is build along y , which is supposed to be the contact normal, with a translation distance of 3.3 × 10 -2 m. Only the nodes depicted in red in fig. 3.16b get a singular entry along y . The material is elastic and set as aluminium. The small deformations are 123 m above the body. The discrete quantities are averaged for the singular nodes between z = 1/3 m and z = 1/6 m. The averaged impulse integrates weighting coefficients depending on the node position: 1/2 on a side, 1/4 at a corner. These averaged quantities are depicted on fig. 3.17.

The displacement on fig. 3.17a and the velocities on fig. 3.17b show that the persistency condition is met with a null velocity during contact. For the impulse on fig. 3.17c, oscillations happens during contact but the stability stays acceptable. But after release, the impulse is not equal to zero. It oscillates around small values. The consequence on the energy balance is visible on fig. 3.17d. Due the non-null impulses, the contact work increases the system energy after release.

Conclusion on singular mass skin formulation

On the contrary to the formulation with normal massless elements of section 3.3, the singular mass skin formulation is less robust. It is limited to small perturbations, and requires a mesh with specific properties for the skin. Moreover the results are less interesting as energy is injected after release, and for a long time. This could lead to stability issues, and decreases the accuracy.

The singular mass skin formulation underlines the real difficulty for a singular mass formulation in the CD-Lagrange scheme. As the velocities are undetermined for the massless DOFs in the dynamics, they have to be fixed by the contact law. And this contact law demonstrates limitations for the non-contact phases. In the normal massless elements formulation, the velocities of massless DOFs are only numerical and local to contact nodes. Their velocities have then less influence on the global solution. But with a singular mass skin, the massless DOFs are integrated in the 3D mesh which links them deeply to the mass DOFs. Their velocities have an influence far more important, which reduces the accuracy of the scheme.

Conclusion: singular mass in the CD-Lagrange

The 1D formulation of a singular CD-Lagrange achieves the energy conservation. If this result is conditioned to the release on a discrete time, it is a new result for non-smooth dynamics particularly for schemes based on velocity contact conditions. If the energy conservation is hardly feasible for multi-DOFs contact, the formulation introduces a local parameter, the skin rigidity. It allows to set the skin response and optimize the energetic behaviour. The undetermined velocities for massless DOF are set thanks to a contact law, which keeps the contact highly stable and close from a persistency condition.

The 1D formulation is extended to 3D meshes in a penalty like formulation, the normal massless elements formulation. If the contact stability is decreased, the energy balance is more accurate than with the consistent CD-Lagrange. And as the formulation is local to the contact nodes, it is fully compatible with large deformations or non-linear materials with an easy implementation. The normal massless elements seems them a promising way to use the singular mass in the CD-Lagrange scheme.

On the contrary, the singular mass skin is less performing. This 3D extension is closer from the existing singular formulations for implicit schemes. But the energetic behaviour is degraded and the formulation seems nonapplicable in large deformations. These extra-difficulties come from the undetermined velocities for the massless DOFs. Indeed the inertia is central in the CD-Lagrange, as the time-integration link between two time-steps. A singular mass formulation directly in the initial mesh seems less suitable to explicit schemes.

Chapter 4

The deformable-deformable contact: mortar methods in the CD-Lagrange scheme This chapter presents a CD-Lagrange scheme with a contact enforced with a mortar method. The goal is to handle a deformable-deformable contact with large sliding. The first part introduces the mortar methods in the CD-Lagrange scheme. It looks at the new contact problem involved, and an algorithm with an acceleration technique i proposed to solve it. The second part numerically validates the preceding formulation with the acceleration technique on a case with small-sliding and frictionless contact. The last part presents numerical cases closer from an industrial case, with large sliding and friction.

Discretization of contact between two deformable bodies

Continuous description of contact

In this section, the contact is enforced between two deformable bodies as described by fig. 4.1. Each body, denoted by the superscript (i), is deformable and described as in section section 1.1. A deformable problem is associated to each body, and each problem is governed by the eqs. (1.9) to (1.12). But now the boundaries ∂Ω (i) 0 of each body are split into three complementary sets:

Γ (i) D ∪ Γ (i) N ∪ Γ (i) C = ∂Ω (i) 0 Γ (i) D ∩ Γ (i) N = Γ (i) D ∩ Γ (i) C = Γ (i) N ∩ Γ (i) C = ∅ (4.1) Γ (1) C and Γ (2)
C denote the parts where contact can potentially happen. The counterparts of C the contact stress is denoted:

Γ (i) C in deformed configuration are γ (i) C . n Γ C,1 Ω 1 (t) Γ D,1 Γ N,1 Γ C,2 (t) Ω 2 (t) Γ D,2 Γ N,2 (a) Initial configuration n γ C,1 (t) Ω 1 (t) Γ D,1 (t) Γ N,1 (t) γ C,2 (t) Ω 2 (t) Γ D,2 (t) Γ N,2 (t) (b) Deformed configuration
λ (i) = σ (i) n (i) (4.2)
with n (i) the inward normal to γ

C , and σ (i) the Cauchy stress tensor. The contact stress (4.2) is decomposed into a normal and a tangential part:

λ (i) = λ (i) N n (i) + λ (i)
T t (i) with λ

(i) N , λ (i) T 0 (4.3)
t (i) is in the orthogonal plane to n (i) . As mentioned in [START_REF] Yang | Two dimensional mortar contact methods for large deformation frictional sliding[END_REF], the actionreaction law links the two contact boundaries by:

∀x (1) ∈ γ (1)
C , λ (1) (x (1) )dγ

C = -λ (2) π 21 (x (1) ) dγ The contact quantities introduced in section 1.2 are still valid. The relative gap and velocities characterize the kinematics of γ C : ) )n (4.7)

g(x (1) ) = x (1) -x (2) (π 21 (x (1) )) • n (4.5) v N (x (1) ) = u(1) (x (1) ) -u(2) π 21 (x (1) ) • n (4.6) v T (x (1) ) = [ u(1) -u(2) ] -v N (x ( 1
Here the contact normal n has still to be chosen, either as the inner normal of γ

C or as the outer one of γ

C . The Signorini's conditions (1.51) describe the normal part of contact quantities, and the Coulomb's law (1.57) the tangential part: ∀(x (1) , t) ∈γ

(1) C × [0, t f ], 0 g(x (1) ) ⊥ λ (1) N 0 (4.8) If λ (1)
N = 0, λ (1)
T = 0 If λ (1)
N > 0, 0 (µλ

(1) N -λ (1)
T ) ⊥ v T (x (1) ) 0 (4.9)

With λ

T = -αv T (x (1) ), α ∈ R +

Here the body (1) is chosen for expressing the contact: the contact quantities are related to γ

C . In the following, the superscript (1) is omitted when the quantity is relative to γ Remark 12. Following [START_REF] Popp | Mortar Methods for Computational Contact Mechanics and General Interface Problems[END_REF][START_REF] Yang | Two dimensional mortar contact methods for large deformation frictional sliding[END_REF], γ

C is the slave or non-mortar surface. In contrast, γ

C is the master or mortar surface. Indeed in this manuscript, the discrete interface for expressing the mortar quantities is chosen between γ As stated in [START_REF] Popp | Mortar Methods for Computational Contact Mechanics and General Interface Problems[END_REF][START_REF] Yang | Two dimensional mortar contact methods for large deformation frictional sliding[END_REF], in order to enforce the contact in a mortar sense, the following term is added in the PVW (1.16) defined in section 1.1: (1) (x (1) ) -δu (2) π 21 (x (1) ) dγ C (4.10) This term is the summation of virtual contact works on both bodies:

δW C = γ C λ • δu
δW C = γ (1) C
λ (1) • δu (1) dγ

(1) C + γ (2) C
λ (2) • δu (2) dγ

(2) C

The integral over γ

C is expressed on γ

C thanks to the action-reaction principle (4.4) and the projection π 21 . If the integrals have different supports γ C , the Signorini's conditions (4.8) ensure that where g(x (1) ) > 0 both λ (1) (x (1) ) and λ (2) π 21 (x (1) ) are equal to zero supposing π 21 is surjective. But to be mathematically precise the integrals should be replaced by duality pairing as mentioned in [START_REF] Popp | Mortar Methods for Computational Contact Mechanics and General Interface Problems[END_REF][START_REF] Wohlmuth | A Mortar Finite Element Method Using Dual Spaces for the Lagrange Multiplier[END_REF].

The functional spaces are those of section 1.1:

• δu ∈ V, with V = {δu ∈ H 1 (Ω)|δu = 0 on Γ D }; • u ∈ U, with U = {u ∈ H 1 (Ω)|δu = û on Γ D };
• δλ ∈ M(λ), with M(λ) the space for contact stresses which enforce the contact conditions of eqs. (4.8) and (4.9).

M(λ) must be included in the dual space M = H -1/2 (γ C ) of the trace space W = H 1/2 (γ C ), the restriction of V at γ C . But in the discrete setting, M is often chosen as W [see 67 , 78]. A distinction may exist between the vectorvalued and scalar-valued spaces M, but it is not assumed here to simplify the notations.

To complete the problem, the contact conditions are stated under the weak form [see 67]:

∀δλ N , δλ T ∈ M(λ), γ C (δλ N -λ N )g N dγ C + γ C (δλ T -λ T )v T dγ C 0 (4.11)
with λ T = λ T t. The equivalence between the weak form (4.11), and the strong contact conditions of eqs. (4.8) and (4.9) is proven in [START_REF] Wohlmuth | Variationally consistent discretization schemes and numerical algorithms for contact problems[END_REF].

Spatial discretization by mortar methods

The mortar methods

The preceding contact terms are discretized thanks to the mortar methods. They are are not specific to contact problems, but relate two non-conforming meshes by a projection in a weak sense. These methods are very robust for non-conforming meshes.

They were introduced in [START_REF] Belgacem | The mortar finite element method for contact problems[END_REF][START_REF] Belgacem | The Mortar finite element method with Lagrange multipliers[END_REF] for contact problems, where they are applied on the Lagrange multipliers for contact enforcement. In [START_REF] Wohlmuth | A Mortar Finite Element Method Using Dual Spaces for the Lagrange Multiplier[END_REF], an other functional space is proposed for the Lagrange multipliers: the dual Lagrange multipliers space. It allows to demonstrate a strict equivalence between the continuous contact constraints and the discrete counterparts. The mortar methods for contact problems were then developed for example in [START_REF] Puso | A mortar segment-to-segment contact method for large deformation solid mechanics[END_REF][START_REF] Yang | Two dimensional mortar contact methods for large deformation frictional sliding[END_REF], or in the thesis works in [START_REF] Hüeber | Discretization techniques and efficient algorithms for contact problems[END_REF][START_REF] Popp | Mortar Methods for Computational Contact Mechanics and General Interface Problems[END_REF]. A clear formalisation of mortar methods is also proposed in [START_REF] Wohlmuth | Variationally consistent discretization schemes and numerical algorithms for contact problems[END_REF]. All these works are devoted to extend the mortar methods to frictional contact, and for 2D or 3D meshes. They all concern implicit time-integration scheme, with as a major concern the linearisation of mortar operators. Indeed it is necessary in the non-linear solver.

Few works exist for explicit time-integration scheme. In [START_REF] Casadei | A mortar spectral/finite element method for complex 2D and 3D elastodynamic problems[END_REF], Casadei et al. use the mortar methods together with the central difference method but not for a contact problem. They relate two meshes, one with finite elements and one with spectral elements.

Discrete-in-space mortar formulation of contact

As done in section 1.1, the continuous functional spaces are approximated by a Galerkin method. The term δW C is approximated with discrete spaces defined on γ C . They arises from those defined on all the domain Ω h , and restricted to γ C . With U h and V h the approximations of U and V, these spaces are:

• X (i) h , the restriction of U h to γ (i) C ; • W (i) h , the restriction of V h to γ (i) C ;
• M h , the space of contact stresses (without any regard to contact constraints).

As explained in [START_REF] Popp | Mortar Methods for Computational Contact Mechanics and General Interface Problems[END_REF], two main choices are reported in literature for M h . The first one is to take M h ≡ W

(1) h

[see [START_REF] Belgacem | The Mortar finite element method with Lagrange multipliers[END_REF][START_REF] Wriggers | Computational Contact Mechanics 2nd[END_REF]. This approach is quite standard in FE software, because the shape functions are those used for the displacement. The second one, called dual Lagrange multipliers and introduced in [START_REF] Wohlmuth | A Mortar Finite Element Method Using Dual Spaces for the Lagrange Multiplier[END_REF], builds M h thanks to a bi-orthogonality condition to W h . They present two major advantage. First they are numerically more efficient because one of the two mortar matrices is diagonal. And they allow a strict equivalence between continuous and discrete mortar conditions.

In the following, the standard approach is selected. The Lagrange multipliers are approximated by the shape functions of displacement restricted to γ (i)

C . The bases of spaces W (i) are denoted φ

(i) j j∈[1..N i ]
with N i the number of interpolation points in γ

(i)
C , leading to the discrete fields for test functions:

δu (1) = i∈[1..N 1 ]
φ (1) (x (1) )δu

(1) i (4.12)

δu (2) = i∈[1..N 2 ] φ (2) (x (2) )δu (2) i (4.13) λ = i∈[1..N 1 ]
φ (1) (x (1) )λ i (4.14)

As in section 1.2, the finite elements are P1. δu

(i)
i and δλ i are then directly the nodal values. Following [START_REF] Yang | Two dimensional mortar contact methods for large deformation frictional sliding[END_REF], the vector valued λ is split into a normal and a tangential parts which are scalars:

λ i = λ N,i n i + λ T,i t i (4.15)
With n i , an inward normal vector at point i; and t i , a tangential one in the tangent plane to n i . Both n i and t i are supposed to be known at all nodes of γ C .

The virtual contact work (4.10) is then split into a normal and a tangential part: (1) )δu (1) (x (1) (1) )δu (2) π 21 (x (1) (1) )δu (1) (x (1) (1) )δu (2) π 21 (x (1) ) dγ C (4.16)

∀δu (1) ∈ W (1) h , δu (2) ∈ W (2) h , δW C = j∈[1..N 1 ] λ N,j n j • γ C φ (1) j (x
) dγ C - γ C φ (1) j (x
) dγ C + j∈[1..N 1 ] λ T,j t j • γ C φ (1) j (x
) dγ C - γ C φ (1) j (x
In order to get the matrix form, the test functions on displacement δu (i) range in the following base of vector-valued space W h :

S (i) = φ (i) i δ ik , i ∈ [1..N (i) ], k ∈ [1..3] with δ ik = [δ 1k δ 2k δ 3k ] t (4.17)
This leads to:

∀ (φ (i) i δ ik ) ∈ S (i) , δW C (φ (i) i δ ik ) = j∈[1..N 1 ] λ N,j γ C φ (1)
j (x (1) )φ

i (x (1) ) dγ C -

γ C φ (1)
j (x (1) )φ

(2) i π 21 (x (1) ) dγ

C n j • δ ik + j∈[1..N 1 ] λ T,j γ C φ (1) j (x (1) )φ (1)
i (x (1) ) dγ

C - γ C φ (1)
j (x (1) )φ

(2) i π 21 (x (1) ) dγ C t j • δ ik (4.18) and with a matrix form: C . B N,1 and B T,1 are then square with a shape of a mass matrix, but B N,12 and B T,12 have a rectangular shape. The expressions of the sub-matrices associated to the interpolation point ij are:

[δW C,i ] i∈[1..N 1 ] = B t N,1 -B t N,12 λ N + B t T,1 -B t T,
[(B N,1 ) ij ] i∈[1..N 1 ] j∈[1..N 1 ] = γ C φ (1)
i (x (1) )φ

(1) j (x (1) ) dγ C n t i (4.20)

[(B N,12 ) ij ] i∈[1..N 1 ] j∈[1..N 2 ] = γ C φ (1)
i (x (1) )φ

(2) j π 21 (x (1) ) dγ C n t i (4.21)

[(B T,1 ) ij ] i∈[1..N 1 ] j∈[1..N 1 ] = γ C φ (1) i (x (1) )φ (1) j (x (1) ) dγ C t t i (4.22) [(B T,12 ) ij ] i∈[1..N 1 ] j∈[1..N 2 ] = γ C φ (1) i (x (1) )φ (2) j π 21 (x (1) ) dγ C t t i (4.23)
Normal and tangential operators are quite similar. The sums of the product φ (2) j are common to B n,12 and B T,12 . The only difference is the multiplication by n i or t i depending on the normal or tangential operator.

The weak Signorini's conditions (4.8) are approximated by the following ones

[see 68] ∀i ∈ [1..N 1 ], 0 (g) i ⊥ (λ N ) i 0 (4.24)
With:

g = [B N,1 -B N,12 ](X) (4.25)
The mortar gap g is a measure of the gap in a weak sense. X gathers the nodal positions x (i) i at times t for both Ω (1) and Ω (2) . These conditions are local to each interpolation points, but enforce the Signorini's conditions in a weak sense on γ C . As proven in [START_REF] Hüeber | Discretization techniques and efficient algorithms for contact problems[END_REF], the point-wise form (4.24) is equivalent to the continuous weak form (4.8) only for dual Lagrange multipliers. If not, it is only an approximation.

These point-wise Signorini's conditions have an equivalent form on velocity [see 61]:

∀i ∈ [1..N 1 ], If (g) i > 0, (λ N ) i = 0 If (g) i 0, 0 (v N ) i ⊥ (λ N ) i 0 (4.26)
v N is the relative normal velocity in a mortar sense. It is computed as for the gap but on the discrete velocities:

v N = [B N,1 -B N,12 ]( U) (4.27)
U gathers the nodal velocities u(i) i at times t for both Ω (1) and Ω (2) . For the weak form of the Coulomb's law (4.9), the discrete point-wise form is:

∀i ∈ [1..N 1 ], If (λ N ) i = 0, (λ T ) i = 0 If (λ N ) i > 0, 0 (µλ N -λ T ) i ⊥ |(v T ) i | 0 (4.28)
With v T defined as v N , but thanks to the tangential mortar operators:

v T = [B T,1 -B T,12 ]( U) (4.29)
The complementary relation involves the absolute value of terms in v T , indeed no assumptions are made on their signs.

Remark 13. Here for non-conforming meshes, g, v N and v T are defined in a weak sense (no more local in space) thanks to the mortar operator. They do not correspond to any physical quantities as in the formulation for conforming meshes in section 1.2.2. As a result, a node geometrically in contact has not necessary a negative mortar gap. And the terms of v N and v T associated to a node do not correspond to the normal and tangential parts of relative velocity ur . During contact the normal relative velocity ur • n is no longer equal to zero and the penetration varies. Indeed the complementary relation does not hold on these quantities but on the mortar ones.

The vectors t i associated to each interpolation point on γ

C are computed thanks to the discrete velocity fields:

u(i) = i∈[1..N i ] φ (i) i u(i) i , u ∈ W (i) h (4.30)
N i is the number of interpolation points in the boundary γ ur (x (1) ) = u(1) x (1) u(2) π 21 (x (1) ) (4.31)

Finally the vectors t i are defined as the opposite tangential direction of the relative velocity at each node of γ

C : ∀i ∈ [1..N 1 ], t i = - ur (x (1)
i ) -ur (x

(1) i ) • n i n i ur (x (1) i ) -ur (x (1) i ) • n i n i (4.32)
with x i the position of the node i.

Finally the space-discrete dynamics, together with the Signorini's conditions under the displacement form, reads:

Md U = F ext -F int (U) dt + B t N dr N + B t T dr T (4.33) ∀i ∈ [1..N 1 ], 0 (g) i ⊥ (λ N ) i 0 (4.34) If (λ N ) i = 0, (λ T ) i = 0 If (λ N ) i > 0, 0 (µλ N -λ T ) i ⊥ |(v T ) i | 0 (4.35)
With dr N and dr T the measures defined in section 1.2 by eq. (1.85).

Using the velocity form of Signorini's conditions, the dynamic equations are:

Md U = F ext -F int (U) dt + B t N dr N + B t T dr T (4.36) ∀i ∈ [1..N 1 ], If (g) i 0, 0 (v N ) i ⊥ (λ N ) i 0 If (g) i > 0, (λ N ) i = 0 (4.37) If (λ N ) i = 0, (λ T ) i = 0 If (λ N ) i > 0, 0 (µλ N -λ T ) i ⊥ |(v T ) i | 0 (4.38)
The mortar operators are:

B N = [B N,1 -B N,12 ] B T = [B T,1 -B T,12 ] (4.39)
They define the gap, and the relative normal and tangential velocities:

g = B N X v N = B N U v T = B T U (4.40)
About the determination of n i . Being the inner normal to γ C at node i, several definitions exist for n i . For a discrete-in-space problem, it is commonly taken as an average on normals of adjacent elements to the node i [see [START_REF] Popp | Mortar Methods for Computational Contact Mechanics and General Interface Problems[END_REF][START_REF] Puso | A mortar segment-to-segment contact method for large deformation solid mechanics[END_REF][START_REF] Yang | Two dimensional mortar contact methods for large deformation frictional sliding[END_REF]. A discrete normal field is then defined thanks to these interpolation values:

n h (x 1 ) = i∈[1..N 1 ] φ i n i (4.41)
According to [START_REF] Popp | Mortar Methods for Computational Contact Mechanics and General Interface Problems[END_REF][START_REF] Puso | A mortar segment-to-segment contact method for large deformation solid mechanics[END_REF][START_REF] Yang | Two dimensional mortar contact methods for large deformation frictional sliding[END_REF], the projection between the discrete γ

C and γ

(2) C

seems more robust if it is based on n h .

Update of n i and t i in the contact solving. In a general framework n i and t i depends on the displacement, which depends also on n i and t i through the contact actions. This causes a space non-linearity in the discrete problem. The discretization of both the virtual contact work (4.10) and the weak contact conditions (4.11) has then to be frame indifferent [START_REF] Yang | Two dimensional mortar contact methods for large deformation frictional sliding[END_REF] in order to integrate this non-linearity in the solving. But here n i and t i are known before the contact solving stage thanks to the time-discretization (see below). This simplifies the spatial discretization, and leads to linear contact problems.

Time-discretization by the CD-Lagrange scheme

The time-integration of the spatial discrete dynamics of eqs. (4.36) to (4.38), is performed with the CD-Lagrange scheme (see the eqs. (1.157) and (1.158)):

U n+1 = U n + h Un+ 1 2 (4.42) M Un+ 3 2 -Un+ 1 2 = h F ext (t n+1 ) -F int (U n+1 ) + B t N,n+1 r N,n+ 3 2 + B t T,n+1 r T,n+ 3 2 (4.43) ∀i ∈ [1..N 1 ], If (g(U n+1 )) i > 0, (r N,n+ 3
2

) i = 0 If (g(U n+1 )) i 0, 0 (v N,n+ 3
2

) i ⊥ (r N,n+ 3
2

) i 0 (4.44) If (r N,n+ 3
2

) i = 0, (r T,n+ 3
2

) i = 0 If (r N,n+ 3
2

) i > 0, 0 µ(r N,n+ 3
2

) i -(r T,n+ 3
2

) i ⊥ |(v T,n+ 3
2

) i | 0 (4.45)
In order to get the reduced dynamics for contact problems, the eq. ( 4.43) is condensed on γ C by a left-multiplication by B N M -1 and B T M -1 . Using B t N B T = 0, this leads to: with the Delassus operators:

H N,n+1 r N,n+ 3 2 = v N,n+ 3 2 -v free N,n+ 3 2 (4.46)
H N = B N,n+1 M -1 B t N,n+1 H T = B T,n+1 M -1 B t T,n+1 (4.48)
The normal and tangential relative free-velocities are defined similarly to v N,n+ 3 The free-velocity is still the velocity without any contact actions:

Ufree n+ 3 2 = Un+ 1 2 + hM -1 F(U n+1 ) (4.50)
As for the rigid-deformable contact in section 1.2.2, the dynamics condensed on γ C together with the contact conditions forms two LCPs. The first one involves only the normal part of contact impulse:

H N,n+1 r N,n+ 3 2 = v N,n+ 3 2 -v free N,n+ 3 2 (4.51) ∀i ∈ [1..N 1 ], If (g(U n+1 )) i > 0, (r N,n+ 3
2

) i = 0 If (g(U n+1 )) i 0, 0 (v N,n+ 3
The second one gives the tangential impulse:

H T,n+1 r T,n+ 3 2 = v T,n+ 3 2 -v free T,n+ 3 2 (4.52) ∀i ∈ [1..N 1 ], If (r N,n+ 3
2

) i = 0, (r T,n+ 3
2

) i = 0 If (r N,n+ 3 2 ) i 0, 0 (µr N,n+ 3 2 -r T,n+ 3 2 ) i ⊥ |(v T,n+ 3 2 ) i | 0
But contrary to the rigid-deformable case, these LCPs are not locally solved at each node. Indeed H N and H T , given by eq. (4.48), are no longer diagonal.

For solving the dynamics (4.43), a crucial point is the computation of B N and B T . B N is determined on the normal to contact boundary in the deformed configuration. The updated displacement is computed once in the time-step due to the explicit feature of the CD-Lagrange. Therefore B N is computed once, which makes the LCP (4.51) linear.

For B T , the crucial quantity is the t i which are the opposite directions of tangential velocities, and then depend on Un+ 3 2 . But, as for the rigiddeformable case, the assumption is made that the directions of relative freevelocities are equal to the directions of the updated relative velocities. As Ufree n+ 3 2 does not depend on contact impulses, B T can be computed once in the time-step, which makes the tangential LCP linear too. In order to compute the {t i } i∈[1..N 1 ] , a free-velocity field is defined upon Ufree (x (1) ) -u(2),free

n+ 3 2 π 21 (x (1) ) (4.54) ∀i ∈ [1..N 1 ], t i = - ufree r,n+ 3 2 (x i ) -ufree r,n+ 3 2 (x i ) • n i n i ufree r,n+ 3 2 (x i ) -ufree r,n+ 3 2 (x i ) • n i n i (4.55)
The energy balance between

[t n+ 1 2 , t n+ 3 2
] comes directly from eq. (1.183) for conforming contact meshes: ensures the same conservation properties as for the rigid-deformable contact: ∆W C,N is conservative everywhere except at impact where it dissipates kinetic energy. When ∆W C,N is null, the contact works on slave and master side are opposite:

1 2 Ut n+ 3 2 M Un+ 3 2 t n+ 3 2 t n+ 1 2 = 1 2 h( Un+ 3 2 + Un+ 1 2 ) t F ext (t n+1 ) -F int (U n+1 , Un+1 ) 1 2 ( Un+ 3 2 + Un+ 1 2 ) t B t N,n+1
∆W (1) C,N = 1 2 ( Un+ 3 2 + Un+ 1 2 ) t B t N,1 r N,n+ 3 2 (4.59) ∆W (2)
C,N = -

1 2 ( Un+ 3 2 + Un+ 1 2 ) t B t N,12 r N,n+ 3 2 (4.60) ∆W C,N = ∆W (1)
C,N + ∆W (2) C,N (4.61)
Algorithm 6 CD-Lagrange for deformable-deformable contact

U 1 2 ← U0 + 1 2 hM -1 F(U 0 ) for (t n+1 , t n+ 3 2 ) ∈ {t n+1 } × {t n+ 3 2 } do U n+1 ← U n + h Un+ 1 2 Displacement update F n+1 ← F(U n+1 , Un+ 1 2 )
Assembly of stresses Ufree

n+ 3 2 ← Un+ 1 2 + hM -1 F n+1 Free-velocity update          B N ← B N (U n+1 ) B T ← B T ( Ufree
← Ufree n+ 3 2 + hM -1 (B t N r N,n+ 3 2 + B t T r T,n+ 3 2
) Contact correction

Un+ 3 2 |Γ D ← UD (t n+ 3 2
)

Dirichlet correction end for

The algorithm 6 is similar to the rigid-deformable case for its main steps. The first step explicitly updates the displacement, which is not corrected after that. Once the displacement is known, the internal stresses can be assembled which allows to compute the free velocity. Then the contact solving step begins with the assembly of mortar operators B N and B T . Unlike in the rigid-deformable problem, a linear solver is required for solving the normal and tangential contact LCPs. The global algorithm 6 is then no more explicit for the contact impulses. Once the normal and tangential impulses are computed, the global velocity is corrected to enforce the contact. The enforcement of Dirichlet's boundary conditions are done then, directly on the velocity.

Solving the contact LCPs

The LCP solver: the non-smooth projected Gauss-Seidel algorithm This section is devoted to the numerical solving of the two LCPs in the CD-Lagrange scheme:

• the LCP (4.51) for the normal contact problem;

• the LCP (4.52) for the tangential one.

Contrary to the rigid-deformable contact of section 1.3.3, the Delassus operators H N and H T are non-diagonal for a deformable-deformable contact. The mortar projection between the two non-conforming meshes connects the contact DOFs. But thanks to CD-Lagrange time-integrator, the directions of the contact problems are known before solving the contact. As a consequence, the contact problems are linear and require only a solver for linear complementary problems.

A first major result is that the LCPs have an unique solution. Indeed, theirs matrices H N and H T are symmetric positive-definite. This places them in P-matrix family, and a LCP defined by such a matrix has a unique solution [see 4].

For the positive-definite feature (here for H N but similar for H T ):

∀x ∈ R n \ 0, x t H N x = x t [B N,1 -B N,12 ]M -1 [B N,1 -B N,12 ] t x = [B N,1 -B N,12 ] t x t M -1 [B N,1 -B N,12 ] t x x t H N x > 0 by positive-definite property of M -1
And for the symmetry: LCP(H, v free ) : Hr = vv free and ∀i, 0 (v) i ⊥ (r) i 0 (4.62)

H t N = [B N,1 -B N,12 ]M -1 [B N,1 -B N,12 ] t t ⇔ = [B N,1 -B N,12 ]M -1 [B N,1 -B N,12 ] t by sym. of M -1 ⇔ H t N = H
The solving is based on a splitting-method [START_REF] Acary | Numerical Methods for Nonsmooth Dynamical Systems 1st ed[END_REF], which decomposes the matrix H in two parts H L and H U . The splitting-method relies on the property:

If r is solution of LCP(H, v free ), r is solution of LCP(H L , v free + H U r)
It leads to an iterative process which gives at convergence the solution of the LCP (4.62):

while r k+1 -r k > tol do r k+1 ← LCP(H L , v free + H U r k ) end while At convergence r k+1 = r k is the solution of LCP(H, v free ). For the splittingmethods, each iterate r k is a closer approximation of the solution. This property brings robustness: the algorithm can be stopped everywhere and gives nevertheless an approximated solution. Several splitting methods exist depending on the split of the matrix [see 4]. They correspond to the standard linear solvers: Jacobi, Gauss-Seidel, overrelaxation methods...

Here a projected Gauss-Seidel [START_REF] Acary | Numerical Methods for Nonsmooth Dynamical Systems 1st ed[END_REF] is used. The matrix H is decomposed into a strictly lower and upper parts:

H L = {H} ij , j < i H U = {H} ij , j > i (4.63)
r k+1 , the solution of LCP(H L , v free + H U r k ), is obtained by (N being the size of r):

for 0 i < N do (r k+1 ) i ← max 0, 1 H ii -v free -H U r k -H L r k+1 i
end for or with a summation form:

for 0 i < N do (r k+1 ) i ← max 0, 1 H ii -(v free ) i -N -1 j=i+1 H ij r j -i-1 j=0 H ij r j
end for This algorithm stores the vector while it updates it: only one vector is necessary in memory, but this algorithm is sequential.

The two contact LCPs, eqs. (4.51) and (4.52), differs in their resolutions because of the different complementary relations. The normal one is solved by algorithm 7 and the tangential by algorithm 8. The tangential solving requires the normal impulse, and must then follow the solving of the normal LCP.

These algorithms require an initial value r 0 . If the algorithm converges with a null r 0 , a clever initialisation is to use the preceding values of r n+ 1 2 for the nodes already in contact. Indeed as the time-step is short, the impulses do not vary a lot between two consecutive times (except at impact). A precise initial value accelerates the convergence.

, v free = v free N,n+ 3 2 , g = g(U n+1 ), H = H N 1: r k+1 ← r 0 2: while r k+1 -r k > tol do 3: r k ← r k+1 4: for 0 i < N do 5: if (g) i 0 then 6: (r k+1 ) i ← 0 7: else 8: (r k+1 ) i ← -1 H ii (v free ) i + N -1 j=i+1 H ij r j + i-1 j=0 H ij r j 9:
(r k+1 ) i ← max 0, (r k+1 ) i 10:

end if 11:

end for 12: end while 13: r ← r k+1 Outputs : r = r N,n+ 3 2 Algorithm 8 LCP solver for tangential problem Inputs : r 0 , v free = v free

T,n+ 3 2 , r N = r N,n+ 3 2 , H = H T 1: r k+1 ← r 0 2: while r k+1 -r k > tol do 3: r k ← r k+1 4: for 0 i < N do 5: if (r N ) i 0 then 6: (r k+1 ) i ← 0 7: else 8: (r k+1 ) i ← -1 H ii (v free ) i + N -1 j=i+1 H ij r j + i-1 j=0 H ij r j 9: (r k+1 ) i ← min µ(r N ) i , (r k+1 ) i 10: end if 11:
end for 12: end while 13: r ← r k+1 Outputs : r = r T,n+ 3 three main reasons: each iterate is an approximate of the final solution, the convergence is fast for predominant diagonal matrices, and can be accelerated by a precise initial guess. They are then more robust for the CD-Lagrange time-integrator.

The choice of a sequential algorithm in the splitting-based methods, as the projected Gauss-Seidel, may seem odd. Indeed the FE softwares use parallel architecture for accelerating the computations. But if the contact problem is small, a sequential solving is not incompatible with a parallel algorithm. Indeed the contact solving can be included in one process of the global algorithm without any loss of accuracy. Moreover the sequential Gauss-Seidel method has a faster convergence rate than the parallel Jacobi methods for symmetric matrices with a predominant diagonal.

The Gauss-Seidel like algorithm seems then more suitable to shock simulations on tire. Indeed the deformable-deformable contact involves few nodes compared to the rigid-deformable one, or to the global problem.

An acceleration technique: lumping mortar operators

In order to accelerate the solving of LCPs, the operators B N,1 and B T,1 can be lumped:

[(B N,1 ) ij ] i∈[1..N 1 ] j∈[1..N 1 ] =      j γ C φ (1)
i φ

(1)

j dγ C n t i if i = j 0 if i = j (4.64) [(B T,1 ) ij ] i∈[1..N 1 ] j∈[1..N 1 ] =      j γ C φ (1)
i φ

(1)

j dγ C t t i if i = j 0 if i = j (4.65)
This lumping technique of mortar operators is proposed in [START_REF] Bussetta | Efficient 3D data transfer operators based on numerical integration[END_REF]. According to Bussetta et al. the error induced by this approximation is acceptable; and it preserves the positivity of nodal values contrary to standard mortar operators. Indeed, the lumping of mortar operators shapes the matrix like the ones get with dual Lagrange multipliers [START_REF] Wohlmuth | A Mortar Finite Element Method Using Dual Spaces for the Lagrange Multiplier[END_REF].

In the following, this technique is validated on a numerical test case. The numerical results come from a code prototype with a major limitation: the intersection between the contact meshes and the assembly of mortar operators are done once in the initial configuration. If this limits theoretically the numerical tests to small deformations, large deformations are yet considered assuming small sliding at contact. The numerical validation is done in two parts:

• with the space-and-time convergence of non-lumped mortar formulation toward a solution;

• with the space-and-time convergence of lumped mortar to non-lumped mortar solution.

The acceleration induced by the lumping technique will be demonstrated too with these numerical results. The test case is described in fig. 4.2. Two bodies impact each other: the dome which is a rectangular cuboid topped by half of a cylinder, and a cube. The upper face of cube is motionless, and the contact involves no friction. The slave side for contact is on the lower side of cube. As the intersection mesh is computed once in initial configuration, the contact normal is y during the whole simulation.

The material is linear and elastic, and set by the values in linear material close from a "compressible" rubber. The dome and the cube are in the same material. A Dirichlet's condition is applied on the upper face of the cube: its velocity is set to zero. At t 0 the dome and the cube are distant of a gap g 0 along y . The initial velocity of the dome v 0 is only along y in order to reduce sliding along x and z directions.

In order to check the numerical convergence, the problem 4.2 is meshed with five characteristic lengths, from a coarse mesh to a refine one. The table 4.2 gives the five lengths with the critical time-step h CFL and the used time-step h. The computation of h CFL is exact thanks to the formula of eq. (1.98), and h is set to h = 0.7 × h CFL .

The fig. 4.3 defines a central zone on the slave contact side on cube. Coloured in green, it forms a square with sides of 0.33 m, located at the center of the face. For each size of mesh, the boundaries of this central zone coincides with edges of elements. It is then possible to average the nodal quantities on it. As the area stays constant, the average values are comparable between meshes. For the nodal impulse, the average integrates coefficients which depend on the position of the node: 1/4 for a node at a corner of the area, 1/2 at an edge, and 1 inside. For other quantities as displacement, the average is standard. The magnitude of the averaged nodal impulse on the cube is denoted r cube .

In order to evaluate the error, the following Hausdorff measure is used

i -f (t i )| (4.66)
where:

• n is the number of discrete times in interval [t 0 , t f ];

• f i is the discrete solution at time t i ;

• f (t i ) is the analytical or reference solution at time t i .

The norm (4.66) evaluates the distance between the point graph from the discrete solution and the continuous graph from the reference. On the problem considered here, the reference solution is the discrete one with the finest mesh and f i and f (t i ) are averaged on central zone.

Numerical validation for non-lumped mortar

This section presents the numerical results obtained with the prototype code and the non-lumped mortar operators. The test case is described in section 4.2.1: the initial velocity of dome has no tangential component in order to limit the sliding, and no friction is considered. The fig. 4.4 gathers the deformation for meshes 1, 3 and 5. Both solids are coloured according to the magnitude of the normal contact impulse. The deformed meshes are represented at four discrete times: The discrete times for the maximum impulse and release are close between the meshes in regard of the different time-steps. But the impact happens clearly later for the mesh 1. This is due to the precision of the mortar gap which depends strongly on the precision of the meshes. The mortar operators is a kind of an average of the geometrical gap on the elements. Then finer is the mesh, closer are the mortar and the geometrical gap. This is clearly visible at impact time: the geometrical penetration is largely greater for mesh 1 on fig. 4.4a than for mesh 3 on fig. 4.4b. The global deformation validates the assumption of small sliding at contact. The fig. 4.5 presents the time-evolution of r cube (an average on the central zone divided by time-step) for meshes 1, 3 and 5. As noticed on the deformations, the space-time refinement changes the times of impact. If the shape is quite different for mesh 1, r cube seems to converge rapidly. Indeed the r cube for mesh 3 and mesh 5 are similar: impacts and releases are close, and reach the same maximal value. In order to estimate the rate of convergence, the norm (4.66) is used to compute the distance between r cube on a mesh and the reference solution on mesh 5. The fig. 4.6a shows a rate of convergence of O(h 2 ). On fig. 4.6b, the convergence rate is found to be the same but for the displacement. The global displacement U is projected along the contact normal, and averaged on the central zone. The error is then estimated thanks to the preceding norm (4.66) as the distance to the solution on mesh 5.

The figure 4.7 describes the energy balance for the solution on mesh 1. The expected behaviour is observed on contact work: the global contact work ∆W C,N of eq. (4.57) is non-null only at impact; and during contact, the contact works on cube and dome are opposite. The only energy loss happens then at impact, as visible on the system energy curve. On fig. 4.8, the system energies are represented for the five meshes. The energy loss at impact decreases with the mesh refinement. If only one impact is visible for mesh 1, multiple impacts occur for the other meshes. Indeed for mesh 1, due to the coarse mesh, all contacting nodes impact at the first impact time. But for other meshes the contacting nodes settle successively, which leads to multiple impact times. This is visible on fig. 4.4: the nodes in contact at the maximum impulse time are not all already in contact at impact. On fig. 4.7 the initial system energy differs from a mesh to an other. As the approximation of geometry becomes better with a finer mesh, the meshed volume increases slightly which results in a higher mass and kinetic energy.

Conclusion.

These results are a first numerical validation for the CD-Lagrange with mortar methods. Indeed the discrete solutions are accurate and stable:

• despite the complementary relation hold on mortar quantities, the nodal impulses on fig.

Numerical validation for lumped mortar

In this section, the lumping technique presented in section 4.1.4 is applied on mortar operators. The analysis of numerical results demonstrates that it does not influence a lot the discrete solution, while accelerating the LCP solving.

Lumped

Lumped and r 0 = r n+ The fig. 4.11 presents the time-evolution of r cube for both lumped and non-lumped mortar operators (weighted value of r cube on the central zone). The discrete solution is slightly influenced by the lumping: the release happens few time-steps later, and the maximal impulse is higher due to the narrower area. But r cube with lumped mortar operators converges to r cube with standard operators. The fig. 4.12 measures the distance between the two graphs (lumped and non-lumped) for the five meshes thanks to the norm (4.66). It shows a convergence rate of O(h).

The table 4.3 gathers the number of iterations required to solve the normal LCP, under a tolerance of 1 × 10 -8 . It shows two way of acceleration:

Conclusion.

The lumping technique is a powerful lever for accelerating the LCP solving:

• the discrete solution is only slightly modified with the nodal impulses condensed on a narrower area, and it converges to the solution with non-lumped mortar operators;

• the projected Gauss-Seidel algorithm converges with ∼90% less iterations.

Results for large sliding and friction

Test case presentation

In this section, the CD-Lagrange together with mortar operators is tested for large sliding and frictional contact. Following the results of section 4. 4.4: Values for numerical cases -Large sliding faces are free. For the initial conditions, the displacement is null everywhere and only the dome has a uniform velocity v 0 . This initial velocity sets up the global motion: a high tangential velocity results in large sliding. The friction is ruled by a Coulomb's law with a friction coefficient µ. This time, the implementation is performed inside the MEF++ software. It notably integrates a library devoted to contact between deformable bodies in large deformations. Thanks to it, the mortar operators are computed at each time-step on deformed configuration.

In order to check the space and time convergence, the computations are performed on the five preceding sizes of mesh. The table 4.6 gathers the characteristic lengths l e , the stable time-steps h CFL and the effective timestep h for the computation. This time, h is chosen at 0.9×f CFL and rounded to 1 × 10 -5 .

By setting v 0 and µ, three numerical cases are defined in table 4.5: large sliding with non-frictional contact, a frictional contact with a small initial tangential velocity, and a frictional contact with a high initial tangential velocity. The small tangential velocity leads to a sticking contact for almost all the nodes, and the high tangential velocity to a sliding contact.

No-friction and large sliding

Friction and sticking

Friction and sliding v 0 (m s -1 ) 5y 0.5x + 5y 5x + 5y µ 0 0.2 0.2

Shape of discrete solutions

No-friction and large sliding

The initial velocity of the dome v 0 is 5 m s -1 along x and y . It generates a large sliding along x for the dome. The global deformation on mesh 3 is represented on fig. 4.15 at four times. The bodies are coloured according to the magnitude of the normal impulse. Being a nodal quantity, this impulse is depicted also by an arrow at each node on contact sides. Each arrows are in the direction of the impulse, and proportional to its magnitude. The green color indicates the cube contact side, and orange the dome. The normal remains close from y . Indeed the slave side is the cube and its deformation remains limited. The energy balance is depicted on fig. 4.16. The system energy decreases surprisingly during the sliding due to the contact work. Indeed as no friction is considered, only the impacts dissipates energy. The fig. 4.17 explains this energy loss. It represents the global deformation and the contact side of cube at two consecutive times with a large energy loss (visible on fig. 4.16 around 3 × 10 -2 s). The contact impulse is denoted by the color shading, and by the arrows. These discrete times are far from the first impact, and close to the time with the maximal impulse. All nodes inside the contact area are then expected to stay in contact. However on figs. 4.17a and 4.17c, a line of nodes leaves the contact in the middle of the contact area. The mortar gap for them is then above the numerical tolerance for detecting an active contact. Indeed the mortar gap differs from the geometrical one, and depends on the mesh precision. These nodes come again in contact just When the discretization becomes finer, these impacts still happen as shown on fig. 4.18. It depicts the system energy for the five meshes. This energy tends to stabilize around 90% of its initial value when the mesh becomes more precise. If impacts are still visible for mesh 4 and 5, the energy jumps decrease in intensity. This is explained by the contact loss which no longer appears inside the contact area but on its edges. This behaviour is expected: with the sliding motion the contact area on cube moves, and the pressure and the penetration are smaller on the edges. The external nodes are then likely to oscillate between contact an non-contact state. The figure 4. [START_REF] Chung | A Time Integration Algorithm for Structural Dynamics With Improved Numerical Dissipation: The Generalizedα Method[END_REF] shows the contact area on the cube for the mesh 5. At the back of contact area, impacts are visible between the two adjacent times but the central nodes stay in contact.

This test case is particularly challenging on this aspect: the tangential velocity are high and not decreased by friction. It shows nevertheless that energy conservation at impact is a key feature. On a discrete system, impacts happen not only when two bodies come in contact but at every nodal shifts between contact and non-contact state.

Despite these impacts, the global resultant seems to converge. It is evaluated as the sum of all nodal resultant on cube contact side. The fig. 4.20 represents the resultants for mesh 1, 3 and 5 along y , as the main direction for contact normals. For mesh 1, the resultant is quite far from the reference resultant (mesh 5): impact time happens later, and the maximum value is higher. But the resultant for mesh 3 is closer from the reference one. The shape of impulse shows no instability despite the challenging conditions of this case.

Friction with a sticking contact

The initial velocity is set here to 0.5 m s -1 along x and 5 m s -1 along y . Together with a friction coefficient µ = 0. The fig. 4.24 shows the contact side of the dome with the tangential impulses for two consecutive times. The change in direction for tangential impulses are clearly visible. This is caused by a small kinetic momentum in tangential direction. At this time the dome is indeed nearly stopped along x due to friction. The part of the free-velocity along x is close to zero and then its direction is not numerically well determined. If the direction of tangential impulse is difficult to determine here, it is not linked with the CD-Lagrange. Indeed, in the Coulomb's law, the direction is determined thanks to the dynamic as the tangential direction of kinetic momentum. A small kinetic momentum leads to a numerical difficulty in determining the direction of friction.

On fig. 4.23 the oscillations on tangential impulses do not disappear with the space-time refinement: both the frequency and amplitude seem close between mesh 3 and mesh 5. This tends to validate a cause in the discrete Coulomb's law, and not an error due to spatial or time discretization. The normal impulse remains smooth and stable.

Friction with a sliding contact

In order to get a sliding frictional contact, the initial velocity is set to 5 m s -1 along x and 5 m s to that with large sliding. This confirms that the energy loss due to normal impulse comes from the sliding movement, which causes the impacts. On this case there are no oscillations on the direction of tangential impulse. Indeed the free-velocity is well defined as the relative tangential velocity remains high during all the movement.

Numerical performances

Accuracy

The fig. 4.28 represents the error with respect to the reference solution on mesh 5 for the normal resultant along y (sum of nodal impulses) on cube contact side. The distance between graphs is evaluated thanks to the norm (4.66). Whatever the sliding or the friction coefficient, the convergence rate is in O(h 2 ).

On fig. 4.29, the same error is depicted for the tangential impulses. The magnitude of tangential impulse in (x , z ) plane is used for the sticking case, as there is no main direction at the end of the movement. And for the sliding cases only the resultant along x is considered, because x is the sliding direction. The convergence rate differs between the two cases: O(h 1 2 ) for sticking, and O(h 2 ) for sliding. The lower convergence rate for sticking case is not clearly explained. Maybe it is due to the strong oscillations in the direction of tangential impulses on some nodes, which pollutes the global estimation. The CD-Lagrange stays nevertheless robust and accurate on these challenging cases: it converges rapidly to a smooth solution.

Efficiency

The efficiency can not be precisely measured in MEF++. Indeed the implementation lacks of optimization for the CD-Lagrange scheme. Some methods are not yet parallelized, and increase the computational time more than the LCP solver. It was then not possible to determine if the sequential LCP solver is limiting for a parallel run.

The efficiency indicator is then again the acceleration providing by the lumped mortar operators. As in the preceding section, it is evaluated through the number of iterations of projected Gauss-Seidel algorithm. Here no guess for r 0 is used, the solver starts with null values everywhere. The table 4.7 gathers the average number of iterations required for solving the LCP of normal contact problem. The preceding gain of ∼90% is observed between lumped and non-lumped mortar operators. For the tangential LCP, the table 4.8 shows the same gain for a sticking contact. For a sliding one, the convergence is very fast for both lumped and non-lumped operators. Indeed the majority of nodes has a tangential impulse equal to the normal one multiplied by the friction coefficient. The algorithm reaches this value in a couple of iterations.

Conclusion

The three cases of this section show both the abilities and the limitations of the CD-Lagrange scheme for deformable-deformable contact with large sliding and friction. The CD-Lagrange scheme is still:

• robust and accurate, the global deformations and the contact quantities are numerically stable, the discrete solutions converge rapidly, and the energy balance is accurate;

• efficient, the prediction of displacement and free-velocity makes the contact problems only linear, and lumping the mortar operators reduces drastically the computation time.

But its major limitation is spotlighted in the challenging case of large sliding.

The energy loss at impact together with the mortar measure decrease the system energy for large sliding motions.

Conclusion: effective mortar methods for tire simulations

The mortar methods are an effective way to handle a deformable-deformable contact for the CD-Lagrange. They do not decrease the contact stability or accuracy compared to a rigid-deformable contact. Their robustness allows to deal with large sliding motions at contact, with or without friction.

The main drawback is that the deformable-deformable contact problems are no longer explicit. As a consequence, the CD-Lagrange loses its fully explicit feature and its ability to be naturally parallelized. Nevertheless the contact solving stays linear thanks to the explicit feature of the CD-Lagrange. Only a linear solver is required, which makes the CD-Lagrange efficient in contact solving. Moreover a lumping technique is applied on the mortar operators. It accelerates the solving without any loss of accuracy.

The choice of the linear solver is lead by the tire simulations framework. Indeed the deformable-deformable contact involves a small number of nodes compared to the rigid-deformable contact, or the global mesh. The contact solving is then not split between several processes in a parallel run, which makes it compatible a sequential solver.

If the CD-Lagrange is ready for shock simulations on tires, no real cases are addressed in this work. A measure of efficiency is required on such a case to verify the assumption on efficiency of contact solving. A direct improvement for the CD-Lagrange is the integration of a parallel solver for the contact problem. Larger contact problems will be then accessible.

Conclusion and outlook

The performances of the CD-Lagrange scheme A suitable explicit integrator for non-smooth dynamics Thanks to an explicit and symplectic time-integrator, and discrete contact conditions based on Moreau-Jean formalism, the CD-Lagrange scheme presents all the properties required in non-smooth contact dynamics.

The central difference method ensures an explicit time-integration for both displacement and velocity. The displacement is updated once in the time-step. The integration of internal stresses is then done with this displacement and the last discrete velocity. As a consequence, the dynamics is linear on velocity even for non-linear materials. For making explicit the dynamics, a lumping technique is applied on the mass matrix. No system solving is then required to get the updated displacement and velocity. This high efficiency in a time-step leads to a low time-cost for the whole simulation despite the stability condition. Moreover this conditional stability creates a space-time relation between the time-step and the mesh precision. This ensures the accuracy as all phenomenons modelled by the mesh are captured by the scheme.

Even in case of contact, the displacement is not corrected in the timestep. Indeed the contact is enforced on velocity. This facilitates the contact solving as the deformed configuration is already known. The explicit feature is conserved in case of contact between a rigid and a deformable body. In this case, the algorithm of the CD-Lagrange scheme is naturally parallel.

Besides its parallel feature, the efficiency of the CD-Lagrange comes also from its time-stepping aspect. Build in the non-smooth contact dynamics framework, the non-smooth events are integrated in time in a weak sense. The time-step remains constant whatever the number of non-smooth events.

The CD-Lagrange inherits of the symplectic properties of the CD method: conservation of the linear momentum, the angular momentum and the energy. For elastic impacts the angular momentum is conserved to its exact value, and the energy balance is conservative. These properties make the scheme accurate in large rotations and long-time simulations.

However the energy conservation is lost with a deformable contact. In-deed at impact, the kinetic energy of contacting nodes is cancelled. The energy balance stays yet accurate as only the impact is not conservative. Moreover this loss of kinetic energy comes directly from the discrete persistency condition, which ensures a high stability for contact quantities.

The CD-Lagrange is then suitable to non-smooth and non-linear dynamics. But it still present two possible improvements:

• the energy conservation through the non-smooth events would made the scheme fully symplectic even for a deformable contact;

• it must be enlarged to deformable-deformable contact.

The mortar methods: an effective deformable-deformable contact formulation

The deformable-deformable contact is addressed with the mortar methods. They link two non-conforming meshes with a robust projection. Even with large sliding at contact, the mortar formulation decreases only slightly the stability and the accuracy compared to the rigid-deformable formulation.

And the contact can involve friction without any extra difficulty.

If the explicit feature is lost for a deformable-deformable contact, the scheme makes the contact problem only linear. As the size of the deformabledeformable problem is small, the solver is chosen sequential. It presents indeed an high robustness, combined with a fast convergence rate. In order to accelerate the solving, a lumping technique is applied on mortar operators. It shows a clear acceleration without any loss of precision. With a small deformable-deformable contact problem, it should provide enough acceleration to not impact the global efficiency even for parallel simulations.

With the mortar methods, the CD-Lagrange is able to simulate shocks on tire with an high efficiency and accuracy. The CD method handles easily the material non-linearities, and ensures an high accuracy in time-integration. The discrete persistency conditions ensure that the contact quantities are stable and precise. Moreover the contact formulation is parallel for the larger contact zones: the contacts between the tire, the ground and the rim are indeed rigid-deformable. If a deformable-deformable contact happens, it stays localized in a smaller area inside the tire. The lumped mortar formulation should then keep the computational cost low.

The singular mass method: a step toward full symplecticity

The adaptation of singular mass formulation in the CD-Lagrange scheme presents a major difficulty. As the inertial term is cancelled, the normal velocities for contacting nodes are no more determined by the dynamics. A contact law substitutes for it. It determines the normal contact velocities in order to enforce a persistency like contact condition.

On the 1D formulation, the singular mass makes the impact conservative. But at release, the energy conservation is conditioned to a release on a discrete time. In this case, the persistency condition is ensured for all discrete times. The energy conservation is a new result for non-smooth time-integrators. Indeed the energy conservation is reached on the discrete system while keeping the contact stable.

The first extension to 3D meshes involves only extra numerical stiffness added at the contact nodes. The implementation is fully compatible with large deformation and non-linear materials. The numerical rigidities allow to set the skin response and improve the energy balance. But the stability at contact is decreased. This extension is yet a promising way for singular mass in explicit schemes.

The second extension to 3D meshes follows the existing singular mass formulations for implicit schemes. The singular modification is done on the initial mesh. The formulation is quite limited: it requires small perturbations and specific properties for the mesh of the skin. The results demonstrates the limits of the contact law. After release, the contact continues to work and degrades widely the energy balance.

Outlook

The CD-Lagrange seems promising for shock simulations on tires. Nevertheless no test on industrial cases have been made during this thesis. Some levers still exist to improve the scheme for industrial applications. The first one is to adapt the scheme to a variable time-step. With non-linearities the time-step is adjusted during the simulation to guaranty the stability. A variable time-step formulation will introduce no error at a time-step switch. The second one is to develop a parallel solver for the deformable-deformable contact problem.

An other direct continuation to this work would be to bind the singular mass formulation and the mortar methods. Indeed in case of large sliding, the numerical results show that the dissipative impact degrades the energy balance. The singular 3D formulation with normal massless elements could improve this.

For the singular mass formulations, maybe it is possible to keep the positivity of contact constraints at release, and then makes the scheme fully symplectic. Indeed it remains to build a time-stepping symplectic integrator for non-smooth dynamics with deformable contact.

Figure 1 . 1 :

 11 Figure 1.1: Initial and deformed configuration of a deformable body

  the total time derivative which defines the velocity u = du dt and the acceleration ü = d 2 u dt 2 .

Figure 1 . 3 :

 13 Figure 1.3: Types of contact

( 1 .

 1 [START_REF] Schindler | Timestepping schemes for nonsmooth dynamics based on discontinuous Galerkin methods: Definition and outlook[END_REF] to the dynamics at a non-smooth time. The eq. (1.72) states a condition on energy (detailed below). The eqs. (1.70) and (1.73) are the Signori's conditions for smooth time intervals, and non-smooth times.By replacing the Lagrangian with its expression of the eq. (1.62) and by taking U ∈ C, the eqs. (1.69) to (1.73) lead to:

( 1 .

 1 42). Simo et al. set this definition by reproducing the continuous one derived from the Hamiltonian [see 73]. For the Newmark's time-integrator, γ = 1/2 and β = 0 is the only set of parameters which corresponds to det(A) = 1.

  1.4).

2 Figure 1 . 5 :

 215 Figure 1.5: A discrete time interval

1 2 . 2 for

 22 The discrete dynamics (1.123) and (1.147), and the timeintegrations relations in eqs. (1.122) and (1.148) are similar. The contact conditions of eqs. (1.118) and (1.149) have the same shape: a complementary relation between the normal velocity and the normal impulse. The difference lies in their time of expression: t n+1 for Moreau-Jean scheme, and t n+ 1 Laursen-Chawla scheme. As v N,n+ 1 2 is a mean value between Un+1 and Un for the Laursen-Chawla scheme, it leads to an oscillating value of U at integer time t n [see 28]. But the energy balance is conservative. Multiplying the dynamics (1.147) at left by Ut n+ 1 2

2 ( 1 1 2or r n+ 3 2 ; 2 . a dissipative impact with 0 e c < 1 where ∆W c 0, as v n+ 1 2 0 and r n+ 3 2 > 0 2 = 0

 2112222020 .168) by v = L N U and the Newton's law.In the following, the contact work of eq. (1.168) is considered locally for each node (or line) i ∈ γ C , but the indices i are dropped for clarity. Two cases have to be distinguished for the value of ∆W C when the contact involves two rigid bodies:1. an elastic impact with e c = 1, where ∆W C = 0 whatever the values of v n+ at impact, and r n+ 3 else.For two rigid bodies, ∆W C,N is then dissipative if e c < 1 and conservative if e c = 1 as expected.

2. during contact, v n+ 1 2 = 0 2 = 0 and r n+ 3 2 = 0 ,

 20220 as the contact was active at the preceding time-step, and then ∆W C = 0; 3. at release, both v n+ 1 and then ∆W C,N = 0.

2 ( 1 . 179 ) 3 2

 211793 In eq. (1.179) the directions of the tangential impulses I T,n+

T,n+ 3 2 = 2 with α > 0

 220 -α ui,T,n+3 which gives in eq. (1.180):

3 2 .to t determined by Ufree n+ 3 2 ,

 22 Thanks the tangential contact problem of eqs. (1.173) to (1.175) is then a LCP: the direction and the sign of (r T,n+ 3 2

≤ 1 Figure 2 . 1 :

 121 Figure 2.1: Bouncing ball

Figure 2 . 2 :

 22 Figure 2.2: Bouncing ball -Position for e c = 1

Figure 2 . 3 :Figure 2 . 4 :Figure 2 . 5 :

 232425 Figure 2.3: Bouncing ball -Position for e c = 0.75

Figure 2 . 6 :Figure 2 . 7 :Figure 2 . 8 :

 262728 Figure 2.6: Van der Pol oscillator -Position

Figure 2 . 9 :

 29 Figure 2.9: Van der Pol oscillator -Convergence on position with no impact

Figure 2 .

 2 Figure 2.10: Rotating spring

Figure 2 . 11 :

 211 Figure 2.11: Rotating spring -Position over 10 s

Figure 2 . 12 :

 212 Figure 2.12: Rotating spring -Angular momentum

Figure 2 . 13 :

 213 Figure 2.13: Rotating spring -Energy Balance

Figure 2 . 15 :

 215 Figure 2.15: Impacting bar -Moreau-Jean and Laursen-Chawla

Figure 2 . 16 :

 216 Figure 2.16: Impacting bar -Paoli-Schatzman and CD-Lagrange

Figure 2 . 17 :

 217 Figure 2.17: Impacting bar -Convergence rate

Figure 2 . 18 :

 218 Figure 2.18: Impacting dome -Initial configuration

Figure 2 . 19 :

 219 Figure 2.19: Impacting dome -Deformation

Figure 2 . 20 :

 220 Figure 2.20: Impacting dome -Impulses and contact displacement

Figure 2 . 21 :

 221 Figure 2.21: Impacting dome -Energy balance

Figure 2 . 22 :

 222 Figure 2.22: Impacting dome -Energy balance

Figure 3 . 1 :

 31 Figure 3.1: Continuous and discrete in space impacting bar problem

˙ u n+ 3 2 2 . 2 and ˙ u n+ 3 2

 3222 where (•) + returns the positive part. The main idea in eqs. (3.10) and(3.11) is to relate the velocities ˙ u of Γ c and the velocity u of Γ b . Two main situations are described by this law:• eq. (3.10) describes the free-of-contact state;• and eq. (3.11) describes the active contact state, where the contact node velocity depends on the skin force.is determined upon the free velocity of the closest node on the bulk un+3 In free-of-contact state, un+ 3

3

 3

2 and ˙ u n+ 3 2

 22

2 + hM - 1 F(U n+1 ) Update of free velocity 3 : r n+ 3 2 = h kδ u n+1 Contact state solving 4 : ˙ u n+ 3 2 2 + M - 1

 21324321 by eq. (3.10) and eq. (3.11) Lt r n+ 3 2

1. If ˙ u n+ 3 2 0, 3 2> 0 ˙ u n+ 3 2 = 0 .n+ 3 2 > 0 ,

 2302020 the beam tends to penetrate into the rigid frontier. r n+ is computed from the dynamic to stop the contact node by setting It compensates both the internal stress of skin and the inertia of contact node. The skin is in compression.2. If ˙u the beam tends to leave the rigid frontier. r n+ 3 2

Figure 3 . 2 :Figure 3 .

 323 Figure 3.2: Singular impacting bar -Position of contact node

Figure 3 . 4 :Figure 3 . 5 :

 3435 Figure 3.4: Singular impacting bar -Contact impulse

Figure 3 . 6 :

 36 Figure 3.6: Singular impacting bar -Comparison of system energies

Figure 3 . 7 :

 37 Figure 3.7: Singular impacting bar -Comparison of impulses

Figure 3 . 8 :Figure 3 . 9 :

 3839 Figure 3.8: Singular impacting bar -Influence of k

Figure 3 . 10 :

 310 Figure 3.10: 3D formulation -Normal massless elements

  (a) Initial configuration (b) Central contact zone

Figure 3 . 11 :

 311 Figure 3.11: Impacting dome with massless elements

Figure 3 . 12 :

 312 Figure 3.12: Imp. dome with massless el. -Deformed mesh

Figure 3 . 13 :Figure 3 . 14 :

 313314 Figure 3.13: Imp. dome with massless el. -Averaged contact quantities

Figure 3 . 15 :

 315 Figure 3.15: 3D formulation -Singular mass skin

Figure 3 . 16 :

 316 Figure 3.16: Impacting dome with singular mass skin

Figure 3 . 17 :

 317 Figure 3.17: Impacting dome with singular mass skin -Contact quantities

Figure 4 . 1 :

 41 Figure 4.1: Contact between two deformable bodies

  π 21 the map defined by eq. (1.45).

C

  . A other way would be to build an extra discrete interface for expressing the contact.

( 1 )

 1 C to themselves, and B N,12 and B T,12 those of γ

j

  are common to B n,1 and B T,1 , and the sum of φ (1) i φ

  nodal values of the discrete field, gathered in the global vector U. The discrete field representing the relative velocity (according to γ (1) C ) is:

H T,n+1 r T,n+ 3 2 = v T,n+ 3 2 -

 22

2 and v T,n+ 3 2 :

 22

  n+

3 2 2 .

 32 the nodal vector valued free-velocity extracted from Ufree n+ 3 This free-velocity is used in eqs. (4.31) and (4.32) and finally the {t i } i∈[1..N 1 ] at t n+

N

  Though the normal and tangential LCPs defined in eqs. (4.51) and (4.52) are different in the complementary relation, they reduce in the following generic LCP (4.62):

  About the choice of the algorithm. The LCP solving methods are classified in two main categories [see 4]: the splitting-based and the pivotingbased methods. The splitting-based methods are preferred in this work for Algorithm 7 LCP solver for normal problem Inputs : r 0

Figure 4 . 2 :

 42 Figure 4.2: Initial configuration -Small sliding

Figure 4 . 3 :

 43 Figure 4.3: Lower side of cube -Central zone

  |f

1 .

 1 at impact for figs. 4.4a to 4.4c; 2. during contact for figs. 4.4d to 4.4f; 3. at the maximum of nodal impulse for figs. 4.4g to 4.4i; 4. at release for figs. 4.4j to 4.4l. 145 (a) Mesh 1 -Impact (b) Mesh 3 -Impact (c) Mesh 5 -Impact (d) Mesh 1 -Contact (e) Mesh 3 -Contact (f) Mesh 5 -Contact (g) Mesh 1 -Max. imp. (h) Mesh 3 -Max. imp. (i) Mesh 5 -Max. imp. (j) Mesh 1 -Release (k) Mesh 3 -Release (l) Mesh 5 -Release

Figure 4 . 4 :Figure 4 . 5 :

 4445 Figure 4.4: Deformed meshes -Non-lumped mortar

  Error on u cube

Figure 4 . 6 :

 46 Figure 4.6: Convergence rates -Non-lumped mortar

Figure 4 . 7 :Figure 4 . 8 :

 4748 Figure 4.7: Energy balance on mesh 1 -Non-lumped mortar

Figure 4 . 9 :

 49 Figure 4.9: Deformation on mesh 3 -Lumped mortar

Figure 4 . 10 :

 410 Figure 4.10: Normal impulse on lower side of cube -Mesh 3

Figure 4 . 11 :

 411 Figure 4.11: Time evolution for r cube on mesh 3

Figure 4 . 12 :

 412 Figure 4.12: Convergence for r cube

Figure 4 . 13 :

 413 Figure 4.13: Localization of terms in H for mesh 3

Figure 4 . 14 : 1 × 10 7

 41417 Figure 4.14: Initial configuration -Large sliding Density ρ 2000 kg m -3 Young's modulus E 1 × 10 7 Pa Poisson's coefficient ν 0.3 Initial time t 0 0 s Final time t f 0.08 s Initial gap g 0 1 × 10 -4 m Table 4.4: Values for numerical cases -Large sliding

Figure 4 . 15 :

 415 Figure 4.15: Deformation on mesh 3 -No-friction and large sliding

Figure 4 . 16 :

 416 Figure 4.16: Energy balance on mesh 3 -No-friction and large sliding

Figure 4 . 17 :

 417 Figure 4.17: Contact loss -No-friction and large sliding

Figure 4 . 18 :

 418 Figure 4.18: System energy -No-friction and large sliding

Figure 4 . 19 :

 419 Figure 4.19: Cube contact side for mesh 5 -No-friction and large sliding

Figure 4 . 20 :

 420 Figure 4.20: Resultants on cube contact side -No-friction and large sliding

2 ,

 2 this leads to a sticking contact. The global deformations are depicted in fig.4.21. The bodies are coloured according to the magnitude of normal impulse. This time, the arrows represent the tangential impulses: scaled proportionally to the magnitude, and along their directions. The green is still for the cube, and orange for the dome. This global deformation shows that the initial kinetic momentum along x is dissipated by the friction. On fig.4.21b, the tangential impulses are globally along x , the direction of initial tangential velocity. But on fig.4.21c, the tangential impulses show that the dome pushes only along y on the cube. A closer look at tangential impulses during the motion shows that no sliding happens. All nodes are in the sticking state of the Coulomb's law.The energy balance for mesh 3 is represented on fig.4.22. The work of tangential impulses is close to zero. The energy loss is almost entirely caused by the work on normal impulse during impacts. Indeed in this case, the tangential impulses are almost negligible regarding to the normal ones.

  160

Figure 4 . 21 :

 421 Figure 4.21: Deformation of mesh 3 -Friction and sticking contact

Figure 4 . 22 :

 422 Figure 4.22: Energy balance on mesh 3 -Friction with sticking contact

Figure 4 . 23 :

 423 Figure 4.23: Resultants on cube contact side -Friction and sticking contact

Figure 4 . 24 :

 424 Figure 4.24: Oscillations in direction of tangential impulses

  -1 along y . The global deformations on fig. 4.25 show indeed a large sliding for the dome on the cube. The bodies are still coloured according to the norm of normal impulse, and the arrows represent

Figure 4 . 25 :

 425 Figure 4.25: Deformation on mesh 3 -Friction and sliding contact

Figure 4 . 26 :

 426 Figure 4.26: Resultants on cube contact side -Friction and sliding contact

Figure 4 . 27 :

 427 Figure 4.27: Energy balance on mesh 3 -Friction with sliding contact

  FOLIO ADMINISTRATIF THÈSE DE L'UNIVERSITÉ DE LYON OPÉRÉE AU SEIN DE L'INSA LYON NOM : DI STASIO DATE de SOUTENANCE : 27/05/2021 Prénoms : Jean, Joseph

3

  The singular mass method and the CD-Lagrange scheme 3.1 The singular mass method in existing schemes . . . . . . . . . 3.2 A 1D formulation for the CD-Lagrange with singular mass matrix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.1 The Impacting bar with a singular mass . . . . . . . . 3.2.2 Numerical results on impacting bar with a singular mass106 3.2.3 Conclusions on 1D formulation . . . . . . . . . . . . . 3.3 Extension to 3D by normal massless elements . . . . . . . . . 3.3.1 The normal massless elements formulation . . . . . . . 3.3.2 Numerical results . . . . . . . . . . . . . . . . . . . . . 3.3.3 Conclusion on normal massless elements . . . . . . . . 3.4 Extension to 3D by the singular mass skin . . . . . . . . . . . 3.4.1 The singular mass skin formulation . . . . . . . . . . . 3.4.2 Numerical results . . . . . . . . . . . . . . . . . . . . . 3.4.3 Conclusion on singular mass skin formulation . . . . . 3.5 Conclusion: singular mass in the CD-Lagrange . . . . . . . . Conclusion: effective mortar methods for tire simulations . . .

	Notations
	Conclusion and outlook
	Bibliography

4 Deformable-deformable contact 4.1 Discretization of contact between two deformable bodies . . . 4.1.1 Continuous description of contact . . . . . . . . . . . . 4.1.2 Spatial discretization by mortar methods . . . . . . . . 4.1.3 Time-discretization by the CD-Lagrange scheme . . . 4.1.4 Solving the contact LCPs . . . . . . . . . . . . . . . . 4.2 Results for small sliding without friction . . . . . . . . . . . . 4.2.1 Test case presentation . . . . . . . . . . . . . . . . . . 4.2.2 Numerical validation for non-lumped mortar . . . . . . 4.2.3 Numerical validation for lumped mortar . . . . . . . . 4.3 Results for large sliding and friction . . . . . . . . . . . . . . 4.3.1 Test case presentation . . . . . . . . . . . . . . . . . . 4.3.2 Shape of discrete solutions . . . . . . . . . . . . . . . . 4.3.3 Numerical performances . . . . . . . . . . . . . . . . . 4.3.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . 4.4

Table 1 :

 1 Notations for a continuous system

	Ω 0 , Ω(t)	:	domain of the initial and deformed configurations
	Γ D , Γ N	:	part of ∂Ω 0 with Dirichlet's and Neumann's con-
			ditions
	x , X	:	3D coordinates vectors in deformed configuration
			and initial configuration
	u, u, ü	:	displacement, velocity and acceleration fields
	F , C , E	:	the deformation gradient, the right Cauchy-Green
			tensor, and the Green-Lagrange strain tensor
	σ	:	the Cauchy stress tensor
	P, S	:	first and second Piola-Kirchoff tensors
	∇ X (•), ∇ X • •	:	the gradient and divergence operators, ∇ being
			the nabla vector differential operator with the ap-
			propriate coordinates indicated by the subscript
	ρ, f v	:	the mass and external forces densities
	U, V	:	the functional spaces for the solutions to dynamics
			and the test functions

Table 2 :

 2 Notations for a discrete system

	M, K	:	the mass and rigidity matrices
	F ext , F int (U, U) :	the integrated contribution of external and in-
			ternal forces
	U, U	:	global vectors gathering the DOFs for displace-
			ment and velocity
	L(U, U)	:	the discrete Lagrangian
	T ( U), V(U)	:	the kinetic and potential energies
	h	:	the time-step
			6

Table 3 :

 3 Notations for a continuous contact

	Γ (i) C , γ	(i) C	:	the contact boundary for the body (i) in initial
				and deformed configurations
	n		:	the inner normal on slave contact boundary
	g(x )		:	the normal gap
	µ		:	the friction coefficient in Coulomb's law
	λ N		:	the normal contact force
	λ T , λ T , t	:	the tangential contact force, its magnitude and its
				direction
	v N , v T	:	the normal and tangential relative contact velo-
				cities
	u r		:	relative velocity according to the slave surface
	e c		:	restitution coefficient at impact

Table 4 :

 4 Notations for a discrete contact

	g(U)	:	a vector containing the normal gap
	L N , L T	:	the projection operators along the contact normal,
			and the tangential contact direction
	λ N , r N	:	vectors with the forces and impulses for normal
			contact actions
	λ T , r T	:	vectors with the forces and impulses for tangential
			contact actions
	v N , v T	:	vectors with the normal and tangential relative
			nodal velocities
	B N , B T	:	the normal and tangential mortar operators
	B N,1 , B T,1	:	the parts of mortar operators concerning the slave
			or non-mortar surface
	B N,12 , B T,12

Table 5 :

 5 Notations for a discrete system with singular mass

	γ C	:	the discrete contact surface with the massless
			nodes
	γ B	:	the interface between the bulk and the skin
	n	:	the inner normal to γ B
	u, ˙ u	:	vectors with the displacements and velocities for
			the massless DOFs of γ C
	ū, u	:	vectors with the displacements and velocities for
			the normal DOFs of γ B
	L, L	:	the projection operators which select the DOFs of
			γ B and γ C
	k	:

  2. to demonstrate the feasibility of robust, efficient and accurate shock simulations on tire. The first objective is purely academic. It aims to identify the weak points of the current explicit time-integrators for impact dynamics, and propose new solutions to bring them further. But the target maturity will not allow an application in an industrial framework, contrary to the second objective. Its goal is to build a time-integrator with mature techniques in order to demonstrate the feasibility of shock simulations on tire. The time-integrator will be implemented in the finite elements solver used by Michelin. It rely on MEF++, a finite elements code co-developed by the GIREF (Groupe Interdisciplinaire de Recherche en Éléments Finis) from the Laval University in Quebec, and industrial partners as Michelin or Hydro-Quebec. MEF++ is a research tool in numerical methods for both academic and industrial projects. It integrates indeed parallel libraries for system solving, which make it efficient enough to address large problems met in industrial framework.

  n+1 and r n+1 are both unknown but linked by the complementary condition. This set of equations forms a Linear Complementary Problem (LCP), which requires a specific solver [see 4].

  The tangential direction t i is known and detailed in the following. The scalar values for normal impulse {r i,N,n+ 3 2 } i and tangential impulses {r i,T,n+ 3

		n+ 3 2	t i	(1.171)
	are gathered in vectors r N,n+ 3 2	and r T,n+ 3 2		2	} i

Table 2 .

 2

	e x

4: Impacting bar -Numerical values momentum is crucial too. In smooth dynamics, an average conservation of angular momentum is enough to ensure a good long-time behaviour. But here with non-smooth events, a drop of angular momentum is observed which leads to a less accurate discrete solution for long time simulations.

Table 2 . 5

 25 -Schatzman scheme achieves the impenetrability condition

	Mesh char. length	l e	0.9 × 1/6 m
	Time-step	h	9.6 × 10 -4 s
	Density	ρ	2000 kg m -3
	Young modulus	E	1 × 10 7 Pa
	Poisson's ratio	ν	0.3
	Initial gap	g 0 1.02 × 10 -1 m
	Initial velocity	v 0	3x + 5y m s -1
	Coulomb's coefficient µ	0.2

: Impacting dome -Numerical values together with a stable contact. But as explained in section 1.3.3, the contact is not necessary dissipative contrary to CD-Lagrange and Moreau-Jean schemes. This can lead to stability issues.

  .[START_REF] Courant | Uber die partiellen Differenzengleichungen der mathematischen Physik[END_REF] shows a convergence rate in O(h) as on the

	10 -3
	10 -4
		0.5			1	1.5	2
						10 -3
			Mesh 1	Mesh 2	Mesh 3
	l e	(m)	0.9 × 1/3	0.9 × 1/6	0.9 × 1/9
	h CFL (s)	2.07 × 10 -3 1.07 × 10 -3 7.44 × 10 -4
	h	(s)	1.86 × 10 -3 9.6 × 10 -4	6.7 × 10 -4
				Mesh 4	Mesh 5
		l e	(m) 0.9 × 1/12	0.9 × 1/15
		h CFL (s)	5.32 × 10 -4 4.13 × 10 -4
		h	(s)	4.8 × 10 -4	3.7 × 10 -4
	Table 2.6: Impacting dome -Characteristics of meshes
				93

Table 3 .

 3 2: Singular impacting bar -Numerical values the validity of the proposed contact law. Indeed the contact node has the same behaviour as with the consistent contact conditions.

	3: Singular impacting bar -Velocity of contact node
	107

Table 3 . 3

 33

: Imp. dome with massless el. -Numerical values

Table 3 . 4

 34

: Impacting dome with singular mass skin -Numerical values assumed. The contact happens on a rigid xz plane positioned 2.67 × 10 -4

  12 λ T (4.19) The scalars {λ N,i } i∈[1..N 1 ] and {λ T,i } i∈[1..N 1 ] are gathered in global vectors: λ N and λ T . The matrices B N,1 , B N,12 , B T,1 , B T,12 form the mortar operators. B N,1 and B T,1 relate the DOFs of γ

table 4

 4

	.1.

Table 4 .

 4

2: Characteristics of meshes -Small sliding cases

  4.5 are stable as the global deformations on fig. 4.4; • the scheme is dissipative only at impact and this energy loss decreases with space-time convergence, as shown by fig. 4.7 and fig. 4.8; • even if established only numerically, the fig. 4.6a and fig. 4.6b indicate the convergence of the discrete contact quantities.

Table 4 .

 4 3: Average number of iterations in LCP solving fig.4.10 represents the normal impulse on lower side of cube. The elements where r cube is null are coloured in blue, and the elements with positive values in a a red shading. On the non-lumped mortar figure 4.10a the area with a positive r cube is larger than that on fig.4.10b for lumped mortar. Indeed the terms of B lump N,1 are more localized around its diagonal which leads to a narrower area for r cube .

	1
	2

  2, the mortar operators are lumped. The test case is close to that of section 4.2: a dome impacts a cube. But the dimensions, described in fig.4.14, are modified for sliding with a larger cube. The material is the same: elastic and linear, it does not represent a real material but gets close to rubber. The table 4.4 gathers the values used in simulations. The boundary conditions are the same: the upper face of cube is held in its initial position, and other

Table 4 .

 4

		5: Values for initial velocity and µ -Large sliding
			Mesh 1	Mesh 2	Mesh 3
	l e	(m)	1/3	1/6	1/9
	h CFL (s)	2.57 × 10 -3 1.19 × 10 -3 7.82 × 10 -4
	h	(s)	2.32 × 10 -3 1.07 × 10 -3 7.00 × 10 -4
				Mesh 4	Mesh 5
		l e	(m)	1/12		1/15
		h CFL (s)	6.13 × 10 -4 4.91 × 10 -4
		h	(s)	5.50 × 10 -4 4.40 × 10 -4

Table 4 . 6

 46

: Characteristics of meshes -Large sliding

Table 4 .

 4 8: Mean number of iterations for solving tangential LCPs

	10 4
	10 3
	10 2	0.6	0.8	1	1.2	1.4 1.6 1.8 2 2.2
						10 -3
		Figure 4.28: Convergence for normal impulses

  TITRE :The CD-Lagrange scheme, a robust explicit time-integrator for impact dynamics: a new singular mass formulation, and an extension to deformable-deformable contact Les pneumatiques sont complexes à simuler car les matériaux y sont hétérogènes, incompressibles et non-linéaires. De plus la géométrie descend jusqu'à l'échelle millimétrique pour les sculptures de la bande de roulement, ce qui requiert un maillage fin. Le modèle éléments finis présente donc un grand nombre de degrés de liberté, reliés par des équations non-linéaires. En dynamique, la simulation est d'autant plus compliquée avec des chocs. Néanmoins elle est cruciale dans le processus de conception pneumatique, où elle apporte une meilleur compréhension de la physique ceci sans tests réels. Les schémas explicites rendent possible les simulations de chocs, car ils résolvent facilement les non-linéarités avec un coup calcul bas. Associés à une formulation de contact précise , ils forment des schémas robustes, précis et efficaces pour la dynamique nonlinéaire avec impacts. Ce travail vise à choisir et un tel schéma, et l'améliorer pour la simulation de chocs sur pneumatiques. La première partie est un benchmark identifiant le schéma CD-Lagrange. L'intégration temporelle est réalisée par le schéma de la différence centrée, et le contact imposé par multiplicateurs de Lagrange sur la vitesse. Deux possibilités d'amélioration sont identifiées. La première est d'atteindre un impact conservatif, seul instant où le schéma n'est pas symplectique. La seconde amélioration est d'étendre la formulation au contact déformable-déformable. La deuxième partie vise à atteindre la conservation de l'énergie à l'impact en adaptant la méthode de la masse singulière au CD-Lagrange. Une première formulation 1D est construite. Elle démontre une amélioration majeure du bilan d'énergie. Deux formulations 3D sont ensuite explorées. La troisième partie introduit les méthodes mortier dans le CD-Lagrange. Elles permettent de traiter un contact déformabledéformable de manière robuste, même en présence de friction et de grands glissements. Une technique d'accélération est proposée pour résoudre le problème de contact, ceci sans perte de précision. MOTS-CLÉS : intégrateurs temporels explicites ; dynamique du contact non-régulière ; intégration symplectique ; masse singulière ; méthodes mortier ; contact déformable-déformable. Laboratoire (s) de recherche : LaMCoS, INSA Lyon, Université de Lyon. Directeur de thèse: David Dureisseix, Anthony Gravouil Président de jury : Pierre Alart Composition du jury : Ludovic Noels, Françoise Krasucki, Barbara Wohlmuth, Gabriel Georges, David Dureisseix, Anthony Gravouil, Thomas Homolle

	NATURE : Doctorat	Numéro d'ordre : 2021LYSEI029
	Ecole doctorale : MEGA (Mécanique, Energétique, Génie civil, Acoustique)
	Spécialité : Génie Mécanique
	RESUME :

Cette thèse est accessible à l'adresse : http://theses.insa-lyon.fr/publication/2021LYSEI029/these.pdf © [J. Di Stasio], [2021], INSA Lyon, tous droits réservés

To go further than this simple model, more complex model were developed. Some of them are described in[START_REF] Wriggers | Computational Contact Mechanics 2nd[END_REF] with friction coefficients depending on the tangential velocity, the contact pressure, the temperature...

Remerciements Acknowledgements

We gratefully acknowledge the French National Association for Research and Technology (ANRT, CIFRE grant number 2017/1555).

This work was supported by the "Manufacture Française de Pneumatiques Michelin".