Alleau Thibaut

Mme Irène

Examinateur M Emmanuel Lefrançois

M Gaël

M Miguel

Rapporteur) Mme A Fernandez

Anne-Virginie Salsac

Directrice De Thèse

) M Laurent Lanquetin

M Eric

Bergoend M Gaetano Burriesci

DEVELOPMENT OF A NUMERICAL PLATFORM TO MODEL THE MITRAL VALVE

Development of a numerical platform to model the mitral valve

Mitral insufficiency is the first valvular disease worldwide, with a 2% prevalence. When open-heart surgery is impossible for the patient, surgeons use percutaneous devices to help the mitral leaflets coapt. However, the only device currently available is based on the edge-to-edge mitral valve repair technique. This type of implant is not adapted for patients suffering from functional mitral insufficiency, where the ventricle is responsible for the lack of coaptation of the leaflets. This thesis aims to provide a numerical platform to help the development of a mitral valve implant adapted for those patients.

Several mitral valve geometries were created from a parametric model using anatomical measurements. Finite element simulations of the mitral valve were performed using ADINA to determine the valve closure under constant pressure. Several material models were developed in large strain and large deformation to model the valve closure accurately. Pathological behaviour such as annulus dilatation and chordae rupture were modelled, and several methods were tested to implement medical devices. Fluid-structure interaction of a 2D mitral valve was obtained using an ALE description and a monolithic coupling approach. Both the systole and the diastole were reproduced and studied, and the hermetic seal of the valve was detailed.

The numerical platform developed is suited to model mitral valve function and can be used to help the development of mitral implants. In addition, the parametric geometry model and the anisotropic material model will be useful to depict with realism the valve function. A 3D fluid-structure interaction of the mitral valve could be developed.

Annexes

6.1 Summary of hyperelastic theory . . . . . . . .

Résumé

Développement d'une plateforme numérique pour modéliser la valve mitrale L'insuffisance mitrale est la valvulopathie mondiale la plus fréquente avec une prévalence de 2%. Lorsque le patient n'est pas en mesure d'être opéré à coeur ouvert, un implant percutané est utilisé pour aider la fermeture des feuillets. Le seul implant actuellement disponible est basé sur la réparation bord à bord de la valve mitrale. Il réduit le reflux vers l'oreillette lors de la systole, mais n'est pas adapté pour les patients souffrant d'insuffisance mitrale fonctionnelle, chez qui la pathologie provient du ventricule et non des éléments de la valve. L'objectif de la thèse est de fournir une plateforme numérique permettant d'aider au développement d'un implant adapté pour ces patients. Plusieurs géométries de valve ont été réalisées au moyen d'un modèle paramétrique, en utilisant des données anatomiques. La dynamique de la valve a été modélisée avec le logiciel ADINA par des simulations éléments finis en grandes déformations. Des modèles structurels de la valve ont permis de représenter la fermeture de la valve sous une pression uniforme. Les lois de comportement de matériaux ont été développé dans le but d'obtenir une fermeture réaliste de la valve. Cela a nécessité la prise en compte de l'hyperélasticité et de l'anisotropie des tissus. Des pathologies valvulaires, telles que la dilatation de l'anneau mitrale ou la rupture des cordages tendineux ont été modélisées, et plusieurs méthodes ont été testées pour y apposer des systèmes médicaux. En utilisant une description ALE et un couplage monolithique, les interactions fluide-structure ont été simulées pour une valve mitrale bi-dimensionnelle. La fermeture hermétique de la valve pendant la systole a pu être reproduite et l'ouverture de la valve étudiée pendant la diastole.

La plateforme numérique développée permet de modéliser la fonction de la valve mitrale et peut être utilisée pour aider au développement d'un implant mitral grâce au modèle paramétrique reproduisant différentes géométries de valve et aux lois matériaux anisotropes. Une perspective reste la création d'un modèle 3D des interactions fluide-structure de la valve mitrale.

Acknowledgments

My PhD thesis would not have been possible without the energy and willpower of both my supervisors, Anne-Virginie Salsac and Laurent Lanquetin. For four years, they believed in me, pushed me forward and helped me reach our objectives. I want to thank you for your trust and your moral support.

I want to thank Gaël Rochefort and Miguel A. Fernández for their review of my work and the correction they provided.

I want to thank Gaetano Burriesci for his overview of my work and the critical points he detailed. He provided helpful insights on mitral valve simulation, particularly the differences between modelled and physiological valves.

I want to thank Emmanuel LeFrançois sincerely for the valuable insight he gave me, especially for discussing the added mass in fluid-structure interaction.

I want to thank Eric Bergoënd, a thoracic surgeon at APHP (Paris Hospital). He always believed in my capacity for this project and provided medical information for my thesis. I would also like to thank Jean-François Deux, a radiologist at APHP, for the medical images he gave me. Finally, I would like to thank Antonio Fiore, who showed me an open-heart surgery.

I want to thank Mme. Irène Vignon-Clémentel for presiding over my PhD defence and her insight on 0D models and their application to my simulation. The models developed could greatly be improved with 0D models and I look forward to implementing them.

I want to thank François-Xavier Gerbeaux, with whom it's always a pleasure to discuss this project and science in general. He is the one who cited Nicolas-Boileau, and it helps me a lot when writing this thesis. Thank You. "Ce que l'on conçoit bien s'énonce clairement, Et les mots pour le dire arrivent aisément"; L'Art Poétique, Nicolas Boileau-Despréaux, 1674.

I want to thank the modelling team at Segula Matra. M. Caudal supported me and provided logistical organization. M. Anoica Vlad for his moral support and M. Lamar Romuald for all the structural and non-linear modelling conversation.

I want to thank the BMBI laboratory, where I did part of my research and PhD, and all the friends I made during the years. Nicolas, Doriane, Baptiste, Xing Yi; Thanks a lot for the time spent.

I want to thank the Research and development team at ADINA and M. Veal Thirom for their support and help during my PhD. In particular, M. Veal Thirom shared his passion for the software ADINA and did his best to help us during these years.

To my family, friends from the engineering school, who believed in me, each day, each time we called, without a doubt. Your unwavering moral support pushed me forward and helped me get this done.

Finally, to my fiancé. She was there when doubt came, when my deception reach summits, when I gave it up. She made me power through, keep it up, continue. I am so grateful she did, and thanks her from the bottom of my heart.

This work was funded by "Association National de la Recherche et de la Technologie" grant CIFRE 2017/1188, Segula Technologies and "Université de Technologie de Compiègne". I sincerely thank all the financial collaborators who made feasible this research.

Introduction

Mitral valve repair and replacement are the second most performed heart valve procedures after aortic valve replacement. Their operative mortality is 1.5% and 5%, respectively, the one of aortic valve replacement being 3% in comparison [1]. The much lower mortality rate compared to mitral valve surgery motivates the recent studies in mitral repair or replacement strategies. Transcatheter valve replacement was a significant breakthrough in cardiac surgery, reducing mortality rate and possible complications. Initially created for the aortic valve, efforts have been made to adapt it to the mitral valve. The mitral valve is more complex than the aortic valve due to its dynamic motion, the possible obstruction of the ventricular aortic output and the diversity in mitral valve geometry between patients. Such difficulties delay the appearance of transcatheter mitral valve replacement devices on the market. Repair procedures have lower mortality rates and complications than the replacement of the mitral valve. They are the favoured pathways for mini-invasive treatment options. A new mini-invasive device must have better performance than the Cardioband or the Mitraclip/Pascal, currently on the market. This PhD thesis, conducted in Segula Technologies (Matra) and UTC University in Compiègne, aims at creating structural and fluid-structural numerical models of the mitral valve motion. The numerical models are to assess healthy and pathological mitral valve behaviours. The ultimate goal is to test the impact of mitral valve repair implants in pathological cases.

The first objective is to obtain a geometrical model of the mitral valve. Ideally, this model could be fitted to patient-specific data to increase its realism. The second objective of this thesis is to identify and represent the mechanical response of the valve's constituents. The third objective is to represent the dynamic and static response of the valve under systolic pressure. The valve shape under systolic loading could be used to estimate the valve realism and closure effectiveness. The fourth and last objective of the thesis is to develop a fluid-structure methodology to model the mitral response to the fluid during cardiac cycles.

Modelling the mitral valve raises several difficulties. First, a numerical representation of the mitral apparatus must be obtained. Numerical representation of human parts can generally be extracted from in-vivo scanning techniques. Nonetheless, the mitral constituents are difficult to observed on 2D in-vivo images, and the 3D reconstruction of the mitral valve requires specific algorithms and tools. Therefore, another approach to reconstruct the mitral valve is to use a geometrical model fitted with anatomical data from in-vivo images or ex-vivo measurements. The first chapter of the thesis discusses this issue and the solution used to obtain a 3D geometry of the mitral valve.

The second difficulty of mitral valve modelling resides in the mechanical response of the mitral tissues. They are mainly composed of fibres surrounded by a protein matrix. Thus, their mechanical responses are hyperelastic and anisotropic. In addition, the mitral valve underwent large deformation during the cardiac cycle. As a result, the tissue response to the deformation requires specific attention.

Hence, the second chapter is dedicated to the mechanical response of the mitral valve.

The third chapter is devoted to the structural simulation of the mitral valve. It details the methods used for the contact, the time scheme, the elements used and the technique applied to fix a hypothetical device to the valve's leaflets. Finally, the mitral valve responses are given and discussed.

The last difficulty of mitral modelling is its interaction with the surrounding blood. The large displacement of the leaflets in the ventricle requires specific attention. In addition, the mitral valve seals the ventricle during systole, which can be challenging to model. Besides, the coupling of the fluid and the structure results in strong added mass and can cause numerical instabilities. The fourth chapter is dedicated to these issues. A simulation of a 2D mitral valve in a ventricle under cyclic loading is used to demonstrate the method used.

Chapter 1

Mitral valve model

To model the valve function in healthy and pathological cases, a representation of the valve is required. This chapter explains how the numerical image of the valve was obtained for the structural simulations. The first section describes the mitral valve and its components, with common patterns and variations. The second section details the methods used in the literature to obtain a numerical representation of the valve. The third section is devoted to the methods used to construct the valve, and the fourth section details the numerical models obtained with the previously detailed methods. Possible improvements and current defaults of the model conclude the chapter.

Anatomy of the mitral valve

This section aims to define what is the mitral valve and to detail anatomical variations. A generic form of the mitral valve is first presented, and variations are then depicted from scans or ex-vivo data.

General definition

The mitral valve separates the left atrium from the left ventricle. It guides the flow of blood from the atrium to the ventricle. It can be opened or closed, depending on the blood movement and the heart contractions. The mitral valve is composed of four key structural components (figure 2.1 & 1.2). The key components are the mitral annulus, the posterior and anterior leaflets, the chordae tendineae, and the papillary muscles.

The annulus is a ring of fibrous tissue that circumscribes the valve perimeter. It has a dynamic saddle shape and is anchored on the fibrous trigones, semi-rigid parts of the heart that is also used by the aortic valve [3]. It shares a fibrous continuity with the aortic valve [4], [START_REF] Veronesi | Quantification of Mitral Apparatus Dynamics in Functional and Ischemic Mitral Regurgitation Using Real-time 3-Dimensional Echocardiography[END_REF]. The right and left trigones separate the valve into anterior and posterior parts, the posteriors being distal to the trigones. The shape of the annulus evolves during the cardiac cycle. During diastole, the annulus is circular and nearly flat while during systole, it has a saddle shape. The annulus provides a fibrous skeleton for the leaflets to anchor into.

Leaflets are a continued band of tissue that extends from the annulus. The anterior leaflet, closer to the aortic valve, is larger than the posterior leaflet, closer to the pericardium [3]. The leaflet surfaces are relatively smooth on the atrium side of the heart. They are more irregular on the ventricular side due to the chordae tendineae implementation. This irregular zone is called the rough zone of the leaflets. Two areas exist, one is close to the free edge of the leaflets and the second is beneath the annulus, in the basal region (figure 1.2). When the mitral valve closes, the leaflets are maintained into the ventricle by the chordae tendineae. The leaflets coapt under the systolic pressure, sealing the [2] gap between the atrium and the ventricle, and preventing blood from crossing the annulus. During diastole, the leaflets move away and the valve opens [4], [6].

Chords, or chordae tendineae, arise from the two papillary muscles and insert into the leaflets in a tree-like manner [START_REF] Lam | Morphology of the human mitral valve. I. Chordae tendineae: a new classification[END_REF], [START_REF] Al-Atabi | Biomechanical assessment of surgical repair of the mitral valve[END_REF]. They are classified depending on their size and attachment area. Marginal chords are the thinnest and insert in the free edge area of the leaflets. Basal chords are thicker and insert in the rough basal zone. Strut chordae arise from the papillary muscle tips and insert in the rough zone of the anterior leaflet (figure 1.2). They are the largest and thickest of all chordae. The different types of chords have different functions. Studying ovine and porcine mitral valves, it was found that marginal chords are primarily responsible for the leaflet's closure during systole, that basal chordae supports and transfer loads to the leaflets, and so protect marginal chords from failure [START_REF] Dagum | Three-dimensional geometric comparison of partial and complete flexible mitral annuloplasty rings[END_REF]- [START_REF] Wilcox | Frequency and diameter dependent viscoelastic properties of mitral valve chordae tendineae[END_REF].

The papillary muscles are muscular extensions of the left ventricular wall connected to the chordae tendineae. Two muscles are usually distinguished. One is anterolateral and the second is posteromedial. They are located underneath the commissures, where the anterior and posterior leaflets merge. Papillary muscles are dynamic structures, alike the annulus, and move through the cardiac cycle toward and from the annulus plane. Their role is not clearly defined in the literature. Papillary muscles contraction would either help opening or closing the mitral valve [START_REF] Oliveira | Geometric description for the anatomy of the mitral valve: A review[END_REF], or even act as shock absorbers during the cardiac cycle [START_REF] Rabbah | Mechanics of healthy and functionally diseased mitral valves: A critical review[END_REF].

Variations and patterns of the mitral valve

The mitral valve differs among individuals. The valve size depends on the size of the ventricle. Its volume is estimated from the individual height. 1 m of a healthy person is equivalent to 32 mL of the left ventricle. Thus, men have on average 58 mL and women have 52 mL. Thus, non-pathological variability is considered to be less than 20%, according to cardiothoracic surgeons. Circular approximation [START_REF] Sakai | Distance between mitral anulus and papillary muscles: Anatomic study in normal human hearts[END_REF] Table 1.1: Size and shape of the mitral valve annulus. * Values approximated from parameters.

Annulus

The annulus can be referred to in size using the anteroposterior length, the inter-commissural length, and the height of the anterior/posterior saddle horn (figure 1.3). The annulus is nearly flat during diastole and has a saddle shape during systole. When flat, the annulus can be described only with a planar shape and the inter-commissural and anteroposterior lengths. Diameter and perimeters can also be used to describe the annulus. Data of the annulus size are given table 1.1.

Leaflets

The leaflets are a continuous band of tissue that extends from the annulus toward the ventricle. Each point of the annulus is linked to a specific point of the leaflet's free edge. The distance between the annulus and the leaflet free edge can thus be monitored using the annular length. The maximal size of the leaflets, their surfaces, and the angle of the leaflets towards the annulus plane are used to describe the leaflets' shapes and sizes (table 1.2a and 1.2b with figure 1.4). The thickness of the leaflet varies on its surface, depending on the anatomical layers and the chordae attachment location. Some anatomical measurements of the thickness and values used in finite element models are given in table 1.2c. Among them, values reported by Sacks and co-workers were obtained using a µ-CT scan from the ovine's mitral valve and are the most accurate [START_REF] Khalighi | Multi-resolution geometric modeling of the mitral heart valve leaflets[END_REF]. 11.4 -9.9 15.6 -12.0 11.4 -9.9 9.0 -6.1 9.0 -6.1 [START_REF] Lau | Numerical simulation of mitral valve function[END_REF] (a) Noticeable leaflets heights (in mm). References explained in figure 1. 4.

leaflets angles (°, ant/post) source 24.6 ± 7 / 34.7 ± 9.6

[14] 55 / 70

[21] 25 / 19

[18] 8.3 / 5.3 [START_REF] Leonardo | Mitral Valve Finite-Element Modelling from Ultrasound Data: a Pilot Study for a New Approach to Understand Mitral Function and Clinical Scenarios[END_REF] (b) Leaflets angle toward annulus plane.

leaflets thickness (mm) source 1.1

[21] 0.75 (iso), 1 (hyper & ortho) [16] 0.5 -3 (non-uniform,scan)

[18] 1 -1.7 (non-uniform)

[23] 0.8 -1.1 (non-uniform)

[24] 1.32 / 1.26 (ant/post) [START_REF] Leonardo | Mitral Valve Finite-Element Modelling from Ultrasound Data: a Pilot Study for a New Approach to Understand Mitral Function and Clinical Scenarios[END_REF] (c) Leaflet thickness used in mitral valve model.

Table 1.3: Heights, angles and thickness variations of the mitral valve

Papillary muscles

Two papillary muscles are commonly described in the human mitral valve. However, one anatomical study shows a large variation of papillary muscles numbers [START_REF] Gunnal | Morphological Variations of Papillary Muscles in the Mitral Valve Complex in Human Cadaveric Hearts[END_REF]. According to the study, only 3% of 116 human valves studied were composed of 2 papillary muscles, 43% of the papillary muscles studies were separated into two groups, 32% were separated into 3 groups, 22% were separated in 4 groups. The papillary muscles actively contract at systole but their tip has a near-constant distance with the annulus during the cardiac cycle [START_REF] Rabbah | Mechanics of healthy and functionally diseased mitral valves: A critical review[END_REF]. Based on this observation, the displacement of the papillary muscle tips can be reconstructed from the movement of the annulus (see [START_REF] Leonardo | Mitral Valve Finite-Element Modelling from Ultrasound Data: a Pilot Study for a New Approach to Understand Mitral Function and Clinical Scenarios[END_REF] for an example using 3D-Tee+t). Anatomical ex-vivo measurements of the PM location on 57 human hearts shows that papillary muscles can be considered symmetrical (figure 1.5). A recent review of the literature lead to the same observation [START_REF] Oliveira | Geometric description for the anatomy of the mitral valve: A review[END_REF] (figure 1.6). The group of Burriesci used the same hypothesis and put the PM tips at 20 -22 mm above the annulus plane for their simulations based on the data table 1.5b, [START_REF] Lau | Numerical simulation of mitral valve function[END_REF]. (a) Representation of the measures used to locate the papillary muscles [START_REF] Oliveira | Geometric description for the anatomy of the mitral valve: A review[END_REF]. Measures of papillary muscles distances toward mitral valve markers [START_REF] Oliveira | Geometric description for the anatomy of the mitral valve: A review[END_REF].

Chordae tendineae

The chordae tendineae are a tree-like structure that connects the papillary muscles to the leaflets. As they divide along the branches, the chord cross-sectional area (CSA) reduces. The number of chords, their branching type, their CSA, and the insertion sites are sources of morphological variability [START_REF] Oliveira | Geometric description for the anatomy of the mitral valve: A review[END_REF]. The high variation of the branching pattern was emphasized in anatomical measurements of 116 [START_REF] Gunnal | Morphological Study of Chordae Tendinae in Human Cadaveric Hearts[END_REF]. f obs = N observed /N valves . The fan-shaped was always observed, mostly for the marginal chordae insertion.

