Keywords: Blunt impact, Dynamic response, porcine ribs, HR-pQCT, Finite Element simulations, LS-Dyna, Split Hopkinson Pressure Bar

Préambule

Le thorax contient le coeur et les poumons qui sont les éléments principaux des systèmes respiratoire et circulatoire. Les fractures costales sont de fréquentes lésions causées par les accidents de la route ou des impacts balistiques à grande vitesse. Les fragments osseux créés par les impacts thoraciques peuvent blesser le coeur, les poumons, la plèvre, et conduire à des pneumothorax ou des hémothorax. Il est donc intéressant de développer une méthode de prédiction des fractures thoraciques sous des chargements de type impact dynamiques. De nos jours, les modèles biomécaniques humains basés sur la méthode des éléments finis (EF) sont de plus en plus employés pour analyser les mécanismes lésionnels. De tels modèles peuvent être utilisés comme substituts aux expérimentations réelles couteuses en temps et en ressources. En particulier, les modèles EF sont devenus des outils intéressants pour étudier la réponses des organes du thorax en conditions tant statiques que dynamiques et cela jusqu'aux impacts balistiques.

Un certain nombres d'études ont été réalisées pour étudier la réponse du thorax dans les crashs automobiles. Cependant, ces travaux ne permettent pas de connaître précisément la localisation de la fracture et le faciès de rupture, malgré les progrès réalisés. En effet, les modèles thoraciques restent simplifiés, car ils ne prennent pas en compte, d'une part, la progressivité de la fracture (propagation) couplée à l'endommagement du matériau, et d'autre part, la géométrie et l'hétérogénéité de la structure et du matériau. Il subsiste donc toujours un besoin de caractériser expérimentalement le comportement osseux en dynamique rapide pour identifier le mécanisme de la rupture et pour élaborer une loi de comportement réaliste implémentée dans un modèle EF innovant.

Objectifs

Dans cette thèse, nous nous intéresserons uniquement aux impacts non pénétrants.

Afin de minimiser les risques d'apparition de lésions létales, il est nécessaire de :

1. Développer une méthodologie expérimentale spécifique pour caractériser la fracture en dynamique rapide sous impact non pénétrant.

iii 2. Développer un modèle numérique tri-dimensionnel réaliste de la structure.

Verrous de la thèse

Pour atteindre les objectifs fixés pour la thèse, nous avons identifié que : (i) les propriétés mécaniques des constituants de la côte, (ii) la loi de comportement et la loi de l'endommagement, (iii) les effets de la vitesse de déformation, doivent être considérées avec une égale importance pour parvenir a une modélisation pertinente de la réponse de la structure osseuse.

Les verrous à lever sont les suivants :

• Propriétés mécaniques -Corrélation de la masse volumique apparente avec l'échelle des niveaux de gris ;

-Définition des relations entre les propriétés (module d'Young, limite élastique, etc.) et la masse volumique ;

• Formuler la loi de comportement et de l'endommagement couplés.

• Implémenter cette formulation dans un code numérique adapté à la dynamique rapide.

Structure du mémoire de thèse

Le présent mémoire s'organise en cinq parties articulées en sept chapitres.

La partie 1 présente l'état de l'art et l'étude bibliographique.

Dans le chapitre 1, nous présentons d'abord une description anatomique de la cage thoracique humaine, les conséquences cliniques des fractures des côtes et différentes échelles empiriques d'évaluation de la gravité des lésions thoraciques. Puis, certains éléments de biologie osseuse sont présentés.

La côte appartient à la catégorie des os plats. Elle a la forme d'un arc dont la ligne moyenne présente une certaine torsion entre ses deux extrémités. Les côtes les plus fréquemment sujettes aux fractures lors d'impacts rapides sont les troisièmes, quatrièmes, cinquièmes et sixièmes côtes. L'évaluation du la sévérité des lésions par les modèles présentés dans ce chapitre ne peuvent pas être utilisés dans cette thèse, car ils s'appliquent au thorax entier.

Concernant la composition biologique des tissus osseux, il existe neuf niveaux hiérarchiques. Parmi ceux-ci, nous n'en considérerons que deux : l'os cortical et l'os trabéculaire (niveau VIII) et la structure côte entière (niveau IX). Nous allons iv supposer que le "cortical" et le "trabéculaire" sont composés d'un matériau de même type. Ainsi, le comportement mécanique sera dû à la variation de l'architecture entre le cortical et le trabéculaire. Les propriétés de ce matériau (en particulier les propriétés élastiques) sont fondamentalement corrélées à la densité apparente.

De plus, nous introduirons la dépendance de certaines propriétés (par exemple la déformation à la rupture) vis-à-vis de la vitesse de déformation.

 portant sur la réponse quasi-statique en flexion trois points de côtes humaine, la partie trabéculaire n'a pas d'influence sur la capacité de résistance à l'effort (chargement) de la côte. pourtant d'autres études soutiennent le contraire. Nous pensons donc que les modèles existants doivent être améliorés tant dans leurs représentation géométrique de la structure que dans leur formulation du comportement du matériau et de leur prédiction de la fracture. Ce triple objectif conduit au développement de modèles de la structure isolée. De plus, il n'existe quasiment pas de modèles permettant d'aborder le régime des vitesses de déformation ( ε) 1 aux valeurs intermédiaires (150 s -1 ≤ ε ≤ 1200 s -1 ) ou élevées ( ε > 1200 s -1 ). Enfin, il nous paraît nécessaire de tenir compte simultanément de l'effet de l'accumulation de l'endommagement.

Dans la partie II, nous présentons l'exploitation d'un ensemble de résultats expérimentaux obtenus antérieurement au laboratoire. Ces données ont été obtenues lors d'une campagne d'essais de flexion trois points avec le montage des barres de Hopkinson. Tous les essais n'avaient pas été dépouillés et leur analyse n'avait été que partielle. Un début d'analyse avait été fait avec un modèle simplifié de géométrie et de comportement. Nous approfondirons l'analyse au moyen du modèle développé dans notre thèse. Parmi les côtes testées, 27 côtes on été considérées pour l'exploitation. Nous avons noté trois réponses différentes (pas de rupture, rupture pendant l'équilibre entre la force entrante et les forces sortantes, et rupture avant l'équilibre), qui sont liées à la vitesse relative des interfaces. Nous présenterons les résultats expérimentaux sous forme de 3 couloirs ou gabarits. Gabarit 1 : pas de fracture ; gabarit 2 : fracture pendant la phase d'équilibre de l'essai aux barres de Hopkinson ; gabarit 3 : fracture avant la phase d'équilibre (ou "hors équilibre"). A l'intérieur de chaque couloir, on note une variation des réponses due aux variations de la géométrie et des propriétés entre les échantillons réels.

Grâce à ces gabarits, nous constatons que le comportement de la côte apparaît comme élastique-viscoplastique, car la vitesse de déformation n'a pas d'influence 1 également appelé taux de déformation v sur la partie linéaire de la réponse. De plus, le pic de force maximale n'influence pas la rupture. La rupture de la côte a toujours comme point de départ une fissure de la face inférieure (sollicitée en tension). Cette fissure se propage ensuite jusqu'au point d'impact.

La partie III concerne les simulations numériques. Elle se compose des chapitres 3 et 4.

Dans le chapitre 3, nous présentons la méthodologie suivie pour construire la structure numérique. Pour construire le maillage 3D, nous avons utilisé le module ScanIP du logiciel Simpleware pour traiter les images obtenues par imagerie à haute résolution (HR-pQCT). Ensuite, la masse volumique apparente ρ app a été corrélée aux niveaux de gris (exprimée en Houndsfiled Units [HU]). Grâce a la masse volumique apparente, nous pouvons relier chaque niveau de gris à la la porosité Po. Les valeurs de porosité, nous permettent de distinguer le type de tissu osseux, puis de déterminer les propriétés mécaniques par des lois puissances. Les propriétés mécaniques des constituants de la côte de porc obtenues par cette méthodologie sont dans le même ordre de grandeur que celles utilisées pour les côtes humaines dans les études publiées antérieurement.

Le chapitre 4 présente la mise en place du modèle numérique du comportement et de l'endommagement. Nous avons constaté dans le chapitre 2 que la côte possède

I dedicate this thesis to my parents, and to his excellency, late. Dr A.P.J. Abdhul Kalam "Sometimes it is the people who no one imagines anything of who do the things that no one can imagine"

-Alan Turing I am grateful to my parents, siblings, friends and acquaintances who wished me the ultimate success. I consider myself nothing without them. They gave me enough moral support, encouragement and motivation to accomplish the personal goals. Most importantly for their help to learn french. Without their help, I could never imagine that I would speak french as fluent as I do today.

un comportement élastique-viscoplastique. Du point de vue numérique, l'effet de la vitesse de déformation au-delà de la limite élastique est décrit par le modèle de Cowper-Symmonds. Cependant, l'écrousissage peut être soit isotrope soit cinématique, et la loi d'endommagement peut être soit formulée de manière intrinsèque au comportement, soit externe à la loi de comportement. Dans le chapitre 4, nous présentons l'étude de trois formulations. Parmi celles-ci, deux utilisent un endommagement externe dépendant de l'état de contrainte et prenant en compte la vitesse de déformation. Ces deux formulations sont capables de prendre en compte l'évolution non-linéaire de l'endommagement (et donc la diminution de la rigidité de l'échantillon). Ensuite, nous présentons aussi l'algorithme de contact utilisé aux interfaces entre l'échantillon et les barres. De plus, l'échantillon contenant des parties ayant des propriétés mécaniques différentes, nous avons utilisé un autre algorithme de contact, grâce auquel, le solveur est capable de tenir compte de l'hétérogénéité de l'échantillon. Au terme de notre étude, nous avons choisi un comportement élastique-viscoplastique avec écrouissage isotrope pour l'os cortical et un comportement élastique-viscoplastique avec écrouissage cinématique pour l'os trabéculaire.

Le chapitre 5 de la partie IV, présente l'application du modèle numérique prévi senté dans le chapitre 4 aux cas expérimentaux.

Nous nous intéresserons aux cas expérimentaux où la fracture se produit pendant la phase hors équilibre de l'essai aux barres de Hopkinson. Pour mettre en évidence l'influence de l'os trabéculaire sur la réponse dynamique, nous avons testé deux échantillons numériques différents (i) sans os trabéculaire et (ii) avec os trabéculaire. L'échantillon (i) montre une réponse bien plus rigide que l'échantillon (ii) pour la même vitesse donnée à l'impacteur. Cela prouve que l'os trabéculaire joue un rôle très important en participant à l'absorption des forces d'impacts.

La réponse numérique de l'échantillon numérique complet (os cortical et trabéculaire) a été calculée pour deux vitesses d'impacts différentes. Dans tous les cas, le modèle numérique est capable de reproduire la forme et l'amplitude moyenne des cas expérimentaux du gabarit 3. Une des réponses expérimentales est particulièrement proche du calcul numérique.

Enfin, les conclusions et les perspectives de ce travail sont présentées dans la partie V.

Mots Clefs : Impacts dynamiques, Côtes, Fracture, Barres de Hopkinson, Flexion trois points, Éléments-Finis, HR-pQCT, Élastique-viscoplastique, Endommagement, LS-Dyna. french vii viii

Prologue

The thorax contains the primary elements of respiratory and circulatory systems.

Rib fractures are frequent thoracic injuries that are caused by Road Traffic Accidents (RTA's) or High Speed Ballistic Impacts (HSBI's). The bone fragments or fractured ribs originating from chest injuries, may injure the heart or lung plural surface or the lung parenchyma, leading to pneumothorax or hemothorax. Therefore, it is important to develop a modeling method to predict the occurrence of thoracic fracture under impact situations. Nowadays, biomechanical human body Finite Element (FE) models play an increasingly important role in modeling injury mechanisms. Such models can be used as a substitute for time consuming and expensive experimental measurements. In particular, FE models of human thorax are becoming an important tool in the investigation of the structural responses of the thorax and the ribs under various static, dynamic and ballistic loading conditions.

A number of studies were performed to investigate the injury mechanisms of the ribs and the entire rib cage under dynamic loading for automotive crash simulation and injury analyses in the past. Nevertheless, they have failed to properly predict the fracture location and fracture profile. Despite the progress performed in this field, finite element models of the thorax and ribs are simplified and don't take into account the damage, the complete progressive nature of fracture, the precise geometrical representation of the ribs considering the cortical bone and trabecular bone regions and the heterogeneity of the material properties. Therefore, there is still need to characterize the dynamic fracture behaviour of the bone, to identify the fracture mechanisms of bone and to develop bone behaviour law and novel validated finite element model to predict the fracture conditions of ribs and thoraxes under high speed dynamic impacts.

Objectives

In order to reduce the risks of lethal injuries, it is necessary to: ix 1. Develop a specific experimental methodology in order to characterize accurately the bone fracture process (and the possible injury mechanisms) induced by dynamic blunt impacts on the chest at high velocities.

2. Develop a three dimensional numerical model of the bone tissue and structure able to simulate the fracture and to predict its occurrence in specific impact conditions (car accidents, falling of the person, ballistic impact). Such a numerical model will be used to guide the experimental studies and to reduce the costs of the experiments.

Problem Constraints

In order to address the objectives of this thesis, we had identified that the mechanical properties of bone constituents of rib, the parameters influencing the mechanical properties, constitutive and damage law, the effects of strain rate have to be given equal importance in such a way that the realistic behavior of rib is simulated as efficient as possible. The constrains addressed in this thesis are:

• Mechanical properties -Calibration of apparent density ρ app with grey scale intensity CT bone -Respective mathematical models to interpolate the mechanical properties with grey scale intensity CT bone

• Simulation of wave propagation in bars.

• A realistic constitutive law coupled with a damage law.

• Contact algorithm.

Thesis structure

This thesis articulates around five main parts, divided into 7 chapters.

Chapter 1, presents an introduction to the human thorax, the clinical consequences of rib fractures and the existing evaluation methods for chest trauma, followed by the biological and mechanical aspects of rib, existing studies on human ribs and imaging methods used in this field of research.

A rib belongs to the flat bone category of the human skeletal system. It is curved and twisted structure that comprises a head, a body, and a tail. A rib articulates with the thoracic vertebrae with its head and a part of its neck at the x anterior end and with the sternum through costal cartilage at the posterior end.

Most common ribs prone to fracture under dynamic environment are the 3 rd , 4 th , 5 th , and 6 th . The injury severity evaluation models presented in this chapter could not be used in this thesis as they consider the whole thorax for a good prediction.

Regarding biological composition of bone tissues, a bone in general is a hierarchical structure with nine different levels of hierarchy. The concerns of this thesis are restricted to the level VIII (cortical and trabecular bone tissues) and level IX (whole structure) only. Several hypothesis have been put forth regarding the differences in the mechanical behavior of cortical and trabecular bone tissues. We embark with an assumption that these tissues are made of same morphological material i.e. same composition, thus the difference in their corresponding mechanical behavior is due to the difference in architectural arrangement. Therefore, the variables influencing the mechanical properties in this thesis are primarily the apparent density and the mechanical environment i.e. strain rate. Yet the variables regarding bone composition are considered to compare and verify the mechanical properties.

The geometrical constraints posed by bone samples have led to studies on characterization of whole bone structures or segments of bone structures. It has been proved that the anatomical location of the segments, the geometrical variables such as area moment of inertia or the cross sectional area of the geometry influence the mechanical behavior of these structures. The most cited study [START_REF] Granik | Human ribs: Static testing as a promising medical application[END_REF] on quasi-static three point bending of human rib segments concluded that the trabecular bone in the rib does not play an important role in load carrying capacity. Yet, recent studies have proved the contrary. The existing numerical models of human thorax do need developments in geometrical representation, mechanical behavior and fracture prediction. These constraints have led to the development of isolated rib models. Despite the existence of different numerical models, there are almost no numerical models in dynamic loading at intermediate (150 s -1 ≤ ε ≤ 1200 s -1 ) and high strain rate ( ε > 1200 s -1 ). Moreover, there is a need for the development of models that consider the effects on strain rate in both mechanical behavior and damage accumulation as well.

Part II presents the processing of existing experimental data base obtained earlier at our laboratory during the master of [START_REF] Aubert | Étude du comportement des côtes sous sollicitations dynamiques[END_REF] and PhD thesis of [START_REF] Pavier | Contribution à la compréhension des phénomènes physiques lors de l'impact d'un corps sur un modèle de structure biologique[END_REF]. This data was obtained from dynamic three point bending (3PB) experiments using the Split Hopkinson Pressure Bar apparatus for dynamic bending at our lab. Due to a lack of time, not every experimental data were processed, yet an initial analysis of the geometry and its approximation is presented. We present xi a deep analysis of this existing database in this thesis.

Among the porcine ribs tested, only 27 were considered to be worthy for post-processing. The experimental results of these ribs are classified into three different categories of response (no fracture, fracture during equilibrium between input and output forces and fracture before equilibrium), which are influenced by the relative interface velocity. All these results are presented in form of corresponding experimental corridors. Thanks, to this way of presentation, we identify that the rib shows an elasticviscoplastic behavior. Moreover, the rupture of the sample is due to a crack originating from the inferior surface (subjected to tension) and its propagation towards the point of impact. We do notice that the peak of the input force does not have an effect of the type of rupture observed, rather it is influenced by the local displacement and velocity of the sample at the interfaces. The methodology adapted for numerical simulations is presented in Part III in two chapters.

Firstly, the generation of numerical rib sample is presented in Chapter 3. The 3D geometry of porcine rib is obtained High Resolution peripheral Computer Tomography (HR-pQCT) imaging technique. The geometry was reconstructed using the ScanIP module of the Simpleware software for the sake of FE mesh. Once the mesh is done, the apparent density ρ app was correlated to the different greyscale values in (expressed in HU). Assuming that both, cortical and trabecular bones are made of same morphological material, we can calculate the porosity P o values corresponding to each grey scale. The P o values were the deciding factor for choosing power law equations for mechanical properties. The properties obtained from apparent density are in the same order of magnitude to that used be previous authors for human ribs. This thereby proves that the porcine rib a reliable biological surrogate and apparent density is a good measure of mechanical properties.

Secondly, the numerical considerations for the developments of the FE porcine rib model are presented in Chapter 4. This speaks about the constitutive law and damage law required for our rib and other aspects such as mesh size and contact algorithm.

The experimental results have shown that the rib shows elasticviscoplastic behavior. Therefore, we had tested three different constitutive laws, capable of considering the effects of strain rate on the post yield behavior through Cowperxii Symmonds model. Out the three laws tested, two necessitates coupling an external damage law for progressive fracture simulation. The so called external damage law is incremental and is dependent on the stress state in the sample. In addition to this, the damage law is also capable of considering non-linear damage accumulation, reduction in stiffness with increase in damage and the effects of strain rate on the fracture behavior of the sample.

The contact algorithm between sample and bar elements is presented at the end of this chapter. Apart for this, another contact algorithm is defined so that the solver considers the heterogeneity of the sample.

In fact, we have chosen elastiviscoplastic behavior with isotropic hardening for cortical bone whilst the trabecular bone is assigned with an elastiviscoplastic behavior with kinematic hardening.

The Chapter 5 of Part IV, presents the application of numerical model put forth in Chapter 4 to experimental cases.

The experimental case 3 is peculiar, since the experimental conditions do not permit the sample to establish an equilibrium. Therefore, we show keen interest in simulating this behavior. In order to understand the role of trabecular bone, we have tested two different samples (i) with no trabecular bone and (ii) with trabecular bone ("Composite sample"). The sample with no trabecular bone is too rigid compared to its counterpart when subjected to sample impact velocities. This proves that the trabecular bone in fact plays an important role in absorbing the impact forces and energies.

The composite sample was then tested at different impact velocities. In all cases, the numerical model is capable of producing results similar to corresponding experimental cases, but with little difference. This difference is in fact due to the difference in the geometry of experimental and numerical sample.

The conclusions and perspectives presented in Part V, conclude this PhD thesis. Keywords: Dynamic Impacts, Ribs, Fracture, Split Hopkinson Pressure Bar, Three point bending, Finite Element simulations, HR-pQCT, Elastic-viscoplastic, Damage, LS-Dyna.

xiii xiv

Contents

Acknowledgments ii

Préambule vii

Prologue xiii

List Of Figures xxi

List Of Tables xxiv

Nomenclature xxix

List of Figures [START_REF] Keller | Predicting the compressive mechanical behavior of bone[END_REF] 1.6 Compressive strength as a function of apparent density, put forth by [START_REF] Rice | On the dependence of the elasticity and strength of cancellous bone on apparent density[END_REF] [START_REF] Keaveny | Strength of Trabecular Bone[END_REF] 1.9 Two variable models for variation of E Currey (1988) Les études sur les côtes humaines traitent soient de côtes entières, soient de segments représentatifs de la variation géométrique des côtes. Aujourd'hui il existe un nombre limité de travaux sur le comportement mécanique des côtes sous charge dynamique aux vitesses de déformation intermédiaires

(1 [s -1 ] ≤ ε ≤ 100 [s -1 ]).
Cela entraîne un manque de compréhension des phénomènes se produisant en régime dynamique, en particulier de la fracture. De plus, les modèles numériques de thorax humain sont limités par les contraintes liées à la représentation géométrique et la description fiable du comportement jusqu'à la rupture complète.

Parmi les techniques disponible aujourd'hui (IRM, DAX, CT, pQCT, HR-pQCT, etc.), le HR-pQCT convient à cette thèse car il donne une distribution volumique de la densité, des paramètres osseux spécifiques aux niveaux de la macro-structure et de la micro-structure. french

Abstract:

This chapter presents the literature study concerning this thesis. First of all, we present the anatomy of thoracic cage, the consequences of rib fractures, the biological and mechanical characteristics of bones. This is followed by a detailed introduction to existing studies on human rib and its biological surrogates and also the existing numerical models on human thorax and ribs. We finish this chapter with an introduction to different imaging techniques used in the field of biomechanics.

There is a huge risk of penetration of the bone fragments and injuries to heart and lungs in the case of rib fractures. Several authors have established different evaluation methods for thoracic injuries (for example the AIS thoracic injury scale). These methods require to take the whole thorax into account. As this thesis deals with isolated ribs, unfortunately, these evaluation methods could not be used.

A human rib is composed of two different bone tissues: cortical and trabecular.

These two tissues show different mechanical behavior, that arises due to the different in architectural arrangement. Moreover, the mechanical properties of bone tissues depend on age of the subject, the mineral composition, and the mechanical environment.

Very few studies could be found today concerning the fracture and behavior of human ribs. These studies concern either whole ribs or segments of ribs representing the geometrical variation of a rib. Yet, there is a void in studies concerning the behavior of human ribs under dynamic loads in the intermediate strain rate regime

(1 [s -1 ] ≤ ε ≤ 100 [s -1 ]
), proving the lack of human knowledge in this regime.

Chapter 1

Several numerical models of human thorax and ribs were developed in recent years.

These numerical models do come with their respective constraints. These so called constraints could be narrowed down to lack of geometrical representation and the use of reliable constitutive and damage law.

Among the imaging techniques used in the biomechanical field today (MRI, DAX, CT, pQCT, etc.), the HR-pQCT imaging technique is the most feasible for this thesis as it provides information on volumetric distribution of density and specific architectural parameters at both macro and micro scale.

Introduction to human thorax

The trunk of a human body is divided by the diaphragm into two parts, the upper part is called the thorax and the lower part is called the abdomen. The human thorax is supported by a skeletal framework called thoracic cage and contains vital respiratory and circulatory organs such as Lungs and Heart [START_REF] Chaurasia | Human Anatomy: Regional & Applied : Upper Limb & Thorax[END_REF]. The thoracic cage is an osseocartilaginous elastic cage (see Fig. 1.1a), which is designed for increasing and decreasing intrathoracic pressure during inspiration and expiration of breathing process. It is formed by:

1. Sternum at the anterior end (Fig. 1.1b) 2. 12 Thoracic vertebrae and intervening intervertebral discs at the posterior end 3. 12 ribs on either sides, that articulate with the thoracic vertebrae of vertebral column (Fig. 1.1a).

(a) Whilst, only 7 rib pairs articulate with the sternum at the anterior end, known as true ribs or vertebrosternal ribs (Fig. 1.1b).

(b) The 8 th ,9 th and 10 th ribs articulate with vertebral column at posterior end and with sternum at anterior end through costal cartilages. Thus they are known as false ribs or vertebrochondral ribs (Fig. 1.1b).

(c) The 11 th and 12 th ribs are known as floating ribs or vertebral ribs since they are free at the anterior end (Fig. 1.1b).

An interesting fact is that, the thoracic wall of infants is highly elastic hence, the rib fracture is rare. On the other hand, the rib fractures in adult may occur at its weakest point. The upper two ribs and the lower two ribs (floating ribs) are the least commonly fractured. [START_REF] Gray | Anatomy of the human body[END_REF] 1. Manubrium, is quadrilateral in shape. It is thick and broad above, narrow and thin at the junction with the Gladiolus. The Manubrium has an articular facet for the clavicle bone and a depression for the 1 st costal cartilage.

2. Gladiolus, is narrow, long and thin than the manubrium. It is the broadest close to its lower end (at facet for 5 th costal cartilage).

3. Xiphoid Process, is the smallest and thinnest of all the three components. It is less oscified when compared to other sternum components in adults, whilst it is a cartilaginous structure amongst the youth.

Thoracic Vertebrae

Section 1.1

The thoracic vertebrae are intermediate in size, placed in between the cervical and lumbar regions. These vertebrae are distinguished by the demi-facets on the sides of body for articulation with the ribs of thoracic cages. They are thinner at the upper extremity and increase in size downwards (see Fig. 1.2b).

Ribs

Rib are bony arches placed one below the other to form a thoracic cage. They can be classified as flat bones. An adult thoracic cage, consists of 12 pairs of ribs (as in Fig. 1.1a) that are classified into three categories:

• True ribs (1 to 7), are those which articulate with the sternum through corresponding coastal cartilage.

• False ribs (8 to 10), are those which articulate with the sternum but with only one single cartilage.

• Floating ribs (11 and 12), are those which do not articulate with the sternum. They a free to move, therefore to "float".

Each rib is supposed to consist three anatomical regions as in Fig. 1.3 [START_REF] Gray | Anatomy of the human body[END_REF].

1. Head, is at the anterior side of the rib. It has two articular facets separated by a wedge. These articular facets help the rib to articulate with corresponding thoracic vertebrae and the one above it.

Chapter 1 2. Neck, in particular has no biological importance. It just connects the head with rib body a point called the "Tubercle". This tubercle has one articular facet, that articulates with the corresponding thoracic vertebrae.

3. Body or Shaft of the rib is flat and curved with external and internal surfaces.

The internal surface has a groove called "Costal Groove", which protects the intercostal vein, artery and nerve from external damage.

Figure 1.3 -Structure of a typical human rib [START_REF] Gray | Anatomy of the human body[END_REF] As seen in Fig. 1.1a, ribs are oriented at a certain angle to the horizontal reference plane. The first rib has peculiar characteristics when compared to its counterparts in terms of orientation as it has a higher angle. In morphological terms, each rib is considered to be different since, the lengths and angle of orientation are not the same. Another common feature of human ribs is the transition of cross-sectional geometry, which is of circular cross-section from the head of the rib to the neck and seems to undergo a transition to elliptical cross-section at the neck and maintain it throughout the body.

Clinical relevance of rib fractures

Rib fractures are most common thoracic injuries, which are principal causative factors in Road Traffic Accidents (RTAs) [START_REF] Cavanaugh | Injury and Response of the Thorax in Side Impact Cadaveric Tests[END_REF] and in High Speed Ballistic Impacts (HSBIs) [START_REF] Shen | Modified Bilston nonlinear viscoelastic model for finite element head injury studies[END_REF]. The Chest Trauma (ChT) is one of the main contributors to the development of complications and mortality amongst poly-trauma patients [START_REF] Horst | Thoracic trauma now and then: A 10 year experience from 16,773 severely injured patients[END_REF]. Among the organs of human thorax, lungs represent the target organ for secondary damage by post traumatic inflammation [START_REF] Mommsen | Productive capacity of alveolar macrophages and pulmonary organ damage after femoral fracture and hemorrhage in IL-6 knockout mice[END_REF]. Moreover, in some cases, lung injuries contribute to Multiple Organ Failure (MOF) and represents around 24% of late deaths after trauma [START_REF] Dewar | Changes in the epidemiology and prediction of multiple-organ failure after injury[END_REF]. For these reasons, ChT is considered to be the second most lethal right next to Head Trauma (HeT) [START_REF] Ludwig | Management of chest trauma[END_REF]. One of the most recent case studies on ChT due to impacts, states that a ChT may arise either from a penetrating injury or a blunt injury [START_REF] Arunan | Chest trauma[END_REF]. They have also investigated and put forth the evolution of critical condition of the victim with respect to time of accident. Evidences show that Blunt Thoracic Trauma (BTT) out ranks Penetrating Thoracic Trauma (PTT) since, it needs more attention and time for recovery [START_REF] Dewar | Changes in the epidemiology and prediction of multiple-organ failure after injury[END_REF], [START_REF] Arunan | Chest trauma[END_REF]. Thus, in this thesis, we concentrate on BTTs only.

