Pablo Rey

Pierre Antolin

Catherine Azzoro-Pantel

Rachel Calvet

Kamila Koziura

Mustapha Julie Hugo

Amandine, David Yuri Gabriel

Chère lectrice, cher lecteur, Ce manuscrit représente les travaux effectués pendant ces trois dernières années passées à Toulouse. Il ne serait pas juste de qualifier ces années de faciles, surtout compte tenu du contexte sanitaire et géopolitique mondial. Néanmoins, cela a été une expérience enrichissante, et je suis heureuse de pouvoir arriver à la fin de ma thèse. Il est maintenant temps d'exprimer ma gratitude envers toutes les personnes qui m'ont accompagnée. Tout d'abord, je remercie mes directeurs de thèse, Robin et Paul, pour leur engagement, leurs conseils scientifiques et leur disponibilité constante. La chance m'a souri en vous ayant comme guides durant cette aventure académique ! Je tiens à exprimer ma gratitude envers David Dureisseix et Gregory Legrain d'avoir accepté d'être les rapporteurs de cette thèse. Je remercie également les examinateur•rice•s, Valentine Rey, Pablo Antolin, Pierre Gosselet, et Eric Florentin, pour avoir consacré du temps à ma soutenance.

Mes remerciements s'étendent également à l'école doctorale MEGEP, représentée par Catherine Azzoro-Pantel, Rachel Calvet, et Kamila Koziura. Leur présence constante, leur capacité à répondre à mes questions et à fournir toutes les explications nécessaires ont été précieuses.

L'idée de faire une thèse trouve ces racines dans ma formation de Master. C'est à ce temps-la que j'ai découvert la beauté de la recherche. Je remercie Alexander Fredin : même si notre collaboration n'a pas duré très longtemps, elle a sûrement laissé des traces dans mon âme.

Tout au long des années de thèse, j'étais entourée par de nombreux•ses collègues au sein du GMM et de l'ICA. Je voudrais mentionner chacun•e, bien que je risque d'en oublier certain•e•s. Je commence par la jeunesse: les doctorant•e•s et les ATER, qu'iel•le•s soient ancienn•e•s ou actuel•le•s. Merci Ali, j'espère tu es en train de surfer et de te profiter des meilleurs vagues possibles. Clément, est-ce que tu écoutes toujours les marches ? Alban, je te souhaite tout le succès dans tes nombreux projets. Anthelme et Ousmane, merci pour les discussions enrichissantes dans le bureau. Hippolyte, tes connaissances en géographie sont vraiment fascinantes ! Julie, Hugo, Mustapha, Amandine, David, Yuri, Gabriel, Dorian, je vous souhaite une bonne continuation. Mouhamad, chaque conversation avec toi est un plaisir, et tes enfants sont vraiment chanceux•ses de t'avoir comme papa. Mes remerciements s'étendent à tous les autres membres du GMM et de l'ICA. Cathy, merci pour ta co-gestion du GMM et ta gentillesse envers moi. Sandrine, merci pour ton bon humeur, toutes les procédures bureaucratiques semblent faciles avec toi. Pierre, merci pour tes conseils scientifiques en mécanique du contact qui m'ont aidée dans mes travaux de thèse. Merci Jean-Charles pour ton intérêt pour mon travail, tes questions ont parfois été difficiles à répondre, mais toujours intéressantes à explorer. Romain, je te souhaite une soutenance de HDR réussie, et j'espère que le gâteau aux épinards t'a plu ! Merci Frédéric, vraiment, si tu n'existais pas, il faudrait t'inventer ! Merci Olivier pour tes connaissances et ta capacité à poser des questions pertinentes qui grandement enrichissaient les discussions du midi. Enfin, Ion, i Dearest gentle reader, This manuscript represents the work conducted over the past three years that I spent in Toulouse. It would not be fair to characterize these years as easy, especially considering the global sanitary and geopolitical context. Nevertheless, it has been a rewarding experience, and I am delighted to reach the end of this work. It is now time to express my gratitude to all the individuals who have accompanied me.

First and foremost, I thank my thesis supervisors, Robin and Paul, for their commitment, scientific guidance, and constant availability. I consider

ii merci pour l'organisation des randonnées : changer le bruit de la ville pour un paysage montagneux fait toujours du bien.

Il est bien connu que la meilleure façon d'apprendre est d'enseigner. En effet, j'ai pu redécouvrir les bases des mathématiques ainsi que la théorie de l'élasticité. Modestement, j'espère que certain•e•s parmi mes élèves ont pu percevoir l'intérêt et l'importance de ces matières. Ici, je remercie Lorick et Anthony, avec qui on s'occupait des cours de maths 2 pour de petits anglophones. Je remercie également toute l'équipe de l'enseignement de la MMC, Eduard, Thomas, Paul, Stéphane, pour leur soutien. Je vous souhaite avoir des étudiants motivés.

Comme vous pouvez le constater, l'intégralité du manuscrit est rédigée en anglais, et pour cela, je remercie Anastasiia Zakharova. C'est avec une chaleur particulière que je me remémore nos cours d'anglais, qui souvent dépassaient le cadre de l'apprentissage de la langue avec des discussions animées et une ambiance amicale.

J'approche tranquillement de la conclusion de ces remerciements. Dans cette dernière section, je tiens à exprimer ma reconnaissance envers celles et ceux qui, dans la vie quotidienne, au-delà du travail, m'ont apporté leur soutien. Je dis merci à ma mère, même si je ne suis pas sûre qu'elle lise ces mots (fais-moi signe si un jour tu y arrives !). C'est grâce à tes sacrifices, ta force pour aller à contre-courant du milieu, et ton amour que je mène aujourd'hui une vie dont je ne pouvais même pas rêver dans mon enfance. Mon cher Thibault, merci pour ton patience, ton amour, ton soutien. Peu importe la langue utilisée, je ne connais pas suffisamment de mots pour exprimer ce que tu es pour moi. Merci également à toute ma belle-famille : Agnès, Alban, Adeline, Camille, Philip. Vous êtes merveilleux•euses, je suis très heureuse de vous connaître. On a partagé de bons moment ensemble et je suis sûre que l'on en partagera encore beaucoup. Enfin, le membre à quatre pattes de la famille : merci, Fortune, d'être toujours d'accord pour accepter les caresses ! Merci à tous•tes! Evgeniia iv It is well known that the best way to learn is to teach. Indeed, I was able to rediscover the fundamentals of mathematics and the theory of elasticity. Modestly, I hope that some of my students have been able to perceive the interest and importance of these subjects. Here, I express my gratitude to Lorick and Anthony, with whom I worked on the Math 2 courses for young English speakers. I also extend my thanks to the entire team of the MMC teaching staff -Eduard, Thomas, Paul, Stéphane -for their support. I wish you to have motivated students.

As you can see, this manuscript is written in English, and for this, I thank Anastasiia Zakharova. It is with a special warmth that I recall our English classes, which often went beyond the scope of language learning with lively discussions and a friendly atmosphere.

I am slowly approaching the end of these acknowledgements. In this final section, I want to express my gratitude to all those who, in everyday life, beyond work, have supported me. I thank my mother, even though I'm not sure if she read these words (let me know if you ever do!). It is thanks to your sacrifices, your strength to go against the current, and your love that I have a life today that I could not even dream of in my childhood. My dear Thibault, thank you for your patience, your love, your support. Regardless of the language used, I don't know enough words to express what you mean to me. Also, thanks to my entire French family: Agnès, Alban, Adeline, Camille, Philip. You are wonderful, and I am very happy to know you. We have shared good moments together, and I am sure we will share many more. Finally, the four-legged member of the family: thank you, Fortune, for always agreeing to receive cuddles! Thank you all!

Evgeniia

Resumé

Avec le développement des techniques d'imagerie volumique, il est désormais possible de construire des modèles mécaniques d'un matériau à l'échelle de la microstructure à partir d'images 3D d'un échantillon. Cependant, le calcul de ces modèles pour des applications réelles reste un défi. L'objectif de ce travail est de présenter une approche numérique pour simuler les matériaux composites à l'échelle de leur microstructure.

Plusieurs défis sont associés à cette tâche. Les modèles sont associés à des données massives, ils peuvent avoir une géométrie complexe et faire appel à des lois de comportement non-linéaires. Pour répondre à ces problèmes, le présent travail propose d'utiliser l'Analyse IsoGéométrique. Cette approche offre une plus grande précision par degré de liberté et une plus grande robustesse que la méthode des éléments finis pour les simulations mécaniques. Cela implique qu'elle est souvent considérée comme un outil de calcul à haute performance. Cependant, ces meilleures propriétés sont limitées au cas où la solution est régulière, ce qui rend l'utilisation de l'Analyse IsoGéométrique discutable pour modéliser les non-linéarités locales. Dans ce contexte, deux stratégies distinctes sont proposées. Premièrement, l'Analyse IsoGéométrique est appliquée à l'échelle globale et couplée de manière non invasive à un modèle local d'éléments finis. Cette approche permet de modéliser divers comportements non linéaires présents dans le modèle local, car la méthode des éléments finis est mieux adaptée à la modélisation de comportements fortement non linéaires et propose une grande variété de procédures déjà mises en oeuvre. Bien qu'initialement proposée pour les matériaux composites, cette stratégie s'est révélée intéressante dans d'autres contextes, tels que la modélisation d'assemblages. Deuxièmement, un cas particulier de modélisation de microstructures hétérogènes avec de nombreuses inclusions linéaires reliées à la matrice linéaire par diverses interfaces non linéaires a été examiné. Pour ce cas spécifique, une méthode de décomposition de domaine mixte appelée méthode LaTIn, combinée avec des techniques de frontière immergée, a été utilisée dans le cadre de l'Analyse IsoGéométrique. Différents exemples numériques montrent les performances et la polyvalence des méthodologies proposées. 1: computed tomography tomography scan of a sample of composite material [START_REF] Wang | Consistent immersed volumetric Nitsche methods for composite analysis[END_REF].

or changes that might not be visible from the outside. Nowadays, these techniques are enough developed to provide highly accurate 3D images that take into account the microstructure of a sample. An example of a 3D image of a composite sample obtained by computed tomography tomography is shown on Figure 1 .

Using such 3D images, one can build mechanical models to compute material's behaviour at the micro-scale (see Figure 2 for an example of such a model). Such models can be used for homogenization analysis [START_REF] Legrain | An X-FEM and level set computational approach for image-based modelling: application to homogenization[END_REF], patient-specific simulations in the biomedical domain [START_REF] Claus | Contact modeling from images using cut finite element solvers[END_REF]. Another application of micromechanical simulation of complex materials lies in the emerging fields of data assimilation and digital twins. Combining the methodology with in-situ experimental tests opens the door for data assimilation at the micro scale [START_REF] Heinze | Experimental and numerical investigation of single pores for identification of effective metal foams properties[END_REF][START_REF] Wu | A combined experimental and numerical method to estimate the elastic modulus of single trabeculae[END_REF], by comparing measured fields coming, e.g., from digital image correlation [START_REF] Bibliography | Digital volume correlation: review of progress and challenges[END_REF][START_REF] Rouwane | Architecturedriven digital volume correlation: application to the analysis of in-situ crushing of a polyurethane foam[END_REF], and simulated fields obtained from computational micromechanics. In other contexts, computational micromechanics offers the opportunity to feed artificial intelligence to create homogenized constitutive laws (Fernández et al., 2021, Masi and[START_REF] Masi | Multiscale modeling of inelastic materials with thermodynamics-based artificial neural networks (tann)[END_REF], or to perform virtual material design [START_REF] Herráez | Computational micromechanics evaluation of the effect of fibre shape on the transverse strength of unidirectional composites: an approach to virtual materials design[END_REF], i.e. to optimize the material microstructure to achieve a desired behaviour.

However, the computation of image-based models seems still out of reach for real mechanical applications. Indeed, this is obstructed due to several reasons. First, these models are associated with massive data since one image can contain 2000×2000×2000 voxels that results in 2000 3 (eight billion) elements [START_REF] Liu | An efficient strategy for large scale 3d simulation of heterogeneous materials to predict effective thermal conductivity[END_REF], Neggers et al., 2018). Thus, a great number of degrees of freedom (DoF) is needed to represent correctly the complicated geometry (see again Figure 1). Secondly, they appeal to complex laws. Indeed, since different subdomains are present (i.e. matrix and fibres), special attention should be paid to modelling the interface between them: there can be, for example, multiple contact zones [START_REF] Claus | A stable and optimally convergent LaTIn-CutFEM algorithm for multiple unilateral contact problems[END_REF], or delamination (Alfano andCrisfield, 2001, Nguyen et al., 2014). Non-linearities can exist not only at the interfaces but also within the subdomains themselves. For example, the matrix material can be elasto-plastic, or fibre may exhibit buckling [START_REF] Qiu | An interface element formulation for the simulation of delamination with buckling[END_REF][START_REF] Saavedra | On a multiscale strategy and its optimization for the simulation of combined delamination and buckling[END_REF]. Finally, there exist a great difference in material properties of the matrix and the fibres. Therefore, significant resources and adapted complex numerical modelling are required to perform numerical simulations.

In line with this background, the objective of this work is to develop an efficient 5

Figure 2: CAD image-based model of fibres [START_REF] Chelaghma | Fonctionnalisation de composites C/PEKK pour application aérospatiale: caractérisation, modélisation et influence sur les propriétés du composite[END_REF].

and accurate numerical method able to compute geometrically complex composite microstructures involving numerous inclusions connected with the matrix and taking into account various non-linear behaviours.

Recent works propose different advanced numerical techniques, such as combining fictitious domain type approaches with shape functions with higher degree and regularity [START_REF] Lian | Image-based computational homogenization and localization: comparison between x-fem/levelset and voxel-based approaches[END_REF][START_REF] Verhoosel | Imagebased goal-oriented adaptive isogeometric analysis with application to the micromechanical modeling of trabecular bone[END_REF][START_REF] Wang | Consistent immersed volumetric Nitsche methods for composite analysis[END_REF], using unfitted Finite Element Analysis with mixed domain decomposition methods (Claus andKerfriden, 2018, Kerfriden et al., 2020), multiscale modelling [START_REF] Saavedra | On a multiscale strategy and its optimization for the simulation of combined delamination and buckling[END_REF], mixed-dimensional coupling [START_REF] Kerfriden | A mixed-dimensional CutFEM methodology for the simulation of fibre-reinforced composites[END_REF], and hybrid isogeometric-finite element approach [START_REF] Maleki-Jebeli | A large deformation hybrid isogeometric-finite element method applied to cohesive interface contact/debonding[END_REF], to name a few.

In the context of this thesis, our focus is on the computational step, assuming that the image processing stage has been completed and we now have a geometric model. This work proposes two strategies to deal with image-based mechanical models of fibre-reinforced composite materials. First, we propose to harness the best of both FEA and IsoGeometric Analysis (IGA) in global/local coupling. This strategies tackles the problem of modelling local non-linearities (such as local contact, local damage, local plasticity) in a global linear structure. FEA is utilized for the simulation of the local possibly non-linear behaviour since its low regular discretization space appears more adapted to capture local, strongly non-linear or even singular phenomena (e.g., contact, delamination, heterogeneities). With FEA on a local level, non-linear phenomena presented within the subdomain (such as non-linear material behaviour) can be efficiently modelled. The global solution is captured by IGA, thus benefiting from its superior behaviour for regular solutions. The coupling scheme is fully non-invasive: the initial global spline model is never modified and the construction of the coupling operators can be performed using conventional FE packages. Secondly, for a particular case of non-linear interfaces in a globally linear structure, a fully IsoGeometric approach is suggested. Whilst meshing is alleviated by immersed boundary-conformal approaches, the non-linear interface behaviour is dealt with mixed domain decomposition method. This approach is naturally parallelizable and flexible to treat any non-linear interface behaviour. Thus, the first approach is more general, since it can be used to deal with any kind of non-linearities, while the second approach is tailored for specific case of non-linear interfaces in globally linear structure.

Outline

In order to present this work, the manuscript is sub-structured into three chapters:

INTRODUCTION

• The first chapter presents a literature review on the problem of multi-scale modelling of scan-based mechanical models. The first part presents IGA from the viewpoint of FEA. The relation between these two methods is discussed, and a FEA-to-IGA bridge is presented. Immersed boundary-conformal approaches are considered in the second part as well as different methods for weak coupling.

The third part addresses the issue of global/local coupling. Monolithic and noninvasive iterative solutions are reviewed. Finally, the fourth part presents a mixed domain decomposition method called LaTIn that is used for the case when the structure is globally linear except for interfaces. Here, this method is presented for a general case, without the precision of the nature of the interface behaviour.

• The second chapter is dedicated to the non-invasive coupling between the global IGA model and local FEA models. The strategy proposed relies on the FEAto-IGA bridge, introduced in Chapter 1 which also enables to arrive at a fully non-invasive strategy in the sense that not only the global/local coupling is noninvasive but also the construction of the coupling operators from only FE resources. A series of numerical experiments that cover 2D and 3D simulations with different local behaviours, such as cracks, contact and delamination are presented using only FEA software, thereby going further than the sole application of composite materials.

• The third chapter addresses the specific case when non-linear behaviour is encountered only on the interfaces. For this case, a fully IsoGeometric practical approach is proposed. A mixed domain decomposition method called LaTIn is used to model interface behaviour, while meshing is alleviated with the use of immersed boundary techniques. The use of the LaTIn method allows to efficiently model different non-linear interface behaviours and also makes possible to perform the computation in parallel. The chapter provides implementation details on different interface behaviours (perfect interface, frictionless contact, Coulomb friction, delamination). The robustness of the approach is illustrated throughout a series of 2D examples.

After that, a general conclusion and perspectives for future work are discussed.

Contributions

This thesis resulted in the following publications and conference talks.