Chordal category [ human mitral valves [START_REF] Gunnal | Morphological Study of Chordae Tendinae in Human Cadaveric Hearts[END_REF] (figure 1.7).

The group of Sacks studied the mitral branching pattern using a generic method [START_REF] Khalighi | Mitral Valve Chordae Tendineae: Topological and Geometrical Characterization[END_REF]. Using µ-CT scan on ex-vivo valves, 3D segmented images were obtained, transformed into skeleton graphs, and then simplified into an isoparametric B-spline model. The branching pattern of the chordae was studied using Reeb graphs of bifurcation and non-planarity angles. Insertion locations of chordae tendineae were also measured. The complex mathematical models were extracted from the data and used to generate realistic mitral valve chordae [START_REF] Khalighi | Development of a Functionally Equivalent Model of the Mitral Valve Chordae Tendineae Through Topology Optimization[END_REF], [START_REF] Sacks | On the Simulation of Mitral Valve Function in Health, Disease, and Treatment[END_REF]. Variations of the cross-sectional area are given in table 1.4 and can be used to generate realistic chords. For finite element simulations, uniform cross-section area of 0.6 mm 2 are used [START_REF] Lau | Numerical simulation of mitral valve function[END_REF], [START_REF] Kunzelman | Fluid-structure interaction models of the mitral valve: function in normal and pathological states[END_REF].

The number of chordae tendineae is also a source of variability. The number of marginal (basal) chordae tendineae is 22.6 ± 7.0 (10.3 ± 2.6 respectively). In general, there are 2 to 4 strut chordae on the mitral valve [START_REF] Gunnal | Morphological Study of Chordae Tendinae in Human Cadaveric Hearts[END_REF]. Finite element model of the mitral valve often contains more chordae tendineae (around 60 chords with 20 basal and 20 marginal) [START_REF] Kunzelman | Fluid-structure interaction models of the mitral valve: function in normal and pathological states[END_REF], [START_REF] Lau | Fluid-structure interaction study of the edge-to-edge repair technique on the mitral valve[END_REF]. A large number of chords aim to model more accurately the function of the chords while neglecting some aspect of the branching pattern.

Mitral valve numerical representation techniques

To model the mitral valve movement, a numerical representation of the valve is required. Two main approaches are described in the literature to obtain a 3D shape of the mitral valve. The first approach consists of imaging the mitral valve using a scanning technique. The second method is to define the mitral valve shape using a parametric model and the data described in section 1.1.2.

Valve geometry obtained from imaging techniques

Several scanning tools are already available in the medical field, and it is possible to obtain a 3D geometric representation of body parts using such tools. Bones and blood vessels are detected in medical routines which often provide accurate 3D models. The 3D geometric representation can then be transformed, cleaned, and used in finite element simulation. The mitral valve can also be detected using imaging techniques. Computed Tomography Scan (CT-scan), Magnetic Resonance Imagery (MRI), and 3D echocardiography transesophageal (3D-TEE) can be used to provide information about the valve function. However, none of these in-vivo medical scanning routines provides a 3D model of the valve that is accurate enough for its uses for finite element modelling.

CT scans

The CT scan image represents several levels of X-rays absorption using a greyscale, which is directly dependent on the tissue. The image of the tissue, obtained in 3 dimensions, is composed of pixels. Imaging software can reconstruct 3D shapes from such images, and extract the geometry in a universal format suited for finite element (IGES). Bones, blood vessels can be obtained using these scans. The method was applied to a coronary CT scan provided by Dr Deux, an experienced vascular radiologist at Henri Mondor Hospital. However, no clear mitral valve tissue delimitation could be obtained. The internal heart wall, the annulus, and the papillary muscles were easily delimited, figure 1.8. The mitral leaflets, the tendineae chordae, and the blood could not be separated. Excluding the blood grey level from the filter results in a torn leaflet and residual blood bubbles (figure 1.8). It seems that the grey level of the leaflet tissue is the same as the grey level of the blood. Besides, the leaflets and the chordae tendineae of the mitral valve are thin structures, about 1 mm and 0.5 mm respectively. This is similar to the precision of the CT scan 0.625 mm, which increases the difficulty of detection. To the author's knowledge, only one group has successfully reconstructed a mitral valve leaflet from computed tomography medical images (figure 1.9a) [START_REF] Wang | Finite Element Modeling of Mitral Valve Dynamic Deformation Using Patient-Specific Multi-Slices Computed Tomography Scans[END_REF], [START_REF] Mao | Fully-Coupled Fluid-Structure Interaction Simulation of the Aortic and Mitral Valves in a Realistic 3D Left Ventricle Model[END_REF]. Although the technique was reused on a tricuspid valve, it is not clear whether it is generally sharp enough to be generalized on any patient valve and be available as a routine [START_REF] Kong | Finite Element Analysis of Tricuspid Valve Deformation from Multi-Slice Computed Tomography Images[END_REF].

MRI

MRI techniques are less frequently used than CT-scan for cardiac imaging but provide more contrasted images. The resolution of MRI varies with the magnetic field used but seems to be limited to 0.5x0.5 mm 2 with slices of 3 mm. A 3D model is reconstructed from the slices. Such scans were used during a cardiac cycle to image the mitral valve leaflets [START_REF] Ma | Image-based fluid-structure interaction model of the human mitral valve[END_REF]- [START_REF] Gao | A Coupled Mitral Valve -Left Ventricle Model with Fluid-Structure Interaction[END_REF]. Points of the interface between the blood and the valve were detected (figure 1.9b). The points were then used for the mitral reconstruc-(a) Multi-slice CT scans [START_REF] Wang | Finite Element Modeling of Mitral Valve Dynamic Deformation Using Patient-Specific Multi-Slices Computed Tomography Scans[END_REF] (b) Leaflets interfaces from MRI scan [START_REF] Ma | Image-based fluid-structure interaction model of the human mitral valve[END_REF] (c) Fibre reconstruction from µ-CT scan [START_REF] Toma | High-Resolution Subject-Specific Mitral Valve Imaging and Modeling: Experimental and Computational Methods[END_REF] Figure 1.9: Imaging techniques used to model the mitral valve tion with a geometric model (using Gambit & SolidWorks). No information on the mitral subvalvular apparatus was obtained.

Ex-vivo scans

Part of the difficulty to scan the mitral valve is linked to in-vivo scanning. Valves from animals or cadavers can be imaged ex-vivo to increase both contrast and resolution. µ-CT scans have an isotropic resolution of 40 µm, and a highly contrasted image of the valve can be obtained when scanned in the open air. High-resolution MRI of the valve in liquid-bath has an anisotropic resolution of 80 µm (7 tesla field), enough to model the valve in 3D with accuracy [START_REF] Stephens | High resolution imaging of the mitral valve in the natural state with 7 Tesla MRI[END_REF]. These ex-vivo scanning techniques are also able to detect the subvalvular apparatus and the chordae tendineae insertion points. A realistic 3D fibre field can be deduced from the geometric definition of these scans [START_REF] Toma | High-Resolution Subject-Specific Mitral Valve Imaging and Modeling: Experimental and Computational Methods[END_REF] (figure 1.9,c). Ovine [START_REF] Drach | Population-Averaged Geometric Model of Mitral Valve From Patient-Specific Imaging Data[END_REF], [START_REF] Toma | High-Resolution Subject-Specific Mitral Valve Imaging and Modeling: Experimental and Computational Methods[END_REF] and porcine [START_REF] Villard | Fast image-based mitral valve simulation from individualized geometry[END_REF] mitral valves have been imaged using these ex-vivo methods. To be scanned, the excised valve must be fixed by the annulus and papillary muscle. A flow is used to maintain a the valve in a static position during the scan (open or close).

Geometry based mitral valve model

Another approach to recreating the mitral valve is to use a geometric shape and to adapt it using mitral valve measurements. The measurements can be obtained from in-vivo scans or anatomical measurements on cadavers. Such geometric models of the mitral valve were the first to be used for finite element purposes, and provide a good compromise between complexity and realism.

Anatomic measurements from cadavers

Kunzelmann and co-workers, in 2007, produced the first complete Fluid-Structure Interaction (FSI) study of the mitral valve based on anatomic measurements from cadavers [START_REF] Kunzelman | Fluid-structure interaction models of the mitral valve: function in normal and pathological states[END_REF]. A similar model was reused by the group of Burriesci to represent the influence of Alfieri's edge-to-edge repair technique on the mitral valve function [START_REF] Lau | Fluid-structure interaction study of the edge-to-edge repair technique on the mitral valve[END_REF]. Due to their appealing simplicity and the accessibility of cadaveric measurements, such models are still being developed and used [START_REF] Shen | The geometric model of the human mitral valve[END_REF].

3D-TEE landmarks

A parametric mitral valve model can be adapted using an in-vivo imaging technique. The 3D-TEE scan can measure the annulus shape, the position of the papillary tips, and the height of the leaflets. The data are used to modify and adapt a geometric mitral valve. Using nine landmarks on a 3D-TEE (left/right trigones, commissural tip locations, anterior/posterior leaflets, papillary muscle tips), the mitral valve movement can be tracked. The movement of the annulus and papillary muscles were implemented based on the in-vivo measurements as evolving boundary conditions [START_REF] Leonardo | Mitral Valve Finite-Element Modelling from Ultrasound Data: a Pilot Study for a New Approach to Understand Mitral Function and Clinical Scenarios[END_REF], [START_REF] Votta | Mitral valve patientspecific finite element modeling from 3-D real time echocardiography : a potential new tool for surgical planning[END_REF]. Using the landmark positions, the mitral valve was modelled in 3D (figure 1.10a) [START_REF] Mansi | An integrated framework for finite-element modeling of mitral valve biomechanics from medical images: Application to MitralClip intervention planning[END_REF], [START_REF] Grbic | Personalized mitral valve closure computation and uncertainty analysis from 3D echocardiography[END_REF]. Information on the chordae tendineae was not available. They were reconstructed from ex-vivo measurements and the literature.

(a) Mitral valve model reconstructed from landmarks set on 3D-TEE [START_REF] Grbic | Personalized mitral valve closure computation and uncertainty analysis from 3D echocardiography[END_REF] (b) Mitral valve from µ-CT scan and leaflet surface deduced from the scan. [START_REF] Khalighi | Multi-resolution geometric modeling of the mitral heart valve leaflets[END_REF] Figure 1.10: Geometric models of the mitral valve µ-CT scan simplification µ-CT scan can be simplified to represent the valve. The group of Sacks has developed mathematical models to extract relevant information from a µ-CT scan and applied it to a parametric model. Leaflet surfaces, describe by superquadric, fit closely the scanned leaflet shape. The technique was applied successfully to fit population-averaged mitral leaflets [START_REF] Drach | Population-Averaged Geometric Model of Mitral Valve From Patient-Specific Imaging Data[END_REF]. The irregularity of the surfaces is decomposed using frequency interpolation and can be reconstructed with several levels of accuracy [START_REF] Sacks | On the Simulation of Mitral Valve Function in Health, Disease, and Treatment[END_REF]. The complex subvalvular apparatus is decomposed using Reeb graphs, bifurcation and non-planarity angles distributions, cross-sectional area, and length distribution. The model also uses insertion sites on the mitral valve leaflet. Depending on the precision desired for the chordae tendineae, several equivalent models can be used [START_REF] Khalighi | Development of a Functionally Equivalent Model of the Mitral Valve Chordae Tendineae Through Topology Optimization[END_REF].

Parametric model of the mitral valve

Several methods are available to obtain the mitral valve shape, as depicted in the previous section. Our main objective is to develop a medical device, so the method chosen must be able to represent several valves with healthy and pathological variations, and be as realistic as possible.

No previous model of the mitral valve was available in our group, so a decision between geometricallybased model and scanned-based model had to be made. The most realistic scanning method is the µ-CT scan. To function properly this method needs ex-vivo mitral valves, as the contrast is higher between the valve component and air than with blood. This kind of model is not used in medical routine, so it is difficult to obtain several healthy and pathological valves. We were not successful in obtaining a coherent mitral valve using a standard multi-slice CT scan, so an automated method to obtain several valves seem out of reach with this method.

We choose to define the mitral valve using a parametric geometry, that can be adapted to in-vivo patient-specific data or anatomical measurement. It is the only current possibility to model the mitral valve using patient-specific data. The first section details a parametric model based on two separated codes, one written in Fortran and another written in Python, used as a FreeCAD macro. The second section explains the methods used in another parametric model, written in Python.

Parametric model using Fortran and Python

The definition of the parametric model is divided into several steps. The annulus is first described. The leaflets are then defined from the annulus, toward the papillary muscle, simplified to two anchorage points. The papillary muscles and the leaflet border are used to generate the chordae tendineae. The shell model of the mitral valve can be thickened using an additional thickness parameter.

In regards to the code used, the key landmarks of the annulus, leaflets and papillary muscles are created using homemade Fortran code. A homemade Fortran library is used to export the data in a Python formalism, that is later read into FreeCAD. FreeCAD geometric functions must be used to define the chordae tendineae, thus a Python module (a library) was written for the task. FreeCAD imports the Python module, interprets the Python code given by Fortran subroutines, and creates geometrical elements such as line, surface, or volumes. The numerical model of the valve is then exported as an IGES file, which can be used in finite element software.

In this study, the parameters control the annulus planar shape, and the length, thickness, and inclination of the valve leaflets. The number of principal and marginal chordae branches and the distances between insertion sites are also controlled by parameters.

The annulus definition

The annulus is assumed to be nearly flat during diastole. In this definition, the annulus is contained in a plane and has a geometrical D-shape. This shape is used in the Carpentier annulus ring. To define the annulus, two parameters are used, the inter-commissural and the anteroposterior lengths in addition to the geometrical relations of the D-shape (figure 1.11b and 1.11a). The annulus is discretized using points at specific locations, where corresponding leaflet heights are known. At these locations, the points are referred to the annulus plane using a radial angle. To define the radial angle, the known leaflet points are projected on the annulus plane, and the projection is decomposed using cylindrical coordinate. Once all the annulus points are defined by the Fortran code, they are written in a Python file. This file will be read by FreeCAD to create the annulus using B-splines (figure 1.11c). A B-spline is a curve defined from a set of ordered points. The same set of points used backwards thus produces a different curve. This is due to the B-spline construction process, in particular to the different boundary conditions for normal directions at starting and ending nodes. At least two B-spline must be used in the model. However, at the junction between the B-splines, a discontinuity in the slope of the edges appeared due to the B-spline construction. This discontinuity was non-physical and led to unrealistic stress concentration in the area.

To correct this default, additional points were added to the B-spline near the junction to reduce the discontinuity to a negligible level. Parabolic and spline interpolations were coded in Fortran to construct these additional points (figure 1.12). These interpolations can be implemented and solved explicitly as they only use 3 points. The increment ∆x between the points was set to 10 -3 cm after some trials. The new points are close enough to the original one so they do not modify where the B-spline passes but balance the slope at the point of interest to ensure tangent continuity between B-splines junctions. Steps of the method used to smooth the B-splines junction using parabolic interpolation.

The leaflet definition

The mitral leaflets are defined using a translation of the mitral annulus. Each annulus point is translated into a leaflet free-edge point using a specific vector. The leaflet free-edge points are then linked together into B-splines, and a surface is created from the annulus and the leaflet curves (figure 1.14). Annulus points are associated with a length, representing the leaflet height for this location on the annulus. For each point representing the annulus, a specific vector is used to translate the annulus point into a leaflet free-edge point. The vector length is set to the measured leaflet height, which can Figure 1.13: Example of leaflet height function of annular length. Data from [START_REF] Lau | Numerical simulation of mitral valve function[END_REF] be represented as a function along the annulus length (figure 1.13). The vector is mainly orthogonal to the annulus plane. It tilts towards the annulus centre using an inclination angle, defined for the anterior and posterior mitral leaflets (figure 1.14a). The angle of the leaflets towards the annulus plane help define the valve in a more realistic position during diastole. Once all the leaflet free-edge points are defined, they are linked together using B-splines. Additional points are used to ensure slope continuity. They are defined using the same method as for the annulus B-splines. Using the leaflets' free-edges and the annulus B-splines, a geometrical loft can be defined (figure 1.14c). A loft is a smooth surface with normal and tangential continuity (C 2 ) which includes the lines given as input. It depicts a spline surface and allows the definition of boat hulls or aircraft's wings with their ribs. In the case of the mitral leaflets, it ensures the leaflet curvature continuity. If desired, additional B-splines can be defined with other inclinations angles and leaflet height profiles. The new B-splines are then added to the loft structure as a new rib, which can increase valve realism.

If need be, the surface obtained from the loft can be uniformly thickened. This feature is useful to model the mitral valve with 3D finite elements instead of shells.

The chordae tendineae definition

The mitral chords are very different from a patient to another. They are complex structures, thus often simplified when modeled. The marginal chordae are responsible for the geometric closure of the valve and are nearly always represented. Basal chordae are mainly responsible for the stress field distribution and have a minor role in the leaflet coaptation.

In this parametric model, only the marginal chordae were model, as they have a primary role in valve closure. Strut chordae are neglected. The papillary muscles were simplified to two anchorage points, as usually done in the literature [START_REF] Lau | Numerical simulation of mitral valve function[END_REF], [START_REF] Kunzelman | Fluid-structure interaction models of the mitral valve: function in normal and pathological states[END_REF], [START_REF] Prot | On modelling and analysis of healthy and pathological human mitral valves: Two case studies[END_REF].

The chordae are modelled using a branching pattern that arises from the papillary muscles and insert into the edges of the leaflet. Principal chords arise from the papillary muscle, split into marginal chords, and merge into the leaflet edge. A regular fan dichotomous shape of the chordae is supposed (as depicted in figure 1.7). The cross-section is uniform on the chordae.

The parameters are the number of main chordae tendineae, the number of sub-branches for each main chord and the distance between the sub-branches. The main chords are uniformly distributed on the edges of the leaflet. The location of insertion points is dependent on both the leaflet geometry and the previous parameter. The mitral leaflet is defined using FreeCAD loft and B-spline functions, thus a specific Python code in FreeCAD had to be made to define insertion sites and the chords.

To create the insertion points, a function was created that splits leaflet free-edge, composed of B-splines. This function was applied twice, one to obtain the position for the main chords and the second time to obtain points for the subbranches. The set of subbranches was then sorted, where points too close or too far from the main chordae were deleted (figure 1.15). Additional "ghost points" were also kept and stored in another set. Once defined, all insertion points are sorted into subgroups. An iterative sorting process balances all the insertion points between each subgroup to minimize the distance between elements of the subgroup (figure 1.16, 1-2). Each subgroup is associated with its closest papillary muscle. Then the chord construction occurs, using the subgroups and the papillary muscles (figure 1. [START_REF] Lau | Numerical simulation of mitral valve function[END_REF][4][START_REF] Veronesi | Quantification of Mitral Apparatus Dynamics in Functional and Ischemic Mitral Regurgitation Using Real-time 3-Dimensional Echocardiography[END_REF][6]. In case the subgroup is relatively close (10%) to several papillary muscles, the subgroup is further subdivided and shared between the several papillary muscles. Additional "ghost points" defined with the insertion points are then added to the subgroups to enable the construction of the chord ramification.

The chord ramification process starts by linking the papillary muscle to the centre of the subgroups (figure 1.16, 4). A ramification area is defined with the maximal distance between elements of the subgroups. Sub-branches are made inside this area using the number of insertion points and the ramification pattern. Only one pattern is currently available in our code, the dichotomous shape pattern (figure 1.7, c). Some options are available to fill the ramification area with cone-shape Bsplines. This can be useful to study mesh variations.

The geometric mitral valve model developed in this study depends on several parameters. The size of the annulus, the height of the leaflets, the number of chordae, and their density can be used to describe several mitral valves in diastole position.

Parametric model based on FreeCAD only

The previous model had several flaws. First, it was limited to the description of the mitral valve in the open configuration (during diastole). Second, it was necessary to compile and execute a Fortran code before interpreting another Python script. Finally, it was complex to work on Python's script input via the Fortran code, as it requires many output calls. As a consequence, a second version of the parametric model was made. It is contained in one Python code called a FreeCAD module. It can be easily used and shared and uses FreeCAD features easily.

In this model, the method used to construct the mitral valve model is straightforward. First, an annulus of any type can be defined. Then the leaflets are created using the annulus and superquadric surfaces, similar to the model of the group of Sacks [START_REF] Sacks | On the Simulation of Mitral Valve Function in Health, Disease, and Treatment[END_REF]. The leaflets' free edges are then detected, and the chords are created.

Annulus definition

The annulus is defined as a list of consecutive points. The points must be described in a clockwise direction from the atrial view. The first point is associated with the middle of the posterior annulus. The parameters used to construct the annulus are the anteroposterior length, the inter-commissural length and the number of points used to discretize the annulus. Alternatively, the diameter and the perimeter can also be used if a shape is given (D-shape or ellipse). Planar D-shape and ellipse annulus are currently coded, but any shape can be easily implemented. The D-shape method was converted from the first parametric model in Fortran, detailed in section 1.3.1.1. Planar shapes can also be deformed along the z-axis using a parabolic function to form a saddle-shaped annulus (figure 1.17). An additional parameter is required to define the anterior and posterior height, in mm.

Leaflet definition

A key feature of this new model is the leaflet definition. The main idea is to use a superquadric surface to draw the leaflet onto, similarly to the model of Sacks [START_REF] Sacks | On the Simulation of Mitral Valve Function in Health, Disease, and Treatment[END_REF]. This surface can be used to represent with realism the mitral valve in several configurations (opened to nearly closed). The equation used to define the superquadric surface of normal z is:

1 = C 1/ 1 + A 1/ 2 + B 1/ 2 2 / 1 (1.1) with A = x-x 0 a 1 2 , B = y-y 0 a 2 2 , C = z-z 0 a 3 2
. Parameters a 1 , a 2 , a 3 are used to modify the shape size in each corresponding axis. The shape is centered in (x 0 , y 0 , z 0 ). 1 and 2 define the squareness factor. 1 influence the z axis, and 2 influence both x and y axis. The squareness factor is always positive. When set to 0, the shape is squared. If set to 1, the shape is round. If set to 2, the shape is a rhombus. If superior to 2, rhombus angles decrease and a star-like shape is observed.

A set of points represent the annulus. Each point of the annulus is associated with its angle in a cylindrical coordinate (figure 1.18a). Using this angle, a height and the superquadric equation, a superquadric translation vector is defined (figure 1.18b). Once all points of the annulus are translated for each desired rib, points are linked with splines. To get coherent splines, the interpolation method used in the previous model is once again used (figure 1.12). An intermediate loft surface is made from the ribs created (figure 1.18b, light grey). The number of ribs is related to the spatial discretization of the superquadric shape along the z-axis. 100 ribs between the annulus and the height a 3 /2 is a good compromise between efficiency and accuracy. This surface is the skeleton on which the mitral leaflets are drawn.

To draw the leaflet edges, new points are added on the surface and their distance to the annulus is compared to the leaflet height definition along the annulus. Once all points are at the correct height, a spline is defined (using the interpolation feature previously developed) to represent the leaflets height (figure 1.18c). This last spline is used as a cutting line to define the leaflet free edges. If desired, the mitral surface can be thicken using an additional parameter.

An additional parameter is the direction of the superquadric. Indeed, the mitral leaflet might not be oriented perpendicularly to the annulus plane. The vector z can be modified to represent other orientations for the mitral valve.

Chordae definition

In this second parametric model we use the same hypothesis as previously (only marginal chords are modelled, uniform cross-section, two papillary muscles are used, and chordae tendineae shape is assumed as a dichotomous shape (figure 1.7).

The chords are defined in consecutive steps. First, the free edge of the leaflet is discretized to obtain insertion sites. The set of points obtained is then separated into groups. Each group is associated with a main chord and a papillary muscle. A chordae is then defined for each group, using a pattern. In this model, the free edge of the leaflet is defined using only a spline. Using this spline directly would result in uniform distribution of the chordae based on the spline length, and would not be realistic. To determine the location of each chord, an angle theta can be given in the cylindrical coordinate of the annulus. The free edge's spline is divided according to the angles given, and each part of the spline is further discretized using the method of the previous model to obtain insertion points. These points are then separated into groups, and coupled to the closest papillary muscle, using again the method developed in the previous model. A chord is created using each pair of papillary muscle and insertions sites. The shape of the chord is a regular fan dichotomous, but other shapes could easily be implemented using the same technique.

The new parametric model is based only on Python, to compare to the use of Fortran and Python. It is easier to modify the valve or the method as only one file is responsible for the function generating the valve. This file is a FreeCAD Macro, and therefore can be easily installed in the software. The main file, which contains the function calls and the data/parameter, can be read directly in the software. This help sharing, using and upgrading this parametric model of the mitral valve. In addition, this new parametric model uses geometric objects and tools develop in FreeCAD. They help the construction of the valve and make it easier to understand or adapt the code. The shapes and structures obtained can be exported in any format, which is helpful for finite element modelling.

Model obtained

The methods defined in the section 1.3 were applied with the data obtained in the section 1.1.2. The rings obtained with both models are detailed with their parameters. Several leaflets are presented using different shapes obtained from the literature and both models. Chordae tendineae are defined using several types of insertion and the number of insertion points.

Annulus

Several annuli were described using the two methods. They both create splines and points in FreeCAD and use the same construction technique. Thus, they produce the same D-shape annulus, and there are no differences between the two models for this shape. The fully Python model also creates ellipses and morphs a flat annulus into a saddle shape. Data of inter-commissural and anteroposterior length from the section 1.1.2 were applied to both models (figure 1. [START_REF] Gunnal | Study of mitral valve in human cadaveric hearts[END_REF]), with several shapes. When not available, diameters and perimeters were converted into inter-commissural and anteroposterior length to construct the annulus (figure 1. [START_REF] Gunnal | Study of mitral valve in human cadaveric hearts[END_REF]). 19: Annuli defined using several techniques. Data from [START_REF] Lau | Numerical simulation of mitral valve function[END_REF] (blue), [START_REF] Khalighi | Development of a Functionally Equivalent Model of the Mitral Valve Chordae Tendineae Through Topology Optimization[END_REF] (green) and [START_REF] Oliveira | Geometric description for the anatomy of the mitral valve: A review[END_REF] (black). 100 points were used for each definition.

Leaflets

To assess only the effects of the leaflets' shapes, the same D-shape annulus was used in all cases, with inter-commissural and anteroposterior length set to 30 and 22.5 mm respectively (blue annulus, figure 1.19). It was morphed into a saddle shape of height 4 mm when used with the second method.

The geometric parameters used for the superquadric are in general unknown. They can be partially deducted from the annulus. a 1 and a 2 are directly set as half of the intercommissural and anteroposterior lengths. a 3 must be related to the maximal height of the leaflet. In this study, a 3 is set to 3 times the maximal height. The squareness factors 1 and 2 were fixed to 3.5 and 0.65 after trials. Ideally, all these parameters should be obtained from a minimization problem using anatomical data.

Leaflet profiles were used to define several mitral leaflets with the two methods (figure 1.20). Using the first model, the leaflet shapes obtain were straight and inclined toward the papillary muscles (transparent leaflets, figure 1.20). With the second model, the leaflets are also bent (opaque leaflet, figure 1.20). This additional curve reduces the total area of the leaflets, resulting in a more compact mitral valve with the same maximal height and input profile. As a result the valve can more easily close while having the same coaptation length. The leaflets obtained from the superquadric are more alike the leaflets from µ-CT scans than initial straight ones obtained.

Using the second model but leaflet profiles initially used for the first one results in non-closing valves. For example, the leaflet anterior shapes obtained in figures 1.20, a & b were too narrow to coapt with the posterior leaflet on a large area. As a result, the valves could only close partially. However, profiles obtain from µ-CT scan (figure 1.20, c) result in larger leaflets, which improve closure when used in finite element simulation. The annulus profile can sometimes be extracted with the leaflet profile [START_REF] Sacks | On the Simulation of Mitral Valve Function in Health, Disease, and Treatment[END_REF]. Using the same method as for the leaflet free edge, the annulus can be drawn on superquadric surfaces and be used to define the second border of the mitral valve (figure 1. 20,c).

Chordae tendineae

To represent the implementation of the chordae tendineae, the same annulus and leaflet are used. The annulus is D-shaped, with a saddle shape of height 4 mm. Inter-commissural and anteroposterior lengths were set to 30 and 22.5 mm respectively. The leaflet and annulus profiles are obtained from µ-CT scan [START_REF] Sacks | On the Simulation of Mitral Valve Function in Health, Disease, and Treatment[END_REF] (figure 1.20, c). The same papillary muscles are used, defined as points, and located using the literature ( [START_REF] Sakai | Distance between mitral anulus and papillary muscles: Anatomic study in normal human hearts[END_REF], see figure 1.5).

Several insertion sites are depicted in figure 1.21. Subfigure a shows uniform distribution for the main chords, while specific distribution is imposed in subfigure b. The distribution favours the free edge border bumps and depression, which is a more common pattern for chord implementation. The number of main chords is used to define the distance between two insertion sites. The number of sub-branches and the distance between sub-branches can define the area of insertion. Figure1.21,a depicts large and dense insertion sites while figure 1.21,b has a smaller inserting area.

Two methods are available to join the papillary muscles to the leaflet's free edge. The number of ramifications is obtained from the elements in the insertion site and the shape chosen for the chord. Setting the shape of the chord to dichotomous produces a ramification pattern using all the insertion points (figure 1.22, a & c). At locations where the free edge is close to both papillary muscles, the ramification is modified so that the insertion points are shared between the two muscles (figure 1.22, c). Using a spline method for the chord shape imposes normal continuity between the sub-branches and the leaflet border (figure 1.22, b). The insertion points previously defined are used to ensure border continuity. No ramification is created, but the insertion area between the spline and the leaflet edge can be filled. This is useful to test insertion area independently of mesh size when running finite element simulations.

Valve created using the first parametric model are shown in figure 1 1.23, c is between systolic and diastolic shapes. It is easier to use this last model for finite element simulation as the referent configuration is likely to be similar to the valve c. In addition, the total leaflet surface of this valve is smaller than the surfaces of valves a and b. This improves the behaviour of the mitral valve during closure and helps numerical convergence, as no leaflet bulges will be produced. It is clear that the valve in figure 1.23, a, will result in unrealistic mitral valve finite element simulation. The implementation of the chords will favor high-stress locations and lead to convergence difficulties. Larger insertion sites, like figure 1.23, b & c, will distribute the force applied by the chord on the leaflet during closure. The stress will be more evenly distributed, and the simulation will be more realistic. Two parametric models were developed in this study. The first one used Fortran and Python scripts. The second one was implemented in a FreeCAD macro, using Python. Both models could model a large variety of mitral valves. The second model has additional features that allow modelling valves that are more realistic, and more fitted for finite element simulation. This FreeCAD only model is also more easily used, implemented, and shared. It is defined per step, so it is possible to work only on a single item (annulus, chord shape, chord location...) and still use other functions of the macro to produce a valve with the new features.

The parametric models developed are dependent on data used to model the valves. The parameters are linked to the annulus shape and size, to the leaflet evolution function of the annulus length, to the number of main marginal chords, and the density of chords at the marginal insertion sites. Data can be obtained from in-vivo imaging techniques or cadaveric measurements. The orientation of the valve and the location of the papillary muscles can be used to shape the valve in a specific position to fit in a left ventricle already modelled in 3D, for example. Pathological variations can be modeled from the modification of a healthy valve using the FreeCAD tools or from specific parameters that lead to a pathological valve.

These codes are efficient tools to create several mitral valve models. A large panel of mitral valves can be obtained using the data from the literature. Such a panel can be quickly generated and implemented in a finite element model. Doing so, medical devices could be tested on pathological valves and optimized to ensure that the same device is fitted to improve several valve functions.

Improvements and discussion

The parametric models developed describe with realism healthy and pathological valves, along with their functions. The second model results in more detailed valves and should be preferred. Possible improvements to the model are discussed in this section. The annulus and papillary muscles are mobile during a cardiac cycle but are considered fix in our model. The chordae tendineae presented could be more realistically implemented on the valve leaflets, and/or have other shapes. Their length, not clearly defined in the models, must be adapted to ensure the coaptation of the leaflets. Finally, the mitral leaflets and chords are composed of an anisotropic material. The geometric definition of the valve should give a map of the fibre preferred direction in the referent configuration.

Annulus and papillary muscles movement

The second parametric model was defined from a static position when the valve is half-opened. As a result, the annulus has a static saddle shape. A finite element simulation of the valve should take into account the dynamic deformation of the annulus, for example by imposing it as a boundary condition. The movement of the papillary muscles, although smaller in magnitude when compared to the annulus, could also be taken into account [START_REF] Leonardo | Mitral Valve Finite-Element Modelling from Ultrasound Data: a Pilot Study for a New Approach to Understand Mitral Function and Clinical Scenarios[END_REF]. The dynamic deformation of the annulus helps to close the valve during systole, according to the same article.

Chords

The parametric model creates only marginal chords, which are a part of the whole subvalvular apparatus. Marginal chords are responsible for the valve closure during systole, but the basal chords help transmit the ventricular load to the leaflet. Modeling only the marginal chords thus results in higher stress near the insertion area, which can be lower to a limit by a larger density of chord on the leaflet border. A major improvement for the model would be to implement the basal chords. This would also alter the valve shapes, reducing the actual curvature of the leaflets during systole. The strut chords should be model for the same reasons.

Structural finite element simulation of the mitral valve tends to focus on the chord's function rather than on chord realism. The chord's function is to maintain the mitral valve closed during systole. If the interaction between the chord and the fluid is neglected, the valve function can be studied accurately while replacing the chord with cables of no apparent thickness. These cables are thus not detected during fluid-structure interaction, and only serve to maintain the valve closed. By doing so, a more generic approach can be applied, like multiple cable implementation with no ramifications. In this case, the cable density on the whole leaflet can be monitored to provide equivalent closure as a full chordae tendineae model [START_REF] Khalighi | Development of a Functionally Equivalent Model of the Mitral Valve Chordae Tendineae Through Topology Optimization[END_REF]. A major drawback of this model is the lack of possible interaction between the chords and the fluid, which is not negligible [START_REF] Toma | Fluid-Structure Interaction and Structural Analyses Using a Comprehensive Mitral Valve Model with 3D Chordal Structure[END_REF].

The current approach aims to represent realistic chords and their branching pattern. While not complete, the model can easily switch between chords of no thickness, invisible to the fluid and thickened ones, that can interact with blood. The methodology developed can be improved by adding a new branching pattern, and new chords implementations. An interesting approach would be to reuse the classification of chords made by the groups of Sacks [START_REF] Khalighi | Mitral Valve Chordae Tendineae: Topological and Geometrical Characterization[END_REF]. This data, coupled with the branching pattern described by Gunnal and co-workers, could improve the realism of the chord generated in our model [START_REF] Gunnal | Morphological Study of Chordae Tendinae in Human Cadaveric Hearts[END_REF]. The ramification of the chord is very dependent on the leaflet geometry, and it is possible to optimize the branching pattern to support the load during systole. As an example, a tree-like generation used in other modelling is biased using a force vector and could be used to represent the load on the mitral valve [START_REF] Runions | Modeling Trees with a Space Colonization Algorithm[END_REF]. This generation has the advantage of using only the geometry of the valve to create a branching pattern. Other optimization shapes, like the one provided in FreeCAD, are interesting candidates for such geometric optimization but are based on finite element resolution (figure 1.24).

A final aspect of the chordae tendineae that was not discussed in our model is the mechanical rest configuration. Independently of the chords modelling, the length at rest is a difficult subject to tackle because it modifies the behavior of the valve, its shape, and stress distribution during closure. In most studies, this length is set to optimize the closure using either trials and errors or some retroactive loop based on simulation results. The length at rest is available in very few studies. This lack of information is mainly linked to the non-detection of the subvalvular apparatus during in-vivo scanning. The studies with ex-vivo µ-CT scans are the only ones ignoring this difficulty as they scan and model the subvalvular apparatus entirely.

Fibre mapping definition

The mitral valve is composed of anisotropic material, but the fibre orientation was not discussed here. They are key components of mitral valve high resistance under blood pressure. Studies of mitral fibre orientation proved that chordae branches locally modify the fibre orientation. 2D fibre mapping are available using SALS technique and are often used as references in the literature [START_REF] Cochran | Nondestructive analysis of mitral valve collagen fiber orientation[END_REF], [START_REF] Sacks | Incorporation of experimentally-derived fiber orientation into a structural constitutive model for planar collagenous tissues[END_REF]. The definition of the fibre map should take into account chords' location on the leaflets and the leaflet geometry. Ideally, once the geometry of the valve and its subvalvular apparatus are defined, an optimization problem should be performed to orientate the fibres. The resulting laplacian problem to solve has boundary conditions defined from the observation of the SALS maps. The mapping obtained is most likely to be realistic, but can only be deduced from the geometry and not directly observed [START_REF] Toma | High-Resolution Subject-Specific Mitral Valve Imaging and Modeling: Experimental and Computational Methods[END_REF].

Chapter 2

Mechanical behaviour of the mitral valve

The mitral valve is a biological material. Most of the biological materials are complex compounds and are difficult to use in finite element simulation. Simplification of the material law is a common process to ease the structural simulation. However, a simplified law can fail to represent the physical problem accurately. Furthermore, the high mobility of the mitral valve requires a complex formulation to be modelled, which is an additional difficulty.

This chapter presents what the mitral valve components are and how they can be modelled effectively. Isotropic material models are then discussed, and their defaults are presented. Finally, a more generic approach to anisotropic materials and their implementation in ADINA is aborded. The mitral valve is composed of 4 main elements (see figure 2.1). The papillary muscles are considered as anchoring points for the chordae tendineae and so are modelled as points. Similarly, the annulus of the mitral valve is incorporated in the fibrous trigone, a part of the fibrous skeleton of the heart. When not considered fixed, the mitral annulus and the papillary muscles are given a cyclic deformation to help the valve function. In general, this movement is not computed from the ventricle contraction but is directly implemented as a boundary condition. Our goal is to simulate the valve closure and not the exact behaviour of the heart. Consequently, only the leaflets and the chordae tendineae of the mitral valve must be modelled to simulate the response of the valve to cyclic pressure variations. Both materials are mainly composed of elastin and collagen and constitute anisotropic material. The chords and the leaflets have a continuous composition of collagen fibres and elastin, but the organization of these two components varies with the valve's geometry. In principle, they could use the same material law. However, in practice, the geometry of the chords is so specific that a simplified material law is used with very little accuracy loss.

Description of the mitral valve components

Chordae tendineae

The chordae tendineae are a compound of biological elements that is difficult to model. They are composed of elastin and collagen fibres that are aligned with the axis of the chords. Their function is to retain the leaflets in the left ventricle while the valve is closed. Because of their cable-like shape and function, they can be model as cables. Doing so allows to simplify the anisotropic material and use an isotropic one instead. This approach is cost-effective and accurate as long as the chord work in tension or compression along their axis. This section describes the chords from a biological point of view, then discusses the isotropic model used in the literature.

Physiological composition

The structural composition of human chords was observed in detail with light and electron microscopy by Millington and Sanders [START_REF] Millington-Sanders | Structure of chordae tendineae in the left ventricle of the human heart[END_REF]. They found that chords have a core made of collagen fibres encircled by the first layer of straight elastic fibres aligned with the chord direction. The second layer of elastin fibres surrounds the first one, with an inclination angle. An endothelium surrounds the whole structure. The collagen fibrils that compose the core of the chord are waving at about 10 µm (figure 2.2). The wavelength of the collagen fibrils differs between chords. The wavelength of the collagen fibrils is uniform through a chord but differ between chords. Their value seems to depend on the chords thickness [START_REF] Liao | A structural basis for the size-related mechanical properties of mitral valve chordae tendineae[END_REF], [START_REF] Liao | Ultrastructure of porcine mitral valve chordae tendineae[END_REF].

Mechanical behaviour

The collagen's fibres and the elastin are together responsible for the mechanical response of the chord. Studies of the mechanical response of the chords show a hyperelastic response to loading with a small viscoelastic effect.

The wavy organization of the collagen's fibre provide a smooth transfer of forces during the chord tension. Once the tension is released, the elastic fibres help the collagen return to its wavy form, helping the smooth tension variation. As a result, the load transmit at low strain is mainly due to the elastin response. When the strain increases, the collagen's fibre gradually tightens and aligns with the deformation, transmitting additional load. This phenomenon is called fibre's recruitment and explains the hyperelastic response of most biological material. The elastin stiffness is about one-tenth of the collagen stiffness; thus, the mechanical response is highly dependent on the material strain. Uniaxial testing of the mitral valve chords characterizes the smooth transition between elastin and collagen fibres (figure 2.3) [START_REF] Kunzelman | Mechanical properties of basal and marginal mitral valve chordae tendineae[END_REF]. Once all the fibres are recruited, the stiffness is the same for both marginal and basal chords, around 3 kg.mm -2 , i.e. 30 MPa. Under physiological range, the non-linear response of the chord is best fitted with an exponential law for the strain energy density function. However, a linear approach can be sufficient and is easier to implement [START_REF] Cai | Some Effects of Different Constitutive Laws on FSI Simulation for the Mitral Valve[END_REF]. The viscoelastic response of the chord was analysed using Dynamic Mechanical Analysis (DMA). Results show that the mitral chords are mostly elastic and slightly viscoelastic in the physiological range [START_REF] Lim | Dynamic elasticity of human mitral valve chordae tendineae[END_REF]. The phase lag induced by the viscoelasticity decreases with higher frequency. A small amount of viscoelasticity is present in the mechanical response and is not equal for each chord type. Thinner ones tend to be stiffer and more viscous than thicker ones. An interesting information is that the storage stiffness is frequently dependent but not significantly different between chords. Even if the mitral chords vary in thickness and stiffness, the actual force transmitted per elongation is nearly the same for all chords. Using uniaxial material testing, Millard and coworkers have proven the same independence of structural stiffness, while the young's modulus ranges from 20 to 80 MPa [START_REF] Millard | Mechanical properties of chordae tendineae of the mitral heart valve: Youngs modulus, structural stiffness, and effects of aging[END_REF]. Results are summarized in table 2.1. No significant mechanical differences exist between anterior and posterior chords [START_REF] Wilcox | Frequency and diameter dependent viscoelastic properties of mitral valve chordae tendineae[END_REF].

The maximum force a chord can transmit before rupture is around 5 N [START_REF] Sedransk | Failure mechanics of mitral valve Chordae tendineae[END_REF], [START_REF] Lodder | Morphological and Mechanical Properties of the Posterior Leaflet Chordae Tendineae in the Mitral Valve[END_REF]. According to an exhaustive literature review, the forces actually transmitted in-vitro and in-vivo range from 0.23 to 1.00 N per chord [START_REF] Rabbah | Mechanics of healthy and functionally diseased mitral valves: A critical review[END_REF]. Marginal chords transmit the lowest loads around 0.2 N [START_REF] Jimenez | Mitral valve function and chordal force distribution using a flexible annulus model: An in vitro study[END_REF], [START_REF] Lomholt | Differential tension between secondary and primary mitral chordae in an acute in-vivo porcine model[END_REF].

Experimental studies show that it is much easier to represent the chord using stiffness rather than the modulus. One model can be used for all the chords independently of their thickness. Furthermore, the viscoelastic response of the chords seems negligible.

young (n=7) old (n=5) Structural stiffness ± SD (kN/m) 2.0 ± 0.4 2.2 ± 0.2

Material law

Material laws used in the literature vary from linear model to specific hyperelastic formulation of the strain energy density function (the use of specific strain energy density function is explained in annexe). Linear models are easy to implement. Complex models are more realistic but also more difficult to solve numerically. Experimental studies show that even if they differ in size, chords have the same stiffness under elongation. The same linear or hyperelastic model can then be used for all the chords. Other approaches use specific material laws and diameters for each type of chord. The specific material law can be implemented using stress-strain relation or strain energy density function. This last approach is closer to the physic of the problem studied but is more difficult to implement.

The uniform simplification was used several times. For example, Watton and co-workers used 300 MPa chords stiff with 0.4 mm 2 cross-section [START_REF] Watton | Dynamic modelling of prosthetic chorded mitral valves using the immersed boundary method[END_REF], while Burriesci's group used 3 MPa chords with the same cross-section.

Non-linear stress-strain law figure 2.3 was used with varying cross area for the chords, from 0.45 mm 2 to 0.79 mm 2 for marginal and basal chords, respectively [START_REF] Kunzelman | Fluid-structure interaction models of the mitral valve: function in normal and pathological states[END_REF]. Burriesci and co-workers reused the same configuration [START_REF] Lau | Numerical simulation of mitral valve function[END_REF].

Several exponential strain energy density function laws were developed to model the chords mechanical response [START_REF] Prot | On modelling and analysis of healthy and pathological human mitral valves: Two case studies[END_REF], [START_REF] Prot | Transversely isotropic membrane shells with application to mitral valve mechanics. Constitutive modelling and finite element implementation[END_REF], [START_REF] Prot | Finite element analysis of the mitral apparatus: Annulus shape effect and chordal force distribution[END_REF]. They contain a linear term and an exponential one to fit both the elastin and the collagen's fibre recruitment. The law used by Prot and co-workers is:

W (I 1 ) = a 1 (I 1 -3) + a 2 [exp(a 3 (I 1 -3) 4 ) -1] + p(det(F ) -1)
with W the strain energy density function, I i the i-th invariant of deformation, a j , j = 1..3 the coefficients obtained from a fit to experimental data (table 2.2), p the quasi-incompressibility pressure enforcement and F the deformation gradient. The cross-sections were set to 0.29 mm 2 , 0.27 mm 2 and 0.61 mm 2 for the anterior marginal, posterior marginal and strut chordae, respectively. [START_REF] Prot | On modelling and analysis of healthy and pathological human mitral valves: Two case studies[END_REF] The transition between chords and leaflets is continuous. Using a realistic and detailed fibre mapping of the mitral valve, it is possible to use the same anisotropic strain energy density function for the leaflets and the chords [START_REF] Toma | Fluid-Structure Interaction and Structural Analyses Using a Comprehensive Mitral Valve Model with 3D Chordal Structure[END_REF], [START_REF] Toma | Fluid-Structure Interaction Analysis of Papillary Muscle Forces Using a Comprehensive Mitral Valve Model with 3D Chordal Structure[END_REF]. This approach requires a µ-CT scan but simplifies the definition of the chord mechanical response (see section 2.1.2.3). The group of Sacks developed another alternative.

a 1 (MPa) a 2 (
After defining the insertion area for each chord, the forces transmitted by the chords were directly adapted to obtain the best structural shape of the leaflet [START_REF] Sacks | On the Simulation of Mitral Valve Function in Health, Disease, and Treatment[END_REF].

The elastic modulus and the thickness varies among the chords, but the stiffness does not significantly differ. Thus a uniform approach can be used to model them. Two methods are mainly used, a linear elastic model or an exponential hyperelastic law. The methods can be implemented using a strain energy density function or a direct stress-strain approach.

Valve leaflets

The leaflet materials are key components of good mitral valve closure. Anterior and posterior leaflets must connect quickly to prevent leakage at early systole, so the material must be easily deformed at low stress. The leaflets also have to maintain a closed seal during the systolic peak. To do so, the material must remain still under high stress. The mitral valve leaflets are thin structures about 1 mm. Yet, they endure high stress 0.5 MPa with low strain, lower than 8% in the radial or circumferential direction according to in-vivo measurements [START_REF] Sacks | In-Vivo Dynamic Deformation of the Mitral Valve Anterior Leaflet[END_REF]. An equivalent elastic isotropic material would have Young's modulus of 3 MPa. Such a valve would be too stiff and wouldn't close at early systole. The composition of the mitral valve is key to understand its mechanical response. Mechanical testing of the leaflets will help characterize their response and eventually simplify their model.

Physiological composition

Histology of the mitral leaflets

The mitral valve, by dry weight, is approximately 70% collagen, 10% elastin and 20% glycosaminoglycans (GAGs), which is similar to the chords. Histological view of the cut valve display several layers (figure 2.4). From the atrium to the ventricle, we found the atrialis, the spongiosa, the fibrosa and the ventricularis, both sides covered by endothelial cells [START_REF] Rabbah | Mechanics of healthy and functionally diseased mitral valves: A critical review[END_REF], [START_REF] Fenoglio | Canine Mitral Complex[END_REF]. The four layers have thicknesses that vary on the leaflets (figure 2.4a). The main portion of the valve is composed of the fibrosa and spongiosa [START_REF] Fenoglio | Canine Mitral Complex[END_REF]. The fibrosa is the thickest layer comprised of dense organized wavy collagen fibres, providing the most load-bearing support to the leaflet. The atrialis and the ventricularis are mostly composed of elastin fibres (figure 2.4b). The elastin fibres facilitate the restoration of the leaflet in the undeformed configuration [START_REF] Kramer | An investigation of layer-specific tissue biomechanics of porcine atrioventricular valve anterior leaflets[END_REF]. This behaviour is similar to the role of the elastin fibres inside the chords. The ventricularis layer also serves as an insertion point for the chords and presents collagen continuity with them [6]. The atrialis contains a layer of loose connective tissue of collagen and elastin fibres supported by a proteoglycan and glycosaminoglycan gel in the spongiosa. The separation between the atrialis and spongiosa is sometimes hard to define, as both contain loosely organized collagen and elastin fibres (figure 2.4b). The hydrated spongiosa layers could help lubricate the atrialis and fibrosa, ensuring better stress distribution between the layer of different stiffness. This lubrication function has yet to be fully understood and detailed [START_REF] Kramer | An investigation of layer-specific tissue biomechanics of porcine atrioventricular valve anterior leaflets[END_REF], [START_REF] Sacks | Heart Valve Function: a Biomechanical Perspective[END_REF].

The crimped fibres unfold, straighten, and elongate as the leaflet undergoes loading during cardiac function, providing the tissue fibre-reinforced stiffening. At low stress, the waviness of the crimped collagen allows the leaflet to undergo substantial deformation even at a low level of stress. This behaviour is supposed to help coaptation and help close the valve.

Collagen fibre orientation

Mitral valve leaflets are mainly a compound of orientated collagen, which orientation seems responsible for the anisotropic response of the valve [START_REF] Grande-Allen | The heterogeneous biomechanics and mechanobiology of the mitral valve: Implications for tissue engineering[END_REF]. Their waviness can be observed using a polarized light microscope [START_REF] Rabbah | Mechanics of healthy and functionally diseased mitral valves: A critical review[END_REF], [START_REF] Grande-Allen | The heterogeneous biomechanics and mechanobiology of the mitral valve: Implications for tissue engineering[END_REF]. The degree of orientation (1 being perfectly oriented, 0 being totally random) ranged from 58% to 65% using small-angle light scattering (SALS) on porcine mitral valves [START_REF] Cochran | Nondestructive analysis of mitral valve collagen fiber orientation[END_REF], [START_REF] Billiar | A method to quantify the fiber kinematics of planar tissues under biaxial stretch[END_REF]. The SALS technique was compared to polarize light microscopy to ensure that the angle of rotation measured by the SALS corresponded to the collagen fibre orientation. The mapping obtained from the SALS technique is given in figure 2.5. The results suggest the fibres be mainly oriented parallel to the annulus in the central portion, then turn orthogonal to the annulus in the commissural region. The crimped collagen fibres are also continuous with the chords fibres [START_REF] Fenoglio | Canine Mitral Complex[END_REF]. Such highly oriented fibrous networks are believed to transmit loads most efficiently along fibre direction and behave as an anisotropic material [START_REF] Cochran | Nondestructive analysis of mitral valve collagen fiber orientation[END_REF].

visco-elastic behaviour

Viscoelasticity is common among biological tissue. It can be characterized by stress relaxation, hysteresis, creep and strain-rate dependence tests. Using uniaxial and biaxial mechanical testing, Grashow, Sacks and co-workers have proven that the mitral valve anterior leaflet have stress-strain responses remarkably independent of strain rates at and below physiological range. In addition, they proved that hysteresis was negligible at physiological rates, thought detectable (around 12%) [START_REF] Grande-Allen | The heterogeneous biomechanics and mechanobiology of the mitral valve: Implications for tissue engineering[END_REF], [START_REF] Grashow | Biaixal stress-stretch behavior of the mitral valve anterior leaflet at physiologic strain rates[END_REF]. This low hysteresis level would dissipate the energy during the cyclic opening and closing of the mitral valve. Another study proved that mitral valve anterior leaflet does not creep, although a significant stress relaxation occurs. The stress relaxation was about 25% in 3h and 6% in 1s [START_REF] Grashow | Planar biaxial creep and stress relaxation of the mitral valve anterior leaflet[END_REF]. In most biological materials, stress relaxation and creep are related, but this is not the case with the mitral valve anterior leaflet. Liao et co-workers have studied the microstructural processes responsible for this behaviour with mechanical testing and small-angle X-ray scattering [START_REF] Liao | The relation between collagen fibril kinematics and mechanical properties in the mitral valve anterior leaflet[END_REF]. The collagen fibres, under stress, reorganize and straighten out toward the principal direction, and once loaded, stay "locked" in this position. However, under imposed strain, the collagen fibre length reduces with decreasing stress.

These studies show that the mitral valve anterior leaflet is not a classic viscoelastic material and can be considered a quasi-elastic anisotropic material [START_REF] Grande-Allen | The heterogeneous biomechanics and mechanobiology of the mitral valve: Implications for tissue engineering[END_REF]. The main viscous characteristic of the mitral valve is stress relaxation, which is about 6% for 1 s.

Biaxial and uniaxial testing of the mitral leaflets

The stress-strain relationships of the leaflets can be separated into 3 zones. Initially, the leaflet undergoes a large strain for small stress, suggesting that the collagen fibres are crimped and not organized. The stress is mostly a result of the elastin fibres and the glycoprotein matrix. With increasing stress, the collagen fibre flattens and orients toward the stress direction, which produces a transitional stress-strain zone. Finally, the fibres are aligned and taut at higher stress, transferring a large amount of stress for a small strain. Results of biaxial mechanical testing revealed a strongly anisotropic and hyperelastic response of the leaflets (figure 2.6). Using equibiaxial and strip biaxial testing, a different behaviour was detected between anterior and posterior leaflets in radial and circumferential directions. The thickness difference between the two leaflets (0.79 and 0.7 mm for anterior and posterior leaflets, respectively) and the presence of collagen fibres in different proportions could explain the differences [START_REF] Kramer | An investigation of layer-specific tissue biomechanics of porcine atrioventricular valve anterior leaflets[END_REF]. The circumferential stress was 5 to 6 times the stress in the radial direction for the same strain during the equibiaxial test ( 2.3: Pre and Post-transitional elastic moduli from three-dimensional analysis for anterior and posterior leaflet during equibiaxial testing [START_REF] May-Newman | Biaxial mechanical behavior of excised porcine mitral valve leaflets[END_REF]. The R ratio represents the ratio of stress between circumferential and radial direction for a given strain.

Kunzelman and Cochran ran similar tests and obtained similar results [START_REF] Kunzelman | Stress/strain characteristics of porcine mitral valve tissue: Parallel versus perpendicular collagen orientation[END_REF]. They proved that the mechanical response is directly modified by the collagen fibre's density and orientation. In addition, biaxial testing of the valve shows that radial and circumferential directions were coupled, suggesting some kind of structural link with the fibres (figure 2.6).

Inverse finite element analysis was used to describe the leaflet property in-vivo under the physiological range [START_REF] Krishnamurthy | Material properties of the ovine mitral valve anterior leaflet in vivo from inverse finite element analysis[END_REF]. Briefly, the stress was obtained from the atrium and ventricle pressures of the sheep model. The strain was measured using radiopaque markers placed on the annulus, the anterior leaflet and the papillary muscles. An equivalent finite element model is used to reverse the strain and stress data to compute information about the material. The study is limited to isovolumic relaxation because it is the only cardiac cycle duration where stress can be computed accurately from pressure. During this period, the leaflets behave linearly with an approximate stiffness of 25 MPa (10 MPa) for the circumferential (radial) direction. These values are coherent with the values obtained ex-vivo for posttransitional stiffness (table 2.3). It proves that the leaflets, during systole, reach the post-transitional state. Therefore to accurately model the leaflet behaviour, a stiffness transition must occur between pre and post-transitional states.

A study of micro-dissected mitral valves in 3 layers (atrialis/spongiosa, fibrosa and ventricular) shows different mechanical responses between each layer. Atrialis/spongiosa layer behave anisotropically, but the fibrosa and ventricularis layers appear to behave isotropically [START_REF] Kramer | An investigation of layer-specific tissue biomechanics of porcine atrioventricular valve anterior leaflets[END_REF]. This new approach offers a better understanding of the valve behaviour and is important for damaging studies of the mitral leaflet.

The leaflets are composed of mainly aligned collagen fibres and elastic fibres. Their mechanical response is hyperelastic, with an order of magnitude of 30 between pre and post-transitional stiffness. Leaflets are anisotropic material, collagen fibres being mostly aligned parallel to the annulus. The leaflets are quasielastic but possess a small stress relaxation.

Numerical modeling

Mitral valve leaflets are anisotropic hyperelastic material and are difficult to model. The complexity of the material increases using hyperelasticity instead of elasticity and anisotropy instead of isotropy. Large deformation and large displacement of the leaflet is supposed. In this context, hyperelastic behaviour is best defined using the strain energy density function W , which relates to the second Piola Kirchhoff stress tensor S and Cauchy Green strain by S = ∂W ∂ . Details on the hyperelasticity formulation are given in the annexe.

Elastic material laws

The simplest material model is isotropic elasticity. Several groups have applied isotropic elastic material to the mitral valves. The standard seems to be 0.75 mm thick leaflets with a Young modulus of 3 MPa [START_REF] Lau | Numerical simulation of mitral valve function[END_REF], [START_REF] Lau | Mitral Valve Dynamics in Structural and Fluid-Structure Interaction Models[END_REF]. The values are coherent with post-transitional values obtained in the literature (table 2.3). A Fluid-Structure Interaction simulation with patient-specific data used the same parameters and model [START_REF] Biffi | A workflow for patient-specific fluid-structure interaction analysis of the mitral valve: A proof of concept on a mitral regurgitation case[END_REF]. A major leak remains due to partial closure.

Watton and co-workers defined mitral (and aortic) valves using a network of elastic springs, with respective modulus 5.4 MPa (and 15 MPa) and thickness 0.1 mm (and 0.125 mm) [START_REF] Watton | Dynamic modelling of prosthetic chorded mitral valves using the immersed boundary method[END_REF], [START_REF] Watton | Effect of ventricle motion on the dynamic behaviour of chorded mitral valves[END_REF]. The mitral valve leaflets are incorporated within a chord structure of stiffness 300 MPa, a disguised method to include anisotropy through fibres.

Ma and co-workers defined a less stiff mitral valve of 800 KPa reconstructed from MRI data [START_REF] Ma | Image-based fluid-structure interaction model of the human mitral valve[END_REF]. Villard and co-workers developed a real-time simulation (15s CPU) of a 400 KPa stiff mitral valve. This low stiffness values increase the coaptation of the valve leaflets but lack realism when the stress increases [START_REF] Villard | Fast image-based mitral valve simulation from individualized geometry[END_REF].

An elastic orthotropic model for shells was developed to capture the complex anisotropic response of the valve more accurately. This model was used in several studies, in particular, to test medical procedures [START_REF] Votta | The Geoform Disease-Specific Annuloplasty System: A Finite Element Study[END_REF], [START_REF] Mansi | An integrated framework for finite-element modeling of mitral valve biomechanics from medical images: Application to MitralClip intervention planning[END_REF], [START_REF] Schievano | Percutaneous Mitral Valve Dilatation: Single Balloon Versus Double Balloon. A Finite Element Study[END_REF], [START_REF] Xu | Augmented mitral valve leaflet area decreases leaflet stress: A finite element simulation[END_REF]. The model's stiffness ranged between 6.23 and 1.89 MPa, coherent with the post-transitional stiffness of the leaflet. The model used a fibre orientation obtained from the SALS map obtained by the group of Kunzelman (figure 2.5).

The major drawback of the linear models is the lack of coaptation obtained. Leaflets tend to fold and resist the coaptation [START_REF] Lau | Numerical simulation of mitral valve function[END_REF], [START_REF] Votta | The Geoform Disease-Specific Annuloplasty System: A Finite Element Study[END_REF], [START_REF] Lau | Mitral Valve Dynamics in Structural and Fluid-Structure Interaction Models[END_REF]. Lowering the leaflet stiffness increases closure to a certain extent but results in unrealistic deformations during systole. Mitral valve reconstructed from imaging technique reach better closure in average but does not present a coaptation area similar to biological response [START_REF] Mansi | An integrated framework for finite-element modeling of mitral valve biomechanics from medical images: Application to MitralClip intervention planning[END_REF], [START_REF] Xu | Augmented mitral valve leaflet area decreases leaflet stress: A finite element simulation[END_REF].

Hyperelastic material laws

Hyperelastic material laws can describe the fibre recruitment of most biological tissues. They help in obtaining coherent leaflet behaviour and the closure of the mitral valve. Adding an anisotropic effect also increases the realism of the leaflets. The best hyperelastic anisotropic material law also takes into account the fibre dispersion around a preferred direction.

Few isotropic hyperelastic materials were used to model the mechanical response of the leaflet. Mooney-Rivlin, 2nd order Ogden and 5th order polynomial hyperelastic models were calibrated from experimental data [START_REF] May-Newman | Biaxial mechanical behavior of excised porcine mitral valve leaflets[END_REF]. The expected increase of the closure area was not observed with Ogden and Mooney-Rivlin model [START_REF] Sharifikia | Structural simulation of human mitral valve behaviour cosidering effects of material nonlinearities[END_REF], [START_REF] Abedini | Numerical simulation of the impact of chorda tendinea rapture on the mitral valve leaflets stress by finite element method[END_REF]. It is not clear if the non-coaptation of the leafets is linked to the simulation configuration or the material law. High deformation of the leaflet and closure of the valve was obtained with the 5th order polynomial strain energy density function [START_REF] Prescott | Computational evaluation of mitral valve repair with MitraClip[END_REF].

May-Newman and Yin created the first anisotropic hyperelastic law for the mitral leaflet to fit specifically their experimental data from equibiaxial testing of porcine valves [START_REF] May-Newman | Biaxial mechanical behavior of excised porcine mitral valve leaflets[END_REF], [START_REF] Newman | A constitutive law for mitral valve tissue[END_REF]. It was reused several times and adapted by Weinberg and Mofrad to help numerical convergence at low strain [START_REF] Weinberg | A finite shell element for heart mitral valve leaflet mechanics, with large deformations and 3D constitutive material model[END_REF]:

W = c 0 exp(c 1 (I 1 -3) 2 + c 2 (I 1/2 4 -1) 4 ) -1 + c 3 (I 1 -3)
where I 1,4 are the invariant of deformations (values in table 2.4). Table 2.4: May-Newman and Yin coefficients fitted from experimental data [START_REF] May-Newman | Biaxial mechanical behavior of excised porcine mitral valve leaflets[END_REF], with the additional Neo-Hookean component [START_REF] Weinberg | A finite shell element for heart mitral valve leaflet mechanics, with large deformations and 3D constitutive material model[END_REF].

Prot and co-workers develop models similar to the one of May-Newmann and Yin. They are supposedly more robust and stable. When tested on a mitral valve using 3D elements, the model presents a large coaptation area [START_REF] Prot | Transversely isotropic membrane shells with application to mitral valve mechanics. Constitutive modelling and finite element implementation[END_REF]. Similar coaptation was observed with leaflets reconstructed from echocardiographic images [START_REF] Leonardo | Mitral Valve Finite-Element Modelling from Ultrasound Data: a Pilot Study for a New Approach to Understand Mitral Function and Clinical Scenarios[END_REF]. Their latest model is [START_REF] Prot | On modelling and analysis of healthy and pathological human mitral valves: Two case studies[END_REF]:

W = µ(I 1 -3) + c 0 [exp(c 1 (I 1 -3) 2 + c 2 (I 4 -1) 4 ) -1]
where I 1,4 are the invariants of the deformation (values table 2 Gao and co-workers developed another strain energy density function for the mitral leaflets [START_REF] Gao | A finite strain nonlinear human mitral valve model with fluid-structure interaction[END_REF]:

W = c 1 (I 1 -3) + a f 2b f (exp[b f (I * f -1) 2 ] -1)
where I * f = max(I 4 , 1),I 1,4 being the invariants of deformation (values in table 2.6). This model is simple, and offer good coaptation area and closure [START_REF] Gao | A Coupled Mitral Valve -Left Ventricle Model with Fluid-Structure Interaction[END_REF]. Acchording to [START_REF] Cai | Some Effects of Different Constitutive Laws on FSI Simulation for the Mitral Valve[END_REF], it is better than the models of May-Newman and Yin and Prot [START_REF] Cai | Some Effects of Different Constitutive Laws on FSI Simulation for the Mitral Valve[END_REF]. [START_REF] Gao | A finite strain nonlinear human mitral valve model with fluid-structure interaction[END_REF], defined from experimental data [START_REF] Wang | Finite element modeling of mitral valve dynamic deformation using patient-specific multi-slices computed tomography scans[END_REF].

c 1 (kPa) a f (
The model of Holzapfel, Gasser and Ogden is known to describe hyperelastic arterial tissues accurately [START_REF] Gasser | Hyperelastic modelling of arterial layers with distributed collagen fibre orientations[END_REF]. Two fibres orientations are modelled with a pseudo-random fibre distribution. It was applied to the mitral leaflets using optimized parameters [START_REF] Wang | Finite Element Modeling of Mitral Valve Dynamic Deformation Using Patient-Specific Multi-Slices Computed Tomography Scans[END_REF]. Finite element simulations with shells and 3D elements using this model result in good coaptation of the leaflet and large closure area [START_REF] Mao | Fully-Coupled Fluid-Structure Interaction Simulation of the Aortic and Mitral Valves in a Realistic 3D Left Ventricle Model[END_REF], [START_REF] Kong | Finite Element Analysis of Tricuspid Valve Deformation from Multi-Slice Computed Tomography Images[END_REF], [START_REF] Rausch | Mechanics of the mitral valve: A critical review, an in vivo parameter identification, and the effect of prestrain[END_REF]- [START_REF] Caballero | Transapical mitral valve repair with neochordae implantation: FSI analysis of neochordae number and complexity of leaflet prolapse[END_REF]. In the model, a parameter describes the angle between the fibres prefered direction. Using an optimization algorithm, a value was obtained but had no physical meaning. A model similar to HGO, but with one fibre alignment instead of two, was proposed by Lee and co-workers [START_REF] Lee | A high-fidelity and micro-anatomically accurate 3D finite element model for simulations of functional mitral valve[END_REF]. This model represents the micro-structure involved and provides good coaptation of the leaflets and closure of the valve. [START_REF] Lee | A high-fidelity and micro-anatomically accurate 3D finite element model for simulations of functional mitral valve[END_REF], [START_REF] Lee | On the effects of leaflet microstructure and constitutive model on the closing behavior of the mitral valve[END_REF] hyperelastic anisotropic with fibre splay Some models aim to define the mechanical response of the leaflet using the micro-structure of the leaflets, with accurate fibre splay description. These models are more compliant, perform better for general deformations but are more computationally expensive [START_REF] Bischoff | Continuous versus discrete (invariant) representations of fibrous structure for modeling non-linear anisotropic soft tissue behavior[END_REF]. On the plus side, they can represent the complex mechanical response of biological material, such as fibre recruitment or the partial viscoelasticity of the mitral leaflets (and fibre locking).

Kunzelman, Einstein and co-workers were the first to develop such a model for the mitral valve [START_REF] Einstein | Nonlinear Acoustic Analysis of the Mitral Valve[END_REF]. The model is defined using isochoric Cauchy strain instead of invariant deformation. It uses each fibre response to the deformation, supposing that the fibres are uniformly distributed around a preferred direction.:

S = pJC -1 + 2J -2/3 DEV ∂ W ∂ C ∂ W ∂ C = αId + π/2 -π/2 S f R(θ)A ⊗ Adθ with DEV[.] = [.] -1 3 ([.] : C) C-1
the deviatoric projection operator, J = det(F ) the Jacobian, the isochoric variable denoted by (.) are computed from the isochoric deformation tensor, F = J -1/3 F . α represent the isotropic matrix stiffness. The second term of the strain energy density is the sum of all fibre's effect. p is the pressure enforcing quasi-incompressibility, defined as p = κ(J -1), κ the bulk modulus (see annexe for details). R(θ) is the distribution of the fibre around the preferred direction:

R(θ) = 1 σ √ 2π exp -(θ -µ) 2 2σ 2
where µ is the mean preferred direction, σ the standard deviation and θ the splay angle of the fibre. Each fibre stiffness is defined using an exponential law:

S f = F 0 exp F 1 2 (A(θ)CA(θ) -1) -1
where F 0 , F 1 are fibre parameters, A is the fibre orientation, C is the Green strain. The model is computationally expensive due to the required integration on θ to compute each fibre contribution. It is 18 times more computationally expensive than standards hyperelastic formulation. An alternative formulation based on modified splay invariants of deformation was developed to reduce this CPU time. Such simplification is possible because the fibre interaction is neglected at the microstructural level. Using this method, the CPU time was reduced by a factor of 2.6 with a negligible loss of accuracy [START_REF] Freed | Invariant formulation for dispersed transverse isotropy in aortic heart valves: An efficient means for modeling fiber splay[END_REF].

F 0 (kPa) F 1 σ (°) α (kPa) κ (
The model was used in several studies [START_REF] Kunzelman | Fluid-structure interaction models of the mitral valve: function in normal and pathological states[END_REF], [START_REF] Toma | High-Resolution Subject-Specific Mitral Valve Imaging and Modeling: Experimental and Computational Methods[END_REF], [START_REF] Toma | Fluid-Structure Interaction and Structural Analyses Using a Comprehensive Mitral Valve Model with 3D Chordal Structure[END_REF], [START_REF] Toma | Fluid-Structure Interaction Analysis of Papillary Muscle Forces Using a Comprehensive Mitral Valve Model with 3D Chordal Structure[END_REF], [START_REF] Einstein | The relationship of normal and abnormal microstructural proliferation to the mitral valve closure sound[END_REF], [START_REF] Einstein | Nonlinear Acoustic Analysis of the Mitral Valve[END_REF], [START_REF] Einstein | Haemodynamic determinants of the mitral valve closure sound: A finite element study[END_REF], [START_REF] Toma | Fluid-Structure Interaction Analysis of Ruptured Mitral Chordae Tendineae[END_REF]. The application range from simple shape leaflet to fluid-structure interaction of a µ-CT scanned valve during a cardiac cycle. It proves that a microstructural approach is doable and provides good results in both structural and fluid-structure interaction simulation. So far, it is the most efficient method to model continuous fibre orientation.

Sacks and co-workers made a similar model [START_REF] Fan | Simulation of planar soft tissues using a structural constitutive model: Finite element implementation and validation[END_REF]. It uses phenomena at the microscale to define the macroscopic behaviour of the leaflet. The splay, but also the progressive recruitment and stretching of the collagen fibres are taken into account. It is to date the most accurate model for the mitral leaflet and is representative of the microstructures involved. In addition, a simplification is available to limit the complexity of the recruitment process and limit the CPU time used [START_REF] Fan | Simulation of planar soft tissues using a structural constitutive model: Finite element implementation and validation[END_REF]. This model was used several times on reconstructed geometries from imaging [START_REF] Khalighi | Development of a Functionally Equivalent Model of the Mitral Valve Chordae Tendineae Through Topology Optimization[END_REF], [START_REF] Sacks | On the Simulation of Mitral Valve Function in Health, Disease, and Treatment[END_REF], [START_REF] Lee | On the effects of leaflet microstructure and constitutive model on the closing behavior of the mitral valve[END_REF], [START_REF] Drach | A comprehensive pipeline for multi-resolution modeling of the mitral valve: Validation, computational efficiency, and predictive capability[END_REF].

Isotropic modelling of the mitral leaflets

Several isotropic models were implemented and tested for the mitral valve leaflets and chords. We will present each model used and their results in terms of equibiaxial response and valve closure.

Definition of the mechanical law used

The laws used are isotropic, but the leaflet presents an anisotropic response to biaxial testing. Therefore, we choose to focus on the circumferential behaviour of the leaflets to fit our model. Hooke's law, Mooney-Rivlin law and Sussman-Bathe law are detailed in this section.

Linear elastic model for the leaflets

To model the mitral leaflets we first used a isotropic Hooke's law. The associated model can be used with large displacements, but is restricted to small strain. For detail concerning this matter, see Section 6.6.1 in [101]. We have, S = C : , where S is the second Piola-Kirchhoff stress tensor, C is the constant elasticity tensor, and the Green-Lagrange strain. Elasticity tensor is defined as:

C ijrs = λδ ij δ rs + µ(δ ir δ js + δ is δ jr )
, δ being Kronecker delta, µ and λ the Lamé parameters. they can be defined from Young's modulus E and Poisson's ratio ν:

λ = Eν (1 + ν)(1 -2ν) µ = E 2(1 + ν)
The Young's modulus and the Poisson's ratio were set 3 MPa and 0.45, respectively. The bulk modulus is defined with K = E/3(1 -2ν), ν being the Poisson's coefficient. Its corresponding value is 1 × 10 7 Pa, which is similar, but about 2 times smaller than the value used by Kunzelman's group [START_REF] Toma | Fluid-Structure Interaction and Structural Analyses Using a Comprehensive Mitral Valve Model with 3D Chordal Structure[END_REF].

Mooney-Rivlin model for the leaflets

In the software ADINA used for the study, the Mooney-Rivlin law has an extended formulation. The corresponding strain energy density function is defined as:

W = i+j<=3 C ij (I 1 -3) i (I 2 -3) j
With I 1..3 the invariant of deformation, and C ij the coefficients. The model is used in large displacement and large strain. The coefficients were fitted to equibiaxial data obtained in [START_REF] May-Newman | Biaxial mechanical behavior of excised porcine mitral valve leaflets[END_REF]. Parameters were obtained for the anterior leaflet in the circumferential direction. Coefficients for the law were obtained using ADINA's fitting method. Briefly, the error between the experimental data and the response from the model is computed and squared. The error is minimized with Gaussian Elimination (GE) or Singular Value Decomposition (SVD). Gaussian Elimination performed better to minimize the error, but the material law obtained produced decreasing stress at low strain. This behaviour is non-physical and leads to numerical difficulties. The SVD algorithm produces a non-negligible error but results in more physical stress-strain laws. The values obtained were similar to a previous study [START_REF] Abedini | Numerical simulation of the impact of chorda tendinea rapture on the mitral valve leaflets stress by finite element method[END_REF]. The same values were used to compare the results (table 2.7). The bulk modulus was set to 4.5 × 10 3 KPa, which is similar to the value used previously in the literature [START_REF] Toma | Fluid-Structure Interaction and Structural Analyses Using a Comprehensive Mitral Valve Model with 3D Chordal Structure[END_REF].

Coefficients

Sussman-Bathe model for the leaflets

The Sussman-Bathe model was created recently to more accurately fit complex uniaxial mechanical test in tensile and compression [START_REF] Sussman | A model of incompressible isotropic hyperelastic material behavior using spline interpolations of tension-compression test data[END_REF]. It is defined as:

W D = w(e 1 ) + w(e 2 ) + w(e 3 )
with e i = ln(λ i ) the principal or logarithmic strain, λ i the stretch, and w the spline interpolation of w obtained from experimental data. Cauchy stress τ along each axis is defined as:

∀i, τ i = ∂W ∂e i + p = w (e i ) + p
Supposing an incompressible deformation of an isotropic material under uniaxial deformation along axis 1, we have e 1 + e 2 + e 3 = 0 and τ 2 = τ 3 = 0 which leads to τ (e) = w (e) -w (-e/2) by combination of the equation on the 3 axis. The derivative of the strain energy density function can be expressed as:

w (e) = ∞ k=0 τ e - 1 2 k
The derivative for tensile strain is defined from stress under tension (e > 0) and compression (e < 0). Therefore the function w needs both uniaxial tests to be unique (i.e. defined). In practice, Cauchy stress and principal strain in compression and tension are interpolated using cubic splines. A cubic spline interpolation representing w is defined from the newly interpolated Cauchy stress. The number of segments to define the splines is chosen according to an estimation of the stress error. Once interpolated, w can be integrated exactly to define w. The constant of integration is fixed using the stress-strain experimental data. For more information see [START_REF] Adina R & D | Theory and Modeling Guide, Volume I: ADINA Solids & Structures[END_REF], section 3.8.

To use this model, uniaxial tensile and compression data must be provided. Both tests from the same set of leaflets were not found in the literature. A conversion from equibiaxial testing was used to determine the leaflet behaviour [START_REF] Adina R & D | Theory and Modeling Guide, Volume I: ADINA Solids & Structures[END_REF]. To obtain compression data, the material was supposed to be incompressible and that e u = -2e b with τ u = -τ b . (.) u refers to the uniaxial compression mechanical test, (.) b refers to the equibiaxial mechanical test. For the uniaxial tensile test, we supposed that

e u = -2e b with τ u = -τ b .
The Sussman-Bathe material law fits closely to experimental data and is sensitive to experimental noise and sampling. Initial equibiaxial data were sampled using Plot Digitizer and smoothed using Tikhonov regularization (or ridge regression, figure 2.9). Once smoothed, the data were converted to equivalent uniaxial compression and tensile test and then used in ADINA to define the spline coefficients.

Equibiaxial test results

A numerical equibiaxial test was implemented in the software ADINA to test the laws used (figure 2.8). A 1 cm square of thickness 0.8 mm was modelled using shells elements. Linear ramp of displacement was imposed from 0 to 0.42 cm, leading to 50% of Green strain ( = 0.5(λ 2 -1)) . Opposed boundaries are fixed in the displacement directions but are allowed to translate in the other directions. 100 timesteps of 1 × 10 -2 s were used to model the equibiaxial mechanical test. 16 shells are used. Large displacement and large strain are used when available. Results obtained are independent of time steps and mesh size. Stress and strain response of the numerical equibiaxial test are presented figure 2.9. The elastic linear law has an equivalent stiffness of the taut leaflets, but the deformation that should have happened at a small strain is totally ignored. The Mooney-Rivlin model fits the experimental data quite well at low strain. However, the polynomial law can not follow the steep rise with increasing strain, resulting in poor realism at higher strain. This observation occurs for all sets of parameters tested. Sussman-Bathe models fit the experimental data well at low strains. They also capture the stress augmentation at higher strain. However, a difference appears between experimental data and the Sussman-Bathe model. It increases with increasing stress. We supposed this difference to be linked with the simplification used to convert equibiaxial data into equivalent compression and tensile data. The Sussman-Bathe model of the posterior leaflet show similar results and has the same default.

Closure test and results

An additional test was made with the material laws on a shell model of the mitral valve-this additional test help to understand and identify the key characteristics required for the valve to close. The geometry of the valve was defined using a Fortran and Python code (chapter 1).

Chordae tendineae structural law

The mitral chords were modelled using an isotropic hyperelastic law obtained from experimental data (figure 2.3, [START_REF] Kunzelman | Mechanical properties of basal and marginal mitral valve chordae tendineae[END_REF]). Only marginal chords are used. The Young modulus after recruitment is around 30 MPa. The chordal cross-section was set to 0.6 mm 2 acchording to similar simulation [START_REF] Lau | Numerical simulation of mitral valve function[END_REF]. When considering uniaxial deformation, the corresponding stiffness modulus can be defined as k = EA/L where k is the stiffness, E Young's modulus, A the cross-section, and L the length. The chord length range from 1 to 1.8 cm in our model. The corresponding stiffness varies from 1.0 kN to 1.8 kN. It is similar to the experimental data obtained in the literature (table 2.1).

Finite element model

The finite element model used comprises a fixed annulus and fixed papillary muscles connected with marginal chords. A systolic pressure ramp uniformly distributed on the leaflets is applied from 0 to 0.2 s. Time step is set to 1 × 10 -3 s. Rayleigh damping is used, with a beta factor of 10 -3 [START_REF] Adina R & D | Theory and Modeling Guide, Volume I: ADINA Solids & Structures[END_REF]. The model is the same for the 3 tests, except for the material laws that differ. The results are observed at 0.2 s when the maximal pressure is applied (figure 2.10a, 2.10b and 2.10c). All the data are presented in CGS (cm, gr, s) units. Nodal gap distance is in cm and stress is in dyn.cm -2 (1 dyn.cm -2 = 0.1 Pa).

Results

Using linear elastic material results in a high-stress area with values above 0.5 MPa. The shape of the leaflets is similar to the ones obtained with hyperelastic materials. Contact between the leaflets occurs, but the area of coaptation is minimal. The leaflets are partially in contact, and some aperture can be seen at the commissural area and on the posterior leaflet. Overall the closure obtained with the linear model is only partial.

The Mooney-Rivlin model reaches lower maximal stress than the linear model, which is more distributed on the leaflet. The middle part of the anterior leaflet undergoes larger deformation and stress when compared to the linear model. Leaflets coapt and the coaptation area is more important than in the linear model. No residual gap can be seen between the leaflets. The anterior leaflet has a major bump above the annulus plane of 1 cm. This bump is called a prolapse and can be seen as a pathology. The Mooney-Rivlin model performs well during closure except for the prolapse, and the resulting mitral valve does not leak.

The Sussman-Bathe model undergoes lower stress than the Mooney-Rivlin model, in particular at the anterior leaflet. The peak of stress is obtained at the insertion site of the chords, as in the Mooney-Rivlin model. The area of coaptation is similar to the one obtained by the Mooney-Rivlin model. This model performs well in terms of closure and valve shape during closure.

The use of hyperelastic material law results in good closure of the valve. The Mooney-Rivlin model has a significant prolapse, but not the Sussman-Bathe nor the linear model. The prolapse relates to the small increase in stiffness under large strain in the Mooney-Rivlin function (figure 2.9). Large displacements must occur at low stress to obtain a coherent leaflet closure, with an important closure area. However, additional strain must be small when larger stress is reached. No linear material can perform well on both small strain and large stress. Mooney-Rivlin is not realistic either at high stress. From a numerical point of view, both hyperelastic cases were difficult to compute. Hyperelasticity introduces additional non-linearity in the model. Non-linearity increase the difficulty for the solver to find equilibrium. This is especially true with the Sussman-Bathe model. The steep stress-strain relationship explains this difficulty partially. The contact between the leaflets adds additional nonlinearity. To capture and compute this contact, an automatic time-stepping algorithm was used. This method decreases the timestep when convergence is not reached and increases it if several attempts at a given timestep are successful. With the Sussman-Bathe model, the automatic time-stepping algorithm reduced the time step from 10 -3 s to an unrealistic 10 -6 s during l simulation. This increased the CPU time from few minutes to an hour in the present model when convergence was obtained. Additional damping effect was added to soften the resolution. It helps slow down the structure that would otherwise oscillate at high frequency. The dynamic of the valve was not realistic, with a damping coefficient large enough to allow convergence. It is not conceivable to have the leaflets material penalizing the timesteps or the dynamic response, particularly when considering the very few numbers of elements (less than a thousand). As a result, and although the Sussman-Bathe law produced a coherent mitral valve closure, other materials were studied to prevent this drastic time step reduction.

Anisotropic modelling of the mitral leaflets

The isotropic models of the valve were presented in the previous section. The Sussman-Bathe material law fits the imposed stress-strain relationship in the circumferential direction but leads to difficulties during the simulation. These could not be explained entirely by the hyperelastic nature of the material. Instead, it suggests a misunderstanding of the physic modelled. This section details the notion of material stability, its definition for practical uses, and its application to our model. We conclude with the presentation of some anisotropic models that could be used for the mitral leaflets.

Definition of the material stability

Material stability refers to the physical stability of a material and can be tested from its constitutive equation. It is not related to numerical discretization. The use of a stable material does not guaranty that the finite element model will be stable. However, a material that does not respect the stability condition can induce numerical instabilities [START_REF] Ericksen | Equilibrium of bars[END_REF]. This condition is different to the buckling instability.

Material stability is commonly a condition linked to the constitutive law. One of such rule is "stress increases with strain" (also called Baker-Ericksen inequality). Another limitation is the one defined by Hadamard and Legendre. Starting from equilibrium and propagation of elastic waves in the material, an inequality can be written ( [START_REF] Hadamard | Leçons sur la propagation des ondes et les équations de l'hydrodynamiques[END_REF], issue 270). Rewritten using Truesdell and Noll formulation [START_REF] Truesdell | The Non-Linear Field Theories of Mechanics[END_REF], it reads :

B A ijkl H ij H kl dv ≥ 0 where A ijkl = ∂ 2 W
∂F kl ∂F ij , W being the strain energy density function (quasiconvex), F the deformation gradient, and H the displacement gradient. This is implied by the standard mathematical definition of strong ellipticity, as defined in [START_REF] Chagnon | Hyperelastic Energy Densities for Soft Biological Tissues: A Review[END_REF]:

1 J F pr F qs ∂ 2 W ∂F ir ∂F js n p n q m i m j > 0
with m, n = 0, F the deformation gradient and J = det(F ) the Jacobian.

The strong ellipticity condition also verifies Baker-Ericksen's rule. A sufficient condition was declined for a certain type of anisotropic material [START_REF] Walton | Sufficient conditions for strong ellipticity for a class of anisotropic materials[END_REF], but strong ellipticity is in general too complex to be used as such. The polyconvexity of the function is often preferred [START_REF] Ball | Convexity conditions and existence theorems in nonlinear elasticity[END_REF].

A practical stability condition can also be defined using stress and strain expression. Drucker's postulate states that the increments of stress due to an incremental deformation can only increase [START_REF] Hill | A general theory of uniqueness and stability in elastic-plastic solids[END_REF], [START_REF] Romanov | The drucker stability of a material[END_REF]: dσ : d ≥ 0 with dσ the Cauchy stress increment corresponding to d the Cauchy Green strain increment. The Drucker criterion for stability looks similar to the deformation term in the Clausius-Duhen inequality, based on the second thermodynamic law. In a constant and uniform temperature field, it reads: -ρ ψ + σ : D ≥ 0 with ψ the Helmholtz energy, D the strain rate tensor, verifying Ė = F T DF . For reversible transformation at constant volume, the variation of Helmholtz energy is 0, which leads to a dynamic condition similar to Drucker's postulate.

Hill-Drucker's postulate is common among plasticity and is the default criterion used in the ADINA software [START_REF] Adina R & D | Theory and Modeling Guide, Volume I: ADINA Solids & Structures[END_REF], [START_REF] Zheng | On the predictive capability and stability of rubber material models[END_REF]. Moreover, this criterion is useful in practice, as material stability can be tested easily and numerically.

Practical definition of the stability criterion

The approach detailed in this section is a generalization of the isotropic stability condition to the anisotropic case [START_REF] Adina R & D | Theory and Modeling Guide, Volume I: ADINA Solids & Structures[END_REF], [START_REF] Zheng | On the predictive capability and stability of rubber material models[END_REF].

Let us consider a unit cube of an hyperelastic incompressible material under applied loads R i , i = 1..3 2.11. The loads are homogeneous and independent of the displacement. The induced displacement is noted u i , i = 1..3. Stretches are noted λ i = 1 + u i , and true strains e i = ln(λ i ) for i = 1..3. The material is hyperelastic, therefor a strain energy density function W exist and can be written as:

Ŵ = W + m(e 1 + e 2 + e 3 )
with m the Lagrange multiplier. The material studied being incompressible, we have e 1 + e 2 + e 3 = 0. It results that the variation of the strain energy density function, also called internal work, are:

δ Ŵ = i ∂W ∂e i + m δe i + δm(e 1 + e 2 + e 3 )
Considering the variations of the external work, we have:

δR = i R i δu i = i R i λ i δe i
because the loads are homogeneous and independent of the displacements. When studying the material at static equilibrium, the internal work balance the external work. This is true for any displacement δe i , thus we have:

∀i, ∂W ∂e i + m = R i λ i e 1 + e 2 + e 3 = 0 (2.1)
The first equation of 2.1 can be differentiate:

∀i, ∆ ∂W ∂e i + m = ∆ (R i λ i ) ∀i, R i λ i ∆e i + λ i ∆R i = ∆m + ∆ ∂W ∂e i ∀i, j ∂ 2 W ∂e i ∂e j -R i λ i ∆e i + ∆m = λ i ∆R i
which can be written in a matrix form:

    KEE 1 1 1 1 1 1 0         ∆e 1 ∆e 2 ∆e 3 ∆m     =     λ 1 R 1 λ 2 R 2 λ 3 R 3 0     (2.2)
with KEE defined as the symmetric matrix:

KEE =     ∂ 2 W ∂e 2 1 -R 1 λ 1 ∂ 2 W ∂e 1 ∂e 2 ∂ 2 W ∂e 1 ∂e 3 sym ∂ 2 W ∂e 2 2 -R 2 λ 2 ∂ 2 W ∂e 2 ∂e 3 sym sym ∂ 2 W ∂e 2 3 -R 3 λ 3    
The Hill's stability condition states that the work produced must be positive:

i ∆u i ∆R i = ∆e 1 ∆e 2 ∆e 3   λ 1 ∆R 1 λ 2 ∆R 2 λ 3 ∆R 3   ≥ 0
With the previous expression of the λ i R i vector, we obtain the following matrix form :

∆e 1 ∆e 2 ∆e 3 ∆m KP P     ∆e 1 ∆e 2 ∆e 3 ∆m     ≥ 0 (2.3)
with:

KP P =     KEE 1 1 1 1 1 1 0    
The inequality 2.3 will be true for any deformation only if the matrix KP P is symmetric positive definite. One equivalent statement is for all the eigenvalues of KP P to be positive, ie:

KP P positive definite ⇐⇒ min λ∈Sp(Kpp) λ > 0
The stability criterion of our material is now based on the positiveness of eigenvalues, which can be computed numerically and plotted as a function of the deformation. To obtain the eigenvalues, the KEE matrix must have numerical values. It is composed of derivatives of the strain energy density function, the stretch and the loads. The derivatives can be obtained for a given hyperelastic material and the deformation. The stretches are easily computed from the deformation. Some hypotheses must be formulated to express the loads, as they are not known in general. Materials are fully characterized using their response to equibiaxial, uniaxial and shear mechanical test. Thus, we supposed three different types of loads for our material.

Under uniaxial testing along the axe i, we have R j = 0, j = i, and e i = -e j -e k j = k = i due to incompressibility. Using combinations of equation 2.1, the Lagrange multiplier can be replaced to obtain an expression for the uniaxial load:

R i λ i = ∂W ∂e i - 1 2 ∂W ∂e j + ∂W ∂e k
Under equibiaxial testing in the plane orthogonal to axe i, we have R i = 0, and e j = e k k = j = i. Using combinations of equation 2.1, the Lagrange multiplier can be replaced to obtain an expression for the equibiaxial load:

∂W ∂e j - ∂W ∂e i = R j λ j j = i
Under shear testing in the plane orthogonal to axe i, we have R i = 0 and e k = 0 k = i. Using combinations of equation 2.1, the Lagrange multiplier can be replaced to obtain an expression for the shear load:

∂W ∂e j - ∂W ∂e i = R j λ j j = i
For any of the three deformation proposed, the matrix KP P can be expressed only with the first and second derivative toward the true strain. However, most of the strain energy density function are defined from the invariant of deformation. In the code proposed, the matrix KP P is defined using the invariant I 1..5 with the derivative composition:

∀i ∂W ∂e i = j ∂W ∂I j ∂I j ∂e i and ∂ 2 W ∂e 2 i = j ∂ ∂I j ∂W ∂e i ∂I j ∂e i
The previous formula is obtained when the invariant are expressed using stretch components. The true strains are defined as e i = ln(λ i ). As a result, the matrix KP P is a function of the deformation and its variation in regards to the invariant:

KP P ∂W ∂I i i , ∂ 2 W ∂I i ∂I j i,j , (e i ) i
Using the strains (e i ) i , it is possible to compute the eigenvalues of the matrix and check if the material is stable according to Hill's stability condition. So far, no hypotheses have been made considering the spatial organisation of the material. If the material were incompressible isotropic, only 3 tests would be necessary (one shear, one equibiaxial and one uniaxial). In addition, the uniaxial test along axis i would verifies e j = e k = -e i /2, j = k = i. Using isotropic assumption, only one strain in a given direction could monitor the three mechanical tests (equibiaxial and shear tests are by definition monitored using only one strain). The minimal eigenvalue of the matrix KP P would be determined using one strain only. The plot of minimal eigenvalues for each mechanical test along principal strain represents the stability of the material. This graphical representation of stability is used in ADINA. Our model was plotted for comparison when the isotropic assumption was used. A comparison of Mooney-Rivlin material is given for illustration (figure 2.12). The material has negatives eigenvalues for biaxial testing for strain larger than 0.25. As a result, a finite element model using this material must overcome additional difficulties to represent this non-physical behaviour.

Additional simplification

We note that the size of the problem 2.3 can be reduced using the equation concerning incompressibility, ∆m = ∆e 1 + ∆e 2 + ∆e 3 = 0, leading too ∆e 3 = -∆e 1 -∆e 2 . The matrix KP P is reduced into a symmetric matrix of size 2. The eigenvalues are then obtained efficiently by solving the roots of the second-order polynomial characteristic. The shear test can be further simplified, leading to a size

Difficulties to obtain stability curves for anisotropic material

The method presented here has two defaults for an anisotropic material. First, anisotropic materials are not fully characterized with only tree mechanical tests because they do not behave identically in any direction. Ideally, 3 material directions should be tested. The fibres in the mitral valve have one preferred direction. This direction can be expressed as a linear composition of 3 vectors aligned with the axis of the Cartesian system. We consider that if the material were stable for the 3 unit vectors, then the material with any fibre orientation would also be stable as a composition. The justification for this hypothesis is that the stability of the material can always be decomposed into the stability of each fibre component. Another way to look at it is to consider the material rotation. Thus the 3 directions are tested under uniaxial, biaxial and shear mechanical tests, which would also prove the stability for any material. To conclude, to test the stability of the anisotropic material in any fibre direction, it is sufficient to prove its stability with 3 fibres aligned on the coordinate system.

The second default is to monitor the minimal eigenvalue under uniaxial test along axis i, the simplification e j = e k = -e i /2, j = k = i was applied. To prevent this, no relation is supposed between the principal strains. The first principal strain e i is imposed and used as abscissa for the uniaxial curve. The second strain e j can vary between the minimum and maximum range given for the first principal strain. The third strain e k is obtained using incompressibility. The lowest eigenvalue obtained for any e j is used.

The anisotropic stability method presented here was implemented in Octave, a free, compatible and complete alternative to Matlab.

Stability results of the material used.

The simulations encountered numerical difficulties with the isotropic hyperelastic material presented in the section 2.13. These have motivated the development of stability criteria to test the mechanical laws. The isotropic version of the stability criterion was applied to the hyperelastic material previously used (figure 2.13). We observed that the minimal eigenvalue obtained could be negative for all materials before the strain reached 0.5. Thus, the hyperelastic material can not be considered stable. This is counterintuitive, especially for the Sussman-Bathe materials that fit closely to experimental data. The physical mitral leaflets used for the experiment are real and thus represent a stable material. The method used must have altered the stability of the leaflets. We suppose that the sampling of the data, along with the conversion method from equibiaxial to uniaxial traction and compression, was a sensitive step that has altered the mechanical response. The Mooney-Rivlin material is known not always to be stable. It appears that the coefficient obtained [START_REF] Abedini | Numerical simulation of the impact of chorda tendinea rapture on the mitral valve leaflets stress by finite element method[END_REF] leads to unstable material, despite the realistic data used to fit the model.

Although the instability limits the study of the mitral leaflets, it also explains the numerical difficulties encountered during the previous simulations (section 2.13). The numerical difficulty could be overcome using a non-realistic but smaller time step. The instability was also softened by introducing high viscosity and a Rayleigh damping. The deformation's speed, increased by the material instability, was lowered by the damping effect. The result was coherent in static but is not a good option for dynamic studies.

To obtain a coherent dynamic simulation, we must focus on a new law for the mitral valve that will respect the stability criterion.

Application to anisotropic strain energy density function

Most of the material made for the mitral leaflets are anisotropic hyperelastic (section 2.1.2.3). To model the mitral leaflets, several strain energy density functions were used. The May-Newman & Yin model [START_REF] Weinberg | A finite shell element for heart mitral valve leaflet mechanics, with large deformations and 3D constitutive material model[END_REF], the model of Prot and co-workers [START_REF] Prot | On modelling and analysis of healthy and pathological human mitral valves: Two case studies[END_REF] and the model of Gao and co-workers [START_REF] Gao | A Coupled Mitral Valve -Left Ventricle Model with Fluid-Structure Interaction[END_REF] were used. The functions and the coefficients used were detailed previously in section 2.1.2.3. The derivatives for each law towards the invariants I 1..5 were implemented in the anisotropic stability code. The three directions N 1 = (1, 0, 0), N 2 = (0, 1, 0) and N 3 = (0, 0, 1) are used as fibre orientations to test the anisotropic response of the material. For each fibre orientation, 3 curves are obtained, corresponding to uniaxial, biaxial and shear stability.

Initial results

Our first implementation of the laws results in negatives eigenvalues. Thus the material was nonphysical (figure 2.14,a). The mechanical test resulting in negative eigenvalues are the equibiaxial test orthogonal to the fibre direction and the uniaxial test orthogonal to the fibre. From a mechanical point of view, these tests compress the fibre (figure 2.14,b). When the fibre is compressed, the term (I 4 -1) 2 contributes to stress. However, the collagen fibres are crimped and fold under compression. Folding offers very little resistance compared to stretching and should be neglected in the material's mechanical response. To do so, the modified invariants I f = max(I4, 1) will be used. This hypothesis can lead to convergence difficulties. Indeed, the strain energy density function should be C 2 to optimize the equilibrium search. Despite this, using the modified invariants is a good method to neglect the fibre behaviour under compression.

The stability curves obtained for the three mechanical laws are presented in the figure 2.15. For each mechanical law, the parameters for the anterior and posterior leaflets were tested. All the material laws respect the stability criterion. As a result of the ignored fibre compression, the stability curves behave differently in compression and extension. The addition of a Neo-Hookean term raises all the curves above the horizontal axis. This observation is distinct for Gao mechanical law and the law of Prot and co-workers for the posterior leaflet (left). It is also true for the model of May-Newman and Yin, although less visible with a Neo-Hookean term of 10 -5 Pa. Null eigenvalues are observed for the model of Prot with the anterior leaflet at 0 strain, but negligible Neo-Hookean term could be added to overcome the difficulty. Identical curves obtained for different directions can be observed. They prove that the material response is the same for all directions orthogonal to the fibre, i.e. the models describe a transversely isotropic material.

Implementation in ADINA and results

Our finite element models were run using the software ADINA. The hyperelastic anisotropic materials chosen were not available in ADINA. Subroutines in Fortran were used to define the material law. After their implementation, the materials were tested to verify the mechanical response and the implementation of our model.

Implementation of the material laws

The Fortran subroutine used to implement the material law is divided into 4 sections. The first and last sections are used to import and print information available in the subroutine. The second section computes stress response to strain. The third section computes the constitutive matrix of the material for a given strain. The section can be run separately. The subroutine is run (at least partially) for each solver iteration, for each element and at each integration point. Thus, it must be efficient from a computing point of view.

The 3 materials laws we would like to implement are anisotropic and hyperelastic and neglect the viscoelastic response of the leaflets. They can be defined using a strain energy density function. The standard formulation for this type of material is a Total Lagrangian formulation (TL). This formulation uses the configuration at time 0 to compute the response at time t + dt. In opposition, Updated Lagrangian formulation such as the Updated Lagrangian Hencky or Jaumann (ULH or ULJ) uses the configuration at time t as a reference to compute the response at time t+dt. Both formulations describe the same behaviour mathematically but result in different numerical implementation and complexity. For example, viscoplastic materials are more easily defined using the updated formulation, whereas hyperelastic materials are more easily defined with TL formulation. Both formulations are available in ADINA and can be used in large strain and large displacements.

For large displacements and large strains, the ULH formulation is used by default in ADINA. ABAQUS, another finite element software used in large displacement and large strains, also uses updated formulation for such simulations. It is not clear why the updated formulation should be preferred for hyperelastic material in large strain and large displacement in general. Using the TL formulation results in easier and computationally efficient implementation of the material in this case. The material response is theoretically the same ULH formulation. Therefore, we choose the TL formulation to implement the material subroutine. To force its use in our software, the option "small strain" must be specified in the material model, while large strain and large displacement are specified as kinematic options (see figure 2.16).

With this option, the Fortran subroutine has as input Green-Lagrange strain and outputs the second Piola-Kirchoff stress tensor and its constitutive matrix defined with reference configuration at time 0. Although the small strain option is specified in the material model, large strain and large displacement are used for the kinematic formulation. To ensure this, an additional test was run with the standard Mooney-Rivlin model in large strain large displacement and compared to the subroutine results. The two simulations produce the same results. The definition of the stress and the constitutive matrix from the strain invariants and the strain energy density function derivatives is given in the annexe. The constitutive matrix construction is used to build the stiffness matrix. The stiffness matrix is used for convergence purposes during equilibrium iterations.

Results of the anisotropic material laws

Once implemented, the material laws were tested on numerical mechanical tests and compared to the experimental data used to fit the parameters. The strain energy density functions chosen were detailed in section 2.1.2.3.

The mechanical law of May-Newman and Yin modified by Weinberg and Mofrad [START_REF] Weinberg | A finite shell element for heart mitral valve leaflet mechanics, with large deformations and 3D constitutive material model[END_REF] was tested under numerical equibiaxial mechanical test (figure 2.17). The corresponding experimental data were extracted from the article using the software Plot Digitizer. They represent the average response of 16 porcines mitral valves leaflets. The stress obtained under numerical simulation fits the data almost perfectly. The exact same differences observed between the model and the data are also present in the original article [START_REF] Newman | A constitutive law for mitral valve tissue[END_REF].

The mechanical law of Prot and co-workers was tested under uniaxial mechanical tests (figure 2.18). The corresponding experimental data were extracted from the corresponding article using Plot Digitizer [START_REF] Prot | On modelling and analysis of healthy and pathological human mitral valves: Two case studies[END_REF]. They represent the mechanical response of a 88 years old mitral valve. The model fits the experimental data well. Unfortunately, no equibiaxial experimental data were available for comparison.

The mechanical law of Gao and co-workers was tested under equibiaxial mechanical tests (figure 2.19). The coefficients for the model were obtained from a fit to experimental data obtained by the same group in a previous study (80 years old mitral valve [START_REF] Wang | Finite Element Modeling of Mitral Valve Dynamic Deformation Using Patient-Specific Multi-Slices Computed Tomography Scans[END_REF]). No comparison between the data and the model are available in the literature. The model response to equibiaxial stretch lack accuracy when compared to the experimental data. This is due to the more humble strain energy density function used. Indeed, the term in the exponential has a lower exponent than in the law of Prot. The quick transition induced by the recruitment of fibres requires higher-order on the invariant I 4 to be captured. In addition, only I 4 invariant is included in the exponential. As a result, the radial direction stress (squares) corresponds to the Neo Hookean term.

Discussion and conclusion

The mitral valve has a complex mechanical response to strain. The chordae tendineae and the leaflets were assumed to be the deforming parts of the valve. According to histological studies and microscopic observations, these parts are mainly composed of collagen fibre embedded in a viscous matrix of glycoproteins. The resulting mechanical response to stretch is viscous, anisotropic and hyperelastic.

The study was simplified to anisotropic hyperelastic response, as the viscous effect has a limited influence on the valve function. This hypothesis permits to use strain energy density function to model both the leaflets and the chordae tendineae.

The mechanical response of all chords under elongation is approximately the same. The most efficient method to model their function is to use the same stress-strain relationship and diameter. The model can be independent of the type of chords (marginal, basal or strut). The chords modelling could be more accurate if a strain energy density function and specific parameters were used.

Several isotropic models for the leaflets were tested under numerical equibiaxial and mitral valve closure tests. These studies have shown that the hyperelasticity of the material is a key feature for the coaptation of the leaflets. However, the stress response to strain is steep, and not all hyperelastic laws seem fitted to model the leaflets. The Sussman-Bathe material provided in ADINA was able to model this response to deformation. However, the use of hyperelastic material increases the model's numerical difficulty, and convergence issues occur.

Physical limitations for the material can be derived from wave propagation, energetic balance or even thermodynamic laws. These inequalities can be applied to the material laws to determine if their behaviour agrees with physics or not. An anisotropic test based on Hill's postulate for hyperelastic strain energy density function was developed. The hyperelastic isotropic materials we used did not agree with Hill's postulate and could lead to instabilities. Several anisotropic hyperelastic laws for the mitral leaflets were tested and verified this postulate. The model developed by Gao and co-workers had the larger margin for this stability criteria. To a lesser extent, the model of Prot and co-workers also had some margin.

Three hyperelastic anisotropic models for the mitral valve leaflets were implemented in ADINA. Fortran subroutines were coded for large displacements and large strains in the software. The models were tested with equibiaxial or uniaxial mechanical tests and compared to experimental data. The model of Mofrad and Weinberg, based on a previous formulation from May-Newman and Yin, fits the porcine experimental data. The model of Prot and co-workers is similar and also fits human experimental data. The model of Gao and co-workers has the lowest accuracy when compared to human experimental data.

Chapter 3

Structural simulation of the mitral valve

The mitral valve act as a hydraulic diode, letting fluid flow in one direction and being hermetic when the flow is going backwards. Fluid-structure interactions govern the opening and closure of the valve, but these are much less pronounced during closure than the opening. Coherent simulations of valve closure can thus be obtained by modelling the valve behaviour under a uniform pressure subjected by the fluid on the structure. Besides, they provide a reasonable estimate of the coaptation surface (i.e. surface area of the leaflets where one is in contact with the other) and possible leakage. The second objective of the thesis is to simulate the mitral valve during closure using a structural numerical model.

An introduction to the mitral valve structural simulation

The main difficulties to face when simulating the dynamics of mitral valves are the large displacements and deformations of the leaflets, the contact of a surface onto itself, the non-linear constitutive laws of the valve material and the coherent closure shape obtained under physiological range. From a numerical point of view, the most difficult is the constitutive law. A chapter of this thesis was dedicated to this issue. Moreover, the shape of the mitral valve can be hard to obtain, and another chapter was dedicated to this issue. The contact algorithm is found in most finite element software and can be handled using several algorithms. The remaining difficulties are the large deformation and large displacement of the leaflets.

This section is divided into two parts-the first part details what we can deduce from the structural simulation of the mitral valve. The second part shows the numerical method used to obtain this information.

Interest for the structural simulation of the mitral valve

The structural simulation of the mitral valve is an efficient method to obtain the static position of the valve during systole. Indeed, the fluid effect on the valve is less crucial when the valve has its closed position during systole. As a result, structural simulations are used to work on the coaptation of the mitral valve, which is only dependent on the static closed position. Several medical procedures were tested using structural only simulation (figure 3.1). [START_REF] Morgan | Neochord Placement Versus Triangular Resection in Mitral Valve Repair: A Finite Element Model[END_REF] used a structural model in LS-DYNA to represent the coaptation increase with several medical repair techniques (resection with/without annulus ring and use of the Neochord procedure with/without ring). A complete finite element heart model in SIMULIA (Dassault systèmes) was employed to model the effect of an infarct on the mitral valve [START_REF] Baillargeon | Human Cardiac Function Simulator for the Optimal Design of a Novel Annuloplasty Ring with a Sub-Valvular Element for Correction of Ischemic Mitral Regurgitation[END_REF]. Ultimately, the model was re-used to restore the coaptation of the leaflets using a new annuloplasty ring (figure 3.1b). Other examples of structural simulation medical procedures assessed using structural simulation of the mitral valve can be found in the literature. For example, the MitraClip procedure was used on several simulations and seems to be efficient in restoring the leaflets coaptation [START_REF] Lau | Numerical simulation of mitral valve function[END_REF], [START_REF] Mansi | An integrated framework for finite-element modeling of mitral valve biomechanics from medical images: Application to MitralClip intervention planning[END_REF]. Another interest in the structural simulation of the mitral valve resides in the fact that internal stress is computed. The internal stress represents how much the leaflet has already worked due to the deformation imposed. This stress is used in linear isotropic material to determine a plastic deformation (where the deformation is not reversible) and or the rupture of the material. High stress is also related to material degradation in the mitral valve and shortens the valve's lifetime.

In figure 3.1a, the group of Ratcliffe have found that sufficient coaptation was obtained with all the procedures tested. However, the average internal stress was the lowest with the annulus ring and 2 Neochord. As a result, this method should have the longest lifetime. According to structural simulation, the MitralClip device increases the maximal stress in the mitral valve, which could lead to a deterioration of the tissue where the device is implanted [START_REF] Lau | Numerical simulation of mitral valve function[END_REF]. There is currently no method to assess the internal stress state of the mitral valve in-vivo. To estimate it, in-vivo deformation can be used with material laws obtained from ex-vivo mitral valve ( [START_REF] Sacks | In-Vivo Dynamic Deformation of the Mitral Valve Anterior Leaflet[END_REF] and [START_REF] Newman | A constitutive law for mitral valve tissue[END_REF] for examples). Inverse Finite element analysis can also be used to estimate this stress [START_REF] Krishnamurthy | Material properties of the ovine mitral valve anterior leaflet in vivo from inverse finite element analysis[END_REF]. The development of complex material laws with the use of accurate µ-CT scan mitral valves enable the use of finite element to determine approximate values for the internal stress in the leaflets (figure 3.2).

The group of Sack have developed an accurate model of the mitral valve. The model was used in several studies to understand mitral dynamic and static behaviour. Under physiological pressure, the stress on the anterior leaflet reaches 0.5 MPa, in accordance with published values [START_REF] Kunzelman | Fluid-structure interaction models of the mitral valve: function in normal and pathological states[END_REF], [START_REF] Toma | Fluid-Structure Interaction and Structural Analyses Using a Comprehensive Mitral Valve Model with 3D Chordal Structure[END_REF], [START_REF] Cai | Some Effects of Different Constitutive Laws on FSI Simulation for the Mitral Valve[END_REF].

Numerical methods used

The finite element used for the mitral valve must be compatible with large deformation and large displacements. The model can be discretized using shells or 3D elements, as the information contained in the thickness of the valve is not a requisite feature of the mitral valve response. Both discretizations were used in the literature. The models using the software LS-DYNA are mainly based on the initial work of Einstein, Kunzelman and co-workers [START_REF] Kunzelman | Fluid-structure interaction models of the mitral valve: function in normal and pathological states[END_REF], [START_REF] Einstein | Dynamic finite element implementation of nonlinear, anisotropic hyperelastic biological membranes[END_REF] (figure 3.3a). They use membrane elements to represent the mitral valve. Shells models of the mitral valve are also implemented in ABAQUS 6, with a similar configuration (explicit time scheme with ∆t = 10 -6 s, S4 or S3 first order shells, quasi-static solution [START_REF] Sacks | On the Simulation of Mitral Valve Function in Health, Disease, and Treatment[END_REF], [START_REF] Lee | On the effects of leaflet microstructure and constitutive model on the closing behavior of the mitral valve[END_REF]).

3D elements were also used to model the mitral leaflets. However, the mitral valve is quasiincompressible, and buckling behaviour can appear. A pressure field can be interpolated separately to prevent instability. This formulation is called hybrid (pressure-displacement). Prot and co-workers model the mitral valve using 3D hybrid brick elements in ABAQUS. Quasi-static simulation of a 3D mitral valve was run with two elements in the thickness of the valve [START_REF] Prot | On modelling and analysis of healthy and pathological human mitral valves: Two case studies[END_REF]. Toma and co-workers model a 3D element mitral valve with tetrahedral elements using µ-CT scan geometry in the software IMPETUS [START_REF] Toma | High-Resolution Subject-Specific Mitral Valve Imaging and Modeling: Experimental and Computational Methods[END_REF], [START_REF] Toma | Fluid-Structure Interaction and Structural Analyses Using a Comprehensive Mitral Valve Model with 3D Chordal Structure[END_REF], [START_REF] Toma | Fluid-Structure Interaction Analysis of Ruptured Mitral Chordae Tendineae[END_REF] (figure 3.3b). A 3D elements mitral valve was modelled during a cardiac cycle with a left ventricle contraction using a formulation based on the immersed boundary method (software IBAMR, ∆t = 10 -4 s to 10 -5 s ) [START_REF] Gao | A Coupled Mitral Valve -Left Ventricle Model with Fluid-Structure Interaction[END_REF]. The free software IBAMR was also used with shell elements to model the mitral valve in a tubular straight channel [START_REF] Ma | Image-based fluid-structure interaction model of the human mitral valve[END_REF].

Alternative formulations for the mitral valve exist. For example, an equivalent tension tissue was obtained from a cross fibre definition of the mitral valve [2]. The structural shape obtained under tension reproduces the pattern observed on anatomical images. The structural model was then used in fluid-structure interaction with an immersed boundary method. Additionally, a one-parameter mitral valve was developed to better capture the effect of the mitral valve on the cardiac flow [START_REF] Celotto | Influence of mitral valve elasticity on flow development in the left ventricle[END_REF]. This valve has a parametric geometric definition. A parameter handles the state of the valve (closes/open) and its elastic response and directly modify its shape.

Methods used to model the mitral valve

The finite element software ADINA was chosen to simulate the mitral valve closure. Details on the finite element formulation in large deformation and large strain can be found in [101]. This section details specific methods used for the simulation of the mitral valve. The temporal resolution used for the simulation is first addressed. The elements used are then detailed, and some information on the buckling effect is given. Several methods to anchor medical devices onto the leaflets are then explained. Finally, a damping technique used to absorb numerical vibration is presented.

Temporal resolution

Estimation of the time-step

A first approach to estimate the time resolution is to consider the characteristic time of the structure. The characteristic velocity of a linear material is:

c s = E ρ
with E Young's modulus and ρ the density and c s the characteristic velocity. The structure's characteristic time t s is deduced as:

t s = L s c s with L s the characteristic length.
In the mitral valve's case, the order of magnitude of the density, the post-stiffening Young's modulus and the characteristic length are given in table 3.1. Thus, we deduce c s 100 m.s -1 and the order of magnitude of the solid response t s = 2 × 10 -4 s. As a result, a reasonable time-step for our simulation is ∆t = 1 × 10 -4 s. Elastic waves also appear when the mitral valve is not under systolic load. When unloaded, the equivalent Young's modulus of the mitral valve is lower; thus, the characteristic time is larger. Pre-stiffening waves will also be detected with a time-step ∆t = 1 × 10 -4 s. The mitral valve is mostly composed of water and thus is nearly incompressible. Artificial compressibility is used to model the mitral valve with 3D elements. The development above only study the elastic response of the material. In case the material is compressible, pressure waves inside the material also exist. Their characteristic speed is defined as:

c κ = κ ρ
with κ the bulk modulus. The pressure waves inside the leaflets are less important for the solid response than the elastic waves but still carry information. The time-step should be smaller than the time required for the pressure wave to go through the domain. Thus ∆t ≤ ∆t κ , ie:

κ ≤ L 2 s ρ ∆t 2
Using the time-step previously obtained and the data table 3.1, κ should be lower than 4 × 10 7 Pa to detect the pressure waves during the simulation. For comparison, the bulk modulus of water is about 2.2 × 10 9 Pa, which correspond to a time-step of ∆t = 1 × 10 -5 s.

Time scheme

Several time-stepping methods are available to represent the temporal evolution of the valve. The explicit method's time-step is limited to the velocity of information through the smallest element in the mesh, but the computational effort for each time-step is relatively small. Implicit methods are unconditionally stable and thus can use larger time-steps. However for each time-step, an inverse problem Ax = b must be solved. Explicit methods are usually preferred in the non-linear dynamic analysis because the small time-step reduces the problem's non-linearity per step. However, implicit methods can also solve non-linear problems efficiently.

In this study, we choose to use an implicit method. Indeed, the software ADINA does not allow the use of user-coded material with an explicit time-stepping. Nevertheless, the non-linearity of the structural problem mainly occurs during the closure of the mitral valve. Once the mitral valve is closed, and the contact is established, the problem is "less" non-linear. In particular, the mitral valve movement is much less pronounced, and a larger time-step could then be used. Also, we aim to solve the interaction between the blood and the valve using a monolithic approach. Thus the same temporal discretization will be used for the fluid and the structure. Implicit methods are generally preferred when modelling incompressible fluids.

Using an implicit formulation that is unconditionally stable allow choosing the time-step regardless of the velocity of pressure waves. As a result, it is not mandatory to lower the bulk modulus of the leaflets or to increase the time-step to detect this phenomenon.

The time scheme used is a second-order implicit time scheme called the Bathe method. It was initially developed for non-linear analysis with large deformation. It is unconditionally stable and naturally provides numerical damping. For each time-step ∆t, the method used two sub-steps. In the first sub-step, the Newmark trapezoidal rule is used. In the second sub-step, the 3-point Euler 59 backward method is employed [101]. The formulation employed is given as:

t+∆t/2 U = t U + ∆t 4 t Ü + t+∆t/2 Ü (3.1) t+∆t/2 U = t U + ∆t 4 t U + t+∆t/2 U (3.2) t+∆t U = 1 ∆t t U - 4 ∆t t+∆t/2 U + 3 ∆t t+∆t U (3.3) t+∆t Ü = 1 ∆t t U - 4 ∆t t+∆t/2 U + 3 ∆t t+∆t U (3.4)
with U nodal displacement, rotation and pressure. The structural dynamics equations at time t + ∆t and t + ∆t/2:

M t+∆t/2 Ü + C t+∆t/2 U + t+∆t/2 F = t+∆t/2 R (3.5) M t+∆t Ü + C t+∆t U + t+∆t F = t+∆t R (3.6)
with F = KU the vector of nodal point forces, M , C and K the mass, damping and tangent stiffness matrix, R the vector of external applied loads and overdot the time derivative. Equations 3.1 to 3.6 can be written in the form:

K1 t+∆t/2 U = R1 (3.7) K2 t+∆t U = R2 (3.8)
with the equivalent tangent stiffness matrix:

K1 = t+∆t/2 K + 16 ∆t 2 M + 4 ∆t C K2 = t+∆t K + 9 ∆t 2 M + 3 ∆t C
and the equivalent external applied loads:

R1 = t+∆t/2 R + M 16 ∆t 2 t U + 8 ∆t t U + t Ü + C 4 ∆t t U + t U R2 = t+∆t R + M 12 ∆t 2 t+∆t/2 U + 4 ∆t t+∆t/2 U - 1 ∆t t U - 3 ∆t 2 t U + C 4 ∆t t+∆t/2 U - 1 ∆t t U
The problem being non-linear, the matrices K1,2 are updated for each sub-time-step, and each iterations of the solver. The resolution of the equations 3.7 and 3.8 are detailed in subsection 3.2.2.

Time-step refinement

When the contact between the leaflets occurs, the nodes of the leaflets' surfaces must be paired to transmit a reaction force. The nodes are paired depending on their relative distances, which vary during the solver iterations for one step. The resulting non-linearity can be too important to solve within the time-step initially provided, especially when the valve largely deformed to obtain its closed configuration. Using smaller time-steps reduces the leaflets' displacements and the non-linearity resulting from the node's connection. ADINA provides a method called ATS (Automatic-Time-Stepping) that can refine the time-step during contact detection. The method uses four parameters: a division factor DIVFACT, a maximal number of sub-increments MAXSUB, a multiplication factor MULFACT and a number of convergent time-steps NCONV. When the solver can not solve the problem for a given time-step ∆t ini , the ATS method divides the time-step with DIVFACT and try to solve the problem again with the new time-step. If the resolution fails, the time-step is further subdivided until convergence is reach or the time-step is below ∆t ini /MAXSUB. If the problem is successfully solved NCONV consecutive times with a time-step, then the time-step is multiplied by MULFACT. The method was configured to use the previous time-step if this criterion was not respected. To summarise:

∆t ini MAXSUB ≤ ∆t ≤

Numerical resolution

The material behaviour and the formulation used are non-linear. Thus,the evolution of nodal forces t F is not linear with the displacement t U , because the tangent stiffness matrix K is also function of the displacement. As a consequence, it is not possible to solve equations 3.7 and 3.8 directly because the tangent stiffness matrices K1 and K2 depends on values defined at t + ∆t/2 and t + ∆t:

t+∆t/2 K1 t+∆t/2 U = R1 (3.9) t+∆t K2 t+∆t U = R2 (3.10)
To be solved, the equations must first be linearized. Newton-Raphson algorithm was used to solved the problem ∀t t f = 0, with ∀t t f = t F -t R.

In this paragraph, only the first sub-time-step is detailled (equation 3.9). The equations can easily be transferred to equation 3.10, as they only differ from the coefficients used to define K1,2 and R1,2 . Considering only the first substep, function f reads:

f (U ) = 0 (3.11) f : U → t+∆t/2 R1 - t+∆t/2 F1 (3.12)
with t+∆t/2 R1 and t+∆t/2 F1 = t+∆t/2 K1 t+∆t/2 U the equivalent external loads and the equivalent nodal point forces defined in section 3.2.1.

Using first order Taylor's development on equation 3.11, we have:

f (U ) = f t+∆t/2 U (i-1) + ∂f ∂ t+∆t/2 U t+∆t/2 U (i-1) ∆U i + O ∆U i 2 (3.13) f (U ) = f t+∆t/2 U (i-1) - ∂ t+∆t/2 F1 ∂ t+∆t/2 U t+∆t/2 U (i-1) ∆U i + O ∆U i 2 (3.14)
with ∆U i = U -t+∆t/2 U (i-1) , and t+∆t/2 R1 being independent of t+∆t/2 U . This assumption supposed that the load t+∆t/2 R in equation 3.5 is not function of the current displacement. If its not the case, and additional term is added in equation 3.14, and can modify the convergence of the algorithm. This can be the case when contact occurs, and is discussed in section 3.2.3. Equation 3.14 use the derivative of the equivalent nodal point forces towards the current displacement. The equivalent nodal point forces represent the forces due to inertia, damping, and elastic response. The mass matrix M is independent of the deformation, thus has no derivative toward the current displacement. Furthermore, no damping is considered in the structural response, and the viscous response of the material was not modelled. As a result, we have:

∂ t+∆t/2 K1 ∂ t+∆t/2 U = ∂ t+∆t/2 K ∂ t+∆t/2 U (3.15)
Using the definition of the equivalent nodal point forces, this leads to:

∂ t+∆t/2 F1 ∂ t+∆t/2 U t+∆t/2 U (i-1) = t+∆t/2 K1 (i-1) + ∂ t+∆t/2 K ∂ t+∆t/2 U (i-1) t+∆t/2 U (i-1) t+∆t/2 U (i-1) (3.16)
For simplicity we note t+∆t/2 K (i-1) the derivative of t+∆t/2 K at t+∆t/2 U (i-1) . Several algorithms are available in the software ADINA to find the solution of the equations 3.9 and 3.10. The equation 3.14 can now be rewritten using equations 3.11, 3.12 and 3.16:

t+∆t/2 K1 (i-1) + t+∆t/2 K (i-1) t+∆t/2 U (i-1) ∆U i = t+∆t/2 R1 (i-1) - t+∆t/2 F1 (i-1) + O ∆U i 2
Considering the iterative process and neglecting higher-order terms, we obtain a linear expression on ∆U i : t+∆t/2 K(i-1)

1 ∆U i = t+∆t/2 R1 (i-1) - t+∆t/2 F1 (i-1) (3.17)
with ∆U i = t+∆t/2 U i -t+∆t/2 U (i-1) and t+∆t/2 K(i-1)

1 = t+∆t/2 K1 (i-1) + t+∆t/2 K (i-1) t+∆t/2 U (i-1) .
The effective solving of the linearized equation 3.17 is dependent on the properties of K(i-1)

1
, which is a combination of the K1 and K matrices. Thus, the mass, damping, elastic and tangent stiffness matrix must be provided specific attention during their constitution.

The mass matrix is banded diagonal and is not modified during our simulation. There is no damping in our first model; thus, the damping matrix is set to 0. As a result, the accurate definition of the matrix K(i-1)

1 is a function of the accuracy of the previous unknown, the tangent stiffness and the constitutive matrix at each iteration. These matrices are computed precisely from the strain energy density function of the material in the user-coded subroutine using a Total Lagrangian formulation. This accurate definition of the matrix K(i-1)

1 leads to a convergence rate of order 2. The solution was considered converged according to relatives residuals. At time t + ∆t/2, the relatives energetic residual is defined as:

E tol = ∆U i T t+∆t/2 K(i-1) 1 ∆U i ∆U 1 T t+∆t/2 K0 1 ∆U 1
And the force residual is defined as:

F tol = t+∆t/2 K(i-1) 1 ∆U i 2 RNORM
The default use of energy relative residual led to poor precision and addition of errors, even when lowered to 10 -6 . We hypothesize that the energy residual uses the displacement computed for the previous time step (∆U 1 ), which can greatly differ significantly from the desired displacement ∆U i . On the other hand, the force residual is less influenced by the initial displacement (RNORM is either fixed or computed from the program) and monitors directly the difference obtained from equation 3.17. It is considered more numerically demanding to solve but leads to more accurate results when compared to convergence with energy residual. The relative force residual is often used in non-linear analysis. Thus, a relative force residual of 10 -3 was used to ensure convergence for each time-step.

A direct solver based on LU Decomposition solves the equations 3.17. The solver is named SPARSE in ADINA and uses sparse storage of the matrices to limit storage usage and CPU time. This direct solver is more efficient than iterative solvers for problems smaller than 500 000 degrees of freedom. For larger problems, iterative solvers scaled linearly with the number of unknowns and can be used to improve CPU time.

Contact method

A contact was implemented to detect the interaction of the posterior leaflet on the anterior leaflet. A contact is an interaction between nodes and an element surface depicted as a target. The gap g is the normal distance between the contactor node and the target surface: g = -(x c -x t ).n with x c,t the position of the contactor or target and n the normal of the target. The overlap h is thus h = -g. The ideal contact configuration is that the normal force of contact λ g is transmitted only when g = 0, and no overlap can be observed (figure 3.4a). This ideal case can be modelled in ADINA using a tolerance N with λ g = N /g (figure 3.4b). N is set to 10 -12 and is nearly independent of the problem studied. A force is always transmitted between the target and the contactor in this configuration, which is not ideal. However, the normal force is difficult to capture in this situation when the gap is close to zero.

A compliant factor is often used to ease the numerical problem and allow convergence. With compliance, the contactor node can go through the segment target, and the normal pressure is a function of the overlap and the compliant factor ĉ (figure 3.4c). The value for ĉ is problem dependant. It must be small enough so that the contact appears as one patch instead of an addition of several areas. If it is too large, the problem is more difficult to solve. The value ĉ = 10 -8 m 3 .N -1 was thus set experimentally. Contact normal response is enforced at each contact node and element using a constraint function based on Lagrange multipliers.

The use of compliance allows the contactor node to go through the target. In order to avoid the geometrical interpenetration of the leaflets, an offset can be applied to the structure. The use of the offset project the surface of contact using the normal of the targets (figure 3.5). Offsets of h off = 0.8 mm was used in our simulations to avoid geometrical interpenetration. in equation 3.17. This force depends on the contactor nodes and surface target and their relative distance. It adds non-linearities to the problem and can lead to oscillations during the Newton-Raphson iterations. To prevent such oscillations, each contact node has its coupling to a surface target fixed after some iterations. The limit number of iterations for coupling the contactor nodes was set to 5.

Elements used

The geometrical mitral valve model obtained in chapter 1 was discretized using shells and 3D elements. Both elements used a hybrid U/P formulation. This formulation allows computing pressure from a distinct degree of freedom instead of the displacement quantities. This feature prevents "locking" behaviour. Thus, fewer elements are required to model the pressure field accurately in an area where compression can be observed.

3D-shells elements

The 3D-shells used are 4-nodes MITC shells (first-order accurate). 3D-shells elements are shells where the mid-fibre position is recomputed as a function of their current deformation. The mid-surface of the element is defined using control vectors (figure 3.6). As the change in thickness is defined with the control vector, the assumption that no stress exists through the shell thickness is unnecessary. The 3D-shells elements also use MITC tying rules to prevent shear locking. As a result, the 3D-shells can be used in large deformation and large strain analyses. The mitral valve was discretized using 3D-shells of size ∆x = 0.1 cm and thickness ∆e = 0.08 cm.

3D elements

The 3D elements used are 11 nodes tetrahedron elements (second-order accurate). The hybrid formulation is used with four pressure degrees of freedom to capt linear pressure variation inside the element. These are the default choices for nearly incompressible material and can be used in large deformation and large strain.

The 3D elements are a natural choice to detect internal stress variations in the mitral leaflets. Medical devices can be apposed on the leaflets and slightly grip them. 3D elements are then natural choices to study the influence of the grasping system.

When considering fluid-structure interaction, the envelope of the structural body is used as a boundary for the fluid. Shell model of the mitral leaflets leads to one surface boundary used for both the atrium and the ventricle leaflet's boundary. From the fluid point of view, the two boundaries can not overlap. Nevertheless, some elements could overlap due to the large deformation and displacement of the structure and the remeshing technique used. Using 3D elements ensure a minimal distance between the two fluid domain is always respected.

Finally, the 3D elements can be used with user-coded material laws in ADINA 9.5, while the shell model is restricted to some material laws already implemented in the software.

In the simulation, a spatial discretization of ∆x=0.1 cm is used unless stated otherwise.

Truss elements

Chordae tendineae are modelled with 2-nodes truss elements. The structural dynamic response of the chords was neglected during the simulation. Only one element is used for each chord to prevent undesired deformations and numerical oscillations. They only contribute to the model during traction. In the simulation, a cross-section of 0.6 mm 2 is used.

Damping methods

When using the Sussman-Bathe material law with 3D elements, numerical oscillations and vibrations appeared early in the simulation, leading to non-convergence. A damping technique called Rayleigh damping was used to ease the resolution of the problem. The Rayleigh damping is an additional damping matrix C R added to the equation 3.5 and 3.6 defined as:

C R = αM + βK (3.18)
where M and K are the mass and the tangent stiffness matrix of the numerical system. α (respectively β) dampen the low (respectively high) frequencies. The corresponding damping function for all frequencies is given in figure 3.7. When α is set to 0, β can be set to dampen the variations of characteristic time higher than T p with T p = πβ. Thus β was initially set to 3.18 × 10 -5 s to prevent noise that should not be detected during the simulation. This value for the Rayleigh parameter was insufficient to reach convergence, and the material instabilities still led to unsuccessful analysis. β parameter was thus increased, along with additional viscous damping on all contact nodes. This additional viscous damper adds a force to the contact nodes proportional to their absolute velocity and a parameter C damp . It was necessary to have β ≥ 1 × 10 -3 s to obtain a converged simulation. When computed successfully, the simulation used time-step from 1 × 10 -4 s to 1 × 10 -6 s because of the automated time-stepping method used. This time-stepping method is necessary to detect and compute the quick contact between the leaflets. The resulting computational time was relatively high (approximately four hours) and was not much influenced by the additional viscous damping (figure 3.8). However, increasing the value of β seems to lower the computational cost. The drawback of this strong added damping is that the dynamical response can only be used as a static result once the valve movement has stopped and closure is reached.

It should be noted that Rayleigh damping should be updated in large displacement simulation. Indeed, the tangent stiffness matrix for our problem is non-linear and depends on the deformation. Thus equation 3.18 should be updated for each time-step, and each increment of the Newton-Raphson t+∆t C (i) r . The update was not possible in our simulations. When updating C r the resolution is even worse and fails to converge even for small damping values. One hypothesis is that the equivalent tangent stiffness matrix used for the Newton-Raphson iteration K(i) 1 in equation 3.17 is penalized by the variation of the damping matrix C r . Oppositely, it accepts its influence when C r is constant.

Results of mitral valve simulations

Several mitral valve simulations were performed to observe the geometrical closure of the valve and the internal stress during the cardiac cycle. Dynamic and static responses were obtained with a shell model and physiological pressure evolution. 3D elements mitral valve model represent static closure at systolic pressure.

Results of the shell model

Dynamic mitral valve closure

The valve geometry used for this model is obtained from the first Fortran & FreeCAD model defined in chapter 1. The annulus is a planar D-shape with an inter-commissural length of 30 mm and anteroposterior length of 22.5 mm. The anterior leaflet height is set to 20.2 mm. The height for the posterior leaflets P1, P2 and P3 are respectively 9.9, 12.0 and 9.9 mm. Commissural height is set to 6.1 mm. The papillary muscles are positioned according to the anatomical measurements (data in figure 1.5). Eighteen chords are used to model the chordae tendineae with five subbranches each. A front view of the mitral valve is depicted in figure 3.9. The mitral valve geometry obtained is also seen figure 1.23, chapter 1.

The annulus and the papillary muscles points have fixed displacements during the whole simulation. Uniform systolic pressure is applied according to the pressure evolution (figure 3.9c). The peak pressure applied is 120 mmHg 16 KPa.

The material used for the leaflets in this simulation is a Sussman-Bathe material previously defined in section 2.13, chapter 2. The material was tested under equibiaxial numerical test and fitted experimental equibiaxial data in the circumferential direction [START_REF] Newman | A constitutive law for mitral valve tissue[END_REF]. The chordae tendineae have a non-linear stress-strain behaviour and follow the physiological response to strain (figure 2.3, [START_REF] Kunzelman | Mechanical properties of basal and marginal mitral valve chordae tendineae[END_REF]).

The valve was discretized using 3D-shells in a hybrid U/P formulation of size ∆x=0.1 cm. Their thickness is set to ∆e = 0.08 cm. The chords are modelled with traction only truss. One truss was used for each chord. Attempts to use smaller meshes resulted in stronger numerical instabilities and failed to reach convergence. This instability could be linked to the material definition and possible incorrect definition, as detailed in chapter 2. The implicit second-order accurate Bathe time scheme was used with a sub-time-step at ∆t/2. The simulation was run for several cycles of period 1 s to observe established mechanical response to pressure variation. The time-step was set to ∆t = 10 -4 s and was divided if convergence was not reached down to 10 -6 s. Unfortunately, due to the material non-linearities and instabilities depicted in chapter 2, the time-step used in practice during the simulation oscillates between 10 -5 s and 10 -6 s. A compliant contact was used with a compliance of 10 -8 m 3 .N -1 to enable the leaflet coaptation. A geometrical offset equal to the shell thickness was added to avoid geometrical interpenetration.

The stress response in the mitral valve closure is given figure 3.11 and 3.12. The maximum stress was obtained at t = 0.5 s, which correspond to the systolic peak pressure. Its value is 1.7 MPa. The maximum stress was located at the connection between the leaflets and the chords. On the inner surfaces of the anterior and posterior leaflets, the stress is relatively small on the mid-surfaces. However, the stress located at the top and bottom shell thickness are more important and similar to the values of 0.5 MPa obtained in the literature (figure 3.10).

The coaptation surface between the leaflets appears at t = 0.283 s and disappears at t = 0.700 s. This corresponds to the systolic phase defined with the pressure variations. The maximal value of the coaptation surface is 2.18 cm 2 . The large coaptation area obtained ensures that the modelled mitral valve close well during systole.

Mitral valve closed configuration

The previous geometry was used with several material laws to test their influence on the close state of the valve. Except for the pressure evolution and the Rayleigh damping, the same discretizations and methods were used. The pressure used in this model was a linear ramp from 0 to the systolic pressure 120 mmHg 16 KPa in 0.2 seconds. The Raleigh damping factors were α = 0 s -1 and β = 10 -3 s to allow convergence of the models with reasonable time-steps. The material laws used were the isotropic law, the Mooney-Rivlin law and the Sussman-Bathe law depicted in chapter 2. The isotropic parameter were E = 3M P a and ν = 0.45. The Mooney-Rivlin coefficients were obtained from the literature and correspond to a fit of experimental data [START_REF] Abedini | Numerical simulation of the impact of chorda tendinea rapture on the mitral valve leaflets stress by finite element method[END_REF]. The Sussman-Bathe coefficients are the same as the dynamic study and correspond to a fit of experimental data in the circumferential direction. Results of the simulation are presented figure 3.13 and where detailed in section 2.2.3. The closure obtained with the isotropic is only partial. A coherent coaptation area is reached using hyperelastic material such as the Mooney-Rivlin or the Sussman-Bathe material. The anterior leaflet of the Mooney-Rivlin model has some kind of protrusion above the annulus plane that is not realistic. Among the isotropic material tested, only the Sussman-Bathe resulted in realistic closure of the mitral valve.

Front view

Results of the 3D-elements model

The shell models of the mitral valve are limited to the material already coded in the software AD-INA. Among the material tested, the Sussman-Bathe law leads to a better closure shape during systole.

However, its simulation required small time-steps and increased the computational cost of the problem. Larger time-steps could be used if the simulation was dampened, but the valve's dynamic response was ignored. Also, it was not possible to use finer meshes. This motivates the development of the 3D-element model, compatible with user-coded material response.

3D-elements model with damping

A first model using 3D elements was developed based on the same mitral geometry as the shell model. Most parameters are the same as the dynamic shell model. A uniform thickness of ∆e = 0.05 cm was applied to thicken the surfaces. The annulus and papillary muscles have fixed displacements. Uniform pressure is applied on the surfaces of the leaflets. The pressure evolution is depicted figure 3.14. The pressure ramp is based on the systolic pressure ramp. A physiological systolic pressure of 120 mmHg is reached at 0.2 s and is then kept constant. The mesh is composed of 11-nodes tetrahedrons using a U/P formulation. 10 000 tetrahedron of referent size ∆x = 0.05 cm were used to model the leaflets. Each chord is modelled with traction only truss of one element. The material used for the leaflets in this simulation is a Sussman-Bathe material previously defined in section 2.13, chapter 2. The material was tested under equibiaxial numerical test and fitted experimental equibiaxial data in the circumferential direction [START_REF] Newman | A constitutive law for mitral valve tissue[END_REF]. The chordae tendineae have a non-linear stress-strain behaviour and follow the physiological response to strain (figure 2.3, [START_REF] Kunzelman | Mechanical properties of basal and marginal mitral valve chordae tendineae[END_REF]).

The implicit second-order accurate Bathe time scheme was used with a sub-time-step at ∆t/2. The simulation was run for several cycles of period 1 s to observe established mechanical response to pressure variation. The time-step was set to ∆t = 10 -4 s and was divided if convergence was not reached down to 10 -6 s.

A compliant contact was used with a compliance of 10 -8 m 3 .N -1 to enable the leaflet coaptation. A geometrical offset of 0.1 mm was added to avoid geometrical interpenetration.

A stiffness proportional Rayleigh damping parameter was used to soften the dynamic response (β = 10 -3 s). The simulation was considered converged when the velocity magnitude of the valve reached 0. Indeed, the damping effect should not affect the static solution obtained.

Results are depicted figure 3.15 and 3.16. The stress obtained reach maximal values on the anterior leaflet and at the insertion sites of the chords. The maximal stress obtained is around 0.6 MPa, in coherence with previous simulation and the literature. The nodal distance shown in figure 3.16 represents a smooth area at the bottom of the anterior leaflet and at the middle of the posterior leaflet. The unequal repartition of the coaptation area questions the resting length of the chords. Looser posterior chords might improve the coaptation area. Overall, the coaptation area consists of one large smooth patch, thus confirming the sealing aspect of the mitral valve when closed.

3D-elements with user-coded material

The user-coded material was used to define the law of Gao et al., the law of Prot et al. and the law of May-Newmann and Yin. Convergence difficulties appeared early in the simulation when the anisotropic material were applied to the previous model. The simulations aborted at the first or second time-steps. Oddly, better results were obtained when no Rayleigh damping was used. We suppose that Rayleigh damping method might perturbs the equivalent tangent stiffness matrix used during the Newton-Raphson iterations. Indeed, the equivalent tangent stiffness matrix is re-evaluated for each iterations, but the Rayleigh-damping matrix C R is only updated at each time-step. However, the matrix C R is defined as a function of the tangent stiffness matrix K, which evolves during the Newton-Raphson iterations and should be updated.

Only the use of the material law from Gao and co-workers leads to the computation of several time-steps and shows relevant information on the valve deformation (figure 3.17 and 3.18). The other material law abort at the first time-steps. It is expected that this material law is the easiest to compute (when compared to the other anisotropic models) because it as the strongest Neo-Hookean term and the smallest anisotropic response to deformation among the anisotropic material laws studied.

The results obtained figures 3.17 and 3.18 are difficult to compare to the literature as the closure state is not obtained. In diastole configuration, the valve is largely opened, and the stress and stretch maximal area are located at the chords insertion. During early systole, the valve closes quickly. However, the maximal stress obtained (0.05 Pa) is not representative of the closed mitral valve maximum stress. Nonetheless, it is small enough to suppose that if the valve would closed, it could be similar to the literature value of around 0.5 MPa.

Conclusion and discussion

Structural simulations of the mitral valve were run successfully using 3D elements and shells model. In both cases, the closed state of the mitral valve was realistic and in agreement with the literature. The complex dynamic response of the valve under physiological pressure was computed successfully with the shell model. The mechanical response of the mitral valve leaflet is hyperelastic and highly non-linear. A standard hyperelastic model like the Mooney-Rivlin model was insufficient to model the closed configuration of the mitral valve. Other hyperelastic material laws such as the Sussman-Bathe material law can accurately model the closed mitral valve. However, this model was related to numerical instabilities and made us use a smaller time-step than required. Non-physical damping was added to the problem using Rayleigh damping. Closure simulations of the mitral valve were then obtained with lower CPU time, but the dynamic response of the valve was ignored.

Other materials that verify a physical stability criterion were developed in chapter 2. These materials do not need the additional Rayleigh damping but have a high non-linear response. The most challenging ones are the model of May-Newman and Yin and the model of Prot and co-workers. The model of Gao et al. is less non-linear, as the fibre stress response to strain is less pronounced. As a result, this model was the only one that could run partially in our models.

The time-step could be further subdivided to improve the convergence rate of the user-coded material. However, attempts to run such simulations with a smaller time-step failed. In addition, using smaller time-steps drastically increase the CPU time required. An alternative would be to use an explicit time scheme, but it is not possible in ADINA using user-coded material.

Another solution would be to use ULH formulation rather than TL formulation to describe the material behaviour. This formulation is successfully used in the literature [START_REF] Sacks | On the Simulation of Mitral Valve Function in Health, Disease, and Treatment[END_REF].

Chapter 4

Fluid-Structure Interaction of the mitral valve

This chapter describes the interaction between the mitral valve and the blood. The first section depicts the advantages and disadvantages of modelling the mitral valve structure and the surrounding fluids. The second section explains the methodology used to model the mitral valve. The third section details the results of the simulations. Finally, the last section contains discussions on the results and possible improvements for the model.

Introduction

This section has two aims. First, to explain why fluid-structure interactions of the mitral valve are interesting. Second, to detail how fluid-structure interactions are run effectively. The motivation toward fluid-structure interaction is first detailed. Next, the critical difficulties to overcome for these simulations are presented. The main types of fluid-structure interactions are then depicted and illustrated, and we conclude by choosing the software adapted to our needs.

Motivations for Fluid-Structure simulations of the mitral valve

The mitral valve acts as a hydraulic diode. When the valve opens, the pressure pushes blood from the atrium toward the ventricle. The valve seals the ventricle when closed, and blood can not cross the annulus to go back into the atrium. Some diseases alter the valve function, and a backward flow can occur when the mitral valve is closed. The volume of reverse flow, called the volume of regurgitation, quantifies the pathology's severity. The mitral regurgitation jet can be measured precisely with echocardiographic imaging tools. The volume of regurgitation is computed from this jet (figure 4.1). The jet speed can reach several meters per second. Regurgitant volume is considered pathological when superior to 45mL [START_REF] Grayburn | Quantitation of mitral regurgitation[END_REF].

Radiologists can use software like Arterys or CMR42 to reconstruct the blood flow. They use patient-specific data from imaging techniques to represent the flow of blood. The software uses invivo data and CFD post-processing technique to plot the flow of blood in a realistic beating heart. In the case of mitral insufficiency, the flow passing through the annulus can be reconstructed during a cardiac pulsation. The maximal jet speed and the volume of regurgitant blood can be computed accurately. The software estimations help to characterise the pathology. The medical details obtained help surgeons and radiologists in making decisions. However, they are based on existing data and can not be used to test the outcome of medical procedures. Structural simulation computes well the mitral valve closed state using pressure evolution and heart contraction and does not require imaging techniques. The fluid simulation uses the Navier-Stokes equation and boundary conditions to compute blood flow. When imaging techniques are unavailable, it is the most accurate method to obtain the regurgitant volume or maximal jet speed. During a cardiac cycle, the mitral leaflets modify the flow of blood. Thus, it is impossible to neglect the structure while computing the fluid equations. The blood and the valve must be modelled simultaneously to obtain realistic fluid flow. This model is called a fluid-structure interaction (FSI). FSI simulation also gives a more accurate pressure distribution in time and space on the leaflets, improving the valve motion and internal stress state. As a result, the dynamic response of the mitral valve is more accurate during fluid-structure simulation than structural simulation. Fluid-structure interaction also captures the vortex resulting from the passage of blood through the valve. The vortex influences the leaflet and helps in closing the valve.

In fluid-structure interactions, blood flow is reconstructed as a function of the heart dynamic evolution and boundary condition. Thus, this type of simulation can represent the effect of a medical procedure on blood flow. For example, a pathological mitral valve can be modelled with a medical procedure such as the MitralClip. The fluid equations reconstruct the blood flow, regurgitant leak and maximal jet before and after the procedure to estimate the procedure efficiency. It is also possible to test several procedures on the same pathological valve and compare the results to determine which one is the best to repair the valve function. With the same method, a medical device can be optimised by modifying its components, shape or implementation.

Fluid-Structure Interaction difficulties

Fluid-structure interaction is complex to model. FSI couples the fluid and the solid mechanical responses and inherits the difficulties from both models. In addition, the fluid and the solid problems must exchange information on their interface, its kinematic and the stress transmitted through it. In the case of the mitral valve, there are additional fluid difficulties. First, the flow modelled during a cardiac cycle can go backward. Second, the fluid contained in the ventricle can be separated from the fluid contained in the atrium during the mitral closure.

This section first discusses the types of FSI. The difficulties related to the mitral valve specific interaction are then detailed.

Types of Fluid-Structure Interaction

Temporal Coupling

Let us consider a domain Ω where the fluid Ω f and the solid Ω s have a common continuous interface Γ. The fluid behaviour is depicted by the Navier-Stokes equation inside Ω f (see annexe eq. 6.13). The equations of solid mechanics are used inside Ω s (see annex eq. 6.12). At their interface Γ, the fluid velocity is equal to the velocity of the solid due to the no-split condition. This condition is referred to as the kinematic coupling condition. In addition, the fluid and solid have a balance of normal stresses at the boundary due to Newton's third of dynamic. This condition is referred to as the dynamic coupling condition. Both fluid and solid depend on the interface and the kinematic and dynamic conditions. The two physics are coupled.

Several approaches can be used to model multi-physic problems. In the partitioned approach, two solvers are used. One is dedicated to fluid (F ). The second (S) compute the structural response. Considering discrete time t n , the variable f n , s n and Ω n represent respectively time fluid response, the solid response and domains and interface at time t n . The simplest coupling method consist in computing the solid response at time t n+1 using previous fluid and interface information. The new interface Ω n+1 can be deduced from s n+1 , and the fluid response computed:

     S(f n , s n , Ω n ) → s n+1 Ω n+1 = func(s n+1 , Ω n ) F (f n , s n+1 , Ω n+1 ) → f n+1
This partitioned coupling using one (or few) iteration(s) between the solver is called weak or loose coupling (figure 4.2). It is fast to obtain coupling results but can lead to instabilities and fail to converge [START_REF] Causin | Added-mass effect in the design of partitioned algorithms for fluid-structure problems[END_REF].

Iterations of the fluid and solid responses can be exchanged several times for a given time step. The successive results ensure that the fluid and solid response are coherent and respect some criteria monitoring the coupling errors: with

∀i ∈ N        S(f (i) n+1 , s (i) n+1 , Ω (i) n+1 ) → s (i+1) n+1 Ω (i+1) n+1 = func(s (i+1) n+1 , Ω n ) F (f (i) n+1 , s (i+1) n+1 , Ω (i+1) n+1 ) → f (i+1) n+1
Ω (0) n+1 = Ω n , s (0) n+1 = s n and f (0) n+1 = f n .
This type of coupling is called strong coupling (figure 4.3). It is the most accurate partitioned coupling method but can be costly as the structure and fluid solvers are computed several times for one timestep. The coupling is unconditionally stable but the number of coupling iterations can be large, in particular when difficulties are encountered and weak coupling fail to converge. In particular, fluid-structure problems where the added-mass is important can be challenging to solve [START_REF] Causin | Added-mass effect in the design of partitioned algorithms for fluid-structure problems[END_REF], [START_REF] Fernández | A projection semi-implicit scheme for the coupling of an elastic structure with an incompressible fluid[END_REF]. The other approach is to directly solve the fluid and solid problem altogether (figure 4.4). This approach is called monolithic coupling or direct approach and does not stagger the fluid response to the structure (and vice-versa). In this method, the fluid and the structural problem are detailed in the same numerical problem. One multi-physic matrix is obtained and solved in one solver M SF , for the time t n . No communication between the solvers is required :

M SF (f n , s n , Ω n ) → s n+1 , f n+1 , Ω n+1

Spatial Coupling

The fluid and the solid interface can be modelled with different techniques. Body-fitted methods, fictitious-domain methods and smooth particle hydrodynamic techniques are discussed in this section. The difficulties inherent to the mitral valve simulation is detailed, with pros and cons for each method. Body-fitted methods define a mesh that follows the interface between the fluid and the solid 4.5b. The Arbitrary Lagrangian-Eulerian (ALE) is a finite element formulation in which the mesh is neither fixed in space (Eulerian mesh) nor fixed to the material (Lagrangian mesh). ALE mesh motion is set to Lagrangian mesh at the boundary to follow the interface accurately, while the core of the mesh can adapt to define the best shape for each element due to the Eulerian formulation. This bodyfitted method can follow a moving interface but has several drawbacks. First, the initial generation of the mesh can be difficult to perform depending on the geometry. Second, if the displacement of the interface is large, the fluid mesh becomes degenerated, and remeshing should occur before continuing the displacement. During remeshing, a new mesh is generated using the deformed geometry, and data from the older mesh is interpolated on the new one, leading to error and diffusivity. Due to meshes generations and interpolations, ALE methods can be computationally expensive. However, they accurately depict the fluid-solid interface and its values, such as the pressure or the wall shear stress. In addition, the formulation is adapted to multi-physic implementation. This method is favoured when remeshing can be avoided. It was used with the software YALES2BIO to model a left ventricle contraction [START_REF] Chnafa | Image-Based Large-Eddy Simulation in a Realistic Left Heart[END_REF]. The ALE formulation was used with Ansys Fluent on a 3D mechanical valve [START_REF] Bavo | Fluid-Structure Interaction Simulation of Prosthetic Aortic Valves: Comparison between Immersed Boundary and Arbitrary Lagrangian-Eulerian Techniques for the Mesh Representation[END_REF]. However, the computation aborted, probably due to remeshing difficulty and the coupling technique used. In [START_REF] Forsythe | Validation of a Fluid-Structure Interaction Model for a Bileaflet Mechanical Heart Valve[END_REF], a free 3D leaflet was modelled using the ALE technique, and a local remeshing technique was used to limit the errors linked to the interpolation.

Fictitious-domain methods use a generic fixed grid for the whole domain, and several approaches are available to detect and enforce the interface response to fluid. The method has a spread interface if the dimension of the interface is the same as the domain. The Immersed Boundary Method (IBM) is the most common among this class of methods. With the IBM, the eulerian variables (pressure p, velocity v...) are computed on the fixed mesh, and the interface is discretised to compute its Lagrangian variables (velocity v Γ , position x Γ , enforcing term f Γ ...). The cells connected to the interface are penalised for enforcing the boundary (green cells in the figure 4.5c). Several methods exist to penalise the cells of the interface. Difficulties with the IBM appear when the interface is moving because the Lagrangian variables must be recomputed at each timestep. Depending on the method used to enforce the interface, small pressure or fluid leak can appear through the interface. Interpolations at the solid interface can lead to less precise wall shear stress and pressure than ALE formulations. However, the mesh and its quality were kept constant throughout the simulation. This appealing method was initially developed for cardiovascular simulation purposes and is promising [2], [START_REF] Peskin | The Immersed Boundary Method[END_REF]. IBM was applied to model prosthetic mitral valve [START_REF] Watton | Dynamic modelling of prosthetic chorded mitral valves using the immersed boundary method[END_REF], aortic valve [START_REF] Bavo | Fluid-Structure Interaction Simulation of Prosthetic Aortic Valves: Comparison between Immersed Boundary and Arbitrary Lagrangian-Eulerian Techniques for the Mesh Representation[END_REF] and [START_REF] Borazjani | Fluid-Structure Interaction, Immersed Boundary-Finite Element Method Simulations of Bio-Prosthetic Heart Valves[END_REF] (with the CURVIB method), and more recently mitral valve with the IBAMR code [START_REF] Gao | A Coupled Mitral Valve -Left Ventricle Model with Fluid-Structure Interaction[END_REF]. The code YALES2BIO was also used with an immersed boundary to model the valve motion [START_REF] Chnafa | Image-Based Large-Eddy Simulation in a Realistic Left Heart[END_REF].

Smooth particle hydrodynamic is a technique where the flow is not directly modelled. Instead, particles are carried with the flow. Hydrodynamic properties are then evaluated at the particles' positions using weighted averages. Thus, the values defined at the particle position are smoothed at the region around the particle position (green area in the figure 4.5d). The advantage of SPH methods is that the relative position of the particles is not modifying the quality of the computation compared to the nodes of a standard mesh. Another advantage is that the particle can be inserted directly in the structure solver, thus neglecting a part of the coupling difficulty. However, to obtain accurate flow, a large number of particles with a small area of influence should be modelled. A trade-off between accuracy and computational efficiency must be found. Fluid-structure interaction of the mitral valve was performed using SPH and geometric valve obtained from µ-CT scans [START_REF] Toma | High-Resolution Subject-Specific Mitral Valve Imaging and Modeling: Experimental and Computational Methods[END_REF], [START_REF] Toma | Fluid-Structure Interaction and Structural Analyses Using a Comprehensive Mitral Valve Model with 3D Chordal Structure[END_REF], [START_REF] Toma | Fluid-Structure Interaction Analysis of Ruptured Mitral Chordae Tendineae[END_REF]. A similar set-up was also used in a realistic left ventricle [START_REF] Mao | Fully-Coupled Fluid-Structure Interaction Simulation of the Aortic and Mitral Valves in a Realistic 3D Left Ventricle Model[END_REF].

Mitral valve FSI difficulties

In general, fluid-structure interactions are more difficult and more computationally expensive to run than the addition of the fluid and the structural problem. The difficulties linked to the structure model were details previously (chapters 3 and 2). The fluid difficulties must be studied with specific discretisation and its numerical resolution. These are discussed in the section 4.2. This section details only the coupling difficulty related to the fluid-structure simulation of the mitral valve.

The coupling strength between the fluid and the structure was quantified using dimensionless analysis of solid and fluid equations. Details of the analysis are given in the annexe 6.4, based on a generic approach for fluid-structure coupling [START_REF] Langre | Fluides et Solides, Ellipses[END_REF].

The kinematic condition coupling the fluid and solid at the interface is ū = D ∂ ξ ∂ t , with u the fluid velocity, ξ the solid displacement, D the displacement ratio. All variables are dimensionless and use the solid as a reference. For the mitral valve, the order of magnitude for the fluid and solid variables are given in the table 4.1. Thus, the kinematic condition can be simplified to:

ū ∂ ξ ∂ t
The dynamic condition obtained from the dimensionless analysis at the fluid-solid interface is M (-pId + 2 Ur Re d).n = σ.n, with M the mass number, p the pressure, U r the reduced velocity, Re the Reynolds number, d the deformation tensor, n the normal at the interface, σ the Cauchy stress. All variables are dimensionless and use the solid reference. Using the value for the mitral valve (table 4.1) we have:

-pn σ.n Fluid variables values Solid variables values L f 2.5 cm L s 2.5 cm ρ f 10 3 kg.m 3 ρ s 10 3 kg.m 3 U f 1 m.s -1 E 100 KPa µ 4 × 10 -3 Pa.s ξ 0 0.5 cm U s = E/ρ s 10 m.s -1 t f = L f /U f 2.5 × 10 -2 s t s = L s /U s 2.5 × 10 -3 s Dimensionless variables values U r = U f /U s 1 × 10 -1 M = ρ f /ρ s 1 D = ξ 0 /L s 0.2 Re = ρ f U f L f /µ 2.5 × 10 3

Table 4.1: Values and numerical application of dimensionless variables for the mitral valve problem

As a result, it is not possible to neglect the dynamic or kinematic coupling conditions. This dimensionless analysis proves that no drastic simplification can occur to couple the physic, and the fluid and solid must be solved together.

The reduced velocity U r often characterised the type of coupling in fluid-structure problems. Suppose U r 1, i.e. the fluid-structure is in the domain of aeroelasticity. Thus, t s t f and from the point of view of the fluid, the structure is static. Suppose U r 1, thus t s t f , and from the point of view of the structure, the fluid does not move. In the case of the mitral valve, U r = 0.1, and it is not clear if the fluid can be considered at rest during diastole or systole. Also, the mitral valve is hyperelastic and is not well described with elastic simplification. Young's modulus of the valve leaflets ranges from 40 kPa to 8 MPa (pre and post stiffening values) (see chapter 2, table 2.3). Thus the solid time should vary between 5 × 10 -3 and 2.5 × 10 -4 s. In this condition, the reduced velocity varies between 0.2 and 0.01. The valve closure begins at the transition phase between diastole and systole. During this moment, the fluid velocity can reach 5 m.s -1 . As a result, the reduced velocity ranges from 0.8 to 0.2, and the hypothesis of fluid at rest no longer hold. However, should a hypothesis be made, it would be that the fluid is at rest from the point of view of the structure. Since the reduced velocity is not large or small, the loose coupling can not be applied efficiently.

The fluid coupling to the structure can be seen as adding several elements (spring, mass, damper) to the solid problem. Suppose the fluid to be at rest; the coupling can be seen as an additional stiffness and mass to the structure. The added mass is an additional inertia term to the solid equation and can penalise the resolution of the equation. It is possible to obtain an estimation of the added mass using the frequency approach [START_REF] Rousseau | Calcul de Matrice de Masse Ajoutée sur Base Modale[END_REF]. A formulation for the added mass can also be obtained if the problem is simplified.

If the mitral leaflets are undeformable structure with one mode displacement that can be decoupled in time and space, we have [START_REF] Langre | Fluides et Solides, Ellipses[END_REF]:

E Cin = 1 2 m A q2
With E Cin the kinetic energy of the fluid during the movement, m A the added mass, and q the speed of the solid movement. Thus, the added mass can be represented as the inertia coefficient representing the fluid's kinetic energy with the solid velocity. During mitral valve closure, the kinetic energy of the fluid can be computed from:

E Cin = 1 2 V ol eject ρ f U 2 eject
with V ol eject the ejected volume during closure and U eject the velocity of blood. Thus, the added mass can be written as:

m A = 2E Cin q2
The ejected volume during closure is about 60% of the ventricular volume ( 36 mL), and the average blood velocity is 0.6 m.s -1 . The leaflet velocity during closure (which is not the caracteristic velocity of the solid) is q = 1 m.s -1 . Thus, we obtain an m A 13 g, to compare to the mass of the leaflets m s 2 g.

The added mass is of particular interest as it binds different modal responses of the system. Coupling modal responses and frequencies can increase the difficulty to solve the numerical problem. In addition, the discretisation term added to communicate between the two models (structural and fluid) increase the latency between the physics and leads to numerical instability. The added mass on incompressible fluid can lead to instability in the loosely partitioned coupling and increase the computational cost of strong coupling [START_REF] Fernández | A projection semi-implicit scheme for the coupling of an elastic structure with an incompressible fluid[END_REF]. Causin, Gerbeau and Nobile proved that loosely temporal coupling formulation could be unconditionally unstable under certain conditions linked to the added mass [START_REF] Causin | Added-mass effect in the design of partitioned algorithms for fluid-structure problems[END_REF]. For incompressible flow under strong added mass, the convergence is no longer guaranteed regardless of the timestep used [START_REF] Van Brummelen | Added mass effects of compressible and incompressible flows in fluidstructure interaction[END_REF], [START_REF] Lefrançois | Strongly Coupling Partitioned Scheme for Enhanced Added Mass Computation in 2D Fluid-Structure Interaction[END_REF]. Some condition on the added mass for stability was obtained and used as a conjecture for convergence [START_REF] Lefrançois | How an added mass matrix estimation may dramatically improve FSI calculations for moving foils[END_REF]

m s > m A .
In our study, the added mass can be six times larger than the solid mass. Thus, some care must be used when modelling fluid-structure coupling. Several solutions exist to ensure convergence under a strong added mass effect. One approach is to use a specific time-stepping for the structure that takes into account the added mass [START_REF] Song | A Partitioned Coupling Scheme Extended to Structures Interacting with High-Density Fluid Flows[END_REF]. In addition to the energetic method, frequential analysis can be used to estimate the added mass [START_REF] Blevins | Formulas for Natural Frequency and Mode Shape[END_REF]. Briefly, the eigenfrequency of the structure in a vacuum and in the liquid is measured and the relation f vaccum /f fluid = 1 + m A /m s can be used for each eigenfrequency [START_REF] Benaouicha | Analysis of Added Mass in Cavitating Flow[END_REF].

The second solution to overcome the coupling difficulty and the added mass effect is to use a direct or monolithic solver. Using a monolithic solver ensures no instability related to the partitioned approach when coupling the physics. Indeed, there is no latency when sharing variables, and no staggered response exists between the fluid and the solid. However, the added mass still couples the modal responses of the problem studied, which increase the problem complexity.

Choosing the software to run Fluid-Structure simulation

The first constraint for the modelling of the mitral valve was to use an industrial code that would be compatible with other industrial studies. This section aims to justify our software choice among industrial modelling codes. SPH models are used in complex fluid simulations. FSI model of the mitral valve was run with specific FE-SPH formulation using Impetus AFEA solver [START_REF] Mao | Fully-Coupled Fluid-Structure Interaction Simulation of the Aortic and Mitral Valves in a Realistic 3D Left Ventricle Model[END_REF], [START_REF] Toma | Fluid-Structure Interaction and Structural Analyses Using a Comprehensive Mitral Valve Model with 3D Chordal Structure[END_REF]. The modelling of the fluid using smooth particles was not used previously in our group and seemed to be a limiting factor for the fluid model. Acceptable closure is obtained for the mitral valve, but fluid information such as vortices and the viscous response seems to be overlooked. SPH methodology seems to be compatible with direct temporal coupling, which would neglect the added mass instability. In addition, it highly benefits from GPU parallelisation but is still currently lacking computational power for fluid realism. Overhaul, the technique is more and more mature and could lead in the future to an easier simulation of the mitral valve.

Immersed boundary methods are heavily used in applied mathematics research teams-the methods trade meshing difficulty for interface detection and interpolation difficulties. Several techniques were developed to enforce the interface in the fluid domain, and some research codes are available. IBAMR, a free code based on the IBM formulation, was used to model the mitral valve [START_REF] Gao | A Coupled Mitral Valve -Left Ventricle Model with Fluid-Structure Interaction[END_REF]. The YALES2BIO code is also compatible with IBM to compute mitral valve function [START_REF] Chnafa | Image-Based Large-Eddy Simulation in a Realistic Left Heart[END_REF]. Apart from research, the application of the method seems to increase. A mitral valve model was run using IBM with ABAQUS CEL (v12) and a specific coupling algorithm called Tango [START_REF] Bavo | Fluid-Structure Interaction Simulation of Prosthetic Aortic Valves: Comparison between Immersed Boundary and Arbitrary Lagrangian-Eulerian Techniques for the Mesh Representation[END_REF]. However, Dassault no longer supports the immersed boundary version of the software. Other studies used fictitious domains through penalty coupling in LS-DYNA [START_REF] Lau | Numerical simulation of mitral valve function[END_REF], [START_REF] Kunzelman | Fluid-structure interaction models of the mitral valve: function in normal and pathological states[END_REF]. The temporal coupling is handled in a partitioned manner. The authors use a smaller bulk modulus equal to one per cent of its physical value to help computational efficiency and convergence difficulty. ALE formulations are widely used in industrial fluid-structure software. STAR-CCM+ can be coupled to ABAQUS using ALE formulation to run fluid-structure interaction. Sloshing was studied using STAR-CCM+ in [START_REF] Brandely | Etude du Ballotement de Fluide Dans les Réservoirs à Carburant. Approches Numérique et Expérimentale[END_REF]. A second mesh around the structure can be used to cover the fluid mesh. When remeshing is required, this covering mesh is used to interpolate precisely data from one mesh to another. The coupling algorithm is partitioned (no direct coupling). Another possibility is to use ANSYS fluent and ABAQUS. It uses ALE formulation and can be coupled in an explicit partitioned fashion. A coupling algorithm called TANGO can be used to obtain an iterative partitioned coupling. However, it was previously shown that ANSYS had meshing difficulty when modelling 3D valve [START_REF] Bavo | Fluid-Structure Interaction Simulation of Prosthetic Aortic Valves: Comparison between Immersed Boundary and Arbitrary Lagrangian-Eulerian Techniques for the Mesh Representation[END_REF].

There is no clear answer to what method is best suited to model the mitral valve fluid structure for industrial applications. SPH methods have good coupling capacities and can be used with excellent structural solvers, but the fluid response and details seem to be limited. Standard fluid solvers such as Fluent or Star-CCM+ are available using a partitioned coupling but require a specific approach. They use ALE formulation, which will need remeshing and increase the difficulty of our problem. Immersed boundary methods have promising results from research teams, but the industrial codes have partitioned coupling to fluid solvers, which results in a computational approximation to compute the problem.

The use of monolithic coupling for industrial code implies that the software was initially created for multi-physic purposes. Very few codes offer direct solving of fluid-structure interaction. To the author's knowledge, only COMSOL and ADINA provide this type of multi-physics solver. ADINA, in opposition to COMSOL, can also model large deformation and large strain of the structure during fluid-structure simulation of the valve. This is a requirement for the structural problem. Also, it is possible to define specific user material using Fortran subroutine. ADINA can handle complex structural problems and couple the study to the fluid flow using a direct solving approach. However, ADINA uses ALE formulation for the fluid-structure interaction and thus will require remeshing. To solve our mitral valve model, we choose ADINA because of its strong capacity for the structure model and the monolithic solver for multi-physics resolution. Thus, we will be able to model the fluid-structure interaction of the mitral valve without the difficulty of the added-mass instability.

Material and methods

This section aims to define the material and methods used in the fluid-structure simulation of the mitral valve. Using ADINA, an ALE formulation for the fluid is used, which will require remeshing due to the large displacement of the structure. When remeshing, interpolation errors might occur. In the first section, an energetic balance applied to the whole fluid domain is detailed. In the second section, techniques applied to the fluid domain are presented. The third section details how the fluid domain can be separated in ADINA when the mitral valve is closed. The last section details the time-stepping

Energetic balance

ADINA uses ALE formulation to model the fluid during a fluid-structure interaction. ALE mesh deforms to follow the interface precisely when it moves. If the displacement is important, the mesh can become too deteriorated to continue the simulation. In this case, a new mesh must be created, and data is interpolated from the old deteriorated mesh to the new one (figure 4.6). The remeshing can be the source of two errors. First, the elements of the deteriorated mesh might already be of bad quality and could have led to an error in the Eulerian fields before the remeshing. The second error can come from the data interpolation of the old mesh on the new one. The interpolation could lead to diffusivity and errors.

Several methods can be used to monitor the deterioration of the fluid domain. It is possible to monitor bad quality elements and assess their number based on a threshold of the Jacobian ratio. However, only the geometry of the element is taken into account, and not its effect on the fluid domain. Thus, no errors related to the interpolation can be detected with this method. Another technique is to compute an energetic balance to ensure that at each timestep, the energy on the fluid domain is conserved. This second technique was applied in our model.

To study energetic conservation, we applied Newtons second law to a continuous media, in a domain Ω function of time. Special care must be taken for the energy dissipated by viscous force, and for the non-fixed domain border. Detailed explanations are given in [START_REF] Garrigues | Equations générales des milieux continus[END_REF]. In Eulerian reference, it reads:

E cin = t φ E cin + E mec int + E mec ext
with E cin the kinetic energy, φ E cin the flow of kinetic energy, E mec int the energy of internal work and E mec ext the energy of external work. The kinetic energy can be expressed as:

E cin = Ω(t)
ρv 2 2 dV with v the vector of velocity and ρ the density. The flow of kinetic energy is:

φ E cin = ∂Ω(t) - ρ.v 2 2 .(v -v ∂Ω(t) ).n(t)dS
with n the normal of the surface oriented towards the outside of domain Ω, ∂Ω(t) the non-fixed border of domain Ω(t) and v ∂Ω(t) the speed of the domain border. The energy of internal work is defined as:

E mec int = t Ω(t) -σ σ σ : DdV dt
where σ is the internal stress tensor and D = (∇v + (∇v) t )/2 the strain tensor. Finally, the energy of external work results from the action of the forces to the domain Ω:

E mec ext = t Ω(t) v.f V dV dt + t ∂Ω(t)
v.σ σ σ.ndSdt with f V the sum of volume forces applying in Ω (such as gravity).

In case of Newtonian fluids, the stress field can be written as σ σ σ = -pId + 2µD + η∇.vId. In addition, the fluid is assumed to be incompressible. The energy of internal work can be rewritten as:

E mec int = t Ω(t)
2µD : DdV dt It represents the energy lost by viscosity. We also neglect all volume forces f V = 0, thus the energy of external work can be simplified:

E mec ext = t ∂Ω(t) v.σ σ σ.ndSdt
It represents the energy produced by external forces on the boundary.

In general, numerical simulations do not totally respect the physical conservation of energy due to the discretisation used to model the fluid flow. The above formula was used to estimate the precision and the realism of the simulations we ran. Remeshing techniques and space and time discretisations were chosen in order to minimise the energy loss.

Remeshing technique used with ADINA

When a mesh contains distorted elements, it is necessary to rebuild a new mesh to continue the simulation. In ADINA, the steered adaptive meshing (SAM) procedure can recreate free form meshing in this situation. SAM procedure works in two distinct steps. The first step consists of detecting the elements that do not respect the criteria. The second step is remeshing the marked elements and then interpolating the variables from the old mesh to the new one. The criterion used in the first step are:

• All elements in a sphere or in a brick • All elements whose size is not between size min and size max • All elements whose quality is not between quality min and quality max • All elements in a zone who does not respect ∇φ × size element = cste, φ being an eulerian field at the time of remeshing

• Any combination of the above criterion

When the elements are marked during the first step, they are also associated with a target element size. When the first step is finished, the new mesh is generated based on the target sizes of the marked elements. Then the variables are interpolated using a second-order scheme. Consecutive interpolations of the meshes might lead to error diffusivity. Therefore, a second-order accurate interpolation method is used for data from the previous mesh to the new one in order to decrease this error. In figure 4.7, two criterion are used. The first criterion is based on the gradient of pressure and vorticity fields and permits to mesh with smaller elements where the gradients are the strongest. The second criterion is based on the belonging of the elements to an area in the fluid domain. In these , the mesh is of bad quality due to deformation (Jacobian Ratio < 0.7). The deformation continues, and at t 2 , the computation abort because of the mesh distortion. At t 3 , a new mesh is built, and data are interpolated from the old mesh to the new one. This is the default behaviour of the SAM procedure. However, the mesh is of bad quality from t 1 to t 2 , and the interpolated data on the new mesh are obtained from a mesh of very bad quality.

The desired behaviour of a remeshing algorithm is to remesh elements when they start to be too much deformed. However, in ADINA, the elements can only be remeshed when the SAM procedure is called and if the elements are marked by the criterion. The SAM procedure is called when convergence is not reached when the mesh is too deformed and need to be repaired. It can also be called at specified timesteps or after a given amount of time T SAM (figure 4.9). Passively waiting for the procedure to be called results in poor meshes quality, and so in additional errors. A more active approach must be used. The SAM procedure needs two previously converged timesteps to compute the new mesh and interpolate the Eulerian field. Thus, the timestep for remeshing is at least two times greater than the timestep used during the simulation. However, the leaflets can move so quickly that two consecutive remeshing are required, and the SAM procedure will fail in remeshing as the previous timestep was also on a new mesh. This can be difficult to deal with as the deformation of the structure (and thus of the mesh) is not known in advance and depends on the simulation. As a result, the timesteps is related to the fluid and solid element deformation, which is also related to their sizes. Thus, the study of mesh size is intrinsically related to the temporal resolution of the simulation and its remeshing.

Discretization of the problem

Time discretisation

Using a monolithic (or direct) coupling method will enforce the same time discretisation for both fluid and solid problems. The same time scheme and timesteps are thus used.

The time scheme used is a composite and implicit method. The solution at time t + ∆t is obtained in two consecutive sub-timesteps:

u(t + γ∆t) = u(t) + γ∆tf (u(t + γ 2 ∆t))
for the first sub-timestep, with γ = 2 -1 α . The second sub-timestep is:

u(t + ∆t) = (1 -β)u(t) + βu(t + γ∆t) + (1 -α)∆tf (u(t + ∆t))
with β = α 2 2α-1 . This method is unconditionnaly L-stable and of second order if 1 2 < α < 1. We used the default configuration, α = 1 √ 2 . This lead to minimal truncation error for linear system.

Spatial discretisation

The weak form of Navier-Stokes equations is obtained using the Galerkin method. The equations are weighted with virtual functions (pressure, velocity or displacement) and integrated over their respective domain. The fluid and solid stress derivatives are then lowered with the divergence theorem, and boundary conditions naturally appear. The virtual functions and the unknowns can be discretised on several bases. If the virtual functions are defined by steps functions, we obtained a discontinuous Galerkin discretisation, called Petrov-Galerkin. This resulting formulation is a hybrid of finite element method and finite volume method and requires using Flow Condition Based Interpolation (FCBI) elements for the fluid domain [START_REF] Kohno | A flow-condition-based interpolation finite element procedure for triangular grids[END_REF]. If the virtual functions are defined on the continuous polynomial Lagrange base, we obtain the finite element method (with Lagrange finite elements). We consider only incompressible fluid. During fluid-structure analysis, a finite element method for the solid and a Petrov-Galerkin formulation for the fluid along with FCBI elements was used. FCBI elements are better than classical FEM elements to ensure local momentum and mass conservation. 2D three-node triangles FCBI elements (first-order accurate) are used to discretise the fluid domain (figure 4.10b). The solid is discretised with 2D nine node quadrilaterals (second-order accurate) (figure 4.10a). Looking now at the transient evolution of energy during the simulation Figure 4.12, we can see that most of the energy is lost during the shift of input velocity, at 0.3 s. By implementing a remeshing based on vorticity and pressure gradients, the velocity variations can be more easily followed if the number of remeshing during a pulsation is sufficient. Using a remeshing criterion based on the gradient of vorticity and pressure, and limiting the size of the elements between 0.2 cm and 0.025 cm, the same fluid only model was launched with repeated calls to the SAM procedure (figure 4.14a). We observe that the relative energetic loss decreases when the number of remeshing increases, down to a minimum reached at 50 remeshing. From 50 calls, the fluctuation of the velocity can be well captured with the remeshing technique, and the energy is well conserved during the simulation. During this simulation, the element size varies from 0.2 cm to 0.02 cm, even if a limit was set to 0.025 cm. This explains the excellent conservation of the total energy, performing even better than the uniform finer mesh with a given (and verified) element size of 0.025 cm. Using 50 calls to the SAM procedure during the fluid simulation resulted in an energetic loss smaller than 2% and was considered acceptable. We note that using more than 50 calls to the SAM procedure requires a timestep smaller than 10 -2 s. This is because the SAM procedure requires for every new mesh a previously converged timestep on the same mesh. The CPU time was also modified positively by the frequent use of the remeshing technique (figure Using regular remeshing with criteria based on vorticity and pressure gradients decreases the energy lost during the simulation from 10% to 3%. This is because the mesh is denser where more variations exist, thus improving the node repartition compared to uniform meshes. Using this remeshing several times during a simulation also decreases the computational time. The minimum is found at 20 calls but was nearly constant until 100 calls.

Definition of the timesteps and the mesh sizes

On the fluid-structure simplified model, the leaflet will move and deform the mesh. Our objective is to determine how many times the remeshing should occur during a simulation to prevent the formation of an element of bad quality. Thus, the fluid-structure simulation was run one time with the SAM procedure called only to repair the elements. During the simulation, the frequency at which the SAM procedure is called was monitored (figure 4.15). Thus to repair the fluid elements, the procedure must be called at least 80 times per second. In the simplified model, we used a regular call to the SAM procedure with a frequency of 100 Hz, ensuring good quality for the elements.

The fluid-structure model was discretised using the SAM procedure on a periodic flow. We are

Simulation of domain closure

Another difficulty when simulating mitral valve dynamics is to model an effective closure of the fluid domain. During systole, a contact is detected on the structural model of an FSI, but this information is hardly ever given to the fluid problem to modify the boundary conditions. Detection of the closure aside, modifying the fluid domain topology during the simulation can be difficult for the solver. The topological change of the fluid domain must be handled properly when the domain is separated and is reunited to prevent leakage or pressure jump.

Gap condition

In ADINA, a "Gap" boundary condition was designed specifically for this task. Defined on a geometry interface (not necessarily physical), it detects closing area depending on values reached by a criterion (area of the interface, time or pressure). Once activated, the nodes at the interface are duplicated for both sides. The created nodes replace the nodes at the interface and are used by the elements instead. When the "Gap" condition is deactivated, the nodes of the two sides of the interface are merged along with the corresponding elements. Simple tests cases were developed to test the GAP condition. The first one is a pressure-driven flow in a closed cavity with an obstacle. A constant pressure difference of ∆P = 180 P a is imposed through normal traction. The simulation methods were detailed previously. The timestep is ∆t = 2 × 10 -2 s, mesh size varies between 0.2 cm and 0.025 cm, and SAM frequency is 10 Hz. At time t = 1s, the domain is closed using the gap condition (figure 4.17). The velocity and pressure field obtained are given in figure 4.18. After the activation of the GAP condition, no residual leak appears in the velocity or pressure interface (figure 4.18). Post-processing the flow after the GAP activation gives no residual flow after activation under constant pressure. and the SAM remeshing technique. In this test, the solid inside the fluid domain is moving periodically up and down and right-to-left. The GAP condition is set on the projection of the solid edge to the closest boundary wall. As the solid moves, the line is recomputed for each timestep to be the shortest distance to the boundary. The GAP condition is set on the line length (figure 4.19). As a result, when the solid approaches the wall, the residual space is closed. The other parameters are the same as the first test.

The second test case shows that the GAP closure can be used with a geometric parameter (figure 4.20, Velocity and Pressure). When the solid moves away, the separated fluid merge together and can cross the GAP line. Once again, the delimitation of pressure and velocity is absolute when the GAP is active, i.e. there is no pressure or fluid leak. This test also shows that the GAP condition is compatible with the SAM remeshing technique before, when active, and when the GAP condition is inactive (figure 4.20).

Pseugo-GAP condition

One default of the GAP technique is that the activation of the GAP is restricted to the boundary condition area, a given time or a pressure difference (only for opening). The pressure difference is restricted to opening the boundary and can not close it. Thus, it is not applicable to the mitral valve. Implementing the GAP condition based on time would neglect all the healthy or pathological valve structural responses and thus is not an option. Thus the GAP condition must be activated using the area of the boundary condition. This is efficient for the mitral valve as long as the shortest distance between the two leaflets can be used as the GAP condition. It is also realistic, as the same distance can be used for the structural offset contact algorithm and the GAP activation. However, this shortest distance is not known prior to the simulation. The method previously used to compute the shortest distance to the wall can only work on a fixed boundary. It is not possible to compute this distance with the mitral leaflets.

An alternative method would be to determine the location of the shortest distance and then implement the GAP condition at this location and rerun the analysis. A structural response under uniform pressure can be used to obtain the deformed shape of the leaflets. However, the structural deformation obtained would lack some information that is available in the fluid-structure simulation.

Another method is to use the elements of the FCBI to naturally enforce the GAP condition (figure 4.21). The idea is to use the formulation of the FCBI elements. The information is only computed at the nodes of the elements. Thus, when the same element is connected to two walls, all the nodes of the element have velocity enforce to 0 because of the no-slip boundary condition. As a result, no flow across the elements is tolerated. As a result, and as long as the element is connected to two different walls, no fluid can cross the element.

This condition is based on the element size and the distance between the walls. In the case of the mitral valve model, this limit the fluid mesh size to the possible structural distance between the leaflet. In other words, the fluid mesh size must be larger than the structural contact offset.

Resolution of equation and solver used

The fluid problem is non-linear due to convection. The fluid-structure simulation also has the nonlinearity of the structural problem, in this case, due to the contact model. The equations must be linearised before being solved. The Newton Raphson algorithm is used to linearise the system and solve the non-linear equations. The linearisation was enforced in a maximum of 50 steps. Relative tolerance was fixed to 10 -3 for all degrees of freedom. When the estimated errors were lower than the tolerance, the solution was considered converged, and the resolution proceeded to the next timestep. To solve the linearised equations, ADINA can use several solvers. Most of them are iterative (GMRES (Generalised Minimal RESidual method), BCG (Bi-Conjugate Gradient method)) and can be preconditioned. There are also two direct solvers, "SKYLINE" and a "sparse" solver, using the Gauss elimination method. This type of solver is less sensitive to the ill-conditioning of matrices than iterative solvers and offers an exact resolution of the discretised equations. Direct solvers are, in general, slower and tend to consume more RAM than iterative ones when correctly used. Indeed, an exact solution is more difficult to obtain than an approached solution. However, using an iterative solver provide negligible benefit when the size of the problem is small and add extra errors to track down.

In this study, the sparse solver was used. It is a direct solver, so there is no need to monitor the error during the resolution of the equations. The "sparse" solver uses a classical LDU decomposition of the matrix A to solve Ax = b, x being the unknown and b the right-hand side. Hence, the problem is solved with x = U -1 D -1 L -1 b. The main point of the "sparse" method is the preservation of the sparsity of the matrix A, reducing required storage and decreasing computer time by order of one to two, in comparison with traditional direct solvers ("SKYLINE" for example). The matrix to invert for FSI problem can be very large, especially when the monolithic approach is used. However, this resolution is available with sparse storage of the matrix and thus drastically reduces the size of the problem.

Backward flow

One of the difficulties of mitral valve simulations is to account for the backward flow. If the mitral valve model is leaky, a flow going from the ventricle to the atrium will appear during systole and could leak to the boundary. In comparison, a healthy mitral valve will direct flow from the atrium to the ventricle, and no backward flow should appear.

If velocity is implemented as the input, the mitral valve behaviour would be imposed, and the health of the valve would not be tested. As a result, a physically leaking valve would not present any leak during the simulation. Thus, the mitral valve behaviour must be tested using a flow driven pressure. In addition, the pressure-driven flow can not have imposed input and output. During diastole, the blood flow is oriented from the atrium to the ventricle. During systole, this orientation should be conserved when healthy valves are used. However, if the valve leaks, the blood flow could be reversed and pushed from the ventricle toward the atrium. Thus, during the same simulation, it is possible that blood flows from the input to the output boundary and then from the output to the input boundary.

To test if the software is able to model pressure-driven flow with no specific direction, a simplified , with e = (∇u + ∇ T u)/2 the velocity strain tensor and ν the viscosity. The mesh is defined using the SAM method based on vorticity and pressure gradients (0.2 cm ≥ ∆x ≥ 0.025 cm). Time scheme is second implicit composite with ∆t = 10 -3 s. The fluid is the same as in the previous study.

The pressure and velocity are monitored during the two seconds and plotted along the atrium and ventricle boundaries (figure 4.23). The average pressure is well defined, which validates the pressure implementation using the normal traction. As a result, fluctuations of the minimum and maximum pressure appear when the average flow is important. The flow of blood is orientated from the atrium to the ventricle during the first second, which corresponds to a diastolic behaviour. From t = 1 s to the end of the simulation, the flow orientation is reversed, and blood goes from the ventricle to the atrium (figure 4.24). This behaviour corresponds to a pathological valve that would not close. In this model, the leaflets are rigid and do not move. Thus they can not close the valve. The leaking volume of fluid can go through the atrium border. In addition, we can see that a recirculation appears on both boundaries due to local pressure variations induced by the flow. This demonstrates that the software is able to accept a recirculation area on the boundaries. Also, the software can model flow independently of the input/output chosen.

2D simulation of the mitral valve

A 2D model of the mitral valve was implemented with the method described in the previous section. To summarize: -ALE formulation is used to compute the moving structure in the fluid. Remeshing is required due to the large displacement of the structure. -The remeshing is carried out with repeated calls to the SAM procedure at each T SAM . The new mesh is computed from the gradients of pressure and vorticity.

-The fluid mesh is subject to error due to the regular remeshing. An energy balance is used to monitor the error linked to the remeshing.

-The flow is driven by pressure, and no input or output is specified. Depending on the structure behaviour, backward flow can appear. Pressure is imposed through normal traction.

-The mitral leaflets can not overlap in the fluid domain. The domain closure is obtained from a GAP boundary condition or a "pseudo" GAP technique.

-The ADINA direct coupling method is used to bypass the multi-physics difficulties.

Model definition

A rigid simplified left heart was used to represent the mitral valve's environment. The initial model was taken from [START_REF] Boilevin-Kayl | Modeling and numerical simulation of implantable cardiovascular devices[END_REF] and modified to represent more accurately the aortic root. The geometry and boundary conditions are depicted in figure 4.25. Maximal dimensions are 8.5 cm and 5 cm. The leaflets are symmetrical, of approximated length 2.3 cm and thickness 0.65 mm. Both leaflets are fixed on the ventricular walls and attached by one chordae tendinea to the papillary muscle tips.

The fluid and solid domains are Ω f and Ω s . Ω s moves inside Ω f . Ω f is fixed except for its internal boundary with Ω s . A no-slip condition is applied on the immobile rigid walls and on the moving deformable valve. A difference of pressure is applied on the domain to mimic the systolic and diastolic phases. The atrium border is imposed 0 Pa pressure, and the aortic root is imposed a physiological aortic pressure variation of 120 mmHG ( 16 KPa) (figure 4.25b).

The physical parameters used for the fluid and solid are depicted in table 4.3. Blood is approximated as a Newtonian, incompressible fluid with constant viscosity µ = 3 mPa.s. The Reynolds number is Re = ρV L µ 2000), so the flow can be considered as laminar. The leaflet is assumed to be a small strain, even if the displacements are large. This hypothesis allows using a classical isotropic elastic constitutive law. Strain reaches 20%; thus, a large strain formulation should be used instead. This assumption does not limit the coaptation of the leaflets in two dimensions. The Young modulus of the valve is set to be E = 6 × 10 5 Pa and the Poisson coefficient to ν = 0.4. Gravity is neglected. The two chordae tendineae have a section area of 6 × 10 -3 cm. They follow the non-linear stress-strain relation defined in ??. The contact between the two leaflets has a 0.02 cm offset and a compliance factor of 10 -8 . These values ensure a smooth area of contact and forbid the leaflet fluid boundary to connect or overlap.

Solid elements are second-order nine nodes bricks in a mixed U/P formulation (where pressure is interpolated separately) with three pressure degrees of freedom and are of an approximate size of 0.1 cm. Chordae tendineae are modelled with one truss element. The fluid elements are first-order three nodes elements called FCBI and can model flow using step functions.

The simulation was run twice to locate the closest nodes between the two leaflets during systole. In the first simulation, a uniform coarse mesh was used (element size fixed to 0.1 cm). The fluid elements connected to both leaflets forbid fluid flow, as the fluid velocity is only computed at the nodes of the element. This "pseudo"-GAP technique was used to virtually seal the ventricle during systole. The simulation was then run again with a finer mesh based on the gradient method (figure 4.27). The seal

Parameter

Value 4.26. The average pressure on the aortic boundary is strictly respected, and no variations can be observed. At the atrium boundary, variations of the imposed pressure appear. The amplitude of the variations is small when compared to the imposed pressure on the aortic root (0.5%). The variations are linked to the high velocity located at the boundary through the use of normal traction conditions. Overall the pressure is well respected.

ρ f = ρ s 1 g.cm -3 µ 3 × 10 -2 g.cm -1 .s -1 E 6 × 10 6 dyne.cm -2 ν 0.4

Results obtained

Meshes, velocity and pressure magnitude variations are depicted at several timesteps during the cardiac cycle (figure 4.27, 4.28 and 4.29 respectively). The results are presented at the 4th pulsation to limit the influence of the initial condition (the fluid is initially at rest). Results are given at specific time (figure 4.25b).

Meshes obtained using SAM procedure

The meshes are generated from the SAM procedure, with T SAM = 10 -2 s and 0.2 cm ≤ ∆x ≤ 0.025 cm. Close to the boundary, an area must be enforced with a given mesh size to post-process the input and output. The meshes follow the variations of pressure and vorticity to refine the mesh where information The elements are generally coherent with the velocity and pressure variations. However, at the area closest to the boundary, an abrupt transition between small and large elements is observed. This is due to the boundary implementation. The boundary elements can not be remeshed. Thus, the area closest to the boundary is given a specific size to move the remeshing area away from the boundaries.

Velocity field

The velocity field on the simplified left ventricle can be observed in figure 4 other, and the blood flows again from the atrium to the ventricle. In our model, the aortic valve is not modelled. Blood flow is not entirely forced in one direction during the whole cardiac cycle. Thus, when the mitral valve opens at the end of the systole, a small flow appears from the aortic root toward the atrium (figure 4.28,g). This flow will be avoided if the aortic valve is closed.

Pressure variations

Pressure variations during the cardiac cycle are presented figure 4.29. A clear delimitation of pressure can be seen at the GAP line. When active, the GAP condition divides the left ventricle into two areas. The atrium has nearly constant pressure. This was also depicted from the pressure measurements at the boundaries (figure 4.26). The ventricle pressure follows the imposed pressure at the boundary. The clear delimitation of pressure ensures that the two zones (atrium and ventricle) are disconnected, and thus no fluid is going through the leaflets. Thus, the leaflets are closed from t = 0.32 s to t = 0.65 s, which correspond to the duration of the systolic pressure variation. The valve is closed during this period of time. The velocity field variation above the leaflets seen in figure 4.28 c,d and f are the results of the movement of the leaflets.

Conclusion and Perspectives

Conclusion

Using ADINA, we have demonstrated that it was possible to obtain full closure of the fluid domain and use pressure-induced flow to model the mitral valve behaviour in a simplified 2D left heart. The ALE formulation was used for the fluid. Remeshing is required to model the large displacement of the valve in the ventricle. ADINA has several types of remeshing available. A gradient-based remeshing technique was used and seemed promising to refine the mesh where information varies the most. Possible errors made during the creation of the new mesh are monitored using an energetic approach.

ADINA can model fluid-structure problems using a direct (or monolithic) coupling method. This approach was chosen to neglect all coupling difficulty that was depicted in the literature. In particular, the fluid is often considered compressible, which might lead to reduced peak stress during systole. However, the direct coupling method can be RAM intensive. In ADINA, the resolution of the multiphysics problem can use sparse storage of the matrix involved, reducing the size of the problem.

The methods developed were applied to a simplified 2D left heart under physiological pressure variations for several cardiac cycles. The velocity and pressure results obtained at the fourth pulsation show that the mitral valve closed well during systole and that no backward flow was observed once the valve was closed. A regurgitant jet can be observed when the valve closes under increasing pressure. Both the inflow velocity and the maximal regurgitant speed are similar to physiological results. The results obtained from this simulation encourage us to follow the study toward a more complex and more realistic multi-physics valve simulation.

Perspectives

The simulation presented in this chapter could be improved in several ways. The domain used for the ventricle could be extended towards the aortic root to move the area of fixed element size away from the ventricle chamber. This would be beneficial to capture more accurately the hemodynamic of the ventricle. Indeed, the area with an element of fixed size neglects the information coming from the pressure and vorticity gradients. The loss of information is important on the aortic root, and efforts toward a correct meshing method in this area should be prioritised.

More realistic boundary conditions would also be beneficial for the ventricle and mitral realism. For example, the aortic valve could be modelled using a GAP function based on a time condition. Thus, the reverse flow observed during the opening of the mitral valve would be reduced or cancelled. If the aortic valve was modelled with the GAP condition, the pressure could not be imposed with normal traction on the aortic valve. Another method must be found to implement the pressure variation inside the ventricle. Physiological pressure variations come from the displacement of the ventricular walls. This displacement can be implemented in the software, but it can be tedious to monitor pressure inside the ventricle using only displacement of the fluid boundary.

The solid used in the simulation is composed of simplified leaflets and the chordae tendineae. The chordae tendineae are modelled as truss and have no interaction with the fluid. They have a non-linear stress-strain relationship. The simplification of the chords is common in the literature. The leaflets, however, were simplified to linear isotropic material. This simplification tends to produce leaflets that resist coaptation in 3D. In 2D, the coaptation can be observed. However, the dynamic response of the leaflets is modified. Some oscillation due to the elastic response appeared when the mitral valve closed. Using hyperelastic response would limit these oscillations and increase the closure realism of the mitral valve.

The present study is limited to a 2D case. A fluid-structure simulation of the 3D valve must be performed to study the valve response and closure accurately. In principle, the same technique could be used to compute the 3D problem. However, the number of elements could increase from the 100 thousand elements used to several million. This can be difficult to deal with, especially when a direct approach is used. Meshes of progressively increasing size should be used to approach this difficulty. The SAM method in ADINA can limit the number of elements used in each criterion. The resolution of large problems can be difficult, even with the sparse L -U decomposition used. In this case, iterative solvers are preferred because they tend to make better use of multi-threading.

Chapter 5

Conclusion and Perspectives

General conclusion

This thesis deals with the modelling of the mitral valve. The ultimate goal is to use the numerical tools developed to help the design and optimization of mitral valve medical devices. To provide this numerical platform, mitral valve geometries were first obtained. The mechanical response of the mitral valve components were then tested, implemented and fitted from experimental data. The structural response of the valve was finally obtained in static and dynamic studies. Fluid-structure of the valve and the surrounding blood was then tested in 2D simulations.

The geometric modelling of the mitral valve was addressed in the first chapter. The complex mitral valve geometry is not easily obtained from in-vivo scanning technique. Ex-vivo scans such as µ -CT scans perform better to obtain the mitral shapes. However, the model can not be used to develop patient specific medical devices. Thus a parametric model was developed and used with anatomical data to create mitral valve shapes. This model can be used with geometrical information obtained in-vivo. It can represent a large variety of pathological and healthy mitral valve.

The second chapter deals with the valve constituent and their mechanical response to deformation. The mitral valve parts are composed of fibers and a surrounding matrix. The mechanical response of the fibers and of the matrix largely differ. As a result, the mitral valve has strong hyperelastic and anisotropic behaviour. Isotropic material laws were first developed and tested to represent the valve. When used to represent the closed state of the valve, these laws were not always realistic. The elastic model failed to obtain a coherent coaptation area. The Mooney-Rivlin model showed a significant protuberance and was not realistic. Only the Sussman-Bathe model lead to a realistic closure of the valve. However, numerical difficulties arise when the Sussman-Bathe material law was used. To detail this difficulty, a stability code was develop for hyperelastic and anisotropic material. The code depicted possible instabilities with the Sussman-Bathe and Mooney-Rivlin model used. Hyperelastic laws were developed for the mitral leaflets that respect the stability criteria. These laws were tested under equibiaxial numerical test and fitted closely experimental data. They were coded in ADINA using user-coded subroutines, hybrid U/P formulation and Total Lagrangian formulation.

The third chapter deals with the structural simulation of the mitral valve. Structural response of the valve is challenging because of the non-linearity of the material law, the large deformation and large strain of the leaflets and the contact of the two leaflets. Implicit time-stepping method is used for compatibility reason and to limit the computational cost of the simulation. Shell and 3D-elements were used to obtained static closure of the mitral valve with the Sussman-Bathe material law. The static closure were obtained using a dynamic method and a strong damping to limit the computation cost. A dynamic response of the valve was obtained with the Sussman-Bathe material with shells. Partial dynamic results was obtained with one anisotropic hyperelastic law.

The fourth chapter deals with the interaction of the valve and its surrounding blood. The coupling between the two physics is complex, and could lead to added-mass effect instabilities. To prevent this, monolithic coupling of the structure and the fluid was used. An ALE method is used in the software ADINA. The structure undergoes large displacement in the fluid domain and it was necessary to remesh. Remeshing technique based on the gradient of the variables was used to provide both correct mesh and local refinement during the simulation. To monitor the solution and possible errors due to the remeshing interpolation, an energetic balance was used. Two techniques were implemented and tested to separate the fluid domain during the valve closure. The methods were applied to a 2D mitral valve immersed in a left heart. Vortices was observed and excellent closure was obtained during systole.

Perspectives

The structural response of the valve lead to interesting results in term of stress and coaptation area. Nonetheless, most structural simulations of the mitral valve were difficult to perform. The non-linearity of the material law is the main issue. It would be of the most interest to perform structural simulation with the laws of Prot, Gao or May-Newman and Yin, as they all represent accurately the hyperelastic and anisotropic response of the leaflet. These laws are highly non-linear due to the exponential form in the strain energy density. To help the convergence and soften the problem, smaller timesteps could be used. However, timesteps down to 10 -6 s were used for the computation, to no avail. If we use the CPU time for the first time steps as reference, and with ∆t = 10 -6 s, the CPU time for one cardiac cycle would be around 1000 days, which is not reasonable. Explicit method could in theory perform this simulation with even smaller timesteps, and smaller CPU time. However, in the software ADINA, it is not possible to use usercoded material model with explicit integration scheme.

The material laws are implemented using a Total Lagrangian, as it seems to be the simpler and more direct method. Nonetheless, it is possible to use Udpated Lagrangian formulation. In this updated formulation the displacements unknwown are increment from a referent configuration at time t, instead of the initial configuration for the Total Lagrangian formulation. This could help the resolution of the model.

A coherent fluid-structure interaction was obtained in two dimensions. This interaction shows realistic displacement, pressure and velocity field, but is of limited interest to test medical devices. The development of a 3D fluid-structure interaction would be key to quantify accurately the effect of medical devices on the mitral regurgitation. This 3D model would help to assess the pathology and the effectiveness of medical treatment.

The parametric geometrical model developed represent realistic mitral valve. The leaflets in particular can be shaped to represent in-vivo images well. Nonetheless, the chords could be more realistic in their implementation by adding strut and basal chordae. The fibers used in the anisotropic models are simplified, and not adapted to the leaflet's geometry. An inverse problem could be solved to orientate the fiber according to the geometry of the leaflets and would increase the realism of the mechanical response [START_REF] Toma | High-Resolution Subject-Specific Mitral Valve Imaging and Modeling: Experimental and Computational Methods[END_REF].

The long term goal is to help the development of medical devices on the mitral valve. The development of the mitral valve dynamic environment would benefit the realism of the valve dynamic. The annulus cyclic deformation along with the ventricle contraction could help in assessing the performance of medical devices more accurately.

109

Chapter 6

Annexes

Hyperelastic material are more difficult to model than linear material, commonly defined with two parameters (Young's Modulus and Poisson's coefficient). The theory used to define hyperelastic material in finite element formulation is summarized in the first section. The second section define material behavior based on displacement and pressure, in opposition to the classic displacement only based formulation. The third section explains the implementation of U/P formulation in most finite element codes using a Total Lagrangian approach. Unless stated otherwise, α i β i means i α i β i and φ , j means ∂φ ∂x j . Derivatives with respect to ordinary matrices are noted ∂φ

∂X = ∂φ ∂X ij ij
, and for symmetrical matrices the convention ∂φ ∂X =

1 2 ∂φ ∂X ij + ∂φ ∂X ji ij
will be used, according to the literature [START_REF] Sussman | A finite element formulation for nonlinear incompressible elastic and inelastic analysis[END_REF], [START_REF] Green | Theoretical Elasticity[END_REF], and [START_REF] Srinivasan | What Is the Gradient of a Scalar Function of a Symmetric Matrix ?[END_REF] (preprint).

Summary of hyperelastic theory

Material definition in large strain and large displacements as well as its implementation is well described in several books on finite element method or hyperelastic deformation [101], [START_REF] Holzapfel | Nonlinear solid mechanics : a continuum approach for engineering[END_REF].

Consider a continuum body Ω in a Cartesian space, figure 6.1. At time 0 and a particule M inside this body, defined at the location 0 x. At time t, the same particule M is at position t x. Ω being continuous, φ( 0 x, t) = t x is a continuous function. We have:

t 0 X = ∇0 x (φ) = ∂φ ∂ 0 x
with t 0 X the deformation gradient, defined as:

∀ i, j t 0 X ij = ∂ t x i ∂ 0 x j = t 0 x i,j
Displacement u is defined as t 0 u( 0 x, t) = t x -0 x The notation t 0 (.) denotes a variable at time t, referred to the original configuration time 0. Any deformation gradient can be decomposed into a rotation and a positive definite stretch tensor, fig 6.2. t 0 U (respectively t 0 V) is the right (respectively left) stretch tensor, and t 0 R is the rotation tensor. We have t 0 X = t 0 R t 0 U = t 0 V t 0 R. The stretch tensors are positive definite and have the same eigenvalues, the principal stretches t 0 λ i . The eigenvectors of t 0 U (respectively t 0 V) are the principal direction of stretches, noted in matrix form t 0 R L (respectively t 0 R E ). Using t 0 Λ = diag t 0 λ i , we have

t 0 U = t 0 R L t 0 Λ t 0 R L T and t 0 V = t 0 R E t 0 Λ t 0 R E T , with t 0 R E = t 0 R t 0 R L .
The rotation t 0 R L represent the rotation of the coordinate system towards the principal directions of deformations.

The rotation of a body does not induced any strain, so rotation independent measures of the deformation are preferred. Cauchy-Green and Green-Lagrange strain tensor are rotation-independent tensor. The right Cauchy-Green tensor C is defined as: t 0 C = t 0 X T t 0 X = t 0 U 2 and does not depend of t 0 R. The Green-Lagrange strain tensor t 0 E is defined as:

t 0 E = 1 2 t 0 C -I
where I is the identity matrix. Hyperelastic material are material that can be defined using a strain energy density function t 0 W . By definition, this function relates to the second Piola-Kirchhoff stress tensor t 0 S:

t 0 S = ∂ t 0 W ∂ t 0 E = 2 ∂ t 0 W ∂ t 0 C ie ∀ij, t 0 S ij = 1 2 ∂ t 0 W ∂ t 0 E ij + ∂ t 0 W ∂ t 0 E ji
The second Piola-Kirchhoff has no physical interpretation. It is preferred in hyperelastic theory, because it is easily defined. It relates to the Cauchy stress (also called true stress) t 0 σ:

t 0 σ = 1 J t 0 X t 0 S t 0 X T 112
where J = det( t 0 X) is the Jacobian of the deformation, and represents the volumetric variation. The constitutive tensor of the material t 0 C will be used to define the mixed formulation in the following section. It is defined as:

t 0 C = ∂ t 0 S ∂ t 0 E ie ∀ ijrs, t 0 C ijrs = 1 2 ∂ t 0 S ij ∂ t 0 E rs + ∂ t 0 S ij ∂ t 0 E sr
In practice, the strain energy density function t 0 W is defined as a function of scalars. The scalars, called invariants of deformation, represents types of deformations. They can be defined from the right Cauchy-Green tensor, or the stretches:

t 0 I 1 = Tr t 0 C = i t 0 λ 2 i t 0 I 2 = 1 2 Tr( t 0 C) 2 -Tr( t 0 C 2 ) t 0 I 3 = det( t 0 C) = i t 0 λ i 2 = J 2
For anisotropic material with a single fiber direction, additional invariants based on fiber orientation are used. We note 0 N the unit vector representing the fiber in the undeformed configuration. The fiber in the deformed configuration t n is obtained from t n = t 0 X 0 N . To describe the anisotropic deformation, two additional invariants are used :

t 0 I 4 = 0 N t 0 C 0 N t 0 I 5 = 0 N t 0 C 2 0 N
In this thesis, the only configuration used as a reference is the initial configuration, defined at time 0. Variables, unless stated otherwise, are defined at time t. The prescripts t 0 (.) are not necessary in this context. They are omitted in the following sections. The hyperelastic formulation obtained is referred as Total Lagrangian, in opposition to Updated Lagrangian formulations. Updated Lagrangian formulation use a referent configuration at time t to define deformation at time t + dt, and are also a relevant choice for large deformations and displacements. The two formulations represent the exact same behavior, but results in different numerical implementation.

Mixed U/P formulation in ADINA

ADINA is a displacement-based finite element code, and uses displacement-based interpolation. Like all displacement-based finite element codes, ADINA meet difficulties in the analysis of almost incompressible materials. Symptoms are ill-conditioning of the stiffness matrix, incorrect stress (pressure), and locking behavior [START_REF] Malkus | Mixed finite element methods -Reduced and selective integration techniques: A unification of concepts[END_REF]. Locking is a phenomenon where the interpolation functions of the element can not describe accurately the field modeled. Locking generally appear as the incompressibility of the material increase. Standard approaches to circumvent incompressibility problems are:

• a mixed formulation where a stress variable linked to hydrostatic pressure is interpolated separately from displacement based variables

• the reduced/selective integration where the volumetric terms in the stiffness matrix are integrated using a lower order of Gauss integration than what is used to integrate the other terms When using mixed formulation, the inf-sup condition (or Brezzi-Babuska condition) prevents the apparition of locking. Mixed U/P formulation is available in ADINA and allows to model nearly incompressible material. To respect the inf-sup condition in ADINA, the 9 nodes element with 3 pressure degrees of freedom can be used in 2D. In 3D problems, the 27 nodes element with 4 degrees of freedom also respect this condition.

In this annex, we describe the U/P formulation used in ADINA with special emphasis and the restrictions implied on the material by the formulation. The interest of reduced invariants is explained, and their definitions are given. For a more detailed approach, see [START_REF] Sussman | A finite element formulation for nonlinear incompressible elastic and inelastic analysis[END_REF].

Definition of the U/P formulation

Mixed formulations use an additional scalar field to describe the pressure. In ADINA, pressure is interpolated separately from displacement variables. We note (.) displacement based variables (obtained from displacement field) and (.) pressure based variables (obtained from pressure field).

Finite element method is based on the principle of virtual work:

0 Ω S : δE d 0 ω = 0 Ω ∂ W ∂E : δE d 0 ω = δ 0 Ω W d 0 ω = t R (6.1)
with R the external virtual work at time t and S : δE = Sij δE ij the matrix double product, Green-Lagrange strain tensor being symmetric. The Second Piola-Kirchhoff stress tensor have been replaced with its definition using the strain energy density function.

In equation 6.1, the strain energy density function is computed solely from the displacement field. To correct the locking behavior, an additional term Q is added to the strain energy density function W , so that eq. 6.1 is now replaced by eq. 6.2:

δ 0 Ω W + Qd 0 ω = t R (6.
2)

The modified strain energy density function W + Q must respect several rules:

• Pressure form the modified strain energy density function must equal the separately interpolated pressure p

• The variation of the modified strain energy density function function of pressure must lead to a constraint relation between the pressure obtained from the displacement and the separately interpolated pressure

• When the pressure obtained from the displacement is equal to the separately interpolated pressure, Q = 0

The hydrostatic pressure (positive in compression) computed from displacement is:

p = -1 3 σ ii = -1 3J X im X in Smn = -1 3J C mn 1 2 ∂ W ∂E mn + ∂ W ∂E nm
Similarly, a pressure function is defined as:

P (.) = -1 3J C mn 1 2
∂(.) ∂E mn + ∂(.) ∂E nm (6.3)

The function respect P ( W ) = p, and the first rule enforced P ( W + Q) = p. The simplest constraint obtained from the second rule, the pressure function obtained previously and the condition of the third rule lead to the system: Q = -1 2P (p) (p -p) 2 P 2 (p) = 0 (6.4)

The function P can be rewritten in terms of volumetric variations. Consider d (C ij ) = C ij dα and d (C ji ) = C ji dα. dα is a variation, and can be set as dα = 2 3J dJ. Thus, using the symmetry of C:

d W = 1 2 ∂ W ∂C ij d(C ij ) + ∂ W ∂C ji d(C ji ) = 1 2 ∂ W ∂C ij + ∂ W ∂C ji C ij dα = 1 2 ∂ W ∂E ij + ∂ W ∂E ji C ij 1 3J dJ
From previous expression we have:

p = P ( W ) = - d W dJ (6.5)
We supposed that pressure is dependent only of J. Using equations 6.4 and 6.5, we deduce that:

p = -κ(J -1) (6.6)
where κ is a constant. We have dp = κdJ, so κ is called bulk modulus for similarity. The hypothesis results in a specific shape for the hydrostatic pressure field. It is coherent with the linear elasticity theory. We note that bulk modulus κ must be independent of volume variation to verify equation 6.4.

Using the hypothesis that pressure is only dependent of the volume, we can obtained a generic expression for the corrected Second Piola Kirchhoff stress tensor:

S kl = Skl + CP P (p -p) 1 2 ∂ p ∂E kl + ∂ p ∂E lk
with CP P = -1 P (p) . From this expression, we deduce that the U/P expression of S depends on:

S = f S, κ, ∂ W ∂E , ∂ 2 W ∂E 2 (6.7)
And that the strain energy density function W must respect in the U/P formulation: p = P ( W ) = -κ(J -1) = -d W dJ (6.8)

Strain energy density function using the reduced invariants

To enforce equation 6.8, we have to isolate the part of the strain energy density that contribute to pressure. For incompressible anisotropic material, the strain energy density function depends on the invariants of deformation W = Wiso (I 1 , I 2 , I 4 , I 5 ), and I 3 = 1. A better assumption is to define: W = Wiso (I 1 , I 2 , I 4 , I 5 ) + Wh (I 3 )

where Wh is a hydrostatic term used to enforced incompressibility. We have:

               P (I 1 ) = -2 3J I 1 P (I 2 ) = -4
3J I 2 P (I 3 ) = -2J P (I 4 ) = -2 3J I 4 P (I 5 ) = -4 3J I 5 (6.9)

We observe from equation 6.9 that invariants of deformation I 1 , I 2 , I 4 and I 5 can contribute to pressure, in addition to the term Wh . This leads to difficulties to implement equation 6.8. The usual trick is to use reduced invariants J 1..5 to separate the terms relative to the hydrostatic dilatation and the ones relative to isochoric deformation. They are defined as: Their influence on the pressure is simplified as ∀i = 3 P (J i ) = 0 and P (J 3 ) = -1. The reduced invariants behave alike the standard invariants, except that they do not influence the pressure. Standard and reduced invariants are equal for totally incompressible material, where I 3 = 1. Using the reduced invariants, we have:      W = Wiso (J 1 , J 2 , J 4 , J 5 ) + Wh (J 3 ) P ( Wiso ) = 0 P ( Wh ) = -κ(J 3 -1) (6.11) with Wiso the isochoric part of the strain energy density function, and Wh the hydrostatic term used to enforced incompressibility, both expressed using reduced invariants. The constraint equation 6.8 is deduced from equation 6.11.

               J 1 = I 1 × I -1/3 3 J 2 = I 2 × I
6.3 Implementation of the hyperelastic material using U/P formulation in the finite element software ADINA

The software ADINA allows the user to define new material laws using Fortran user-coded subroutine.

In the case of a mixed formulation with the Total Lagrangian formulation, Green Lagrange strain tensor is given as input and Second Piola-Kirchhoff stress tensor must be given as output. In addition, the constitutive matrix must be given to complete the stiffness matrix and the U/P formulation. A quick overview of the ADINA subroutine steps is given below. In ADINA, the material subroutine is called for at each element, for each non-linear iteration, for each time steps. Computational efficiency is key for practical uses of the material subroutine.

In the subroutine, the strain and stress tensor are symmetric. The values of each tensor are stored as a vector of size 6, instead of a matrix of size 9. The relation between the vector and the matrix index is:

1 → 11 2 → 22 3 → 33 4 → 12&21 5 → 13&31 6 → 23&32 116
For example, the Green Lagrange strain tensor can thus be written as:

  E 11 E 12 E 13 E 21 E 22 E 23 E 31 E 32 E 33   = (E ij ) ij = E = (E r ) r =         E 1 E 2 E 3 E 4 E 5 E 6         1

2-) Compute first derivatives of invariants with respect to Green Lagrange strain

Using invariants definition, compute: In the U/P formulation, the strain energy function is defined as:

∀k ∈ 1..5, ∂I k ∂E = 1 2
∂I k ∂E ij + ∂I k ∂E ji ij = ∂I k ∂E r r 3
W = W + Q
However, the computation of the corrective term Q is handle automatically by ADINA. Thus it is only needed to compute: W = Wiso (J 1 , J 2 , J 4 , J 5 ) + Wh (J 3 )

Wiso must be any of any kind, but must respect P ( Wiso ) = 0, with P (.) defined in equations 6.3 and 6.5. For equation 6.8 to be respected, we also implement:

Wh = κ 2 (J 3 -1) 2
Thus, we can compute the first derivatives of the strain energy density function based on the displacement field:

∀k ∈ 1..5 ∂ W ∂J k = ∂ Wiso ∂J k + κ(J 3 -1)δ k,3 5

-) Compute Second Piola-Kirchhoff tensor

The corrective term in the strain energy density function is handled by ADINA automatically. The resulting second Piola-Kirchhoff stress tensor, S Q will be added to the displacement part of the second Piola-Kirchhoff stress tensor, S. The displacement part of this stress tensor is defined as:

S = ∂ W ∂E = ∂ W ∂J k ∂J k ∂E = ∂ Wiso ∂J k ∂J k ∂E + κ(J 3 -1) ∂J 3 ∂E Sr r = ∂ Wiso ∂J k ∂J k ∂E r + κ(J 3 -1) ∂J 3 ∂E r r

6-) Compute variables used for mixed U/P formulation correction

The values for the mixed U/P corrections are provided to compute the corrective term of the strain energy density function Q: p = -κ(J 3 -1)

P (p) = κ ∂J 3 ∂E = 1 2J 3 ∂I 3 ∂E = 1 2J 3 ∂I 3 ∂E r r ∂ 2 J 3 ∂E 2 = 1 2J 3 ∂ 2 I 3 ∂E 2 - 1 4J 3 3 ∂I 3 ∂E ⊗ ∂I 3 ∂E = 1 2J 3 ∂ 2 I 3 ∂E r ∂E s - 1 4J 3 3 ∂I 3 ∂E r ∂I 3 ∂E s rs
We note that I 3 is twice differentiable, and so ∂ 2 I 3 ∂Er∂Es is always symmetric. 7-) Compute second derivatives of invariants with respect to Green Lagrange strain ∀k ∈ 1..5,

∂ 2 I k ∂E 2 = 1 2 ∂ ∂E ij ∂I k ∂E r + ∂ ∂E ji ∂I k ∂E r ij r = ∂ 2 I k ∂E r ∂E s rs
We note that ∀k, I k and f k (I 3 ) are twice differentiable, and so their hessian matrix towards deformation are always symmetric.

8-) Compute second derivatives of reduced invariants with respect to Green Lagrange strain Using the decomposition with second order tensor: The corrective part of the constitutive matrix will be added automatically by ADINA. Inside the subroutine, only the the displacement part of the constitutive matrix must be computed. Using the definition of the constitutive matrix and the vector notation for the index we have:

∀k ∈ 1..5 ∂ 2 J k ∂E 2 = f k (I 3 ) ∂ 2 I k ∂E 2 + ∂I k ∂E ⊗ ∂f k (I 3 ) ∂E + ∂f k (I 3 ) ∂E ⊗ ∂I k ∂E + I k ∂ 2 f k (I 3 ) ∂E 2 = f k (I 3 ) ∂ 2 I k ∂E r ∂E s + ∂I k ∂E r ∂f k (I 3 ) ∂E s + ∂I k ∂E s ∂f k (I 3 ) ∂E r + I k ∂ 2 f k (I 3 ) ∂E r ∂E s rs
D = ∂ S ∂E = Drs rs = ∂ Sr ∂E s rs = ∂ 2 W ∂J l ∂J k ∂J k ∂E ⊗ ∂J l ∂E + ∂ W ∂J k ∂ 2 J k ∂E 2 = ∂ 2 W ∂J l ∂J k ∂J k ∂E r ∂J l ∂E s + ∂ W ∂J k ∂ 2 J k ∂E r ∂E s rs
The expressions W = Wiso + Wh and Wh = κ/2(J 3 -1) 2 can be used to obtain:

D = ∂ 2 Wiso ∂J l ∂J k ∂J k ∂E ⊗ ∂J l ∂E + ∂ Wiso ∂J k ∂ 2 J k ∂E 2 + κ(J 3 -1) ∂ 2 J 3 ∂E 2 + κ ∂J 3 ∂E ⊗ ∂J 3 ∂E = ∂ 2 Wiso ∂J l ∂J k ∂J k ∂E r ∂J l ∂E s + ∂ Wiso ∂J k ∂ 2 J k ∂E r ∂E s + κ(J 3 -1) ∂ 2 J 3 ∂E r ∂E s + κ ∂J 3 ∂E r ∂J 3 ∂E s rs
The general form of Wiso is unknown, so the Hessian of Wiso is in general not symmetric. This leads to non-symmetrical constitutive matrix. To impose the symmetry of the constitutive matrix, it is sufficient to require Wiso ∈ C 2 with respect to invariants (Young's theorem). In ADINA, it is possible to use non-symmetric constitutive matrix. The use of non-symmetric matrix will increase the difficulty and the computational cost of the resolution, and will be avoided if possible. Thus, we will also require the condition: Wiso ∈ C 2 in addition to equation 6.11 for the strain energy density function.

Dimensionless analysis of Fluid-Structure coupling

In this section a dimensionless analysis of fluid-structure problem is performed. The fluid, solid and coupling equation are used to obtain information on the coupling strength of the problem studied. This development is obtained from [START_REF] Langre | Fluides et Solides, Ellipses[END_REF]. For a common variable φ belonging to both fluid and solid, φ s refer to solid and φ f to fluid. Vector x will be refered as x and matrix M will be refered as M . Undimensional variable φ for the fluid is referred as φ. Undimensional variable φ for the solid is referred as φ. Ω s is the solid domain, Ω f is the fluid domain and Γ is the interface between Ω s and Ω f . Gravity is neglected.

To simplify the study, isotropic elastic material response is supposed. Thus, the equations governing the structure mechanical response are:

in Ω s : Blood is simplify to a Newtonian and incompressible fluid, so the deviatoric stress tensor can be written as: where n is the normal vector to the surface Γ.

       ρ s ∂ 2 ξ ∂t 2 = div(σ) = (1 + ν) E σ - ν E tr (σ)Id
τ f = 2µd (
Equations (6.13) to (6.15) can be written in a dimensionless form using the fluid characteristic length L f and peak velocity U f . Equation (6.12) can also be written in a dimensionless form using the length L s and the solid wave speed U s = E/ρ s . The corresponding characteristic times are t f = L f U f for the fluid and t s = Ls Us for the solid. To express all the equations in the same temporal reference, we use the dimensionless solid time t s for all the equations. The resulting equations are:

in Ω f : Ls . ξ 0 is the displacement of reference.

List of Figures [START_REF] Lau | Numerical simulation of mitral valve function[END_REF] (blue), [START_REF] Khalighi | Development of a Functionally Equivalent Model of the Mitral Valve Chordae Tendineae Through Topology Optimization[END_REF] (green) and [START_REF] Oliveira | Geometric description for the anatomy of the mitral valve: A review[END_REF] (black).

Figure 1 . 1 :

 11 Figure 1.1: Scheme of the mitral valve and structural components.

Figure 1 . 2 :

 12 Figure 1.2: Ventricular side of porcine mitral valve [2]

Figure 1 . 3 :

 13 Figure 1.3: Measures of the mitral annulus sizes during systole. Int. com. L.: inter-commissural length; Ant. L.: Anteroposterior length; Ant./Pos. H.: anterior/posterior horn height; adapted from[START_REF] Oliveira | Geometric description for the anatomy of the mitral valve: A review[END_REF]

Figure

  Figure 1.4: Noticeable leaflets heights. See table 1.2a for values. H. Ant.: anterior leaflet height; H. P1-3: posterior leaflet height of scallop 1-3; H. PM. com.: posteromedial commissural height; H. AL. Com.: Anterolateral commissural height.

Figure 1 . 5 :

 15 Figure 1.5: Papillary muscle locations from anatomical markers on 57 human mitral valve [20].

Figure 1 . 6 :

 16 Figure 1.6: Papillary muscle locations from recent literature review [14].

Figure 1 . 7 :

 17 Figure 1.7: CT branching types and frequency of observation [26]. f obs = N observed /N valves . The fan-shaped was always observed, mostly for the marginal chordae insertion.

Figure 1 . 8 :

 18 Figure 1.8: Partial reconstruction of the mitral valve from a CT scan with 3D-Slicer. Holes in the leaflet and blood bubbles have the same grey level, leading to a difficult reconstruction of the leaflets.

  (a) Parameters used to construct the annulus (b) D-shape annulus geometric construction (c) Points created for the annulus definition

Figure 1 . 11 :

 111 Figure 1.11: Definition of the plane D-shape mitral valve annulus

Figure 1 .

 1 Figure 1.12: Steps of the method used to smooth the B-splines junction using parabolic interpolation.

Figure 1 . 14 :

 114 Figure 1.14: Construction of the mitral leaflets using translation vectors, B-splines, and a loft.

Figure 1 . 15 :

 115 Figure 1.15: Steps to define the insertions points. 1) Initial spline. 2) First subdivision for the main chordae. 3) Second subdivision for sub-branches. 4) Select closest points 5) Unselect points below tolerance 6) Clear the set of nodes.