Several authors have investigated BTT through case studies, most of them Chapter 1 state that RTAs is the common causative factor for BTTs [START_REF] Roth | Chest injuries in children: An analysis of 100 cases of blunt chest trauma from motor vehicle accidents[END_REF], [START_REF] Hamid | Blunt thoracic trauma[END_REF] and [START_REF] Blau | 4 -High-Energy, Blunt-Force Trauma: Motor Vehicle (MV) Incidents[END_REF]. Apart from RTAs, Industrial Accidents (IAs) and Falls are other causative factors. In fact, the injury mechanism plays an important role in the severity of injury and critical condition of victim [START_REF] Burnside | Blunt thoracic trauma[END_REF] and [START_REF] Fallouh | Blunt thoracic trauma[END_REF].

Assessment of an injury can be done using a scoring system, and there are a variety of scoring systems available today. These so called scoring systems can be classified on the basis of their predicting nature (i.e.) systems based on anatomical description of injury, systems based on physiological parameters and systems combining both anatomical and physiological factors. The later being the most commonly used in predicting of outcome after trauma [START_REF] Chawda | Predicting outcome after multiple trauma: which scoring system?[END_REF]. The Tab. 1.1 presents the different trauma scoring systems available such as • Combined scoring systems (a) Trauma and Injury Severity Score (TRISS).

• Other scoring systems (a) Glasgow Coma Score (GCS) (b) Thoracic Trauma Severity Score (TTSS.)

Bone biology

Bone is one of the building block of the human musculoskeletal system, that takes part in several roles such as protection, transmission of mechanical forces and serves as a mineral reservoir. Moreover, it is a determinant of the global structural stiffness.

Bone tissue is composed of organic materials, inorganic materials and water. The following minerals are the subcategories of above mentioned materials:

1. Organic materials (a) Type I collagen molecules in high concentrations (b) Type III, VI and V collagen molecules in less amounts [START_REF] Keaveny | Bone Mechanics[END_REF][START_REF] Reznikov | Bone hierarchical structure in three dimensions[END_REF] (c) Noncollagenous proteins such as osteocalcin, osteonectin, osteopontin and sialopretein.

2. Inorganic materials (a) Ceramic crystalline material in impure form, known as hydroxyapetite (Ca 10 (PO 4 ) 6 (OH) 2 ). It is considered to be impure because of the existence of potassium, magnesium, strontium, carbonate and chloride/fluoride. [START_REF] Lowenstam | On biomineralization[END_REF] 3. Water in different forms. It is bound to crystal surface at microscopic levels, in between the triple helix collagen molecules and also in fluid form. Water present in these three forms is known as bound water. Moreover, water can also be found in unbound form in canaliculi, lacunae and blood vessels.

The interesting properties of bone is due to its hierarchical structure, which is of nine levels according to [START_REF] Reznikov | Bone hierarchical structure in three dimensions[END_REF] as in Tab. 1.2. The materials and tissue array at each level is responsible for forming different tissues. For example, the lamellar packets of ordered and disordered units form the trabeculae, which is the structural element of the trabecular bone. Whereas, the osteon and fibrolamellar units form the cortical bone [START_REF] Reznikov | Bone hierarchical structure in three dimensions[END_REF]. Thereby, forming a composite structure seen at IX hierarchical level, referred as 'bone' . Thus, at a macroscopic level, bone is non-homogeneous anisotropic and porous with two different tissues namely, cortical and trabecular.

It is clearly seen in Fig. 1.4 that a rib is composed by two different tissues: cortical and trabecular bone. Therefore, we are interested in biological and mechanical characteristics of the level VIII tissues (Tab. 1.2) in this thesis.

During their lifetime, cortical and cacellous bones are formed two different building blocks, namely woven and lamellar bone (see level VI of tab. 1.2 and Fig. 1.5).

Both bones are composed only of woven bone during the embryo stage and initial stages of development (age < 5 yrs) as well as during the fracture healing process. Later, the pockets of woven bones are replaced by lamellar bone. A woven bone is formed rapidly, poorly organized with random arrangement of collagen fibers and mineral crystals and is considered to be a transient tissue. [START_REF] Reznikov | Bone hierarchical structure in three dimensions[END_REF] Figure 1.5 -Woven (to the left) and Lamellar (to the right) bone