• Articles in scientific journals 1. Lapina, E., [START_REF] Lapina | A fully non-invasive hybrid IGA/FEM scheme for the analysis of localized non-linear phenomena[END_REF] This chapter aims to provide an overview of advanced simulation strategies for modelling complex structures. It sequentially covers the following topics: IsoGeometric Analysis in relation to Finite Element Analysis, Immersed Boundary approach as a viable solution to model geometrically complex structures, Global/Local coupling as a means to model localized non-linear behaviours, and the Large Time Increment method, which is tailored to scenarios where non-linearities occur only at interfaces. et al. (2005) and seeks to unify the fields of Computer Aided Design (CAD) and Finite Element Analysis (FEA). The core idea is to resort to the same higher-order and smooth bases, such as made of B-Spline and Non-Uniform-Rational-B-Spline (NURBS) functions (Cohen et al., 1980, Piegl and[START_REF] Piegl | The NURBS book. 2Nd[END_REF], for the representation of the geometry in CAD as well as for the approximation of solution fields in FEA. Basically, this is the main difference between FEA and IGA. In a quite simplified way, this idea is shown on Figure 1.1. If FEA takes Lagrange polynomials used to discretize the solution and employs them to create model's geometry, IGA goes in the opposite direction: starting with CAD geometry, it takes functions used in this geometry (such as B-Splines) to discretize the solution.

Figure 1.1: Schematic representation of the main idea of IGA and FEA: FEA uses polynomials to approximate the geometry, whilst IGA utilizes functions from CAD to discretize the solution.

The use of such families of functions quickly made IGA highly attractive for engineers since a common geometrical model can be used by both designers and analysts, thereby facilitating the dialogue between their two worlds. Indeed, in a classical FEA procedure the step of creating a Finite Element (FE) model from a CAD geometry is often extremely tedious and time-consuming. It is known to be the major bottleneck in engineering analysis procedures [START_REF] Cottrell | Isogeometric Analysis: Toward Integration of CAD and FEA[END_REF]. More of it, each modification or update to the initial model geometry requires restarting this step, leading to significant inefficiencies and delays. However, utilizing CAD models directly for engineering analysis is not always a seamless process. Note that CAD models are primarily designed for visualization and geometric accuracy, without consideration for the specific requirements of engineering analysis. Analysis-aware model requires a boundary fitted discretization while there is no such a need for geometry models in CAD: in this case, the geometry needs only to be rendered correctly. Therefore, numerous approaches are being developed these days to truly arrive at a CAD-analysis link. In this work, we will consider fictitious domain type methods as a solution to this problem. More details 1.1. IGA FROM A FEA PERSPECTIVE 11 regarding this point will be given further with an emphasis on local enrichment (see Figure 1. 3 and related paragraph).

Even though the original aim was to streamline the time-consuming process of creating mechanical models from CAD programs, the use of splines offers its own advantages from an analysis point of view. Indeed, spline functions can be C (p-1) regular between elements for a polynomial degree p, while Lagrange polynomials, which are used in the standard FEA, attain only a C 0 regularity at these locations. As a result, IGA exhibits higher per-degree-of-freedom accuracy and robustness compared to the standard FEA for mechanical simulations [START_REF] Evans | n-Widths, sup-infs, and optimality ratios for the k-version of the isogeometric finite element method[END_REF]. Owing to this superior behaviour, IGA has found its application in various branches of mechanics. One can name contact mechanics (Antolin et al., 2019a[START_REF] Seitz | Isogeometric dual Mortar methods for computational contact mechanics[END_REF][START_REF] Temizer | Contact treatment in isogeometric analysis with nurbs[END_REF][START_REF] Wriggers | Smooth cĺ1-interpolations for twodimensional frictional contact problems[END_REF], non-linear solid mechanics [START_REF] Ambati | Isogeometric Kirchhoff-Love shell formulation for elasto-plasticity[END_REF], Bouclier et al., 2015a[START_REF] Elguedj | Isogeometric analysis of nearly incompressible large strain plasticity[END_REF], shell analysis (Bouclier et al., 2015b[START_REF] Echter | A hierarchic family of isogeometric shell finite elements[END_REF][START_REF] Kiendl | Isogeometric shell analysis with Kirchhoff-Love elements[END_REF], fluid-structure interaction [START_REF] Apostolatos | Systematic evaluation of the interface description for fluid-structure interaction simulations using the isogeometric mortar-based mapping[END_REF][START_REF] Kamensky | An immersogeometric variational framework for fluid-structure interaction: Application to bioprosthetic heart valves[END_REF], modelling elastic solids undergoing chemical transformations [START_REF] Morozov | Modeling chemical reaction front propagation by using an isogeometric analysis[END_REF], Digital Image Correlation [START_REF] Dufour | CAD-based displacement measurements with stereo-dic: principle and first validations[END_REF][START_REF] Réthoré | On the use of NURBS functions for displacement derivatives measurement by digital image correlation[END_REF], and shape optimization [START_REF] Kiendl | Isogeometric shape optimization of shells using semi-analytical sensitivity analysis and sensitivity weighting[END_REF][START_REF] Nagy | Isogeometric design of anisotropic shells: optimal form and material distribution[END_REF][START_REF] Wall | Isogeometric structural shape optimization[END_REF].

The relationship between IGA and FEA can be viewed from different angles: considering technologies or the resulting approximation spaces. From a technological perspective, IGA provides the flexibility to employ basis functions with up to C p-1 regularity, where p represents the polynomial degree, in contrast to only C 0 continuous polynomials in FEA. Thus, IGA seems to encompass FEA. Concerning discretization, the potentially more regular approximation spaces offered by IGA are included within the C 0 spaces of classic FEA, provided that both IGA and FEA come with the same number of elements and polynomial degree. As a result, IGA can also be interpreted as a projection of FEA onto a specific reduced basis.

This section undertakes to present IGA from the second point of view (Bouclier andPassieux, 2023, Tirvaudey et al., 2020). First, fundamental elements of IGA are revisited. Then, the existing link between IGA and FEA is reviewed through the presentation of the initial Bézier and more recent Lagrange (Schillinger et al., 2016b) extractions. Finally, advantages and limitations of both approaches are summarized.

Remark 1.1.1

It is important to note that in this work, we restrict ourselves to B-Splines and NURBS. These are the technologies widely used in CAD. However, other geometry descriptions exist, such as T-Splines [START_REF] Bazilevs | Isogeometric analysis using T-splines[END_REF] or U-Splines [START_REF] Thomas | U-splines: Splines over unstructured meshes[END_REF]. Their advantage in comparison with B-Splines and NURBS is the possibility to perform local refinement. Indeed, while IGA offers efficient strategy for global refinement (which is discussed further, see Section 1.1.2), local refinement is not straightforward due to the rigid tensorproduct structure of multi-variate basis functions (see Eq. (1.5) and related paragraph). However, in this work, we do not proceed to the local refinement; thus, standard B-Splines and NURBS are employed.

Basics of NURBS-based IGA

We begin with key elements of the B-spline and NURBS geometric modelling techniques. As IGA is now mature and relatively well-known in the scientific computing community, we do not enter into the details of the technology here. For more information, the interested reader is referred to the seminal contributions [START_REF] Cottrell | Isogeometric Analysis: Toward Integration of CAD and FEA[END_REF], 2005) and e.g., [START_REF] Bouclier | IGA: Non-conforming Coupling and Shape Optimization of Complex Multipatch Structures[END_REF] for a recent review and to the works cited hereafter. Widely used in CAD for geometry modelling, NURBS allow the accurate representation of many shapes, such as conic sections (circles, cylinders, spheres, ellipsoids, etc.). NURBS can be considered as a generalisation of B-Splines. Since NURBS are based on B-Splines, we start the presentation below with B-Splines.

B-Splines

A set of non-decreasing real numbers where B-Spline basis functions are defined is called a knot vector, Ξ

1 = ξ 1 1 , ξ 2 1 , ...ξ n+p+1 1
, where ξ i 1 is the i-th knot, n is the number of basis functions and p is the polynomial order. Interval [ξ 1 1 , ξ n+p+1 1 ] forms an isogeometic patch, while coordinates ξ i 1 , called knots, divide the patch into knot-span elements [ξ i 1 , ξ i+1 1 ]. Whilst in standard FEA each element has its own parametrization, B-Splines possess a parametric space that is confined to the patch level, akin to the notion of a macro-element.

A knot vector with the first and the last knots which appear p + 1 times makes the basis functions interpolatory at the ends of the parametric space interval. This knot vector is called an open knot vector and is commonly used in CAD packages.

The univariate B-spline basis functions N p i (ξ) of order p are defined recursively by the Cox-de-Boor formula (Cohen et al., 1980, Piegl and[START_REF] Piegl | The NURBS book. 2Nd[END_REF]:

       N 0 i (ξ 1 ) = 1 if ξ i 1 ≤ ξ 1 < ξ i+1 1 , 0 otherwise, p = 0 N p i (ξ 1 ) = ξ1-ξ i 1 ξ i+p 1 -ξ i 1 N p-1 i (ξ 1 ) + ξ i+p+1 1 -ξ1 ξ i+p+1 1 -ξ i+1 1 N p-1 i+1 (ξ 1 ), p ≥ 1.
For p = 0 and p = 1 the IsoGeometric basis functions coincide with those for standard piecewise constant and linear FE. However, for p ≥ 2 they are different.

It is important to note that unlike FEA basis functions, which always lie on the geometry, splines are not (in general) interpolatory. If it is necessary to ensure that an interior knot is interpolatory, its multiplicity must be increased. When the multiplicity of a knot value is precisely p, the basis function is interpolatory at that knot.

Figure 1.2(right) depicts quadratic splines formed from an uniform open knotvector. It can be seen that the functions are always non-negative, N p i (ξ 1 ) ≥ 0, ∀ξ 1 and that the basis constitutes partition of unity:

n i=1 N p i (ξ 1 ) = 1, ∀ξ 1 .
They have larger support than FEA basis functions do (see Figure 1.2(left) for FEA quadratic basis functions). Indeed, for a B-spline of order p, the support is p + 1 elements. However, for 1.1. IGA FROM A FEA PERSPECTIVE 13 both IGA and FEA, each basis function shares its support with 2p functions. Thanks to this, the use of B-Splines does not increase bandwidth in numerical methods.

B-Spline curves

B-spline curves in R d (where d is the dimension of the problem) are constructed by taking a linear combination of B-Spline basis functions. The vector-valued coefficients of the basis functions are referred to as control points. This is similar to FEA, where nodal coordinates are the coefficients of the basis functions. However, due to the noninterpolatory nature of the basis, there is no direct interpretation of the control point values. Given n basis functions, N p i (ξ 1 ) , i = 1, 2, ..., n, and corresponding control points

P i ∈ R d ; i = 1, 2, ..., n, a piecewise-polynomial B-Spline curve C BS (ξ 1
) is given by:

C BS (ξ 1 ) = n i=1 P i N p i (ξ 1 ). (1.1)
Note here that P i can live either in a 1D (

P i = x i ), 2D (P i = (x i , y i )) or 3D (P i = (x i , y i , z i )).
Eq. (1.1) can be rewritten in matrix form:

C BS (ξ 1 ) = N T x BS , (1.2)
where N is the matrix of B-Spline basis functions and x BS is the vector collecting the locations of the associated control points. For a 2D example, N and x BS formally read:

N =            N 1 0 . . . . . . N n 0 0 N 1 . . . . . . 0 N n            , x =            x 1 . . .
x n y 1 . . .

y n           
.

Non-Uniform Rational B-Splines

In addition to B-Spline geometries, NURBS offer a precise means of representing conics. NURBS functions, being a generalization of B-Splines, share several fundamental properties with B-Splines, including enhanced smoothness.

Starting with a B-Spline curve, each control point P i is associated with a weight ω i . The NURBS basis functions are defined as follows:

R i (ξ 1 ) = N p i (ξ 1 )ω i n j=1 N p j (ξ 1 )ω i . . (1.3)
With these basis functions, the general formula of a NURBS curve has the same form as for a B-Spline curve :

C N U RBS (ξ 1 ) = n i=1 P i R i,p (ξ 1 ). (1.4)
Notably, when all weights are equal, the NURBS entity simplifies to a B-Spline entity. Eq. (1.4) can be rewritten in matrix form as follows:

C N U RBS (ξ 1 ) = R T x N U RBS ,
where R is the matrix of NURBS basis functions and x N U RBS is the vector collecting the locations of the associated control points.

Multivariate basis functions

Of course, B-Splines and NURBS are not restricted to parametric curves. Multivariate entities such as surfaces and volumes can also be described in the parametric form by using these technologies.

Bivariate basis functions are obtained by a tensor-product of univariate ones:

N pq i (ξ 1 , ξ 2 ) = N p i1 (ξ 1 ) × N q i2 (ξ 2 ).
(1.5)

Basis functions N p i1 (ξ 1 ) and N q i2 (ξ 2 ) are univariate B-Spline basis functions of degree p and q. They are formed from two different knot vectors, Ξ 1 of size n 1 + p + 1 and Ξ 2 of size n 2 + q + 1, n 1 and n 2 being the number of basis functions associated to directions ξ 1 and ξ 2 , respectively. Thus, the number of bivariate B-Spline functions is

n 2D = n 1 × n 2 .
With the use of bivariate basis functions, a B-Spline surface if defined as:

S BS (ξ 1 , ξ 2 ) = n 2D i=1 N pq i (ξ 1 , ξ 2 )P i .
The tensor nature of these entities has significant implications on the regularity of surfaces (C 0 points in one direction propagate along a line into the other direction and form C 0 lines) and on local mesh refinement.

As for the univariate case, the construction of a NURBS surface requires the additional definition of one weight per control point. More specifically, a NURBS surface reads:

S N U RBS (ξ 1 , ξ 2 ) = n 2D i=1 R pq i (ξ 1 , ξ 2 )P i ,
with R pq i (ξ 1 , ξ 2 ) a bivariate NURBS basis function, defined in a similar manner as in Eq. (1.3) but from bivariate B-spline functions:

R pq i (ξ 1 , ξ 2 ) = N pq i (ξ 1 , ξ 2 )ω i n 2D j=1 N pq j (ξ 1 , ξ 2 )ω j .
. Therefore, the tensor product operation is also involved for NURBS surfaces. As for NURBS curves, the weights give more shape control on the surface and allow it to describe exactly, for example, cylinders and spheres.

Increased continuity and its benefits

One of the compelling advantages of splines is their higher degree of regularity. Specifically, a B-Spline function with a degree of p can attain a C p-1 regularity at knot ξ i 1 if that knot is singular in Ξ 1 . As a consequence, when comparing meshes with the same polynomial degree and number of elements, a B-Spline mesh with C p-1 continuity requires fewer degrees of freedom than the corresponding C 0 FE mesh. This is entirely logical, as the space of C p-1 functions is included into the space of C 0 functions. For instance, in the case of a 3D solid mesh composed of 100 elements in each of the three parametric directions, we achieve a reduction in number of degrees of freedom by a factor of approximately 8, 25, and 57 for p = 2, 3, and 4, respectively [START_REF] Bouclier | IGA: Non-Invasive Coupling with FEM and Regularization of Digital Image Correlation Problems[END_REF]. This distinctive feature of IGA results in significantly enhanced accuracy per degree of freedom when dealing with smooth solutions, compared to FEA.

Note that there can be more than one knot at a given location of the parametric space. More precisely, if m is the multiplicity of a given knot, the functions are C p-m continuous at that location which is in contrast with standard FE where only a C 0 regularity is encountered on the element boundaries. Thus, if m = p, functions are C 0 at the given knot and they are interpolatory.

Global refinement

Spline functions present efficient global refinement procedures which allow enhancement of the approximation space without changing the geometry. To perform a refinement, two techniques exist: knot insertion and order elevation.

Knot insertion in IGA is similar to the h-refinement in FEA. New knots may introduce new values in the existing knot vector Ξ or duplicate already presented ones. In the latter case, their multiplicity increases and the continuity of basis functions at this point is reduced. This procedure coincides with the h-refinement in FEA if each new knot is inserted p times (so that basis functions are C 0 at elements boundaries). The new refined version of the B-Spline curve takes the same degree p but is built from the new knot-vector which concatenates the initial knots with the inserted ones in ascending order.

Order elevation is the analogous of the p-refinement in FEA. During this procedure, the polynomial degree of the basis functions p is increased. However, the initial regularity of the basis functions is preserved: all values in the knot-vector are repeated. Thus, if one starts with a spline mesh with p = 1 and C 0 at the knots, order elevation techniques coincides with p-refinement in FEA (the basis functions with augmented order p will still be only C 0 . ) However, refinement procedure in IGA are not restricted to these two. There exist a third procedure, k-refinement, which is specific to IGA. In this case, degree-elevation is performed first, and then the new knot values are inserted (with multiplicity one) so that elements are added while ensuring the maximum available regularity of the basis functions at the knots, namely C p-1 . An interested reader can find more details about refinement procedures, e.g., in Cottrell et al. (2007b), [START_REF] Piegl | The NURBS book. 2Nd[END_REF]. In this work, unless otherwise stated, we will perform k-refinement to take advantage of the superior properties of splines Remark 1.1.2

In practice, one can use matrix representation of the spline refinement procedures [START_REF] Piegl | The NURBS book. 2Nd[END_REF]). Thus, denoting by R c and R f (resp. x c and x f ), the matrices (resp. vectors) collecting the coarse and fine spline functions (resp. control points), we can build the refinement operator D IG cf such that:

R c = D IG cf R f and x f = D IG cf T x c .
For more details on refinement strategies of splines and their matrix representations, reference is made to [START_REF] Bouclier | IGA: Non-conforming Coupling and Shape Optimization of Complex Multipatch Structures[END_REF], Cottrell et al. (2007b), [START_REF] Piegl | The NURBS book. 2Nd[END_REF].

IsoGeometric principle (splines in FEA)

IGA uses the same variational formulations as FEA. From the analysis point of view, IGA can be seen as a generalization of FEA that considers smooth and higher-degree functions instead of typical Lagrange polynomials. With the same idea as in FEA, the weak form is converted to systems of linear equations with the use of the spline basis. Thus, the displacement field is approximated as:

u h (ξ 1 ) = n i=1 R i (ξ 1 )u i ,
where R i (ξ 1 ) are the IsoGeometric basis functions and u i represent control points values. As it stems from the core idea of IGA, it is the geometry of the problem that indicates the choice of the basis functions. Therefore, we end up with the linear system: Just like in FEA, K is called the stiffness matrix, u is the displacement vector (displacement values at degrees of freedom -in this case, control points, not nodes), f is the force vector (also defined at control points).

Remark 1.1.3

Numerical integration is important for building the stiffness matrix. In this study, for the sake of simplicity and non-invasiveness with respect to current practices in FEA, we employ standard Gauss quadrature, utilizing p + 1 Gauss points per element. Nevertheless, it's important to highlight that this approach may not be the most optimal choice for spline functions. There exist advanced integration rules, see, e.g., [START_REF] Adam | Selective and reduced numerical integrations for NURBS-based isogeometric analysis[END_REF], [START_REF] Auricchio | A simple algorithm for obtaining nearly optimal quadrature rules for NURBS-based isogeometric analysis[END_REF], [START_REF] Bouclier | Locking free isogeometric formulations of curved thick beams[END_REF], [START_REF] Hughes | Efficient quadrature for NURBS-based isogeometric analysis[END_REF], [START_REF] Schillinger | Reduced Bézier element quadrature rules for quadratic and cubic splines in isogeometric analysis[END_REF].

Enforcement of boundary conditions in IGA

As previously mentioned, control points in IGA are not interpolatory, which stands in contrast to FEA nodes that consistently align with the geometry of the model. In FEA, Dirichlet boundary conditions (both homogeneous and non-homogeneous) are prescribed directly at the nodes defining the boundaries. In IGA, only homogeneous boundary conditions can be applied in the same way, i.e. directly at the control points. This is due to the fact that if based on an open knot vector, spline basis functions are interpolatory and C 0 -continuous on its end. Hovewer, to apply non-homogeneous boundary conditions, a special treatment is requiered. One can consider variational forms for handling boundary conditions in IGA, such as penalty or Lagrange multipliers methods.

Spline discretization and challenges: illustration for local enrichment

Even if based on the same basis functions, the practical realization of the CADanalysis connection, as initially advocated by IGA, presents non-trivial challenges. Basically, as stated above, the principle of IGA is to use B-spline and NURBS functions to build the approximation spaces when applying Galerkin's method. Although this can provide higher accuracy and robustness compared to the standard FEA, this also, unfortunately, makes the local modelling of geometrically complex inclusions within a spline patch highly challenging. This point seems to be closely related to what is called the analysis-suitable model issue in the field. Indeed, standard IGA requires a boundary-fitted discretization for the analysis while in CAD programs, where the only matter is the rendering of the geometry, entities are described as collections of their boundary surfaces.

As an illustration, let us take the simple example of Fig. 1.3. In CAD, this object may consist of: (i) a one-patch B-spline surface for the whole plate (linear, 1 element) and (ii) a NURBS trimming curve (quadratic, 2 elements) that forms the boundary between the global and local domains Ω m and Ω f , respectively (see Fig. 1.3(a)). As a consequence, the underlying spline surface is unaffected by the trimming object and preserves its topology. Conversely, using standard IGA for the analysis of such a problem would require a delicate spline re-parametrization of the whole model to explicitly define the two domains. In practice, this would inevitably lead to the splitting of the geometry into several (tensor-product) patches with C 0 continuity at the boundaries (see Fig. 1.3(b) for an example of a boundary fitted NURBS parametrization of the considered problem). In the case of multiple inclusions, the situation is getting even trickier since the resulting geometry of the matrix strongly differs topologically from a square. The new spline model of Fig. 1.3(b) is commonly referred to as an analysis-suitable model in the sense that it can be easily enhanced using classic spline refinement [START_REF] Bouclier | IGA: Non-conforming Coupling and Shape Optimization of Complex Multipatch Structures[END_REF], Cottrell et al., 2007a[START_REF] Lee | Degree elevation of nurbs curves by weighted blossom[END_REF] to compute the solution of a corresponding mechanical problem. This simple example underlines the difficulties of generating analysis-suitable models for multi-scale IGA. One can perform as in Fig. 1. 3(b); that is, one can strive to remove all trimmed regions by invoking spline re-parametrization strategies, which generally appears cumbersome in practice. Alternatively, one can implement specific numerical schemes adapted to the models with the mentioned defects coming from the geometric modelling [START_REF] Marussig | A review of trimming in isogeometric analysis: Challenges, data exchange and simulation aspects[END_REF]. The current works goes in the latter direction, resorting to fictitious domain methods (see Figure 1.3(c) for an illustration and Section 1.2 for a further discussion).

Establishing the Link between IGA and FEA

In this section, we aim to bridge the gap between IGA and FEA, helping readers who are already familiar with FEA gain a better understanding of IGA. Furthermore, this 18 CHAPTER 1. STATE OF THE ART: ON THE MULTI-SCALE MODELLING AND COMPUTATION OF COMPLEX STRUCTURES WITH IGA connection will play a pivotal role in Chapter 2, where we will establish a global-IGA/local-FEA coupling which is non-intrusive with respect to FEA.

To establish a connection between IGA and FEA, we resort to a global algebraic bridge between IGA and FEA [START_REF] Colantonio | Noninvasive multilevel geometric regularization of mesh-based three-dimensional shape measurement[END_REF][START_REF] Tirvaudey | Non-invasive implementation of nonlinear isogeometric analysis in an industrial FE software[END_REF]. This bridge relies on the utilization of Bézier-based operators [START_REF] Borden | Isogeometric finite element data structures based on Bézier extraction of nurbs[END_REF][START_REF] Kamensky | tIGAr: Automating isogeometric analysis with FEniCS[END_REF], Schillinger et al., 2016b). Moreover, we show the connection between IsoGeometric operators (stiffness matrix and load vector) and their FE counterparts. It allows the implementation of IGA even without access to a global IsoGeometric code.

Viewpoint of the technologies

Bézier extraction

Bézier extraction allows to formulate a smooth B-spline (resp. NURBS) function in terms of C 0 polynomial Bernstein (resp. rational Bernstein) functions. Considering that the space generated by smooth spline functions is included into the one generated by C 0 functions, the idea here is to extract the smooth part from C 0 functions.

Let us first consider B-Spline functions. Thus, the goal is to map a Bernstein polynomial basis into a smooth B-Spline polynomial basis. Starting with B-Spline functions, it is easy to build C 0 polynomials: it suffices to repeat all the interior knots of the knot vectors until they reach a p multiplicity. The advantage of Bernstein functions is that they exhibit an elementary structure which is similar to FEA.

This procedure can be written in a matrix form. Denoting N (resp. B) matrix collecting B-Spline (resp. Bernstein) functions, it reads:

N = D Béz B, (1.6)
where D Béz is the corresponding Bézier extraction operator. New control points should also be computed. To do so, we remind that the initial B-Spline geometry (defined by Eq. (1.2)) remains unchanged during the procedure. Consequentially, the expression for control points can be obtained by inserting Eq. (1.6) in Eq. (1.2):

G(ξ) = N T x BS = B T D T Béz x BS = B T
x Ber , with x Ber = D T Béz x BS . In the equation above, x BS (resp. x Ber ) stands for location of the B-Spline (resp. Bernstein) control points.

Lagrange extraction

Here, the idea is to extend the concept of Bézier extraction to standard nodal FE functions (Schillinger et al., 2016b), effectively building the bridge between B-Splines and Lagrange polynomials. With Bézier extraction in hand, we first present a passage from Lagrange to Bernstein polynomials, and then construct the entire Lagrange-to-B-Splines operator by combining these two technologies.

To construct Lagrange-to-Bernstein operator, one simply need to express the Bernstein functions as a linear combination of the Lagrange functions at some interpolation points. Making use of the interpolatory property of the Lagrange functions, this operator can be efficiently generated by evaluating Bernstein functions at nodal points corresponding to the Lagrange basis. Denoting by D LB such an operator and by L the classical FE basis functions, the connection is expressed as follows:

B = D LB L.
Then, with the same reasoning as for Bézier extraction, FE nodes can be derived from the Bernstein mesh. The FE nodes, obviously, interpolate the geometry and are obtained as follows:

x

F E = (D LB ) T x Ber .
With this two ingredients in hand (Bézier extraction operator B Béz and Lagrangeto-Bernstein operator D LB ), the Lagrange extraction operator can be viewed as a combination of them:

N = DL with D = D Béz D LB .
(1.7)

The FE nodes are obtained as:

x F E = D T x BS .
Lagrange extraction (Schillinger et al., 2016b[START_REF] Tirvaudey | Non-invasive implementation of nonlinear isogeometric analysis in an industrial FE software[END_REF] allows to formulate a smooth polynomial B-Spline discretization in terms of a standard FE discretizations. In other words, with the Lagrange extraction, we somehow extract the smooth part of C 0 polynomials to derive the C (p-1) B-Spline basis functions.

Remark 1.1.4

While Eq. (1.7) provides a clear understanding of the Lagrange extraction operator D LB , it does not reflect the way how it is built in practice. In practice, constructing D LB does not involve going through Bernstein polynomials. Instead, it simply requires evaluating the B-Spline basis functions at the nodal points associated with the Lagrange polynomials Schillinger et al. (2016b).

The extraction in case of NURBS

The Lagrange extraction can be extended to the case of NURBS. Following Eq. (1.3), NURBS basis functions matrix R can be found from B-Spline functions matrix as:

R = WN W , (1.8)
where W is the diagonal matrix of NURBS weights, and W = n i=1 N p i ω i . With Eq. (1.7) B-Spline functions can be replaced with Lagrange polynomials, thus:

R = WDL W .
(1.9)

A link between NURBS weights and rational Lagrange functions weights can be made. To do so, let us express function W in a matrix form, with w as a vector, collecting NURBS weights: (1.10) where w L = D T w is a vector collecting rational Lagrange functions weights. Similarly to NURBS in Eq. (1.8), rational Lagrange polynomials L R can be obtained from L as follows:

W = n i=1 N p i ω i = w T N = w T DL = (D T w) T L = (w L ) T L = W L ,
L R = W L L W L , ( 1
R = WD(W L ) -1 L R = D R L R .
This new extraction operator D R = WD(W L ) -1 links NURBS with rational Lagrange functions. Finally, the rational Lagrange control points depend on the NURBS control points: x RL = D T R x N U RBS . We emphasize at this stage that this link from rational Lagrange functions to NURBS is exact since both bases are rational and NURBS are of higher-order smoothness. To truly involve Lagrange polynomials, a projection allowing to formulate rational functions in terms of polynomials is required. To do so, operator D LL is introduced as:

L R = D LL L.
(1.12)

Details about construction of this operator can be found in [START_REF] Tirvaudey | Non-invasive implementation of nonlinear isogeometric analysis in an industrial FE software[END_REF].

The FE nodes can be obtained as:

x F E = (D LL ) T x RL . (1.13) Remark 1.1.5
Let us underline that equalities (1.12) and (1.13) do not strictly hold in case of NURBS. Indeed, moving from a rational to a polynomial geometry necessarily leads to some approximations. However, as demonstrated in [START_REF] Tirvaudey | Non-invasive implementation of nonlinear isogeometric analysis in an industrial FE software[END_REF], the error related to such approximations is largely insignificant when considering refined geometries compared to the associated NURBS discretization error.

Figure 1.4: Link between IsoGeometric (B-Splines or NURBS) mesh and a FE mesh with the use of different global operators. The left side of the figure illustrates the process of deriving the IsoGeometric mesh from an FE mesh. Then, the same operators are used to recover IsoGeometric stiffness operator from the FE one computed using a classical FE software, taken as a black-box (depicted on the right side). The strategy for B-Splines is outlined in the first row with black arrows, while the strategy for NURBS is illustrated in the second row with blue arrows. Solid lines represent an exact link, while dashed lines represent an approximation.

Figure 1.4(left) illustrated the process of creating the FE mesh from initial IsoGeometric(B-Splines or NURBS) mesh. Once the FE mesh is obtained, classical FE software can be used to compute FE operator (stiffness matrix and load vector). Then, IsoGeometric operators can be derived from the FE operators using the same 1.1. IGA FROM A FEA PERSPECTIVE 21 operators that were initially used to create the FE mesh from the IsoGeometric one. For B-Splines functions, the passage can be done as:

K BS = DK F E D T , f BS = Df F E ,
where K BS (resp. K F E ) and f BS (resp. f F E ) are the stiffness matrix and the load vector associated with B-Spline (resp. FE) mesh.

In case of NURBS the procedure to recover IsoGeometric operators from FE ones involves a two-step process. First, the stiffness matrix K RL and the load vector f RL associated with the rational Lagrange mesh are computed:

K RL = D LL K F E (D LL ) T , f RL = D LL f F E .
Then, the stiffness matrix K N U RBS and the load vector f N U RBS associated with the NURBS mesh are derived as follows:

K N U RBS = D R K RL (D R ) T , f N U RBS = D R f RL .
Thus, FEA software can be used to build FE operators, and then they can be transformed to B-Splines or NURBS operators. This is illustrated on Figure 1.4(right). Then, the corresponding linear system (K BS u BS = f BS for B-Splines and K N U RBS u N U RBS = f N U RBS ) can be solved to obtain the displacement vector.

Finally, let us note that the FE nodal displacement vector u F E can be derived from the displacement vector associated to the IsoGeometric model. Thus, in case of B-Splines:

u F E = D T u BS ,
and in case of NURBS:

u F E = (D LL ) T (D R ) T u N U RBS .
This implies that once u F E obtained, classical FE software can be used for postprocessing.

Finally, let us generalize the link for any case of IsoGeometric basis functions:

u F E = (D F E ) T u IG , (1.14a) K IG = D F E K F E (D F E ) T , f IG = D F E f F E ,
where u IG is the displacement vector associated to the IsoGeometric mesh and D F E is an appropriate passage operator with regard to the type of the IsoGeometric basis functions.

Remark 1.1.6

While this work focuses on B-Splines and NURBS, FEA-to-IGA bridge and its use for achieving a non-invasive implementation with respect to FEA could be straightforwardly extended to many other spline technologies. These may include T-Splines [START_REF] Kamensky | tIGAr: Automating isogeometric analysis with FEniCS[END_REF], Scott et al., 2011), hierarchical B-Splines and NURBS (D'Angella and Reali, 2020[START_REF] Hennig | Bézier extraction and adaptive refinement of truncated hierarchical NURBS[END_REF], and hierarchical T-Splines (Chen andDe Borst, 2018, Evans et al., 2015), since the operators of Bézier extraction are available for these technologies as well.

Advantages and limitations of IGA and FEA

After this overview on IGA with respect to FEA, , let us summarize advantages and limitations of both IGA and FEA. FEA benefits from more that 50 years of developments and a myriad of softwares based on this approach is presented today, making FEA at hand for engineers. IGA, on its side, is not as easily accessible to the general public yet. Hitherto, there exist possibilities to perform IGA simulation with the use of existing FEA software (e.g., , LS-Dyna [START_REF] Dj Benson | Blended isogeometric shells[END_REF][START_REF] Chen | Current status of lsdyna isogeometric analysis in crash simulation[END_REF][START_REF] Hartmann | About isogeometric analysis and the new NURBS-based finite elements in LS-DYNA[END_REF][START_REF] Lozinski | A new fictitious domain method: Optimal convergence without cut elements[END_REF], Abaqus [START_REF] Duval | abqNURBS: implémentation d'éléments isogéométriques dans abaqus et outils de pré-et post-traitement dédiés[END_REF], Elguedj et al., 2012[START_REF] Lai | Integrating CAD with Abaqus: a practical isogeometric analysis software platform for industrial applications[END_REF], Radioss [START_REF] Occelli | LR B-splines implementation in the Altair RadiossTM solver for explicit dynamics IsoGeometric Analysis[END_REF]), but they often require additional efforts, such as user-defined elements [START_REF] Khakalo | Isogeometric analysis of higher-order gradient elasticity by user elements of a commercial finite element software[END_REF]. As for research code, different codes exist, such as IGAFEM for MatLab [START_REF] Nguyen | Isogeometric analysis: an overview and computer implementation aspects[END_REF], libraries PyIGA (https://github.com/c-f-h/pyiga) and Nutils (http://www.nutils.org) for Pyfthon. The use of splines for analysis offer several benefits in comparison with Lagrange polynomials:

• the exact representation of the geometry minimizes the geometrical discretization errors.

• the increased continuity of basis functions allows to achieve better accuracy with the same number of degrees of freedom compared to FEA for regular solutions.

• in contact mechanics, it allows a smooth representation of contact surfaces, while in FEA, surfaces have faceted representations which results in discontinuities in the normal vector. FEA can suffer from high jumps and artificial oscillations in the value of contact tractions (De Lorenzis et al., 2014). In contrast, IGA, thanks to the higher-order continuity of spline basis functions, exhibit more accurate results in contact problem modelling when compared to the traditional FE approach [START_REF] Corbett | NURBS-enriched contact finite elements[END_REF]Sauer, 2014, De Lorenzis et al., 2014).

However, there are some considerations to be made. First, tackling complex geometries within IGA can be a non-trivial task (see again Figure 1.3). It often requires such as multipatch coupling or immersed boundary techniques. In contrast, FEA demands only a sufficiently fine mesh to adequately represent the geometry, making it a more straightforward process in this regard. Secondly, the increased continuity of basis functions may not necessarily be beneficial when dealing with sharp discontinuities or singularities in the solution. In such cases, FEA might be a more suitable choice.

Immersed Boundary Methods for Modelling Complex Geometries

In this section, we address the issue of modelling geometrically intricate structures. As previously discussed in Section 1.1.2 (see Figure 1.3), the modelling of complex geometries within the IGA framework poses significant challenges. In this section, we review the Immersed Boundary approach as a viable strategy to address these challenges. Given that geometrically complex models often require the use of multiple patches and, consequently, the coupling of these patches, we also present here weak coupling techniques.

Introduction

Conventional Partial Differential Equations (PDE) system resolution techniques use meshes, for which the computational domain is based on the physical domain under study. These meshes aim at accurately capturing the real geometry of the physical domain. Hence, when the geometry of the computational domain is complex, a highly detailed mesh is needed to accurately represent it. Thus, the primary drawback of such methods lies in their substantial computational cost both in terms of computational time and storage. To address these limitations, the immersed boundary approach was proposed. Initially, it was introduced within the context of the finite difference method to model blood flow around the human heart valve [START_REF] Peskin | Flow patterns around heart valves: a numerical method[END_REF]. One distinguishing feature of this approach is the ability to perform the entire simulation on a fixed Cartesian mesh (see Figure 1.5 for an illustration). The core idea consists in immersing any complex geometry in a more simple domain (usually a rectilinear one) and representing the solution of the physical problem using the degrees of freedom related to the embedding mesh. The main advantage of this method is that it alleviates the meshing step. This facilitates the analysis of geometrically complex and time-evolving problems. However, to obtain a satisfactory approximate solution over the physical domain, it is essential to accurately account for the physical boundaries. To do so, specific integration procedures are developed. Thus, the step of taking into account the physical boundaries of the considered domain is shifted from meshing to integration. Following Peskin's pioneering work, various adaptations and enhancements have emerged, leading to a number of variants of this approach. However, it was not until the early 2000s that the immersed approach was applied to FEA, notably through the concept of cutting finite elements as an unfitted meshing technique pioneered by Hansbo [START_REF] Hansbo | An unfitted finite element method, based on Nitsche's method, for elliptic interface problems[END_REF]. The pace in the development and impact of immersed finite element methods increased significantly with the introduction of the Finite Cell Method (FCM) [START_REF] Divi | Topology-preserving scan-based immersed isogeometric analysis[END_REF], Rank et al., 2012[START_REF] Schillinger | The finite cell method: A review in the context of higherorder structural analysis of CAD and image-based geometric models[END_REF][START_REF] Verhoosel | Imagebased goal-oriented adaptive isogeometric analysis with application to the micromechanical modeling of trabecular bone[END_REF], which combines the cut element concept with higher-order basis functions, and CutFEM (Burman andHansbo, 2012, Burman et al., 2015).

Because representing complex geometries in IGA can be challenging (as it was shown on Figure 1.3), immersed boundary approaches have emerged as an area of significant interest within this field. Recent years have witnessed the integration of immersed analysis concepts with IGA [START_REF] Hughes | Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement[END_REF][START_REF] Marussig | A review of trimming in isogeometric analysis: Challenges, data exchange and simulation aspects[END_REF], Rank et al., 2012[START_REF] Schillinger | The finite cell method: A review in the context of higherorder structural analysis of CAD and image-based geometric models[END_REF], often referred to as IGA-FCM [START_REF] Rank | Geometric modeling, isogeometric analysis and the finite cell method[END_REF] or immersogeometric analysis [START_REF] Hoang | Skeleton-stabilized immersogeometric analysis for incompressible viscous flow problems[END_REF][START_REF] Kamensky | An immersogeometric variational framework for fluid-structure interaction: Application to bioprosthetic heart valves[END_REF]. This approach has demonstrated the ability to harness the favourable properties of spline basis functions used in IGA, while enhancing the versatility of the simulation workflow, particularly in scenarios where boundary-fitting spline geometries are not readily accessible, such as in scan-based analyses.

In this thesis, a focus is performed on the FCM because its IsoGeometric version based on B-splines [START_REF] De Prenter | Condition number analysis and preconditioning of the finite cell method[END_REF][START_REF] Schillinger | The finite cell method: A review in the context of higherorder structural analysis of CAD and image-based geometric models[END_REF], Schillinger et al., 2012) opens multiple attractive applications thanks to the structured aspect of the embedding mesh. Examples of using FCM combined with IGA include analysis of flow in porous media [START_REF] Hoang | Skeleton-stabilized immersogeometric analysis for incompressible viscous flow problems[END_REF], fluid-structure interactions in vesicles [START_REF] Casquero | The divergence-conforming immersed boundary method: Application to vesicle and capsule dynamics[END_REF], and composite kirchhoff plates [START_REF] Patton | Cost-effective and accurate interlaminar stress modeling of composite Kirchhoff plates via immersed isogeometric analysis and equilibrium[END_REF], to name a few. In the context of image-based modelling, the FCM has been successfully applied to the modelling of trabecular bones [START_REF] Divi | Topology-preserving scan-based immersed isogeometric analysis[END_REF][START_REF] Verhoosel | Imagebased goal-oriented adaptive isogeometric analysis with application to the micromechanical modeling of trabecular bone[END_REF], composite materials [START_REF] Wang | Consistent immersed volumetric Nitsche methods for composite analysis[END_REF], cellular materials [START_REF] Rouwane | Adjusting fictitious domain parameters for fairly priced image-based modeling: Application to the regularization of digital image correlation[END_REF], microarchitectured structures obtained with additive manufacturing [START_REF] Korshunova | Image-based material characterization of complex microarchitectured additively manufactured structures[END_REF].

This chapter starts with a description of the immersed boundary method, its funda- (2015). mentals and challenges. Specific attention is paid to the integration procedure and to addressing the issue of possible ill-conditioning of the linear system. When dealing with complex models, such as scan-based models of composite materials, coupling of two or more domains is necessary. In the context of the immersed boundary approach, this becomes more complicated because the physical boundaries can go through elements. Thus, the coupling should be performed on trimmed elements. A reasonable response to this challenge is weak coupling strategies. The two of them, Mortar (Lagrange multipliers) and Nitsche's, are discussed in this chapter after the immersed boundary method. Even though in this work, we consider only linear elasticity, the approach is not limited to this case and can be applied to a wide range of problems.

Fundamentals and Challenges

Problem description

The idea behind the immersed boundary approach is to avoid the cumbersome process of meshing geometrically complex objects by relaxing the requirement for meshes to be boundary conformal. Thus, simple unfitted structured meshes can be used for the interpolation of the mechanical fields. This concept is shown on Figure 1.5: starting with the initial domain Ω phys , an additional fictitious part Ω f ict is introduced, such that the final domain Ω is a simple rectangular domain comprising both Ω phys and Ω f ict . We consider domain Ω phys . Its boundary is denoted as Γ phys . Displacements u g are imposed on boundary Γ u g , and surface forces F g are applied on boundary Γ F g . We denote f g the body forces that act on domain Ω phys (see Figure 1.5 (left)). This notation for boundary conditions is valid throughout this work. We suppose linear elastic behaviour and small deformations for the sake of simplicity.

Following the immersed boundary approach, domain Ω phys is embedded in a simple rectangular domain Ω. Now, a simple Cartesian mesh can be built on the entire domain Ω (see again Figure 1.5(right)).

We consider the classical linear elasticity problem. To present it in a weak form, we start by defining the functional spaces U i and V i that contain the displacement solution and test functions respectively:

U = {u ∈ [H 1 (Ω i )] d , u| Γu = u g }; V = {v ∈ [H 1 (Ω i )] d , v| Γu = 0}.
The well-known weak form reads: find u ∈ U such that,

a(u, v) = l(v), ∀v ∈ V,
where a(u, v) and l(v) are standard bilinear and linear forms:

a(u, v) = Ω phys ε(v) : Cε(u)dΩ, l(v) = Ω phys v • f g dΩ + Γ F v • F g dΓ F .
σ denotes the Cauchy stress tensor, ε(u) the infinitesimal strain tensor, C the Hooke tensor. Note that even if in the equation above, the integrals are performed over the physical domain Ω phys , it may not be exactly the case. The technical aspects related to the method are discussed further.

We aim to solve the problem within an 'unfitted' setting, where Ω phys is not meshed exactly but instead represented implicitly by a scalar level-set function φ(x):

φ(x) =        > 0, if x ∈ Ω phys , 0, if x ∈ Γ phys , < 0, if x ∈ Ω f ict .
(1.15)

Finally, all one has to do is to solve the well-known system of linear equations Ku = f , where the stiffness matrix K is correctly built taking into account the physical frontiers of the problem.

Thus, the immersed boundary approach shifts the initial meshing challenge to the integration step. In this context, three main challenges arise:

1. integration scheme: special integration rules should be applied to construct the operators accurately;

2. possible ill-conditioning of the system: elements trimmed arbitrarily may be situated in the fictitious domain except for a small part. This gives small terms in the stiffness matrix and thus leads to the ill-conditioning of the linear system;

3. application of Dirichlet boundary conditions and coupling through the immersed boundary: since physical boundary can go through elements, enforcement of Dirichlet boundary conditions in the strong form is not always possible. Similar to domain coupling, appropriate strategies should be applied to deal with boundaries which cut elements.

The upcoming sections are dedicated to addressing these challenges.

Numerical integration

One of the primary challenges in simulating irregular geometries is how to handle elements that are partially inside and partially outside the physical domain of interest.

Clearly, the standard Gauss quadrature method becomes invalid when dealing with the arbitrary shape of the cutting interface. Various strategies have been devised to address the challenge of integrating functions over arbitrarily shaped domains within the context of fictitious domain methods. Moment fitting techniques, e.g., [START_REF] Abedian | Performance of different integration schemes in facing discontinuities in the finite cell method[END_REF], [START_REF] Garhuom | Non-negative moment fitting quadrature for cut finite elements and cells undergoing large deformations[END_REF], [START_REF] Legrain | Non-negative moment fitting quadrature rules for fictitious domain methods[END_REF], have been developed to enable integration over such domains by selecting arbitrary integration points and optimizing integration weights. Numerous efforts have also been dedicated to the creation of efficient quadrature rules specifically designed for trimmed elements with arbitrary shapes and topologies, with contributions from studies like [START_REF] Kudela | Smart octrees: Accurately integrating discontinuous functions in 3d[END_REF], Legrain (2013), Nagy and Benson (2015b), [START_REF] Stavrev | Geometrically accurate, efficient, and flexible quadrature techniques for the tetrahedral finite cell method[END_REF]. A way of avoiding integrating cut mesh elements was proposed in φ-FEM, where basis functions are combined with the given level-set function [START_REF] Duprez | φ-fem: a finite element method on domains defined by level-sets[END_REF]Lozinski, 2020, Lozinski, 2016). Additionally, a novel approach presented in [START_REF] Wei | Immersed boundary-conformal isogeometric method for linear elliptic problems[END_REF] involves the decomposition of trimmed elements into reparametrization cells.

CHAPTER 1. STATE OF THE ART: ON THE MULTI-SCALE MODELLING AND COMPUTATION OF COMPLEX STRUCTURES WITH IGA

In this work, following FCM, we apply adaptive integration, that allows to approximate the physical domain [START_REF] Verhoosel | Imagebased goal-oriented adaptive isogeometric analysis with application to the micromechanical modeling of trabecular bone[END_REF]. While integrating basis functions, the position of each element concerning physical boundary Γ phys is verified. Thus, three options are possible:

1. the element considered is entirely inside physical domain Ω phys . In this case, a full Gauss integration is performed.

2. the element considered is entirely inside fictitious domain Ω f ict . In this case, there is no contribution to global operators from the element.

3. the element considered is situated both in physical domain Ω phys and in fictitious domain Ω f ict . In this case, the element (also named a cell in the FCM terminology) is divided into four sub-cells. With the use of sampled points uniformly distributed along its boundaries and values of the level-set function Eq. ( 1.15) at these points, the position of each sub-cell is evaluated. If the sub-cell intersects the material interface the algorithm subdivides it into four new sub-cells. This decomposition is repeated until a predefined maximum level is reached (denoted l max ). For the last cut sub-cell, a Delaunay tessellation is constructed (see [START_REF] Rouwane | Adjusting fictitious domain parameters for fairly priced image-based modeling: Application to the regularization of digital image correlation[END_REF], [START_REF] Verhoosel | Imagebased goal-oriented adaptive isogeometric analysis with application to the micromechanical modeling of trabecular bone[END_REF] for more details). On the contrary, if the level-set function has values greater or smaller than the threshold value on all sampling points, the sub-cell is subsequently either integrated with a full Gauss quadrature or removed from the integration procedure.

This procedure is illustrated on Figure 1.6. The initial discretization (the uniform Cartesian grid used for simulation) is shown in green. The trimming curve (the physical boundary Γ) is in red. The maximum level of decomposition equals to three in order to keep the example illustrative. Thus, in the final step, sub-cells are three times smaller than the initial cells (green elements). In this step, the physical boundary Γ is interpolated and a tessellation-based scheme is applied.

Ill-conditioning of the linear system

When the computational domain Ω is trimmed arbitrarily by a physical frontier, certain basis functions have their support both in physical and fictitious parts. It can occur, that for some of them, the support belongs mostly to the fictitious part. As a result, the operators are ill-conditioned and, therefore, solver accuracy is affected. This occurs because by reducing the active support of these basis functions, the minimum eigenvalues of the linear system matrix are also reduced, while the maximum eigenvalue remains unchanged.

There are different ways of addressing this issue. One of the ways to deal with this drawback is to associate a virtual stiffness to the fictitious part so that the elasticity tensor is multiplied by a factor α [START_REF] Rank | Geometric modeling, isogeometric analysis and the finite cell method[END_REF][START_REF] Schillinger | The finite cell method: A review in the context of higherorder structural analysis of CAD and image-based geometric models[END_REF], Schillinger et al., 2012):