Figure 1 . 16 :

 116 Figure 1.16: Steps to define the chord ramification. 1) Initial set of insertion points. 2) Definition of subgroups. 3) Subgroup centre and coupling with papillary muscle. 4) Define the main chord. 5) Define the ramification area. 6) Make sub-branches.

Figure 1 . 17 :

 117 Figure 1.17: Mitral annulus shapes defined using the FreeCAD model. 100 points are used.

Figure 1 . 18 :

 118 Figure 1.18: Definition of mitral leaflets using superquadric function. Axis z is the direction of the superquadric.

Figure 1 .

 1 Figure1.19: Annuli defined using several techniques. Data from[START_REF] Lau | Numerical simulation of mitral valve function[END_REF] (blue),[START_REF] Khalighi | Development of a Functionally Equivalent Model of the Mitral Valve Chordae Tendineae Through Topology Optimization[END_REF] (green) and[START_REF] Oliveira | Geometric description for the anatomy of the mitral valve: A review[END_REF] (black). 100 points were used for each definition.

  (a) First profile (isotropic model from[START_REF] Lau | Numerical simulation of mitral valve function[END_REF]) (b) Second profile (hyperelastic model from[START_REF] Lau | Numerical simulation of mitral valve function[END_REF]) (c) Third profile (µ-CT scan from[START_REF] Sacks | On the Simulation of Mitral Valve Function in Health, Disease, and Treatment[END_REF]

Figure 1 . 20 :

 120 Figure 1.20: Definition of several leaflets using both parametric models. Opaque leaflets are obtained from the FreeCAD only model. Transparent leaflets are obtained from the Fortran & FreeCAD model.

  .23, a & b. A valve from the second parametric model is shown in figure 1.23, c. When compared to µ-CT scans, we observed that

5 Figure 1 . 21 :

 5121 Figure 1.21: Insertion sites use for chord implementation. Parameters are N , the number of principal chord, N sub the number of sub-branches per chord and d sub the distance between each sub-branches

Figure 1 . 22 :

 122 Figure 1.22: Chordae tendineae obtained with automated construction methods.

  (a) Planar D-shape annulus, straight and perpendicular leaflets, no chord ramification (b) Planar D-shape annulus, straight and tilted leaflets, chord ramification (c) Elliptic saddle-shaped annulus, curved leaflets, chord ramification

Figure 1 . 23 :

 123 Figure 1.23: Examples of mitral valves model created with the two parametric models

Figure 1 . 24 :

 124 Figure 1.24: Possible improvements for the chord branching pattern

Figure 2 . 1 :

 21 Figure 2.1: Scheme of the mitral valve main components[START_REF] Kunzelman | Fluid-structure interaction models of the mitral valve: function in normal and pathological states[END_REF]

  (a) Collagen core observed using light microscopy, illumined in the direction of the chord. The waviness of the collagen fibrils can be seen as long as the straight elastic fibril (E) (b) Scheme of the chord internal structure

Figure 2 . 2 :

 22 Figure 2.2: Internal chord structure[START_REF] Millington-Sanders | Structure of chordae tendineae in the left ventricle of the human heart[END_REF]

Figure 2 . 3 :

 23 Figure 2.3: Chordae tendineae mechanical response to uniaxial tests. Values are ± SEM. [52].

Figure 2 . 4 :

 24 Figure 2.4: Several histological views of the mitral valve leaflet. (a) Region and layers of the anterior leaflet microstructure (Movat's pentachrome). Blood flow direction is indicated with the thick arrow. The scale bar is 2 mm [64]. (b) Porcine mitral leaflet composition (Movat's pentachrome). Black is elastin, pink is collagen and blue is ground substances [65].

Figure 2 . 5 :

 25 Figure 2.5: Normalized mitral valve collagen orientation of four anterior and posterior leaflets[START_REF] Einstein | The relationship of normal and abnormal microstructural proliferation to the mitral valve closure sound[END_REF]

Figure 2 . 6 :

 26 Figure 2.6: Membrane stress-strain relation from equibiaxial (open symbol) and strip (filled symbols) of a porcine mitral valve. Circles, circumferential direction; triangles, radial direction [73].

Figure 2 . 8 :

 28 Figure 2.8: Equibiaxial numerical test (ADINA)

Figure 2 . 9 :

 29 Figure 2.9: Stress-strain relationships during an equibiaxial numerical test. Blue and red are for anterior and posterior leaflets, respectively. Markers + and x are experimental data smoothed using the Tikhonov algorithm. Linear elastic model and Mooney-Rivlin parameter model are obtained from the literature. [16], [82].

Figure 2 . 10 :

 210 Figure 2.10: Mitral valve closure at systolic pressure, with several materials. (a) Stress in the principal direction (b) Nodal gap between the two leaflets

Figure 2 . 11 :

 211 Figure 2.11: Unit cube under applied load. R i the load, u i the displacement induced, i = 1..3 [103].

Figure 2 . 12 :

 212 Figure 2.12: Stability criterion applied for Mooney-Rivlin material W = c 1 (I 1 -3) + c 2 (I 2 -3) with c 1 = 1 P a and c 2 = 1P a. 500 different stretch were used for each deformation type.

Figure 2 . 13 :

 213 Figure 2.13: Stability curves of hyperelastic material used. Green uniaxial, blue equibiaxial, red shear.

  (a) Stability curves for anisotropic laws. Green is uniaxial, blue is equibiaxial, red is shear.

  Cauchy cube at rest (1) and under equibiaxial mechanical test (2).

Figure 2 . 14 :

 214 Figure 2.14: Instability under fibre compression. Equibiaxial and uniaxial tests compress the fibres and result in negatives eigenvalues.

Figure 2 . 16 :

 216 Figure 2.16: Formulation used in ADINA with the user coded material subroutine

Figure 2 . 17 :

 217 Figure 2.17: Stress response of the Mofrad & Weinberg model to equibiaxial stretch. Lines are obtained from the subroutine. Circles (respectively squares) are experimental data in the circumferential (respectively radial) directions. Blue (respectively red) corresponds to the anterior (respectively posterior) leaflet.

Figure 2 . 18 :

 218 Figure 2.18: Stress response of the Prot et al. model to uniaxial stretch. Lines are obtained from the subroutine. Circles (respectively squares) are experimental data in the circumferential (respectively radial) directions. Blue (respectively red) corresponds to the anterior (respectively posterior) leaflet.

Figure 2 . 19 :

 219 Figure 2.19: Stress response of the Gao et al. model to equibiaxial stretch. Lines are obtained from the subroutine. Circles (respectively squares) are experimental data in the circumferential (respectively radial) directions. Blue (respectively red) corresponds to the anterior (respectively posterior) leaflet.

  (a) Coaptation modification using Neochord placement or resection and annulus ring[START_REF] Morgan | Neochord Placement Versus Triangular Resection in Mitral Valve Repair: A Finite Element Model[END_REF] (b) Coaptation of the mitral valve at systole after infarctus without (left) and with (right) an annuloplasty ring[START_REF] Baillargeon | Human Cardiac Function Simulator for the Optimal Design of a Novel Annuloplasty Ring with a Sub-Valvular Element for Correction of Ischemic Mitral Regurgitation[END_REF]

Figure 3 . 1 :

 31 Figure 3.1: Coaptation modification observed using structural simulation.

Figure 3 . 2 :

 32 Figure 3.2: Stress values obtained in structural simulation of the mitral valve

  (a) Mitral valve meshed with Belytschko-Lin-Tsay shell elements in LS-DYNA[START_REF] Lau | Numerical simulation of mitral valve function[END_REF] (b) Mitral valve meshed with 3D tetrahedral elements in IMPETUS[START_REF] Lau | Numerical simulation of mitral valve function[END_REF]

Figure 3 . 3 :

 33 Figure 3.3: Mitral valve meshed with shell and 3D elements

1 :

 1 Characteristic values of the structural problem

Figure 3 . 4 :

 34 Figure 3.4: Configuration of the normal response to gap in ADINA (CAT 1 contact). λ g the normal force, g the gap, p g the normal pressure and A the contactor area.

Figure 3 . 5 : 2 F

 352 Figure 3.5: Surfaces representation when the offset is used

Figure 3 . 6 :

 36 Figure 3.6: 3D shell control vectors in undeformed (a) and deformed (b) configuration

Figure 3 . 7 :

 37 Figure 3.7: Damping amplitude ξ n function of pulsation ω n obtained with Rayleigh damping method.

Figure 3 . 8 :

 38 Figure 3.8: CPU time function of Rayleigh parameter β for different values of contact damping C damp .

Figure 3 . 9 :

 39 Figure 3.9: Boundary condition, mesh and pressure variation of the shell model

Figure 3 . 10 :

 310 Figure 3.10: Principal stress on the mitral valve with 3D shell elements (top/bottom shell)

Figure 3 . 12 :

 312 Figure 3.12: Closure of the mitral valve (systole to early diastole). Left and middle: Smoothed stress in the principal direction. Right: isosurface of coaptation of contact (coaptation surface)

Figure 3 . 13 :

 313 Figure 3.13: Smoothed principal stress of the valve leaflets in closed configuration

Figure 3 . 14 :

 314 Figure 3.14: Mesh, boundary condition and pressure variation of the 3D-element model

Figure 3 . 15 :

 315 Figure 3.15: Smoothed principal stress of the closed 3D element mitral valve.

Figure 3 . 16 :

 316 Figure 3.16: Nodal contact gap of the closed 3D element mitral valve.

Figure 3 . 17 :

 317 Figure 3.17: Fiber stretch, Principal stretch and Principal stress from the 3D-elements mitral valve with Gao material law during diastole.

Figure 3 . 18 :

 318 Figure 3.18: Fiber stretch, Principal stretch and Principal stress from the 3D-elements mitral valve with Gao material law at early systole.

Figure 4 . 1 :

 41 Figure 4.1:Mitral regurgitation jet at late systole imaged using colour Dopler[START_REF] Grayburn | Quantitation of mitral regurgitation[END_REF]

Figure 4 . 2 :

 42 Figure 4.2: Weak coupling of fluid and structure solvers.

Figure 4 . 3 :

 43 Figure 4.3: Strong coupling of fluid and structure solvers.

Figure 4 . 4 :

 44 Figure 4.4: Monolithic coupling or direct approach for fluid-structure problems.

Figure 4 . 5 :

 45 Figure 4.5: Several formulations applied to the same domain. Ω s represent the solid domain, Ω f is the fluid domain, and Γ is the interface.

Figure 4 . 6 :

 46 Figure 4.6: Remeshing of the fluid domain after large displacement

Figure 4 . 7 :

 47 Figure 4.7: Example of mesh obtained after using the SAM method. Two criterion are used, one based on the gradient of Eulerian variables and the second based on the belonging to the orange squares.

(a) t 1 (b) t 2 (c) t 3

 123

Figure 4 . 8 :

 48 Figure 4.8: Fluid mesh reconstruction function of mesh deformation obtained in ADINA

Figure 4 .

 4 [START_REF] Al-Atabi | Biomechanical assessment of surgical repair of the mitral valve[END_REF] gives an example of an automatic remeshing technique. At t 1

Figure 4 . 9 :

 49 Figure 4.9: Calls of the Steered Adaptive Meshing (SAM) procedure. Purple is repair only, green is at specified timesteps.

Figure 4 .

 4 Figure 4.10: 2D solid and fluid elements and nodes.

  To define the spatial and temporal resolution, we choose to monitor the conservation of energy. The formulation defined in section 4.2.1 is used. Simplified fluid and fluid-structure models were used to test the method (figure 4.11). The fluid-structure model is a straight and rigid rectangle with half a disk laid on top of the domain. The boundary of the domain is rigid. The structure part is a moving leaflet, of thickness T and length L, fixed to the border of the half-disk. A velocity is imposed as input using an established Poiseuille profile and the velocity flow displayed in figure 4.11a. Output is free outflow. Flow is laminar with incompressible material of density ρ f and viscosity µ. The solid material is Hookean in large deformation of Young's modulus E, density ρ s and Poisson's ratio ν. Values for the problem are given in table 4.2. The fluid model is the same as the fluid-structure model, except for the moving leaflet of length L that was taken out.

Figure 4 . 11 : 6 ×

 4116 Figure 4.11: Simplified model used to define spatial and temporal resolution

Figure 4 . 12 :

 412 Figure 4.12: Example of energetic balance on uniform mesh

  (a) Error function of mesh size (b) CPU time function of mesh size

Figure 4 . 13 :

 413 Figure 4.13: Error and CPU time function of discretization used for the fluid model

  (a) Error obtained function of number of remeshing during the simulation (b) CPU time function of number of remeshing during the simulation

Figure 4 . 14 :

 414 Figure 4.14: Effect of the number of remeshing on the fluid model. Remeshing based on vorticity and pressure gradient with 0.2 cm > ∆x > 0.025 cm

Figure 4 . 15 :

 415 Figure 4.15: Evolution of remeshing frequency function of the simulation time. Only repair is used. Criterion is pressure and vorticity gradient with 0.5 cm > ∆x > 0.05 cm.

  (a) Fluid mesh at t = 0.19 s (b) Fluid mesh at t = 0.35 s (c) Fluid mesh at t = 0.53 s (d) Fluid mesh at t = 0.84 s

Figure 4 . 16 :

 416 Figure 4.16: Fluid mesh of the fluid-structure simulation at the fifth pulsation. Boundary layers and vortex grid mesh are directly linked to the vorticity gradient condition ∇Ω x .∆ x const, with ∆ x the spatial discretisation and Ω x the vorticity field. This condition leads to 0.025 cm ≤ ∆ x ≤ 0.2 cm.

Figure 4 . 17 :

 417 Figure 4.17: Implementation of the GAP function based on a time criteria

  (a) Pressure, t = t GAP -∆t (b) Pressure, t = t GAP (c) Velocity, t = t GAP -∆t (d) Velocity, t = t GAP

Figure 4 . 18 :

 418 Figure 4.18: Pressure and velocity field before and after the GAP activation. Pressure ranges from 0 (purple) to 180 Pa (red). Velocity ranges from 0 (purple) to 50 cm.s -1 (red).

Figure 4 . 19 :

 419 Figure 4.19: Boundary condition of the second test case

  (a) Mesh, t = t GAP (b) Mesh, t = t GAP + ∆t (c) Mesh, t = t GAP + 4∆t (d) Pressure, t = t GAP (e) Pressure, t = t GAP + ∆t (f) Pressure, t = t GAP + 4∆t(g) Velocity, t = t GAP (h) Velocity, t = t GAP + ∆t (i) Velocity, t = t GAP + 4∆t

Figure 4 . 20 :

 420 Figure 4.20: Results of the second case with a GAP condition based on the distance with the wall. Pressure ranges from 0 (purple) to 180 Pa (red). Velocity ranges from 0 (purple) to 50 cm.s -1 (red).

  (a) 2D FCBI elements (b) Forcing no flow condition with FCBI elements

Figure 4 .

 4 Figure 4.21: "Pseudo" GAP condition using FCBI elements and the no-slip boundary condition

  (a) Ventricle model (b) Mesh at t = 0.917 s (c) Mesh at t = 1.987 s

Figure 4 . 22 :

 422 Figure 4.22: Domain of the simplified 2D ventricle model and examples of meshes used. Pressure is imposed at the atrium and ventricle border (dashed in (a)).

  (a) Domain used and boundary condition (b) Normalized imposed pressure and time postprocessed

Figure 4 . 25 :

 425 Figure 4.25: Domain and boundary condition of the 2D left heart. Dashed boundary are used to impose pressure.

Figure 4 . 26 :

 426 Figure 4.26: Pressure and average velocity variation during a cardiac cycle (4-th pulsation).

  (a) t = 0.20 s (b) t = 0.30 s (c) t = 0.32 s (d) t = 0.35 s (e) t = 0.50 s (f) t = 0.65 s (g) t = 0.70 s (h) t = 0.90 s

Figure 4 . 27 :

 427 Figure 4.27: Mesh evolution during a cardiac cycle

Figure 4 . 28 :

 428 Figure 4.28: Velocity magnitude evolution during a cardiac cycle. Red is 60 cm.s -1 and purple is 0 cm.s -1 .

  (a) t = 0.20 s (b) t = 0.30 s (c) t = 0.32 s (d) t = 0.35 s (e) t = 0.50 s (f) t = 0.65 s (g) t = 0.70 s (h) t = 0.90 s

Figure 4 . 29 :

 429 Figure 4.29: Nodal pressure evolution during a cardiac cycle. Red is 16 KPa (120 mmHg) and purple is 0 Pa.

Figure 4 . 30 :

 430 Figure 4.30: Cardiac output of the simplified 2D left ventricle

Figure 6 . 1 :

 61 Figure 6.1: Deformation of a continuum body Ω between time 0 and t.

Figure 6 . 2 :

 62 Figure 6.2: Decomposition of the deformation using the polar theorem (Truesdell & Noll, 1965).

  -) Compute Right Cauchy Strain tensor C and its invariants Implement C = 2E+I. From C and N , the fiber in the undeformed configuration, define invariants I 1 , I 2 , I 3 , I 4 and I 5 .