Another aspect of bone is its self repairing through a process called remodeling. Bone remodeling occurs on surface of trabeculae of cancellous bone and at haversian system of cortical bone. The modeling and remodeling process of bones is carried out by a group of bone cells (osteoblasts, osteoclasts, bone lining cells and osteocytes), known as "Basic Multicellular Units" (BMU's).

Level VIII tissues: cortical and trabecular

Cortical bone is in general more mineralized and less porous structure (5% ≤ P o ≤ 30%) . In the case of an healthy adult human, cortical bone exists in form of osteonal or Haversian bone i.e. an arrangement of 10-15 lamellae in concentric cylinders around the Haversian canal. These canals carry several blood vessels, nerves and blood cells which aide in communication with bone cells (see Fig. 1.6a). Since the bone tissue has the capacity of self repair through bone remodeling, the newly formed osteons are separated from the interstitial bone by a thin layer, known as cement line. This cement line forms a weak interface and is known to improve fatigue properties of cortical bone through microcraks, which server as reserves for energy dissipation during crack propagation [START_REF] Keaveny | Bone Mechanics[END_REF]. The most distinguished characteristic of trabecular bone is its porosity resulting from voids external to mineralized trabeculae, known as Inter-Trabecular space porosity (ITP), represented by voids in Fig. 1.7. These so called 'trabeculae' exist either in forms of rods or plates, depending on the need for force absorption. For example, the trabecuale in the human femoral neck (as in Fig. 1.7c) are in form of plates where as they exist in form of both rods and plates in the human vertebrae (as in

Mechanical insights on bone

Mechanical behavior of cortical and trabecular bone

We had seen the variables influencing mechanical properties of bone in Sec. 1.3.

Most of the research data presented were either on isolated cortical or on trabecular bone samples. These data suggests that there is an essential difference between the mechanical behavior of these two tissues. Quite a few authors suggest that these are two different tissues and thus the difference in mechanical behavior. On the other hand, these two tissues are based on the same morphological tissue, but the difference in mechanical behavior is due to the difference in architectural arrangement.

For example, cortical bone shows a behavior corresponding to elastic-plastic where as the trabecular bone shows behavior similar to that of polymeric foams.

Chapter 1 (a) Cortical bone [START_REF] Crowninshield | The response of compact bone in tension at various strain rates[END_REF] (b) Trabecular bone [START_REF] Carter | The compressive behavior of bone as a two-phase porous structure[END_REF] Figure 1.8 -Mechanical behavior of cortical and trabecular bone Another main difference between the mechanical behavior of cortical and trabecular bone is in their post yield response to loading-reloading tests.

Cortical bone

1. It exhibits mechanical properties degradation due to damage accumulation, thus it's post yield response is damage related.

2. It has a relatively low yield offset point and can be determined by standard yield offset technique (0.2%).

3. When reloaded beyond yield point, the initial modulus is not the same as its intact modulus (see Fig. 1.9a) .

4. According to [START_REF] Keaveny | Bone Mechanics[END_REF], the microdamage patterns include microcracks in longitudinal and transversal directions, cross-hatched shear band patterns and damage diffusion.

2. It starts to yield at ≈ 1% and can sustain large deformations upto 50%.

3. It reloading modulus is close to intact modulus but loses its stiffness rapidly (see Fig. 1.9b) .

4. This does not depend on volume fraction since similar response can be seen in bovine trabecular bone (as in Fig. 1.7a [START_REF] Keaveny | Bone Mechanics[END_REF].

5. Therefore, there is an evidence of residual strain and back stress.

These data suggest that cortical bone can be modeled as transverse isotropic material with isotropic hardening and trabecular bone can be modeled as an anisotropic material with kinematic hardening.

In either cases simple yield criteria such as von-Mises yield cannot be applied for cortical and trabecular bone materials, since they show evidences of more strength in compression with respect to tension. Moreover, damage seems to play an important role in post yield behavior, thereby suggesting the use of damage dependent yield and progressive failure criteria.

Variables influencing mechanical properties of bone

In order to understand the behavior of bone under certain mechanical conditions, it is important to quantify its mechanical properties and corresponding variables influencing it. Since bone can adapt it self by remodeling, its properties change with Chapter 1 respect to age, diseases and the mechanical environment it is subjected to. Therefore, it is important to quantify the mechanical properties before understanding its mechanical behavior.

Age

As said above, certain authors have investigated the influence of age on the variation of mechanical properties (Young's modulus (E) , ultimate strength (σ ult ) , energy absorbed at failure (W f ) and strain at failure (ε f )). For example, [START_REF] Currey | Changes in the impact energy absorption of bone with age[END_REF] had postulated that the energy absorbed at failure decrease with increase in age. This is since the ash content of bone increases with increase in age. Later, [START_REF] Currey | The effects of ageing and changes in mineral content in degrading the toughness of human femora[END_REF] had stated that "there is a strong decrease in properties with increase in ash content of bone, but other factors such as age are probably also responsible for deterioration of bone properties". In contrast, [START_REF] Leng | Effect of age on mechanical properties of Collagen phase in different orientations of human cortical bone[END_REF] have shown that E and σ ult (tensile) are age and orientation dependent i.e. there was a decrease in these properties in longitudinal and transversal directions with ageing as opposed to [START_REF] Mccalden | Age-related changes in the tensile properties of cortical bone. The relative importance of changes in porosity, mineralization, and microstructure[END_REF] and [START_REF] Wang | Age-related changes in the collagen network and toughness of bone[END_REF]. This peculiar behavior was correlated to the collagen matrix in cortical bone. However, the effects of ageing on mechanical properties of bone are still ambiguous.

Bone composition

The next parameter supposed to influence the mechanical properties of bone is its composition. A question often posed is which density is the appropriate independent variable characterizing mechanical properties of bone ?. There are different measures of densities, which are either measured directly or indirectly correlated with QCT images.

1. Apparent dry density (ρ dry )1 , is the ratio of dry specimen mass to specimen volume ρ dry = Dry specimen mass specimen volume (1.1)

2. Apparent wet density (ρ wet ), is the ratio of hydrated specimen mass to specimen volume ρ wet = hydrated specimen mass specimen volume (1.2) 3. Apparent ash density (ρ α ), is expressed in terms of ash fraction and apparent dry density as ρ α = αρ dry (1.3) where α is the ash content of a bone, defined as the ratio of ash weight to dry specimen weight.

The studies concerning the influence of bone composition on it mechanical properties could be grouped as per the variables used.

Group 1: ash fraction (α)

Two most cited studies in this group are [START_REF] Vose | Bone strength-its relationship to X-raydetermined ash content[END_REF] Quite a few other authors have had investigated the influence of apparent ash density ρ α on E and σ ult . Among the authors of this group, [START_REF] Keller | Predicting the compressive mechanical behavior of bone[END_REF] and [START_REF] Keyak | Correlations between orthogonal mechanical properties and density of trabecular bone: use of different densitometric measures[END_REF] [START_REF] Keller | Predicting the compressive mechanical behavior of bone[END_REF] used pooled data of cortical and trabecular bone samples from human vertebrae and femur, leading to a conclusion that ρ α explained over 96% of variations in mechanical behavior. Whereas [START_REF] Keyak | Correlations between orthogonal mechanical properties and density of trabecular bone: use of different densitometric measures[END_REF] used trabecular bone samples from human proximal tibia, leading to the conclusion that ρ α had the strongest correlation.

Group 3: apparent density (ρ app )

Chapter 1

Apart from mineral content, apparent density was also considered as control variable by several authors. [START_REF] Galante | Physical properties of trabecular bone[END_REF], had investigated the effect of apparent density ρ wet on compressive strength σ ult of cancellous bone sample extracted from human vertebrae. They concluded that "porosity P o is a major factor in determining ρ wet , which in turn influences σ ult ". Moreover, their tests proved that σ ult showed a positive linear correlation with ρ wet . This was backed up by [START_REF] Carter | The compressive behavior of bone as a two-phase porous material[END_REF], who had proved that E ∝ ρ wet 3 and

σ ult ∝ ρ wet 2 .
Rice et al. (1988) put forth a pooled data in concern of the dependence of Elastic Modulii and strength on the apparent density. For their Young's Modulus estimates, the authors had used a pooled data from [START_REF] Carter | The compressive behavior of bone as a two-phase porous material[END_REF][START_REF] Bensusan | Tensile, Compressive and torsional testing of Cancellous bone[END_REF][START_REF] Williams | Properties and an anisotropic model of cancellous bone from the proximal tibial epiphysis[END_REF] [START_REF] Rice | On the dependence of the elasticity and strength of cancellous bone on apparent density[END_REF] In the case of compressive strength, the data were taken from [START_REF] Galante | Physical properties of trabecular bone[END_REF], [START_REF] Behrens | Variations in strength and structure of cancellous bone at the knee[END_REF], [START_REF] Carter | The compressive behavior of bone as a two-phase porous structure[END_REF], [START_REF] Williams | Properties and an anisotropic model of cancellous bone from the proximal tibial epiphysis[END_REF] and [START_REF] Bensusan | Tensile, Compressive and torsional testing of Cancellous bone[END_REF] [START_REF] Rice | On the dependence of the elasticity and strength of cancellous bone on apparent density[END_REF] As both Young's Modulus and Strength (Compressive & Tensile) are dependent on the square of apparent density, they are proportional to each other. This proportionality was also substantiated by [START_REF] Currey | The effects of ageing and changes in mineral content in degrading the toughness of human femora[END_REF], [START_REF] Brown | Mechanical property distribution in the cancellous bone of Human Proxima Femur[END_REF] and [START_REF] Bensusan | Tensile, Compressive and torsional testing of Cancellous bone[END_REF].

Apart from this, they also proposed extrapolated equations of dependence of mechanical properties of Compact bone on apparent density in contradiction to Wolff's statement: "Compact bone is simply more dense than Cancellous Bone" Whilst Galante et al. (1970), Carter.D.R and[START_REF] Carter | The compressive behavior of bone as a two-phase porous material[END_REF], [START_REF] Gibson | The mechanical behaviour of cancellous bone[END_REF] and [START_REF] Rice | On the dependence of the elasticity and strength of cancellous bone on apparent density[END_REF] have investigated the isotropic behavior of cortical and cacellous bones, some authors have also presented evidences of transverse isotropic behavior of bone. For example, [START_REF] Lotz | Mechanical properties of metaphyseal bone in the proximal femur[END_REF] [START_REF] Lotz | Mechanical properties of metaphyseal bone in the proximal femur[END_REF] They had concluded that "The difference between global average values of longitudinal and transverse properties is due to the difference in density and decrease in transverse anisotropy." Similar to [START_REF] Lotz | Mechanical properties of metaphyseal bone in the proximal femur[END_REF], [START_REF] Keaveny | Strength of Trabecular Bone[END_REF] had presented the variation of properties in longitudinal and transverse directions for trabecular bone as follows:

Property

Orientation Equation [START_REF] Keaveny | Strength of Trabecular Bone[END_REF] On the other hand, [START_REF] Gibson | The mechanical behaviour of cancellous bone[END_REF] stated that cancellous bone is a cellular structure connected with a network of rods or plates. She had quantified the dependence of mechanical properties to relative density ρ rel , expressed as follows:

ρ rel = ρ bone /ρ solid (1.4)
where ρ bone is the density of the cancellous bone sample and ρ solid is the density of solid bone in the cell wall. Based on ρ rel , she had classified the type of cellular structure exhibited by cancellous bone.

Cellular structure

   Open cell (rod like) if ρ rel < 0.2 Closed cell (plate like) if ρ rel > 0.2
Each cellular structure appears to exhibit either symmetry or no symmetry, therefore, resulting in 4 structures (asymmetric cells with open cell structure, asymmetric cells with closed cell structure, columnar cells with open cell structure and columnar cells with closed cell structure).

She had then performed compression tests on samples in order to correlate mechanical property, mechanism of deformation with ρ rel .

• Asymmetric cells:

-with open cell structure showed that the young's modulus E and ultimate stress in compression σ ult are proportional to the square of relative density

(∝ ρ rel 2 ).
-with closed cell structure showed that the young's modulus E and ultimate stress in compression σ ult are proportional to the cube of relative density (∝ ρ rel 3 ).

• Columnar cells:.

-For both open and closed cells, the young's modulus and ultimate stress in compression along longitudinal direction were proportional to relative density (∝ ρ rel ).

Group 4: Two variables

The models presented above are single parameter models i.e. they consider either mineral composition (α) or apparent density (ρ app ) as control variables. On the other hand, several authors have proposed two parameter models. For example, Currey (1988) had proposed both power law and linear regression models to explain the variation of Young's modulus (E) with respect to calcium content (Ca mg.g -1 ) and Section 1.3 volume fraction (V f %) of cortical bone samples extracted from various mammals and birds. The volume fraction (V f ) is expressed in terms of ρ app and ρ tissue as

V f = ρ app ρ tissue (1.5)
He had concluded that "Young's modulus shows a positive correlation with both calcium content and volume fraction" and "power law regression seems to fit the data better than linear equations". His results are given in Tab. 1.9

Model Tension

Bending

Linear E = -71.2 + 0.197Ca + 42.5V f E = -52.3 + 0.143Ca + 35.3V f Power law log 10 E = -6.30 + 3.17log 10 Ca + 3.52log 10 V f log 10 E = -4.84 + 2.55log 10 Ca + 3.13log 10 V f
Table 1.9 -Two variable models for variation of E Currey (1988) Although the models presented by [START_REF] Currey | The effect of porosity and mineral content on the Young's modulus of elasticity of compact bone[END_REF], [START_REF] Keller | Predicting the compressive mechanical behavior of bone[END_REF] and [START_REF] Keyak | Correlations between orthogonal mechanical properties and density of trabecular bone: use of different densitometric measures[END_REF] were capable of explaining the fact that the mechanical properties of bone depend both on ρ α and ρ app , their models did not separate the effects of each individual variable. This so called drawback was brought into light by [START_REF] Hernandez | The influence of bone volume fraction and ash fraction on bone strength and modulus[END_REF], who had formed an equation for true tissue density ρ tissue based on data of [START_REF] Keller | Predicting the compressive mechanical behavior of bone[END_REF] as follows

ρ tissue ≈ 1.41 + 1.29α (1.6)
V f can be expressed in terms of α as

V f = ρ app 1.41 + 1.29α (1.7)
Based on α and v f from Eq. (1.7), they had proposed a two variable power law regression models for E and σ ult

E [GPa] = 84.37V f 2.58 α 2.74 (1.8) σ ult [MPa] = 794.33V f 1.92 α 2.79 (1.9)
We can notice that most of the regression models presented above are either based on uniaxial compression or three point bending, as the tensile properties are usually expressed in percentages of compressional properties. Several authors have investigated and presented different values for the ratio of tensile to compressive properties. In fact, this ratio depends on the origin of the sample, test direction and type of Chapter 1 bone tissue tested [START_REF] Reilly | The elastic and ultimate properties of compact bone tissue[END_REF]. A list of studies in concern with the ratio between tensile and compressive properties is given in Tab. 1.10. Even though this ratio seems to vary between 0.33 and 1.27 depending on the anatomical location, type of bone tissue and origin of test sample, [START_REF] Keyak | Prediction of femoral fracture load using finite element models: an examination of stress-and strain-based failure theories[END_REF] had considered this ratio to be a constant and is best when considered between 0.7 and 1.0.

Cortical bone architecture

Cortical bone in general, is a solid structure with osteons and canals immersed in its matrix. The orientation of these osteons describes the anisotropic nature of cortical bone. Since the principle direction of these osteons is along the longitudinal axis, cortical bone is more strong along longitudinal direction than transverse direction. Apart from osteons, the circumferential lamallae (as in Fig. 1.5) do exhibit different material properties in radial and tangential directions. According to [START_REF] Reilly | The elastic and ultimate properties of compact bone tissue[END_REF], this structure is influenced by secondary osteons, therefore cortical bone has transverse isotropic nature.

Trabecular bone architecture

In contradiction to cortical bone, trabecular bone is dominated by voids with little space occupied by trabeculae. The anisotropic nature of trabecular bone is better described by the orientation of these trabeculae. Several authors have presented different hypothesis regarding the influence of trabecular orientation on its mechanical properties. [START_REF] Whitehouse | A stereological method for calculating internal surface areas in structures which have become anisotropic as the result of linear expansions or contractions[END_REF] was the first to propose a methodology to measure the microscopic anisotropic structure of trabecular bone in terms of distance between bone marrow and trabeculae, known as 'Mean Intercept Length (MIL)' . Later, [START_REF] Cowin | The relationship between the elasticity tensor and the fabric tensor[END_REF] proposed a relation between the elastic tensor and a second order tensor called fabric tensor. The Fabric tensor is expressed as the inverse square root of the MIL tensor. The major assumption considered by [START_REF] Cowin | The relationship between the elasticity tensor and the fabric tensor[END_REF] is the matrix material is isotropic in order to calculate the anistropic nature arising from different orientations of the matrix materials (trabeculae). This hypothesis was used to put forth an anisotropic strength criterion depending on the fabric tensor by [START_REF] Cowin | Fabric dependence of an anisotropic strength criterion[END_REF]. The fabric dependence anisotropic nature of trabecular bone was then validated by [START_REF] Turner | The fabric dependence of the orthotropic elastic constants of cancellous bone[END_REF].

Mechanical environment

Other parameters influencing the mechanical properties of bone is the mechanical environment in which it is supposed to function. For example, bone behaves differently in static, quasi-static and dynamic environments, therefore, there is an influence of strain rate ( ε).

McElhaney (1966) was probably the first and the most cited article in this concern. He had studied these effects through compression test on bone samples form a 3 year old bovine femur and a 24 years old male femur using an air gun under ε ∈ [0.001 s -1 , 1500 s -1 ]. Properties like ultimate stress (σ ult ), energy absorption capacity (e c ) , Young's modulus (E), strain at failure (ε f ) were studied. The variation in these properties were studied in terms of a critical velocity (v c )3 . In case of bone, v c occurred when ε ∈ [0.1 s -1 , 1 s -1 ]. He had also pointed out that ε also had an influence on fracture mechanism. The average values of these properties are given in Tab. 1.11. It can be noticed that as ε increased, e c begins to increase in the begining and shows a drastic decrease for ε

∈ [0.1 s -1 , 1 s -1 ]. On the other hand, ε f tends to decrease but shows a drastic decrease for ε ∈ [0.1 s -1 , 1 s -1 ]
. This shows that the critical velocity v c for bovine femur bone exists for ε

∈ [0.1 s -1 , 1 s -1 ].
Another peculiar results to be noted from the results of [START_REF] Mcelhaney | Dynamic response of bone and muscle tissue[END_REF] is that both σ ult and E continued to increase with increasing ε. Moreover, McElhaney (1966) had stated that stress at a particular strain level, can be related to the corresponding strain rate with an exponential as in Eq. 1.10.

σ = ALn( ε) + B (1.10)
where, A and B are constants determined from experimental data. A = 2.96 × 10 7 Pa.s -1 and A = 2.68 × 10 8 Pa for bovine femur bone. 4

Strain rate ε, [s [START_REF] Mcelhaney | Dynamic response of bone and muscle tissue[END_REF] Increase in young's modulus with increasing strain rate was also reported by [START_REF] Crowninshield | The response of compact bone in tension at various strain rates[END_REF], [START_REF] Currey | The effects of strain rate, reconstruction and mineral content on some mechanical properties of bovine bone[END_REF], [START_REF] Wright | Tensile testing of bone over a wide range of strain rates: effects of strain rate, microstructure and density[END_REF] and [START_REF] Hansen | The Effect of Strain Rate on the Mechanical Properties of Human Cortical Bone[END_REF]. On the other hand, a decreasing effect was reported by [START_REF] Ferreira | Mechanical properties of bovine cortical bone at high strain rate[END_REF] whilst [START_REF] Pithioux | Comparison of compact bone failure under two different loadings rates: experimental and modelling approaches[END_REF] had reported no correlation.

-1 ] σ ult , [MPa] e c , [J.m -3 ] E, [GPa] ε f , [%] 0.
Similar to Young's modulus, several authors have reported the effect of strain rate on post yield behavior (yield stress and yield strain). For example, [START_REF] Crowninshield | The response of compact bone in tension at various strain rates[END_REF] and [START_REF] Currey | The effects of strain rate, reconstruction and mineral content on some mechanical properties of bovine bone[END_REF] have reported an increasing correlation with increasing strain rate as opposed to [START_REF] Hansen | The Effect of Strain Rate on the Mechanical Properties of Human Cortical Bone[END_REF], who had observed a decreasing correlation. [START_REF] Mcelhaney | Dynamic response of bone and muscle tissue[END_REF], [START_REF] Crowninshield | The response of compact bone in tension at various strain rates[END_REF], [START_REF] Pithioux | Comparison of compact bone failure under two different loadings rates: experimental and modelling approaches[END_REF] have reported a decreasing trend in fracture strain with respect to increasing strain rate for cortical bone. [START_REF] Guedes | Viscoelastic behaviour and failure of bovine cancellous bone under constant strain rate[END_REF] had reported an increasing trend in fracture strain with respect to increasing strain rate for trabecular bone.

Whilst the authors mentioned above have reported the effects of strain rate on isolated cortical and trabecular sample, [START_REF] Wood | Dynamic response of human cranial bone[END_REF] had investigated the behavior of human cranial bone samples extracted from male subjects with age between 25 4 note: Ln denotes natural logarithm Section 1.4 yrs to 95 yrs (with bith cortical and trabecular). These samples were subjected to a tensile load for ε ∈ [0.005 s -1 , 150 s -1 ]. Properties like ultimate stress (σ ult ), energy absorption capacity (e c ), Young's modulus (E), strain at failure (ε f ) were studied. The regression models proposed by [START_REF] Wood | Dynamic response of human cranial bone[END_REF] , 1971) He concludes that:

1. human cranial bone exhibits transverse isotropic behavior.

2. young's modulus E and ultimate tensile strength σ ult tend to increase with increasing strain rate.

3. failure strain ε f decreases with increasing strain rate.

4. energy absorbed at failure e c , on the other hand seems not to be influenced by strain rate.

Numerous studies were presented in Sec. 1.3 concerning the mechanical behavior of bones. It is noted that these studies were based on small sized bone samples extracted from bone organs. The reason for samples of such sizes is the nature and aim of research studies mentioned above i.e. to characterize and to apprehend the mechanical behavior under specific mechanical environments.

Despite these studies, quite a few authors have also studied the behavior of whole bone organs to characterize the response of the structure. In such cases, small sized samples do not give good insights. Therefore, we are obliged to use samples of decent size and different materials. This chapter, presents existing studies on ribs at Level IX of Tab. 1.2 or samples with both cortical and trabecular bone.

Chapter 1

Experimental studies on structural response of ribs

In case of structural response characterization, it is logical to use either whole rib samples or an anatomical segment. Therefore, the tests on structural response can be categorized as:

1. Three point bending tests on whole or partial rib samples 2. compression or tensile tests on mutli-material cylindrical samples of decent size.

Quasi-Static tests

Three point bending [START_REF] Granik | Human ribs: Static testing as a promising medical application[END_REF] is one of the most cited articles regarding the fracture of ribs and their response in static three point bending. They had tested segments from 6 th and 7 th human ribs (see Fig. 1.10a). The segments of ribs with least apparent curvature were selected (4 in long). Their data suggest that the thickness distribution of cortical bone is non-uniform and the average dimensions of the cross section at point of loading were 0.650inches × 0.290inches as in Fig. 1.10b. [START_REF] Granik | Human ribs: Static testing as a promising medical application[END_REF] They have concluded that the skew and taper are not important and did not influence the mechanical properties of the rib segment. Moreover, the trabecular bone of a rib shows negligible contribution to the mechanical properties of rib, since its stiffness and strength are much less when compared to that of whole structure.

The stiffness (eq. 1.11) and strength (eq. 1.12) of rib sections were calculated as follows

E f RS = ∆P L 3 ∆δ48I (1.11) σ f RS = P max Lc 4I (1.12) E ′ f RS = ∆P L 3 ∆δ768 10 I + 3 I b + 3 I c (1.13)
where, ∆P is the change in load (expressed in pounds lbs), ∆δ is the change in deflection (expressed in inches), L is the segment length (expressed in inches), I, I b and I c are area moments of inertia at point of loading, L/4 and 3L/4 respectively (inches 4 ). [START_REF] Got | Les Caractéristiques Morphologiques, Chimiques et Physiques des côtes et leurs relations avec la déflexion provoquée du thorax[END_REF] have tested 80mm long straight bone segments from 4 th ,5 th and 6 th human ribs. They have presented correlations between ash fraction per unit mass of the segment (A s /M s ), ash content per unit length of the segment (A s /L s ), age and sex.

• Their results shows that the slope of elastic region in quasi-static regime is strongly correlated with the ash content per unit length (A s /L s ), thickness of cortical bone and cortical bone surface area.

• Where as the dynamic three point bending tests prove that the maximum force in force-displacement response is strongly correlated to (A s /L s )

In contrast to [START_REF] Granik | Human ribs: Static testing as a promising medical application[END_REF] and [START_REF] Got | Les Caractéristiques Morphologiques, Chimiques et Physiques des côtes et leurs relations avec la déflexion provoquée du thorax[END_REF], [START_REF] Stein | Rib structure and bending strength: An autopsy study[END_REF] have studied the response of human rib segments from 79 male cadavers with respect to age and possible disease. They concluded that "apart from myeloma and arthritis, none of the other diseases had an effect on the response."

Another commonly cited work on quasistatic three point bending behavior of human ribs is [START_REF] Charpail | Analyse du comportement mécanique des côtes humaines en dynamique[END_REF], who had studied the behavior of 6 th ,7 th and 8 th human ribs. The experimental setup used is shown in Fig. 1.11. The behavior of rib was characterized through the following variables

• t max , is the time corresponding to maximum force observed • t end , is the time at which a rapid decrease in force is observed, indicating rupture

• F max , is the maximum force observed

• δ max , is the rib displacement corresponding to maximum force observed

• K, is the slope of linear part of the force displacement curve Figure 1.11 -Experimental setup [START_REF] Charpail | Analyse du comportement mécanique des côtes humaines en dynamique[END_REF] According to [START_REF] Charpail | Analyse du comportement mécanique des côtes humaines en dynamique[END_REF], the mechanism of rupture under quasi-static three point bending can be explained in 5 stages: The tests carried out by [START_REF] Charpail | Analyse du comportement mécanique des côtes humaines en dynamique[END_REF] have led to the following conclusions :

• There is a difference between the response of female and male ribs. This difference is due to the variation in geometry of corresponding rib. A similar conclusion was also put forth by [START_REF] Kimpara | Biomechanical properties of the Male and Female Chest Subjected to Frontal and Lateral Impacts[END_REF].

• The modulus and ultimate strength found by three point bending tests are ≈ 1.5 times greater than those found in compression tests. This can be explained by the difference in the behavior of bone subjected to uniaxial compression and tension respectively.

• The thickness of cortical bone at a particular cross-section plays an important role in bending stiffness and strength. "more the thickness, less are bending stiffness and strength".

The data for modulus and ultimate strength subjected to three point bending tests found in literature are presented in Fig. 1.13a and Fig. 1.13b respectively. This data is grouped with respect to the loading direction, LM Dy signifies Dynamic loading in Lateral Medial direction, LM QS signifies Quasi Static loading in Lateral Medial direction, ML QS signifies Quasi Static loading in Medial Lateral direction and NR signifies Not Reported by the author or ambiguous. 1.13 -Existing data on human ribs under three point bending [START_REF] Mayeur | Personnalisation géométrique et mécanique multiéchelles du thorax humain[END_REF] The studies presented above have used human ribs for experimental characterization. In case of constraints to obtain human ribs, biological surrogates can be used. One of the recent studies concerning quasi static three point bending had used porcine ribs as substitute [START_REF] Aubert | Étude du comportement des côtes sous sollicitations dynamiques[END_REF].

The porcine rib was represented as a curve beam with an uniform circular cross section in order to calculate the Young's modulus through a curved beam subjected to small displacements hypothesis. They found that the average Young's modulus Section 1.4 of 6 porcine rib sample was around 14 GPa, though the Eq. 1.14

K pr = 2EI R 3 δ(α) (1.14)
The results states that on loading the same rib beyond a displacement of 4 mm, the initial elastic modulus is not the same. For example, the Fig. 1.14, shows the response of porcine rib to loading-unloading and reloading response for the case δ = 3.5mm in red and δ = 4mm in green respectively. It is clear that once the same rib is loaded at different displacements, the initial elastic response is almost the same since the curves overlap.

Figure 1.14δ ≤ 4mm [START_REF] Aubert | Étude du comportement des côtes sous sollicitations dynamiques[END_REF] On the other hand, if the displacement is greater than 4 mm, there is a shift in the initial elastic response as in Fig. 1.15. This change in response with respect to displacement is due to the accumulation of damage and plastification of the rib. In contradiction to [START_REF] Granik | Human ribs: Static testing as a promising medical application[END_REF], the response in Fig. 1.15 proves that there is an influence of trabecular bone on the response of whole rib.

Multimaterial cylindrical samples

The three point bending data presented above lead to a conclusion that the geometry of the rib and the segment has an influence on the structural response of a rib. Since it is difficult to obtain a straight and uniform sample to study the compressive and tensile behavior of rib. Thus, quite a few authors have suggested to conduct tests on cylindrical segments of ribs containing both cortical and trabecular bone. Quite a few studies can be found in the literature regarding uniaxial tests on human rib samples are shown in Fig. 1.16a and Fig. 1.16b Yet this proves that the sample location does play an important role for the response of human ribs. One of the studies concerning the dispersion of properties with respect to anatomical location of the sample was proposed by [START_REF] Charpail | Analyse du comportement mécanique des côtes humaines en dynamique[END_REF]. She concludes that the sample extracted from posterior end were more stronger than that from the anterior end. Therefore, showing the influence of porosity.

Dynamic tests

Most of the studies presented above, were based on quasi static test or dynamic tests using pendulum. In order to characterize the structural response of a rib under similar dynamic environments similar to automotive impacts, it is necessary to subject the rib to a three point bending load at high strain rates. This can be possible through the Split Hopkinson Pressure bar method.

The inconvenience in split hopkinson pressure bar is the choice of material for the bars so that the strain signals can be treatable for results. This depends on the mechanical impedance of the bar and sample materials. As we know that bone has very low mechanical impedance compared to metals, polymeric materials must be considered. The use of polymeric materials for testing low impedance materials in dynamic environment is not new, it was used by [START_REF] Zhao | On the use of SHPB techniques to determine the dynamic behavior of materials in the range of small strains[END_REF] and most recently by [START_REF] Pavier | Experimental study of the coupling parameters influencing the terminal effects of thoracic blunt ballistic impacts[END_REF] for low dense polymers. Unfortunately, there are only a few studies concerning three point bending tests on whole ribs using SHPB [START_REF] Aubert | On the use of Hopkinson bar-bending apparatus to study soft impact on porcine ribs[END_REF][START_REF] Pavier | Contribution à la compréhension des phénomènes physiques lors de l'impact d'un corps sur un modèle de structure biologique[END_REF]. Section 1.5 The experimental setup used by [START_REF] Pavier | Contribution à la compréhension des phénomènes physiques lors de l'impact d'un corps sur un modèle de structure biologique[END_REF] is considered as reference for this thesis and is presented in the next chapter.

Numerical models of human thorax and ribs

Several models and proposals were put forth for the human thorax. These models can be classified as either Isolated Thorax models or Full body human models. On the other hand, quite a few authors have also presented isolated rib models with both cortical and trabecular bones.

Human Thorax Models

The first finite element model of the human thorax was proposed by [START_REF] Roberts | On Some Geometric Properties of Human Ribs[END_REF], [START_REF] Roberts | Global geometric characteristics of typical human ribs[END_REF] and [START_REF] Sundaram | Finite element analysis of the human thorax[END_REF]. In these models, the ribs were represented as beam elements. [START_REF] Roberts | On Some Geometric Properties of Human Ribs[END_REF] is one of the most commonly cited models, which was based on postmortem studies on first 8 ribs from one female cadaver. Even though this model was a break through for its time in terms of numerical representation of a human thorax, there was no correlation with cadaver experimental database. Another drawback of this model is that it did not consist the intercostal muscle tissues and spine. This loophole was filled by [START_REF] Andriacchi | A model for studies of mechanical interactions between the human spine and rib cage[END_REF], who added the tissues and spine. The bones were modeled as rigid bodies and tissues of the thoracic cage as deformable beams. Later, [START_REF] Sundaram | Finite element analysis of the human thorax[END_REF] proposed a FE model for the thorax by considering the entire thoracic cage as deformable body. The thorax model proposed by [START_REF] Sundaram | Finite element analysis of the human thorax[END_REF], consisted of almost all the crucial components of the thoracic cage including lungs and heart. A major draw back of this model was that, it was not capable of simulating large displacements, non-linear behavior and elements were not fine enough.

Later with the improvement of computer graphics and competencies, several well refined models were proposed. Some of the most commonly used models in the industry are the HUman MOdel for Safety, HUMOS [START_REF] Robin | HUMOS: HUman MOdel for Safety -A joint effort towards the development of refined human like car occupant models[END_REF], the Total HUman Model for Safety THUMS [START_REF] Maeno | Development of a Finite Element Model of the Total Human Model for Safety (thums) and Application to Car-Pedestrian Impacts[END_REF], the RADIOSS human model [START_REF] Arnoux | Radioss finite element model of the Thor dummy[END_REF], the H-model [START_REF] Haug | Human Models for Crash and Impact Simulation[END_REF] and the Isolated thorax model [START_REF] Schoell | Age-and sexspecific thorax finite element model development and simulation[END_REF]. [START_REF] Mayeur | Personnalisation géométrique et mécanique multiéchelles du thorax humain[END_REF], had pointed out a few drawbacks and simplifications of most frequently used human models concerning the human ribs of these models as follows:

1. there is need for accurate geometrical representation of ribs 2. need for integrity between bones and other tissue 3. different properties are to be allocated to different ribs

Isolated Rib Models

Isolated rib models are the ones aimed in improving either the geometrical representation or mechanical behavior of the ribs in the human body models. Some of the most recent isolated human rib models are [START_REF] Charpail | Analyse du comportement mécanique des côtes humaines en dynamique[END_REF], [START_REF] Li | Rib fractures under anterior-posterior dynamic loads: Experimental and finite-element study[END_REF], [START_REF] Yates | Subject-Specific Modeling of Human Ribs: Finite Element Simulations of Rib Bending Tests, Mesh Sensitivity, Model Prediction with Data Derived From Coupon Tests[END_REF] and [START_REF] Lundin | Detailed FE rib modelling for fracture prediction[END_REF] 2018) were adapted by a homogenization method proposed by Amine (2017), which was not reported in their document. The authors presented above, have tested the mesh sensitivity and have proposed all-hex / hex-shell models for ribs. The most feasible model was supposed to be the shell-solid model i.e. the cortical bone was represented by shell elements and trabecular with solid elements. Among them, [START_REF] Charpail | Analyse du comportement mécanique des côtes humaines en dynamique[END_REF] tested the human rib with no trabecular bone under lateral medial loads (Fig. 1.17a). These results were used to optimize the properties of cortical bone in anterior-posterior loads with trabecular bone (Fig. Several imaging techniques are available and in use in the filed of clinical densitometry for analysis of bone strength and the parameters influencing it.

Speaking of which, the type of technique differs with the requirements of different applications. For example, techniques such as radiography, dual-energy X-ray absorptiometry (DXA), peripheral quantitative computed tomography (pQCT), computed tomography (CT) and the magnetic resonance imaging (MRI) are sufficient enough for analysis of macro-architectural properties. Micro-architectural properties in contrast to their counter part, require imaging techniques with more resolution such as multidetector CT, MRI, and high-resolution peripheral quantitative computed tomography (HR-pQCT). The above mentioned points were adapted from [START_REF] Griffith | Bone mass and architecture determination: state of the art[END_REF].

However, as said by [START_REF] Wong | A comparison of peripheral imaging technologies for bone and muscle quantification: a technical review of image acquisition[END_REF], the pQCT, the HR-pQCT and MRI techniques excel in depiction of bone architecture. He had also pointed out the limitation of MRI technique used for bone architecture that "one must be cautioned as to the accuracy of MRI measures of bone structure, especially when the in-plane pixel size is larger than 0.2 mm". The effects of different imaging techniques and their dosage levels on age of the subject were brought into light by [START_REF] Irogi | Update on bone density measurements and their interpretation in childrens and adolescents[END_REF] in their article concerning the use of these techniques for children and adolescents.

All together, we can narrow down to HR-pQCT imaging technique for its reliable and efficient depiction of both bone mineral density and its architecture. To understand the working principle of HR-pQCT, we shall first consider QCT. A typical QCT technique works on the principles of photoelectric absorption. AS said by [START_REF] Wong | A comparison of peripheral imaging technologies for bone and muscle quantification: a technical review of image acquisition[END_REF], thicker and denser the object, fewer the photons are detected. In fact the intensities of incident (I in ) and transmitted (I out ) photons are related to the thickness (th) and linear attenuation coefficient (µ) of the object examined as follows

I out I in = e µth (1.15)
As we know that the bone has the capacity to attenuate more photons than its surrounding tissues (muscles and fat for example), QCT scans are capable of quantifying theses materials separately. The HR-pQCT on the other hand, is capable of providing high resolution images of both trabecular and cortical bone architecture, whilst capable of given volumetric density distributions separately.

Summary

In a nutshell, despite the existence of huge research in concern to thoracic trauma assessment tools, the most common conclusion encountered is "None of the existing trauma scoring systems considers all the parameters such as variation in population, age , physiological variability,., etc and continuous validation is needed. Therefore, we must be cautious in the use of currently existing trauma scoring systems".

All these studies show that rib fractures due to impact are often associated with serious morbidity and mortality, which tend to increase with number of ribs fractured. Key medical complications associated with rib fractures are hemothorax, pneumothorax, pulmonary laceration, pulmonary contusion and vascular injury. The complication endangering a victim is fail chest, that arises due to multiple ribs / adjacent ribs fractured leading to free-floating and individual movements of ribs. Therefore it is important to either improve thoracic trauma prediction or at least to mitigate the effects in case of blunt chest trauma. A good understanding on the fracture of ribs under these specific mechanical loading conditions can lead us to address the aforementioned improvements. In order to have a good apprehension of rib fractures, it is necessary to give equal importance to the geometry of rib, the mechanical properties of its tissues (cortical and trabecular) and the mechanism Section 1.7 of fracture and the variables influencing it. The following chapters address these constraints from bio-medical as well as mechanical point of view.

Bone is a self repairing and evolving tissue. The complicated mechanical properties and behavior is due to the constant erosion of under developed tissues by mature tissues during initial stages of development that occurs at microscopic scale and lower levels of hierarchy.

A human rib consists of two different tissues, namely cortical and trabecular bone. The difference in mechanical behavior of these tissues is essentially due to the difference in their architecture and structure. Since both bone tissues are strong in compression when compared to tension, common yield criteria such as von-Mises yield criterion is not suitable for numerical models. Moreover, the behavior reported by [START_REF] Keaveny | Bone Mechanics[END_REF], suggests that damage accumulation has an important role in the determination of post yield behavior of these tissues.

Another interesting fact is that the mechanical properties are influenced by many variables such as age, bone composition (density and mineral content), architectural arrangement and also the mechanical environment. These data suggest that corresponding variables of each category must be given equal importance for a valid simulation. Since this thesis deals with high speed dynamic loading, the strain rate plays a key role in the behavior of numerical rib sample.

Therefore, we propose that elastic-viscoplastic behavior with isotropic hardening is suitable for cortical bone and an elastic-viscoplastic behavior with kinematic hardening for trabecular bone. Moreover, as damage accumulation is also involved in the post yield behavior, a strain rate sensitive external damage law based on the continuum damage mechanics is the one and only suitable law for a progressive fracture propagation is suggested.

The response of human rib as a structure is done through three point bending tests on either whole human ribs or segments of certain length. These tests prove that the geometry of the rib had a strong influence on the mechanical response. It was initially concluded that the trabecular bone did not change the response of a rib, since its rigidity and stiffness is low when compared to that of whole rib. Moreover, the thickness of the cortical bone played a key role in characterizing the response. The results of quasi static bending tests on porcine ribs show that the trabecular bone has an influence if the displacement of the rib exceeds a certain amount.

Chapter 1

In order to make sure that the geometry had an influence, studies on mutlimaterial cylindrical samples were carried out by certain authors. Their results prove that the anatomical location of the sample in a rib did influence the rigidity and stiffness of the sample. The fact that there is a huge dispersion in the properties of ribs is due to the difference in load direction, sample size and location and sex. Very few studies can be found regarding the response of human rib or its biological surrogates to dynamic three point lateral medial bending. When it comes to numerical modes of human thorax, there is always a question regarding the resolution and representation of its constituents. Even though these models are used quite often in the industry, there are still improvements to be made notably in 1. the geometrical representation 2. the mechanical constitutive law

the consideration of mechanical environment i.e. strain rate effects

To confront the problems of geometrical representation and mechanical behavior, quite a few authors have presented isolated rib models. Their results have proved that use of elastic plastic constitutive law for bone is feasible. Yet unfortunately, no studies can be found to consider the strain rate effects on damage accumulation and progressive crack propagation.

Therefore, the numerical model of a human rib or its surrogate must have a good geometrical representation, a constitutive law with strain rate effects and most importantly a damage law that is capable of considering the rate effects and progressive damage accumulation and crack propagation.

Part II

Analysis of experimental data obtained previously at the laboratory

Chapter 2

Processing of experimental data

Résumé :

Ce chapitre présente la caractérisation expérimentale des côtes de porc en flexion dynamique "trois points". La base de données expérimentale brute a été obtenue par [START_REF] Aubert | Étude du comportement des côtes sous sollicitations dynamiques[END_REF]. Les données ont été analysées pendant cette thèse. This chapter presents the mechanical characterization of porcine ribs subjected to dynamic three-point bending loads. The raw experimental data base was obtained prviously by [START_REF] Aubert | Étude du comportement des côtes sous sollicitations dynamiques[END_REF] A total of 20 test were analyzed. Three categories of responses were identified based on the force at interfaces: (i) tests without fracture, (ii) tests with fracture during dynamic equilibrium between input and output forces and (iii) fracture before dynamic equilibrium. We had noticed that the relative interface velocity, plays an important role in determining the type of rupture observed. As the relative velocity increases, the fracture tend to occur before dynamic equilibrium. Moreover, we had also noticed a decrease in time at rupture and displacement at rupture. This permitted us to establish a threshold in terms of energy at rupture, time at rupture and displacement rupture. If the sample is rigid enough to overcome this threshold, it shall establish the aforementioned dynamic equilibrium, which is characterized by a plateau kind response in the forces at interface. The complete rupture of the sample is characterized by a rapid decrease in input and output forces, and is due to the propagation of a crack from tensile surface towards the impact point (the zone in compression). The maximum values of force does not have an influence on the rupture. The ruptures observed are controlled by the local displacement of the sample.

Le banc expérimental est constitué par des barres de

Introduction to Hopkinson Bar

The Hopkinson pressure is one of the common methods to study the dynamic behavior of materials. It was first studied by John Hopkinson in the year 1872, when he investigated propagation of stress waves in iron wire. Later, [START_REF] Hopkinson | A method of measuring the pressure produced in the detonation of high, explosives or by the impact of bullets[END_REF] introduced a methodology to save pressure vs time graphs. The accuracy of this methodology was then improved by [START_REF] Davies | A critical study of the Hopkinson pressure bar[END_REF]. The Split Hopkinson pressure bar known as of today, was introduced by [START_REF] Kolsky | An Investigation of the Mechanical Properties of Materials at very High Rates of Loading[END_REF], he placed the specimen of interest in between an input bar and an output bar. The original SHPB setup is used for comprehension of material behavior subjected to dynamic compression. The SHPB setup had been modified for tests in tension and three-point bending. In general, the compression and tension test bench are used to test materials, whilst, the three-point bending test bench is used to understand the dynamic behavior of structures.

Principle of SHPB

The modified three bar SHPB setup used in dynamic testing is based on the Kolsky bar technique. The sample is placed in between input and output bar(s). A striker is launched at the input bar through a pressurized air cannon. Once the striker impacts the input bar, longitudinal compression waves travel though out the test bench. The waves involved in SHPB are as follows (see Fig.

Z b = A b E b ρ b (2.1)
The characterization of tests carried out on the SHPB is based on the propagation of elastic waves within the input bar and output bar(s) [START_REF] Zhao | A three dimensional analytical solution of the longitudinal wave propagation in an infinite linear viscoelastic cylindrical bar. Application to experimental techniques[END_REF], [START_REF] Zhao | On the use of a viscoelastic split hopkinson pressure bar[END_REF]. Another important characteristic that is mandatory to obtain decent results with SHPB is the slenderness ratio of the bars. This slenderness ratio helps us to approximate 1D wave propagation in bars, thereby avoiding the non plane wave effects.

To understand the physics involved in SHPB, let us consider a volume element in the bar. The dynamic equilibrium of this volume element is

∂σ(x, t) ∂x = ρ b ∂ 2 u(x, t) ∂t 2 (2.2)
where σ(x, t) is the stress in the volume element, that is represented as

σ(x, t) = E b ε(x, t) = E b ∂u(x, t) ∂x (2.3)
where, ε(x, t) represents the strain in the element and u(x, t) represents the displacement along x direction in the element. Substituting Eq. 2.3 in Eq. 2.2, we get the 'D'Alembert Equation' (Eq. 2.4) A general solution can be written to the D'Alembert equation (Eq. 2.4) using the principle of superposition as follows:

ρ b E b ∂ 2 u(x, t) ∂t 2 = ∂ 2 u(x, t) ∂x 2 (2.
u(x, t) = f (x -c b t) + g (x + c b t) (2.6)
where, f (xc b t) represents a wave propagating in positive x direction and g (x + c b t) represents a wave propagating in the opposite direction. In fact the strain ε(x, t), the stress σ(x, t) and the particle velocity v(x, t) in the element can be expressed in terms of first order spatial derivative of f (xc b t) and g (x + c b t) in spatial domain.

The expression of ε(x, t) shows that the strain is a linear combination of spatial derivative of the two waves (f and g) that propagate in opposite direction. Thus the strain in a section of the bar can be expressed in terms of incident wave ε inc and reflected wave ε ref as

ε(x, t) = ε inc + ε ref (2.7)
Similarly, we can write the input and output forces in terms of ε inc and ε ref as in Eq. 2.8 and Eq. 2.9 respectively.

F in = A b E b [ε inc + ε ref ]
(2.8)

F out = A b E b [ε tr ] (2.9)
On the other hand, the velocities of the input and output bar interfaces can be written in terms of c b , ε inc and ε ref as in Eq. 2.10 and Eq. 2.11 respectively.

v in = -c b [ε inc -ε ref ]
(2.10)

v out = -c b [ε tr ] (2.11)

Experimental Difficulties

Impedance adaptation

The characterization of the material response using SHPB is based on evaluation of force and displacement imposed on the sample at input and output interfaces at every time. Since the forces (Eq. 2.8, Eq. 2.9) and velocities (Eq. 2.10, Eq. 2.11) can be calculated though the strain signals, for an acceptable characterization, the mechanical impedance must be adapted between the bars and the sample. The mechanical impedance of the bar material (Z b ) and that of sample (Z s ) , controls the amplitude of reflected wave ε ref through Eq. 2.12 and transmitted wave ε tr through Eq. 2.13 respectively.

R = Z b -Z s Z b + Z s (2.12) T = 2Z s Z b + Z s (2.13)
where, Z b is the mechanical impedance of bar material and Z s is the mechanical impedance of sample.

Therefore, there is a strong influence of material used as bars in the characterization of sample response. Since, the testing material is bone and it not a uniform material i.e. some areas of bone are made of soft tissue whilst other are of hard tissue. Moreover, the mechanical impedance of bone is feeble. Thus the use of bars with high mechanical impedance such as steel or aluminum leads to incomprehension of different waves involved in SHPB tests and also high noise-to-signal ratio. Several authors had investigated the use of polymeric bars such as PMMA [START_REF] Zhao | On the use of a viscoelastic split hopkinson pressure bar[END_REF], [START_REF] Sawas | Dynamic characterization of compliant materials using an all-polymeric split Hopkinson bar[END_REF] and [START_REF] Chen | A split Hopkinson bar technique for low-impedance materials[END_REF]. Therefore, we'd chosen bars of Nylon to characterize the dynamic response of bone.

Shifting of waves in spatial and temporal domain

As said above, the forces and velocities are calculated at the interface, this poses another difficulty. The reason is, it is impossible to differentiate between the incident and reflected wave at any given instant at the interface due to complex wave interaction phenomena and the change in sign of these both waves. Therefore in practice, we capture the signal at a certain distance from the interface as in Fig. 2.2. Later these waves are shifted in such a way that they were captured at the interface. This methodology is perfect for bars made of elastic materials.

For example, let us consider that the strain signals acquired are in terms of absolute time frame (t ′ ). The shifting procedure consists of defining an origin time t ′ 0 which is the time of arrival for the incident pulse at the input bar sample interface. Then a relative time t is defined in such a way that the strain signals were acquired at the interfaces of bars and sample.

In case of polymeric bars, wave attenuation and wave dispersion phenomena due to the viscous property of these bars. According to [START_REF] Bacon | An experimental method for considering dispersion and attenuation in a viscoelastic Hopkinson bar[END_REF], in order to shift these signals in time, it is necessary to reconstruct these signal by taking the dispersion and attenuation into consideration. Moreover, [START_REF] Gary | Correction de dispersion pour l'analyse des petites déformation aux barres d'Hopkinson[END_REF] had proposed to introduce a viscoelastic correction factor based on the Pochammer-Chree model for an infinite cylinder. He had considered the polymeric bar as a linear isotropic visco-elastic material, that lead to an equation in frequency domain with a complex solution. [START_REF] Zhao | On the use of a viscoelastic split hopkinson pressure bar[END_REF] put forth two keen relations that were used in the reconstruction of these waves, (i) a relation between the frequency and phase velocity and (ii) a relation between frequency and damping of waves. We follow the same procedure of [START_REF] Zhao | On the use of a viscoelastic split hopkinson pressure bar[END_REF] i.e. firstly the measured signal is decomposed into a Fourier series, then the phases in the fourier series is modified by taking the dephasing in the wave due to transport into consideration. Lastly, the amplitude of the wave is modifies by taking damping into account.

Non-plane wave

Another well known issue in a SHPB setup is non plane wave effect. A non plane wave effect is nothing but a non uniform distribution of strain through out the cross section normal to the direction of propagation of waves. In other words, it is the possibility of having zero strain at the outer surfaces and a non-zero strain within the cross-section. This issue of critical importance was studied by several authors and one of the note works was proposed by [START_REF] Zhao | On the errors associated with the use of large diameter shpb, corrections for radially non-uniform distribution of stress and particle velocity in shpb testing[END_REF]. They have presented a detailed discussions on the non plane wave effects by taking the radius of the bar, the spectrum of pulse wave and the material of bar into consideration. According to them, the assumptions of uniform strain distribution in the bar for the experimental conditions in this thesis, leads to a maximum error of less than 5%. Therefore, the corrections applied in this study are only for dispersion but not for three dimensional Chapter 2 effects.

Experimental Setup at our Laboratory

The setup used in this thesis is shown in Fig. 2.2. It consists of a nylon striker of 15 mm in diameter and either 500 mm or 1300 mm long. The striker is launched at the input bar through an air cannon. The pressure in the air cannon varies between 1 bar and 5 bar. The velocity of the striker is measured using an optical method through reflective bands that are placed at regular intervals on the striker (silver bands on striker in Fig. 2.3).

Measurement of signals

It is established that the striker has reflective bands placed at certain intervals to reflect the laser so that the striker velocity can be measured. Once the striker passes the optical sensors, it triggers the measurement chain used for strain signals. The time is set to t ′ = 0, which is an absolute quantity. The strain gage on the input bar measures the incident and reflected strains vs time t ′ : ε inc (t ′ ) and ε ref (t ′ ) respectively. On the other hand, the two strain gages on the two output bars, measure the transmitted waves in absolute time t ′ : ε tr1 (t ′ ) and ε tr2 (t ′ ) respectively.

Figure 2.5 -Strain gage for signal measures

These gages are in a Wheatstone bridge arrangement so that one gage measures the longitudinal strain (ε) whilst the other measures the transverse strain signals (-νε). Another advantage of Wheatstone bridge arrangements is that it is sensible and it eventually eliminates the effects of bending in bar. Once the bridge is at equilibrium, we modify the gain of the system using an oscilloscope. This is done be calibrating the variation in resistance and the variation in the measured tension through:

∆Em = 1 4 ∆R 1 R 1 - ∆R 2 R 2 + ∆R 3 R 3 - ∆R 4 R 4 × V (2.14)
We know that the strain and the variation in resistance are related as

∆R 4 R 4 = k × (ε) (2.15)
where k is strain gage factor and for strain gages used in testing, k ≈ 2.

The signals are saved using National Instrument PXI-5105 (as in Fig. 2.6), with an acquisition memory of 12 bits and 8 ports. This helps us to save input strain signal, output strain signal and striker velocity data simultaneously. All these signals are saved with a sampling time interval of 1 × 10 -6 s and are triggered when

Preparation of rib

Porcine ribs were used as human surrogates for the tests carried out, it is quite common to use porcine ribs for studies in thoracic response as it shares geometrical and anatomical similarities with human ribs.

Fresh porcine ribs brought from a local butcher. The muscles and membranes surrounding the rib were removed. The extremities of the rib were then sanded to a planar surface in order to establish a good contact surface. Every rib specimen were then numbered and their characteristics such as mass, global dimensions were then noted. All the ribs were stored in ambient temperature before testing. The ribs were placed in between input and output bar(s) with zero initial stress condition and the distance between these interfaces was always 120 mm so that we have uniform results.

Geometrical statistics of bone sample

A rib in general can be approximated to a curved beam. The most simple way to do this is to assume that the rib has a constant cross-section as shown in Fig. 2.8. Let us say that we have a set of data and we know is mean µ and its standard deviation σ, then the variance is σ 2 . Knowing these values, the probability density function pdf of normal distribution is expressed as

pdf = 1 √ 2πσ 2 e -(x-µ) 2 /2σ 2 (2.16)
This pdf signifies the relative likely hood when the value of a sample equals to that of the sample in a given space.

Experimental Results

The experimental results presented in this section are based on the waves shifted in time and space i.e. ε inc (t), ε ref (t), ε tr1 (t), and ε tr2 (t). We have two transmitted waves as the setup consists of two transmitted bars. If the stress distribution is below the yield limits of the bar material, the dynamic variables used for characterization of material / structure could be calculated as in Eq. 2.17a to Eq. 2.17f through the plane wave propagation.

F in = A b E b [ε inc + ε ref ] (2.17a) F out1 = A b E b [ε tr1 ] (2.17b) F out2 = A b E b [ε tr2 ] (2.17c) v in = -c b [ε inc -ε ref ] (2.17d) v out1 = -c b [ε tr1 ] (2.17e) v out2 = -c b [ε tr2 ] (2.17f)
Apart form the equations Eq. 2.17a to Eq. 2.17f, two supplementary equations were also used to characterize the response of rib, they are the relative interface velocity V and relative interface displacement U and are expressed as follows

V = v in - v out1 + v out2 2 (2.18) U = ˆt 0 V dt (2.19)
The tests carried out had permitted us to identify an ensemble of 20 treatable tests, with and without fracture of the specimen. With the help of air canon, we were capable of subjecting the specimens to relative velocities V between 5 m.s -1 and 13 m.s -1 . The ensemble of data had led us to the identification of transition phenomenon in the dynamic response of rib structure. The results are presented below.

Tests without fracture

The name itself says that no visible fracture was observed in this test. Screen shots of high speed camera are shown in figure below. (Note: The red mark indicates the initial position of the input bar-sample interface). This peculiar case represents the response of rib at V = 6 m.s -1 , v str = 9.27 m.s -1 .

Through Fig. 2.10, we can quantitatively say that the input bar imposes a displacement of 7.5 mm during 1.5 × 10 -3 s from the point time of impact. We do see negligible displacement of output bars. Even though Fig. 2.10 provides us with some quantitative information, the analysis of forces involved is important for a better comprehension. 2.11a. These forces permits us to identify and analyze the equilibrium state involved in this test. We were able to identify three different phases in the response of rib:

• Phase A, which is seen during the first 200 × 10 -6 s of the response. In this phase, we observe a drastic increase in the F in , followed by a sharp shoot. On the other hand, the ΣF out is zero during 150 × 10 -6 s. This proves that the waves generated at the point of impact in the sample, propagate in the sample and are yet to reach the output bar-sample interface. At 200 × 10 -6 s, the waves start to agitate the output bars, they start to exert reaction force on the sample.

• Phase B, this phase seems to exist between 200 × 10 -6 s and 600 × 10 -6 s. It is during this phase, we observe the reflection of waves from the output barsample interface towards the input bar-sample interface. This phenomenon is characterized by the oscillation of ΣF out around F in .

• Phase C, this phase seems to come into existence beyond 600 × 10 -6 s. The equilibrium is established between F in and ΣF out during this phase and the sample it subjected three-point bending.

Force vs. Displacement response

The F in vs. U aids in the comprehension of rib structural response. The U is calculated from Eq. 2.19.

Tests with Fracture During Dynamic Equilibrium

The Fig. 2.12, represents a case of fracture during the dynamic equilibrium phase, corresponding to an experimental case with V = 9 ms -1 and v str = 17.39 ms -1 . Similar to the response in the case 2.4.2, we can differentiate three phases of response. Firstly, a transition zone during first 100 × 10 -6 s, followed by a zone with increase in output forces between 100 × 10 -6 s and 600 × 10 -6 s. Beyond this point, the F in and ΣF out are almost equal, which helps us to identify the equilibrium is established at ≈ 625 × 10 -6 s. In addition to this, the Fig. 2.13a alongside Fig. 2.12, also helps us to identify that a macro crack stats to develop at 985 × 10 -6 s and leads to total rupture of the specimen at around 1100 × 10 -6 s from the time of impact.

(a) t = 0µs (b) t = 1090µs (c) t = 1450µs

Force vs. displacement response

The F in vs u rib for this case is presented in Fig. 2.13b. Based on the data, we can identify that the rib was subjected to a displacement of 9 mm, the debut of crack noted at u rib = 7.9 mm leading to rupture at u rib = 8.75 mm.

The analysis of all the experiments showing evidences of the presence of equilibrium phase, leads us to a conclusion that "the dynamic equilibrium is established between 600 × 10 -6 s and 800 × 10 -6 s". The time at which equilibrium is established, depends on properties of rib such as celerity of waves and geometry. In fact, in this case, we seen only the first two phases i.e. Phase A and Phase B. Moreover, this case does show evidence of a time difference between F out1 and F out2 , which is assumed to appear since the tested rib sample had some geometrical asymmetry.

The reason why, this experimental case is presented in this section is that we do not observe the Phase-C in force vs time response. A keen observation to Fig. 2.15a, shows that the rupture occurs at 520 × 10 -6 s, which is < 600 × 10 -6 s. Thereby proving that equilibrium is never achieved.

Force vs. displacement response

The F in vs U for this case is presented in Fig. 2.15b. Based on the data, we can identify that the rib was subjected to a displacement of 6 mm, the debut of crack noted at U ≈ 4 mm leading to rupture at U ≈ 5.6 mm.

Analysis of Experimental Results

Effect of striker velocity on time for fracture

The variation of time at fracture t f rac from the instance t ′ 0 is presented in Fig. 2.16. The o points signify cases with fracture during equilibrium whilst the stars represent fracture before equilibrium. We can notice that as the striker is launched at velocities greater than 15 [m.s -1 ], t f rac decreases to below 0.8 × 10 -3 [s]. Another interesting point to notice is that the experimental cases with fracture before equilibrium occur for t ≤ 0.8 × 10 -3 [s]. To be precise, out of 5 such experimental cases, 4 of them are around 0.5 × 10 -3 [s]. Thus, we can deduce that the characteristic time for fracture t f rac to have a fracture before the dynamic equilibrium is established is around 0.5 × 10 -3 [s].

Energy at fracture

Similar to t f rac , the energy at rupture W r as a function of relative interface displacement at rupture U r is presented in Fig. 2.17. We notice that around 75% of the experimental results are grouped in displacements between 0.005 [m] and 0.01 [m] whose corresponding energies at rupture vary between 2 [J] and 6 [J].

The black markers present the mean of each category.

Comparing W r with t f rac , we can notice this dispersion of energy with respect to fracture before and during equilibrium in a more significant manner (see Fig. All out together, we can deduce that higher the striker velocity, less is the time and energy required for fracture. Thus, we can establish a barrier for transition between fracture before equilibrium and fracture during equilibrium in terms of t f rac ≈ 0.5 × 10 -3 [s] and W r ≈ 4.2 [J]. If the sample is rigid enough to surpass this barrier, it shall be capable of establishing an equilibrium thereby producing a plateau kind response in the force response.

Fracture Mode in ribs subjected to dynamic impacts

The question posed quite often in three-point bending tests, is the mode of fracture observed in the sample whilst subjected to dynamic impact loads. In this case, we find three possible modes of fracture:

1. A rupture controlled by shear stress without global bending of the sample.

2. A rupture due to crack originating from the tensile surface due to global bending of the sample.

A combination of two preceding modes.

The ensemble of 20 treatable experimental tests carried out had permitted us to conclude that

• The rupture of the sample is due to a crack originating from the tensile surface of the sample.

• This so called rupture was observed to occur either during the dynamic equilibrium or before the equilibrium was established.

• The maximum force imposed on the sample by the input bar i.e. F in does not influence the rupture. Since, the maximum peak in F in is seen during the transition phase (Phase A) but neither in Phase B nor in Phase C, where the rupture is supposed to occur.

• Rather, the rupture is dependent on the local displacement of the rib. We had observed that the displacement attained by the sample during rupture was in between 4 mm and 12 mm, whose lengths varied in between 100 mm and 135 mm.

Characterization of Structural response

On the other hand, the experimental protocol presented in this thesis had provided us the useful data insights for understanding the structural response of porcine rib. The insights in this concern are given through the impedance adaptation between the bars and different samples tested, there by helping us to identify the parameters influencing the mechanical behavior of bone subjected to dynamic loads.

Since the geometry of rib plays an important role in its mechanical behavior, it is therefore impossible to characterize the structural response using fixed geometry of rib. To avoid these difficulties, we had decided to characterize the response in terms of F in vs. U curves and the corridors obtained from maxima and minima of these curves. The F in vs. U curves were analyzed with respect to rate of displacement V at the input bar -sample interface. This led us to identify three different corridors:

1 Corridor 1, 5 m.s -1 ≤ V ≤ 7 m.s -1 . (see Fig. 2.19a) 2 Corridor 2, 8 m.s -1 ≤ V ≤ 10 m.s -1 . (see Fig. 2.19b) 3 Corridor 3, 11 m.s -1 ≤ V ≤ 13 m.s -1 . (see Fig.2.19c) This methodology of presenting the results in form of experimental corridors, helps us to understand the influence of V on F in and the bending displacement of ribs. We were able to identify the following:

• On the peak of F in seen during Phase A (of sec.2.4.1 and sec.2.4.2). At low rates of displacement (i.e. V ≤ 7), the peak of F in varies between 400 N and 800 N. On the other hand, for high rates of displacement (i.e. V ≥ 8), the peak of F in varies between 800 N and 1400 N.

• As V increases, the fracture of rib tends to occur before dynamic equilibrium. The time at which the rupture occurs, depends on boundary conditions imposed on the specimen and also on the rate of loading. In other words, more the rate of displacement, less time for fracture to occur.

• Finally, we do remark an influence of V on the shape of F in vs. U. We can see that at low rates of displacement, the force tends to increase continuously till rupture (see Fig. 2.19a). Whereas, at relatively high rates of displacement, the force decreases from its peak, followed by a subtle increase and stays constant (plateau like response) as in Fig. 2.19b and Fig. 2.19c. This plateau like response corresponds either to plastification or damage of structure.

Figure 2.19 -Experimental corridors

Note that the experimental response case of fracture before dynamic equilibrium is an interesting case as the energies and displacements are high enough that the Section 2.5 bone fractures too early. Moreover, in reality, we cannot assume the establishment of equilibrium all the time. Therefore, we have decided to try to replicate two experimental cases in the numerical simulations at experimental scale namely case 1 and case 2. Case 1 represents to an experimental case with a striker velocity v str = 17 [m.s -1 ] and Case 2 represents to the experimental case with a striker velocity v str = 19.21 [m.s -1 ], that is an average of striker velocities corresponding to all experimental cases with fracture before dynamic equilibrium. The experimental force displacement responses in corridor 3 are presented below

Influence of geometry on structural response

In order to understand the effect of geometry on the type of response seen i.e.

corridor type, we had grouped the corresponding samples with respect to corridors (see Fig. 2.21). We do see that the sample corresponding to corridor 2 (C-2) are concentrated near the global mean of all samples. Whereas, the samples corresponding to that of corridor 1 and corridor 3 have means less than that of the global mean. For example, let us consider the property l of the samples. The global mean presented in 2.3 µ l is 113.59 mm, the means of lengths corresponding to corridors 1,2 and 3 are 109.68, 113.59 and 111.78 respectively.

A similar tendency is noticed among all other parameters. Unfortunately, we are unable to express the distinct nature of geometrical effects. In fact, to clearly understand the effects of geometry, it is advisable to carry out tests on a much larger sample size with considerable variations in their geometries.

Summary

The difficulties involved in a typical Split Hopkinson Pressure Bar system were over come by choosing nylon bars. Since these bars have a mechanical impedance similar to bone material, we were able to differentiate and shift the different waves of SHPB. Moreover, the sanding of porcine rib anterior and posterior extremities had ensured the establishment of decent contact with output bars.

These altogether, led to the comprehension of rib's response to dynamic three-point bending loads. It was found that the geometry of the specimen plays an important role similar to the human rib.

In the cases of fracture during dynamic equilibrium and fracture before dynamic equilibrium, the rupture of rib was due to a crack originating at the tensile surface of the rib. On the other hand, the maximum force imposed on the sample did not influence the rupture, since the rupture was observed during plastification phase of the response. Rather, the rupture was highly influenced by the displacement of rib, which varied between 4mm and 12mm. This variation in displacement was once again due the difference in geometry of the sample.

It was also noted that the velocity at the sample and bar(s) interfaces influenced rib displacement. Which in turn influenced the response of rib in terms of increase in peak force, time for fracture and shape of the response. As the interface displacement rate increased, 1. there was an increase in the magnitude of peak force.

2. the fracture tends to occur before the dynamic equilibrium is established.

3. the force tend to decrease from its peak, followed by a subtle increase and a plateau response.

These results show that the rib as a whole structure is sensible to strain rate.

Chapter 3

Generation of numerical rib sample

Résumé :

L'échantillon numérique a été généré à partir des images DICOM obtenus par imagerie tomographique à haute résolution (HR-pQCT). La géométrie 3D a été ensuite reconstituée à partir des contours des images DICOM en utilisant le module de numérisation IP du logiciel Simpleware.

Les propriétés mécaniques (module d'Young, contraintes limites,...) sont reliées à la masse volumique, laquelle est elle même corrélée au niveaux de gris (échelle HU). Comme il y a plusieurs définitions de la masse volumique, nous en avons testé deux : la masse volumique apparente ρ app et la masse volumique apparente de contenus de cendre ρ α . Dans un premier temps la densité est corrélée avec la variable HU des images tomographiques. Ensuite il est possible de connaître les propriétés mécaniques associées à chaque zone de niveaux de gris par l'intermédiaire de la relation ρ app (HU) ou ρ α (HU). Nous avons utilisé les relations entre les propriétés et la densité publiées dans la littérature.

Nous constatons que les propriétés mécaniques en fonction de ρ app sont les plus pertinentes dans le cadre de cette thèse.

Comme ρ app nous permet de connaître la fraction volumique et la porosité P o de chaque constituant de la côte numérique. Connaissant la porosité, nous pourrions utiliser les équations de loi de puissance pour connaitre les propriétés mécaniques de ces constituants. french

Abstract:

The numerical sample was generated from the DICOM images obtained through high resolution peripheral computed tomography (HR-pQCT). The 3D geometry and mesh were constructed from the contours of DICOM images using the Scan IP module of Simpleware.

The mechanical properties of bone constituents are linked to density values. We have tested two different mass density definitions in order to ensure which density definition is a reliable measure for mechanical properties for predicting the rib response observed experimentally. First of all, the apparent density ρ app was correlated with grey scale values expressed in HU through a linear equation. Secondly, the apparent ash density ρ α was tested. The mechanical properties of each rib constituent were then found through power law regression equations available in literature.

It is found that the mechanical properties from ρ app are the most reliable for this thesis. As ρ app , gives us access to corresponding volume fraction V f and porosity P o values, respective power law equations were used to differentiate between cortical and trabecular bone.

HR-pQCT Scan procedure

A batch of porcine ribs were bought from a local meat processing plant in order to ensure homogeneity of samples. These samples were wrapped in a damp cloth, soaked in saline (9% NaCl) to avoid dehydration. The damp samples were then stored at stored in refrigerator at 4 • C whilst not in use. They were defatted in an aqueous solution of Sodium Hypochlorite (NaClO) and Salt (NaCl) prior to the CT scan. Later, they were scanned using High Resolution peripheral Quantitative Computed Tomography (HR-pQCT)(XtremeCT,Scanco,Switzerland). The chosen to determine the heterogeneous mechanical properties of eight parts within the rib. A grey intensity (or level) is associated with a particular value of the linear attenuation coefficient µ. Five-sample calibration phantoms (Mindways R ) with known equivalent K 2 HPO 4 and H 2 O densities were placed below the rib during the scan. Using an algorithm implemented into Simpleware R code, we obtained the linear calibration coefficients µ to convert the raw data into normalized data (GS), expressed in Hounsfield Units (HU) using the following linear expression [START_REF] Taylor | Determination of orthotropic bone elastic constants using FEA and modal analysis[END_REF][START_REF] Teo | Relationship between ct intensity, micro-architecture and mechanical properties of porcine vertebral cancellous bone[END_REF]:

GS = 1000 (µ b -µ w ) (µ w -µ a ) ; µ w = 0; µ a = -1000. (3.1)
where: µ b , µ w , and µ a are the values for bone, water and air respectively. This was done by calculating the average GS within each of the five calibration samples along the length of the scan and calculating the slope and intercept vs the known equivalent K

Interpolation of Apparent Density

To quantify the mechanical properties of the constituents of porcine rib, the apparent density was first interpolated using the hypothesis put forth by [START_REF] Taylor | Determination of orthotropic bone elastic constants using FEA and modal analysis[END_REF]. The hypothesis is as follows 1. Having verified the units of CT values, first we calculate the effective density ρ ef f based on the relation proposed by [START_REF] Taylor | Determination of orthotropic bone elastic constants using FEA and modal analysis[END_REF]. According to them, ρ ef f includes the mass of bone, bone marrow and bone fluids. In other words it is equivalent to ρ wet

ρ ef f [kg.m -3 ] = 0.523GS + 1000 (3.2)
2. As we know that the dry density of water is zero and the apparent density of the hard bone segment in the rib corresponds to max(ρ ef f ). Using this hypothesis, the apparent density is correlated to HU with a line of slope a and zero intercept. This leads to the following relation A close look to the posterior end of the rib in Fig. 3.3 shall reveal that there is almost zero distribution of high density value elements at the posterior end. This is due to the defatting process carried out before CT scans. The corresponding region of the real porcine rib shall also be noted in Fig. 3.1. This shows that the defatting process does modify the structure a little bit. Unfortunately, this was the best and optimized result that we could obtain for this thesis.

Mechanical properties for numerical rib sample

As said in Sec. 1.3, two different regression models were tested i.e. one variable model (either ρ app or ρ α ) as control variable and two variables model with volume fraction and ash fraction (V f and α) as control variables.

One variable model

Three different models were tested for mechanical properties based on single variable power law regression. Out of which, one model considered ρ app as control variable (referred as Model-1a) and two models considered ρ α as control variable (referred as Model-1b and Model-1c)

Apparent density ρ app

Model-1a: As said above, this model considers only apparent density in Fig. 3.2. The constituents of porcine rib were then identified by porosity Po app from volume faction V f app through the Eq.

V f app = ρ app ρ tissue (3.4) Po app = 1 -V f app × 100 (3.5)
Based on the hypothesis proposed by [START_REF] Keaveny | Bone Mechanics[END_REF], cortical and trabecular bone constituents were identified as in Tab. 3.2. Based on the porosity values, corresponding power law regression equations were used. The equations for cortical and trabecular bone are given in tab. 3.3 respectively.

Apparent ash density ρ α

As said above, two different models were tested for ash fraction (α) and apparent ash density (ρ α ). Model 1b is based on the study proposed by [START_REF] Keller | Predicting the compressive mechanical behavior of bone[END_REF] whilst the ρ α in Model 1c was interpolated through Equivalent Mineral Density (ρ EQM ) from CT P R .

Model 1b:

The ρ α in this model was interpolated to ρ app of porcine rib data using the following relation proposed by [START_REF] Keller | Predicting the compressive mechanical behavior of bone[END_REF] ρ α = 0.522ρ app + 0.007; f or ρ app < 1000 kg.m -3

(3.6) Since the power law regression presented by [START_REF] Keller | Predicting the compressive mechanical behavior of bone[END_REF] is based on pooled data from human femur and vertebrae, the equations relating to the later were considered for interpolation of mechanical properties and are given below. In order to compare the two models (1b and 1c), the ash fractions (from Eq.3.13) and the mechanical properties through Eq. 3.8, Eq. 3.9 and Eq. 3.10 were considered. The range of α from Eq. 3.6 and Eq. 3.7 is 0.795 < α lit < 0.8 whilst the range of α from ρ EQM was 0.672 < α EQM < 0.676. According to [START_REF] Martin | Porosity and specific surface of bone[END_REF], the ash fraction of a fully mineralized bone is ≈ 0.7 and that of an osteon is 0.0.

ρ α = 0.
Based on this hypothesis and comparing the α values from model 1b and model 1c, suggest that α EQM is more reliable for porcine rib than that of α lit . It is clearly seen that the variation in strength of bone is more pronounce due to change in ash content rather than in apparent density or volume fraction. On the other hand, models based on ash content seem to under estimate the elastic property of bone.

Two variables model

In order to have a better understanding of the influence of volume fraction and ash fraction, [START_REF] Hernandez | The influence of bone volume fraction and ash fraction on bone strength and modulus[END_REF] had put forth two variables based models that were supposed to give a better insight. They had considered both V f and α as control variable for the pooled data produced by [START_REF] Keller | Predicting the compressive mechanical behavior of bone[END_REF]. The application of this model to numerical porcine rib is referred as model 2 in this thesis. Model 2: the tissue density ρ tissue is considered as a function of α as in Eq. 3.14, Section 3.3 which is based on the the data of [START_REF] Martin | Porosity and specific surface of bone[END_REF].

ρ tissue =
ρ app 1.41 + 1.29α (3.14) Since the range of α EQM is more relevant to data produced by [START_REF] Hernandez | The influence of bone volume fraction and ash fraction on bone strength and modulus[END_REF], it was considered as reference for two variables model. The interpolated V f and ρ tissue are presented in Fig. 3.9 It is clearly seen that the ρ tissue density determined from this hypothesis is surely over estimated and is much higher than that considered in Sec. 3.3.1. Thereby the V f is under estimated, thus suggesting that the porcine rib consists of only one cortical bone constituent in contrast to Tab. 3.2. Yet this proves that there is an influence of V f but is less comparing to that of α.

Comparison of mechanical properties

In order to make sure which hypothesis shall be considered for the material properties of bone, the results of single variable model and two variable model were compared and verified with the element distribution in the numerical rib sample. • The ρ α from model 1b is slightly higher than that determined from model 1c. This difference in ash density can be explained by the origin of sample. For example, the relations proposed by [START_REF] Keller | Predicting the compressive mechanical behavior of bone[END_REF] were based on pooled from human vertebrae and femur. Thus showing a difference in mineralization for different species.

• On the other hand, the ρ tissue density from model 2 is also higher than that used in model 1a. which can be explained due to the consideration of α in Eq. 3.14.

The ash fraction does not only influence density measures, it rather has an effect on almost all the mechanical and architectural properties. For example, let us consider the V f and Po depicted by blue and orange colored markers respectively in Fig. 3.12.

Regarding the mechanical properties, following conclusions can be drawn

• The young's modulus obtained thorough apparent density and two variable model are of similar range and are largely high to those obtained from apparent ash density based models (see Fig. 3.13a).

• The ultimate properties are highly influenced by α or ρ α as it can be seen in Fig. 3.13b. Even though the model 2 considers α, it does produce results similar to those obtained from apparent density. These data prove that either V f or ρ app based regression models are efficient and reliable for having a decent insight on mechanical properties of a numerical model.

Summary

The apparent density ρ app follows a linear relationship with grey scale intensity of porcine rib CT P R . Regarding the interpolation of mechanical properties, two one variable models and one two variable models were tested. Results show that the considering either apparent ash density ρ α or ash fraction α as independent variable, overestimates the ultimate and post yield properties. Therefore, mechanical 95 Section 3.4

properties based on ρ app are in good correlation and also are in the same range of properties used in previous numerical studies of human ribs (see Tab. 3.4). Yet this proves that porcine ribs are biologically reliable as human surrogates.

Chapter 4

Definition of numerical model

Résumé :

Dans ce chapitre, nous présentons les choix principaux concernant : la loi de comportement et la loi d'endommagement, le maillage, l'algorithme numérique du contact.

Les résultats expérimentaux nous montrent que la vitesse relative des interfaces

n'influence pas la partie linéaire de la réponse de la côte. Par contre, elle a une forte influence sur la partie plastique. Donc, une loi de comportement élasto-viscoplastique a été retenue. Trois formulations de lois de comportement élasto-viscoplastique ont été testées.

Une loi d'endommagement incrémentale et dépendante de l'état de contrainte a été couplée avec les lois de comportement testées. Cette loi d'endommagement est capable de tenir compte de l'accumulation non-linéaire de l'endommagement et ses conséquences (perte de rigidité...). De plus, l'influence de la vitesse de déformation ( ε) sur la déformation à la rupture (ε f ) a été prise en compte.

Un algorithme de contact entre les surfaces, basé sur la méthode des pénalités a été choisi pour les interfaces entre les barres et l'échantillon. Pour tenir compte de la cohésion entre les éléments internes de la côte (zones dotées de propriétés mécaniques différentes), un contact intérieur a été défini. Ainsi, le solveur numérique considère la côte comme une structure composite. french

Abstract:

This chapter presents the important choices considered for simulation of porcine rib at experimental case, such as the constitutive law, the damage law, mesh, and contact between the sample and bar interfaces.

It is well known that the bone is strain rate sensitive, thus use of constitutive laws such as viscoelastic-viscoplastic or elasticviscoplastic could be noted in the literature available today. The experimental results obtained in this thesis proves that the relative interface velocity does not affect the initial linear response of the sample, rather it has a pronounced effect on the post yield behavior. Therefore, an elasticviscoplastic constitutive law is the suitable law in this thesis. A total of three elasticviscoplastic constitutive laws were tested.

As all three constitutive laws do not offer coupled damage law, an incremental and stress-state dependent damage law was chosen. This damage law is dependent on stress-state of the sample with a capacity to consider non-linear damage accumulation and reduction in force carrying capacity with increase in damage. Moreover, the effects of strain rate ε on fracture strain ε f was taken into consideration through a logarithmic equation.

A penalty based surface-to-surface contact was chosen for contact at bar-sample interfaces. As the numerical sample is a distribution of different densities, thus different mechanical properties, a contact ("Contact Interior") was defined internally between each elements of the sample in order to represent the internal cohesion of the structure. By doing so, the solver is capable to apprehend that these constituents form a single structure.

The bibliographical research and experimental results have shown that the numerical sample must have a decent geometrical representation, the material model for bone must consider the effects of strain rate both in mechanical behavior and damage. Therefore, we present a numerical model generated through CT scan images of porcine rib for a good geometrical representation and an elasticviscoplastic constitutive law, coupled with an external damage law that takes the effects of strain rate on damage accumulation.

Apart from these considerations, the numerical model does involve one dimensional wave propagation in bars and contact between each constituent of the numerical model. The numerical constraints addressed are:

• Wave propagation in bars

• Constitutive and damage law for rib and its constituents

• Contact algorithm for bar and sample.

Wave propagation in a bar

First of all, the choice of constitutive law for bars is important. The bars of SHPB experimental setup are made of nylon in order to have mechanical impedance of same order of magnitude as bones. Even though nylon is a viscoelastic material, the viscosity of nylon induces only small variations to the local linear relations between particle velocity, stress and strain [START_REF] Zhao | On the use of a viscoelastic split hopkinson pressure bar[END_REF]. This led to the use of classic elastic constitutive law for nylon bars in our simulations. Another classical problem encountered in a numerical model with wave propagation in a non-uniform mesh is that it produces wave reflection and there is a difference in wave celerity. Thus, the accuracy of a numerical simulation depends on the mesh size and mesh density [START_REF] Roth | Chest injuries in children: An analysis of 100 cases of blunt chest trauma from motor vehicle accidents[END_REF]. Moreover, the mesh size and its variation does affect the uniform propagation of waves in a numerical simulation [START_REF] Bazant | Spurious reflection of elastic waves in nonuniform finite element grids[END_REF]. The main assumption considered in SHPB setup is one dimensional propagation of waves. Therefore, in order to avoid spurious wave reflection, we had Section 4.1 adapted the principle put forth by [START_REF] Roth | Chest injuries in children: An analysis of 100 cases of blunt chest trauma from motor vehicle accidents[END_REF]. This relation is used to determine the mesh density (Γ ) and mesh size (∆x) for a desired wave length (λ) so that the celerity of numerical wave (c F EM ) is almost equal to theoretical celerity of waves (c th ) as in Eq. 2.5.

According to [START_REF] Bazant | Spurious reflection of elastic waves in nonuniform finite element grids[END_REF], the relation between c F EM and c th is given as

c F EM c th = 2 Φ (sin Φ 2 ) -2 - 2 3 m -1 2 (4.1)
where, Φ = 2π∆x

λ (4.2) Γ = λ ∆x (4.3) ∴ Φ = 2π Γ (4.4)
where, m in Eq. 4.1 represents type of mass distribution m = 0 for lumped mass matrix, m = 1.0 for consistent mass matrix and m = 1.5 for hybrid mass matrix. The variation of celerity ratio for respective mass matrices is presented in As the bar is considered a deformable body and deformable solid elements in LS-Dyna are made of consistent mass matrix i.e. 'the mass of a element is distributed evenly at all the eight nodes of the element'. Therefore, the variation of CR , resembles to that of Fig. 4.1b.

In order to test the desired mesh size for a uniform wave propagation in bars, a samll scaled model was used as reference. This model consists a striker 0.1 m long and 3 × 10 -3 m in diameter and a bar 0.6 m long and 4 × 10 -3 m in diameter. Three different values of mesh density were tested for a specific wavelength (Γ = 12, 16 and 20). Three elements at l The change in mesh size did show a negligible affect on the strain signal response (see Fig. 4.2). On the other hand, the contact force was not affected. Even though Γ = 20 takes more time per cycle of the calculation, it is more accurate in terms of wave celerity. Thus, the corresponding mesh size was considered for bars. Another common problem encountered in numerical simulation of wave propagation is the element formulation with respect to number of Gaussian integral points. Previous studies that investigated the model convergence recommended using an element of size less than 3 mm to accurately capture the heterogeneous variations in the mechanical properties [START_REF] Keyak | Correlations between orthogonal mechanical properties and density of trabecular bone: use of different densitometric measures[END_REF][START_REF] Bessho | Prediction of strength and strain of the proximal femur by a CT-based finite element method[END_REF] Similar to cortical bone tissue, several constitutive laws were proposed for trabecular bone. The main difference is the architecture of trabecular bone, which makes it anisotropic. Moreover, the high porosity of trabecular bone leads to a foam like behavior. Some constitutive laws considering the trabecular bone architecture can be seen in [START_REF] Cowin | Fabric dependence of an anisotropic strength criterion[END_REF], [START_REF] Pietruszczak | A fabric-dependent fracture criterion for bone[END_REF] and so on. On the other hand, [START_REF] Gibson | Biomechanics of cellular solids[END_REF], [START_REF] Karkar | Personnalisation morpho-mécanique de la voûte crânienne humaine à différentes vitesses de sollicitations[END_REF] have considered foam like models for trabecular bone.

′ b /4, l ′ b /2
For simulations concerning the fracture of bone tissue, damage accumulation and evolution plays and important role. Numerous studies have pointed out the limitations of basic fracture criteria models such as von-Mises equivalent stress criteria, Hill's criteria to model the damage of bone. For example [START_REF] Malik | Equine cortical bone exhibits rising R-curve fracture mechanics[END_REF], [START_REF] Taylor | Microdamage and mechanical behaviour: predicting failure and remodelling in compact bone[END_REF], [START_REF] Vashishth | Experimental validation of a microcracking-based toughening mechanism for cortical bone[END_REF], [START_REF] Ural | Effects of intracortical porosity on fracture toughness in aging human bone: a microCT-based cohesive finite element study[END_REF] and Abdel-Wahab and Silberschmidt (2011) had used Fracture mechanics approach but failed to predict the complete fracture patterns since these models are restricted to single dominant idealized planar crack. Whilst, Hambli (2011a), [START_REF] Hambli | Multiscale prediction of crack density and crack length accumulation in trabecular bone based on neural networks and finite element simulation[END_REF] and [START_REF] Dall'ara | A nonlinear QCT-based finite element model validation study for the human femur tested in two configurations in vitro[END_REF] have proposed damage models based on continuum damage mechanics (CDM).

Note that the porcine rib model presented in this thesis is a homogenized model i.e. it does not consider the architecture of trabecular bone. Therefore, considering of foam like constitutive law for trabecular bone components may not be the priority here. Thus, we concentrate on considering rate effects in constitutive law and coupling it to a credible damage law.

Before choosing a reliable constitutive law, let us take a keen look at the experimental corridors presented in Fig. 2.19a,Fig. 2.19b and Fig. 2.19c. We notice that the slope of initial elastic response stays almost unaffected by the change in strain rate ε. This lays the foundation for an elastic-viscoplastic behavior of rib. The difference in behavior of cortical and trabecular bone is considered through hardening. Cortical bone was assumed to show isotropic hardening behavior in contrast to trabecular bone, which shows kinematic hardening behavior.

On the other hand, the experimental results have proved that the crack initiation occurred in the tensile surface and propagated towards compressive surface in the vicinity of impact load. Moreover, we could see some material erosion at the point of contact with input bar. Therefore, a credible damage law must consider this phenomenon.

Constitutive laws in LS-Dyna

As the experimental results prove that an elastic-viscoplastic constitutive law is suitable for the porcine rib, we must now choose a good candidate based on the data in our possession. In this case, we can narrow down to three possible candidates that are capable of considering elastic-viscoplastic behavior: This material model, offers elastic-viscoplastic behavior through Cowper-Symmonds model (i.e. post yield stress scaling) and is internally coupled with a damage law based on the hypothesis of effective stress proposed by [START_REF] Lemaitre | A course on damage mechanics[END_REF]. The yield function (φ)

φ = 1 2 S ij S ij - σ 2 y 3 ≤ 0 (4.5)
The yield function in Eq. 4.5, is expressed in terms of Deviatoric stress components (S ij ) and the dynamic yield stress (σ y ) . The effects of ε on post yield stress are considered through the Eq. 4.6.

σ y = 1 + ε C 1 p σ y 0 + f h (ε P ef f ) (4.6)
Here, C and p are Cowper-Symmonds parameters, σ y 0 is the initial yield stress and

f h (ε P ef f ) is the hardening function .
The hardening function is treated in two different ways:

1. The user can give different curves in a tabular form for f h (ε P ef f ) and LS-Dyna interpolates the hardening with respect to the curves provided in the table.

A default linear treatment, that is expressed as

f h (ε P ef f ) = E P (ε P ef f ) (4.7)
where, E P is the plastic hardening modulus, expressed in terms of tangent modulus (E T AN ) and Young's modulus (E) as in Eq. 4.8 and ε P ef f is the 104 Chapter 4

effective plastic strain, as in Eq. 4.9.

E P = E T AN E E T AN + E (4.8) ε P ef f = ˆt 0 2 3 εP ij ε ij P 1/2 dt (4.9)
The damage law is a continuum damage mechanics model, that considers effective stress concept. The effective stress for a damaging material is expressed as

σ = σ (1 -D) (4.10)
here, D is the damage variable. The evolution of D with time is considered as below

Ḋ =      Υ ς(1-D) ṙ for r > r D and σ I > 0 0 otherwise (4.11)
where, Υ is the strain energy density release rate, ς is the strain energy release rate, r is the damage effective plastic strain, r D is the damage threshold and σ I is the maximal principal stress. The strain energy density release rate in Eq. 4.11 is expressed in terms of triaxiality function R v (η) as below 2 (4.12)

Υ = σ 2 vm R v (η) 2E(1 -D)
R v = 2 3 (1 + ν) + 3(1 -2ν)η 2 (4.13)
where η is the stress triaxiality expressed in terms of hydrostatic stress (1/3tr σ)

and von-Mises stress (σ vm ) as

η = 1 3 tr σ σ vm (4.14)
The rate of damage accumulated plastic strain ṙ is expressed in terms of effective plastic strain rate ( εP ef f ) and damage (D) as in Eq. 4.11 is represented as

ṙ = εP ef f (1 -D) (4.15)

Modified elastic-viscoplastic material model (*MAT-003)

This material model has the ability to consider either kinematic hardening or isotropic hardening or even a combination of both kinematic and isotropic hardening. This can be achieved by varying a parameter β h between 0.0 and 1.0. If β h = 0.0, the material is considered as kinematic hardening and for β h = 1.0, the material is isotropic hardening. As we know that the center of the yield surface shifts according to the behavior for kinematic hardening, the yield condition can be expressed as

φ = 1 2 ξ ij ξ ij - σ 2 y 3 = 0 (4.16
)

ξ ij = S ij -α ij (4.17) σ y = 1 + ε C 1 p σ y 0 + β h E P ε p ef f (4.18)
S ij and α ij of Eq. 4.17 are deviatoric stress and the center of yield surface respectively. The variable σ y of Eq. 4.18 is the rate dependent yield stress, σ y 0 is the initial yield stress defined in the material parameters and the variables C and p are Cowper-Symmonds parameters. The other variables of Eq. 4.18 are calculated as follows:

ε = εij εij (4.19) E P = E T AN E E T AN + E (4.20) ε P ef f = ˆt 0 2 3 εP ij ε ij P 1/2 dt (4.21)
The plastic strain rate εP ij in Eq. 4.21 is difference between total strain rate εij and elastic strain rate εe

ij εP ij = εij -εe ij (4.22)
Numerical implementation of *MAT-003:

• First the stress are updated elastically through the Eq. 4.23

σ * ij = σ n ij + C ijkl ∆ε kl (4.23)
here, σ * ij is the trial stress tensor, σ n ij is the stress tensor from previous iteration, C ijkl is the elastic tangent modulus matrix and ∆ε kl is the increment in strain.

• The yield function is checked.

(a) If the yield function of Eq. 4.16 is satisfied, the solution advances to next iteration.

(b) If the yield function is violated, an increment in plastic strain is calcu-lated, the stresses are scaled back to yield surface and the yield surface center is updated.

Let us assume the trial deviatoric stress S * ij and it's trial invariant ξ * ij as of the state n + 1 as in Eq. 4.24 and Eq. 4.25 respectively

S * ij = σ * ij - 1 2 σ * kk (4.24) ξ * ij = S * ij -α ij (4.25)
Then the yield function is defined as follows

φ = 3 2 ξ * ij ξ * ij -σ 2 y = Λ 2 -σ 2 y    ≥ 0 for elastic response < 0 for plastic hardening (4.26)
where,

Λ = 3 2 ξ * ij ξ * ij
In the case of plastic hardening (i.e. for φ < 0), the following steps are executed

• updating of plastic strains

ε P n+1 ef f = ε P n ef f + Λ -σ y 3G + E P (4.27)
• scaling back of stress deviators

σ n+1 ij = σ n ij - 3G∆ε P ef f Λ ξ * ij (4.28) here, G = E/[2(1 + ν)]
• updating the center of yield surface

α n+1 ij = α n ij + (1 -β h ) E P ∆ε P ef f Λ ξ * ij (4.29)

Piecewise Linear Plasticity model (*MAT-024)

This model is a basic elastic-plastic constitutive law, that gives the user to choose yield stress scaling with respect to strain rate ε, thereby producing an elastic-viscoplastic law. The yield function in this model is

φ = 1 2 S ij S ij - σ 2 y 3 ≤ 0 (4.30) σ y = 1 + ε C 1 p σ y 0 + f h (ε p ef f ) (4.31)
The yield function in Eq. 4.30 is expressed in terms of deviatoric stress tensor S ij and dynamic yield σ y as opposed to that in Eq. 4.16, which considers the invariant of deviatoric stress component ξ ij . Moreover, the difference between *MAT-003 and this model is that it considers isotropic hardening only and the hardening function

f h (ε p ef f
) is treated differently. The *MAT-024 offers two different ways of treatment for the hardening function:

1. The user can give different curves in a tabular form for f h (ε P ef f ) and LS-Dyna interpolates the hardening with respect to the curves provided in the table.

A default linear treatment, that is expressed as

f h (ε P ef f ) = E P (ε p ef f ) (4.32)
where, E P and ε P ef f are expressed in similar manner as in Eq. 4.20 and Eq. 4.21 respectively. This material model has another advantage over *MAT-003. The rate effects can be considered either through a user defined table consisting the stress-strain curves of the material for different strain rates or by the for the variable VP. The later is considered in LS-Dyna as follows:

VP =            Cowper -Symmonds if VP = -1.0 Yield stress scaling if VP = 0.0 Entirely Viscoplastic if VP = 1.0
Therefore, this material model has an upper hand to *MAT-003 when it comes to consideration of strain rate effects. Numerical implementation of *MAT-024:

• the stress are updated in similar fashion to *MAT-003.

• the yield function in Eq. 4.30 is checked.

• If satisfied, the deviatoric stresses are accepted.

• if the yield function is not satisfied, an increment in plastic strain is computed through

∆ε P ef f = 3 2 S * ij S * ij 1/2 -σ y 3G + E P (4.33)
then the trial deviatoric stresses S * ij are scaled back to the yield surface

S n+1 ij = σ y 3 2 S * ij S * ij 1/2 S * ij (4.34)
Summary of feature available in all three constitutive laws is presented in Tab. 4.3.

Model yield function φ strain rate effects σ y ( ε) Damage Law MAT-105 1 2 S ij S ij - σ 2 y 3 ≤ 0 1 + ε C 1 p σ y 0 + f h (ε P ef f ) coupled (CDM) MAT-003 ξ ij ξ ij - σ 2 y 3 = 0 1 + ε C 1 p σ y 0 + β h E P ε P ef f none MAT-024 1 2 S ij S ij - σ 2 y 3 ≤ 0 1 + ε C 1 p σ y 0 + f h (ε P ef f ) none
Table 4.3 -Key factors of elastic-vicoplastic models in LS-Dyna Note that S in Tab. 4.3, denotes the devitoric stress components whereas ξ denotes the invariant of deviatoric stress components and in expressed as in terms of deviatoric stress components and center of the hardening curves. The material models MAT-105 and MAT-024 are the same in terms of numerical formulations of the constitutive law, whereas the MAT-003 approaches the problem in a different manner as it includes an additional parameter β h for taking the effects of kinematic hardening into consideration.

In contrast to similarities in constitutive behavior of these models, the material models MAT-003 and MAT-024 need an external damage law for an effective simulation of fracture and damage propagation. This makes them less efficient in terms of calculation time.

Understanding the constitutive laws

Elastic-viscoplastic formulation

Since the material model MAT-024 the only providing two different types of yield stress scaling for considering the effects of strain rate on plastic part, simulations were carried out on single element model as in Fig. 4.4a.

To understand the difference between VP = 0 and VP = -1.0, a force boundary condition at different loading rates (as in Fig. 4.4b) was imposed. The parameters of Cowper-Symmonds model C = 2.5 and P = 7.0 were adapted from [START_REF] Li | Rib fractures under anterior-posterior dynamic loads: Experimental and finite-element study[END_REF]. The stress strain response corresponding to each loading rate for VP = 0 (solid curves) and VP = -1.0 (dashed curves) is presented in Fig. 4.4c.

It is clear that the MAT-024 model with VP = -1.0 is the one more effective and does consider the effects of strain rate as it is supposed to. The yield stress scaling behavior of MAT-003 is also executed in a similar fashion. Moreover, as MAT-105 is just an extension of MAT-024 for constitutive law, it does behave similar to that of MAT-024 and MAT-003. Therefore, we can say that the Cowper-Symmonds model is the suitable candidate.

Consideration of hardening

The parameter β h is the one responsible for hardening in MAT-003. In order to understand its influence, a random cyclic displacement imposed on 4 nodes of one extremity whilst the other 4 nodes were constrained with 0 displacement in all directions as in Fig. 4.5. The response of material for the displacement curve is shown in Fig. 4.6 and Fig. 4.7. The effect of β h is shown in Fig. 4.6. We can notice that for β h < 0.5, the material response is similar to that of kinematic hardening. Whereas for β h > 0.5, it tends to be more isotropic.

The Fig. 4.7 shows the difference between the material response for kinematic hardening (β h = 0.0 in red) and isotropic hardening (β h = 1.0 in blue). We can notice that the material response depicted by red curve shows a shift in the stress strain curve, yet stress of similar order of magnitude to that of initial loading. This is characterized by a shift in α ij of the yield surface and almost constant radius. On the other hand, the response for isotropic hardening, depicted by blue curve shows an increase in magnitude and a negligible shift in α ij .

To understand, the difference between the formulation of isotropic hardening of *MAT-024 and that of *MAT-003, the same one element model was carried out but this time with *MAT-024 for the same load as in Fig. 4

External damage law for bone tissue

As said above, a credible damage law must be able to differentiate the tensile and compression surfaces in the rib sample and also consider the effects of strain rate in the reduction of strain at fracture.

Even though material models *MAT-003 and *MAT-024 provide failure option with strain at failure though element erosion, they are not sufficient enough to consider the progressive nature of fracture and strain rate effects, thus they must be coupled with an external damage law. In this case, the only external damage law satisfying our needs is *MAT-Add-Erosion of LS Dyna material library. This law is based on element deletion / erosion criteria. The *MAT-Add-Erosion, offers three different approaches for damage evolution: Erosion is a general approach, that considers the failure strain and elements are removed through erosion once this criteria is satisfied. On the other hand GISSMO and DIEM are more sophisticated damage models that work on the basis of element deletion depending on the criteria chosen by the user. For example, GISSMO is a phenomenological model that considers incremental damage accumulation, softening and failure as opposed to DIEM which relies on criteria given by user for damage initiation and evolution. A total of 5 criteria can be chosen for DIEM model. As we know that the stiffness and rigidity of bone decreases with increase in damage and the objective of this thesis is to capture progressive fracture process, we had chosen GISSMO as a potential candidate for an external damage law.

*MAT-Add-Erosion with GISSMO

The GISSMO model offers an incremental and stress-state dependent criteria, moreover it considers progressive evolution of damage, reduction in stiffness and strain rate effects as well. It was proposed by [START_REF] Neukamm | Considering damage history in crashworthiness simulations[END_REF] and [START_REF] Haufe | A Comparison of recent Damage and Failure Models for Steel Materials in Crashworthiness Application in LS-DYNA[END_REF]. The stress-state of the material is considered in terms of triaxiality η as follows

η = 1 3 tr σ σ vm (4.35)
here, 1 3 tr σ represents the hydrostatic stress or mean stresses and σ vm represents the von-Mises / equivalent stress.

Damage Evolution

The damage evolution is in fact a generalization of Johnson-Cook model. The path dependent criteria is taken into consideration through a user defined curve of static fracture strain as a function of ε f (η). Moreover, it also offers non linear damage evolution through the variable n. The damage evolution is expressed as follows

∆D = n ε f (η) D (1 -1 n ) ∆ε p (4.36)

Post Critical Behavior

Post critical behavior signifies the effects of damage on the behavior of a material such as softening, reduction in stiffness and rigidity until entire rupture. There are three different ways to invoke the post critical behavior, which is controlled by ECRIT variable in the LS-Dyna Add Erosion card.

1. By giving a constant value for ε crit . Once ε p reaches this value, the current damage value is noted as critical damage value i.e. D | (εp = ε crit ) = D crit . This point is considered as the onset of reduction in stress carrying capacity and the damage is coupled to stress through a modified equation of Lemaitre's effective stress concept:

σ = σ 1 - D -D crit 1 -D crit m (4.37)
where σ is the effective stress (as in Lemaitre concept), m is the fading exponent. The Eq. 4.37 continues to be couple until the D value reaches 1 and the element is removed from the calculation.

2. Second way is to provide a curve for ε crit with respect to η. In this case, an incremental instability measure is calculated in the similar manner to ∆D as in Eq. 4.36.

∆F = n ε crit (η) F ( 1 -1 n ) ∆ε p (4.38)
the variable F represents the instability and is supposed to vary as 0 ≤ F ≤ 1.

Once ∆F reaches unity, the damage is coupled to stress with Eq. 4.37.

3. By declaring ECRIT = 0 and giving D crit values. In this case, coupling of stress with damage starts as soon as D reaches D crit . This method gives us flexibility to enter any arbitrary value. If the value for D crit is unknown and the value of the fading exponent m is unity, the Eq. 4.37 is rewritten as follows

σ = σ (1 -D) (4.39)
The Eq. 4.39 is the same as that proposed by Lemaitre, found in MAT-105

In our case, the third method was chosen with values of D crit taken from literature and the fading exponent is unity.

Effects of strain rate

Since we deal with ribs subjected to dynamic impacts, literature proves that ε effects the damage behavior of bones [START_REF] Mcelhaney | Dynamic response of bone and muscle tissue[END_REF], [START_REF] Wood | Dynamic response of human cranial bone[END_REF] and [START_REF] Johnson | A viscoelastic, viscoplastic model of cortical bone valid at low and high strain rates[END_REF]. Thus, we had provided *MAT-Add-EROSION with an additional curve for the effects of ε, which is controlled by the variable LCSRS. The curves were based on the equation given by [START_REF] Wood | Dynamic response of human cranial bone[END_REF] and the static ε f values were changed in order to study its effects on the response of rib. The ratio of dynamic strain dependent fracture strain ε f d to that of static fracture strain ε f s varies with respect to strain rate ε as follows

ε f d /ε f s = A -BLn ε (4.40)
The difference in damage law considered in

Working principle of damage law

The damage law is incremental, stress state dependent and is capable of considering non-linear damage accumulation, and reduction in force carrying capacity with increase in strain rate. In order to have coherent results, this damage model needs a curve for considering the variation of fracture strain ε f with respect to stress state η, which is presented in Fig. 4.9a.

The lack of data for the variation of ε crit with respect to triaxiality η, forced us to use a critical damage value D crit was specified in the damage law. When the yield specified in the material card of the sample is reached, the damage D starts to couple with the stress field. Meanwhile, the fading exponent data is also unavailable for porcine rib, thus forcing us to set it to unity. The parameters required for post yield behavior such as critical damage D crit , damage exponent for non linear damage n and fading exponent m for coupling with effective stress were adapted from literature and are presented in Tab. 4.5. m/m, which corresponds to a quasi-static strain-rate of 3.18 × 10 -3 s -1 . Similarly, in the present study Eq. 4.40 is not applicable for ε lower than 3.00 × 10 -3 s -1 , which are quasi static. For all values 0 s -1 ≤ ε ≤ 3.00 × 10 -3 s -1 , a linear approximation is applied so that

Property

ε f d /ε f s = 1 at ε = 0 s -1 .
To understand the effects of damage law parameters on the behavior of sample, a simple three-point bending tests was carried out using rigid walls and a rectangular sample of 80 mm long and 100 mm 2 square cross-section Fig. 4.10a. The rigid walls were provided with velocities from one of the experimental cases (see Fig. 4.10b).

The velocities are negative since they are in the opposite Y direction. A total of three simulations were carried out, (i) with no *MAT-Add-Erosion i.e. the failure criteria was provided as a constant value in the material model card

ε f s = 2.0×10 -2 [m/m],
(ii) with *MAT-Add-Erosion but no effects of strain rate i.e. ε f (η), and (iii) with *MAT-Add-Erosion with effects of strain rate ε f d (η, ε). The input force F in is obtained from the contact algorithm between input rigid Section 4.3 wall and the sample, where as the relative displacement U is calculated through the formula in Eq. 2.19. The impulse was calculated through the integral of F in with respect to time. The energy data presented was obtained through part energy history in LS Dyna.

We can notice that the sample without damage law (represented by solid lines in Fig. 4.11) is subjected to more displacement and more energy before entire rupture. This is due to the absence of coupling between the damage D and flow stresses. On the other hand, the sample under the influence of damage law (ε f (η) in dash-dot line and ε f d (η, ε) in dashed line) do eventually rupture way before their counterpart. For example, let us consider the energy absorbed by the sample before rupture (Fig. 4.11c). The curves corresponding to the cases ε f η and ε f d (η, ε) have values much less than that of ε f s . This could be explained by the coupling of damage with flow stresses as soon as the stresses in the sample exceed its static yield limit σ yT . Another peculiar response is that of the sample with ε f d (η, ε) passes into plateau region for an energy of ≈ 2.8[J]. This is due to the reduction in fracture strain with increasing strain rate.

Numerical Model Considerations

The numerical considerations such as decisions for which constitutive law to use, coupling of constitutive law with an external damage law if required and contact parameters were done using a single constituent porcine rib model (Fig. The difference between the real porcine rib sample and single part sample is that, the later is made of only one bone constituent instead of 8.

The rigid walls were modeled using the *MAT-020-Rigid of LS-Dyna material library. This material model furnishes the ability to convert any geometrical mesh of either shell or solid elements to convert into a rigid body. As the elements of "rigid body" are bypassed to store history variables, this model is cost effective. In order to replicate the sample behavior as close as possible to the original SHPB model, input and output interface velocities of one of the experimental cases were given to respective input and output rigid walls (see Fig. 4.10b). Since thee definition of contact between porcine rib sample and rigid walls is mandatory, the rigid wall were attributed with same properties of nylon as in Tab. 4.1.

To finalize the constitutive law, damage law and contact algorithm, all three constitutive laws (Sec. 4. Some researchers have studied the type of fracture in bones and have found that there is a ductile to brittle transition in bones [START_REF] Hasson | A ductile-to-brittle transition in bone[END_REF] and in case of dynamics, [START_REF] Kirchner | Ductility and brittleness of bone[END_REF] reported that this transition (not sharp) from ductile to brittle but occurs when the strain rate is around 10 s -1 . Thus the fracture we are supposed to see in our simulations is either brittle fracture or quasi-brittle fracture. This lead to the decision to approximate the damage law parameters for cortical bone as in Tab. 4.9.

NOTE:

The damage parameter ς is an approximation based on that for concrete. Concrete was chosen since, it shows quasi-brittle fracture.

The simulation seemed to be credible at initial stages till the penetration was encountered between the elements of input rigid wall (Master segment) and single part rib sample (Slave segment). It seemed that the penetration was due to parameters chosen in contact algorithm. The problem still persisted upon changing the contact parameters. Upon further investigation, we noticed that the elements were eroded and removed from the calculation at the compression surface of rib. Even though this material model is complex and efficient, it is unable to differentiate between the compression and tension surfaces when applied to three point bending. Moreover, the erosion of element are not only due to this reason, yet they might occur due to the choice of damage parameters. Nevertheless we have no choice to let go off this constitutive law.

Single Part Model-24

A total of 3 different models were carried out i.e. Model-24a [START_REF] Pavier | Contribution à la compréhension des phénomènes physiques lors de l'impact d'un corps sur un modèle de structure biologique[END_REF], Model-24b [START_REF] Pavier | Contribution à la compréhension des phénomènes physiques lors de l'impact d'un corps sur un modèle de structure biologique[END_REF] plus VP = -1.0 and Model-24c [START_REF] Li | Rib fractures under anterior-posterior dynamic loads: Experimental and finite-element study[END_REF]

Model-24a Results

Crack Propagation (Model-24a)

On applying the phenomena specified of dynamic crack propagation [START_REF] Miannay | Time-Dependent Fracture Mechanics[END_REF] to the results of this simulation, we see that the damage is still at the microscopic level till 0.46 ms. At 0.48 ms, we have a crack that surfaces at macroscopic scale and is initiated in the sample. At 0.8 ms, we see the arrest of the crack. Apart from this, we also notice along-branching or kincking. (see figure: The forces were recovered at the contact interface of Rigid walls and Sample.

Initially, we see a great amount of increase in Input Contact force till it reaches its maximum value of 4.15kN and starts to decrease. The decrease in contact force signifies the increase in microscopic damage and reduction in force carrying capacity between t = 0.132 ms and t = 0.46 ms. Once the microscopic voids/cracks surface at macroscopic level, we see the crack initiation at t = 0.48 ms and propagation between t = 0.48 ms and t = 0.80 ms. After which, the sample breaks entirely. (see figure:4.15)

Model-24b Results

In Model-24b, the numerical rib was give the same properties an in table:4.10. In addition to that, this model also take into account the strain rate effects through Cowper-Symmonds model i.e. VP = -1.0. The Cowper-Symmonds parameters were adapted from [START_REF] Li | Rib fractures under anterior-posterior dynamic loads: Experimental and finite-element study[END_REF], C = 2.5 and P = 7.0.

Crack Propagation (Model-24b)

Whereas, on including the strain rate formulation through Cowper-Symmonds model, we see that the crack initiation is delayed by 0.04 ms, which is not a significant difference. The crack is initiated at 0.5 ms and at 0.8 ms we see the arrest of the crack. In contrast to model-2a, we do not see knicking. Initially, we see a great amount of increase in Input Contact force till it reaches its maximum value of 4.95kN and starts to decrease. The decrease in contact force signifies the increase in microscopic damage and reduction in force carrying capacity between t = 0.132 ms and t = 0.50 ms. Once the microscopic voids/cracks surface at macroscopic level, we see the crack initiation at t = 0.50 ms and propagation between t = 0.50 ms andt = 0.80 ms. After which, the sample breaks entirely.(see figure:4.17

Model-24c Results

Model-24c, is based on the properties adapted from [START_REF] Li | Rib fractures under anterior-posterior dynamic loads: Experimental and finite-element study[END_REF] and it does consider the strain rate effects through Cowper-Symmonds model with same values of C and P. In contrast to [START_REF] Li | Rib fractures under anterior-posterior dynamic loads: Experimental and finite-element study[END_REF], we have also added a damage law based on GISSMO model with triaxiality curve.

Crack Propagation (Model-24c)

On applying the phenomena specified above to the results of this simulation, we see that the crack is initiated at time between t = 0.32 ms and t = 0.36 ms. At t = 0.36 ms, we have a crack that is already initiated and propagating in the sample.

At t = 0.40 ms, we see the arrest of the crack. Apart from this, we also notice that there is no along-branching or kincking. (see figure:4.18). Once the crack comes to arrest (in this case the sample is completely broken), the fragments start rotating.

Apparently, the Input rigid wall re-establishes contact and we see Contact force signals again. Initially we see increase in contact force till t = 0.1179 ms, at this point Microscopic damage is in progress till we have a crack that surfaces at macroscopic level (t = 0.36 ms). This crack then propagates till the sample is entirely broken t = 0.545 ms. But we do not see any kinking. Alongside, we do notice that the sum of output forces is greater than Input force for most of the simulation i.e. the system is not in equilibrium. The contact force at the Input wall-sample interface does not drop down to zero as the wall reestablishes contact with the broken sample and the elements stay in contact.

Model-3a Results

Crack Propagation in Model-3a

On applying the phenomena specified above to the results of this simulation, we see that the crack is initiated at time between 0.32 ms and 0.34 ms. At 0.34 ms, we have a crack that is already initiated and propagating in the sample. At 0.4 ms, we see the arrest of the crack. Apart from this, we also notice that there is no along-branching or kincking. (see figure: Moreover, we also notice that the sum of contact forces between Output wall and the sample is greater than force between input wall and sample, this signifies that the system is not in equilibrium. This phenomenon coincides with the experimental results of [START_REF] Pavier | Contribution à la compréhension des phénomènes physiques lors de l'impact d'un corps sur un modèle de structure biologique[END_REF] as the experimental velocity were taken from a case that was not in equilibrium.

Crack Propagation (Model-3b)

In contrast to model-24b (*MAT-024), this model behaves a bit differently. We see a huge delay in crack initiation i.e. instead of initiation at t = 0.5 ms, the crack in model-3b is initiated in between t = 0.66 ms and t = 0.68 ms. Moreover, two cracks are initiated instead of one but only one crack propagates from t = 0.68 ms and the crack is arrested at t = 0.8 ms.

Force Analysis (Model-3b)

Like all other models, initially the input force increases till t = 0.1635 ms but has the same magnitude as in Model-24b. The decrease in input force is backed-up by increasing microdamage thereby a decrease in Force Carrying capacity. The Input force does not descend to zero as the Input wall always see some elements of the rib.

Contact Algorithm

An automatic surface to surface contact was used for contact between bars and porcine rib. This contact algorithm uses a segment based penalty approach. On the other hand, an interior contact was defined between elements of each constituents of porcine rib in order to assist transmission of forces and to avoid negative elements.

The segment based automatic surface to surface contact is not dependent on segment orientation i.e. the algorithm looks on both salve and master segments, thereby producing a symmetric contact. Thanks to the algorithm, we can obtain the forces involved in the contact. In fact the contact force depends on the shape of the geometry. The contact algorithms and their differences are presented in AppendixB For example, let us consider porcine rib in three point bending (to your left in The two geometrical configurations were modeled using the same contact algorithm and material properties. The resultant forces at the input bar-sample interface are presented in Fig. 4.25. "3PB" represents three-point bending and "PC" represents pure compression.

Analysis of forces with respect to geometrical configuration Three-point bending (3PB)

1. As the sample starts to bend under load, the contact forces are written. Since, only a set of nodes are supposed to be in contact with the bars at corresponding interfaces, we see a difference between F in and F out .

2. In this case, both F in and F out drop down to zero, which signifies the rupture and loss of contact at the sample-bar(s) interface.

Pure compression (PC)

1. As opposed to rib sample, the input and output bars are in contact with all salve nodes in the salve segment throughout the simulation. Thus F in and F out superimpose.

2. In contrast to 3PB sample, the F in and F out do not drop to zero as the sample is always in contact even though it is subjected to damage and elements being eroded.

Discussions

The mesh density and mesh size have an influence on the numerical propagation of longitudinal waves and wave celerity c F EM . A mesh density of Γ = 20 per wavelength gives the best possible results yet not cost efficient in terms of calculation time. Therefore, for a better propagation of longitudinal waves, the element size in the direction of wave propagation was considered to be 1 mm.

The choice of constitutive law and damage law plays an important role in the results. We had pointed out some hurdles in experiments regarding behavior of porcine rib cortical and trabecular bone sample under dynamic environment. For example, due to miniature geometrical nature of porcine ribs, it is difficult to extract uniform rectangular or dog-bone samples for dynamic tension or compression tests. On the other hand, cylindrical samples of trabecular bone extracted from porcine rib might not capable of producing comprehensible experimental results. Moreover, as said in Sec. 4.2, the numerical porcine rib model used in this thesis is homogenized i.e. we do not consider the trabecular bone architecture. Therefore, we had no other choice rather to let go off using foam like behavior for trabecular bone constituents of porcine rib.

These all together , have led us to use a elastic-viscoplastic constitutive law coupled with an external damage law in order to simulate the time and path dependent progressive crack propagation until entire fracture. In this case, three material models were tested (*MAT-105,*MAT-024 and *MAT-003).

1. Out of which the first model was ruled out as there were penetrations of the Input rigid wall into the sample. It seems that the damage material parameter ς of material model *MAT-105 was the reason.

2. The response of ribs with properties from Pavier (2013) (Models-24a,24b and 3b) seem to respond well even though there is a significant amount of difference in the force magnitudes.

a Model-24a is a replication of its first version but the numerical sample is solid bone instead of a hallow tube as in [START_REF] Aubert | Étude du comportement des côtes sous sollicitations dynamiques[END_REF]. We strongly believe that the rigidity of the solid bone increases the force alongside a decrease in displacement. a Model-3a is an upgrade to its first version through a damage law that helps to differentiate between tensile and compressive failure instead of a simple Failure strain criteria as in [START_REF] Li | Rib fractures under anterior-posterior dynamic loads: Experimental and finite-element study[END_REF]. The force magnitudes are less than that of Models-24a,24b and 3c due to the difference in material properties of [START_REF] Pavier | Contribution à la compréhension des phénomènes physiques lors de l'impact d'un corps sur un modèle de structure biologique[END_REF] and [START_REF] Li | Rib fractures under anterior-posterior dynamic loads: Experimental and finite-element study[END_REF].

b Model-24c was executed to compare the material models *MAT-003 and *MAT-024. The crack propagation in this model is an evidence for the influence of kinematic hardening parameter β h of *MAT-003.

Part IV

Simulations

Chapter 5

Simulation at experimental scale

Résumé:

Ce chapitre présente l'application du modèle numérique mis en place dans le chapitre ( 4 Even though the rigid wall model produces results resembling to experimental cases, it is always better to carry simulations at experimental scale to have reliable and accurate comparisons. Thus the simulations at experimental scale were carried out with exact dimensions. Two numerical samples were considered, (i) single part porcine rib, and (ii) multi-part porcine rib. In both cases, the rib sample was meshed with under-integrated eight node solid elements of 0.15 mm in size and were provided with a stiffness based Flanagan-Belytschko hourglass was provided to avoid hourglass modes.

As concluded in §4.5, the bars and striker were meshed with under-integrated eight node solid brick elements of 1 mm in size along the direction of propagation of longitudinal waves. A stiffness based Flanagan-Belytschko hourglass formulation was used to avoid hourglassing modes.

An automatic-surface-to-surface contact with segment based penalty formulation was considered between every master and slave surface in the model. In addition to this, an interior contact algorithm was specified between each constituent of porcine rib if required. This chapter presents the different numerical tests carried out on single and multi part porcine rib samples as well as a generalized 5th human rib of 80 mm in length.

Single-part porcine rib sample

The numerical model consists of single part porcine rib sample, a short striker, an The next step is to calculate the dynamic variables involved at the interface for characterization of the sample. This is done using Eq. 2.17a to Eq. 2.17f, Eq. 2.18, and Eq. 2.19. Note that the curves presented in red, green and blue represent the variables at input bar-sample interface, output bar -sample interface 1, and output bar-sample interface 2 respectively. The curves in black represent either the sum or relative variables.

On the other hand, the energy absorbed by the sample was calculated through the methodology proposed by [START_REF] Lundberg | A split hopkinson bar study of energy absorption in dynamic rock fragmentation[END_REF]. The energy absorbed by the sample during fracture in three point bending can be expressed in terms of the energy of incident, reflected and transmitted stress pulses as follows

W P R = A b c b E b ˆt 0 σ 2 inc -σ 2 ref -σ 2 tr1 -σ 2 tr2 dt (5.1)
The Eq. 5.1 is modified to calculate the energy in three bars SHPB test.

Let us consider the time history of W P R and in Fig. 5.4d for further analysis.

Initially, W P R exhibits a linear evolution till t = 0.298×10 -3 [s], which corresponds to the energy for crack initiation W ci . From this point of time, we see fluctuation till entire rupture of the sample at t = 0.4998 × 10 -3 [s], corresponding to the energy for entire rupture W r .

The relative interface velocity V for this corresponding case is ≈ 12.5 [m.s -1 ],

thus it is logic to compare the F in vs U curve to the established experimental corridor-3.

Even though the energy and displacement values correspond to that of experimental values, we can see that the force amplitude is more than that of experimental. This proves that the trabecular bone plays an important role in determining the overall response and flexural strength of the rib. This is in opposition with that [START_REF] Granik | Human ribs: Static testing as a promising medical application[END_REF] had proposed.

Multi-part porcine rib sample

As seen above, a singe part porcine rib (i.e. with only cortical bone) is more rigid when compared to tested real porcine ribs. Therefore, a multi part porcine rib

Chapter 5 was tested. The trabecular bone elements were modeled using elastic-viscoplastic with kinematic hardening and cortical bone using elastic-viscoplastic with isotropic hardening. Since a model with different mechanical constants would be complicated to understand at this stage, the cortical and trabecular bone elements were provided with average material properties an in Tab. 3.4. The numerical setup with multi part rib resembles that as shown in Fig. 5.6. interior contact is added for different parts of porcine rib so that the solver can recognize that these parts form a single structure.

Case 1

The case 1 is with v str = 17 [m.s -1 ], so that the effect of trabecular bone can be pointed out. The raw strain signals in terms of absolute time t ′ for this case are presented in Fig. 5.7. The same methodology was followed for post processing of strain signals and the respective results are presented below. As the bar dimensions and material are the same, the arrival time t ′ 0 for multi-part sample is also 1.80 × 10 -3 [s].

Comparing the shifted strain signals of multi-part and single part porcine rib samples, we can notice that the reflected wave signal of multi-part (Fig. 5.8a) has changed a bit whilst the overall shape remains similar. On the other hand, the force responses at the interface show a pronounced change. For example, the Max(F in ) of single part was around 2000 [N] whereas, the Max(F in ) of multi-part is around 1200 [N]. This proves that the single part sample more rigid and takes less time for fracture (t f rc ≈ 3.46 × 10 -4 [s]) when compared to multi-part whose t f rc is prolonged to ≈ 5.80 × 10 -4 [s], thus proving that the trabecular bone has an effect on fracture of the structure.

Case 2

This case corresponds to a numerical striker velocity of 19.21 [m.s -1 ], which is the average velocity of all experimental cases in corridor 3 with fracture before dynamic equilibrium.

It is logic that the increase in striker velocity v str must increase the amplitude of strain in the bar. As the difference in v str is not huge, the increase in strain must not be large. This increase in v str will in turn increase the interface velocities. Comparing Fig. 5.9b and Fig. 5.10b, we can observe that the mean value of V for case 1 is ≈ 11.41 [m.s -1 ] where as that for case 2 is ≈ 12.45 [m.s -1 ]. Thus, the sample has further less time to establish an equilibrium with input and output forces. This is justified by a decrease in time for fracture t f rac (see Fig. 5.10a).

Similarly, the energy of incident wave W inc (t) is also increased due to the increase Section 5.2 in v str . Whereas, the energies required for crack initiation and rupture (W ci and W r ) have a negative correlation with striker velocity. The values presented in Fig. 5.9d and Fig. 5.10d for cases 1 and 2 respectively, speak for themselves.

Analysis of numerical response

Consistency of force response

As said in sec. 4.4, the contact algorithm does provide us the ability to have force at the interface. The forces obtained from shifted strain signal and contact force shall be similar, given that the contact algorithm has a significant representation of its physical counterpart. The Fig. 5.12a shows F inp (t) from contact algorithm of LS-Dyna and F inp (t) obtained from shifted strain signals through Eq.2.17a.

We can see that the initial linear response is the same for force obtained through both methods. After a certain time limit, fluctuations are observed in F inp from strain signals (as pointed by an arrow in Fig. 5.12a). This fluctuation is due to the accuracy of selection of starting point of incident and reflected waves of raw strain signal in shifting process. This accuracy issue is in fact a well known problem of SHPB experiments. Moreover, we do notice a time lag of ≈ 0.19 × 10 -3 [s] in the force obtained from strain signals. It is important to recall that the forces obtained from strain signals of SHPB is based on the assumption: "there is total contact between sample and bars at corresponding interfaces". Rather in reality, only a certain area of the sample is in contact and the area of contact is not constant throughout the test.

In simple terms, the difference between contact force and force from strain signals depends on the geometry of sample. A more detailed explication is given in sec. 4.4 and the results are presented in Fig. 4.25.

To make sure that the contact force is accurate than that from strain signals, the F inp vs U curves were compared to corresponding experimental results as in Fig. 5.12b.

As we can see that the linear elastic response of contact force superimposes with that of experimental case, this proves that the contact force is more physical than that of force from strain signals. Yet it is always better to compare the later with experimental response in order to have similar post processing methods. Before seeing the influence of scale factor curves on force response, we need to first understand the trend seen in all numerical and experimental responses corresponding to fracture before equilibrium. the Fig. 5.13b, shows F inp vs t response with all five scale factor curves.

The structure initially behaves as an elastic body (response to the left of the vertical dotted line) in Fig. 5.13b and continues to increase till it reaches maximum value of F inp . Till this point of time, the response obtained with each scale factor curve is almost similar. A second peak is observed, as the structure tries to establish an equilibrium with the output forces. This is where the damage starts to couple with flow stresses. As the interface velocity is high when compared to other experimental cases, the structure does not have sufficient time to establish the equilibrium. Thus, we see a rapid decrease in F inp till rupture of the sample. Therefore, a good candidate for the scale factor curve must produce the above characteristics (i.e. two peaks).

The curve 5 of Fig. 5.13a is unable to reproduce the second peak, which means that the rate of decrease in scale factor with strain rate is so brutal that the flow stresses couple to damage before hand and the rupture occurs way before the output forces come back to the impact point. Hence, the suitable curve for porcine rib is among curves 1 to 4.

Response obtained from each scale factor curve was then compared to corresponding experimental results to choose the suitable curve. The curve capable of reproducing similar results to experimental case was chosen see. The Fig. 5.13c

shows F inp vs U numerical response (solid lines with corresponding colors to scale factor curves used) and experimental result (represented with red pentagons).

It is clear that the curves 5,4 and 3 could be ruled out as either they failed to produced the second peak or a second peak of less magnitude. In contrast to curves 4 and 3, the curve 1 had a tendency to overestimate the second peak magnitude.

The best response is obtained with scale factor curve 2. Even though the magnitude of peak is overestimated, the displacement and time at fracture U f rac and t f rac is almost similar to that observed in experimental case.

Similar to force, a decreasing response can also be noted in the energy at fracture W r . The energy at rupture for the sample with curve 5 is very low that it is not in the established energy corridor. On the other hand, the energy with curve 3 and curve 4 could be barely accepted. Rather, it is curve 2 that gives acceptable The experimental data are represented by solid makers whilst, the numerical data are represented by red markers with no fill. We do notice that the numerical energy at rupture clusters alongside that of experimental case. Moreover, the numerical values are close to that of experimental mean (diamond and star black markers respectively).

We can notice two distinctive responses corresponding to single part porcine rib sample (inverted triangle red marker) and multi-part porcine rib sample case 2 (square red marker) in Fig. 5.14b. We see that these two samples needed less t f rac with respect to samples with similar W r . This proves that the composition and striker velocity play an important role in determining the time at fracture. We had already established that more the relative interface velocity V , the more fracture tends to down shift towards transition phase or the more it tends to occur before dynamic equilibrium. Moreover, the neither of the samples are of similar geometries. Therefore, it would be better to compare the correlation coefficient (r X,Y ) of displacement at rupture U r and energy at rupture W r for both experiments and numerical values. The correlation coefficient for experimental values (r U,W Exp.) is 0.8608 where as that of numerical values (r U,W Num.) is 0.7689.

Summary

The single part porcine rib results show that the sample is too rigid with respect to experimental results. Therefore, trabecular bone does play an important role in absorbing impact forces under these mechanical conditions.

The ensemble of mutli part simulations proves that the structure behaves in an elastic manner initially till it reaches the peak force and it depends on the distribution of mechanical properties and stiffness in the vicinity of impact. Later, the force starts to decrease till the point where the reaction forces arrive at the point of impact. This is the reason why we see an increase in input force. The second peak represents the maximum deflection the sample can sustain before crack initiation. The decrease in force after this point represents the propagation of crack and loss of structural rigidity.

The numerical simulation was carried out with a striker velocity of 17.39m/s, leading to a numerical interface velocity V of 10.46m/s. The numerical results are compared to an experimental case of similar conditions i.e. Corridor-3. The best numerical result corresponds to that obtained with input curve 2 of the damage law and is shown in Fig. 5.13a.

The differences between experimental and numerical response arise due to the difference in the geometry of samples tested.

Chapter 6 Discussions

The mechanical properties of porcine rib obtained through one variable and two variable models were compared. It is clear that the properties obtained from ρ app as independent variable are more representative than that obtained from either ash fraction α or properties as a function of ash fraction and volume fraction (α, V f ). Moreover, the constituents of porcine rib from apparent porosity values match to that of element distribution obtained through grey scale finite element mesh. It is rather interesting to see that the average mechanical apparent properties of rib constituents is of similar order of magnitudes to that used in past for human rib. This shows that the apparent mechanical properties for porcine rib are biologically reliable, given that the apparent density is correctly correlated with grey scale intensities. Yet proving that the porcine rib as human surrogate is an adequate choice.

The experimental set up used in this thesis had given important insights on the behavior of rib under dynamic loads at intermediate strain rates. The rupture of sample occurred either post equilibrium or pre equilibrium. In fact the rupture of the sample is dependent on rate of displacement V and so is the peak input force F inp . At low rates of displacement (i.e.) V ≤ 7 [m.s -1 ], the sample was capable to resist the force for large displacements. Yet as the rate of displacement increased, the peak of input force F inp increased and the displacement before fracture decreased. Another interesting point is that as V increases, the sample ruptured before establishing a dynamic equilibrium. Moreover, for high rates of displacements (i.e.) V ≥ 8 [m.s -1 ], the sample attains its peak, followed by a subtle increase in force and a constant plateau like response. This plateau response corresponds either to plastification or damage of the structure or both. Upon further increasing the rate of displacement (i.e.) V ≥ 11 [m.s -1 ], the sample attains its peak, followed by a subtle increase in force and a drastic decrease in force. This proves that the response of rib depends on rate of displacement V and therefore, strain rate ε. These results have led to establish the parameters such as time at fracture t f rac , energy at fracture W r and displacement at fracture U r . For rates of displacements between 8 [m.s -1 ] and 10 [m.s -1 ], the mean values of the aforementioned variables are t f rac ≈ 1.15 × 10 -3 [s], W r ≈ 5.12 [J] and t f rac ≈ 8 × 10 -3 [m]. On the other hand, for rates of displacements between 11 [m.s -1 ] and 13 [m.s -1 ], the mean values of the aforementioned variables are t f rac ≈ 0.55 × 10 -3 [s], W r ≈ 4.89 [J] and t f rac ≈ 6 × 10 -3 [m]. The fracture modes observed in experimental tests is due to a crack originating from tensile surface due to global bending of the sample, which tend to propagate towards the compression surface. The rupture is not influenced by the maximum force imposed on the sample. It rather depends on the local displacement of the rib.

A close look to the experimental data reveals that the rate of loading has a more pronounced influence on post yield behavior than on the linear elastic part. This suggests that the rib exhibits an elastic-viscoplastic behavior. Three different material models with elastic-viscoplastic capability were chosen initially. The material model MAT-105, which has damage law coupled internally is more efficient in calculation time than the other two. Unfortunately, this model was dropped as it could not consider the effects of strain rate on fracture strain.

Meanwhile, the literature suggests that the trabecular bone exhibits kinematic hardening behavior that is more pronounced than that exhibited by its counterpart cortical bone. Therefore, the material model *MAT-003 with kinematic hardening and MAT-024 with isotropic hardening was chose for trabecular and cortical bone respectively. These two material models were coupled to an incremental and stress-state dependent external damage law. Thanks to the curves provided by user, it can consider the variation of fracture strain with both stress-state and strain rate.

The parameters used for damage law permits to work in such a way that, when static yield stress of the material is reached, the damage D starts to accumulate. At the instance when D = D crit , damage starts to couple with the flow stresses.

The dynamic fracture strain ε f d , influences the accumulation of D. The more ε f d decreases with increase in ε, the more drastic is the increase in damage accumulation and plastic strain. The ultimate damage value D = 1 is attained through the non-linear coupling between damage D, dynamic fracture strain ε f d (η, ε), the plastic strain ε p and damage coupled flow stress. When this value is reached (D = 1), the element does not exist in the calculation, thereby initiating a crack. This continues until entire rupture of the sample.

Regarding numerical simulations at experimental scale, two different numerical samples were tested (i) the porcine rib with an assumption of only cortical bone, and (ii) the porcine rib was assumed to consist a distribution of only one cortical and one trabecular constituents with average of mechanical properties. In either case, the static tensile yield given in the material card of rib constituent is reached at point corresponding to the peak of the contact force at the input bar and sample interface. Yet this proves that the interface displacement rate or the velocity of the corresponding interface does not influence the initial linear response.

The fracture of samples is a pure local phenomenon. The fracture process is related to wave propagation and forces. From the point of impact till the peak, the response is elastic, which depends on the local distribution of the mechanical properties. From this instance till the second peak, we see a complex three dimensional wave propagation in the sample. The increase in F in before the second peak is due to the arrival of reflected waves from the sample-output bar interfaces. This increase in force signifies the onset of coupling between damage and flow stresses. It is at this point, the waves reflected from anterior end-output bar interface travel towards the posterior end-output bar interface and vice-versa. These waves traveling in opposite directions are the ones responsible for subjecting the sample to further tension and fracture.

Part V

Conclusions and Perspectives

Chapter 7

Conclusions and perspectives

Conclusions

Rib fractures are important when it comes to chest trauma. The consequences of rib fracture shall lead to complex medical conditions and eventual death. The existing trauma evaluation scores do give critical insights to the condition of victim but they are not applicable prior to a radio-logical examination. Moreover, these so called evaluation scores are not applicable for single ribs. Thereby driving to a need for comprehension of rib fracture under dynamic loads.

Most of the blunt trauma are due to road accidents. Other driving factors are ballistic impacts, sport accidents and industrial accidents. The advances in computational capabilities have led to the development of several numerical models of the human thorax. Even though these models are being used quite often, there is a lack of geometrical representation, mechanical behavior and fracture prediction to be precise. The goal of this thesis was to carry numerical simulations on isolated porcine ribs with a feasible geometrical representation and realistic constitutive law coupled to damage law.

As this thesis deals with isolated porcine, it is always better to characterize the different bone tissues from anatomy point of view. Speaking of rib anatomy, a rib is composed of two different bone tissues at macroscopic level: cortical and trabecular bone. These two tissues are different in terms of mechanical behavior and properties. Therefore, care must be taken so that a reliable constitutive law.

Literature says that the trabecular bone does not play a significant load in the load carrying capacity. Speaking of considering the rib as a structure has an inconvenience as we cannot differentiate the role played by two bone in precise manner. Moreover, there is a lack of experimental database for ribs subjected to dynamic loads at intermediate and high strain rates.

The experimental setup used in this study had led us some key insights on the dynamic behavior in three point bending. In order to differentiate the waves involved in split hopkinson pressure bar and for impedance adaptation, bars made of nylon were used. Even though the split hopkinson pressure bar setup is widely used in dynamic characterization of materials and structures, it imposes three main difficulties: (i) impedance adaptation, (ii) time and space shifting of waves and (iii) non-plane wave effects. Theses three issues of critical importance must be verified prior to experimental characterization, if not they lead to either poor characterization of sample of incoherent results. Since the ribs used during the experimental tests are not the same, there will be variations in the experimental results. Despite the difficulties, we were able to present three different kinds of force response.

The force displacement curves obtained from experimental characterization proved that the rate of displacement had a more pronounced influence on the post yield behavior of the sample. Therefore, an elastic-viscoplastic constitutive law was chosen. Out of the material models presented in this thesis, the models MAT-024 and MAT-003 were chosen as potential candidates. Based on the literature study, we were able to decide that the trabecular bone exhibits more kinematic hardening than that of cortical bone. Therefore, isotropic hardening was chosen for cortical bone (MAT-024) and kinematic hardening was chosen for trabecular bone (MAT-003). A yield stress scaling based on Cowper-Symmonds model was used to consider the effect of strain rate on the post yield behavior of both cortical and trabecular bone.

Even though, these material models are capable of simulating failure, they are not sufficient enough for considering crack propagation till entire rupture. Therefore, we have decided to couple these models with an external damage law.

The chosen damage law is stress state dependent and does consider the effects of strain rate on the fracture strain. The lack of possibility to obtain uniform and reliable samples of cortical and trabecular bone from porcine rib is a major constrain to characterize their behavior in dynamic tensile loading. Therefore, we had chosen a relation put forth for cranial bone for the reduction in fracture strain with increase in strain rate. This choice of critical importance was made as the structure at macroscopic levels of cranial bone and rib resemble. A total of 5 variations of this curve were tested. It is evident that the reference curve of cranial bone is not suitable for rib as it drove to a much less rigid response of the rib with respect to corresponding experimental data. Therefore, a less decreasing curve was selected and the one that gave the best response was ε f d = 0.63 -1.25 × 10 -2 Ln ε.

Regarding the contact between bars and bone, a segment based penalty contact with a static friction coefficient of 0.1 was chosen. The segment based penalty contact is applicable here in order to adapt between the difference in material properties of bone and bar. The segment based contact is more precise for determining the contact force as it removes the initial penetration of each time step iteration before calculating the contact force of the current iteration. This is why it is expensive in terms of calculation time. As the porcine rib sample consists of elements with different mechanical properties and behavior (i.e. a composite structure), another contact is necessary to represent the internal cohesion of the structure, which helps the solve to understand that the sample is composed of different material properties. This in turn avoids a negative volume elements.

Concerning the numerical simulations at experimental scale, the proposed constitutive law coupled to damage law, with the help of contact is able to consider the effects of strain rate in both behavior and in fracture. First of all, the numerical force response are similar to that of experimental (both in terms of time and displacement). The mean fracture time, displacement at rupture and energy at rupture are approximately the same for cases with fracture before equilibrium. Thanks to the numerical simulations, we can deduce that the energy for rupture before equilibrium is around 3.8 [J]. Given that the rate of displacement imposed on the sample at the bar interfaces is not brutal enough, the sample might need more energy , time and displacement for rupture. Therefore, the rate of displacement or velocity of impact do play an important role.

Therefore, the current study indicates successful FE implementation of complicated constitutive and damage law with the consideration of strain rate effects and non-linear damage accumulation.

Limitations of the proposed methodology

Despite the consistent performance of the mode presented in this thesis, there are quite a few limitations. For example,

• It is isotropic i.e. it does not consider the effects along other axis even though it is three dimensional.

• It does not consider the difference of behavior in compression and tension.

• The parameters used of material properties and constitutive laws are optimized form of existing mathematical models.

In order to overcome these limitations, I suppose that the following experimental steps could give us important insights

• Tests on cortical and trabecular samples in different stress state modes and directions

• Implementation of a user defined model based on triaxiality both in constitutive law and damage law i.e. damage variable shall be considered as a vector rather than scalar.

• Dynamic tensile and compression tests on porcine bone coupons for the effects of ε f ( ε).

Perspectives

Even though the correlation between the experimental and numerical fracture pattern is adequate, an accurate representation of a more realistic fracture behavior requires model improvements such as: bone heterogeneity, transmission of forces between the different constituents of the rib. Moreover, optimization of the damage law parameters through respective experimental data, will further improve the efficiency of the proposed FE model.

The improved porcine rib model shall be tested on a numerical human rib so that we could understand the critical parameters for human ribs under similar mechanical environments. Once this is done, it could be implements to numerical thorax models to give insights on injury severity, thereby contributing to betterment of passenger safety or to a link between prediction of critical condition and treatment to be provided at trauma centers. Thus, the virtual work over n elements can be written as follows where, N is the interpolation matrix, B is the strain displacement matrix, σ is the stress tensor, a is the acceleration vector, b is the body force vector and t is the applied traction load.

node hexahedral element

Let us consider an eight node solid hexahedral element as in Fig.

A.2. The Eq.

A.16 can be written as ements with an explicit solver for dynamic applications or applications involving large deformations. Yet the use of under integrated elements requires use of hourglass properties in order to avoid the Hourglass modes.

x i (X α , t) = x i (X α (ξ, η, ζ), t) =

Hourglass

The use of under integrated elements such as one point integration, leads to the presence of zero energy modes commonly known as "Hourglassing Modes". These modes are oscillatory and tend to have periods less than that of the structure's response. On the other hand, the hourglass modes with time periods similar to structural response are known to have negligible effects on the response. Thus, the hourglassing modes with less time periods have to be omitted. This could be either

Type I

In order to understand the difference between these formulations and why hourglassing modes exist, let us consider a 8 node hexahedral element with only one Gaussian integral point. This element can be related to its spatial coordinates in terms of nodal coordinates using the isoparametric shape functions as follows (Flanagan and In other words, the one point integration is based on the assumption that the mean stress depend only on mean strain but not the nonlinear terms. The mean kinematic terms uij in Eq. A.34 are expressed in terms of the displacement field over the element as uij = 1

V ˆV ui,j dV (A.35) This leads to the definition of B matrix as follows:

B ik = ˆV φ ik dV (A.36)
Thorough mathematical manipulations, [START_REF] Flanagan | A uniform strain hexahedron and quadrilateral with orthogonal hourglass control[END_REF] had concluded that in a one point integral framework, only the base vectors lead to stresses or nodal forces. These base vectors are given in Tab. A.1.

Node Π = 1 Π = 2 Π = 3 Π = 4

Γ Π1 1 1 1 1 Γ Π2 1 -1 -1 -1 Γ Π3 -1 -1 1 1 Γ Π4 -1 1 -1 -1 Γ Π5 -1 -1 1 -1 Γ Π6 -1 1 -1 1 Γ Π7 1 1 1 -1 Γ Π8 1 -1 -1 1
Table A.1 -Hourglass base vectors and modes for an eight node hexahedral element [START_REF] Flanagan | A uniform strain hexahedron and quadrilateral with orthogonal hourglass control[END_REF] For a solid element with one point integration, when the opposite nodes have identical velocities, the strain rate εij = 0 due to the anti-symmetric property of strain matrix. The shape functions of a hexahedral solid element with one point

Types IV and V

The types II and III are executed by adding certain viscosity to the elements. This unfortunately absorb more energies thereby, leading to discrepancies in the total energy involved in the model. This can be avoided by replacing the additional viscosity with springs of certain stiffness. Therefore, we have the hourglass types IV and V that are mathematically identical to types II and III represented by equations A.42, A.43 and A.44.

In general, a hourglass with stiffness form is used in problems involving large deformations (automotive crash environments for example) where as, a viscous from hourglass is chosen in problems involving deformations at high velocities. Thus, the a hourglass of stiffness from was chosen for both bone and bar materials in this thesis.

Thus, we can write the intersection of boundaries at t 0 and their intersection at t i as in Eq. B.1 and Eq. B.2 respectively. δB 1 ∩ δB 2 = ∅ (B.1)

δb 1 ∩ δb 1 = ∅ (B.2) (B.3)
Therefore, the expression of C at t i becomes

C = δb 1 ∩ δb 2 (B.4)
A general rule of thumb for slave and master surfaces involved in contact is that a coarse mesh surface is always declared as master (say the B 1 ) and the other as salve surface (B 2 ). Therefore, all the nodes on the discretized boundary δb 1 are treated as master nodes and the others as slave nodes (on the boundary δb 2 ). This is the basis of contact algorithms available on LS-Dyna.

Contact between sliding or impacting interfaces can be implemented by three different methods in LS-Dyna :

1. Kinematic Constraint Method (KCM): The constraints are imposed on the global equations of the system through transformation of nodal displacement of slave surface on the contact interface. This transformation, eliminates the normal degrees of freedom on the slave nodes. In order to retain the efficiency of the explicit time integration, the mass is lumped in such a way that only nodes of master surface are coupled. The disadvantage of this method is lack of accuracy when the master surface has a mesh finer than slave surface.

This leads to unnecessary penetration of master nodes into the slave surface

Penalty Formulations

Standard Penalty Formulation

This is most commonly used penalty contact formulation. On implementation, each and every slave node n s in contact is searched for penetration through master surfaces. Once identified, two cases can show up 1. In the case of no penetration of n s , nothing is done and the calculation is advanced to next iteration 2. If penetration occurs, an interface force between n s and contact point is applied which is proportional to the amount of penetration as follows.

• the penetration is identified mathematically as

P = n i × [t -r(ξ c , η c )] < 0 (B.10)
here, n i is normal to the master segment at contact point.

• an interface force f s is added to n s and a force f i m to each node of the corresponding master segment S i . These forces are expressed in equations below f s = -Pk i n i if P < 0 (B.11)

f i m = φ i ((ξ c , η c ))f s if P < 0 (B.12)
here, k i is the stiffness factor for master segment and is expressed as

k i = f si K i A 2 i V i (B.13)
where, f si is the stiffness factor for the interface (0.01 by default). K i , A i and V i are the bulk modulus, area and volume of the corresponding master segment respectively.

Whist used for contact between material with different bulk modulus, the contact stiffness f si can be under estimated thereby making a weak contact. Therefore, the standard penalty formulation might require manual adjustments in order to avoid instabilities.

Soft Constraint Penalty Formulation

The soft based penalty formulation was developed as a counter action to standard penalty formulation. This formulation is well suited for application to materials with huge difference in bulk modulus.

An additional stiffness based on Courant's stability criterion of the local masses (master and slave) connected by a spring is used. This additional stiffness is ex-

Segment based Penalty Formulation

The contact stiffness in this approach is calculated similar to that of Soft constraint.

The only difference is that the mass function m is the mass of master and slave segments as opposed to mass of master and slave nodes.

k cs (t) = 0.5 • SLSFAC         SFS or SMS         • m 1 m 2 m 1 + m 2 1 ∆t c (t) 2 (B.16)
where,m 1 and m 2 are masses of master and slave segments in contact. In case of solid elements, the segment masses are one half of element mass. SFS and SMS are penalty scale factors for slave and master segments respectively, SLSFAC is the scale factor for sliding interfaces in contact.

Another difference is that ∆t c is updated only when the time step of solution is greater than 5%∆t c .

  from the École National Supérieur d'Arts et Métiers for having accepted and granted their valuable time to review my PhD thesis. I would equally like to express my sincere gratitude to the members of the jury Pr. Jean-François GANGHOFFER of Laboratoire d'étide des Microstructures et de Mécanique des Matériaux (LEM3) from the École Nationale Supérieur d'Electricité et de Mécanique (ENSEM), to Pr. Salah NAILI of Laboratoire Modélisation et Simulation Multi Echelle (MSME) from the Université Paris Est Créteil (UPEC), to Dr. Julien PAVIER of Nexter Munitions and to Dr. Nicolas PRAT of Institut de Recherche Biomécanique des Armées (IRBA) for having granted their time to exterminate my thesis and for having participated in the jury. I would like to than the French Région Centre Val de Loire for having provided the funding required for this thesis. Meanwhile, I would like to express my sincere gratitude to Dr. Eric LESPESSAILLE (MD,PhD,PAST) of the Centre Hospitalier Régional Orléans for his valuable assistance in the High Resolution peripheral i Quantitative Computer Tomography (HR-pQCT) procedure. I would also like to express my sincere gratitude to former head of the M2RI research masters program at polytech Tours, retired Pr. Narayanaswami Mohan Ranganathan , to Pr. Mame Williams-Louis of Université d'Orléans from Laboratoire Pluridisciplinaire de Recherche Ingénierie des systèmes, Mécanique, Énergétique, to Pr. S.H. Venkatasubramanian and to Pr. Ramesh D. Kandadai of SRM University, India for having encouraged me to pursue my career in research and to do my PhD.

1

 1 . . . . . . . . . . . . 3.7 Comparison of α from model 1b and 1c . . . . . . . . . . . . . . . . . 3.8 Mechanical properties of porcine rib with ρ α . . . . . . . . . . . . . . 3.9 variation of V f and ρ tissue with α . . . . . . . . . . . . . . . . . . . . 3.10 Mechanical properties of porcine rib with two variable model . . . . . 3.11 Different density measures for porcine rib . . . . . . . . . . . . . . . . 3.12 Comparison of architectural properties . . . . . . . . . . . . . . . . . 3.13 Comparison of mechanical properties for porcine rib . . . . . . . . . . 3.14 Element distribution in porcine rib sample . . . . . . . . . . . . . . . 4.1 Variation of celerity ratio with mass matrices . . . . . . . . . . . . . . 4.2 Incident strain signals at mid point of the bar . . . . . . . . . . . . . 4.3 Finite element mesh of porcine rib and bars (a) and zoom of the left rib extremity showing the hexahedral elements (b) -dimensions in mm102 4.4 Yield stress scaling MAT-024 . . . . . . . . . . . . . . . . . . . . . . 4.5 Displacement boundary condition . . . . . . . . . . . . . . . . . . . . 4.6 Effect of β h *MAT-003 . . . . . . . . . . . . . . . . . . . . . . . . . . 4.7 Kinematic and Isotropic hardening *MAT-003 . . . . . . . . . . . . . 4.8 Comparison of *MAT-003 and *MAT-024 . . . . . . . . . . . . . . . . 4.9 Input curves of damage law . . . . . . . . . . . . . . . . . . . . . . . 4.10 Model for damage law validation . . . . . . . . . . . . . . . . . . . . 4.11 Effects of damage law parameters . . . . . . . . . . . . . . . . . . . . 4.12 Numerical Porcine Rib Sample used for model considerations . . . . . 4.13 Penetration of Input wall into the sample . . . . . . . . . . . . . . . . 4.14 Crack growth in Model-24a . . . . . . . . . . . . . . . . . . . . . . . 4.15 Contact Forces in Model-24a . . . . . . . . . . . . . . . . . . . . . . . 4.16 Crack growth in Model-24b . . . . . . . . . . . . . . . . . . . . . . . 4.17 Contact Forces in Model-24b . . . . . . . . . . . . . . . . . . . . . . . 4.18 Crack propagation in Model-24c . . . . . . . . . . . . . . . . . . . . . 4.19 Contact Forces in Model-24c . . . . . . . . . . . . . . . . . . . . . . . 4.20 Crack propagation Model-3a . . . . . . . . . . . . . . . . . . . . . . . 4.21 Contact forces Model-3a . . . . . . . . . . . . . . . . . . . . . . . . . 4.22 Crack propagation Model-3b . . . . . . . . . . . . . . . . . . . . . . . 4.23 Contact forces Model-3b . . . . . . . . . . . . . . . . . . . . . . . . . 4.24 Influence of contact surface on force . . . . . . . . . . . . . . . . . . . 4.25 F vs t form contact force data . . . . . . . . . . . . . . . . . . . . . . 5.1 Single part porcine rib with SHPB at experimental scale . . . . . . . 5.2 Raw numerical strain signals . . . . . . . . . . . . . . . . . . . . . . . 5.3 Shifted strain signals of single part porcine rib . . . . . . . . . . . . . 5.4 Dynamic variables and energy for single part porcine rib . . . . . . . xx 5.5 Comparison of single part to experimental corridor-3 . . . . . . . . . 5.6 Multi-part porcine rib with SHPB at experimental scale . . . . . . . . 5.7 Raw numerical strain signals for case 1 . . . . . . . . . . . . . . . . . 5.8 Shifted strain signals of multi-part porcine rib case 1 . . . . . . . . . 5.9 Dynamic variables and energy at interfaces for case 1 . . . . . . . . . 5.10 Dynamic variables and energy at interfaces for case 2 . . . . . . . . . 5.11 Comparison of case 1 (left) and case 2 (right) with experimental data 5.12 Consistency of force response . . . . . . . . . . . . . . . . . . . . . . 5.13 Effect of scale factor curves . . . . . . . . . . . . . . . . . . . . . . . 5.14 Consistency with experimental data . . . . . . . . . . . . . . . . . . . A.1 Body in a Cartesian system . . . . . . . . . . . . . . . . . . . . . . . A.2 8 node hexahedral element . . . . . . . . . . . . . . . . . . . . . . . . B.1 Two bodies in R 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B.2 Surface treatment in kinematic constraint method . . . . . . . . . . . B.3 Slave and Master nodes for slave search . . . . . . . . . . . . . . . . . B.4 Unique normal of a master segment S 1 . . . . . . . . . . . . . . . . . Different scoring systems . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Hierarchical Organization of Bone . . . . . . . . . . . . . . . . . . . . 1.3 α based equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4 ρ α based equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.5 Young's modulus E [GPa] as a function of apparent density as put forth by Rice et al. (1988) . . . . . . . . . . . . . . . . . . . . . . . .

  . . . . . . . . . 1.10 Ratio of tensile to compressive properties of bone . . . . . . . . . . . 1.11 Average data for bovine femur bone . . . . . . . . . . . . . . . . . . . 1.12 Variation of tensile properties of cranial bone with ε, [s -1 ] . . . . . . 1.13 Existing data on human ribs under three point bending . . . . . . . . 1.14 Data of existing isolated human rib FE model . . . . . . . . . . . . . 1.15 Existing numerical models details for human rib in LS Dyna . . . . . 1.16 Mechanical properties of numerical porcine rib . . . . . . . . . . . . . 2.1 Mean and standard deviation values for all properties . . . . . . . . . 3.1 Grey scale values for porcine rib . . . . . . . . . . . . . . . . . . . . . 3.2 Different constituents of Porcine rib based on ρ app . . . . . . . . . . . 3.3 Power law regression equations used for mechanical properties (unit for stresses is here MPa) . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Comparison of average porcine rib apparent properties with human rib 4.1 Elastic properties of Nylon for bars in SHPB . . . . . . . . . . . . . . 4.2 Effect of Γ on accuracy and calculation time . . . . . . . . . . . . . . 4.3 Key factors of elastic-vicoplastic models in LS-Dyna . . . . . . . . . . 4.4 Summary of simulations with single part porcine rib and rigid walls . 4.5 Damage law parameters of porcine rib . . . . . . . . . . . . . . . . . 4.6 List of different ε f d /ε f s vs ε relations tested . . . . . . . . . . . . . . 4.7 Summary of simulations with single part porcine rib and rigid walls . 4.8 Properties of cortical bone for Model-105 . . . . . . . . . . . . . . . . 4.9 Damage parameters of Cortical Bone . . . . . . . . . . . . . . . . . . 4.10 Mechanical properties for single part model-24 . . . . . . . . . . . . . 4.11 Tensile fracture strain for single part model-24 . . . . . . . . . . . . . xxiii A.1 Hourglass base vectors and modes for an eight node hexahedral High Resolution -peripheral Quantitative Computed Tomography

  Figure 1.1 -Lateral and anterior perspectives of thoracic cage[START_REF] Gray | Anatomy of the human body[END_REF]

Figure 1 . 2 -

 12 Figure 1.2 -Sternum and Thoracic Vertebrae[START_REF] Gray | Anatomy of the human body[END_REF]

Figure 1 .

 1 Figure 1.4 -Evolution of rib cross-section

•

  Anatomical scoring systems (a) Abbreviated Injury Scale (AIS) (b) Injury Severity Score (ISS) (c) New Injury Severity Score (NISS). • Physiological scoring systems (a) Revised Trauma Score (RTS) (b) Acute Physiology And Chronic Health Evaluation (APACHE).

  The porosity associated with cortical bone are Vascular Porosity (VP) formed from the canals carrying bone fluids, nerves (osteon canals and Volkmann canals), Lacunar-Canalicular Porosity (LCP) arising from the spaces in canals such as lacunae and canaliculi and Collagen-Apatite Porosity (CAP) associated with the spaces between collagen and hydroxyapetite (S.C.Cowin, 1999) (see Fig.1.6b). Cortical bone exists at the shafts of long bone such as femur or as an external shell around the trabecular bone in flat or short bones such as ribs or vertebrae.

  Figure 1.6 -Cortical bone

Fig. 1 .

 1 Fig. 1.7d).

Figure 1 . 7 -

 17 Figure 1.7 -Trabecular bone at different anatomical sites[START_REF] Keaveny | Bone Mechanics[END_REF]

Figure 1 .

 1 Figure1.9 -Loading-reloading response of cortical and trabecular bone[START_REF] Keaveny | Bone Mechanics[END_REF]

Figure 1 .

 1 Figure 1.10 -Data from[START_REF] Granik | Human ribs: Static testing as a promising medical application[END_REF]

Stage 1

 1 the striker comes into contact with the external cortical bone in rib and an equilibrium is established with supports (Fig. 1.12a) Stage 2 the cortical bone under the striker, caves in due to local damage (Fig. 1.12b) Stage 3 due to further displacement, the caved in cortical bone clamps itself with the striker (Fig. 1.12c) Stage 4 further buckling occurs (Fig. 1.12d) Stage 5 rupture occurs from the interior cortical bone end (Fig. 1.12e).

Figure 1 .

 1 Figure1.12 -Different stages of response[START_REF] Charpail | Analyse du comportement mécanique des côtes humaines en dynamique[END_REF]

Figure 1 .

 1 Figure 1.13 -Characteristics of human rib in three point bending

Figure 1 .

 1 Figure 1.15δ > 4mm Aubert (2012)

Figure 1 .

 1 Figure 1.16 -Compressive properties of cylindrical human rib samples

  1.17b) in Radioss solver.

  (a) Lateral Medial bending (b) Anterior-Posterior bending

Figure 1 .

 1 Figure1.17 -Three point bending models[START_REF] Charpail | Analyse du comportement mécanique des côtes humaines en dynamique[END_REF]

  Figure 1.18 -Other Anterior-Posterior bending models

  and the results were analyzed in the present thesis work. The experimental setup is based on Hopkinson pressure bars, which permits us to study the response of sample in three point bending. It consists of three bars, an input bar and two output bars. The sample is in fact placed between input and output bars. The input bar subjects the sample to mechanical loads by transmitting a compression wave to the sample. The bars are made of nylon in order to have mechanical impedance of similar orders of magnitude as that of the sample. The analysis of the sample response is done through dynamic variables such as force at interfaces and strain rate etc.,. These dynamic variables are based on the elastic strain waves measured in the input and output bars.

  2.1): 1 Waves in striker (a) a compression wave propagating towards the free end of the striker (b) a traction wave reflected from the free end, that cancels the stress imposed by the compression wave. 2 Waves in input bar (a) Incident Wave (ε inc ): A wave traveling in the input bar towards the input bar -sample interface. The difference in the mechanical impedance of the sample and bars gives rise to the separation of (ε inc ) into two waves. (represented in solid red line in Fig. 2.1). (b) Reflected Wave (ε ref ): A part of the incident wave, that propagates in the opposite direction to the incident wave. (represented in dotted red line in Fig. 2.1). 3 Waves in output bar(s) (a) Transmitted wave (ε tr ): The other part of incident wave that passes through the sample into the output bar (represented in solid green line in Fig. 2.1). The bars in SHPB are characterized by their density (ρ b ), Young's modulus (E b ) and cross-sectional area (A b ). All these variables together, influence the mechanical impedance (Z b ).

Figure 2 . 1 -

 21 Figure 2.1 -The Lagrange diagram for a typical SHPB setup

Figure 2 . 2 -

 22 Figure 2.2 -The modified three-point bending SHPB setup

Figure 2 .

 2 Figure 2.3 -Reflective bands placed on striker and optical sensors

Figure 2

 2 Figure 2.4 -P can vs v str for long striker[START_REF] Pavier | Contribution à la compréhension des phénomènes physiques lors de l'impact d'un corps sur un modèle de structure biologique[END_REF]

  the optical sensor capture any motion of striker. The acquired signals are then treated using a software DAVID, developed at the "Laboratoire de Mécanique des Solides de Palaiseau" (LMS). The DAVID software, permits us to reconstruct the strain signals (shifting procedure) by taking the dispersion correlation and attenuation of the bars used in the SHPB. The user must execute a supplementary test with the bars (i.e. without sample) in order to calibrate the dispersive and dissipative behavior of SHPB bars.

Figure 2 .

 2 Figure 2.6 -Signal acquisition setup

Figure 2 . 7 -

 27 Figure 2.7 -Initial Position of the sample

Figure 2 .

 2 Figure 2.8 -Curved beam approximation of a rib

Figure 2 . 9 -

 29 Figure 2.9 -Normal distribution of all samples tested 61

Figure 2 .

 2 Figure 2.10 -Test without fracture

Figure 2 .

 2 Figure 2.11 -Experimental data for v str = 9.27 m.s -1

Figure 2 .

 2 Figure 2.12 -Test with fracture during dynamic equilibrium

Figure 2 .

 2 Figure 2.13 -Experimental data for v str = 17.39 m.s -1

Figure 2 .

 2 Figure 2.14 -Test with fracture before dynamic equilibrium

Figure 2 .

 2 Figure 2.15 -Experimental data for v str = 19.21 m.s -1

Figure 2 .

 2 Figure 2.16t f rac vs v str

Figure 2 .Figure 2 .

 22 Figure2.17 -Energy at fracture W r vs displacement U r

  3, 11 m.s -1 ≤ V ≤ 13 m.s -1

Figure 2 .

 2 Figure 2.20 -Experimental force displacement for Case 1 and Case

Figure 2 .

 2 Figure 2.21 -Samples grouped with respect to experimental corridors 75

  DICOM image slices were obtained using a multi-sclice device at 120 kVp, 160 mAs with an 512 × 512 image matrix and 150 × 150 mm field view. Each image sclice was 7 mm thick with a spatial resolution of 0.25 mm × 25 mm as in Fig. 3.1.

Figure 3 . 1 -

 31 Figure 3.1 -Treated porcine rib and CT scan images

ρFigure 3 . 2 -

 32 Figure 3.2 -Correlation of ρ app with respect to CT P R

Figure 3

 3 Figure 3.3 -Cross-sectional views of the rib FE model showing the density distribution.

  Figure 3.4 -V f app and Po app vs CT P R

  Figure 3.5 -Mechanical properties of porcine rib with ρ app

Figure 3

 3 Figure3.6ρ α vs ρ app of porcine rib based on[START_REF] Keller | Predicting the compressive mechanical behavior of bone[END_REF]

  σ ult,T en [MPa] = 3.079 × 10 -5 ρ 2.27 α (3.10) Model 1c: Similar to Model 1b, this model also considers ρ α as control variable. The only difference is that the ρ α is interpolated from the Equivalent Mineral density of porcine rib, which is determined from CT P R .ρ EQM [kg.m -3 ] = 0.6452CT P R + 0.487 (3.11)ρ αEQM [kg.m -3 ] = 1.27ρ EQM + 0.0452 (3.12)

Figure 3 . 7 -

 37 Figure 3.7 -Comparison of α from model 1b and 1c

  Figure 3.8 -Mechanical properties of porcine rib with ρ α

Figure 3

 3 Figure 3.9 -variation of V f and ρ tissue with α

  Figure 3.10 -Mechanical properties of porcine rib with two variable model

FirstlyFigure 3

 3 Figure 3.11 -Different density measures for porcine rib

Figure 3

 3 Figure 3.13 -Comparison of mechanical properties for porcine rib

1 -

 1 Elastic properties of Nylon for bars in SHPB

Fig

  Fig. 4.1

Figure 4 . 1 -

 41 Figure 4.1 -Variation of celerity ratio with mass matrices

Figure 4 . 2 -

 42 Figure 4.2 -Incident strain signals at mid point of the bar

1 .

 1 Elastic-viscoplastic constitutive law with internally coupled damage law (*MAT-105) 2. Modified elastic-viscoplastic constitutive law with kinematic and isotropic hardening (*MAT-003) 3. Modified elastic-viscoplastic constitutive law with isotropic hardening only (*MAT-024) Elastic-Viscoplastic with Continuum Damage Mechanics (*MAT-105)

Figure 4 . 4 -

 44 Figure 4.4 -Yield stress scaling MAT-024

  Figure 4.5 -Displacement boundary condition

Figure 4

 4 Figure 4.6 -Effect of β h *MAT-003

2 .

 2 GISSMO (Generalized Incremental Stress-State dependent damage MOdel) 3. DIEM (Damage Initiation and Evolution Model).

  Figure 4.9 -Input curves of damage law

Figure 4 .

 4 Figure 4.10 -Model for damage law validation

Figure 4 .

 4 Figure 4.11 -Effects of damage law parameters

  4.12a) with rigid walls instead of bars to economize time as in Fig. 4.12b (a) Single Part (b) Numerical model

Figure 4 .

 4 Figure 4.12 -Numerical Porcine Rib Sample used for model considerations

7 -

 7 2.1), coupled with damage law if necessary were tested with the contact algorithm of Sec. 4.4. The different models tested are presented below. Summary of simulations with single part porcine rib and rigid walls 4.3.1 Single Part Model-105 Since the material model *MAT-105 has damage law readily coupled to constitutive law, it would be interesting to test it. Therefore, we had embarked with this model. The material properties of cortical bone were considered from ρ app and are shown in Tab. 4.8. In addition to these mechanical properties, three additional damage parameters are required for complete implementation of *MAT-105. 1 Damage threshold(r D ): The value of effective plastic strain, at which material softening begins. parameter D c : That signifies the damage variable value at which element(s) are deleted from the calculation.

Figure 4 .

 4 Figure 4.13 -Penetration of Input wall into the sample

  Figure 4.14 -Crack growth in Model-24a

  Figure 4.15 -Contact Forces in Model-24a

  Figure 4.16 -Crack growth in Model-24b

  )

  Figure 4.17 -Contact Forces in Model-24b

Figure 4 .

 4 Figure 4.18 -Crack propagation in Model-24c

  Figure 4.19 -Contact Forces in Model-24c

  Figure 4.20 -Crack propagation Model-3a

  Figure 4.21 -Contact forces Model-3a

Figure

  Figure 4.22 -Crack propagation Model-3b

Fig. 4 .

 4 Fig. 4.24) and an assumed cylindrical bone sample in pure compression (to your right in Fig. 4.24)

Figure 4 .

 4 Figure 4.25 -F vs t form contact force data

2 Figure 5 . 1 -Figure 5 . 2 -

 25152 Figure 5.1 -Single part porcine rib with SHPB at experimental scale

Figure 5 . 5 -

 55 Figure 5.4 -Dynamic variables and energy for single part porcine rib

Figure 5

 5 Figure 5.6 -Multi-part porcine rib with SHPB at experimental scale

Figure 5

 5 Figure 5.7 -Raw numerical strain signals for case 1

Figure 5

 5 Figure 5.10 -Dynamic variables and energy at interfaces for case 2

  Figure 5.12 -Consistency of force response

  Figure 5.14 -Consistency with experimental data

(

  m is the contribution of each element and is expressed asδπ m = ˆvm ρẍ i δx i dv m + ˆvm σ ij δx i, j dv m -ˆvm ρf i δx i dv m -ˆ∂κ∩∂vm t i δx i ds = 0

  Figure A.2 -8 node hexahedral element

  achieved by adding additional stiffness or additional viscosity, yet having no effects on the global stiffness in the mean time. LS Dyna offers 5 different formulation of Hourglass for solid elements, as enumerated below Type 1 Standard Viscous form Type 2 Flanagan-Belytschko integration in viscous from Type 3 Flanagan-Belytschko with exact volume integration in viscous from Type 4 Flanagan-Belytschko integration in stiffness from Type 5 Flanagan-Belytschko with exact volume integration in stiffness from

  the definition of another vector in terms of base vectors, that is known as hourglass shape vectors γ Πk as followsγ Πk = Γ Πkφ k,i 8 k=1 x k i Γ Πk (A.42)therefore, the analogue for the hourglass components of Eq. A.38 and hourglass resisting forces are rewritten as equations A.42, A.43 and A.44 are the basis for both types II and III. The only difference between type II and type III is that the shape function derivatives are evaluated at the centroid of the element for type III instead of calculating them at origin of nodal coordinate system. By doing so the anti-symmetry property of strain displacement tensor vanishes, thereby the tyep III requires few more number of evaluations for each time step. Yet the advantage of a type III hourglass formulation is that it produces an exact volume integration.

(

  represented by X in Fig. B.2).

Figure B. 2 -

 2 Figure B.2 -Surface treatment in kinematic constraint method

  Figure B.3 -Slave and Master nodes for slave search (John.O,[START_REF] Halloquist | LS-Dyna Theory Manual[END_REF]

  is the scale factor used to constrain forces in Soft based contact, m is a function of masses of slave and master nodes and ∆t c is time step for contact, that is updated at every iteration to avoid unstable contact behavior. The final contact stiffness used for contact is expressed as followsk soft = max (k cs , k i ) (B.15)NOTE the ∆t c is always set to the current solution time step value and is updated every time step increment.

  1.1 Lateral and anterior perspectives of thoracic cage . . . . . . . . . . . 1.2 Sternum and Thoracic Vertebrae . . . . . . . . . . . . . . . . . . . . 1.3 Structure of a typical human rib . . . . . . . . . . . . . . . . . . . . . 1.4 Evolution of rib cross-section . . . . . . . . . . . . . . . . . . . . . . . 1.5 Woven (to the left) and Lamellar (to the right) bone . . . . . . . . . 1.6 Cortical bone . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.7 Trabecular bone at different anatomical sites . . . . . . . . . . . . . . 1.8 Mechanical behavior of cortical and trabecular bone . . . . . . . . . . 1.9 Loading-reloading response of cortical and trabecular bone . . . . . . Initial Position of the sample . . . . . . . . . . . . . . . . . . . . . . . 2.8 Curved beam approximation of a rib . . . . . . . . . . . . . . . . . . 2.9 Normal distribution of all samples tested . . . . . . . . . . . . . . . . 2.10 Test without fracture . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.11 Experimental data for v str = 9.27 m.s -1 . . . . . . . . . . . . . . . . 2.12 Test with fracture during dynamic equilibrium . . . . . . . . . . . . . 2.13 Experimental data for v str = 17.39 m.s -1 . . . . . . . . . . . . . . . . 2.14 Test with fracture before dynamic equilibrium . . . . . . . . . . . . . 2.15 Experimental data for v str = 19.21 m.s -1 . . . . . . . . . . . . . . . . 2.16 t f rac vs v str . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Samples grouped with respect to experimental corridors . . . . . . . . 3.1 Treated porcine rib and CT scan images . . . . . . . . . . . . . . . . 3.2 Correlation of ρ app with respect to CT P R . . . . . . . . . . . . . . . . 3.3 Cross-sectional views of the rib FE model showing the density distribution. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 V f app and Po app vs CT P R . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Mechanical properties of porcine rib with ρ app . . . . . . . . . . . . . 3.6 ρ α vs ρ app of porcine rib based on

1.10

Data from

[START_REF] Granik | Human ribs: Static testing as a promising medical application[END_REF]

. . . . . . . . . . . . . . . . . . . 1.11 Experimental setup . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.12 Different stages of response . . . . . . . . . . . . . . . . . . . . . . . . 1.13 Characteristics of human rib in three point bending . . . . . . . . . . 1.14 δ ≤ 4mm (Aubert, 2012) . . . . . . . . . . . . . . . . . . . . . . . . 1.15 δ > 4mm Aubert (2012) . . . . . . . . . . . . . . . . . . . . . . . . 1.16 Compressive properties of cylindrical human rib samples . . . . . . . 1.17 Three point bending models . . . . . . . . . . . . . . . . . . . . . . . 1.18 Other Anterior-Posterior bending models . . . . . . . . . . . . . . . . 2.1 The Lagrange diagram for a typical SHPB setup . . . . . . . . . . . . 2.2 The modified three-point bending SHPB setup . . . . . . . . . . . . . 2.3 Reflective bands placed on striker and optical sensors . . . . . . . . . 2.4 P can vs v str for long striker . . . . . . . . . . . . . . . . . . . . . . . . 2.5 Strain gage for signal measures . . . . . . . . . . . . . . . . . . . . . 2.6 Signal acquisition setup . . . . . . . . . . . . . . . . . . . . . . . . . . 2.7 2.17 Energy at fracture W r vs displacement U r . . . . . . . . . . . . . . . 2.18 Energy at fracture W r vs t f rac . . . . . . . . . . . . . . . . . . . . . . 2.19 Experimental corridors . . . . . . . . . . . . . . . . . . . . . . . . . . 2.20 Experimental force displacement for Case 1 and Case . . . . . . . . . xix 2.21

Table 1 . 2

 12

	Lamellar bone on the

-Hierarchical Organization of Bone

Table 1 .3 -α based equations Group 2: ash density (ρ α )

 1

	and Currey J.D

  . The Tab.1.4.

	Author	Origin	variable	E		σ ult
	Keller (1994)	human	ρ α	10500ρ α	2.37	117ρ α	1.93
	Keyak et al. (1994) wild rabbit	ρ α	33900ρ α	2.2	137ρ α	1.88
		Table 1.4 -ρ α based equations
	Note:

Table 1 .7 -Mechanical properties as a function of apparent density as put forth by

 1 have carried out three point bending tests on rectangle sample extracted from metaphyseal shell of proximal femur and femoral diaphysis. They had considered ρ app as a control variable to study its influence on longitudinal, transversal Young's modulus (E long , E trans ) and longitudinal, transversal compressive strength (σ ult,l , σ ult,t ).

	Property	Orientation	Equation
	E, MPa	Longitudinal 2065ρ app	3.09
	σ ult,l , MPa		72.4ρ app	1.88
	E, MPa	Transverse 2314ρ app	1.57
	σ ult,t , MPa		37.0ρ app	1.51

Table 1 .

 1

	Author	Origin	Anatomical	Bone type Load	Ratio 2
			location	direction
	Reilly and Burstein (1975) human Carter et al. (1980) human	cortex proximal	cortical trabecular longitudinal longitudinal ≈ 0.7 transverse ≈ 0.5 1.0
			femur
	Stone et al. (1983) Neil et al. (1983)	bovine human	humerus lumbar	trabecular longitudinal ≈ 0.33 trabecular longitudinal 1.0
			spine
	Kaplan et al. (1985)	bovine	proximal humerus	trabecular longitudinal ≈ 0.62
	Røhl et al. (1991)	human	proximal tibia	trabecular longitudinal ≈ 1.27
	Keaveny et al. (1994)	bovine	tibia	trabecular -	≈ 0.72

10 -Ratio of tensile to compressive properties of bone

Bone Organization

Bone samples have same bone composition i.e. similar density or porosity and mineralization can show different material properties, this difference can be backed up by a difference in organization of bone tissues at different levels of Tab. 1.2. These organizational factors contribute to the anisotropic behavior of bone.

Table 1

 1

	001 0.01 0.1 1 300 1500	175.81 206.84 230.97 251.65 282.68 365.42	9.68 × 10 6 1.15 × 10 7 1.16 × 10 7 7.20 × 10 6 6.36 × 10 6 6.09 × 10 6	18.61 19.99 24.13 27.57 33.09 42.05	1.88 1.82 1.75 1.25 1.00 0.90

.11 -Average data for bovine femur bone

Table 1 .

 1 are given in Tab. 1.12 12 -Variation of tensile properties of cranial bone with ε, [s-1 ]

	Property	regression equation
	E, GPa	15.99 + 1.93Ln	ε
	σ ult , MPa	0.0826 + 6.94Ln	ε
	ε f , % e c , kg.m -3	ε 1.179 × 10 6 0.63 -0.04Ln

(Wood

Table 1 .

 1 .AuthorBone ρ, kg.m -3 E, GPa σ y , MPa ε f , % 14 -Data of existing isolated human rib FE modelThe mechanical properties of trabecular bone used by Lundin and Storm (

	Section 1.5

Table 1 .

 1 15 -Existing numerical models details for human rib in LS DynaOn the other hand,[START_REF] Pavier | Contribution à la compréhension des phénomènes physiques lors de l'impact d'un corps sur un modèle de structure biologique[END_REF] had presented a numerical model of porcine rib in dynamic three point bending to lateral medial load. The porcine rib was represented by a circular arc with no trabecular bone. He had used an elasticplastic constitutive law with no strain rate effects but with an external damage law. The properties of the numerical rib are given in Tab 1.16.

	Bone	Li et al. (2010)	Lundin and Storm (2018) Yates and Untaroiu (2018)
		Elastic-Plastic	Elastic-Plastic	Elastic-Plastic
	Cortical	Isotropic hardening Strain rates effects	Isotropic hardening -	Isotropic hardening -
		failure strain	failure strain	failure strain
		Elastic-Plastic	NR	Elastic-Plastic
	Trabecular	Kinematic hardening Strain rates effects	--	Isotropic hardening -
		failure strain	NR	failure strain

1.6 Image acquisition techniques and densitome- try Bone

density alone is not sufficient for a detailed description of bone strength. It needs several other parameters such as architectural parameters. These architectural parameters are of two categories, (i) Macro-architectural parameters such as bone size and geometry and (ii) Micro-architectural parameters such as cortical and trabecular bone thickness or structural detail.

  Si l'échantillon est suffisamment rigide pour dépasser ce seuil, il s'établit un équilibre entre la force entrante et les forces sortantes, ce qui correspond à un plateau dans la réponse. La rupture complète de l'échantillon est due à une fissure apparaissant dans la zone sollicitée en traction et qui se propage vers la zone

	comprimée (située au voisinage du point d'impact). La force maximale appliquée sur
	l'échantillon n'influence pas la rupture. On observe plutôt que la rupture dépend du
	déplacement local de la côte. french
	Abstract:

Hopkinson permettant

l'étude de la réponse en flexion. Il se compose de trois barres : une barre dite "entrante" qui va transmettre une onde de compression à l'échantillon, lequel est en appui simple sur les deux barres "sortantes" à ses deux extrémités. Les barres sont en nylon pour que leur impédance mécanique soit en accord avec celle de l'échantillon. L'analyse de l'essai et l'obtention des grandeurs intéressantes (forces aux interfaces, vitesses de déformations,...) sont basées sur la mesure des ondes de déformations élastiques sur les barres.

Au total, 20 essais ont été analysés. A partir des historiques des forces aux in-

terfaces entrantes et sortantes, trois catégories de réponses ont été identifiées : (i) essais sans rupture, (ii) essais avec rupture pendant l'équilibre dynamique entre la force entrante et les forces sortantes

1 , et (iii)

essais de rupture hors équilibre. Nous pouvons observer que la vitesse relative de l'interface joue un rôle très important dans le type de réponse avec fracture. Si cette vitesse augmente, la fracture tend à se produire avant que l'équilibre dynamique ne soit établi. De plus, une diminution du temps de rupture et du déplacement à la rupture sont observées. Un seuil critique est détecté en termes d'énergie à la rupture, de temps à la rupture et de déplace-ment à la rupture.

Table 2 .

 2 Out of 27 ribs tested, only 19 were considered to be worth for further analysis. A general normal distribution with 95% confidential bounds is used for calculating mean µ and standard deviation s.d. for each variable is given in tab. 2.1 and is also represented graphically with corresponding probability density function in Fig.2.9. 1 -Mean and standard deviation values for all properties

	• Local geometry :
	-Major axis 2a
	-Minor axis 2b
	-superior layer thickness h s
	-interior layer thickness h i

Table 3 .

 3 2 HPO 4 and H 2 O densities.

	Part Number Greyscale Intensity
	1	184.314
	2	364.706
	3	552.941
	4	741.176
	5	921.569
	6	1109.800
	7	1298.040
	8	1400.000

1 -Grey scale values for porcine rib

  V f app and Po app vs CT P R

	Part Number CT P R ρ app [kg.m -3 ] V f app	Po app [%]	Classification
	1	181.31 222.61	0.1203	87.96	Trabecular
	2	374.70 440.48	0.2381	76.18	Trabecular
	3	552.90 667.83	0.3610	63.90	Trabecular
	4	741.17 895.18	0.4839	51.51	Trabecular
	5	921.56 1113.10	0.6017	39.83	Trabecular
	6	1109.80 1340.40	0.7245	27.54	Trabecular
	7	1298.00 1567.80	0.8474	15.26	Cortical
	8	1400.00 1690.90	0.9140	8.60	Cortical
		Table 3.2 -Different constituents of Porcine rib based on ρ app
		Property	Cortical	Trabecular
		Young's Modulus E Compressive Strength σ uC 1.659 × 10 -4 ρ 1.88 app 1.109 × 10 -6 ρ 3.09 app Tensile Strength σ uT 1.61 × 10 -4 ρ 1.88 app Compressive Yield σ yC 1.422 × 10 -4 ρ 1.85 app Tensile Yield σ yT 3.064 × 10 -4 ρ 1.67 app	0.5903 ρ 1.20 app 8.534 × 10 -4 ρ 1.53 app app 8.201 × 10 -4 ρ 1.07 app 5.167 × 10 -4 ρ 1.60 7.5686 × 10 -3 ρ 1.04 app
	Table 3.3 -Power law regression equations used for mechanical properties (unit for
	stresses is here MPa)

  and 3l ′ b /4 respectively were used to record the strain signal to verify c F EM .Table 4.2 -Effect of Γ on accuracy and calculation time

	Γ	c F EM , m.s -1	Number of cycles	Time (min s)
	12 1635.98	15200	3 min 33 s
	16 1648.26	15200	4 min 48 s
	20 1659.50	15200	6 min 56 s

Table 4 .

 4 4 -Summary of simulations with single part porcine rib and rigid walls

	and that of GISSMO model

Table 4 .

 4 11 -Tensile fracture strain for single part model-24 123 As explained in Sec. 4.2.2, this damage law Every model was compared in terms of fracture propagation and contact force response. These contact force responses were then compared with experimental results.

  . Similar to Model-24b, we have effects of strain rate in this model. Moreover, this model is based on Kinematic Hardening principle instead of Isotropic. The Kinematic Hardening parameter β h delays the time of crack initiation. In all the three cases i.e. (Models-24a,24b and 3b), the magnitude of Input force is consistent.3. The response of ribs with properties fromLi et al. (2010) (Models-24c and 3a) too correspond well qualitatively with experimental responses.

	b Model-24b considers the strain rate effects through Cowper-Symmonds
	formulation. The values considered for Cowper-Symmonds parameters C
	and P do not correspond to bone material. As a reason, this has a strong
	influence in further increase in contact forces.
	c Model-3b was executed in first place for comparing the material model
	*MAT-024 and *

  ) aux cas expérimentaux. Les signaux de déformation numériques ont été enregistrés aux même points (sur les barres) que les signaux expérimentaux mesurés

	sur les barres de Hopkinson.
	Dans un premier temps, nous avons testé un échantillon numérique fictif en
	le supposant uniquement formé d'os cortical. Les résultats numériques montrent
	que cet échantillon est trop rigide par rapport aux résultats expérimentaux. Ceci
	prouve que l'os trabéculaire joue un rôle très important dans l'absorption des forces
	d'impacts.
	Deuxièmement, l'échantillon numérique "composite" (cortical + trabéculaire) a
	été testé avec les conditions expérimentales correspondant à la fracture hors équilibre.
	Les résultats de ces simulations sont conformes aux gabarits expérimentaux définis
	dans le chapitre 2.
	Les résultats numériques prouvent que le temps, l'énergie et le déplacement à la
	rupture diminuent avec l'augmentation de la vitesse donnée à l'impacteur.
	Si on le compare au cas expérimental retenu, le modèle numérique proposé
	montre une grande efficacité dans la prévision du déplacement, des forces aux
	interfaces, des seuils d'énergie et des fractures.

represented as ρ app in this thesis.

critical velocity occurs when a material property exhibits large variations over a small range of ε,[START_REF] Mcelhaney | Dynamic response of bone and muscle tissue[END_REF]

L'équilibre dynamique est établi lorsque F in ≈ F out1 + F out2

ACKNOWLEDGEMENTS

I would like to sincerely thank my PhD supervisors Dr. André LANGLET and Pr.

Ridha HAMBLI for their valuable suggestions, insights and for having supported me during the time period of my thesis. I would also like to express my sincere gratitude to Pr. Patrick CHABRAND of l'Institut des Sciences du Mouvement from the Aix-Marseille Université and to Pr. Sébastien LAPORTE of l'Institut de Biomécanique Humaine Georges Charpak

Bibliography

We can conclude that the ash content in bone has a strong positive correlation on strength of bone and a negative correlation on the initial elastic behavior of bone. The reason for this influence in particular for human ribs or its biological surrogates is yet to be explained.

Unfortunately, there is a lack of data on the influence of α, ρ α or two variable models on σ y . Again, this difference in V f is due to the consideration of α in model 2.

Author

Section 4.2 and to ensure accurate prediction of the fracture path growth [START_REF] Hambli | Micro-CT finite element model and experimental validation of trabecular bone damage and fracture[END_REF].

Therefore, an element size of about 0.15 mm was retained to mesh the ribs leading to a total of 226520 eight nodes hexahedral (solid brick) elements.

Constitutive and Damage Law for bone

Several authors have presented different constitutive laws for cortical bone tissue.

However, the damage accumulation and the considerations its accumulation into numerical models have been a point of arguments among these authors. We can find many articles using an elastic-plastic constitutive law, for example [START_REF] Garcia | Elastic plastic damage laws for cortical bone[END_REF], [START_REF] Pavier | Contribution à la compréhension des phénomènes physiques lors de l'impact d'un corps sur un modèle de structure biologique[END_REF]. A realistic constitutive law for bone must consider the effects of strain rate. Regarding the effects of strain rates, quite a few authors have suggested either viscoelastic-viscoplastic constitutive model [START_REF] Johnson | A viscoelastic, viscoplastic model of cortical bone valid at low and high strain rates[END_REF] or elastic-viscoplastic constitutive model [START_REF] Li | Rib fractures under anterior-posterior dynamic loads: Experimental and finite-element study[END_REF].

Abstract:

This chapter presents an application of the numerical model established in the preceding chapter at experimental scale. In order to ensure similar conditions of experimental setup, numerical strain signals were extracted from elements at same distances from the sample interfaces.

First of all, a numerical sample was assumed to be made of only cortical bone was tested (namely single part). The results show that the single part sample is too rigid with respect to experimental cases. This proved that the trabecular bone plays a crucial role in absorbing the impact forces of a rib.

Secondly, a numerical sample with both cortical and trabecular was tested (namely multi part sample). Experimental cases corresponding to fracture before dynamic equilibrium were the center of interest. Therefore, striker velocities of corresponding cases were simulated. The results of multi part sample are in accordance with the established experimental corridors.

The numerical results show that as striker velocity increases, the time , energy and displacement at fracture decreases. The numerical threshold values of time, energy and displacement at fractures are similar to that of experimental threshold values.

The proposed FE model shows high efficiency in predicting force-displacement, energy thresholds and similar fracture patters to that of experimental cases.

Comparison with experimental data

The relative interface velocity V for this case varies around 11.14 [m.s -1 ]. Therefore, the results were compared to the experimental corridor 3.

Comparison of numerical and experimental energies

As we have already established the experimental energy at rupture W r in sec. 2.5.2, we can now compare the consistency of numerical data. Since the aim of numerical simulations is to have insights on fracture before equilibrium, The numerical data were compared to corresponding experimental values.

Appendices

Appendix A

Solid Element Formulation Preliminary

Let us consider a body denoted by κ in a Cartesian coordinate system initially at a point represented by X α , (α = 1, 2, 3) which displaces to a point x i , (i = 1, 2, 3) that is denoted by κ in the same coordinate system with time t as in The time dependent displacement of this body using the Lagrangian approach can be written in terms of its initial coordinates as follows:

Thus, the initial condition at time t = 0 can be expressed as

where, v i is the initial velocity.

Similar to the initial conditions, the boundary condition on the bodies, such as traction boundary condition of the system, displacement boundary condition on the boundary κ and the contact discontinuity boundary condition on the boundary κ can be expressed as in Eq. A.4, Eq. A.5 and Eq. A.6 respectively.

here, σ ij is the Cauchy's stress tensor, n ij is unit normal of an element on the boundary of body configuration κ and D i is the displacement matrix of the body κ. The notations + and -represent the quantities exterior and interior to the boundary of κ.

Having established the initial and boundary conditions, the governing equations for the system in Fig.

A.1 are written as follows:

Conservation of mass:

where, ρ 0 is the reference or initial density and V is the relative volume and is expressed in terms of determinant of deformation gradient matrix F ij .

Conservation of momentum:

where σ ij, j represents the co-variant differential of stress tensor and ẍi represents the acceleration. The solution of momentum equation Eq. A.8 in its weak form must satisfy the boundary conditions in Eq. A.4, Eq. A.5 of κ and Eq. A.6 on κ respectively.

Conservation of energy:

the variables s ij and p represent the deviatoric stress tensor and pressure respectively and are expressed as follows

here, δ ij represents the Kronecker delta and q represents the bulk viscosity.

The displacement δx i , that satisfies all the boundary conditions on κ can be written using the divergence theorem as:

as we know that

The Eq. A.12 can be rewritten as

This leads to the statement of virtual work in its weak form as in Eq. A.15

Principle of virtual work

Using the principle of virtual work, the change in external work done in a system can be expressed in an integral form as follows:

Having established the principle of virtual work for our system, let us consider that the system is discretized into n elements with k nodal points. Then, the relation representing the displacement of these nodal point is written as

here, φ j represents the shape or interpolation functions in natural coordinate system (ξ, η, ζ) and x i j represents the nodal coordinate of j th node in i th direction.

Thus, all the terms of strain displacement matrix can be calculated by the Jacobian J

Volume integration

The volume integration in LS-Dyna is executed through Gaussian quadrature method. For this case, let us assume a function 'g' defined over a volume with N number of integral points, then the integral of g over the volume is:

The integrals in Eq. A.26, can be approximated as summation over integral points in terms of weighting functions w j , w k , w l as follows

These fully integral elements have the following disadvantages:

• The computation effort required to compute the matrices is huge and is expensive in time. For example, (a) time required for computing the strain matrix is around 25 times more than that required for a one point integration.

(b) the number of multiplies for strain and nodal force matrices is 16 times more than that for one point integration.

(c) even the constitutive evaluation needed for a fully integrated element is substantially more.

• in addition to cost, the other common inconvenience with fully integrated elements is element locking whist used in plasticity problems or problems involving failure.

These disadvantages of fully integrated elements can be overcome using one point integration solid elements. Moreover, it is not new to use one point integration el-Belytschko, 1981)

Similarly, the displacement field and velocity fields can be expressed for a hexahedral element with one point integration can be expressed as follows

The velocity gradient with respect to the direction and nodal coordinate can be expressed as ui,j = uik φ k,j (A.30)

where D ij is the velocity strain tensor that is symmetric and W ij is the instantaneous rotation tensor that is anti-symmetric. On the other hand, the subscript (i, j), represents the differentiation of the quantity with respect to nodal coordinate.

The principle of virtual work can be expressed in terms of nodal force as :

In fact the one point Gaussian integral neglects the nonlinear terms, thereby considering a uniform stress-strain state through out the element. Thereby, leading to the velocity gradient in the following form .34) integration in terms of hourglassing modes can be expressed as 2, 3, 4 (A.37) where i denotes the directions and Π denotes the hourglassing modes. When the nodal velocities have no hourglass modes, the product of base vectors is zero. Therefore, the hourglass component can be written as

Hence, the hourglass resisting forces f k iΠ can be expressed as It is evident that the base vectors of Type I are not orthogonal to the linear velocity fields. Therefore, LS Dyna offers another formulation, in which the base vectors Γ Πk are orthogonal to the linear vector filed by adapting an other manipulation proposed by [START_REF] Flanagan | A uniform strain hexahedron and quadrilateral with orthogonal hourglass control[END_REF], which is referred to as Type II.

Types II and III

The hourglass type II is based on the hypothesis of isolating the linear modes from the total velocity or displacement field gives rise to hourglass fields. This is mathe-

Appendix B Contact Algorithms

To understand the function of contact algorithms, let us consider two bodies B 1 and interact according to the principles of mechanics of continuous media".

This comparison of different contact algorithms have led us to the conclusion that the Penalty Constraint Method (PCM), has an upper hand over all three and is chosen as potential candidate for the contact analysis in this thesis.

Slave Search

The In case these two nodes do not coincide, then n s is assumed to be in any of the segments (for example say S 1 ). This can be verified by using the following equations

here, c i and c i+1 are vectors along the adjacent edges of S 1 pointing outward from m s . s is the projection of vector g on the segment S 1 (see Fig. B.3b), that is expressed as

Logically speaking, the contact surfaces are not always uniformly meshed. In case the contact surfaces contain badly shaped elements, it is not easy to locate the slave nodes as there is a probability that the salve node n s is not located in the identified master segment (as in L'objectif de cette thèse est de caractériser le comportement mécanique des côtes isolées sous impact en dynamique rapide par l'expérimentation et les calculs EF. Dans le cadre de cette étude, les côtes de porc on été choisies comme substituts de la côte humaine. Les essais de flexion trois points en dynamique rapide on été réalisés avec les barres de Hopkinson. Les résultats de campagne expérimentale ont mis en évidence l'influence de la vitesse de déformation sur le temps de fracture, le faciès de rupture et également, sur l'amplitude des forces d'impact. Un modèle numérique 3D d'une côte de porc a été développé, grâce aux images obtenues par micro-CT en Haute Résolution (HR-pQCT). Ensuite, les propriétés élastiques des constituants on été attribués à partir des niveaux de gris des images de HR-pQCT. Une loi de comportement élastique-viscoplastique, modifiée pour prendre en compte les effets de la vitesse de déformation couplée avec une loi d'endommagement a été élaborée. La loi d'endommagement présentée, nous permet de prendre en compte l'endommagement non-linaire, la réduction de la rigidité de la structure et les effets de la vitesse de déformation. Le modèle EF présenté dans cette thèse est capable de prédire la fracture des côtes sous impact en dynamique rapide.

Mots clés : Chocs, Impactes, Réponse dynamique, côte de porc, HR-pQCT, simulations EF, LS-Dyna, SHPB

Characterization of rib fractures subjected to impact loads: numerical simulations and experimental validation

Abstract :

The objective of this study is to characterize the mechanical behavior of isolated porcine ribs subjected to dynamic impact through experiments and numerical simulations. The porcine rib was used as a human surrogate rib in order to develop an adapted methodology for the human rib. A three Split Hopkinson Pressure Bar (SHPB) setup for three-point bending tests was used. An ensemble of 20 tests data was considered to be comprehensible for experimental characterization, thereby, showing an influence of strain rate on both time for fracture and amplitudes of force response. A three-dimensional porcine rib model was generated from the DICOM images of High-Resolution peripheral Quantitative Computed Tomography (HR-pQCT) scans. Material properties interpolated using power law regression equations based on apparent density were assigned to the numerical rib. A modified elastic-viscoplastic constitutive law, capable of considering the effects of strain rate was elaborated. An incremental and stress-state dependent damage law, capable of considering effects of strain rate on fracture propagation, non-linear damage accumulation and instabilities was coupled to the constitutive law. The proposed Finite Element model is able to predict satisfactory force-displacement curve and fracture patterns of the tested ribs indicating that the developed numerical model may be used to investigate the fracture behavior of human ribs under dynamic loads.