σ α = ασ = αC : ε, α(x) = 1 ∀x ∈ Ω phys 10 -q ∀x ∈ Ω f ict
Thus, the contribution of the fictitious domain is penalized by α = 10 -q with typical values of q between 5 and 10 ( [START_REF] Schillinger | The finite cell method: A review in the context of higherorder structural analysis of CAD and image-based geometric models[END_REF]. In this case, a compromise between the accuracy of the integration and the conditioning of the system should be found.

Figure 1.6: Integration procedure for l max = 3: discretization of the problem is shown in green, discretization used for the integration -in black, along with the corresponding Gauss points, the trimming curve -in red.

For the case of scan-based geometry, an alternative way was proposed in [START_REF] Verhoosel | Imagebased goal-oriented adaptive isogeometric analysis with application to the micromechanical modeling of trabecular bone[END_REF]. According to this idea, there is no penalization of the stress tensor, but basis functions are removed depending on criteria that characterize their contributions to the operators.

While dealing with the integration of trimmed elements, the basis functions with negligible contribution should be removed. The following criteria (see [START_REF] Rouwane | Adjusting fictitious domain parameters for fairly priced image-based modeling: Application to the regularization of digital image correlation[END_REF], [START_REF] Verhoosel | Imagebased goal-oriented adaptive isogeometric analysis with application to the micromechanical modeling of trabecular bone[END_REF] for more details) may be used:

1. Criteria based on the basis functions value on the support:

Supp(N p i )∩Ω phys N p i (x, y)dx Supp(N p i ) N p i (x, y)dx ≤ ε (1.16)
2. Criteria based on the area fracture of the support:

A(supp(N p i ) ∩ Ω phys ) A(supp(N p i )) 1/d ≤ ε,
where

N p i is the basis function in question, A(supp(N p i ))
is the area of the support of a basis function N p i and d corresponds to the dimension of a problem. When the support of a basis function N p i is completely inside the physical domain, both fractions are equal to one, thus, the function will never be removed.

The selection of ε represents a trade-off between solution accuracy-where smaller values of ε yield more precise integration-and the conditioning of the stiffness matrix. Smaller ε values lead to smaller terms in the operator. A value of ε = 0 implies that all basis functions are considered, regardless of how minor their contributions may be.

Linear system preconditioning

As stated before, the accuracy of the integration can lead to poor conditioning of the system and therefore may be sacrificed to avoid numerical problems. However, 28 CHAPTER 1. STATE OF THE ART: ON THE MULTI-SCALE MODELLING AND COMPUTATION OF COMPLEX STRUCTURES WITH IGA the integration may be performed accurately, provided a preconditioner is applied to the linear system. In this work, we apply a diagonal scaling (Antolin et al., 2019b[START_REF] De Prenter | Condition number analysis and preconditioning of the finite cell method[END_REF]. For the linear system Ku = f it results in the following preconditioning:

PKPx = Pf , u = Px,
where P is the diagonal scaling preconditioner:

P =     1 √ K11 . . . 1 √ Knn .    
This preconditioner allows to renormalized the basis function with respect to their contribution to the stiffness matrix K. The interested reader is referred to Antolin et al. (2019b), [START_REF] Buffa | A minimal stabilization procedure for isogeometric methods on trimmed geometries[END_REF][START_REF] De Prenter | Condition number analysis and preconditioning of the finite cell method[END_REF] for a more detailed discussion on the conditioning of the stiffness matrix in trimmed geometries.

Boundary conditions

In this work, constant Dirichlet boundary conditions are imposed only on complete exterior rectilinear boundaries. In this simple case, control points describing these boundaries lie on the geometry and the applied displacement is constant, thus, Dirichlet boundary conditions are directly applied to these degrees of freedom. In the general case, this strategy is not possible and Dirichlet boundary conditions are imposed in a weak sense [START_REF] Ruess | Weakly enforced essential boundary conditions for NURBS-embedded and trimmed NURBS geometries on the basis of the finite cell method[END_REF][START_REF] Schillinger | The finite cell method: A review in the context of higherorder structural analysis of CAD and image-based geometric models[END_REF].

Weak coupling strategies

Complex geometric models often comprise several patches, naturally raising the challenge of patch coupling. Each patch is typically defined independently, leading to non-conformal interfaces between them. This complexity is further compounded when physical boundaries intersect trimmed elements within a domain (as an example, see Figure 1.3(c)).

To address the issue of coupling in IGA, established strategies from FEA have been adapted and can be categorized into three main approaches:

• Mortar Method: The Mortar method seeks to enforce continuity and compatibility conditions across non-conformal interfaces by introducing Lagrange multipliers.

• Penalty Method: The Penalty method introduces penalty terms into the governing equations to penalize deviations from the continuity and compatibility conditions at the interface.

• Nitsche's Method: Nitsche's method reformulates the problem by adding stabilization terms to the weak form of the governing equations. It provides a more robust approach for handling non-conformal interfaces and avoids the need for Lagrange multipliers or penalty parameters.

In the following sections, we present Mortar and Nitsche's approaches, leaving penalty coupling behind. This is due to the fact that the current work utilizes these two methods and not the penalty one.

Problem definition

In the present section, we consider two non-overlapping subdomains, Ω 1 and Ω 2 , which are to be joined along an interface Γ, as in Figure 1.7. For i = 1, 2, Ω i ∈ R d with d = 2 or 3 being the dimension of the problem. We assume that Ω i is subjected to body forces f g i , surface force F g i is applied on the boundary Γ F g i and displacement u g i is imposed over boundary Γ u g i (see again Figure 1.7). We assume linear elastic behaviour and small deformations.

In each subdomain, the equilibrium equations, the kinematic constraints, and the constitutive relations should be verified. The corresponding governing equations read; for i ∈ {1, 2} :

div(σ i ) + f g i = 0 in Ω i ; σ i n i = F g i on Γ F g i ; u i = u g i on Γ u g i ; σ i = C i ε(u i ) in Ω i .
(1.17)

Here ε(u i ) denotes the infinitesimal strain tensors, σ i the Cauchy stress tensors, C i the Hooke tensors, n i the outward unit normals.

For the sake of simplicity, here, we consider only the perfect interface. Thus, kinematic compatibility and equilibrium of the tractions should be added to the coupling interface, which reads as follows:

u 1 -u 2 = 0 on Γ;
(1.18a)

σ 1 n 1 + σ 2 n 2 = 0 on Γ. (1.18b)
Combining equations inside the subdomains (1.17) with the interface conditions (1.18) is equivalent to solving the linear elasticity problem inside the whole domain Ω = Ω 1 ∪ Ω 2 ∪ Γ. The interface conditions (1.18) can not be imposed in a strong way when the interface is not conformal. Thus, these conditions should be imposed in a weak form.

The Mortar Method

The Mortar Method, also known as Lagrange Multiplier Method, is a numerical technique used to handle non-conformal interfaces or contact conditions between different 30 CHAPTER 1. STATE OF THE ART: ON THE MULTI-SCALE MODELLING AND COMPUTATION OF COMPLEX STRUCTURES WITH IGA subdomains. Mortar-based approaches have been successfully applied within the framework of IGA (e.g., [START_REF] Bouclier | Development of a new, more regular, Mortar method for the coupling of NURBS subdomains within a NURBS patch: Application to a non-intrusive local enrichment of NURBS patches[END_REF], [START_REF] Brivadis | Isogeometric Mortar methods[END_REF], [START_REF] Dornisch | The weak substitution method-an application of the Mortar method for patch coupling in NURBS-based isogeometric analysis[END_REF], [START_REF] Miao | Isogeometric Bézier dual mortaring: The enriched Bézier dual basis with application to second-and fourthorder problems[END_REF], [START_REF] Wunderlich | Biorthogonal splines for optimal weak patch-coupling in isogeometric analysis with applications to finite deformation elasticity[END_REF] to name a few).

Here we also define the functional spaces U i and V i that contain the displacement solution and test functions respectively for i ∈ {1, 2}:

U i = {u i ∈ [H 1 (Ω i )] d , u m | Γu i = u g i }; V i = {v i ∈ [H 1 (Ω i )] d , v m | Γu i = 0}
. The idea of the Mortar Method is to introduce an additional interface field which is the Lagrange multiplier field. A Lagrange multiplier λ ∈ M (where M is an appropriate space) represents the interface traction forces in (1.18b), -σ 1 n 1 = σ 2 n 2 = λ. These are additional unknowns introduced at the interface or contact regions. The compatibility condition (1.18a) is imposed weakly over Γ. Lagrange multipliers are solved along with the primary unknowns of the problem.

The Lagrangian of the coupled problem is:

L ((u 1 , u 2 ), λ) = 1 2 a 1 (u 1 , u 1 ) + 1 2 a 2 (u 2 , u 2 ) -l 1 (u 1 ) -l 2 (u 2 ) + b(λ, u 1 -u 2 ), (1.19)
where for i ∈ {1, 2}:

a i (v i , u i ) = Ωi ε(v i ) : C i ε(u i )dΩ i and l i (v i ) = Ωi v i • f g i dΩ i + Γ F i v i • F g m dΓ Fi .
The bilinear form b is defined as an integral over the interface Γ :

b(µ, u) = Γ µ • udΓ.
With Eq. (1.19), we can obtain the following formulation of the coupled problem:

find u 1 ∈ U 1 , u 2 ∈ U 2 , λ ∈ M such that:        a 1 (u 1 , v 1 ) + b(λ, v 1 ) = l 1 (v 1 ), ∀v 1 ∈ V 1 ; a 2 (u 2 , v 2 ) -b(λ, v 2 ) = l 2 (v 2 ), ∀v 2 ∈ V 2 ; b(µ, u 1 -u 2 ) = 0, ∀µ ∈ M.
Thus, the Mortar Method is implemented by solving the augmented problem, which includes the primary unknowns (displacements) and the Lagrange multipliers. This results in the following discrete version in the linear equations system:

  K 1 0 C T 1 0 K 2 -C T 2 C 1 -C 2 0      u 1 u 2 λ    =    f 1 f 2 0    .
In the equations above, K 1 (resp. K 2 ) and f 1 (resp. f 2 ) are the classical stiffness matrix and the vector forces associated to domain Ω 1 (resp. Ω 2 ). For i ∈ {1, 2} C i is the Mortar coupling operator and is defined as follows:

C i = Γ N λ N T i dΓ,
where N i (resp. N λ ) are standard shape function matrices related to subspaces U i (resp.M). The compelling advantage of this approach is that the stiffness operators of the subdomains stay unmodified. Indeed, the subdomains communicate with each other only by means of the Lagrange multiplier.

While the Mortar Method is a powerful technique for handling non-conforming interfaces and contact conditions, there are some considerations to be made:

1. mesh locking: when coupling non-conformal meshes with substantial discrepancies in material parameters through the Mortar coupling method, the phenomenon of mesh locking, as discussed in [START_REF] Sanders | A Nitsche embedded mesh method[END_REF], becomes a notable concern. Mesh locking arises when the number of constraints imposed on the interface exceeds the available degrees of freedom, rendering these constraints unable to fully accommodate the physical deformations that the elements should permit. This situation is particularly pronounced when coupling a domain which is both finer and stiffer with a coarser grid. Concurrently, oscillations may manifest in the Lagrange multipliers, utilized to enforce these constraints. The consequences of mesh locking are twofold: a notable loss of overall solution accuracy and a less-than-optimal convergence rate in affected problems. Conversely, an insufficient imposition of constraints may lead to the underrepresentation of compatibility conditions and inaccurate modelling of surface fluxes. Hence, special care should be paid while performing the Mortar approach on a non-conformal interface.

2. suboptimal convergence rate for high order basis functions: the Lagrange multiplier approach appears to fail to reproduce the optimal convergence rates, especially when increasing the polynomial order of the IsoGeometric shape functions [START_REF] Bouclier | A Nitsche-based non-intrusive coupling strategy for global/local isogeometric structural analysis[END_REF]. This may have origin in the non-conformal nature of coupling and, thus, a complicated choice of the appropriate Lagrange multipliers field.

3. choice of Lagrange multipliers field and computational cost: for severely non-conformal meshes the question of the choice of Lagrange multipliers naturally arises. Special attention should be paid to this question since the performance of the method highly depends on the compatibility between approximation spaces of the subdomains and the interface.

Nitsche's method

Nitsche's method is an alternative to the Mortar Method for handling contact, interface, or Dirichlet boundary conditions in numerical simulations. It was originally proposed in [START_REF] Nitsche | Über ein variationsprinzip zur lösung von Dirichlet-problemen bei verwendung von teilräumen, die keinen randbedingungen unterworfen sind[END_REF] for week reinforcement of essential boundary conditions in elliptic PDEs. It has since become an essential technique in numerical simulations, particularly for problems involving contact, fluid-structure interaction, and various types of interfaces. In recent years, the method has been successfully applied to different problems. It has been used for modelling contact and friction problems (see [START_REF] Chouly | An overview of recent results on Nitsche's method for contact problems[END_REF] for a review on this topic) and for thermomechanical contact [START_REF] Seitz | Nitsche's method for finite deformation thermomechanical contact problems[END_REF].

It was combined with CutFEM to model the interaction of the fluid with poroelastic solids [START_REF] Ager | A Nitsche-based cut finite element method for the coupling of incompressible fluid flow with poroelasticity[END_REF], and with XFEM for modelling of the interaction between incompressible viscous fluid and elastic structures [START_REF] Alauzet | Nitsche-XFEM for the coupling of an incompressible fluid with immersed thin-walled structures[END_REF]. The core idea behind Nitsche's method lies in the weak enforcement of boundary or interface conditions. Unlike the Mortar method, which introduces additional degrees of freedom at the interface, Nitsche's method embeds the continuity and equilibrium conditions directly into the initial variational formulation of the problem.

Nitsche's method starts with the standard variational formulation of the PDE problem. To account for the coupling, both the interface conditions (1.18a) and (1.18b) are imposed weakly and a penalty-like stabilization term is added to ensure the ellipticity of the boundary value problem. The corresponding weak formulation can be written

: find (u 1 , u 2 ) ∈ U 1 × U 2 such that, i∈{1,2} a i (u i , v i ) - Γ {Cε(v) n} • [[u]]dΓ - Γ [[v]] • {σ n}dΓ +ζ Γ [[v]] • [[u]]dΓ = i∈{1,2} l i (v i ), ∀v i ∈ V i ; (1.20)
where a i is the standard bilinear form l i is the linear form associated to subdomain Ω i , [[u]] = u 1 -u 2 and {σ n} denote the displacement jump and the stress flux across Γ, respectively, and ζ is a stabilization parameter that shall depend on the mesh discretizations and material properties of the two subdomains to be coupled. Following [START_REF] Antolin | Overlapping multipatch isogeometric method with minimal stabilization[END_REF], [START_REF] Wei | Immersed boundary-conformal isogeometric method for linear elliptic problems[END_REF], we adopt the one-sided flux, i.e. that we consider only the term for the conformal subdomain:

{σn} = σ 1 n 1 = C 1 ε(u 1 )n 1 .
As for the stabilization parameter ζ, which is necessary to ensure the ellipticity of the boundary value problem, there exist different expressions (see, e.g., [START_REF] Hansbo | Nitsche's method for interface problems in computational mechanics[END_REF], [START_REF] Johansson | Multimesh finite element methods: Solving PDEs on multiple intersecting meshes[END_REF], [START_REF] Wei | Immersed boundary-conformal isogeometric method for linear elliptic problems[END_REF]). In de Prenter et al. ( 2018) the stabilization parameter is discussed within the framework of an unfitted mesh. It is worth noting, that even though the last term in Eq. (1.20) remind us of the Penalty coupling approach, the choice of the parameter ζ exerts less influence on the solution compared to the penalty parameter in the Penalty coupling approach. This term is essential to guarantee problem coercivity but not to enforce the interface conditions. Finally, let us recast formulation (1.20) by writing it as follows:

find u ∈ U such that, a(u, v) = l(v), ∀v ∈ V, (1.21)
where and a(u, v) and l(v) are the left-hand side and righthand side of Eq. (1.20), respectively. We recall that Eq. (1.21) enforces (1.18)-(1.17).

U = U 1 × U 2 , V = V 1 × V 2 ,
In the discrete form, Eq. (1.21) gives:

K 1 + ζN 11 -2(K N 21 ) T -ζN N 12 -2K N 21 -ζ(N N 12 ) T K 2 -2K N 22 -2(K N 22 ) T + N N 22 u 1 u 2 = f 1 f 2 ,
where K 1 (resp. K 2 ) and f 1 (resp. f 2 ) are the classical stiffness matrix and the vector forces associated to domain Ω 1 (resp. Ω 1 ). The operators K N ii read:

K N 11 = Γ - 1 2 (n 1 D 1 B 1 ) T N 1 dΓ ; K N 12 = Γ 1 2 (n 1 D 1 B 1 ) T N 2 dΓ ; K N 21 = Γ - 1 2 (n 1 D 2 B 2 ) T N 1 dΓ ; K N 22 = Γ 1 2 (n 1 D 2 B 2 ) T N 2 dΓ .
N N ii are defined as follows:

N N 11 = Γ (N 1 ) T N 1 dΓ ; N N 12 = Γ (N 1 ) T N 2 dΓ ; N N 22 = Γ (N 2 ) T N 2 dΓ ; 1.3. GLOBAL/LOCAL COUPLING 33
In the above equations, N 1 and N 2 are the standard shape function matrices, B 1 and B 2 are the standard strain-displacement operators and D 1 and D 2 represent Hooke constitutive law in Voigt notation.

Nitsche's approach has several compelling advantages. The coupled stiffness matrix is positively defined and the formulation is variationally consistent. Furthermore, unlike its Mortar counterpart, Nitsche's method does not introduce additional fields. Thus, we avoid choosing an approximation space for Lagrange multiplies.

Remark 1.2.1

There is a non-symmetric variant of Nitsche's approach available. This formulation maintains variational stability and eliminates the need for any additional stabilization terms. But as its names suggests, it breaks the symmetry of the stiffness operator. For more details on this topic, the interested reader is referred to [START_REF] Bouclier | A Nitsche-based non-intrusive coupling strategy for global/local isogeometric structural analysis[END_REF], [START_REF] Chouly | Symmetric and non-symmetric variants of Nitsche's method for contact problems in elasticity: theory and numerical experiments[END_REF], Schillinger et al. (2016a).

Remark 1.2.2

Even though in this work we use Mortar and Nitsche's methods to perform coupling between subdomains, both of them can be used to apply the Dirichlet boundary conditions. In the context of the immersed boundary approach, this is the way to ensure Dirichlet boundary conditions, since it is not always possible to apply them in the strong form.

Global/Local Coupling for Modelling Localized

Non-Linearities

Introduction

Most composite materials inherently possess a multi-scale nature, where the scale of the individual constituents is of a considerably lower order than that of the resulting material and structure. This inherent complexity necessitates a careful consideration of how to analyse these materials effectively.

Traditionally, many analyses of composite structures resort to the use of effective or homogenized material properties. Rather than accounting for the individual component properties and geometrical arrangements, these analyses simplify the material behaviour. However, sometimes such analyses are not accurate enough. Ideally, one could turn to the microscopic scale for a more precise understanding of composite behaviour. However, in practice, microscopic models are often overly complex and computationally intensive, making them impractical for analysing large structures.

A compelling solution to this issue is the concept of global/local modelling. It involves the coupling of macroscopic and microscopic models to harness the efficiency of the former and the precision of the latter. This approach allows for a more comprehensive understanding of composite materials, bridging the gap between the macro and micro worlds of mechanics. Global models typically make simplifications in geometry, kinematics, and material properties. In contrast, local models tackle greater complexity, with their primary objective being the accurate representation of structural behaviour.

CHAPTER 1. STATE OF THE ART: ON THE MULTI-SCALE MODELLING AND COMPUTATION OF COMPLEX STRUCTURES WITH IGA

Global and local models should be correctly coupled to provide accurate results. The most commonly used method is submodelling. It has been used in different contexts and is often available with commercial FEA software. The global model, encompassing the whole structure, has a coarse mesh, while the local one describes only a zone of interest with a finer discretization. According to this approach, first, the solution for the global model is computed. Then, it is used to impose boundary conditions of the local model, in terms of displacement [START_REF] Voleti | Global-local analysis of large-scale composite structures using finite element methods[END_REF] or stress [START_REF] Jara-Almonte | The specified boundary stiffness/force sbsf method for finite element subregion analysis[END_REF]. Its compelling advantage is the possibility of avoiding detailed mesh on the whole structure while accurately modelling the local zone of interest. However, there is no real dialogue between global and local models (there is no information from the local model communicated to the global one). As a result, this approach is applicable only if the local model has no real influence on the global one.

When submodelling is not sufficient, a strong coupling should be applied. In this case, the influence of localized nonlinearities on the global response is taken into account. To do so, the models can be coupled with different techniques (such as Mortar and Nitsche's methods, discussed in 1.2.3). They result in a coupled linear system, which is solved using a direct solver. Even if this allows us to obtain accurate solutions, there are significant drawbacks to this technique. First, the size of the coupled linear system can be large, especially if several detailed local models are present. Second, even though the non-linearities are present only locally, the whole system should be solved with a non-linear solver, which is far from optimal. Moreover, the global model should be modified specifically for the considered local models. Finally, when the local model changes in size, for example, because of the growth of a crack or a damaged area, it's necessary to update both the local and global models in the case of strong coupling. Changing geometry, and, consequently, mesh, is a time-consuming process.

As an answer to the above-mentioned issues, a new class of method has emerged: the non-intrusive coupling [START_REF] Duval | Non-intrusive coupling: recent advances and scalable nonlinear domain decomposition[END_REF][START_REF] Gendre | Non-intrusive and exact global/local techniques for structural problems with local plasticity[END_REF][START_REF] Whitcomb | Iterative global/local finite element analysis[END_REF]. Within this method, the existing numerical model can be modified, without actually modifying its stiffness matrix. The non-intrusive global/local coupling has recently been applied to different problems (see [START_REF] Gendre | Non-intrusive and exact global/local techniques for structural problems with local plasticity[END_REF] 2020) for fracture modelling with the phase-field approach, [START_REF] Guinard | Multiscale analysis of complex aeronautical structures using robust non-intrusive coupling[END_REF] for real aeronautical structures, and [START_REF] Wangermez | Interface coupling method for the global-local analysis of heterogeneous models: A second-order homogenization-based strategy[END_REF] for multi-scale periodic heterogeneous materials, [START_REF] Guguin | Nonintrusive coupling of 3d and 2d laminated composite models based on finite element 3d recovery[END_REF] for 2D/3D coupling, to name a few). In addition, within the non-intrusive coupling, research codes can be seamlessly integrated into commercial software, since there is no need for software modifications.

In this section, we introduce the non-invasive global/local coupling approach through the analysis of a simple problem involving the presence of a hole within a global structure (see Figure 1.8 for illustration). It's important to note that while we currently discuss this within the context of composite materials, the non-invasive coupling technique holds value for modelling various multi-scale structures, such as aeronautical structures or subsea drilling equipment. For more information, please refer to [START_REF] Bouclier | IGA: Non-Invasive Coupling with FEM and Regularization of Digital Image Correlation Problems[END_REF].

Non-Invasive Approach

In the following sections, we will outline the development of the non-invasive coupling strategy. We start with the formulation of the reference global/local coupling problem. Then, the non-invasive algorithm is presented as an alternative to the conventional monolithic resolution approach.

Governing equations

We start by introducing the mechanical multi-scale global/local problem that we seek to solve. We consider a global (coarse) model of a structure that is characterized by a physical domain Ω 1 (see Figure 1.8(a)). This domain is divided into two disjoint, open and bounded subsets Ω 11 and Ω 12 . These two non-overlapping subdomains share a common interface denoted by Γ such that Ω 1 = Ω 11 ∪ Ω 12 ∪ Γ and Ω 11 ∩ Ω 12 = ∅. The global structure is considered linear elastic. We assume that such behaviour and the coarse discretization are sufficient to accurately capture the solution except in the small region Ω 12 where a local (possibly non-smooth, singular, or even discontinuous) phenomenon is to be introduced. As a consequence, a local (more detailed) model characterized by domain Ω 2 is constructed to replace the global model in Ω 12 (see Figure 1.8(a)). The substitution of the local model within the global one is achieved through interface Γ. The resulting global/local problem to be solved is a classical multi-domain problem in Ω 11 ∪ Ω 2 ∪ Γ. Hence, we consider here the same coupling problem as in Section 1.2.3.

Subdomain Ω 11 (resp. Ω 2 ) is subjected to body force f g 11 (resp. f g 2 ). Surfaces forces

F g 11 (resp. F g 2 ) are applied on boundary Γ F g 11 (resp. Γ F g 2 ).
Additionally, displacements u g 11 (resp.u g 2 ) are imposed on boundary Γ u g 11 (resp. Γ u g 2 ) (see again Figure 1.8(a)). For the sake of simplicity, we consider the local model to be linear elastic. However, the non-invasive global/local coupling method can be applied for any (possibly nonlinear) local behaviour. Interface Γ is a perfect interface.

Thus, following Section 1.2.3, in each subdomain, the equilibrium equations, the kinematic constraints, and the constitutive relations should be verified. Therefore, for i ∈ {11, 2} Eqs. (1.17

u 11 -u 2 = 0 on Γ; σ 11 n 11 + σ 2 n 2 = 0 on Γ.
Conversely to Section 1.2.3, where the coupled problem was solved directly, here we will present an iterative strategy of solving this problem.

Monolithic resolution

We start by formulating the coupling with a Lagrange multipliers approach. We don't care about the non-invasive formulation yet and we apply the Mortar coupling presented in Section 1.2.3.

The formulation of the coupled problem being introduced in Section 1.2.3, here we present only its discrete version. Following the notation of Section 1.2.3, we denote the matrices of shape functions associated with the model over Ω 11 and local model over Ω 2 by N 11 and N 2 , respectively. Matrix N λ is the matrix of shape functions associated with interface Γ. The displacement vectors u 11 , u 2 corresponds to Ω 11 and Ω 2 , respectively. The Lagrange multipliers vector is denoted by λ. With this notation, we obtain the following linear equations system:

  K 11 0 C T 11 0 K 2 -C T 2 C 11 -C 2 0      u 11 u 2 λ    =    f 11 f 2 0    . (1.25)
Following the previous notation, K 11 (respectively f 11 ) and K 2 (resp. f 2 ) are the classical stiffness matrices (resp. load vectors) associated with subdomains Ω 11 and Ω 2 . The operators C 11 and C 2 are the Mortar coupling operators.

The resolution of the global/local problem (1.25) constitutes the classical monolithic approach: the coupled model of Figure 1.8(b) is computed directly using a single direct solver. This strategy is invasive in the sense that it requires (i) to modify the initial global model (and thus its operators) to remove some of its elements, or possibly pieces of elements (see Figure 1.8 again) which may lead to ill-conditioned stiffness operators [START_REF] De Prenter | Condition number analysis and preconditioning of the finite cell method[END_REF], Wei et al., 2021), and (ii) to set up an additional solver that merges the contributions of the two models. In case the local detail grows up (during crack propagation, or expansion of damage or plasticity for instance), the situation gets even worse since not only the local operator K 2 but also the global operator K 11 have to be fully re-built, and the augmented system (1.25) re-factorized during the simulation.

Iterative resolution

In contrast to directly solving system (1.25), the non-invasive strategy adopts an iterative exchange approach that alternates between solving global problems over Ω 1 and local problems over Ω 2 . The derivation of this strategy involves two steps. Firstly, we split the original system (1.25) to identify two different in terms of boundary conditions problems: Neumann problem on Ω 11 , and Dirichlet problem on Ω 2 . Then, we extend the global solution from Ω 11 to Ω 12 . This allows us to apply the additivity of the integral with respect to domain Ω 1 = Ω 11 ∪ Ω 12 ∪ Γ to recover the initial whole global model.

Splitting of the coupled problem.

The idea is to split system (1.25) into two and apply an iterative algorithm. We obtain the following asymmetric algorithm in the sense that Dirichlet and Neumann problems with respect to Γ are alternatively solved until convergence. More precisely, for the nth iteration, starting with initial guesses λ (0) and u (0) 11 , we look for u

(n) 11 , u (n) 2
and λ (n) such that:

1. Resolution of a Neumann problem (with respect to Γ) over Ω 11 :

K 11 u (n) 11 = f 11 -C 11
T λ (n-1) .

(1.26)

2. Resolution of a Dirichlet problem (with respect to Γ) over Ω 2 :

  K 2 -C 2 T -C 2 0     u 2 (n) λ (n)   =   f 2 -C 11 u (n) 11   . ( 1
K 1 u 1 = K 11 u 1 + K 2 u 1 ,
where K 1 and u 1 are the stiffness matrix and the displacement vector associated to domain Ω 1 , K 11 and K 2 are prolongations of operators K 11 and K 2 into domain Ω 1 .

To obtain K 11 and K 2 , K 11 and K 2 are extended by inserting zeros to ensure they have the same dimensions as u 1 . In the same manner, the load vector associated to domain Ω 1 is:

f 1 = f 11 + f 2 ,
with f 11 and f 2 extensions of f 11 and f 2 to Ω 1 . Finally, we obtain the following algorithm: for the nth iteration, starting with initial guesses λ (0) and u (0)

1 , we look for u

(n) 1 , u (n) 2
and λ (n) such that:

1. Resolution of a Neumann problem (with respect to Γ) over Ω 1 :

K 1 u (n) 1 = f 1 -C 1 T λ (n-1) + λ (n-1) 12
.

(1.28)

2. Resolution of a Dirichlet problem (with respect to Γ) over Ω 2 : The procedure can be interpreted as a fixed point strategy aiming at ensuring the equilibrium of the interface reaction forces (see Eq. (1.28)) provided that the displacement is transferred at each iteration between the two models (see Eq. (1.29)). An illustration of the algorithm is provided in Fig. 1.9. For more information, we advise the interested reader to consult the following reviews on the subject Allix and Gosselet (2020), [START_REF] Duval | Non-intrusive coupling: recent advances and scalable nonlinear domain decomposition[END_REF]. Since the initial global model is now unmodified, its stiffness operator remains well-conditioned regardless of the shape of the local region, and it can be assembled and factorized only once during the pre-processing step (Bouclier andPassieux, 2018, Bouclier et al., 2016).

  K 2 -C 2 T -C 2 0     u 2 (n) λ (n)   =   f 2 -C 1 u (n) 1   . ( 1

Remark 1.3.1

We recall that the method is by no means limited to a linear elastic local model. Indeed, as long as we can apply Dirichlet boundary conditions to the local problem and compute (directly or in a post-processing step) the corresponding reaction forces, any local behaviour can be considered, as demonstrated in the large literature on the topic (see [START_REF] Duval | Non-intrusive coupling: recent advances and scalable nonlinear domain decomposition[END_REF]

Large Time Increment Method for Addressing

Non-Linear Interfaces

Motivation

In previous sections, different coupling approaches for perfect interfaces were discussed. However, in some problems, the overall mechanical structure is composed of linear subdomains interconnected by potentially non-linear interfaces. On these interfaces, behaviours like sliding or decohesion can happen. In this context, the main challenge lies in the accurate simulation of mechanical problems involving this interface behaviour. Addressing this challenge is at the core of this section. To accurately represent local interface behaviour, we propose the use of the Large Time Incremental (LaTIn) method, a mixed domain decomposition approach. This method efficiently captures local behaviours and simplifies the parallelization of linear system solving processes. We propose to utilize the LaTIn method within the framework of IGA. Thus, not only we reduce the number of degrees of freedom due to IGA, but we also allow the parallelization of the problem. This strategy addresses the issue of simulating models associated with big data.

In this section, we present the LaTIn method in a general context without specifying interface behaviour. Detailed implementations of specific cases, such as perfect interfaces, frictionless contact, Coulomb friction, and decohesion, are discussed in Chapter 3.

The LaTIn method was proposed in [START_REF] Ladevèze | Sur une famille d'algorithmes en mécanique des structures[END_REF], [START_REF] Ladevèze | Nonlinear computational structural mechanics: new approaches and nonincremental methods of calculation[END_REF] to model problems, where local non-linear equations can be separated from global linear ones. Among its first applications, one can name such non-linear problems as plasticity and (thermo-) viscoplasticity [START_REF] Boisse | A large time increment method for elastoplastic problems[END_REF][START_REF] Whitcomb | Iterative global/local finite element analysis[END_REF][START_REF] Cognard | A large time increment approach for cyclic viscoplasticity[END_REF], Cognard et al., 1999). Over the past three decades, the method has witnessed continuous development. The growing research body, that emerged over the past three decades, highlights the versatility of the LaTIn method.

There are the following examples of its applications. It was applied to model tensegrity structures [START_REF] Alart | A scalable multiscale LaTIn method adapted to nonsmooth discrete media[END_REF], to diverse frictional contact problems [START_REF] Giacoma | Toward an optimal a priori reduced basis strategy for frictional contact problems with LaTIn solver[END_REF], also with complex loading [START_REF] Blanzé | Contact problems in the design of a superconducting quadrupole prototype[END_REF]. In [START_REF] Oumaziz | A parallel noninvasive multiscale strategy for a mixed domain decomposition method with frictional contact[END_REF], authors proposed a non-invasive formulation and applied it to 3D structures involving many interfaces with frictional contact. The quasi-static formulation was extended to dynamics: to model heterogeneous materials [START_REF] Odièvre | A parallel, multiscale domain decomposition method for the transient dynamic analysis of assemblies with friction[END_REF], complex 3D structural assemblies with frictional contact [START_REF] Boucard | A parallel and multiscale strategy for the parametric study of transient dynamic problems with friction[END_REF].

The idea of separating equations found its application not only for non-linear mechanical problems but also in the case of multi-physical problems (see e.g., [START_REF] Bibliography | A LaTIn computational strategy for multiphysics problems: application to poroelasticity[END_REF] for a coupled fluid-solid problem in the context of poroelasticity). Industrial applications list assemblies modelling [START_REF] Champaney | Modular analysis of assemblages of three-dimensional structures with unilateral contact conditions[END_REF][START_REF] Oumaziz | A parallel non-invasive mixed domain decomposition-implementation and applications to mechanical assemblies[END_REF]. In the FEA context, the LaTIn method was combined with advanced approaches. Along with X-FEM, the LaTIn method was used for modelling crack propagation [START_REF] Guidault | A two-scale approach with homogenization for the computation of cracked structures[END_REF] while taking into account frictional contact between the crack faces [START_REF] Ribeaucourt | A new fatigue frictional contact crack propagation model with the coupled X-FEM/LaTIn method[END_REF]. Together with CutFEM, it was used for frictionless contact computation between a femur and a hip bone in image-based models [START_REF] Claus | Contact modeling from images using cut finite element solvers[END_REF], braided composite material [START_REF] Claus | A stable and optimally convergent LaTIn-CutFEM algorithm for multiple unilateral contact problems[END_REF], composite materials with quasi-1D fibrous reinforcements [START_REF] Kerfriden | A mixed-dimensional CutFEM methodology for the simulation of fibre-reinforced composites[END_REF], to model delamination (Allix andLadevèze, 1992, Guinard et al., 2002), delamination combined with buckling [START_REF] Saavedra | On a multiscale strategy and its optimization for the simulation of combined delamination and buckling[END_REF], and to predict the damage in [START_REF] Aubard | Damage modeling at two scales for 4d carbon/carbon composites[END_REF].

However, hitherto and to the author's knowledge, the LaTIn method has not been employed within the context of the IGA. The current work introduces a novel approach by combining the LaTIn method together with immersed boundary techniques in the framework of IGA.

General idea

Consider a model decomposed in subdomains connected by an interface. We suppose the linear elastic behaviour inside the subdomains, whilst the interface can exhibit non-linear behaviour. In the context of (quasi-) static mechanical problems, the LaTIn method relies on the two following principles:

1. The equations are separated into two groups: local non-linear equations that describe the possibly non-linear behaviour on the interface and global linear equations relative to a subdomain. This results in two spaces of partial solutions: space A of the solutions that verify the equilibrium of subdomains, space L of the solutions that verify the interface behaviour.

2. A two-stage iterative approach is applied, where the first stage is the local nonlinear and the second one is the global linear. The scheme is characterized by two parameters called search directions: k + to go from A to L, k -to go from L to A. At each iteration, two stages are performed consecutively:

(a) Local non-linear stage: for a known solution S (n) ∈ A find a solution Ŝ(n) ∈ L following the search direction k + .

(b) Global linear stage: for a known solution Ŝ(n) ∈ L find a solution S (n+1) ∈ A following the search direction k -. At the end of this step, a relaxation is applied to ensure convergence.

This idea is schematically illustrated in Figure 1.10.

General LaTIn numerical scheme

Governing equations

Once again, we undertake to solve a multi-domain problem characterized by a physical domain Ω ⊂ R d , d = 2 or 3 being the topological dimension of the domain. More precisely, since our target applications are composite materials such as fibre-reinforced ones, we consider the coupling of a global model (say the matrix) with multiple local models (say the fibres) through non-linear interfaces. For the sake of simplicity, let us take here only two subdomains. Thus, this problem corresponds to the one considered in Section 1.2.3, where two elastic domains were coupled by a perfect interface (see Figure 1.7). However, in this case, the interface behaviour may be more complicated than simply perfect (it may be non-linear, see Chapter 3). Thus, the equations corresponding to subdomains Ω 1 and Ω 2 are exactly as in Section 1.2.3 (Eqs. (1.17)). However, the interface condition should take into account possibly non-linear behaviour. From a general point of view, the latter can be written as:

σ 1 n 1 + σ 2 n 2 = 0 on Γ ;
(1.31a)

g(u 1 , u 2 , σ 1 n 1 , σ 1 n 1 ) = 0 on Γ ; (1.31b)
where Eq. (1.31a) ensures the equilibrium of the two subdomains along Γ (the same as for a perfect interface, Eq. (1.18b)), and (1.31b) expresses the possible non-linear law between the interface displacements and tractions through a given function g. In the Chapter 3, g will be detailed considering the case of contact and delamination.

Separation of the equations

To formulate the method, let us start by introducing the following interface fields for subdomains Ω i : ∀i ∈ {1, 2}, λ i represents a surface force applied to Ω i through boundary Γ, and w i is the trace of u i over Γ. With these new variables, we can define the first partial solution spaces A i that group the solutions of the linear equations set on the subdomains: ∀i ∈ {1, 2},

A i : (λ i , w i ) such that         
(1.17)

σ l i n l i = λ i on Γ w i = u l i on Γ
.

(1.32)

The spaces A i are affine spaces often called spaces of admissible fields. We denote their union:

A = i∈{1,2} A i .
Then, let us define additional force and displacement interface fields ( λi , ŵi ) i∈{1,2} that belong to the second partial solution space L, i.e. that verify the behaviour of the interface:

L : ( λi , ŵi ) i∈{m,f } such that    λ1 + λ2 = 0 on Γ g( ŵ1 , ŵ2 , λ1 , λ2 ) = 0 on Γ . (1.33)
L is a manifold referred to as the local space since it is defined by pointwise independent equations. Finally, it is required to add the search directions k + and k -to communicate between spaces A and L and thus close the problem. This is performed in a mixed way as follows: ∀i ∈ {1, 2},

k + : λi -λ i -k + i ( ŵi -w i ) = 0 on Γ ; (1.34a) k -: λi -λ i + k - i ( ŵi -w i ) = 0 on Γ ; (1.34b)
where k + i and k - i represent interface stiffnesses. In practice, we choose k + i = k - i = k i > 0 which is the classical setting of the LaTIn algorithm [START_REF] Ladevèze | Nonlinear computational structural mechanics: new approaches and nonincremental methods of calculation[END_REF].

Iterative algorithm

With the above separation of equations in hand, we can then perform a fixed point numerical scheme to solve the problem. More precisely, for the nth iteration, starting

i ) ∈ A i , ∀i ∈ {1, 2}, we subsequently perform two steps:

1. Local (non-linear) stage:

∀i ∈ {1, 2} , given S (n) i = (λ (n) i , w (n) i ) ∈ A i , find Ŝ(n) i = ( λ(n) i , ŵ(n) i ) ∈ L ∩ k + .
2. Global (linear) stage:

∀i ∈ {2, 2} , given Ŝ(n) i = ( λ(n) i , ŵ(n) i ) ∈ L, find S (n+1) i = (λ (n+1) i , w (n+1) i ) ∈ A i ∩ k -.
To ensure the convergence of the algorithm, a relaxation step is also performed at the end of the linear stage:

∀i ∈ {1, 2} , S (n) i ← θS (n+1) i + (1 -θ)S (n) i ,
where we choose θ = 0.5 following again [START_REF] Ladevèze | Nonlinear computational structural mechanics: new approaches and nonincremental methods of calculation[END_REF]. The overall procedure is illustrated in Fig. 1.10.

Finally, an indicator of error is used to quantify the distance between A and L for two successive partial solutions. It is written in an energy norm such that:

η = i∈{1,2} Γ k i w (n) i - ŵ(n-1) i 2 + 1 k i λ (n) i - λ(n-1) i 2 dΓ i∈{1,2} Γ k i w (n) i 2 + ŵ(n-1) i 2 + 1 k i λ (n) i 2 + λ(n-1) i 2 dΓ .
While this error indicator is commonly employed, it should be used with caution. First, it uses the search directions, which values therefore can potentially interfere in convergence. Secondly, in some problems, such as a frictional contact problem, tangential and normal parts of displacements and forces might be of a different order. Therefore, the standard error indicator is no longer suitable for such problems. Other convergence indicators should be considered, see, e.g., [START_REF] Ribeaucourt | A new fatigue frictional contact crack propagation model with the coupled X-FEM/LaTIn method[END_REF]. Finally, it provides only global information about the solutions, potentially missing localized phenomena.

Remark 1.4.1: Local stage resolution

The local stage equations are solved at the points belonging to the interface. Whilst FE nodes always lie on the geometry, it is not the case for control points in IGA. Thus, especially for non-linear interface behaviour and higher-order splines, the problem of solving the local stage arises. To overcome this issue, in this work (and so it is in some others, e.g., [START_REF] Claus | A stable and optimally convergent LaTIn-CutFEM algorithm for multiple unilateral contact problems[END_REF], [START_REF] Ribeaucourt | A new fatigue frictional contact crack propagation model with the coupled X-FEM/LaTIn method[END_REF]) local equations are satisfied at the integration points (see Chapter 3 for more details).

Implementation: Static case

Let us present the weak form of the linear stage equations to make it easier to understand the essence of the method. Regarding the linear stage, the problems to be solved are actually subdomain-wise elastic problems subjected to generalized Robin boundary conditions. Replacing λ i in Eq. (1.32) by its expression from Eq. (1.34b), and returning to u i by applying the trace over Γ, the problems reads: ∀i ∈ {1, 2}, find u i ∈ U i such that,

a i (u (n) i , v i )+ Γ k - i u (n) i •v i dΓ = l i (v i )+ Γ ( λ(n-1) i +k - i ŵ(n-1) i )•v i dΓ, ∀v i ∈ V i , (1.35)
where the functional spaces U i and V i contain the displacement solution and test functions, respectively; a i and l i being the standard bilinear and linear forms, respectively (see again Section 1.2.3 for their definition).

The equation of the linear stage (Eq. (1.35)) is written independently for each subdomain. Subdomains receive information through mixed boundary conditions applied on interface Γ. Also, it can be noticed, that during iterations, only the right hand side of the Eq. (1.35) changes. That means, that the matrix of the linear system can be computed only once, before the iterative procedure. However, this requires that the search direction k -is constant through the whole computation. There is no such a requirement for k + .

Considering the local stage, we aim at finding unknown interface hat values ŵ1 , ŵ2 , λ1 , λ2 and equations to be solved are:

λ1 + λ2 = 0; (Mechanical equilibrium) g( ŵ1 , ŵ2 , λ1 , λ2 ) = 0; (Interface behaviour) λ1 -λ 1 -k + 1 ( ŵ1 -w 1 ) = 0; (Search directions) λ2 -λ 2 -k - 2 ( ŵ2 -w 2 ) = 0.
(1.36)

These equations are defined at the integration points on interface Γ, meaning that at each Gauss point one should solve four equations to obtain four unknown hat values.

Remark 1.4.2

Here, we presented the LaTIn method in the static framework. However, nonlinear interface behaviour (e.g., frictional contact) demands to take into account the loading history. To do so, the quasi-static formulation exists, which is presented in Section 3.3.2.

On the parallelization of the LaTIn method

With the LaTIn method, we somehow naturally end up with a non-linear parallel domain decomposition algorithm. At each iteration, two stages (the local and the global ones) are performed consecutively. However, each of this stages can be efficiently parallelized:

Local stage: as the name suggests and as it can be seen from Eq. (1.36), all equations of the local stage are solved locally, i.e. at the interface points. This implies that the local stage can be computed in parallel.

Global stage: each subdomain communicate with the rest of the domain only by its interface. Indeed, as it can be seen in Eq. (1.35), subdomains receive forces and displacement in the form of mixed boundary conditions on the interface. Consequently, subdomains can be independently processed in parallel.

Conclusion

In this chapter, an overview of existing advanced simulation methods was provided.

In the first two sections, we addressed the question of discretization of a computa-44 CHAPTER 1. STATE OF THE ART: ON THE MULTI-SCALE MODELLING AND COMPUTATION OF COMPLEX STRUCTURES WITH IGA tional domain. IGA, utilizing smooth basis functions both for geometry modelling and solution approximation, reveals higher per-degree-of-freedom accuracy and robustness compared to the standard FEA for mechanical simulations (provided the solution is smooth), which makes this approach often seen as a HPC tool. Nevertheless, it is rather challenging to model geometrically complex structures within the framework of IGA. It requires performing multi-patch coupling, which leads to reduced regularity at the coupling interfaces. Thus, the compelling advantage of increased regularity in IGA is not fully harnessed in this case. For this reason, we appeal to immersed boundary methods. According to this approach, the mesh is no longer required to be geometrically conformal to the physical boundaries of the domain considered. The tedious meshing process is thus alleviated. The combination of IGA with the immersed boundary approach appears to be a promising blend, capitalizing on the superior behaviour of IGA and the geometric flexibility of immersed boundary methods. Using immersed boundary approaches results in elements cut by physical boundaries. Hence, the necessity for non-conformal coupling becomes apparent. The standard Lagrange multipliers (Mortar) coupling strategy was discussed, followed by the presentation of Nitsche's approach, with an emphasis on the advantages of the latter in the context of linear elasticity.

In the last two sections, different approaches to solving the mechanical problem were discussed. When dealing with localized non-linearities, the non-invasive local/global coupling approach becomes particularly interesting. Indeed, it allows for incorporating complex local behaviours within a global model without modifying global operators. Finally, the LaTIn method was introduced as a useful tool for modelling multiple an various non-linear interface behaviours. The principle of the LaTIn approach for domain decomposition is to separate the non-linear interface equations from those related to the subdomains, with the latter being linear and subdomain-wise independent. Therefore, the advantage is that the subdomain equations can be solved efficiently in parallel, and the non-linear interface equations can be written locally in an easy manner, regardless of the treated non-linearity.

With these methods in hand, we propose to use IGA at least at the global scale. The use of IGA at a local scale, if the local solution is expected non-smooth, e.g., with sharp discontinuities or singularities, may not be reasonable. The next chapters address in more details these two problems. Firstly, we will consider a non-invasive global IGA/local FEA coupling to deal with any non-linear local models in a global one. Following that, we will introduce a fully IGA method tailored for scenarios involving non-linear interfaces within a globally linear structure. To achieve this, we will leverage the LaTIn method to model a range of non-linear behaviours effectively.

A fully non-invasive hybrid IGA/FEA scheme for the global/local simulation of structures

In this chapter, a new hybrid modelling approach is proposed. It leverages the advantages of both IGA and FEA, offering efficiency in global analysis and versatility in local behaviour modelling while preserving the non-invasive character of the coupling process. The non-invasiveness is offered not only in terms of global/local coupling, but also in the way coupling operators are constructed. The isogeometric operators are built from their FE counterparts in a straightforward way using FEA software. Finally, the approach is illustrated throughout a series of numerical examples.

Contents

Introduction

While dealing with modelling mechanical structures, it occurs that some particular phenomena characterized by a complex behaviour is present locally. For instance, cracks, local contact, delamination or local heterogeneities can be present in a global structure. These localized complexities introduce challenges that necessitate special modelling and simulation techniques to accurately capture and analyse the behaviour of these intricate regions. In the context of IGA, the main implication of the presence of such phemonema is that in these local zones of interest, the solution sought is not highly regular any more. As stated before, the one of the main advantages of using the higher-order and smooth basis functions is the increased per-Degree-of-Freedom (DoF) accuracy. However, this is applicable only for regular solutions. Thus, IGA does not appear suitable when the solution is not regular any more. As a result, numerous sophisticated methods have been developed over the years to make possible local simulations within IGA. For example, regarding for instance fracture and/or delamination, one may refer to the IG version of XFEA, namely XIGA [START_REF] Luycker | X-FEM in isogeometric analysis for linear fracture mechanics[END_REF][START_REF] Fathi | Geometrically nonlinear extended isogeometric analysis for cohesive fracture with applications to delamination in composites[END_REF][START_REF] Yuan | Multi-patch local mesh refinement XIGA based on LR NURBS and Nitsche's method for crack growth in complex cracked plates[END_REF], or to IG cohesive elements [START_REF] Dimitri | NURBS-and T-spline-based isogeometric cohesive zone modeling of interface debonding[END_REF], Verhoosel et al., 2011), or even to phase-field approaches [START_REF] Borden | A higher-order phase-field model for brittle fracture: Formulation and analysis within the isogeometric analysis framework[END_REF][START_REF] Paul | An adaptive space-time phase field formulation for dynamic fracture of brittle shells based on lr nurbs[END_REF][START_REF] Proserpio | A framework for efficient isogeometric computations of phase-field brittle fracture in multipatch shell structures[END_REF]. However, all these methods seem to have a very high level of complexity and therefore may potentially demand considerable effort for implementation within a conventional IGA framework. Contrary to IGA, FEA appears adapted to simulate local, strongly non-linear or even singular behaviours due to its reduced regularity and its meshing flexibility. As stated previously (see Section 1.1.4), FEA benefits from more that 50 years of developments and practices so numerous enhanced FE implementations, both efficient and robust, exist to simulate various local behaviours.

In this context, we introduce a hybrid global-IGA/local-FEA modelling approach, where we couple a global IsoGeometric (IG) model with a local Finite Element (FE) model. This coupling is achieved in a non-invasive manner, as discussed in Section 1.3. Thus, this approach benefits from both analysis technologies: efficiency of IGA for geometric description and for capturing global, regular response; and, ability of FEA to compute local, strongly non-linear or even singular behaviours. Moreover, the noninvasiveness of the approach appears attractive when coupling a global IG model with a local FE one. Firstly, there is no need to incorporate a local region within the initial IG model, thus, it naturally avoids costly spline re-parametrization procedures. Secondly, the global IG stiffness operator can be assembled and factorized only once and the IG system to be solved remains well-conditioned regardless of the shape of the local region. Finally, it offers the opportunity to simply couple an IG code with any existing robust FE code suitable for the modelling of complex local behaviours, using the FEA-to-IGA bridge discussed in Section 1.1.3.

This introduction is followed by presenting the proposed strategy in Section 2.2, starting by the reference mechanical problem in Section 2.2.1, then the construction of FE mesh with a conformal interface is presented Section 2.2.2, and, finally, the fully noninvasive approach is shown in Section 2.2.3. Then, in Section 2.3, the performance of the proposed implementation is demonstrated through a series of benchmarks involving complex local behaviours. To conclude, the chapter wraps up with a summary in Section 2.4.

Non-invasive global-IGA/local-FEA

We introduce here our hybrid global-IGA/local-FEA scheme that seeks to combine the interests of both technologies for global/local simulations. Let us underline at this stage that our approach is generic in terms of programming environments: the users may have in hand an IG code (performing standard elasticity) and wish to couple it with a specific FE software to model complex local phenomena, or the users only have at their disposal FE packages. The aim is to arrive at an automatic coupling between global-IGA and local-FEA.

To couple global IGA and local FEA models, the non-invasive global/local algorithm (discussed in Section 1.3) is applied. Moreover, specific construction of the models along with the dedicated implementation is used to achieve a fully non-invasive hybrid IGA/FEA procedure. Our approach relies on the complete algebraic bridge that directly goes from Lagrange nodal polynomials to B-Spline and NURBS functions, discussed in Section 1.1.3. Thus, IG operators (stiffness matrix and load vector) are constructed in a straightforward manner from their FE counterparts computed using a standard FE code (see Eqs. (1.14) for the link between IG and FE operators). This allows to perform the implementation even without having a global IG code in hand, as will be shown in Section 2.2 (see Remark 2.2.3).

The global IGA/local FEA problem

In our approach to use non-invasive global/local coupling, we focus on the mechanical coupling problem discussed in Section 1.

Construction of the FE model to reach a conformal global/local interface

Incorporating a specific local region in an IG patch without care may result in the overlap of some global knot-span elements due to the rigid tensor product structure of (standard) multivariate spline bases (see Figures. 2.1(a) and 1.9 again). More precisely, the difficulty relies on (i) the evaluation of integrals over pieces of knot-span elements (to get the interface reaction force λ IG 12c , see Eq. (1.30), here, subscript c emphasizing that it is associated with the coarse discretization), and (ii) the formulation of a coupling method adapted to an immersed interface. As a remedy, the idea here is to call upon efficient (classic) FE meshing procedures to reach a conformal interface (similar pragmatic approach as in [START_REF] Guinard | Multiscale analysis of complex aeronautical structures using robust non-intrusive coupling[END_REF] in the context of global/local FEA). In order to do so in a simple, automatic and consistent way, we make use of the FEAto-IGA bridge of Section 1.1.3, which will also enable to arrive at a fully non-invasive strategy in the sense that not only the global/local coupling is non-invasive but also the construction of the coupling operators (see Eq. (1.29)) from only FE resources.

Figure 2.2: Illustration of the proposed procedure to build in a simple, automatic and consistent way a conforming global-IGA/local-FEA discretization. The strategy strongly relies on the FEA-to-IGA bridge, which also offers the opportunity to simplify the implementation of the mechanical solver by involving only FE operators.

The proposed procedure for the construction of the conforming global-IGA/local-FEA modeling is illustrated in Figure 2.2. This figure also presents the notations used in this chapter. This notation follow the notation of Chapter 1 (Sections 1.3 and 1.1 more) with some additional details: superscripts IG and F E specify the type of discretization associated while subscripts c and f precises if the mesh is coarse or refined. Starting with a global coarse IG model of the whole structure (a), a specific local FE model (e), meant to replace the global IG model in an area of interest, is built. The area of interest is obviously chosen as a set of initial global knot-span elements so that we are able to reach a conformal coupling interface. More precisely, from the initial coarse global IG model (a), we apply standard spline refinement procedures (associated operator D IG 1cf , see Remark 1.1.2) to obtain the refined global IG model (b). Then we can obtain the corresponding refined global FE model (c) using the FEAto-IGA bridge presented in Section 1.1.3, operator D F E 1f performing this passage from refined global IG model to refined global FE one (see operator D F E in Eq. (1.14)). It is therefore possible to extract the FE description of the interface (d) by calling upon a trace operator (T 1f ) that selects only the nodes (or the DoF) concerned with the interface. Finally, by applying existing optimized FE meshing procedures (e.g., GMSH [START_REF] Geuzaine | Gmsh: A 3-D finite element mesh generator with builtin pre-and post-processing facilities[END_REF], or Salome-Meca [START_REF] Ribes | Salome platform component model for numerical simulation[END_REF] which will be used for the numerical illustrations), we can build a local FE mesh that is conformal to the interface. In order to do so, it may be noticed at this stage that the constructed FE mesh must have the same polynomial degree as the interface.

As is standard practice in conformal global/local FEA, we choose to take the trace along the interface of the functions of the local model (L T ) to discretize the Lagrange multiplier field. Thus, Mortar coupling operator from system (1.29) becomes:

C F E 2 = Γ L T L T 2 dΓ,
where L 2 refers to the shape functions associated with the local FE model. Thus, C F E 2 is the mass matrix associated with the local FE model of the FE interface. We never encountered instabilities in our numerical experiments with such a choice.

Fully non-invasive implementation of global-IGA/local-FEA

In general, the computation of the Mortar coupling operators may not appear trivial.

It is necessary to build an integration technique on the interface and to have access to the values of the shape functions on both sides at each integration point, which is not a classical output of industrial codes and breaks the concept of non-invasiveness. The meshing strategy presented in Figure 2.2 offers a solution to circumvent these challenges. By utilizing the FEA-to-IGA bridge, the hybrid IGA/FEA coupling operators become explicit and can be computed more readily. Indeed, operator C 1 in Eq. 1.28, which is denoted by C IG 1c in the present context (since it is associated with the coarse IG model) can be computed as follows:

C IG 1c = Γ L T R T 1c dΓ = Γ L T L T 1f dΓ D IGF E 1cf T = C F E 1f D IGF E 1cf T ,
where

D IGF E 1cf = D IG 1cf D F E 1f
is an operator composed of the FEA-to-IGA operator D IGF E and the spline refinement operator D IG cf . Operator C F E 1f is the mass matrix associated with the refined global FE mesh of the FE interface. With these computations and following the notation of this chapter, Eqs. (1.28)-(1.29) become:

1. Resolution of a Neumann problem (with respect to Γ) over the initial (coarse) IG model:

K IG 1c u IG 1c (n) = f IG 1c -D IGF E 1cf C F E 1f T λ (n-1) + λ IG 12c (n-
  K F E 2 -C F E 2 T -C F E 2 0     u F E 2 (n) λ (n)   =   f F E 2 -C F E 1f D IGF E 1cf T u IG 1 (n)   . (2.2)
Thus, there is no longer operators that merge basis functions from IGA and FEA. Nevertheless, the computation of the FE interface mass matrices C F E 1f and C F E 2 may still not be straightforward in practice using commercial FE codes. However, as often performed in a transparent manner when coupling domains in FEA, these operators are not truly required here; the trace operators T 1f and T 2 , which are more readily available in FE codes, are actually sufficient. To highlight this, we first introduce the FE mass matrix of the interface:

C F E T = Γ L T L T T dΓ,
which is invertible (it symmetric definite positive by construction). Then, we take λ = C F E T λ which has the dimension of a load vector, and we multiply the second part of Eq. (2.2) by C F E T -1 . With these manipulations and making use of equalities:

T 1f = C F E T -1 C F E 1f ; T 2 = C F E T -1 C F E 2 ,
Eqs. (2.1) and (2.2) are then modified to use these trace operators, simplifying the coupling process (modifications highlighted in grey colour):

K IG 1c u IG 1c (n) = f IG 1c -D IGF E 1cf T T 1f λ(n-1) + λ IG 12c (n-1) ; (2.3)   K F E 2 -T T 2 -T 2 0     u F E 2 (n) λ(n)   =   f F E 2 -T 1 D IGF E 1cf T u IG 1 (n)   .
(2.4)

The trace operators being most of time available in FE codes, the coupling can be performed using only FE industrial packages, making the methodology fully non-invasive.

The passage is further illustrated by

Remark 2.2.1

In case if the users have only FE codes at their disposal, it is possible to reformulate the iterative process in terms of only FE operators. To do so, firstly, we note that the initial coarse IG solution u IG c can be obtained from the refined 2.3. NUMERICAL RESULTS

51

FE solutions with the use of operator D IGF E cf ad follows:

K IG c u IG c = f IG c ⇔ D IGF E cf K F E f D IGF E cf T u IG c = D IGF E cf f IG f .
This leads to the following link between operators related to the initial coarse IG mesh and the refined FE mesh:

K IG 1c = D IGF E 1cf K F E 1f D IGF E 1cf T and f IG 1c = D IGF E 1cf f F E 1f ;
(2.5)

Reaction forces at interface Γ and displacement vectors can also be connected:

λ IG 12c = D IGF E 1cf λ F E 12f with λ F E 12f = K F E 12f u F E 1f -f F E 12f , u F E 1f = D IGF E 1cf T u IG 1c .
This leads to the following iterative process:

D IGF E 1cf K F E 1f D IGF E 1cf T u IG 1c (n) = D IGF E 1cf f F E 1f -D IGF E 1cf T T 1f λ(n-1) + D IGF E 1cf λ F E 12f (n-1) ; (2.6)   K F E 2 -T T 2 -T 2 0     u F E 2 (n) λ(n)   =   f F E 2 -T 1 D IGF E 1cf T u IG 1 (n)   , (2.7)
where only FE operators are required (see Figure 2.2 for the notations).

Remark 2.2.2

With the proposed method, it is possible to take a higher-order IG model (i.e., p > 1) with a classic low-order FE mesh (i.e., p = 1). In this case, we take as many nodes at the interface for the refined global FE model (see Figure 2.2(c)) as for the local FE model (see Figure 2.2(e)). This means that, for instance, for quadratic IGA versus linear FEA, there are two FE elements in front of one IG element along the interface. Two approximations are thus performed implicitly: (i) the (possibly curved) interface of the IG model is facetted between each of the interface nodes for the local model, (ii)

C F E T -1 C F E
1f is not strictly equal to T 1f any more (the mass matrix of a quadratic FE interface is approximated by the mass matrix of a linear FE interface with twice more elements). However, our results seem to indicate that the error related to such approximations is very low compared to that associated with the discretization (see section 2.3.1 and Figure 2.7).

Numerical results

To assess the performance of the developed non-invasive hybrid global-IGA/local-FEA algorithm, we now present a series of numerical experiments that cover 2D and 3D simulations with different local behaviours, such as cracks, contact and delamination. All the implementations have been carried out using the open-source FE industrial software package code-aster (ÉDF, 1989(ÉDF, -2017) ) developed by the EDF R&D company. No IG codes have been used; we consider the case where we have at our disposal only the FE code code-aster. We thus more precisely implemented algorithm (2.6)-(2.7) and limited ourselves to quadratic spline functions since code-aster does not go 52 CHAPTER 2. A FULLY NON-INVASIVE HYBRID IGA/FEA SCHEME FOR THE GLOBAL/LOCAL SIMULATION OF STRUCTURES beyond second-order Lagrange finite elements. Yet, we underline that the proposed implementation schemes (2.1)-(2.2) and (2.3)-(2.4) could be applied to higher-order splines if one has an IG code in hand to compute the IG stiffness and load vector operators. The automatic procedure described in Figure 2.2 for the construction of the conforming global/local discretization was performed using the mesh generator Salome-Meca [START_REF] Ribes | Salome platform component model for numerical simulation[END_REF]

Linear elastic 2D curved beam

The first example consists of a 2D linear elastic curved beam subjected to end shear adapted from [START_REF] Zienkiewicz | The finite element method: its basis and fundamentals[END_REF]. Such an example has been widely used in IGA to assess the performance a method. The global geometry was perfectly generated using a single NURBS patch composed of only one quadratic element. and 2.4c in terms of displacement and of Von Mises stress, respectively. We note that it is the converged solution in Ω 11 ∪ Ω 2 that is mapped (the fictitious prolongation of the global solution over Ω 12 is not represented). For all examples, we will perform this way for the visualization. On this simple case, the iterative non-invasive algorithm converges very quickly without acceleration techniques: 3 iterations were needed with a stopping criterion based on the equilibrium of the interface global and local reaction forces (tolerance of 10 -8 here). The solution appears smooth and in a good agreement with [START_REF] Zienkiewicz | The finite element method: its basis and fundamentals[END_REF].

To go further, the convergence of the method with the refinement of the mesh was studied. In order to do so, the computation depicted in Figure 2.4a was repeated for several global/local discretizations. Starting with a global quadratic NURBS mesh of 6 (circumferential direction) × 4 (radial direction) elements, the refinement was increased to reach 12 × 8, 24 × 16, 48 × 32 and 96 × 64 elements. We kept the same region Ω 12 and the replacement of 1 NURBS element by 4 T6 triangles for the local model for each discretization. We proceeded in the same way as in [START_REF] Zienkiewicz | The finite element method: its basis and fundamentals[END_REF]; that is, the convergence behaviour of the strain energy is considered through to computation of the relative energy error: tion when considering the global model everywhere and "Standard FEA" corresponds to the solution when all the NURBS elements are replaced by 4 T6 triangles. The results show that the same rate of convergence was achieved with the proposed hybrid IGA/FEA scheme as with the reference solutions, which accounts for the accuracy of our method. More specifically, it can be noticed that for a given mesh refinement, the errors are about the same for the three solutions. Only the number of DoF changes: it decreases when IGA is used. This illustrates the increased per-DoF accuracy of IGA and is totally consistent with our interpretation on IGA as a projection of FEA onto a reduced, regular basis. The solution here being smooth, it is well captured with IGA as with FEA, but IGA comes with less DoF due to its higher regularity. Our hybrid IGA/FEA solution obviously appears in between the two reference curves since the DoF of the global and local models are summed. Let us underline here that our way of counting the DoF has not a concrete meaning from a computational cost point of view since the IG and FE problems are solved separately in our non-invasive strategy. One could have chosen to take the maximum number of DoF between the global and local models, i.e. the same number of DoF as for the standard IGA solution which would have led to the superposition of the "Standard IGA" and "Non-invasive global-IGA/local-FEA" curves.

Err h = |E ref -E h | E ref where E ref
For completeness on this test case, we finally carried out the same numerical experiments but with a local FE model composed of linear triangles (i.e., T3 triangles). For this purpose, the strategy described in Remark 2.2.3 was applied: two FE elements were put in front of one NURBS element along the interface. Figure 2.6a shows a zoom on the local region when considering a global model made of 24 × 16 elements. This time, 16 T3 triangles replaced 1 NURBS element in the local region. The results for this discretization are given in Figures 2.6b and 2.6c in terms of displacement and of Von Mises stress, respectively. Of course, some discontinuities for the stress can be observed in the local region since this field is now piecewise constant in this area, but the solution still appears ro be in a good agreement with [START_REF] Zienkiewicz | The finite element method: its basis and fundamentals[END_REF]. The convergence curve in terms of relative energy error is then plotted in Figure 2.7 along with the equivalent single-model solutions; that are, the solution when considering the global quadratic IG model everywhere and the solution when all the NURBS elements are replaced by 16 T3 triangles. As expected, the convergence rate of our hybrid quadratic-IGA/linear-FEA scheme is now driven by the linear-FEA part of the solution (same convergence rate for the coupled solution as for the linear-FEA case). However, the coupled solution appears much more accurate than the standard linear-FEA one (drastic reduction of the constant factor) which is due to the higher order and higher regularity of the global model. This numerical experiment validates the proposed procedure to couple higher-order IGA with standard linear FEA (see again Remark 2.2.3) and further confirms the interest of making use of IGA for the global response to reach an increased per-DoF accuracy.

2D curved beam with holes, cracks and contact

With the second example, we illustrate the potential of our methodology to include geometrical details along with non-linear behaviors within a global NURBS model. More precisely, the global linear elastic 2D curved beam model of the previous test case was recycled. This time, it was clamped on its lower part and subject to a vertical distributed load on its left edge, see Figure 2.8(left). The global mesh was composed of 24 × 16 quadratic NURBS elements. This model was locally enriched by a specific FE mesh made of quadratic triangles (i.e., T6 triangles) incorporating several holes and cracks. In addition, frictional contact was modeled between the lips of the cracks. The cracks were initially open (initial gap between the crack lips of about 0.004 mm). A zoom on the local model is performed in Figure 2.8(right). Let us underline that solving such a problem within the sole IG framework would be delicate due to the geometric complexity of the local region (perhaps, some advanced IG immersed technologies [START_REF] Ruess | Weak coupling for isogeometric analysis of non-matching and trimmed multi-patch geometries[END_REF][START_REF] Wang | Consistent immersed volumetric Nitsche methods for composite analysis[END_REF], Wei et al., 2021) or extended IGA [START_REF] Luycker | X-FEM in isogeometric analysis for linear fracture mechanics[END_REF][START_REF] Fathi | Geometrically nonlinear extended isogeometric analysis for cohesive fracture with applications to delamination in composites[END_REF][START_REF] Yuan | Multi-patch local mesh refinement XIGA based on LR NURBS and Nitsche's method for crack growth in complex cracked plates[END_REF] or phasefield procedures [START_REF] Borden | A higher-order phase-field model for brittle fracture: Formulation and analysis within the isogeometric analysis framework[END_REF][START_REF] Paul | An adaptive space-time phase field formulation for dynamic fracture of brittle shells based on lr nurbs[END_REF][START_REF] Proserpio | A framework for efficient isogeometric computations of phase-field brittle fracture in multipatch shell structures[END_REF]) should be required). On the contrary, realizing a locally boundary-fitted mesh is straightforward with standard FEA. Besides, the displacement has to be discontinuous on either side of the cracks which is natural with a FE mesh that fits the cracks.

The Von Mises stress map obtained once the non-invasive algorithm has converged is given in the whole top crack closes while only subparts of the other cracks are in contact. In addition, slight sliding can be observed on the large bottom crack. For completeness, the global-IGA/local-FEA displacement obtained through our non-invasive strategy was compared to the displacement field obtained with a full FE discretization of the whole problem (the FE mesh that allows to recover the IG global solution by projection (see Eq. (2.5)) was used in Ω 11 ). Since the solutions are very close, the relative discrepancy between the hybrid solution and the FE one is plotted in Figure 2.9c. Less than 2% of local mismatch can be observed while the two solutions come from different approximation subspaces in Ω 11 . This result confirms that IGA is sufficient to accurately capture the global response, even when the local one exhibits some discontinuities.

Finally, the convergence of the non-invasive algorithm was investigated in Figure 2.10. For this example, the Aitken's dynamic relaxation acceleration [START_REF] Duval | Non-intrusive coupling: recent advances and scalable nonlinear domain decomposition[END_REF], Gosselet et al., 2018) seems to be necessary, which was expected since the stiffness gap between domains Ω 12 and Ω 2 was significant. Taking advantage of the Aitken's update, we were able to make the number of iterations relatively low: a residual of 10 -4 was obtained in 23 iterations.

2D plate with multiple inclusions and delamination

To illustrate our approach in the context of modelling of composite materials, the third example constitutes a 2D illustration of a matrix and several inclusions inside. At the matrix/inclusion interface delamination was modelled using cohesive zone law. Moreover, this test case will allow to highlight another attractive property of our noninvasive algorithm: it results in an efficient non-linear domain decomposition solver when several local models are considered.

More precisely, the structure was composed of a square shape matrix containing 1, 4 or 16 (uniformly distributed) circular inclusions, as depicted in Figure 2.11. Symmetry boundary conditions were applied on the left and bottom edge of the structure, while a constant horizontal displacement was prescribed over the right edge, which results in a structure globally subjected to traction. The global model, which obviously describe here the square shape matrix (without inclusions), was discretized with 16 × 16 quadratic IG elements. This IG model was then enriched by as many local FE models as there are inclusions. The associated local FE meshes were the same for all the inclusions and they can be seen on Figure 2.12. They all included one inclusion along with 1cf . In addition, we followed Remark 2.2.3 to adopt a global quadratic-IGA/local linear-FE modelling. We first refined the IG global model 4 times and applied the FEA-to-IGA bridge to obtain a refined global FE model made of 64×64 quadratic FE elements. We then extracted the FE nodes on the global/local interface and built the FE mesh of Figure 2.12 by putting T 3 elements between each interface nodes. We eventually considered the same material for the inclusions and the matrix: a linear elastic isotropic material with the Young modulus E = 10000 MPa and Poisson ratio ν = 0.3. The cohesive elements at the matrix/inclusion interfaces followed a standard bilinear law, see Figure 2.13. This law is described by 3 parameters: σ c the critical stress, G c the density of the critical energy of the material and ρ a penalisation coefficient to control the stiffness of the undamaged domain. In our case, we chose: σ c = 10 MPa, G c = 0.005 mm.MPa and ρ = 0.1. The displacement jumps δ 0 and δ c , corresponding to the initiation and completion of damage are derived from these parameters as δ 0 = ρ • 2 Gc σc and δ c = 2 Gc σc . At the beginning, the interface is characterized by the initial stiffness k 0 . As the displacement jump increases, the critical stress σ c is reached and the interface starts to be damaged, the surface force decreases and in the case of unloading, the elastic return follows a damaged stiffness k. structure. The problem was solved for each loading step using our non-invasive algorithm (2.6)-(2.7) accelerated with the Aitken dynamic relaxation. In terms of results, we first show in Figure 2.14 that the convergence of the non-invasive global/local algorithm does not depend on the number of inclusions. Underlining that the different non-linear local problems can be naturally solved in parallel here, this accounts for the scalability of the algorithm. In other terms, our approach can be used as an efficient non-linear domain decomposition solver, as proposed in [START_REF] Duval | Non-intrusive coupling: recent advances and scalable nonlinear domain decomposition[END_REF] in standard FEA. The global model actually plays the role of the coarse problem of domain decomposition approaches, which allows to transmit directly the information all over the sub-domains. Furthermore, it may be noticed that the convergence of the non-invasive algorithm is similar for each loading step: between 10 and 15 iterations are sufficient to obtain a residual of 10 -4 . Then, Figure 2.15 shows the stress distribution for the 16-inclusion problem at the last loading step. Stress concentrations can be observed around the inclusions which means that parts of the interfaces were severely damaged. This is further confirmed in Figure 2.16 where a zoom over one inclusion is performed. The stress distribution can be better appreciated in Figure 2.16a and in Figure 2.16b we plot the status of the damage in the cohesive elements. We see that almost all the interface were damaged and moreover the right part was almost completely broken. Finally, in Figure 2.17, we plot the reaction forces according to the prescribed displacement. The impact of the damage on the cohesive zones is noticed as the reaction force (blue curve) decreases compared to the curve with undamaged cohesive zones (dashed red curve). Let us note here that the load steps were quite large and this avoided numerical difficulties to capture some snap-back phenomena. Indeed with potential brutal loss of stiffness due to quick rupture of the cohesive zone, the precedent curve could involve snapbacks. With a finer load discretization, the Newton solver could fail and continuation techniques should be used.

3D mechanical assembly example with preload and frictional contact

In this example, we apply our approach in the context of assembly modelling to show its versatility. In this section, a 3D frictional contact problem of a bolted assembly . The purpose of this example is to show that the developed non-invasive hybrid procedure can be easily applied to solve general 3D contact problems. The assembly was composed of two quarter cylinders and one bolt that acts on two perforated plates which extend each of the pieces of cylinders. For the modeling, the global model was only made of a continuous half cylinder with quadratic NURBS elements (see blue region in Figure 2.18). This half cylinder was cut and bolted through a complex local FE model (see gray region in Figure 2.18). More precisely, the local model consisted in five solids: two perforated plates linked to the global model across the global/local interface, and a screw and two nuts that enforced the contact between the two plates thanks to the application of a preload in a pre-processing step. The dark blue region in Figure 2.18 thus concerns domain Ω 12 , the NURBS solution is replaced by the local FE solution in this part. Overall, the local model included 5 interfaces (see again Figure 2.18): a frictional contact interface between the two plates (see red line), two preload interfaces between the nuts and the screw that were used to impose a relative displacement in the bolt (see green lines), and two perfect interfaces between the plates and the nuts (see blue lines). In the end, these five interfaces produced tension in the screw and compression in the plates. As for the discretization, quadratic tetrahedrons and pyramids were used for the local model and the strategy depicted in Figure 2.2 was applied in 3D to recover a conforming global/local interface (see Figure 2.19). Eventually, displacement boundary conditions were prescribed at the bottom of the global model (half cylinder) to extend the bolted assembly.

The convergence of the non-invasive algorithm is provided in Figure 2.20. An Aitken's acceleration was necessary to obtain a residual of 10 -4 in a few tens of iterations as in the previous non-linear test cases. Then, the deformed shape along with the Von Mises stress is depicted in Figure 2.21. Once again, a smooth transition of the stress is observed at the global/local interface which confirms the performance of our non-invasive coupling scheme in 3D. Moreover, due to the applied Dirichlet boundary conditions, a slight detachment can be observed between the two quarter cylinders near the interior radius whereas the bolt maintains the contact around the screw. Obviously, stress concentrations can also be reported close to the nuts. Finally, we show in Figure 2.22 the parts of the assembly that are in compression (i.e. where σ xx is negative). As expected, the screw is in tension due to the preload and leads to a large compression Half of the structure is shown (a cut is performed in the symmetry plan of the structure). zone under the nuts, which helps keeping the plates in contact. This good behavior of this last case demonstrates the potential of our method to treat more representative applications.

Remark 2.3.1

Let us notice that several bolts could be easily considered using exactly the same strategy but with several local models, in the previous test case. The superiority of the non-invasive approach would be even clearer from a modelling and a solution point of view. Indeed, solving a multi-contact problem in a monolithic way is not an easy task as convergence properties of the non-linear solvers are worsened by the increase of the number of contact surfaces. Here, the different contact problems associated with each local models could be naturally solved in parallel, therefore we would end up with a non-invasive, non-linear Domain Decomposition strategy (see again [START_REF] Duval | Non-intrusive coupling: recent advances and scalable nonlinear domain decomposition[END_REF] for more information on this topic).

Conclusion

In this chapter, we tackled the problem of modelling presence of local non-linearities in a global structure. To end up with an efficient approach, we combined different technologies: global (more smooth) solutions are captured by the means of IGA, while FEA is employed to simulate local complex behaviours. This allows to benefit from both modelling techniques.

Our approach can be viewed from different perspectives. On one hand, it enhances traditional full global/local FEA simulations by establishing a direct link with CAD for global geometry representation and by significantly reducing the DoF while maintaining accuracy in the global region. On the other hand, it enriches global/local IGA by enabling the modelling of complex local behaviours within IG patches through the utilization of FEA routines tailored for this purpose.

The global/local coupling is performed in the non-invasive way. The method relies on an iterative Dirichlet-to-Neumann process where the initial global IG model to be enriched is never modified. This constitutes the first key point in IGA to reach a noninvasive strategy since it avoids costly spline re-parametrization procedures that may have been necessary otherwise to incorporate a free local region. In addition, it has to be noted that since only the initial global IG operator is involved, the IG system to be solved remains well conditioned regardless of the shape of the local region.

Moreover, we resorted to the existing FEA-to-IGA bridge, based on Bézier or Lagrange extraction, to transform the interface within the initial global IG model into a FE interface on which the local FE mesh can be constructed by calling upon efficient (classic) meshing procedures. The situation of a conforming coupling interface was thus reached and we were able to express the hybrid IGA/FEA coupling by means of only the FEA-to-IGA operator (easily built from any spline libraries) and standard FE trace operators that are available in (possibly industrial) FE codes. It results that any robust FE code suitable for the modelling of complex local behaviours can be used in a plug-and-play manner. Finally, depending on how deeply the FEA-to-IGA operator is used, our implementation is generic in terms of programming environments: the users may have in hand an IG code (performing standard elasticity) and wish to couple it with a specific FE software to model complex local phenomena, or the users only have at their disposal FE packages. To validate the strategy, a series of numerical experiments was conducted. The implementation was performed using the open-source FE industrial software package code-aster ÉDF (1989ÉDF ( -2017) ) developed by the EDF R&D company. Firstly, an academic problem was examined, and the results were promising in terms of the obtained solution quality and mesh convergence. Then, it was enhanced with presence of holes, cracks and frictional contact at the local scale, to show the potential of the methodology. Further, we shifted our focus to composite materials, and a matrix with several inclusion was considered. Each local model comprised one inclusion with adjacent matrix material and delamination was considered at the matrix/inclusion interface. It allowed to highlight that our strategy results in an efficient non-linear domain decomposition approach. Finally, we shifted our focus to 3D assembly modelling, where frictional contact was considered within the local model. Thus, even in this thesis we focus on composite microstructures, the proposed strategy can be successfully used in other contexts.

Immersed boundary-conformal IsoGeometric LaTIn method for multiple non-linear interfaces

The objective of this chapter is to present an efficient solver for the simulation of geometrically complex composite microstructures involving numerous inclusions connected with the matrix through various non-linear interfaces. The proposed approach is fully-IGA, it leverages the recently introduced immersed boundary-conformal method and is enhanced with the Large Time INcremental method. This results in an immersed hybrid mixed higher-order numerical scheme that is flexible to treat any non-linear interface behaviour. The accuracy and efficiency of the developed algorithm are demonstrated by solving a range of non-linear examples in 2D, including different numbers of inclusions in unilateral contact, frictional contact, and delamination with the matrix of the composite microstructure.

Contents

Introduction

In many engineering applications, mechanical structures are often composed of multiple elastic domains interconnected by potentially non-linear interfaces. The presence of such interfaces makes the simulation much more complicated. The accurate representation requires detailed mesh as well as implementation of complex behaviour laws. These reasons make the final simulation both time-and resources-consuming.

Modelling such structures can be considered as a particular case of modelling localized non-linearities, which was considered in Chapter 2. Here, non-linearities are presented only at the interfaces and not within the bulk. In the context of fibre-reinforced composites, it is a problem of such complex behaviour as contact or delamination on the interface between the matrix and fibres, while the matrix and the fibres follow linear elasticity.

To tackle the problem of simulation of linear structure with possibly non-linear interfaces, this chapter proposes a pragmatic approach. As a first step towards efficiency, we use IGA both for the geometry representation of interfaces and the discretization of the domain. Thus, conversely to the non-invasive global-IGA/local-FEA strategy presented in Chapter 2, here we use IGA on all levels. This allows to fully benefit from higher per-degree-of-freedom accuracy and robustness compared to the standard FEA. Nevertheless, as it was shown in Section 1.1.2 (see Figure 1.3) modelling multi-scale and geometrically complex objects, such as composite microstructures involving a great number of local inclusions, is challenging in IGA due to the rigid tensor-product structure of the basis functions. To tackle this issue, we resort to the family of immersed boundary methods (see Section 1.2), thus allowing the mesh not to follow the geometry of the model.

In this case, the main challenge consists in accurate coupling formulation and implementation between the matrix and the multiple inclusions. This appears far from trivial since the coupling interfaces are expected to be both non-conformal (the interface can cut arbitrary elements of the meshes) and non-linear (such as incorporating contact or delamination). As a remedy, the first ingredient of our approach is to consider the recently introduced Immersed Boundary-Conformal Method (IBCM) [START_REF] Wei | Immersed boundary-conformal isogeometric method for linear elliptic problems[END_REF] that has proved to be accurate and robust for perfect (linear) interfaces. Our approach consists in a pragmatic strategy that transforms the initial interface, through the construction of conformal layers from it, into three different interfaces: the initial one between the matrix and the inclusion that becomes conformal, and two non-conformal ones that now lie within the matrix and inclusion. Therefore, the strategy leverages the geometric flexibility of the immersed methods with the advantages of conformal discretizations, which seems to us all the more relevant in our context. Indeed, (i) the solution and especially the stress fields can be properly described around the interfaces with controlled conformal discretizations, (ii) complex non-linear laws can be enforced in standard manners through the conformal interfaces, and (iii) the non-conformal aspect can be treated within the framework of linear elasticity, which we perform by Nitsche's approach, presented in Section 1.2.3. Moreover, let us note that the difficulty of building conformal layers is generally mitigated when the method is applied to composite materials, as the local inclusions may have rather simple geometries (for instance, they merely consist of circular cylinders in case of fibre-reinforced composites).

In this context, we present in this chapter an immersed boundary-conformal IsoGeometric LaTIn method that aims in accurate and efficient computation of composite microstructures involving multiple inclusions. In our framework, the LaTIn method is therefore applied across conformal interfaces without any difficulty, which results in the local formulation of the non-linear interface equations at each interface Gauss point, and in the parallel solution of subdomain-wise Robin-type linear elastic problems. The stabilization of the formulation is moved within the bulk equations related to the matrix and the inclusions, where Nitsche couplings with minimal stabilization are performed [START_REF] Antolin | Overlapping multipatch isogeometric method with minimal stabilization[END_REF], Wei et al., 2021). The implemented algorithm consists in an immersed, hybrid (bulk/interface), mixed (displacement/force,) higher-order, and smoother numerical scheme that appears accurate and scalable to treat multiple inclusions in frictionless contact, frictional contact, and delamination with the matrix of the composite microstructure.

This chapter is organized as follows: after this introduction, the proposed strategy for dealing with both non-linear and non-conformal interfaces is presented in Section 3.2, as well as the implementation of Nitsche's method within the immersed boundary approach. The LaTIn approach is discussed in Section 3.3, with details on different non-linear behaviours. Next, a range of numerical experiments are carried out in 2D to assess the performance of our method with respect to accuracy and efficiency in Section 3.4, followed by conclusions and discussion in Section 3.5.

Isogeometric immersed LaTIn method

As stated in the introduction, we propose to draw inspiration from the IBCM recently introduced in [START_REF] Wei | Immersed boundary-conformal isogeometric method for linear elliptic problems[END_REF]. IBCM proposes to construct a layer of discretization conformal to the boundary whilst a simple background mesh is used for the remaining domain. Thus, it leverages the geometric flexibility of the immersed boundary method and the advantages of a conformal discretization. Notably, in our case, the conformal discretization allows to capture accurately the solution near the non-linear interface. Thus, in the present framework, the idea is to separate the non-linear and nonconformal aspects of the interfaces by making appear two different types of interface (similar mindset as in [START_REF] Guinard | Multiscale analysis of complex aeronautical structures using robust non-intrusive coupling[END_REF] for instance): non-conformal but perfect interfaces (i.e., that ensure the equilibrium and the standard kinematic compatibility between the subdomains) and non-linear but conformal ones (i.e., where the boundaries of the subdomains are aligned between each other). The method is described in this section, starting with the principle, and then focusing on the numerical schemes to address the non-conformal (but perfect) interfaces and finally the non-linear (but conformal) ones.

Reference problem

For the sake of simplicity, let us take here only two subdomains. The reference problem is shown on Figure 3.1. Since our target applications are composite materials such as fibre-reinforced ones, we consider the coupling of a global model (say the matrix) with multiple local models (say the fibres) through non-linear interfaces. Subscripts m and f may be viewed as referring to the matrix and fibre, respectively. Linear elasticity is assumed for the two bodies while the interface is expected to exhibit non- linear behaviour, in particular unilateral or frictional contact, or even delamination. Boundary conditions and applied loads are shown on Figure 3.1. Given that at this step we do not precise the interface behaviour, this problem is the coupling problem discussed in Section 1.4. Thus, ∀i ∈ {m, f }, Eqs. (1.17) should be verified and interface conditions (1.31) read with the use of the current notation:

σ m n m + σ f n f = 0 on Γ ; g(u m , u f , σ m n m , σ f n f ) = 0 on Γ ;
In the following (see Section 3.3), g will be detailed considering the case of contact (frictionless or frictional) and delamination. Let us note that in this step, interface Γ is expected to be both non-conformal and non-linear.

Principle

Starting from interface Γ depicted in Figure 3.1, the idea is to extrude it both in the matrix and fibre directions to obtain two conformal layers. The situation is illustrated in Figure 3.2. As a result, the possibly non-linear interface Γ becomes conformal, while the non-conformal aspect is moved inside subdomains. As for notations, we introduce the superscript j ∈ {l, b} to precise whether we are in the layers (l) or in the bulks (b). Consequently, we denote Ω l m (resp. Ω l f ) the conformal layer within the matrix (resp. fibre) and Ω b m (resp. Ω b f ) the complementary matrix subdomain (resp. fibre subdomain). Finally, the process makes appear two new interfaces, denoted Γ m and Γ f , that allow to couple the layer and the bulk within the matrix and the fibre, respectively. These interfaces, while being non-conformal, are perfect in the sense that there is no presence of non-linear behaviour. Eventually, given the conformal nature of Γ, the non-linear interface behaviour can be enforced in a standard manner. In this work, we utilize the LaTIn method (Section 1.4 to enforce non-linear interface behaviour locally, especially pointwise at each interface Gauss point. Regarding the perfect but non-conformal couplings, we resort to Nitsche's method (see Section 1.4) following [START_REF] Wei | Immersed boundary-conformal isogeometric method for linear elliptic problems[END_REF].

Returning to the formulation, the subdomain-wise Eqs. (1.17) become (see again Figure 3.2: Introducing conformal layers into the initial problem. Now, each subdomain consists of two, Ω l i and Ω b i (i ∈ {m, f }) for the layer and the bulk, respectively. The possible non-linear coupling between the two layers is performed through interface Γ (still in red) that becomes conformal. The non-conformal aspect is moved inside the subdomains: for i ∈ {m, f }, non-conformal but perfect interfaces Γ i appear (depicted in green) to connect Ω l i and Ω b i . The LaTIn method will be applied for the coupling at interface Γ, while a Nitsche coupling scheme will be used at interfaces Γ i , thereby arriving at a robust, simple, and efficient strategy.

j i ) + f g i = 0 in Ω j i ; σ j i = C i ε(u j i ) in Ω j i ; u j i = u g i on Γ u j i ; σ j i n j i = F g i on Γ F j i ; (3.2)
with the following perfect interface conditions: ∀i ∈ {m, f },

σ l i n l i + σ b i n b i = 0 on Γ i ; u l i -u b i = 0 on Γ i . ( 3.3)
And lastly, Eqs. (3.1) now concern only the layers, which yields:

σ l m n l m + σ l f n l f = 0 on Γ ; (3.4a) g(u l m , u l f , σ l m n l m , σ l f n l f ) = 0 on Γ. (3.4b)
The object in what follows is to derive robust and efficient numerical schemes to solve problem (3.2)-( 3.3)-(3.4a)-(3.4b).

Non-conformal coupling with Nitsche's approach

As stated before, we use Nirsche's method to enforce Eqs. (3.2)- (3.3), i.e. the equations related to the perfect and non-conformal couplings across Γ i , ∀i ∈ {m, f }.

CHAPTER 3. IMMERSED BOUNDARY-CONFORMAL ISOGEOMETRIC LATIN METHOD FOR MULTIPLE NON-LINEAR INTERFACES

Nitsche's coupling formulation

Nitsche's approach was discussed in Section 1.2.3. Thus, for the problem in question, the weak formulation (1.20) becomes: ∀i ∈ {m, f }, j∈{l,b}

a j i (u j i , v j i ) - Γi {C i ε(v i )} n l i • [[u i ]]dΓ - Γi [[v i ]] • {σ i } n l i dΓ +ζ i Γi [[v i ]] • [[u i ]]dΓ = j∈{l,b} l j i (v j i ), ∀(v l i , v b i ) ∈ V l i × V b i ; (3.5)
This formulation ensures Eqs. (3.2)- (3.3), thus insuring coupling between the bulk and the layer inside of each model. We remind that in Eq. (3.5), ζ i is a stabilization parameter. Following [START_REF] Antolin | Overlapping multipatch isogeometric method with minimal stabilization[END_REF],we take ζ i as:

ζ i = β i ((h l i ) -1 + (h b i ) -1 ) with β i = 6(p max i ) 2 × 8(E max i /(1 -2ν max i )) ;
where h j i is the maximum element size in Ω j i , and p max i , E max i and ν max i are the maximum polynomial degree, Young modulus, and Poisson ratio, respectively, of the subdomains (Ω j i , ∀j ∈ {l, b}) to be coupled. With such choices, the formulation is rather simple and has proved to be stable and accurate (see again [START_REF] Antolin | Overlapping multipatch isogeometric method with minimal stabilization[END_REF], [START_REF] Wei | Immersed boundary-conformal isogeometric method for linear elliptic problems[END_REF]). Here, given that the layer consists of the same material as the corresponding bulk, we obviously have

E max i = E l i = E b i = E i and ν max i = ν l i = ν b i = ν i grouped in the Hooke tensor C i , ∀i ∈ {m, f }.
Finally, let us condense formulation (3.5) by writing it as follows: ∀i ∈ {m, f }, find

u i ∈ U i such that, a i (u i , v i ) = l i (v i ), ∀v i ∈ V i , (3.6)
where

U i = U l i × U b i , V i = V l i × V b i
, and a i (u i , v i ) and l i (v i ) are the left-hand side and right-hand side of Eq. (3.5), respectively. We recall that Eq. (3.6) enforces Eqs. (3.2)- (3.3). It thus remains to ensure Eqs. (3.4), which is the object of the next Section 3.3.

Implementation care with immersed methods

The implementation of formulation (3.6) requires special care to (i) evaluate integrals over pieces of d-variate elements (for subdomains Ω b i ), (ii) evaluate integrals over (d -1)-variate interfaces that cut the meshes (for computing the Nitsche coupling operators over Γ i ), and (iii) manage the possible ill-conditioning of the stiffness matrices over Ω b i (i.e., the operators associated to bilinear forms a b i ). For point (i), for simplicity and robustness, we use the recursive, quad-tree based quadrature approach employed in FCM [START_REF] Ruess | Weak coupling for isogeometric analysis of non-matching and trimmed multi-patch geometries[END_REF][START_REF] Schillinger | The finite cell method: A review in the context of higherorder structural analysis of CAD and image-based geometric models[END_REF]. This approach was covered in Section 1.2.2. For sure, more geometrically faithful quadrature rules (Antolin et al., 2019b[START_REF] Garhuom | Non-negative moment fitting quadrature for cut finite elements and cells undergoing large deformations[END_REF][START_REF] Kudela | Efficient and accurate numerical quadrature for immersed boundary method[END_REF][START_REF] Legrain | Non-negative moment fitting quadrature rules for fictitious domain methods[END_REF], Nagy and Benson, 2015a) could also be considered since our approach is generic in terms of integration schemes.

Regarding point (ii), we pay attention to discretize interfaces Γ i as the intersection of the two meshes on both sides of it (see Figure 3.3). This ensures that on each interface element, all the involved basis functions from the two subdomains to couple are polynomials rather than piecewise polynomials.

Finally, for point (iii), we start by removing from the corresponding stiffness matrices all the degrees of freedom associated with basis functions whose support does not intersect with the integration domain. We use the criteria based on the basis functions value on the support (see Section 1.2.2). We take into account all contributions Figure 3.3: Interface discretization for coupling non-conformal meshes using the Nitsche approach: to compute integrals over Γ m , the interface discretization is chosen as the intersection of two non-conformal meshes on both sides of it. Interface Γ m is shown, Γ f is of course discretized in a similar manner.

to the physical domain, no mater how small they are (which means to take ε = 0 in Eq. (1.16)). Then, we make use of a diagonal scaling preconditioner presented in Section 1.2.2. With the latter, we never encountered instabilities in our numerical experiments.

Enforcement of non-linear interface behaviours with the LaTIn method

Let us now move to the non-linear interface Γ. The object here is to add Eq. (3.4) to formulation (3.6). In order to do so, we recall that we consider the LaTIn method (see Section 1.4). In what follows, we discuss the implementation of the linear stage for both static and quasi-static frameworks, as well as modelling different non-linear interface behaviour during the local stage.

Implementation of the linear stage

Regarding the linear stage, at iteration n, ∀i ∈ {m, f } Eq. (1.35) should be verified. Since we consider here a conformal interface Γ, this problem should be rather simple to solve and do not require additional stabilization terms as in the case of non-conformal LaTIn interfaces (see, e.g., [START_REF] Claus | A stable and optimally convergent LaTIn-CutFEM algorithm for multiple unilateral contact problems[END_REF]). The stabilization is actually moved to the Nitsche coupling in Section 3.2.3 in our approach. In the discrete version, Eq. (1.35) gives: ∀i ∈ {m, f }, find u

(n) i ∈ U i such that, K * i u i = f * i . (3.7)
where ∀i ∈ {m, f }, K * i is composed of the initial stiffness matrix of subdomain Ω i and additional interface stiffness:

K * i = K i + k - i Γ (N l i ) T N l i dΓ N i Ω→Γ . (3.8)
In the equation above, the operators K i denote the stiffness associated Ω i , thus referring to the stiffness of the matrix or the fibre. This stiffness operator is derived through the discretization of Eq. (3.5), and it encompasses both the conformal layer and the bulk, with their coupling achieved through Nitsche's approach. The matrices N l i represent the standard shape function matrices. Since the integrals are computed on interface Γ, these basis functions correspond to the conformal layer. Finally, operators N i Ω→Γ

CHAPTER 3. IMMERSED BOUNDARY-CONFORMAL ISOGEOMETRIC LATIN METHOD FOR MULTIPLE NON-LINEAR INTERFACES

are matrices containing elements of 0 and 1, designed such that w i = N i Ω→Γ u i , where w i is the vector of displacement values at integration points. Thus, we simply take the existing (d -1)-variate (boundary) basis functions associated to subdomain Ω l i that generate Γ to discretize w i . Let us note that the last term in Eq. (3.8) is actually a mass matrix of interface Γ, and, unlike the stiffness matrix, consists of an integral of basis functions and not their derivatives.

The right hand side in Eq. (3.7) reads:

f * i = f i + Γ (N l i ) T λi dΓ + k + i Γ (N l i ) T ŵi dΓ, ( 3.9)
where f i is the load vector associated to domain Ω i . For integrals in Eq.( 3.9), we actually only need the evaluation of λi and ŵi at the interface Gauss points, which is the case at the end of the local stage. As for the linear step variables, displacement values w i and forces λ i , evaluated at Gauss points are approximated with the corresponding to domain Ω l i basis functions. Remark 3.

3.1

It may be noticed at this stage that interface Γ in the proposed LaTIn strategy actually needs to be geometrically conformal but it can be non-matching. In other words, interface Γ needs to be aligned with the boundary edges of the subdomains to be coupled but the discretizations of these subdomains over Γ (mesh refinements and/or polynomial degrees) can differ (see, e.g., [START_REF] Bouclier | Development of a new, more regular, Mortar method for the coupling of NURBS subdomains within a NURBS patch: Application to a non-intrusive local enrichment of NURBS patches[END_REF] for the nomenclature). With a non-matching interface, only additional care is required to define the integration points over Γ for the local equations: in this case, interface Γ is discretized as the intersection of the two meshes on both sides of it and:

n gp = max i∈{m,f } (p i ) + 1
Gauss points are taken per (d -1)-variate elements of the latter. That being said, let us underline that taking different discretizations for the two layers may appear strange since there is no reason to prioritize the accuracy of the interface fields of one of the subdomains to correctly represent the non-linear interface behaviours of interest.

Treatment of the non-linearity in specific cases

Now let us give a few more details regarding the local stage (1.33). For simplicity, we omit in what follows the superscript (n). Thus, ∀i ∈ {m, f }, starting with quantities (λ i , w i ), the aim is to find the hat quantities ( λi , ŵi ), at each interface Gauss point, such that Eqs. (1.33) and (1.34a) are satisfied. In this work, we consider contact (with and without friction) and delamination (using a cohesive interface) for the non-linear behaviours.

Perfect interface

Before providing details on these non-linear behaviours, let us present a simple case of a perfect interface. Thus, Eq. (1.33) should impose continuity of the displacement and equilibrium of the forces, i.e. force and displacement interface fields ( λi , ŵi ) i∈{m,f } should verify the following behaviour of the interface:

L : ( λi , ŵi ) i∈{m,f } such that    λm + λf = 0 on Γ ŵm -ŵf = 0 on Γ . (3.10)
Then, Eq. (1.34) links the interface quantities with the bulk quantities through search directions.

Thus, Eqs. (3.10) and (1.34a) form a linear system of four equations for four unknowns ( λm , λf , ŵf , ŵm ). This system is solved explicitly at each interface point.

Remark 3.3.2

In the case of the perfect interface, the problem becomes a standard coupling problem considered in Section 1.2.3. Actually, solving local stage equations (3.10) and (1.34a) and linear stage equation (1.35) is equivalent to solving Eqs. (1.17) (subdomain equations) with the interface conditions (1.18).

Quasi-Static formulation

Non-linear behaviour should take into account the loading history, which was not the case for the perfect interface. In order to do so, we actually need to introduce a pseudotime and then write the LaTIn method in a quasi-static framework. Consequently, the above LaTIn algorithm is slightly modified as the search directions now link, at each time step t, the surface forces λ

(t) i , λ (t)
i and the pseudo-velocities ẇ(t) i , ŵ(t) i such that:

k + : λ(t) i -λ (t) i -k + v i ( ŵ(t) i - ẇ(t) i ) = 0 on Γ ; (3.11a) k -: λ(t) i -λ (t) i + k - v i ( ŵ(t) i - ẇ(t) i ) = 0 on Γ. (3.11b)
Since in practice, we take k + = k -, in what follows the superscripts + andare omitted.

Next, a standard implicit scheme is used to express the velocities as a function of the displacements:

ẇ(t) i = w (t) i -w (t-1) i ∆t ; ŵ(t) i = ŵ(t) i - ŵ(t-1) i ∆t . (3.12a)
Finally, contact and cohesive behaviours require the use of a local basis, at each interface Gauss point, defined by the tangent and normal vectors associated with the interface. In the following, the subscript τ (resp. n) refers to the tangential (resp. normal) component of a quantity expressed in this basis. The normal vector is arbitrarily defined from the matrix to the fibre. Here, we also benefit from isogeometric description of the interfaces to be in contact. Indeed, they are intrinsically smooth enough so that the local basis {n, τ } can be defined in any point.

Remark 3.3.3

Note that k v i in Eqs. (3.11a) and (3.11b) does not strictly have the same dimension as k i in Eqs. (1.34a) and (1.34b), due to the introduction of the time increment ∆t. However, let us emphasize that we can take in practice ∆t = 1 since we are in the case of a quasi-static transformation (i.e., without real dynamic effects). Therefore, we take for our numerical experiments the same values for k v i as for k i , ∀i ∈ {m, f }. The specific values of the search directions employed in this work can be found in Section 3.3.3.

CHAPTER 3. IMMERSED BOUNDARY-CONFORMAL ISOGEOMETRIC LATIN METHOD FOR MULTIPLE NON-LINEAR INTERFACES

Frictionless contact

As our first example of non-linear interface behaviour, we will consider frictionless contact on the interface. We introduce j n as the initial normal gap between two subdomains. At each Gauss point and at each time step, the interface quantities should satisfy the following equations:

λ(t) m + λ(t) f = 0, ( Mechanical
C (t) n = ŵ(t-1) mn - ŵ(t-1) fn + j n + ∆t ẇ(t) mn - ẇ(t) fn -∆t k v -1 m λ (t) mn -k v -1 f λ (t)
fn .

Thus,

C (t)
n is a point-wise indicator, which is derived from normal components of the interface quantities. It refers to the contact state of the corresponding interface Gauss point at time step t:

Contact ⇔ C (t) n ≤ 0 No contact ⇔ C (t)
n > 0 The demonstration of these equalities is beyond the scope of this work but can be found, e.g., in [START_REF] Oumaziz | Une méthode de décomposition de domaine mixte non-intrusive pour le calcul parallèle d'assemblages[END_REF].

Once indicator

C (t)
n has been computed, two situations are possible for each Gauss point. In the first scenario, where C (t) n > 0, it signifies the absence of contact at the point. Consequently, there is no contact forces. The equations governing the local stage in this context are as follows (for each time step t):

λ(t) mn = λ(t) fn = 0, λ(t) mτ = λ(t) fτ = 0, ŵ(t) m = ẇ(t) m -k -1 vm λ (t) m , ŵ(t) f = ẇ(t) f -k -1 vf λ (t) f , ŵ(t) m = ∆t ŵ(t) m + ŵ(t-1) m , ŵ(t) f = ∆t ŵ(t) f + ŵ(t-1) f . (3.14)
In the second scenario, when C n ≤ 0, it indicates that the specific Gauss point is in contact, necessitating the computation of contact forces. In this context, the equations governing the local stage for each time step t read:

λ(t) mτ = λ(t) fτ = 0, λ(t) mn = (k -1 m + k -1 f ) -1 C (t) n , ŵ(t) m = ẇ(t) m + k -1 vm ( λ(t) m -λ (t) m ), ŵ(t) f = ẇ(t) f + k -1 vf ( λ(t) f -λ (t) f ), ŵ(t) m = ∆t ŵ(t) m + ŵ(t-1) m , ŵ(t) f = ∆t ŵ(t) f + ŵ(t-1) f .
(3.15) Thus, at each time step, six unknowns are computed for each Gauss point:

( ŵ(t) m , ŵ(t) f , ŵ(t) m , ŵ(t) f , λ(t) m , λ(t) f ).
It is done by explicitly solving the corresponding system of linear equations (Eqs. (3.14) in the case of no contact, Eqs. (3.15) in the case of contact).

Frictional contact

Now, let us expand our discussion to include the case of frictionless contact for a contact interface that follows the Coulomb friction law. We introduce µ, which represents Coulomb's friction coefficient. The interface quantities need to satisfy in the case of friction contact the following equations at each interface Gauss point and time step:

λ(t) m + λ(t) f = 0, (Mechanical equilibrium) ŵ(t) mn - ŵ(t) fn + j n ≥ 0, (Non interpenetration) λ(t) mn - λ(t) fn ≥ 0, (Positive reaction forces) ( λ(t) mn - λ(t) fn )( ŵ(t) mn - ŵ(t) fn + j n ) = 0, (Signorini conditions)    If ∥ λ(t) iτ ∥ < µ| λ(t) in | then ŵ(t) mτ - ŵ(t) fτ = 0, If ∥ λ(t) iτ ∥ = µ| λ(t) in | then ∃α > 0, ŵ(t) mτ - ŵ(t) fτ = -α λ(t) mτ - λ(t) fτ .
(Sticking) (Sliding) (3.16) Thus, the problem of the local stage consists in solving Eqs. (3.16), (3.11a) and (3.12a). Incorporating the frictional aspect is achieved by introducing new conditions in cases involving contact. Initially, the indicator C (t) n is computed. Subsequently, the determination of normal components, as explained in the preceding subsection, is carried out. In situations where contact occurs, an additional indicator G (t) τ is then calculated to evaluate the possibility of sliding. Similar to

C (t) n , this new indicator G (t) τ
is explicitly determined from known fields at time step t:

G (t) τ = λ (t) mτ -k v m ẇ(t) m - k v f k v m + k v f λ (t) mτ + λ (t) fτ -k v m ẇ(t) mτ -k v f ẇ(t) fτ .
The following equalities allow to determine the status of a Gauss point in contact [START_REF] Oumaziz | Une méthode de décomposition de domaine mixte non-intrusive pour le calcul parallèle d'assemblages[END_REF] in the same simple way as for the previous case:

Sticking ⇔ |G τ | ≤ µ | λt+1 mn | Sliding ⇔ |G τ | > µ | λt+1 mn |
In summary, the meaning of these indicators is as follows:

C (t)
n signifies the contact state of the corresponding interface Gauss point at time step t, which can be either in contact (C

(t) n ⩽ 0) or not (C (t)
n > 0). In the case of contact, the second indicator G

λ(t) mτ = G (t) τ , λ(t) mτ + λ(t) fτ = 0, ŵ(t) mτ = ẇ(t) mτ + k -1 vm ( λ(t) mτ -λ (t) mτ ), ŵ(t) fτ = ẇ(t) fτ + k -1 vf ( λ(t) fτ -λ (t) 2τ ), ŵ(t) m = ∆t ŵ(t) m + ŵ(t-1) m , ŵ(t) f = ∆t ŵ(t) f + ŵ(t-1) f .
(3.17 In the sliding case, the tangential forces should be modified and the following equations have to be solved:

λ(t) mτ = µ λ(t) mn G (t) τ |G (t) τ | , λ(t) fτ = µ λ(t) fn G (t) τ |G (t) τ | , ŵ(t) mτ = ẇ(t) mτ + k -1 vm ( λ(t) mτ -λ (t) mτ ), ŵ(t) fτ = ẇ(t) fτ + k -1 vf ( λ(t) fτ -λ (t) fτ ), ŵ(t) m = ∆t ŵ(t) m + ŵ(t-1) m , ŵ(t) f = ∆t ŵ(t) f + ŵ(t-1) f . (3.18)
The frictionless contact interface can be regarded as a specific case of the frictional contact interface. Eqs. (3.15) can be derived from Eqs. (3.18) by setting µ = 0. Thus, there is no tangent forces and Eqs. (3.18) are simplified.

Remark 3.3.4

It should be noted, that the above equations are derived for a 2D case, meaning that there is only one tangential direction. Hence, the computation of tangential forces is simplified in this case. For a 3D case, the expressions are slightly more complicated, see, e.g., [START_REF] Oumaziz | A non-invasive implementation of a mixed domain decomposition method for frictional contact problems[END_REF].

The overall procedure of the implementation of a contact interface (both frictionless and frictional) in the LaTIn method is summarized in Algorithm 1.

Remark 3.3.5

Two frameworks can be used to take into account the quasi-static assumption. The first one is the incremental approach, where the LaTIn method can be viewed as a non-linear solver used to find the solution for one time step before determining the next time step. The second framework is the non-incremental approach. The LaTIn method is seen as a more sophisticated solver since all time steps are embedded in the LaTIn iterations. In this variant, each local stage and linear stage involves finding an approximation of the solution for all the time steps. In contrast to Algorithm 1, this means that the order of the time loop (for each t do) and the iteration loop (while Not converged do) is inverted. This non-incremental approach is particularly adapted to multiparametric studies [START_REF] Roulet | A parallel strategy for the multiparametric analysis of structures with large contact and friction surfaces[END_REF] or time/space decomposition with model order reduction [START_REF] Giacoma | Toward an optimal a priori reduced basis strategy for frictional contact problems with LaTIn solver[END_REF], which is not our case here so we carry out the incremental approach.

Cohesive interface

Let us now specify the equations regarding the cohesive model for the interface and how it is managed within the LaTIn method.The fundamental concepts of cohesive zone models were initially introduced in [START_REF] Barenblatt | The formation of equilibrium cracks during brittle fracture. General ideas and hypotheses. Axially-symmetric cracks[END_REF], [START_REF] Dugdale | Yielding of steel sheets containing slits[END_REF] and [START_REF] Barenblatt | The mathematical theory of equilibrium cracks in brittle fracture[END_REF]. Cohesive models are based on the idea that the interface is defined by a stiffness that links the surface force to the normal displacement jump. This traction-separation relationship is called the cohesive law. Various cohesive laws have been proposed in the literature, see e.g., [START_REF] Shet | Analysis of energy balance when using cohesive zone models to simulate fracture processes[END_REF] for a compilation of some of the 2015) for a review in the context of the numerical simulation of matrix cracks and delamination. Delamination can be viewed as the gradual degradation of the interface's stiffness. Initially, the cohesive interface is undamaged, but during loading, the interface can be progressively damaged until it breaks completely. A fully damaged interface can be considered as a contact interface.

In this part, we assume that the interface is opening and the normal displacement jump is positive. The case of a closed interface and a null displacement jump is considered as contact behaviour in the normal direction and is enforced as we described previously. Therefore, to avoid cumbersome notations, we do not specify the positive part of the displacement jump in the following.

In general terms, the constitutive law can be written as follows:

λ(t) fn = k [[ ŵ(s) n ]] [[ ŵ(t) n ]], with s ≤ t,
where k is the stiffness of the interface which depends on the normal displacement jump [[ ŵ(t) n ]] = ŵ(t) fn -ŵ(t) mn . In this case, the complete problem to be solved at the local stage consists in finding a solution that verifies the above constitutive law, the equilibrium of the interface, and the search direction equations, which read: f . As a result, the computation of the surface forces are fully explicit and there is no need for a non-linear solver. We advise the interested reader to consult [START_REF] Fernandez | The effect of the mesh refinement on a multiscale domain decomposition method for the non-linear simulation of composite structures[END_REF], [START_REF] Kerfriden | A three-scale domain decomposition method for the 3D analysis of debonding in laminates[END_REF], [START_REF] Violeau | Micromodel-based simulations for laminated composites[END_REF] for more details on this topic.

λ(t) fn = k [[ ŵ(s) n ]] [[ ŵ(t) n ]], s ≤ t, (Constitutive law) λ(t) fn + λ(t) mn = 0, (Equilibrium of surface forces) λ(t) m -λ (t) m -k v m ( ŵ(t) m -ẇ(t) m ) = 0, (Search direction for the matrix) λ(t) f -λ (t) f -k v f ( ŵ(t) f - ẇ(t) f ) = 0, (

Remark 3.3.6

So far we have written the problem for the normal component of the interface quantities. We can define the same kind of equation for the tangential ones. The main difference is that compressive behaviour does not make much sense for the tangential component, shear damage will occur in any direction in the tangential plane. Therefore, the tangential displacement jump is considered as an absolute value.

Bilinear law

To describe the cohesive model, we consider a simple bilinear law [START_REF] Alfano | Finite element interface models for the delamination analysis of laminated composites: mechanical and computational issues[END_REF]), as we applied in Section 2. 3.3 (see Figure 2.13). However, this time it links interface hat variables λn and [[ ŵn ]]. As mentioned in Section 2.3.3 this law requires three parameters (energy restitution rate G c , critical stress σ c , penalisation coefficient p) and δ 0 = p • 2 Gc σc and δ c = 2 Gc σc (displacement jumps corresponding to the initiation of the damage and to the complete damage, respectively).

To introduce the implementation of the bilinear cohesive law, we start by defining the maximal displacement jump both for normal and tangential components as follows:

[[ ŵ(t) n ]] max = max s≤t {[[ ŵ(s) n ]], 0}; [[ ŵ(t) τ ]] max = max s≤t { [[ ŵ(s) τ ]] , 0}.
Conformally to what was written above (see Remark 3.3.2), only positive normal displacement jumps are considered in the model, while for tangential displacement, the absolute value is utilized. This approach accounts for shear damage which can occur in any direction.

Then, depending on the values of [[ ŵ(t) n ]] max and [[ ŵ(t) τ ]] max , three scenarios are possible leading to three different expressions for the interface forces λ:

∀ i ∈ {n, τ } :        if [[ ŵ(t) i ]] max ≤ δ i0 , λi = k i0 [[ ŵ(t) i ]]; if δ i0 ≤ [[ ŵ(t) i ]] max < δ ic , λi = k i0 (1 -d i )[[ ŵ(t) i ]]; if [[ ŵ(t) i ]] max ≥ δ ic , λi = 0.
In the equation above, k n0 (resp. k τ0 ) is the initial stiffness associated to the normal (resp. tangential) direction. Parameters δ n0 (resp. δ τ0 ) and δ nc (resp. δ τc ) are the displacement corresponding to the initiation and completion of damage in normal (resp. tangential) direction. Parameters d n and d τ are the damage parameter and are defined as follows:

∀ i ∈ {n, τ } : d i = δ ic δ ic -δ i0 [[ ŵ(t) i ]] max -δ i0 [[ ŵ(t) i ]] max .
The minimum value of d, d = 0, means that the interface is undamaged. The maximum value of d, d = 1 corresponds to a fully damaged interface.

Linear stage under a quasi-static assumption

Finally, let us underline that the linear stage Eq. (1.35) is slightly modified to take into account the quasi-static assumption. To do so, we simply need to repeat the procedure to obtain Eq. (1.35), but with Eq. (3.11) as search directions. Combining this with the time integration scheme Eq. (3.12), the problem at the linear stage now reads: ∀i ∈ {m, f }, at time step t find u (t) i ∈ U i such that:

a i (u (t) i , v i ) + Γ 1 ∆t k i u (t) i • v i dΓ = l (t) i (v i ) + Γ λ(t) i + k i ŵ(t) i + 1 ∆t k i w (t-1) i • v i dΓ, ∀v i ∈ V i .
(3.20)

Search directions

The choice of search directions significantly impacts the convergence of the method. Ideally, the optimal search direction would precisely represent the behaviour of the complement of the subdomain. However, constructing such an operator is unfeasible in real applications. First, the complement of the subdomain is of large size, and therefore, the computation of this operator comes with a huge cost. Second, in a parallel context, subdomains generally do not have access to information about other subdomains. Indeed, as at can be seen in Eq. (1.35) and Eq. (3.20), the subdomains exchange information only be the means of Robin boundary conditions on the interface.

CHAPTER 3. IMMERSED BOUNDARY-CONFORMAL ISOGEOMETRIC LATIN METHOD FOR MULTIPLE NON-LINEAR INTERFACES

Unfortunately, the computation of such an operator would require significant global communication between the subdomains. An approximation must therefore be made. Traditionally, various works on the LaTIn method propose a straightforward choice of search direction operator. In 2D, we take

k m = E f L f and k f = Em
Lm , where E f (resp. E m ) and L f (resp. L m ) are the Young modulus and characteristic length of subdomain Ω f (resp. Ω m ). This choice represents a stiffness of the structure and aims to simply represent the overall behaviour of the structure.

Remark 3.3.7

Even though this classical choice of search direction is rarely questioned, it is noteworthy that alternative selections have been documented in the literature. For example, in [START_REF] Saavedra | An enhanced nonlinear multi-scale strategy for the simulation of buckling and delamination on 3d composite plates[END_REF] a study of the search direction for buckling phenomena in slender structures is proposed.

Numerical results

To assess the performance of the developed isogeometric IBCM LaTIn scheme, we now present a series of numerical experiments in 2D that covers different geometries, discretizations and non-linear interface scenarios. First of all, we start by considering well-known Kirsch problem to investigate the influence of the thickness of conformal layers in IBCM. Then, we investigate 4 test cases: (i) one inclusion in unilateral contact with the matrix, (iii) several inclusions in unilateral contact with the matrix, (ii) a twobody assembly with frictional contact inside, and (iv) one inclusion in delamination with the matrix.

On the thickness of conformal layers in IBCM

While dealing with Immersed Boundary-Conformal methods the question of the right choice for the shape of the conformal layers arises naturally. In a related study [START_REF] Wei | Immersed boundary-conformal isogeometric method for linear elliptic problems[END_REF], the authors, while acknowledging that this question remains open, provide readers with two important ideas: (i) IBCM reveals more accurate results compared to pure immersed boundary approach for all conformal layers configurations considered, (ii) among these configurations, the thinnest conformal layer tends to yield the most accurate results, assuming the same discretization. The second point has a simple explication: indeed, the thinnest layer has the finest mesh, therefore, the solution inside the layer is captured with higher precision and accuracy.

However, we would like to note that these results are relevant to the conformal layer rather than to the whole model. Assuming the conformal layer is built around a zone of interest, where the solution can be less regular or even singular. When reducing the thickness of the conformal layer (keeping the same discretization), two components come into play. First, elements in the conformal layer become smaller, which leads to increased accuracy within the layer. Second, the geometrical domain of the background mesh extends closer to the region of interest. Consequently, this potentially singular solution near the region of interest is captured by the pure immersed boundary method, which results in an increased error in the background mesh.

To demonstrate this idea, we considered the well-known Kirsch problem, which consists of an infinite plate with a circular hole under in-plane tension. The problem details are presented on Figure 3.4(a). Since the problem is symmetric, we considered only a quarter of the region with the appropriate symmetry boundary conditions on the lower and left borders. The exact traction was applied on the upper and right borders. The material in consideration was linear elastic, and we assumed small deformations. The analytical solution is well-known and can be found, for example, in [START_REF] Sadd | Elasticity: theory, applications, and numerics[END_REF]:

σ r = T x 2 (1 - a 2 r 2 ) + T x 2 (1 + 3a 4 r 4 - 4a 2 r 2 ) cos(2θ), σ θ = T x 2 (1 + a 2 r 2 ) - T x 2 (1 + 3a 4 r 4 ) cos(2θ), σ rθ = - T x 2 (1 - 3a 4 r 4 +
2a 2 r 2 ) sin(2θ), where T x is the applied traction, a is the radius of the hole, r and θ are polar coordinates (see again Figure 3.4(a)).

In this numerical example, we considered IBCM without applying the LaTIn method. Thus, the problem consisted in coupling the conformal layer with the bulk by Nitsche's approach (Eq. (3.5)). The background mesh contained the whole square domain (a square L × L), whilst the conformal layer's thickness t differed to show its influence on the error. For the background mesh integration within the immersed boundary approach, the predefined maximum level of decomposition l max = 8 was used. We considered three cases: t = 0.2a, t = 0.6a, and t = a. See Figure 3.4(b) for a discretization example.

For each of this models we considered six different discretizations: the coarsest one contained 5 × 5 elements in the background mesh and 2 × 2 elements in the conformal layer (as shown in Figure 3.4(b)), each following mesh was obtained from the previous one by doubling the number of elements in each direction. We used the same discretization for all geometries of the conformal. Splines of degree 2 were used.

As for error estimator, we took the energy norm, which is defined as: that this technique allows us to circumvent the difficulty of dealing with a coupling interface both non-conformal and incorporating contact.

Err = ∥u h -u ref ∥ ∥u ref ∥ = Ω ε(u h -u ref ) : Cε(u h -u ref ) dΩ Ω ε(u ref ) : Cε(u ref ) dΩ
The following results were obtained for a rather fine mesh containing 32×32 NURBS quadratic elements in each layer and 80 × 80 B-spline quadratic elements in the matrix background mesh. We chose the same element size for the inclusion's background mesh. The search directions for the LaTIn method, due to the significant contrast between the properties of the materials, were chosen as the Young modulus of the complement divided by its characteristic size (see Section 3.3). The LaTIn iterations were stopped once indicator η reaches 10 -5 . All the shown results correspond to converged LaTIn solutions. Figure 3.8 shows the obtained displacement on a deformed geometry with a scale factor of 10. Figure 3.9 depicts stresses on the deformed geometry. We observe that since σ yy is continuous on the top of the inclusion, materials are in contact in this area, while there is no contact on the right of the inclusion. In addition, Figure 3.10 (left) shows the status of the integration points which is defined by the sign of the contact indicator C (t) n . Normal tractions on the interface Γ are depicted on Figure 3.10 (right) where angle α is the angle between the position vector of a point and the positive direction of x axis. The absence of oscillations in these results numerically accounts for the stability of the numerical scheme.

Then, the convergence of our method with the mesh refinement was studied. For this purpose, five meshes were considered, see Figure 3.11. In addition, NURBS of degrees 2 and 3 were used. Finer meshes were obtained from the initial one in a hierarchical manner by consequently doubling the number of elements in each direction. As a reference solution, we took a very refined one following our approach (320 × 320 cubic elements in background meshes and 128 × 128 cubic elements in the annuli). We investigated the convergence rate in the energy norm, which is defined as:

Err = ∥u h -u ref ∥ ∥u ref ∥ = i∈{m,f } Ωi ε(u h i -u ref i ) : C i ε(u h i -u ref i ) dΩ i∈{m,f } Ωi ε(u ref i ) : C i ε(u ref i ) dΩ
, where u h and u ref denote the solution of interest and the reference solution, respectively. For the evaluation of the error, the coarser solutions were interpolated at the integration points of the reference solution. We show the convergence plots for 2and 3-degree splines on Figure 3.12a with respect to size indicator h. It can be observed that when evaluating the error over the whole domain the rate of convergence for the quadratic solution is almost the one encountered in linear elasticity (h 2 ). This demonstrates the superior accuracy of our method with respect to the state-of-the-art of the LaTIn technique that usually restricts to linear FEA. Nevertheless, it can also be noticed that going to cubic splines does not improve much the convergence. This is expected due to the lack of regularity of the contact problem [START_REF] Kikuchi | Contact problems in elasticity: a study of variational inequalities and finite element methods[END_REF]. For completeness on this point, we finally computed the error only in the matrix background mesh. The associated convergence curves are given in Figure 3.12b. This time, we did obtain the convergence rates usually reached in linear elasticity for both the quadratic and cubic solutions. The solution in the matrix background mesh being quite distant from the contact interface is actually less affected by the lack of regularity of the problem.

Remark 3.4.1

This mesh refinement study suggests that it might be worthwhile to discretize the layers with splines of relatively low polynomial degrees (say 2) and to use splines of higher polynomial degrees (say 3 or more) for the bulk. This would allow to correctly catch the non-linear behaviour at the interface while increasing the per-degree-of-freedom accuracy in the bulk. Note that this treatment is possible with our method without any effort, since it relies on a stabilized Nitsche scheme that allows to couple spline immersed subdomains of different degrees.

Matrix with multiple inclusions

We are now interested in the convergence of the LaTIn method for different number of inclusions. To this end, we considered models consisting of an 8mm × 8mm squareshaped matrix with 1, 2, 4, or 16 (uniformly distributed) inclusions, see Figure 3.13. Each inclusion was of a 1mm diameter, each annulus thickness was 0.1mm. We apply symmetry boundary conditions on the left and bottom edges of the matrix and a horizontal displacement on the right edge. As in the previous case, both matrix and inclusions materials were linear elastic and we assumed again plane stress state and small deformations. We considered the same material for all inclusions. Material parameters were the same as in Section 3.4.2, so that there was a significant contrast in the matrix and inclusions Young moduli. Furthermore, we considered once again frictionless contact behaviour on all matrix/inclusion interfaces. For the numerical simulations, we took rather fine meshes. Matrix background mesh contained 32 × 32 elements, each of its conformal layers being composed of 16 × 4 elements. All inclusions had identical meshes but, contrary to the previous test case, the inclusions' background meshes were eight times finer than the matrix one (see again Figure 3.13 for meshes representation).

Figure 3.14 shows stress distribution for the 4-inclusion problem on a deformed geometry with a scale factor of 10. As can be seen from the continuity of stresses σ yy the highest and the lowest parts of inclusions are in contact with the matrix while some detachments appear on the left and right parts of the inclusions. This seems consistent with the Poisson effect of the plate. For completeness on this point, Figure 3.15 depicts the contact status of interface integration points for this case.

Eventually, the convergence of the LaTIn method for 1, 2, 4, and 16 inclusions problems was studied. Figure 3.16 shows the corresponding convergence curves. It can be observed that the convergence of the developed algorithm is almost independent of the number of inclusions. More precisely, the algorithm needs only about 30 iterations for different numbers of inclusions to attain an indicator η of 10 -5 . Recalling that the different interface and subdomain computations during the local and global stages, respectively, can be carried out in parallel, our approach can thus be used as an efficient non-linear domain decomposition solver. The matrix actually plays the role of the coarse problem of parallel domain decomposition approaches [START_REF] Bosy | A domain decomposition method for isogeometric multi-patch problems with inexact local solvers[END_REF][START_REF] Hirschler | A dual domain decomposition algorithm for the analysis of non-conforming isogeometric Kirchhoff-Love shells[END_REF][START_REF] Widlund | BDDC deluxe algorithms for twodimensional H (curl) isogeometric analysis[END_REF]. It allows to transmit directly the information all over the inclusions, which makes the algorithm scalable.

Frictional contact problem

In the next step, we applied our method to a model with frictional contact at interfaces. To this end, a model with two curvilinear interfaces was created (see Figure 3.17). This model was inspired by [START_REF] Boucard | A suitable computational strategy for the parametric analysis of problems with multiple contact[END_REF] where interfaces were straight.

The model comprised two domains, Ω 1 and Ω 2 , where Ω 1 , analogous to the matrix model in the previous examples, was itself composed of 2 parts, physical and fictitious ones. Two conformal layers were built on curvilinear non-linear interfaces 1 and 2 to treat the non-conformal aspect of domain Ω 1 . Domain Ω 2 , in its turn, was a simple solely physical domain that had a conformal interface discretization with interfaces 1 and 2. Frictional contact between Ω 1 and Ω 2 was modeled by the Coulomb friction law with friction coefficients µ 1 and µ 2 on interfaces 1 and 2, respectively. In addition, interface 3 imposed boundary conditions on the right boundary of Ω 2 . It was a frictionless contact interface with initial gap j between the interface and domain Ω 2 . During the linear stage, the displacement w corresponding to interface 3 was forced to be zero as if there were a rigid body. Then, interface forces λ were calculated from w with the use of the search direction. Here, both materials are linear elastic with equal properties. The discretization used to solve the problem involved 16 × 48 quadratic B-spline elements for Ω 1 , 32 × 2 quadratic B-spline elements for the conformal layers, and 32 × 32 quadratic B-spline elements for Ω 2 (see again Figure 3.17).

The problem studied in this section is quasi-static, and it involved two load steps. In the first step, at time t 0 , a vertical load F 1 was applied on the top boundary of domain Ω 1 , causing interfaces 1 and 2 to come into contact and thus resulting in non-zero contact forces on the whole of these interfaces. During this step, there was no external load on Ω 2 , so F 2 = 0. In step two, at t = t 1 , load F 2 = F max 2 was applied on domain Ω 2 while F 1 is being kept constant. The aim was to investigate the effect of friction coefficients on the reaction forces at interface 3. We varied µ 1 and µ 2 from 0 to 0.6 with 20 values for each coefficient. The model parameters are listed in Table 3.1.

In this situation of frictional contact, we considered that the algorithm had converged when the criterion η < 10 -7 was achieved. The required number of iterations varies depending on the time step and on the values of µ 1 and µ 2 . However, it is of the order of a hundred iterations.

The contact zone at interface 3 can vary based on the chosen friction coefficients. It could be fully in contact, partially in contact, or not in contact. We computed the resulting reaction force R 3 on this interface to quantify the contact status. Results are depicted in Figure 3.18, where we can see 2 main zones. A zone for which the friction coefficients are too high and prevent contact. It results in a null reaction force (see zone close to point (d) in Figure 3.18). In the second zone, a partial or full contact is reached on interface 3: the size of the contact area depends on the friction coefficient, the smaller the friction coefficients, the higher the reaction force (see zone close to point (a) in Figure 3.18).

Remark 3.4.2

It can be seen that contrary to the model with straight interfaces from [START_REF] Boucard | A suitable computational strategy for the parametric analysis of problems with multiple contact[END_REF], the reaction force for µ 1 = µ 2 = 0 is bigger than the applied horizontal load F 2 . Indeed, when an interface is in contact, it results in nonzero normal forces. For the curvilinear interfaces, these normal forces have a nonzero horizontal component, which gives additional horizontal load on domain Ω 2 .

To give a closer look at the results, Figure 3.19 shows horizontal displacement on the deformed configuration for four particular pairs of (µ 1 , µ 2 ). It can be seen that interface 3 becomes completely in contact if there is no friction on interfaces 1 and 2 (Figure 3.19a). Then, Figure 3.19b and Figure 3.19c depict two cases where interface 3 is only partially in contact. Finally, for µ 1 = µ 2 = 0.6 (Figure 3.19d) there is no contact on interface 3.

Matrix with one inclusion and delamination

Finally, we investigated the case of a cohesive inclusion/matrix interface with the last example. The problem description is provided in Figure 3.20(a), and we took the same material properties as in the first example (see Section 3.4.2) for the matrix and inclusion. Regarding discretization, the matrix contained 32 × 16 NURBS elements in the conformal layer and 32 × 32 B-spline elements in the background mesh. For the inclusion, we opted for a mesh two times finer than for the matrix, both in the bulk and in the layer. Additionally, we employed different splines degrees: p = 2 for the inclusion and p = 4 for the matrix (see Figure 3.20(b)). Therefore, in view of assessing the robustness of the developed approach (and especially regarding the LaTIn strategy), we placed ourselves in the case of a geometrically conformal but non-matching interface here (see remark 3.3.1). Dirichlet boundary conditions were imposed on the top and bottom of the region. The simulation involved 40 time steps with linear increments, resulting in a final displacement of u g = 0.2mm. To simulate the delamination of the matrix/fibre interface, we employed the standard bi-linear law described in sec. 3.3.2. In our computations, we set G c = 0.5 mm•MPa, and σ c = 9 MPa, ρ = 0.1.

The reaction force at the bottom of the domain is shown in Figure 3.21 at each time step. More precisely, Figure 3.21 shows three configurations: (a) the case of an undamaged structure, for which the behaviour of the interface remains in the first linear part (see again Figure 2.13); (b) the structure with the cohesive interface starting from an undamaged one and (c) the structure with a fully damaged interface. As expected, the behaviour of the standard structure falls between the two extreme configurations. It starts with an undamaged interface to reach the behaviour of a fully damaged interface as the load increases. In fact, damage to the interface progressively reduces the overall stiffness of the structure.

Finally, Figure 3.22 presents the damage parameter at the integration points of the interface at the last time step. A parameter value of zero indicates an undamaged interface, while a value of one corresponds to complete interface damage. As it can be seen, the interface is fully damaged on the top and at the bottom, where jumps in normal displacement are maximum, while it stays intact on the left and on the right, where these jumps are insignificant. This can also be seen in Figure 3.23, as the stress at the top and bottom of the inclusion is zero. Since the interface is mostly broken, the inclusion is hardly loaded in the y-direction.

Conclusion

In this chapter, we addressed the challenge of modelling complex structures characterized by the presence of non-linear interfaces. The proposed approach fully relies on IGA, thus benefiting from its superior behaviour. IGA was used for the geometric representation of the different interfaces and for the problem discretization. To offer geometric flexibility, the immersed boundary-conformal method was applied with Nitsche's approach used to couple conformal layer and the bulk. Finally, LaTIn technique was implemented in this context to obtain a parallel domain decomposition algorithm capable of treating any non-linear interface behaviour. In other words, the key point was to separate the difficulty of treating interfaces both non-conformal and non-linear, so as to obtain a simple numerical scheme which stability is ensured within the bulk equations through Nitsche's coupling.

Regarding the non-conformal coupling by Nitsche's approach, the issue of building coupling operators for non-conformal meshes was addressed by special interface discretization. Furthermore, we detailed the implementation of the LaTIn method for addressing three distinct cases of non-linear interface behaviour: frictionless contact, frictional contact, and delamination.

A series of 2D numerical examples was considered to validate the strategy. Firstly, before applying the full strategy, we addressed the question of construction of the conformal layer for immersed boundary-conformal methods and some insights were provided on this challenge, investigating the error independently in the conformal layer and the bulk. Then, a matrix with one central inclusion and frictionless contact on the matrix/inclusion interface was considered. Mesh convergence study was performed for this case. Furthermore, to show the scalability of our approach, convergence in terms of iterations of the LaTIn method was investigated for models with different numbers of inclusions. Finally, numerical examples with frictional contact and delamination were presented.

General conclusions

Final Overview

This thesis represents the search for robust and efficient methods for modelling composite structures at the microstructure level, paving the way for more advanced and accurate simulations. The main focus is on the computation of image-based models of such materials. This models account for microstructure of a material and simulations represent different challenges: the models are associated with big data, they call for complex behaviour laws, there is a great difference in material properties of the matrix and fibres, which may affect the numerical stability. The general objective of the thesis is to respond to the above scientific challenge in order to end up with an efficient digital environment to simulate, in non-linear, heterogeneous microstructures obtained by imaging. We suppose that image processing step is completed, and therefore, we have a geometric description of the model at our disposal. In this context, two strategies were proposed: a fully non-invasive global-IGA/local-FEA coupling and a fully IsoGeometric LaTIn method.

The first strategy, presented in Chapter 2, is suitable for a general case of localized non-linearities presented in a global model. In the context of fibre-reinforced composite materials, non-invasive global/local coupling allows to easily and efficiently simulate the presence of local fibres in the global matrix. The use of IGA on the global scale and FEA on the local scale allows to benefit from superior behaviour of IGA for modelling regular solutions and geometric flexibility and versatility of FEA for modelling local non-linear phenomena. Importantly, the proposed strategy is non-invasive not only in terms of coupling approach, but also in terms of computation of stiffness operators. To achieve it, we resorted to the existing FEA-to-IGA bridge, based on Bézier or Lagrange extraction. This strategy offers great flexibility because it enables the use of any robust FE code capable of modelling complex local behaviours without major modifications. Depending on the level of integration required, our implementation is adaptable to various programming environments. Users might already possess an IG code and wish to couple it with a specific FE software to tackle complex local phenomena, or they may have access only to FE packages. Finally, it's important to highlight that this strategy, although primarily focused on composite structures, is not limited to this domain. We conducted numerical examples related to assembly modelling to demonstrate the versatility of the approache. This strategy can find its application in different domains where accurate modelling of local non-linearities within the bulk is essential.

The second strategy, as discussed in Chapter 3, focuses on a fully IGA approach tailored for simulating composite microstructures comprising multiple inclusions connected to the matrix through non-linear interfaces. This problem can be seen as a particular case of the previous one: here, non-linearities are localized specifically on interfaces and not within the bulk. This strategy addresses the core challenge of dealing with non-linear and non-conformal interfaces. The strategy separates these two challenges: non-conforming coupling is shifted in the bulk using the immersed boundaryconformal method, while the non-linear aspect is handled through the LaTIn method. 100 GENERAL CONCLUSIONS Thus, the non-conforming coupling aspect is achieved using Nitsche's approach, coupling specially constructed conformal layers with the bulk. On the other hand, the LaTIn method provides a framework to model various non-linear behaviours at conformal interfaces, including contact (both frictionless and frictional) and delamination. The implementation details of the LaTIn method for these cases are thoroughly explained. An advantage of the LaTIn method is its suitability for parallel computation, which is especially valuable for large image-based composite models. Consequently, this strategy leverages the superior behaviour of IGA, offers geometric flexibility through the immersed boundary approach, efficiently implements non-linear behaviour, and allows for potential parallelization, thanks to the LaTIn method.

Prospects

Several prospects for this work can be considered. Firstly, let us note that the immersed boundary-conformal IsoGeometric LaTIn method proposed in Chapter 3 may be extended to the 3D case. In this respect, the main issue may concern the surface integration over the non-conformal interfaces for Nitsche's coupling, which would require computing the intersection of a spline surface mesh with the interior of a spline volume mesh [START_REF] Antolin | Region extraction in mesh intersection[END_REF]. Moreover, one will likely need to use a more efficient quadrature rule to build stiffness matrix K, such as moment fitting [START_REF] Abedian | Performance of different integration schemes in facing discontinuities in the finite cell method[END_REF][START_REF] Garhuom | Non-negative moment fitting quadrature for cut finite elements and cells undergoing large deformations[END_REF][START_REF] Legrain | Non-negative moment fitting quadrature rules for fictitious domain methods[END_REF] or decomposition-based techniques [START_REF] Antolin | Robust numerical integration on curved polyhedra based on folded decompositions[END_REF], Wei et al., 2021). Finally, in the 3D scenario, given the size of the problem, parallelization becomes a necessity. In this regard, one can consider utilizing MPI with PETSc [START_REF] Balay | PETSc, the portable, extensible toolkit for scientific computation[END_REF].

Secondly, it might be worthwhile to extend the approach to non-linear material behaviour. Indeed, the matrix material in composites is typically elastoplastic. Therefore, we could consider domain decomposition in the context of non-linear problems using the global/local non-intrusive coupling algorithm. The concept would involve dividing the matrix into a set of non-overlapping local models [START_REF] Duval | Non-intrusive coupling: recent advances and scalable nonlinear domain decomposition[END_REF][START_REF] Kerim | Asynchronous global-local non-invasive coupling for linear elliptic problems[END_REF]. Also, fictitious domain decomposition may be interesting, with the use of the global/local non-intrusive coupling based on Nitsche's technique [START_REF] Bouclier | A Nitsche-based non-intrusive coupling strategy for global/local isogeometric structural analysis[END_REF].

Another prospect consist in extending the strategies to mixed-dimensional problems. In the context of composites, beam models may be used to represent fibres in 3D matrix [START_REF] Kerfriden | A mixed-dimensional CutFEM methodology for the simulation of fibre-reinforced composites[END_REF]. Indeed, a 3D model is not always required, since different 1D models are adopted in the framework of IGA [START_REF] Gan | An isogeometric approach to beam structures[END_REF]. This expansion could also broaden the applicability of the approach, especially in civil engineering and related fields.

Finally, it would be interesting to apply the proposed strategies to real image-based models. Starting from images, e.g. obtained with X-ray computed tomography, and making use of specific image-based spline procedures [START_REF] Passieux | Image-based isogeometric twins of lattices with virtual image correlation for varying cross-section beams[END_REF], it would be attractive to take into account the real inner geometries of materials. This idea can find its application also in the biomedical domain, where patient-specific geometry should be considered [START_REF] Claus | Contact modeling from images using cut finite element solvers[END_REF].

In summary, the potential for extending and applying these methods to various domains is substantial, offering opportunities for further research and application.

A. Résumé des travaux (en Français)

Les matériaux composites offrent une large gamme d'applications dans diverses industries en raison de leurs propriétés avantageuses. Avec le progrès des techniques d'imagerie 3D, la création de modèles mécaniques à micro-échelle devient possible, ouvrant la voie à des simulations plus avancées et plus précises ainsi qu'à des applications telles que l'analyse d'homogénéisation et l'assimilation de données pour les jumeaux numériques. Cependant, le calcul de modèles basés sur des images pour de véritables applications mécaniques présente d'importants défis. Ces modèles impliquent une quantité massive de données en raison du grand nombre de voxels dans les images 3D, ce qui conduit à un grand nombre de degrés de liberté (DDL) nécessaires pour une représentation précise. De plus, il est nécessaire de prendre en compte un comportement de matériau complexe, en particulier aux interfaces entre différentes sous-domaines, telles que les zones de contact ou le délaminage. Des non-linéarités peuvent également exister au sein des sous-domaines, ce qui complique davantage le processus de modélisation et de calcul. Enfin, la différence substantielle de propriétés méchaniques entre la matrice et les fibres nécessite des ressources considérables et des techniques numériques avancées pour des simulations numériques réalistes.

Dans ce contexte spécifique, l'objectif général de cette thèse est de répondre aux défis scientifiques ci-dessus afin d'aboutir à un cadre de travail numérique efficace pour simuler le comportement, y compris non-linéaire, des microstructures hétérogènes obtenues par imagerie. Dans ce travail, nous supposons que l'étape de traitement de l'image est terminée et que nous disposons donc d'une description géométrique du modèle.

Le premier chapitre, à travers une revue de la littérature, présente des stratégies de simulation avancées pour la modélisation de structures complexes qui sont utilisés dans ces travaux. D'abord, nous avons abordé la question de la discrétisation d'un domaine de calcul. L'analyse IsoGéométrique (IsoGeometric Analysis, IGA) vise à unifier les domaines de la Conception Assistée par Ordinateur (CAO) et de l'Analyse par Éléments Finis (Finite Elements Analysis, FEA). L'idée centrale est d'utiliser des fonction de forme lisses et de plus d'ordre plus élevé, telles que les fonctions B-Splines et B-Splines Rationnelles Non Uniformes (Non-Uniform Rational B-Splines, NURBS), à la fois pour la représentation de la géométrie en CAO et pour l'approximation des champs de solution en FEA. Bien que l'objectif initial était de simplifier le processus chronophage de création de modèles mécaniques à partir de modèles de CAO, l'utilisation de splines présente ses propres avantages du point de vue de l'analyse. En effet, les fonctions splines peuvent être C (p-1) régulières entre les éléments pour un degré polynomial p, tandis que les polynômes de Lagrange, qui sont utilisés dans la FEA standard, n'atteignent qu'une régularité C 0 à ces emplacements. Par conséquent, l'IGA présente une précision supérieure par DDL par rapport à la FEA standard pour les simulations mécaniques. Cela rend cette approche particulièrement adapté au calcul à haute performance (High-Performance Computational, HPC). l'IGA peut être interprétée comme une projection de la FEA sur une base réduite spécifique. En effet, les espaces d'approximation potentiellement plus réguliers issus de l'IGA sont inclus dans les espaces de la FEA classique de C 0 , à condition que l'IGA et la FEA aient le même nombre d'éléments et le même degré de polynôme. Dans le contexte de cette relation, il existe une approche de transition entre la FEA et l'IGA, qui repose sur l'extraction de Bézier et l'extraction de Lagrange. Cette technique permet de générer un maillage d'Éléments Finis à partir d'un maillage IsoGéométrique en utilisant différents opérateurs globaux. De même, un lien peut être établi pour les matrices de rigidité, les vecteurs de force et, en fin de compte, le champ de déplacement.

Pour modéliser des géométrie complexe avec IGA, l'approche de la frontière immergée est utilisée. Contrairement aux techniques classiques de résolution des équations aux dérivées partielles (EDP) qui utilisent des maillages pour représenter suffisamment précisément la géométrie du domaine, l'approche de la frontière immergée ne demande pas au maillage d'être conforme à la géométrie du domaine. L'idée centrale consiste à immerger toute géométrie complexe dans un domaine plus simple (généralement rectiligne) et à représenter la solution du problème physique à l'aide des DDL liés au maillage cartésien. L'avantage principal de cette méthode est qu'elle simplifie l'étape de maillage, ce qui facilite l'analyse de problèmes géométriquement complexes et évolutifs. Cependant, pour obtenir une solution approximative satisfaisante dans le domaine physique, il est essentiel de tenir compte avec précision des frontières physiques. Pour ce faire, des procédures d'intégration spécifiques sont développées. Ainsi, l'étape de prise en compte des limites physiques du domaine considéré passe du maillage à l'intégration. Vu que plusieurs sous-domaines doivent être couplé, les approches de couplage faibles sont également présentés.

Après avoir abordé la question de la discrétisation, différentes approches pour résoudre le problème mécanique ont été discutées. Premièrement, l'approche du couplage local/global non-invasif a été présentée. Cette approche permet des modifications du modèle globale existant sans altérer sa matrice de rigidité. La méthode repose sur un processus itératif. En suite, la méthode LaTIn a été présentée. Elle permet de modéliser les comportements multiples et variés des interfaces non linéaires. Le principe de la méthode LaTIn pour la décomposition du domaine est de séparer les équations non linéaires de l'interface de celles liées aux sous-domaines, ces dernières étant linéaires et indépendantes du sous-domaine. L'avantage est donc que les équations des sous-domaines peuvent être résolues efficacement en parallèle et que les équations non linéaires de l'interface peuvent être écrites localement de manière simple, quelle que soit la nonlinéarité traitée.

Le second chapitre aborde le problème de la modélisation des non-linéarités locales dans un modèle linéaire global. Ce chapitre est basé sur la publication [START_REF] Lapina | A fully non-invasive hybrid IGA/FEM scheme for the analysis of localized non-linear phenomena[END_REF]. La présence des non-linéarités locales impliquent notamment que dans ces zones d'intérêt locales, la solution recherchée n'est pas régulière. Comme mentionné précédemment, l'IGA utilise les fonctions de forme d'une régularité supérieure et, par conséquent, ne semble pas adaptée à la modélisation de ces zones locales. Contrairement à l'IGA, la FEA semble adaptée à la simulation de comportements locaux, fortement non linéaires, voire singuliers, en raison de sa régularité réduite et de la flexibilité géométrique. De plus, de nombreuses implémentations améliorées de 105 la FEA, à la fois efficaces et robustes, existent pour simuler divers comportements locaux.

Dans ce contexte, nous avons porposé une approche de modélisation hybride globale-IGA/locale-FEA, dans laquelle nous avons couplé un modèle global IsoGéométrique (IG) avec un modèle local d'Éléments Finis (FE). Cette approche bénéficie donc des deux technologies d'analyse : l'efficacité de l'IGA pour la description géométrique et la représentation d'une réponse globale et régulière, et la capacité de la FEA à calculer les comportements locaux, fortement non linéaires, voire singuliers. Le couplage entre le modèle globale et le modèle locale est réalisé de manière noninvasive. Cela permet de ne pas incorporer une région locale dans le modèle IG initial globale, ce qui évite naturellement les procédures coûteuses de re-paramétrisation des splines. De plus, l'opérateur de rigidité IG global peut être assemblé et factorisé une seule fois et le système IG à résoudre reste bien conditionné quelle que soit la forme de la région locale. En outre, nous avons utilisé le lien FEA-to-IGA existant, basé sur l'extraction de Bézier ou de Lagrange, pour transformer l'interface au sein du modèle IG global initial en une interface FE sur laquelle le maillage FE local peut être construit. Une interface de couplage conforme a ainsi été construite et nous avons pu exprimer le couplage hybride IGA/FEA en utilisant uniquement l'opérateur FEA-to-IGA et les opérateurs de traçage FE standards qui sont disponibles dans les codes FE classiques. Ainsi, la stratégie proposée peut être utilisée d'une manière prête à l'emploi avec n'importe quel code FE robuste adapté à la modélisation de comportements locaux complexes.

Notre stratégie est non-invasive dans deux sens : le couplage global/local est effectué en suivant l'algorithme non-invasif et tous les opérateurs IGA sont calculés dans un logiciel FEA sans qu'il soit nécessaire de modifier ce logiciel.

Pour valider la stratégie, une série de simulations numériques a été réalisée. La mise en oeuvre a été effectuée à l'aide du progiciel industriel FE à source ouverte code-aster. Tout d'abord, un problème académique a été examiné et les résultats étaient prometteurs en termes de qualité de la solution obtenue et de convergence du maillage. Ensuite, il a été enrichi avec la présence de trous, de fissures et de contacts par friction à l'échelle locale, afin de montrer le potentiel de la méthodologie. Nous nous sommes ensuite concentrés sur les matériaux composites, et une matrice avec plusieurs inclusions a été considérée. Chaque modèle local comprenait une inclusion avec un matériau matriciel adjacent et le délaminage a été considéré à l'interface matrice/inclusion. Cela a permis de mettre en évidence que notre stratégie aboutit à une approche efficace de décomposition de domaine non linéaire. Enfin, nous nous sommes concentrés sur la modélisation d'assemblages en 3D, où le contact par friction a été pris en compte dans le modèle local. Ainsi, même si cette thèse se concentre sur les microstructures composites, la stratégie proposée peut être utilisée avec succès dans d'autres contextes.

Le troisième chapitre, basé sur sur la publication Lapina et al., aborde le problème de la simulation d'une structure linéaire avec des interfaces éventuellement non linéaires. Ce problème peut être vu comme un cas particulière des non-linéarités localisée, consideré dans le second chapitre. Ici, nous avons proposé une approche pragmatique. Tout d'abord, nous avons utilisé l'IGA à la fois pour la représentation géométrique des interfaces et pour la discrétisation du domaine. Ainsi, à l'inverse de la stratégie globale-IGA/locale-FEA non-invasive présentée dans le second chapitre, nous avons utilisé l'IGA à la fois à l'échelle globale et à l'échelle locale. Cela permet de bénéficier pleinement d'une précision supérieure par DDL présentée par l'IGA par rapport à celle présentée la FEA standard. Néanmoins, la modélisation d'objets multi-106 APPENDIX A. RÉSUMÉ DES TRAVAUX (EN FRANÇAIS) échelles et géométriquement complexes, tels que les microstructures composites avec un grand nombre d'inclusions locales, est un défi dans l'IGA en raison de la structure rigide du produit tensoriel des fonctions de forme. Pour répondre à ce problème, nous avons utilisé des méthodes de frontières immergées, ce qui permet au maillage de ne pas suivre la géométrie du modèle. Dans ce cas, le principal défi consiste à formuler et à mettre en oeuvre un couplage entre la matrice et les inclusions multiples. Cette tâche est loin d'être triviale puisque les interfaces de couplage sont censées être à la fois non conformes (l'interface peut couper des éléments arbitraires des mailles) et non linéaires (comme l'incorporation d'un contact ou d'un décollement). Le premier ingrédient de notre approche consiste à prendre en compte la méthode des frontières immergées conformes (Immersed Boundary-Conformal Method , IBCM) récemment introduite, qui s'est avérée précise et robuste pour les interfaces parfaites. Notre approche consiste en une stratégie pragmatique qui transforme l'interface initiale, par la construction de couches conformes, en trois interfaces différentes : l'interface initiale entre la matrice et l'inclusion qui devient conforme, et deux interfaces non conformes qui se trouvent maintenant à l'intérieur de la matrice et de l'inclusion.

Ainsi, dans le présent cadre, l'idée est de séparer les aspects non linéaires et non conformes des interfaces en faisant apparaître deux types d'interface différents: les interfaces non conformes mais parfaites (i.e., qui assurent l'équilibre et la compatibilité cinématique standard entre les sous-domaines) et les interfaces non linéaires mais conformes (i.e., où les limites des sous-domaines sont alignées l'une sur l'autre). Le couplage non conformes est fait par la méthode de Nitsche. La construction d'opérateurs de couplage pour les maillages non conformes a été traitée par une discrétisation spéciale de l'interface : les interfaces ont été discrétisées comme l'intersection des deux maillages de part et d'autre. Cela garantit que sur chaque élément de l'interface , toutes les fonctions de base impliquées des deux sous-domaines à coupler sont des polynômes plutôt que des polynômes par morceaux. Pour modéliser le comportement non linéaire de l'interface, la méthode LaTIn a été mise en oeuvre. Nous l'avons détaillée pour trois cas distincts de comportement d'interface non linéaire : le contact sans frottement, le contact avec frottement et le délaminage. De plus, la méthode LaTIn nous a permit d'obtenir un algorithme de décomposition de domaine parallèle.

Une série de simulations numériques en 2D a été étudiée pour valider la stratégie. Tout d'abord, avant d'appliquer la stratégie complète, nous avons abordé la question de la construction de la couche conforme pour les méthodes des frontière immergées conformes. Ensuite, une matrice avec une inclusion centrale et un contact sans frottement sur l'interface matrice/inclusion a été considérée. Une étude de convergence du maillage a été réalisée pour ce cas. En outre, pour montrer l'extensibilité de notre approche, la convergence en termes d'itérations de la méthode LaTIn a été étudiée pour des modèles avec différents nombres d'inclusions. Enfin, des exemples numériques avec contact frictionnel et délaminage ont été présentés.

  from a FEA perspective . . . . . . . . . . . . . . . . . . . . . 10 1.1.1 Introduction to IGA . . . . . . . . . . . . . . . . . . . . . 10 1.1.2 Basics of NURBS-based IGA . . . . . . . . . . . . . . . . 11 1.1.

  Figure 1.2: Quadratic (p = 2) basis functions: Lagrange polynomials (on the left), B-Splines (on the right) for open, uniform knot vector Ξ = {0, 0, 0, 0.25, 0.5, 0.75, 1, 1, 1}.

  16 CHAPTER 1. STATE OF THE ART: ON THE MULTI-SCALE MODELLING AND COMPUTATION OF COMPLEX STRUCTURES WITH IGA Ku = f .

Figure 1 . 3 :

 13 Figure 1.3: Illustration of the analysis-suitable model issue in IGA. (a) Initial CAD parametrization; (b) Spline re-parametrization to obtain a boundary-fitted analysissuitable model; (c) Discretizations considered in the case of immersed methods. In Fig. (c) the meshes are associated with the basis functions of the mechanical fields, while the grey areas correspond to the integration domains.

  STATE OF THE ART: ON THE MULTI-SCALE MODELLING AND COMPUTATION OF COMPLEX STRUCTURES WITH IGA

  Figure 1.7: Schematic representation of the coupling problem: two domains, Ω 1 and Ω 2 , should be coupled through interface Γ.

  32 CHAPTER 1. STATE OF THE ART: ON THE MULTI-SCALE MODELLING AND COMPUTATION OF COMPLEX STRUCTURES WITH IGA as

  for local plasticity, Fuenzalida-Henriquez et al. (2022), Li et al. (2021), Meray et al. (2022), Passieux et al. (2013) for crack propagation, Aldakheel et al. (2021), Gerasimov et al. (2018), Noii et al. (

Figure 1

 1 Figure 1.8: Example of a global/local problem. The global model over subdomain Ω 12 is replaced by the finer local model of domain Ω 2 through interface Γ, which enables to integrate geometrical details (e.g., a hole) along with possible non-linearities around (e.g., contact, cracks, plasticity) within the initial regular coarse model.

  ) should be verified. Interface conditions Eqs. (1.18) with the 36 CHAPTER 1. STATE OF THE ART: ON THE MULTI-SCALE MODELLING AND COMPUTATION OF COMPLEX STRUCTURES WITH IGA current notation become:

Figure 1

 1 Figure 1.9: The iterative non-invasive exchange procedure. Starting with a global Neumann resolution (with respect to Γ), the interface displacement is transferred from the global to the local model. Then, a local Dirichlet problem is solved and the interface traction force coming from the local model is applied to the global model along with the interface traction force produced by the covered part of the global model at previous iteration.

  , Gendre et al. (2009), Li et al. (2021) to name a few). Examples of non-linear local behaviour are also shown in this work, see Chapter 2.

  Figure1.10: Latin method: a schematic depiction of the iterative process ( S(i) is an approximation before relaxation).

  3.2. However, here, we use an IG model for the global representation while employing FEA to construct the local model. The global IG model has a coarse discretization and represents the entire structure, characterized by the physical domain Ω 1 (see Figure 2.1(a)(left)). As in Section 1.3.2, Ω 1 is divided in Ω 11 and Ω 12 , then a local, more detailed FE model characterized by domain Ω 2 is constructed to replace the global model in Ω 12 (see Figure 2.1(a)(right)). The substitution of the FE local model within the IG global one is achieved through interface Γ. Thus, we seek to solve a hybrid IGA/FEA multi-domain problem in Ω 11 ∪ Ω 2 ∪ Γ, the global solution in Ω 12 being discarded.

Figure 2 . 1 :

 21 Figure 2.1: Example of a global-IGA/local-FEA problem. The global IG model over sub-domain Ω 12 is replaced by the finer local FE model of domain Ω 2 through interface Γ, which enables to integrate geometrical details (holes) along with possible non-linearities around (e.g., contact, cracks, plasticity) within the initial regular coarse model.

  Figure 2.3, which visualizes the communications between the global IG and local FE models in line with Eqs. (2.3)-(2.4).

Figure 2 . 3 :

 23 Figure 2.3: Communications between the global IG and local FE model through the fully non-invasive strategy.

  4b

  denotes the reference exact strain energy and E h the strain energy of the discrete model. The convergence curve is given versus the number of DoF in Figure2.5 (see the green curve). The number of DoF was computed as the sum of the global IG DoF and of the local FE DoF. For comparison purpose, the convergence curves of the equivalent single-model solutions are also plotted: "Standard IGA" represents the solu-

Figure 2 . 4 :

 24 Figure 2.4: Global/local non-invasive analysis of the linear curved beam problem (NURBS mesh composed of quadratic 24 (circumferential direction) × 16 (radial direction) elements for Ω 1 , and standard FE mesh composed of 128 T6 elements for Ω 2 .

  Figure 2.5: Convergence of the relative energy error for the linear elastic 2D circular beam (quadratic case).

  (a) Local T3 triangles (in black) compared to the initial quadratic IG mesh over region Ω 12 (in red). (b) Displacement field. (c) Von Mises stress.

Figure 2

 2 Figure 2.6: Global/local non-invasive analysis of the linear curved beam problem (NURBS mesh composed of quadratic 24 (circumferential direction) × 12 (radial direction) elements for Ω 1 , and standard FE mesh composed of 512 T3 elements for Ω 2 .

  Figure 2.7: Convergence of the relative energy error for the linear elastic 2D circular beam with global quadratic NURBS elements and local T3 triangles.

Figure 2

 2 Figure 2.8: Non-invasive introduction of holes, cracks and frictional contact between the lips of the cracks in an initial 2D NURBS beam with the proposed hybrid global-IGA/local-FEA methodology (the cracks, where a contact model is applied, are colored in blue).

  (a) Von Mises stress. (b) Zoom on contact zones of the deformed configuration (scale factor 1)contact is reached on the red zones. (c) Relative displacement discrepancy with respect to an equivalent full FE solution.

Figure 2 . 9 :

 29 Figure 2.9: Solution obtained for the 2D curved beam with holes, cracks and contact with the proposed hybrid global-IGA/local-FEA algorithm.

Figure 2

 2 Figure 2.10: Convergence of the non-invasive global/local algorithm for the 2D curved beam with holes, cracks and contact.

Figure 2 .

 2 Figure 2.11: Description of the 2D plate problem with multiple inclusions and delamination. 1, 4 or 16 inclusions are considered. The global/local interfaces are underlined in blue while the cohesive interfaces are plotted in green. One local FE model includes one inclusion plus a matrix region surrounding the inclusion, and thus one cohesive interface. The different model parameters are indicated on the right.

Figure 2 .

 2 Figure 2.12: Local FE meshes considered for the plate with multiple inclusions problem. It is made of T 3 triangles. The inclusion is in red and cohesive elements are incorporated at the inclusion-to-matrix interface. On the boundary of the complete local model, i.e. at the global/local interface, 2 T 3 triangles are put in front of one element of the intermediary refined global model, following Remark 2.2.3.

Figure 2 .

 2 Figure 2.13: Bilinear law for the cohesive elements at the inclusion-to-matrix interfaces.

Figure 2

 2 Figure 2.14: Convergence of the non-invasive global/local algorithm for the plate problem with multiple inclusions and delamination. One graph corresponds to one loading step.

  (a) Stress σxx. (b) Stress σyy.

Figure 2 .

 2 Figure 2.15: Stress distribution for the plate problem with 16 inclusions at the last loading step.

Figure 2 .

 2 Figure 2.16: Zoom on a local model for the plate problem with 16 inclusions at the last loading step.

  Figure2.17: Reaction force versus prescribed displacement at each loading step for the plate problem with 16 inclusions.

Figure 2 .

 2 Figure 2.18: 3D mechanical assembly example: description and data of the problem.

  (a) Overview and zoomed window around the global/local interface. (b) Discretization of the contact zone (a cut is performed in between the two perforated plates in contact).

Figure 2 .Figure 2

 22 Figure 2.19: Meshes for the 3D mechanical assembly example.

Figure 2 .

 2 Figure 2.21: Obtained Von Mises stress and deformed configuration (scale factor 200). Half of the structure is shown (a cut is performed in the symmetry plan of the structure).

Figure 2 .

 2 Figure 2.22: Highlighting of the compression zones (i.e. where σ xx is negative).

70CHAPTER 3 .

 3 Figure 3.1: Schematic representation of the reference problem: two domains, Ω m and Ω f , are coupled through interface Γ.

Figure 3 .

 3 Figure 3.2 for the notations): ∀i ∈ {m, f } and ∀j ∈ {l, b}, div(σj i ) + f g i = 0 in Ω j i ; σ j i = C i ε(u j i ) in Ω j i ; u j i = u g

  whether the solution involves sticking or sliding. Now, let us present the equations describing the local stage in each case. Sticking case consists in solving the following equations:

) 78 CHAPTER 3 .

 783 IMMERSED BOUNDARY-CONFORMAL ISOGEOMETRIC LATIN METHOD FOR MULTIPLE NON-LINEAR INTERFACES

80CHAPTER 3 .

 3 IMMERSED BOUNDARY-CONFORMAL ISOGEOMETRIC LATIN METHOD FOR MULTIPLE NON-LINEAR INTERFACES popular cohesive laws, and Abrate et al. (

Figure 3 . 4 :

 34 Figure 3.4: Infinite plate with a circular hole: (a) problem definition; (b) discretization example: the conformal layer is shown in green, the background mesh in black; Γ m is a non-conformal perfect interface, t is the thickness of the conformal layer.

,

  where u h and u ref denote the solution of interest and the reference solution, respectively. Reference solution u ref is the exact solution to the Kirsch problem.

Figure 3 . 6 :

 36 Figure 3.6: Matrix with a central inclusion in frictionless contact: description and data of the problem.

Figure 3 . 7 :

 37 Figure 3.7: Model decomposition: each part (matrix and inclusion) is composed of two subparts, a conformal layer, and a background mesh. The meshes used for the solution approximation in the matrix are colored in black while the meshes used for the solution in the inclusion are colored in red. The grey color corresponds to the integration area for stiffness operators.

86CHAPTER 3 .

 3 Figure 3.8: Horizontal (left) and vertical (right) displacements plotted over the deformed geometry with a scale factor of 10.

Figure 3 .Figure 3 .

 33 Figure 3.9: σ xx stress (on the left) and σ yy stress (on the right). Deformed geometry with a scale factor of 10.

88CHAPTER 3 .

 3 Figure 3.11: Meshes considered for convergence study: matrix model (black) and inclusion model (red).

( a )

 a Energy norm error for splines of second and third degrees computed in the whole domain.(b) Energy norm error for splines of second and third degrees computed only in the matrix background mesh (matrix without its conformal layer).

Figure 3 .

 3 Figure 3.12: Convergence plots for 2-and 3-degree splines.

Figure 3 .

 3 Figure 3.13: Description of the matrix with multiple inclusions problem: 1, 2, 4, and 16 inclusions. Black color corresponds to the matrix mesh, orange -to the inclusions meshes. Model parameters are given at the bottom.

90CHAPTER 3 .

 3 Figure 3.14: σ xx stress (on the left) and σ yy stress (on the right) for the model with 4 inclusions. Deformed geometry with a scale factor of 10.

Figure 3 .Figure 3 .

 33 Figure 3.15: Integration points on the interfaces: in contact (orange) and not in contact (grey). Model with 4 inclusions.

Figure 3 .

 3 Figure 3.18: Reaction force (in N) on interface 3 for variation in friction coefficients µ 1 and µ 2 .

94CHAPTER 3 .

 3 Figure 3.19: Horizontal displacement at the end of the last time step. Deformed geometry with a scale factor of 50. Four pairs of coefficients considered. Note that for the two central cases (pairs (µ 1 = 0, µ 2 = 0.6) and (µ 1 = 0.6, µ 2 = 0)) the displacement varies throughout interface 3 (with the maximum equal to gap j) which means that only a part of interface 3 is in contact.

Figure 3 .

 3 Figure 3.20: Delamination problem description.

Figure 3 .

 3 Figure 3.21: Reaction force versus applied displacement at each loading step.

Figure 3 .

 3 Figure 3.23: σ yy stress in MPa for the last time step. Deformed geometry with a scale factor of 2.

conclusions Appendices A Résumé des travaux (en Français) Bibliography Introduction INTRODUCTION Figure

  Reference problem . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.2 Principle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.3 Non-conformal coupling with Nitsche's approach . . . . . . . . . 3.3 Enforcement of non-linear interface behaviours . . . . . . . . . . . . . . 3.3.1 Implementation of the linear stage . . . . . . . . . . . . . . . . . 3.3.2 Treatment of the non-linearity in specific cases . . . . . . . . . . 3.3.3 Search directions . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Numerical results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4.1 On the thickness of conformal layers in IBCM . . . . . . . . . . . 3.4.2 Matrix with a central inclusion in frictionless contact . . . . . . . 3.4.3 Matrix with multiple inclusions . . . . . . . . . . . . . . . . . . . 3.4.4 Frictional contact problem . . . . . . . . . . . . . . . . . . . . . . 3.4.5 Matrix with one inclusion and delamination . . . . . . . . . . . . 3.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

	2	TABLE OF CONTENTS
	Table of contents 3.2.1 General

State of the art: on the multi-scale modelling and computation of complex structures with IGA

  . Immersed boundaryconformal isogeometric LaTIn method for multiple non-linear interfaces. Submitted for publication. https://hal.science/hal-04224392 2. Lapina, E., Oumaziz, P., Bouclier, R. and Passieux, J.C., 2023. A fully noninvasive hybrid IGA/FEM scheme for the analysis of localized non-linear phenomena. Computational Mechanics, 71(2), pp.213-235.

• Conference presentations 1. Oumaziz, P., Lapina, E., Bouclier, R., Passieux, J.C. "Couplage hybride IGA/FEM non-intrusif pour la simulation globale/locale non-linéaire de structures". 15ème Colloque National en Calcul des Structures, 2022, Giens, France. INTRODUCTION 1.

Isogeometric Analysis from a Finite Element Analysis Perspective 1.1.1 Introduction to IGA

  General idea . . . . . . . . . . . . . . . . . . . . . . . . . . 39 1.4.3 General LaTIn numerical scheme . . . . . . . . . . . . . . 39 1.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

		CHAPTER 1. STATE OF THE ART: ON THE MULTI-SCALE MODELLING
	10	AND COMPUTATION OF COMPLEX STRUCTURES WITH IGA
	1.1

[START_REF] Lapina | IsoGeometric LaTIn method for non-conformal coupling with non-linear interface behaviour[END_REF]

Establishing the Link between IGA and FEA . . . . . . . 17 1.1.4 Advantages and limitations of IGA and FEA . . . . . . . 21 1.2 Immersed Boundary Methods . . . . . . . . . . . . . . . . . . . . . 22 1.2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . 22 1.2.2 Fundamentals and Challenges . . . . . . . . . . . . . . . . 23 1.2.3 Weak coupling strategies . . . . . . . . . . . . . . . . . . . 27 1.3 Global/local coupling . . . . . . . . . . . . . . . . . . . . . . . . . 32 1.3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . 32 1.3.2 Non-Invasive Approach . . . . . . . . . . . . . . . . . . . . 34 1.4 LaTIn method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 1.4.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . 38 1.4.2 IsoGeometric analysis (IGA) was originally introduced in Cottrell et al. (2009), Hughes

  .29) The coupling operator C 1 is a prolongation of C 11 into Ω 1 . λ 12 is introduced to denote the discrete reaction forces at Γ produced by the covered part Ω 12 of the global model. STATE OF THE ART: ON THE MULTI-SCALE MODELLING AND COMPUTATION OF COMPLEX STRUCTURES WITH IGA Once the algorithm converged and the solution obtained, the global solution in Ω 12 is discarded (see Figure1.8(b)).

	It emerges to counterbalance the effect of this covered region since this one is not
	present in the reference coupling problem (see Fig. 1.8(b) again). It reads at iteration
	n -1:	λ	(n-1) 12	= K 12 u

(n-1) 11

f 12 ,

(1.30)

where K 12 and f 12 are the extensions into Ω 1 of the classical stiffness matrix K 12 and load vector f 12 of Ω 12 , respectively.

38

CHAPTER 1.

  Linear elastic 2D curved beam . . . . . . . . . . . . . . . . 52 2.3.2 2D curved beam with holes, cracks and contact . . . . . . 54 2.3.3 2D plate with multiple inclusions and delamination . . . . 57 2.3.4 3D mechanical assembly example with preload and frictional contact . . . . . . . . . . . . . . . . . . . . . . . . . 61 2.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65

	This chapter is based on the publication Lapina et al. (2023):
	Lapina, E., Oumaziz, P., Bouclier, R. et al. A fully non-invasive hy-
	brid IGA/FEA scheme for the analysis of localized non-linear phenomena.
	Computational Mechanics, 71(2):213-235, 2023. https://doi.org/10.1007/
	s00466-022-02234-2

2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 2.2 Non-invasive global-IGA/local-FEA . . . . . . . . . . . . . . . . . 46 2.2.1 The global IGA/local FEA problem . . . . . . . . . . . . . 47 2.2.2 Construction of the FE model to reach a conformal global/local interface . . . . . . . . . . . . . . . . . . . . . . . 48 2.2.3 Fully non-invasive implementation of global-IGA/local-FEA 49 2.3 Numerical results . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 2.3.1

  included in code-aster. Finally, every computed fields (displacement, stress) are expressed in terms of FE quantities (see Eq. (1.14a)) so that the standard code-aster post-processing functionalities are used to visualize the results. In the illustrations, we keep the notations introduced previously in the paper; in particular, domain Ω 1 = Ω 11 ∪ Ω 12 ∪ Γ characterizes the global IG model, and the local FE model of domain Ω 2 is expected to replace the global IG model in sub-domain Ω 12 .

  FULLY NON-INVASIVE HYBRID IGA/FEA SCHEME FOR THE GLOBAL/LOCAL SIMULATION OF STRUCTURES following the non-invasive algorithm and all IGA operators are computed within a FEA software without need for any modifications in this software.

Thus, our strategy is non-invasive in two senses: global/local coupling is performed 66 CHAPTER 2. A

  Principle . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 3.2.3 Non-conformal coupling with Nitsche's approach . . . . . 71 3.3 Enforcement of non-linear interface behaviours . . . . . . . . . . . 72 3.3.1 Implementation of the linear stage . . . . . . . . . . . . . 73 3.3.2 Treatment of the non-linearity in specific cases . . . . . . 74 3.3.3 Search directions . . . . . . . . . . . . . . . . . . . . . . . 81 3.4 Numerical results . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 IMMERSED BOUNDARY-CONFORMAL ISOGEOMETRIC LATIN METHOD FOR MULTIPLE NON-LINEAR INTERFACES This chapter is based on the publication Lapina et al.: Lapina, E., Oumaziz, P., Bouclier, R. et al. Immersed boundary-conformal isogeometric LaTIn method for multiple non-linear interfaces. Submitted for publication. https://hal.science/hal-04224392

	3.4.1 On the thickness of conformal layers in IBCM . . . . . . . 82
	3.4.2 Matrix with a central inclusion in frictionless contact . . . 84
	3.4.3 Matrix with multiple inclusions . . . . . . . . . . . . . . . 89

3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 3.2 Isogeometric immersed LaTIn method . . . . . . . . . . . . . . . . 69 3.2.1 Reference problem . . . . . . . . . . . . . . . . . . . . . . 69 3.2.2 3.4.4 Frictional contact problem . . . . . . . . . . . . . . . . . . 90 3.4.5 Matrix with one inclusion and delamination . . . . . . . . 93 3.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95 68 CHAPTER 3.

  Model description for the frictional contact problem. Black colour corresponds to domain Ω 1 , orange -to domain Ω 2 . Note that due to the geometric simplicity of domain Ω 2 , the latter is directly made conformal to the frictional interfaces, so only one layer is built for Ω 1 at each of these interfaces.

	CHAPTER 3. IMMERSED BOUNDARY-CONFORMAL ISOGEOMETRIC LATIN
	92	METHOD FOR MULTIPLE NON-LINEAR INTERFACES
	Figure 3.17: Parameters	Values
		E	210 GPa
		ν	0.3
		h	50 mm
		e	10 mm
		j	0.04 mm
		d	5 mm
		F max 1	30 MPa
		F max 2	30 MPa
		µ 1 , µ 2	∈ [0, 0.6]
		∆t	1 s
		η	10 -7

Table 3 .

 3 1: Model parameters for the frictional contact problem.

(a) problem description and geometrical parameters.(b) spline discretizations with a zoom on the non-matching interface.

Remerciements Acknowledgements

(b) Energy norm error for splines of second degree computed only in the background mesh. The error was computed separately for the conformal layer and the background mesh. Figure 3.5a shows the convergence plot for the conformal layer with respect to mesh size indicator h. As expected, and conformally to the results revealed in [START_REF] Wei | Immersed boundary-conformal isogeometric method for linear elliptic problems[END_REF], thinner the conformal layer is, smaller its error is. As for the background mesh, the inverse behaviour is observed: the model with the thinnest conformal layer has the biggest error (see Figure 3.5).

These results imply that the solution in the whole domain Ω may not be sensitive to the thickness of the conformal layer. Nevertheless, further research is necessary to draw definitive conclusions in this regard.

Matrix with a central inclusion in frictionless contact

To start with, we considered a single circular inclusion in a rectangular matrix, as described in Figure 3.6. Both matrix and inclusion materials were linear elastic, and we assumed plane stress state and small deformations. However, a large elastic contrast between the matrix and inclusion materials was prescribed: in particular, the Young modulus was chosen 80 times larger for the inclusion than for the matrix. Contact behaviour without friction is considered on the matrix/inclusion interface. Due to the symmetry of the problem, only one quarter of the area was actually modeled with corresponding symmetrical boundary conditions (see again Figure 3.6).

We decomposed the model into two independent subdomains corresponding to the matrix and to the inclusion. Then, following our approach, in each model we added two conformal layers on the matrix/inclusion interface (see Figure 3.7). The thickness of each of the conformal layers was 0.2mm. Thus, each subdomain consists of two subparts: a conformal layer and a regular background mesh. Each layer was coupled by Nitsche method with its corresponding model (through a non-conformal but perfect interface). The stiffness operators of these corresponding regular background meshes, following FCM, were built by integrating solely in the physical domain which is colored in grey on Figure 3.7. Together, the subdomains were connected by the developed LaTIn method through interface Γ (a non-linear but conformal interface). We recall Appendices