3 )∂f k (I 3 )

 33 -) Compute reduced invariants and their first derivatives with respect to Green Lagrange strainWe can write: ∀k ∈ 1..5,J k = I k × f k (IThus the derivatives of the reduced invariants are defined from the decomposition rule:∂E ij + ∂f k (I 3 )∂E ji ij 4-) Compute Strain Energy Density Function and its first derivative with respect to reduced invariants

with ∂ 2

 2 f k (I 3 )

2 ∇

 2 t ξ + ∇ξ is the strain tensor, σ the Cauchy stress, ν Poisson's ratio, E Young's modulus, ρ the density, ξ the displacement.Navier-Stokes simplified equation for incompressible fluids are used:in Ω f : u.∇) (), τ f the deviatoric (or viscous) stress tensor, u the velocity vector and p the pressure.

6 . 14 ) where d = 1 2 ∇

 6142 t u + ∇u is the deformation tensor and µ is the dynamic viscosity. The kinetic and dynamic boundary conditions implemented at the interface between the fluid and solid domains are: + τ f ).n = σ.n(6.15)

  Us is the reduced velocity, Re =ρ f U f L f µ the Reynolds number, M = ρ fρs the mass number and D = ξ 0

  Com.: Anterolateral commissural height. . . . . . . . . . . . . . . . . . 1.5 Papillary muscle locations from anatomical markers on 57 human mitral valve [20]. . . 1.6 Papillary muscle locations from recent literature review [14]. . . . . . . . . . . . . . . . 1.7 CT branching types and frequency of observation [26]. f obs = N observed /N valves . The fan-shaped was always observed, mostly for the marginal chordae insertion. . . . . . . 1.8 Partial reconstruction of the mitral valve from a CT scan with 3D-Slicer. Holes in the leaflet and blood bubbles have the same grey level, leading to a difficult reconstruction of the leaflets. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.9 Imaging techniques used to model the mitral valve . . . . . . . . . . . . . . . . . . . . 1.10 Geometric models of the mitral valve . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.11 Definition of the plane D-shape mitral valve annulus . . . . . . . . . . . . . . . . . . . 1.12 Steps of the method used to smooth the B-splines junction using parabolic interpolation. 1.13 Example of leaflet height function of annular length. Data from [16] . . . . . . . . . . 1.14 Construction of the mitral leaflets using translation vectors, B-splines, and a loft. . . . 1.15 Steps to define the insertions points. 1) Initial spline. 2) First subdivision for the main chordae. 3) Second subdivision for sub-branches. 4) Select closest points 5) Unselect points below tolerance 6) Clear the set of nodes. . . . . . . . . . . . . . . . . . . . . . 1.16 Steps to define the chord ramification. 1) Initial set of insertion points. 2) Definition of subgroups. 3) Subgroup centre and coupling with papillary muscle. 4) Define the main chord. 5) Define the ramification area. 6) Make sub-branches. . . . . . . . . . . . . . . 1.17 Mitral annulus shapes defined using the FreeCAD model. 100 points are used. . . . . . 1.18 Definition of mitral leaflets using superquadric function. Axis z is the direction of the superquadric. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.19 Annuli defined using several techniques. Data from

  100 points were used for each definition. . . . . . . . . . . . . . . . . . . . . . 1.20 Definition of several leaflets using both parametric models. Opaque leaflets are obtained from the FreeCAD only model. Transparent leaflets are obtained from the Fortran & FreeCAD model. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.21 Insertion sites use for chord implementation. Parameters are N , the number of principal chord, N sub the number of sub-branches per chord and d sub the distance between each sub-branches . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  1.4: Noticeable leaflets heights. See table 1.2a for values. H. Ant.: anterior leaflet height; H. P1-3: posterior leaflet height of scallop 1-3; H. PM. com.: posteromedial commissural height; H. AL. Com.: Anterolateral commissural height.

	H. Ant.	H. P1	H. P2	H. P3	H. AL.	H. PM.	source
					com.	com.
	21.2 ± 3.0	5 -17	9.8 ± 2.0	5 -17	5 -13	4 -12	[14]
	23.4 ± 2.9	11.0 ± 2.3	13.8 ± 2.9	11.0 ± 2.3	8.7 ± 2.4	8.3 ± 2.0	[20]
	24	11	14	11	0	0	[17]
	23.3 -20.2

Table 1 .

 1

		13]	[30]	[27]	[31]
	Marginal	0.08	0.38 (±0.18)
	Basal	0.19	0.71 (±0.25)	0.50	0.21(±0.15)
	Strut	0.39	2.05 (±0.40)

4: Cross Sectional Area of chordae tendineae (CSA (± SD)) in mm 2 obtained in the literature.

Table 2 .

 2

1: Structural stiffness evolution between young (18-26 weeks) and old (over 2 years) porcine mitral valve

[START_REF] Millard | Mechanical properties of chordae tendineae of the mitral heart valve: Youngs modulus, structural stiffness, and effects of aging[END_REF]

Table 2 .

 2

	MPa)	a 3

2: Material parameters for the strain energy density function defined in

table 2 .

 2 3).

	Leaflet	Direction	E pre kPa	E post kPa	R
	Anterior	Circ.	89.3 ± 91.7	8960 ± 8830	5.71 ± 4.38
		Rad.	60.5 ± 43.8	2400 ± 1480
	Posterior	Circ.	93.8 ± 52.1	3430 ± 2090	4.27 ± 4.15
		Rad.	42.3 ± 23.9	1920 ± 1310
	Table

Table 2 .

 2 5: Coefficients of the Prot model fitted from experimental data[START_REF] Prot | On modelling and analysis of healthy and pathological human mitral valves: Two case studies[END_REF]

	.5

Table 2 .

 2

	kPa)	b f

6: Coefficients for the model of Gao and co-workers

Table 2 .

 2

	kPa)

7: Material value defined for the fibre splay model

[START_REF] Einstein | Haemodynamic determinants of the mitral valve closure sound: A finite element study[END_REF]

Table 3 .

 3 ∆t iniThe parameters used for the simulation of the mitral valve are given in table 3.2. Thus, the time-step varies between 10 -4 s and 10 -6 s.

	If there is no convergence then
		∆t =	∆t prev DIVFAC
	If there are NCONV consecutive convergences then
		∆t = MULFACT × ∆t prev
	Else
		∆t = ∆t prev
	DIVFAC	MAXSUB	MULFACT	NCONV
	2	100	2	5

2: Parameters used for the ATS method.

Table 4 .

 4 2: Values used for the simplified valve model

Table 4 .

 4

	3: Physical values of the 2D mitral valve model
	4.3.2 Results at the boundary

Pressure and average velocity at the boundaries are given in figure

  .28. Initially, the blood flows from the atrium to the ventricle. The inflow velocity reaches 1 m.s -1 . The valve quickly closed when the pressure increased to the systolic value. During closure (from 0.3 s to 0.35 s), part of the ventricle blood is ejected toward the atrium. This regurgitant jet is a result of the leaflet movement and can not be avoided. The closing leaflets help propel the jet toward the atrium (figure4.28,c & d). The regurgitant jet has an average velocity of 50 cm.s -1 at the atrium border and a local maximal value of 3 m.s -1 (figure 4.26). Small vortices appear in the atrium as a result of this regurgitant jet. In the ventricle, a large vortex appears as the valve closes and seems to help maintain the valve in a closed state. The intensity of the vortex decreases over time. Apart from the initial jet, no leaking is observed when the valve is closed (figure 4.28,e). When systole ends, the leaflets move away from each

  1.1 Scheme of the mitral valve and structural components. . . . . . . . . . . . . . . . . . . 1.2 Ventricular side of porcine mitral valve [2] . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Measures of the mitral annulus sizes during systole. Int. com. L.: inter-commissural length; Ant. L.: Anteroposterior length; Ant./Pos. H.: anterior/posterior horn height; adapted from [14] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4 Noticeable leaflets heights. See table 1.2a for values. H. Ant.: anterior leaflet height; H. P1-3: posterior leaflet height of scallop 1-3; H. PM. com.: posteromedial commissural height; H. AL.