Keywords: models, 18 2, 2, 3 Discretization methods

The microstructural modelling of multiple microcracking with competing crack paths and complex crack intersections, branching and interactions in heterogeneous quasibrittle cohesive grain-based materials is a complex task with many far-reaching applications. For instance, transfer properties in grain-based rocks are strongly impacted by microcracking, which requires its modelling. The numerical modelling of such phenomena requires three main ingredients: (i) a proper account for the microstructural heterogeneous geometry representation, i.e. a description of grain boundaries, (ii) a discretization of the governing equations, allowing a kinematical description of crack propagation, and (iii) the formulation of an efficient crack propagation model.

In this thesis a flexible and general stable displacement-Lagrange multiplier mixed formulation is developed to model distributed cracking in cohesive grain-based materials in the framework of the cut finite element method using a non-conforming background mesh. The displacement field is discretized on each grain separately, and the continuity of the displacement and traction fields across the interfaces between grains is enforced by Lagrange multipliers. The design of the discrete Lagrange multiplier space is detailed for bilinear quadrangular elements with the potential presence of multiple interfaces/discontinuities within an element. We give numerical evidence that the designed Lagrange multiplier space is stable and provide examples demonstrating the robustness of the method. Relying on the stable discretization, a cohesive zone formulation equipped with a damage constitutive formulation expressed in terms of the traction is used to model the propagation of multiple cracks at the interfaces between grains. The damage propagation is governed by an energetic formulation. To prevent the crack faces from self-penetrating during unloading, a contact condition is enforced. The solutions for the mechanical fields and the damage field are separately obtained and an explicit damage update algorithm allows using a non-iterative approach. The damage formulation couples the normal and tangential failure modes, accounts for different tension and compression behaviours and takes into account the compression-dependent fracture energy in mixed mode. The framework is applied to complex 2D problems inspired by indirect tension tests and compression tests on heterogeneous rock-like materials.

Tools and state of the art Summary

First, I would like to thank my supervisors, Prof. Thierry J. Massart and Prof. Nicolas Moës for their continuous help during my PhD. I am very grateful to Thierry because of the regular meetings and his prompt replies even at the weekends. He patiently reminded me several times of the physical nature of the problem and the global picture when I went off the track. The additional funding he provided helped me finishing the PhD. I express my appreciation to Nicolas as well. In the first few months, I had to explain my progress rigorously in front of a whiteboard. Although inconvenient at that time, I truly benefited a lot from his teaching technique later. His tendency to use few assumptions when solving a problem and striving for simple algorithms greatly influenced how I solve a task today.

Benoît Lé explained me the original version of the explicit damage update algorithm. I also thank him for keeping answering my questions after I had left ECN. My regular discussions at the ECN with Abdullah Waseem, Baptiste Reyne, Erwan Grelier and Quentin Ayoul-Guilmard about science, mathematics, programming and L A T E X typesetting increased my knowledge and gave me a wonderful time in Nantes. Additional thanks to Quentin for answering my questions in functional analysis. At the ULB, I would highlight Héctor Hernández with whom I had several interesting talks about discretization methods.

The members of the jury made interesting remarks and asked for more details, leading to better explanations in the final version of the thesis.

The everyday chats and daily interactions with the other members of the reasearch groups at the ULB and at the ECN also contributed to my well-being.

My work was supported by the Erasmus Mundus Joint Doctorate SEED project (European Commission, 2013-1436/001-001-EMJD). Their financial help is gratefully acknowledged.

The aim of this chapter is three-fold. A short introduction to the challenges of rock mechanics modelling is given. Then, with the application of the developed computational framework in mind, a series of objectives and major modelling steps are formulated, indicating which sections of the work elaborate on them. Finally, the major contributions of the present work are emphasized.

This introductory chapter starts with the motivation of the thesis, explaining why it is relevant to simulate cracking in heterogeneous grain-based materials. It is well-known that heterogeneous materials exhibit complex behaviours which are difficult to capture at the macroscopic level with phenomenological constitutive equations. Therefore, the modelling is carried out at the meso-scale where the shape, size distribution and material properties of the constituents can be taken into account. These microstructural properties are considered as parameters and the developed framework allows carrying out simulations on different realizations of aggregates of grains. We focus our attention on quasi-static crack propagation in rocks as a practical application. It turns out from the available literature that in many practically relevant cases, it is a good approximation to assume that the cracks propagate along grain boundaries (interfaces).

Restricting the possible crack paths to grain boundaries allows relatively cheap simulations. For representing cracks, models can be based on the incorporation of degradation of those interfaces according to a failure criterion. Alternatively, one can assume that the cracks are possibly present on all the interfaces from the beginning, but they are only activated at a position when a physical criterion predicts so. This way, maintaining closed, a priori positioned would-be cracks leads to constraint equations. Their mathematical treatment is classical in the continuous setting but can cause instability when the governing equations are discretized.

Experiments on rocks show that they can fail under compression if the accompanying shear is large enough. Another experimental observation is that the shearing fracture energy increases with normal compression. These phenomena, that have to be accounted for, make the modelling challenging.

In view of such challenges, the major contributions of this work include the formulation of a stable mixed cut finite element method with non-conforming Cartesian meshes of Q1 elements arbitrarily cut by interfaces, choosing the crack path for multiple competing cracks, and a cohesive model which takes into account the variable fracture energy. The solution scheme allows an embarrassingly parallel implementation.

Background and motivation

Rocks are geological materials with complex behaviour due to many factors [START_REF] Jing | A review of techniques, advances and outstanding issues in numerical modelling for rock mechanics and rock engineering[END_REF]. Some rock types are porous and in certain applications, the thermo-or hydromechanical coupling cannot be neglected. Rocks, in most cases, contain pre-existing fractures both at the macro-and at the micro-scale. Moreover, they are elastically inhomogeneous and anisotropic. In "weak" rocks, creep and relaxation may be relevant mechanisms [START_REF] Jing | A review of techniques, advances and outstanding issues in numerical modelling for rock mechanics and rock engineering[END_REF]. Also, rocks cannot be considered as virgin materials and due to tectonic movements, rock masses are subjected to a pre-existing stress state. Figure 1.1: Microstructure of a Lac du Bonnet granite. Mineral types: quartz (Qtz), plagioclase (Plag), potassium feldspar (Kfsp), biotite (Bt). Microcracking types: grain boundary (Grb), intragranular (Intr), transgranular (Tr). Picture taken from [START_REF] Lim | In-situ stress and microcracking in granite cores with depth[END_REF].

Making decisions in rock engineering problems is made more difficult due to size effects. Commonly used rock mechanics tests typically operate on samples of sizes ranging from 10 mm to 500 mm. By measuring their stress-strain response, macroscopic data (Young's modulus, Poisson's ratio, etc.) can be obtained. However, it is known that the variation in the type and size of the constituents may influence the macroscopic properties. It is important to take the size effect into account when making predictions on rock masses based on experiments. The surface roughness of rock joints depends on the specimen size, leading to size effect as well. As the hydromechanical properties also depend on the surface roughness, it is relevant to characterize it somehow. The authors in [START_REF] Fardin | The scale dependence of rock joint surface roughness[END_REF] carried out experiments on different sample sizes and realized that a minimum size is necessary so that the surface roughness in a rock mass can be characterized by measurements on laboratory samples.

Various failure mechanisms act in rocks. Experiments showed that failure of granite in compression is governed by the following phenomena [START_REF] Martin | The progressive fracture of Lac du Bonnet granite[END_REF]: (i) closure of existing microcracks, (ii) linear elastic deformation, (iii) crack initiation and stable crack growth, (iv) unstable crack growth (v) post-peak behaviour. Microcracking leads to the progressive degradation of the macroscopic properties.

We will consider here cohesive grain-based materials that are made up of tight packing of constituents without significant void parts. To separate these grains, non-zero traction is needed. Working with such cohesive grain-based materials, including rocks, is important in many geomechanical applications. Nuclear waste disposal in geological formations is nowadays considered. Although the type of the neighbouring rock is selected to have low permeability, the excavation required for disposal alters the stress state in the surrounding of the underground opening [START_REF] Van Den Eijnden | Multiscale modelling of the hydromechanical behaviour of argillaceous rocks[END_REF]. The stress redistribution in the excavation disturbed zone initiates microcracking in the rock. The microcracks grow and nucleate to form larger scale fractures which act as preferential pathways for fluid. The infiltration of nuclear waste to the surface through these channels has serious environmental consequences. One of the major sources of material nonlinear behaviour in brittle rocks such as granite is microcracking [START_REF] Souley | Damage-induced permeability changes in granite: a case example at the URL in Canada[END_REF]. It is therefore of significant importance to study the permeability increase with respect to microcracking, which requires robust computational tools representing distributed cracking.

Several experimental techniques exist to measure the material parameters at the micro-scale and at the macro-scale. The ultrasonic pulse velocity testing was used in [START_REF] Vasconcelos | Ultrasonic evaluation of the physical and mechanical properties of granites[END_REF] to determine the tensile and compressive strengths and Young's modulus in granites. Macroscopic laboratory tests such as the direct tension test, uniaxial and triaxial compression tests, or the Brazilian test can provide the global strength and Young's modulus based on the measured stress-strain curve. The individual phases (mineral species) often have different material properties that may be difficult to evaluate. Micro-indentation is a test procedure used to determine the micro-mechanical properties. For heterogeneous materials, grid indentation is better suited [START_REF] Mahabadi | A novel approach for micro-scale characterization and modeling of geomaterials incorporating actual material heterogeneity[END_REF]. The fracture toughness can also be extracted from scratch test results. Similarly to the characterization of the microstructure, numerical methods can also be used to investigate which microscopic material property values allow reproducing an experimentally observed probability distribution.

Newer modelling strategies and the improving computational capabilities made numerical simulations a useful technique, complementing experimental approaches to solve engineering problems in rock mechanics. Due to the lack of available data and the complexity of rock behaviour, numerical simulations can often only provide qualitative results [START_REF] Jing | A review of techniques, advances and outstanding issues in numerical modelling for rock mechanics and rock engineering[END_REF]. They are however important because they can complement experiments that have their own difficulties:

• long sample preparation time • in heterogeneous materials, such as rock, the experimentally observed crack path substantially differs for each sample. It is therefore almost impossible to obtain the local characteristic material properties and the repeatability of the experiments is also an issue.

• the response of rock samples is often brittle. Hence the loading device must have a very careful displacement control or mixed control facility; suitable loading speed must be applied to achieve a quasi-static response and prevent abrupt failure [START_REF] Vasconcelos | Experimental characterization of the tensile behaviour of granites[END_REF] • due to size effect, the properties measured on the sample may not represent the value of parameters in the large scale Numerical methods allow the user developing a qualitative understanding of rock fracturing and allows quickly performing parametric studies and sensitivity analyses [START_REF] Jing | A review of techniques, advances and outstanding issues in numerical modelling for rock mechanics and rock engineering[END_REF]. Depending on a specific engineering problem at hand, only few of the previously mentioned mechanisms are considered in the developed model to manage the complexity. Comparing the numerical and the experimental results is made harder by the fact that real-world materials always contain heterogeneities and existing microcracks. If measurements are available only at the macro-scale, careful calibration of the fine-scale quantities is required in the numerical model. Rocks are heterogeneous both at the macroscopic and at the microscopic scales. Macroscopic heterogeneity comes from the different rock types and from the preexisting fractures. The micro-level heterogeneity is partly caused by pre-existing microcracks and the different mineral species with distinct material properties present in the rock [START_REF] Mahabadi | Investigating the Influence of Micro-scale Heterogeneity and Microstructure on the Failure and Mechanical Behaviour of Geomaterials[END_REF]. Accounting for such heterogeneities is crucial for an accurate description of such geomaterials. The distribution and shape of the grains can be determined by µCT, followed by image processing to map the colours of the CT image to the aggregate types [START_REF] Mahabadi | A novel approach for micro-scale characterization and modeling of geomaterials incorporating actual material heterogeneity[END_REF]. If the exact microstructure of a particular sample is not considered important, computer-generated microstructures can be cheap alternatives. They allow analysing the effect of specific features such as the size distribution, the convexity and the clustering of species. Among the many microstructure generation frameworks, we mention the one in [START_REF] Sonon | A unified level set based methodology for fast generation of complex microstructural multi-phase RVEs[END_REF], which is able to generate a wide variety of, possibly densely packed, microstructures with great control. For cohesive grainbased materials, Voronoi tessellations have long been applied [START_REF] Alonso-Marroquín | Role of anisotropy in the elastoplastic response of a polygonal packing[END_REF][START_REF] Massart | Stress-induced permeability evolution in a quasi-brittle geomaterial[END_REF][START_REF] Van Den Eijnden | A FE 2 modelling approach to hydromechanical coupling in cracking-induced localization problems[END_REF]. Assigning different elastic properties (E, ν) to the Voronoi cells was demonstrated to be able to capture the effect of elastic heterogeneity [START_REF] Lan | Effect of heterogeneity of brittle rock on micromechanical extensile behavior during compression loading[END_REF].

The challenges of distributed microcracking modelling and the many important applications motivated the topic of this thesis.

Objectives

From the above introduction it becomes clear that the proper characterization of the microstructure and the microcracks is crucial to understand the failure of rock-like materials. This thesis is concerned with the numerical modelling of diffuse or multiple microcracking in cohesive grain-based materials, with an emphasis on rocks such as granite. By diffuse or multiple cracking, we mean a configuration in which cracks can propagate along multiple, a priori known, paths with potential related crack interactions, intersections and branching. This is not to be confused with completely arbitrary, a priori unknown, microcracks modelled as damage growth in the bulk. Although there exist numerous numerical techniques, some of them being reviewed in Section 2.1, there is still room for new models and solution methods which are accurate, consisting of physical parameters, and that require little user intervention. It is not our intention to replace the existing modelling methodologies, rather we wish to develop an efficient framework endowed with the following capabilities.

P.

The physical model must capture some relevant mechanisms in the failure process of rocks. It includes P. [START_REF] Abdelaziz | Grain based modelling of rocks using the combined finite-discrete element method[END_REF] an intergranular crack propagation framework. In grain-based rock microstructure the preferential locations for microcrack formation are often the grain boundaries. This mechanism is called intergranular cracking. Cracks going through the grain are referred to as transgranular crack propagation. Several factors have an influence on the nature of the fracture of polycrystalline materials. When studying piezoelectric polycrystals, it was found that intergranular cracks are mainly present for small grains, while transgranular cracks emerge as the average grain size increases [START_REF] Verhoosel | Modelling inter-and transgranular fracture in piezoelectric polycrystals[END_REF]. Intergranular fracture is favourable in polycrystalline materials if the fracture toughness of the interfaces is less than that of the grains, as reported in [START_REF] Nguyen | Multi-phase-field modeling of anisotropic crack propagation for polycrystalline materials[END_REF]. The competition between inter-and transgranular fracture was thoroughly analysed in [START_REF] Shabir | The role of cohesive properties on intergranular crack propagation in brittle polycrystals[END_REF], where the authors identified that the magnitude of the fracture toughness of the grains and the grain boundaries have a significant impact whether intergranular or transgranular fracture is realized in quasi-brittle materials.

P.2 a failure criterion which takes into account the compression for damage initiation. The compression-dependence of the mechanical response is significant in rocks, especially considering heterogeneity and deep configurations. It was observed that rocks do not only fail in tension but also in shear combined with compression. In rock mechanics, two of the most commonly used failure criteria are the Mohr-Coulomb and the Hoek-Brown [START_REF] Hoek | Empirical Strength Criterion for Rock Masses[END_REF] criteria. Although well-suited for the shear and compression regime, these two empirical failure envelopes are not successful in fitting experimental results in tension. Therefore, both criteria were combined with cut-offs in the tension part of the stress state.

P.3 a non-negligible fracture process zone. Initially, two types of approaches were used in rock mechanics: the strength of materials and the linear elastic fracture mechanics (LEFM). The first one predicts the onset of damage initiation when a certain failure criterion is met, but does not model the post-peak response. LEFM can be used to investigate the conditions for the propagation of a crack, but the validity requirements of LEFM do not generally hold in rock masses. Indeed, cracks have a considerable characteristic length compared to the whole sample size. Furthermore, [START_REF] Rubin | Tensile fracture of rock at high confining pressure: Implications for dike propagation[END_REF] noticed that in a deep underground configuration at large confining pressure, the process zone is significant. This is why rocks should be modelled as quasi-brittle and not completely brittle materials.

P.4 non-uniform fracture energy. Material parameters, including fracture properties, of rocks are not uniform in general. The dependence of the fracture energy on the compressive stress is indicated in Fig. 3.9. in [START_REF] Goodman | Introduction to Rock Mechanics[END_REF].

P.5 a continuum-based model. Discrete modelling methods are a natural choice for granular materials such as sand, but computations on initially cohesive rocks can efficiently be carried out by discretization methods that rely on a continuum description. The finite element method is widely used because of its mature mathematical basis, its accuracy, and the accumulated experience with it.

M. The mechanical model should M. [START_REF] Abdelaziz | Grain based modelling of rocks using the combined finite-discrete element method[END_REF] give rise to simple equations, the discretization of which does not lead to a too costly numerical scheme. The grains of the rock are considered to undergo small deformation and rotation, and are therefore described by the equations of linear elasticity. All the complex phenomena are assumed to take place at the grain boundaries. Recalling what we wrote in P.3, this modelling choice begs for the use of the cohesive zone concept.

M.2 make use of the a priori known potential crack paths. Owing to criteria P. [START_REF] Abdelaziz | Grain based modelling of rocks using the combined finite-discrete element method[END_REF] and M.1, all the interfaces are assumed to be completely open. To suppress these openings before the crack actually appears there, constraint equations are needed, which are detailed in Section 2.4. In addition to the displacement field defined on the individual grains, the Lagrange multiplier field is also present, responsible for enforcing the constraints.

M.3 couple well with the cohesive model. As the cohesive tractions of the cohesive zone formulation are related to the Lagrange multipliers (shown in Appendix A), the Lagrange multipliers are not eliminated from the two-field formulation.

D. The discretization method

D.1 must lead to a stable discrete mixed method for reliability and optimal rate of convergence. The background for stability is discussed in Sections 2.4.2.2 and 2.4.3, and is applied to our problem in Section 3.6.

D.2 should require no intervention from the user. Mesh generation for complicated grain configurations can be challenging for an automatic mesh generator. Human help in the meshing is tiring and time-consuming, especially when many grain realizations are tested in statistical studies. To decouple the function approximation from the mesh, the cut finite element method will be used. The idea behind mesh-independent methods is explained in Section 2.5 and will be used for our problem in Sections 3.4-3.5 of Chapter 3. By mesh-independent in this thesis, we point to a method that allows solving the equilibrium problem based on a non-conforming mesh, i.e. the mesh is decoupled from the domain. This is not to be confused with meshless methods, which only contain nodes without a mesh data structure.

D.3 should be scalable. The computational geometric algorithms ideally have linear or close to linear complexity so that many grains can be handled in a reasonable amount of time. Large grain count is necessary to represent the real microstructure. Parallelization of the main steps of the complete framework is also desirable.

Major contributions

To complete the objectives worded above, the following main originalities are presented in this thesis.

• A stable discrete Lagrange multiplier space is constructed for the displacement space discretized with bilinear quadrilateral elements. The inf-sup conditions are verified computationally to demonstrate that the constructed Lagrange multiplier space is stable.

• An unfitted structured background mesh is sufficient to treat complex heterogeneous aggregates of grains with arbitrary polygonal shapes.

• An explicit interfacial damage formulation is developed for crack propagation combining a Mohr-Coulomb criterion with tension cut-off and compression cap. Contact is handled in this formalism and different behaviours for tension and compression are taken into account with normal compressive stress-dependent mixed mode fracture energy.

• The above features and the careful implementation enable to conduct fast computations on highly heterogeneous structures.

• The use of the developed framework is illustrated based on indirect tension test-like configurations (e.g. Brazilian test) and on uniaxial compression tests, as frequently used on rock specimens in the literature.

This chapter provides a literature review on research areas the subsequent chapters heavily rely on.

The first topic is about crack propagation in rocks and is structured around three concepts: the main phenomena underlying it, the existing models to describe it and the numerical methods developed so far. All of these aspects are discussed in the view of the objective to represent multiple, concurrent cracking in cohesive grain-based materials.

11

The continuity of the displacement and traction fields in the primal formulation of the elasticity equations for disjoint bodies can be enforced by constraints. The various techniques to achieve this are discussed in this chapter, and their strengths and weaknesses are compared. Imposing constraint equations at the discrete level for non-matching meshes requires special care and robust algorithms borrowed from computational geometry.

When constraints are enforced with Lagrange multipliers or when mixed finite elements are used, the discretized weak form gives rise to an algebraic saddle point problem. The existence and uniqueness of such linear systems is important when designing stable formulations. The solution of saddle point systems is more intricate than linear systems arising from the single-field finite element methods, partly because they are indefinite and often ill-conditioned.

Partition of unity based discretization methods became very popular over the last twenty years. After giving the rationale of why they are useful in certain situations, we compare some of these techniques, the way they generalize the finite element method and the corresponding emerging challenges. We provide some applications of these discretization methods, mainly concentrating on cracking. Our focus is on the extended finite element method (X-FEM) and the cut finite element method (CutFEM).

Some modelling techniques for rocks

A very detailed overview of the use of computational methods in rock mechanics can be found in [START_REF] Jing | A review of techniques, advances and outstanding issues in numerical modelling for rock mechanics and rock engineering[END_REF]. Here we only highlight microstructure-based and phenomenological models, homogenization and the discrete element method as the main class of models.

As mentioned in Section 1.1, the complex behaviour of rocks can be partly attributed to their microstructure. It is rather difficult to construct phenomenological constitutive laws at the macro-scale accounting for microstructural features. An accurate macroscopic stress-strain relation therefore usually requires constructing empirical expressions which fit the experimental data. However, except for a few cases, experimental information is scarce. Furthermore, repeating experiments for each loading case and material remains very time-consuming. An alternative for modelling is to take into account the heterogeneity of the microstructure, on a scale at which the behaviour of the constituent materials can be characterized by simpler and well-understood constitutive models. When dealing with cohesive grain-based geomaterials, the grains can be assumed to behave in a linear elastic way, and the complex phenomena related to cracking can be assumed to remain restricted to the grain boundaries (interfaces). Some publications use macroscopic variables motivated by phenomena at the micro/meso-scale. In [START_REF] Halm | A Model of Anisotropic Damage by Mesocrack Growth -Unilateral Effect[END_REF], a macroscopic damage tensor -a symmetric second-order tensor -was introduced. Its spectral decomposition shows that three mutually orthogonal meso-cracks can be described by such a tensor. Their model also handled meso-crack closure under compressive loading, resulting in the recovery of the elastic moduli. This macroscopic damage model was extended in [START_REF] Halm | An anisotropic model of damage and frictional sliding for brittle materials[END_REF] by considering frictional sliding on the surfaces of the meso-cracks as another dissipative mechanism. This improvement allowed the partial recovery of the shear moduli as well. Frictional sliding has a blocking effect on crack growth and is an important ingredient in the microstructural modelling of rocks. The authors of [START_REF] Halm | A modular damage model for quasibrittle solids -interaction between initial and induced anisotropy[END_REF] further modified their previous model, keeping the macroscopic description. In that work, damage-induced anisotropy was complemented with initial anisotropy by incorporating additional tensors (so-called fabric tensors). In some rocks, irreversible strains can be important. In [START_REF] Shao | Modelling of induced anisotropic damage in granites[END_REF], the authors added a term to the free energy function of [START_REF] Halm | A Model of Anisotropic Damage by Mesocrack Growth -Unilateral Effect[END_REF] to account for such irreversible strains without introducing plasticity. They tested their model on the degradation of a Lac du Bonnet granite around a tunnel.

Modelling the whole rock mass at the micro-level is, of course, impossible with the current numerical methods and computers. This is why homogenization techniques became popular. Among the several homogenization methods, computational homogenization is the most successful for complex (evolving) microstructures, accompanied by possibly large rotations (e.g. fragmentation of rocks). It was first used by [START_REF] Smit | Prediction of the mechanical behavior of nonlinear heterogeneous systems by multi-level finite element modeling[END_REF] and its main features are described in [START_REF] Geers | Multi-scale computational homogenization: Trends and challenges[END_REF].

Homogenization is a powerful approach, allowing the incorporation of inter-and transgranular cracking and multiphysics at the micro-scale. The effect of microcracking on the permeability increase in granite was analysed in [START_REF] Massart | Stress-induced permeability evolution in a quasi-brittle geomaterial[END_REF] using first-order computational homogenization for the mechanical problem. The RVE consisted of an assembly of grains generated by Voronoi tessellation and intrinsic cohesive elements were inserted between each pair of grains. An exponential cohesive law was parametrized with the interfacial damage, from which the permeability was computed. Loading the RVE causes the microcracks to open, acting as fluid channels. The phenomenological coupling of the permeability and the microcracking in [START_REF] Massart | Stress-induced permeability evolution in a quasi-brittle geomaterial[END_REF] was replaced by the parallel plate flow model in [START_REF] Massart | Computational modelling of crackinduced permeability evolution in granite with dilatant cracks[END_REF], relying on the computed crack opening. The PhD thesis [START_REF] Van Den Eijnden | Multiscale modelling of the hydromechanical behaviour of argillaceous rocks[END_REF] and the related papers [START_REF] Van Den Eijnden | A FE 2 modelling approach to hydromechanical coupling in cracking-induced localization problems[END_REF][START_REF] Van Den Eijnden | Modeling the strain localization around an underground gallery with a hydro-mechanical double scale model; effect of anisotropy[END_REF] were dedicated to the hydro-mechanical coupling at the micro-scale, combined with a macroscopic second-order gradient model, to determine the localization in claystone.

Besides continuum models, rocks are often described by discrete particles at the microstructural scale. These methods are preferred to continuum-based approaches in case of fragmentation. The discrete element method (DEM) models granular materials as an assembly of rigid or deformable parts, tracking the motion of the constituents and their contact. It is well-suited for naturally fragmented materials or when the complete understanding of the model is not known. On the other hand, DEM simulations usually require many parameters to be calibrated based on available experimental test data [START_REF] Potyondy | A bonded-particle model for rock[END_REF]. The DEM usually contains more parameters than a continuum-based model. Those parameters can be classified as material properties and numerical parameters [START_REF] Coetzee | Review: Calibration of the discrete element method[END_REF]. The material properties at the scale of the particles are either calibrated through macroscopic tests (Brazilian test, direct shear test, compression tests) or directly measured (which is experimentally complicated). The material properties in the DEM are similar to the that of the continuum models at the scale of the grains (e.g. bond strength, Young's modulus, bond stiffness). On the other hand, the DEM requires several parameters that are based on purely numerical grounds, such as the particle shape (often taken spherical to ease the contact detection), damping parameters, resistance to particle rotation, etc. Among these parameters, the particle shape has a major effect on the tensile to compressive strength ratio [START_REF] Cho | A clumped particle model for rock[END_REF]. The problem with spherical particles is alleviated by the clumped particle model in [START_REF] Cho | A clumped particle model for rock[END_REF], grouping neighbouring particles, handled as a single rigid body, to achieve irregular shapes. To ease the process of calibration, calibration charts were introduced in [START_REF] Fakhimi | Application of Dimensional Analysis in Calibration of a Discrete Element Model for Rock Deformation and Fracture[END_REF] to form dimensionless quantities. The combination of the FEM and DEM, the so-called FEM-DEM, was introduced in [START_REF] Munjiza | A combined finite-discrete element method in transient dynamics of fracturing solids[END_REF] in which a finite element mesh is associated with each discrete element which allows their deformation. An intrinsic cohesive zone with a large initial stiffness is used to detect the onset and propagation of damage. When the cohesive zone fails (complete separation), the discrete elements are handled as individual bodies. Their positions are tracked by explicit time integration and various contact algorithms ensure the non-penetration of the fragments. For loading with high strain-rate, FEM-DEM is especially advantageous. However, FEM-DEM models generally contain more parameters than the pure FEM method, resulting in longer calibration time. Moreover, FEM-DEM simulations are computationally more demanding.

Modelling the degradation of quasi-brittle materials

The fracture process zone is the region ahead of the crack tip where nonlinearities take place. Quasi-brittle materials are characterized by a relatively large fracture process zone. Materials such as concrete, rocks, ice, masonry and some ceramics are considered quasi-brittle. Their degradation is characterized by strain localization at the continuum level. During this intensive deformation, caused by the initiation and coalescence of microcracks, the gradient of the (equivalent) strain becomes large in a narrow band. As the load increases, this band becomes narrower and narrower, and finally, a macroscopic crack emerges. The macroscopic crack can be described by the discontinuity of the displacement field across the crack surface. This representation is referred to as discrete crack model. An energetically equivalent characterisation is the smeared crack approach in which the crack is not represented as a surface of discontinuity, but rather as a damaged region of finite width. Both approaches have their advantages as detailed below.

Discrete crack models

Discrete crack approaches idealize the strain localization band as a sharp interface.

In continuum mechanics, the displacement solution to the boundary value problem (BVP) describing an elastic body under loading is the minimizer of the potential energy Π = U -F , where U is the strain energy and F is the work done by external forces. Griffith's theory of fracture [START_REF] Griffith | The Phenomena of Rupture and Flow in Solids[END_REF] laid the foundation of fracture mechanics.

The amount of energy required to open a crack of surface area A is called the surface energy and is computed as S = G c A, where the material constant G c is the fracture energy. On the other hand, the formation of a new crack dissipates a portion of the so-far accumulated elastic energy of the body. This is characterized by the energy release rate giving the decrease of the potential energy upon the formation of crack faces:

G = - ∂Π ∂A (2.1)
If the released energy can cover the surface energy, the appearance of the crack is favourable. Therefore, the energy balance according to Griffith can be written, reformulated into a variational setting [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF], as

A Ȧ ≥ 0 the crack can only grow (irreversibility) (2.2) G ≤ G c
the energy release rate is bounded

(2.3) (G -G c )A Ȧ = 0 crack propagates when the energy release rate is critical (2.4)
In Griffith's model the crack path is assumed to be known beforehand. Then, the computed crack increment determines the position of the crack tip. In many real-world situations, the crack path is not known a priori. In those cases, additional techniques have to be used to determine the direction of the crack growth (see [START_REF] Bouchard | Numerical modelling of crack propagation: automatic remeshing and comparison of different criteria[END_REF] for a comparison of some of those techniques in conjunction with the finite element method), which is particularly difficult in three-dimensional crack propagation. Another drawback of the Griffith model is that only existing cracks can be propagated, it cannot initiate a crack.

Linear elastic fracture mechanics (LEFM) builds on Griffith's model, assuming that the body is linear elastic everywhere except along the crack. In LEFM the displacement and stress fields have singularities at the crack tip. The LEFM is valid for brittle fracture in which the inelastic deformations are confined to the neighbourhood of the crack tip.

The variational approach to brittle fracture was developed to compensate for the deficiencies of the Griffith model. After reformulating Griffith's original theory with equations (2.2)-(2.4), the authors in [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF] developed a variational framework which contains the crack path as an unknown. The quasi-static displacement solution and the energetically optimal crack path are obtained at the same time by finding the global minimum of the total energy Π + S. In [START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF] an elliptic regularization was proposed. A detailed description and the summary of its numerical implementation can be found in the monograph [START_REF] Bourdin | The Variational Approach to Fracture[END_REF].

In the previous models, it was assumed that the size of the fracture process zone is negligible compared to the size of the body. In ductile and quasi-brittle materials this assumption is not valid. The cohesive zone model (CZM) was introduced by Dugdale [START_REF] Dugdale | Yielding of steel sheets containing slits[END_REF] for ductile and by Barenblatt [16] for quasi-brittle materials. According to this model, the nonlinearities (coalescence of microcracks, crack bridging, etc.) in the fracture process zone are lumped to a single surface and are expressed in the form of a traction-separation law. In other words, the CZM replaces the stress-strain constitutive model of the localization band with a traction-separation relation. Let us presume that the crack path for the whole fracturing process is known and is denoted by Γ. Then the CZM distinguishes three parts of Γ as shown in Fig. 2.1. In the fully developed crack region the tractions are zero. In the cohesive zone the tractions change as a function of the displacement jump. Further away from the crack the material is intact. The length of the cohesive zone is inversely proportional to the brittleness of the material and is defined as [START_REF] Hillerborg | Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements[END_REF]

ℓ ch = EG c f 2 t , (2.5)
where f t is the tensile strength. The mathematical crack tip separates the cohesive zone from the intact part of the body, while the physical crack tip marks the point which still has zero traction. Originally written for mode I crack opening, the CZM was later modified to work under local mixed mode loading. It is customary to define equivalent, sometimes called effective, quantities in the cohesive zone: the equivalent traction t eq (t n , t t ) and the equivalent displacement jump u eq ( u n , u t ), where the subscripts denote the normal and tangential components of the traction vector t and the displacement jump vector u . It is important to note that the tractions remain finite as opposed to LEFM in which the stress is unbounded at the crack tip.

The CZM is a simple and powerful concept, being able to represent the continuous degradation of the material through the softening law (marked in red in Fig. 2.1), and can also predict the initiation of the crack due to the presence of the strength in the model. Moreover, the fracture energy of the Griffith model is also present in the CZM:

G c = u eq,c ∫ 0 t eq d u eq (2.6)
Another salient feature of the CZM is that it is unnecessary to track the position of the crack tip as the traction-separation relation automatically gives the mathematical crack tip when the strength t eq,c is reached. The shape of the traction-separation curve ideally comes from experiments. It does not necessarily have to be a monotonically decreasing function but its maximum value in a homogeneous body must be the strength, as explained in [START_REF] Elices | The cohesive zone model: advantages, limitations and challenges[END_REF]. In the work of Dugdale, the traction-separation function was constant. A linear relation was first used in [START_REF] Hillerborg | Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements[END_REF] and later polynomial and exponential functions were constructed. For some quasi-brittle materials, the shape of the softening function has a significant effect [START_REF] Volokh | Comparison between cohesive zone models[END_REF][START_REF] Alfano | On the influence of the shape of the interface law on the application of cohesive-zone models[END_REF]. For granite-type rocks, the strength and the fracture energy are the two main parameters. The cohesive zone models can be derived from a potential function, making them thermodynamically consistent in case of mixed mode loading [START_REF] Park | Cohesive Zone Models: A Critical Review of Traction-Separation Relationships Across Fracture Surfaces[END_REF].

Cohesive zone models were initially formulated with an initial elastic region [START_REF] Xu | Void nucleation by inclusion debonding in a crystal matrix[END_REF], as illustrated in Fig. 2.2a with linear softening. In contrast to these intrinsic cohesive models, the extrinsic CZM of [START_REF] Camacho | Computational modelling of impact damage in brittle materials[END_REF] is initially rigid with linear softening (see Fig. 2.2b). It allows the adaptive insertion of the cohesive zones in the bulk when a certain material failure criterion is met. The need for this additional failure criterion, not contained in the cohesive zone formulation, led to the term extrinsic. The cohesive zone is initiated when the equivalent traction t eq at a point reaches the critical value t eq,c . This critical value is the strength and is given by some failure criterion. Common in the models reviewed so-far is that they represent the crack as a sharp discontinuity. This allows handling contact and friction between the crack lips, modelling fluid transport in open cracks and avoiding distorted meshes in the FEM implementation. On the other hand, discrete crack models are valid only after extreme localization. Therefore, they cannot describe precisely the microcrack formation. Furthermore, their treatment in FEM needs special care as we will see in Section 2.2.3.

Smeared crack models

Continuum damage mechanics introduces a progressive degradation of continuum properties of the material and is able to initiate and coalesce diffuse microcracks, thereby naturally handling crack branching and crack interactions. The stress-strain relation of isotropic damage in a linear elastic body is given by

σ = (1 -D)C : ε, (2.7)
where D ∈ [0, 1] is the bulk damage. Standard continuum models using local stressstrain relations are inadequate for softening because the width of the localization band tends to zero upon mesh refinement, and therefore the dissipated energy tends to zero as well. This way, failure can happen without energy dissipation, which is clearly not physical. The reason for this behaviour of the model is that the governing differential equations lose their ellipticity, indicated by the singularity of the acoustic tensor [START_REF] Jirásek | Mathematical analysis of strain localization[END_REF].

To extend the validity of continuum damage models for highly localized strain, they are regularized by introducing an internal length scale, which is a parameter of the model and sets the size of the damaged zone [START_REF] Lasry | Localization limiters in transient problems[END_REF]. Since this length is not zero, the dissipated energy remains finite. In non-local constitutive models the stress at a point depends not only on the strain at that point but also on the strain in the neighbourhood [START_REF] Pijaudier-Cabot | Nonlocal Damage Theory[END_REF]. Strain gradient models introduce higher-order gradients of the strain field [START_REF] Peerlings | Gradient enhanced damage for quasi-brittle materials[END_REF]. These higher-order strains appear with coefficients involving the length scale. As they include higher-order derivatives, additional boundary conditions need to be imposed. Others keep the Cauchy continuum description and introduce the gradients of internal variables, such as the gradient of the damage. Introduced in [START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF], the phase-field model is an elliptic regularization of the Francfort-Marigo functional of [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF] to tackle its two main difficulties: strong discontinuity and unknown crack path, both difficult to implement in the conventional FEM. The solution proposed in [START_REF] Bourdin | Numerical experiments in revisited brittle fracture[END_REF] consists of including a function -later called phase-field -varying between 0 (undamaged) and 1 (fractured), defined on the whole domain. An additional parameter ℓ 0 , the length scale, was also introduced. It was proved that this modified functional tends to the original one in [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF] as ℓ 0 → 0.

The crack band model was put forward in [START_REF] Bažant | Crack band theory for fracture of concrete[END_REF] as a simple alternative to the previously mentioned regularization techniques. In this model, the stress-strain relation is not fixed for a given material, but it is scaled according to the width of the localization band. The crack band model is equivalent to the CZM in the sense that the same strength is used and the fracture energy is also set to be the same by changing the constitutive law according to the width of the band. Its implementation in FEM is easy because only the material stiffness needs to be reduced. Since the correspondence of the fracture energy is derived for uniaxial tension, and the stress state in the crack front is not uniaxial, not the same amount of energy is dissipated for different band widths [START_REF] Jirásek | Numerical aspects of the crack band approach[END_REF]. For skewed finite element meshes, an equivalent crack band width must be determined, and simple criteria based on element area do not give proper results [START_REF] Jirásek | Numerical aspects of the crack band approach[END_REF]. Therefore, the crack band model is an easy-to-use but inaccurate method to circumvent the pathological mesh-sensitivity of softening materials in the FEM.

Given that continuum damage models are better at predicting the initial phase of degradation and discrete crack models have advantages for fully developed cracks, some attempts were made to implement a transition from continuum damage to sharp discontinuity in an automated way. Reference [START_REF] Cazes | A thermodynamic method for the construction of a cohesive law from a nonlocal damage model[END_REF] describes the transition from a non-local damage model to CZM. The thick level set approach [START_REF] Moës | A level set based model for damage growth: The thick level set approach[END_REF] is not only a regularization technique but it also automatically inserts sharp discontinuities in the location where the material is completely damaged. We refer to [START_REF] Tamayo-Mas | A continuous-discontinuous model for crack branching[END_REF] for an up-to-date review on the existing continuous-discontinuous models.

Discretization methods

The finite element discretization of continuum damage models requires sufficiently many elements in the width of the band of localization to be able to capture the large strain gradient. If the band is narrow, i.e. the characteristic length scale is small, a fine mesh is needed. If it is initially not known where the localization takes place, a very fine overall mesh is necessary, leading to a significant computational cost. Moreover, occasional remeshing is required to avoid the overly distorted elements. In phase-field models, the additional difficulty lies in finding the global optimum of the phase-field functional which is numerically very challenging. From now on, we concentrate on the discretization of discrete crack models.

The FEM working with C 0 shape functions is not capable of representing strong discontinuities within the element. The classical solution is to fit the mesh in such a way that the crack lies on element edges. Assigning distinct degrees of freedom to the nodes of the elements neighbouring the crack discontinuity allows the crack to open. Every time the crack propagates, the mesh needs to be modified in the vicinity of the new crack increment. Even though better and better meshing tools are available, the solution from the previous mesh needs to be transferred to the current mesh. Projecting functions between unrelated (i.e. unstructured) meshes is costly and introduces errors. This motivated researchers to design finite element methods which can capture the discontinuity within the element. The basic assumption of the enhanced assumed strain (EAS) formulation [START_REF] Simo | A class of mixed assumed strain methods and the method of incompatible modes[END_REF] is that the strain tensor can be decomposed into the sum of a compatible part, ∇ s u, and an incompatible part. The incompatible part is again written as the sum of a bounded and an unbounded term. Substituting these decomposed strain tensors into the Hu-Washizu weak form and choosing the space of enhanced strain fields to be orthogonal to the space of stress fields, the stress unknowns and then the incompatible strain unknowns can be eliminated at the element level. The advantage of EAS is that the extra degrees of freedom (discretized enhanced strain field) are defined at the element level. However, the EAS only provides a piecewise continuous displacement jump (i.e. the crack path is not continuous) across the element boundaries. In [START_REF] Roubin | Multi-scale failure of heterogeneous materials: A double kinematics enhancement for Embedded Finite Element Method[END_REF] the EAS formulation was extended to handle both strong (emerging microcracks) and weak (matrix-inclusion boundary) discontinuities. A two-scale approach was used to simulate crack propagation in concrete, where the aggregates were assumed to be linear elastic bodies and the cohesive model was an extrinsic exponential CZM. The paper demonstrated that such simple constitutive relations at the meso-scale are able to represent the experimentally observed anisotropy at the macro-scale. The extended finite element method (X-FEM) is a newer method, which is not an element-based but a node-based approach. Therefore, the displacement jump can be interpolated continuously across the element edges/faces. Moreover, the crack tip enrichment further improves the accuracy of X-FEM. Also, the element type is not restricted as opposed to the EAS method. On the other hand, the computational complexity is higher than that of the EAS because the enriched degrees of freedom are coupled. The X-FEM will be presented in details in Section 2.5.3, therefore only its applications to cracking is discussed in this section.

Implementing the CZM into a FEM framework was first done in [START_REF] Hillerborg | Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements[END_REF]. In the conventional FEM, cohesive zones are zero-thickness interface elements. If the crack path is known (e.g. debonding of composites), these interface elements can be inserted along the crack path. In the intrinsic CZM, the initial elastic part of the constitutive model permits to represent some elastic deformation. In case of an unknown crack path, these cohesive elements are inserted between each pair of bulk elements. To suppress the early activation of the intrinsic cohesive zones, a high penalty (also called dummy) stiffness enforces the closure of the cohesive zone. This method is easy to implement but a reasonably fine mesh is required to capture the true crack path. Moreover, the crack path in the computation is restricted to element edges/faces, which renders the result mesh-dependent. An equally unpleasant consequence of the large penalty stiffness is the large condition number of the system matrix. Furthermore, one has to choose a proper quadrature scheme for interface elements. The use of Gauss-Legendre quadrature rules results in oscillating traction profiles. As the traction appears in the cohesive relation (2.6), the computed opening becomes inaccurate. The first explanation of this behaviour was published in [START_REF] Schellekens | On the numerical integration of interface elements[END_REF] for interface elements. It turned out that the exact Gauss-Legendre integration couples unrelated degrees of freedom of the interface element. If the quadrature points are placed at the nodes of the interface element, the spurious coupling is not present because the shape functions are zero at the unrelated quadrature points. This is why low-order Gauss-Lobatto or Newton-Cotes quadrature formulae solve the oscillation problem. On the other hand, if the dummy stiffness approach is used in X-FEM, the coupling of the unrelated degrees of freedom cannot be eliminated with this technique [START_REF] Simone | Partition of unity-based discontinuous elements for interface phenomena: computational issues[END_REF]. An alternative to the a priori insertion of interface elements is the application of an extrinsic CZM for which the complications mentioned before do not exist. The adaptive insertion of extrinsic cohesive zones requires flexible data structures to rearrange the mesh connectivity every time the crack is extended. Based on the extrinsic CZM of [START_REF] Camacho | Computational modelling of impact damage in brittle materials[END_REF], the authors in [START_REF] Pandolfi | An Efficient Adaptive Procedure for Three-Dimensional Fragmentation Simulations[END_REF] described a reference implementation in C. The X-FEM was also used with cohesive zones, first in [START_REF] Wells | A new method for modelling cohesive cracks using finite elements[END_REF], who inserted extrinsic cohesive zones on-the-fly. A similar approach in [START_REF] Mergheim | A finite element method for the computational modelling of cohesive cracks[END_REF] used the phantom node method to represent the discontinuity across the crack faces. In Ref. [START_REF] Moës | Extended finite element method for cohesive crack growth[END_REF] instead of the four singular crack tip functions used in LEFM, a single non-singular tip enrichment function was introduced to account for the finite displacement at the cohesive crack tip. The crack length control could tackle the snap-back of the load-displacement response and the use of the stress intensity factors improved the accuracy of the crack growth condition compared to the local stress evaluation. The CZM for 3D problems with contact and unloading detection was considered in [START_REF] Ferté | 3D crack propagation with cohesive elements in the extended finite element method[END_REF]. A three-field formulation inspired by [START_REF] Wall | An embedded Dirichlet formulation for 3D continua[END_REF] was the starting point in [START_REF] Sabada | A stable X-FEM in cohesive transition from closed to open crack[END_REF] to create a stable X-FEM approximation for a priori closed cracks.

A challenge in the discrete crack methods is how to locate the crack tip. In the CZM the progressive evolution of the interfacial traction automatically separates the traction-free (fully open) crack from the partially open cracks (fracture process zone). However, in case of brittle fracture, the crack tip must additionally be located. In [START_REF] Sukumar | Extended finite element method for three-dimensional crack modelling[END_REF] the same crack tip enrichments were used for 3D planar cracks as what is conventionally utilized for 2D in LEFM. An alternative in 3D is to use the Heaviside enrichment only but with h -p refinement [START_REF] Duarte | A high-order generalized FEM for through-the-thickness branched cracks[END_REF]. If the crack tip is not represented, it must lie on the element boundary. Therefore, the nodes of the edge containing the crack tip are not enriched with the Heaviside function [START_REF] Wells | A new method for modelling cohesive cracks using finite elements[END_REF]. Similarly, the crack tip element is not doubled in CutFEM [START_REF] Song | A method for dynamic crack and shear band propagation with phantom nodes[END_REF].

The other difficulty in discrete crack models is crack branching. The authors of [START_REF] Daux | Arbitrary branched and intersecting cracks with the extended finite element method[END_REF] introduced special enrichment functions to enable crack branching, which applied the Heaviside enrichment for each crack segment with additional enrichment functions at the junctions. In [START_REF] Budyn | A method for multiple crack growth in brittle materials without remeshing[END_REF] the growth of multiple cracks was modelled using the theory of linear elastic fracture mechanics (LEFM), requiring the handling of multiple crack faces and crack tips, and a crack growth condition. For an elaborate overview of X-FEM applied to cracking problems, [START_REF] Sukumar | Extended finite element method in computational fracture mechanics: a retrospective examination[END_REF] is suggested for reading.

Finally, we mention the element deletion technique, an old and very simple method, which does not need to represent the strong discontinuity, therefore it is applicable in the conventional FEM without any remeshing. In the element deletion method the crack is represented by a set of deleted elements. Technically, the elements are not removed, but their stress and stiffness are set to zero. The element deletion technique allows the standard FEM to model fracture with minimal modification of the code. Although the energy dissipation on an element can be scaled so that it agrees with the surface energy needed to open the crack segment, the method is inherently meshdependent and produces unphysical crack paths especially for structured meshes. The element deletion, the inter-element crack and the X-FEM were compared on several benchmark problems for dynamic crack propagation in [START_REF] Song | A comparative study on finite element methods for dynamic fracture[END_REF].

Intergranular cracking of polygonal grains has been investigated in many works using different numerical techniques. Interpolation defined on polygons is useful when a domain made up with polycrystalline grains need to be discretized. The Voronoi cell finite element method (VCFEM) [START_REF] Ghosh | Voronoi cell finite elements[END_REF], which treats each Voronoi cell as a finite element, was coupled with a three-level homogenization procedure in [START_REF] Ghosh | A multi-level computational model for multi-scale damage analysis in composite and porous materials[END_REF] to simulate damage propagation in polycrystalline and porous materials. The quite recent virtual element method (VEM) found an application in cracking of cohesive grain-based materials in [START_REF] Benedetto | Applications of the Virtual Element Method for Cracking Analysis of Cement-Based Composites Using Interface Elements[END_REF]. Brittle fracture in realistic polycrystalline microstructures was examined using the X-FEM in [START_REF] Sukumar | Brittle fracture in polycrystalline microstructures with the extended finite element method[END_REF] and in the related paper [START_REF] Sukumar | Finite element-based model for crack propagation in polycrystalline materials[END_REF]. The competition between intergranular and transgranular fracture was determined by the ratio of the grain internal and grain boundary fracture toughness values. The authors in [START_REF] Simone | A Generalized Finite Element Method for polycrystals with discontinuous grain boundaries[END_REF] introduced linearly independent generalized Heaviside enrichment functions along with a cohesive zone formulation to allow intergranular crack opening. This framework was utilized in [START_REF] Shabir | The role of cohesive properties on intergranular crack propagation in brittle polycrystals[END_REF] to thoroughly study the effect of the cohesive strength, fracture energy and the resulting cohesive length on the crack path and the load-displacement response. The boundary element method (BEM) with a linear cohesive law on the grain boundaries was used in [START_REF] Sfantos | A boundary cohesive grain element formulation for modelling intergranular microfracture in polycrystalline brittle materials[END_REF] to study intergranular cracking. Only the grain boundaries need to be discretized in this formulation. Inter-and transgranular cracking in 2D and 3D was analysed by BEM in [START_REF] Gulizzi | Modelling intergranular and transgranular micro-cracking in polycrystalline materials[END_REF].

Variational formulation and discretization

To understand the subsequent techniques, let us consider the following abstract variational problem

{ Find u ∈ U such that a(u, v) = ℓ(v), ∀v ∈ V (2.8)
where U and V are vector spaces endowed with norms ∥ • ∥ U and ∥ • ∥ V , respectively. Furthermore, a(•, •) is a continuous bilinear form on U × V and ℓ(•) is a continuous linear form on V . By ℓ(v) the duality pairing ⟨ℓ, v⟩ V ′ ,V is meant, where V ′ is the topological dual space of V , i.e. the space of continuous linear functionals on V . The Lax-Milgram lemma, and under more general conditions, the Banach-Nečas-Babuška theorem ensures the well-posedness of (2.8) [START_REF] Ern | Theory and Practice of Finite Elements[END_REF].

In general, the variational problem (2.8) can only be solved numerically. The Galerkin method is a discretization framework which replaces U and V with the finite-dimensional vector spaces U h (trial space) and V h (test space). Then the discrete version of (2.8) becomes:

{ Find u h ∈ U h such that a h (u h , v h ) = ℓ h (v h ), ∀v h ∈ V h (2.9)
A discretization method is said to be consistent if the exact solution u ∈ U satisfies the discrete problem, i.e. a h (u,

v h ) = ℓ h (v h ), ∀v h ∈ V h . The method is conforming if U h ⊂ U and V h ⊂ V . The bilinear form a(•, •) is called coercive if ∃α > 0 such that a(u, u) ≥ α∥u∥ 2 U , ∀u ∈ U .
The conditions for the well-posedness of the discrete problem (2.9) and the error estimates depend on the type of the discretization.

In a conforming and consistent discretization, the well-posedness of the continuous problem (2.8) directly implies the well-posedness of its discrete version (2.9) and Céa's lemma gives the best approximation property. Additional symmetry and coercivity assumptions on a(•, •) further sharpen the error bound. When the discretization is conforming but inconsistent, the first Strang lemma is applicable, while in case of an inconsistent and non-conforming case, the second Strang lemma can be used. For the detailed mathematical background on these topics, see e.g. [START_REF] Ern | Theory and Practice of Finite Elements[END_REF]. A special case of the Galerkin method is the Bubnov-Galerkin method, which uses the same space as trial and test function spaces: V h = U h . The finite dimensional space V h is given by its basis {ψ i } dim V h i=1 , with dim denoting the dimension of a vector space. The finite element method (FEM) is one of the most common discretization methods to systematically create the approximation functions ψ i with the help of a mesh.

The discretization is demonstrated on the Poisson equation with inhomogeneous Dirichlet boundary condition. Although this is the simplest elliptic PDE, it has all the characteristics we will need later in the thesis. Moreover, this boundary value problem (BVP) has a well-known treatment with weak constraints, making the methods of Section 2.4.2 comparable. We note that the linear elasticity equation, used in this thesis, could have been considered as well, but that would not add more information about the core of the techniques and would make the presentation obfuscated. The model problem is the following.

-∆u = f on Ω u = g on Γ (2.10)
where Γ is the boundary of the domain Ω. We define the spaces V = H 1 (Ω),

V g := {v ∈ V | v| Γ = g} and V 0 := {v ∈ V | v| Γ = 0}
, where v| Γ = g is the shorthand notation for γv = g with γ : H 1 (Ω) → H 1/2 (Γ) being the trace operator. The source term f is assumed to be f ∈ L 2 (Ω). Then the weak form of (2.10) reads:

⎧ ⎪ ⎨ ⎪ ⎩ Find u ∈ V g such that a(u, v) = ℓ(v), ∀v ∈ V 0 where a(u, v) = (∇u, ∇v) Ω , ℓ(v) = (f, v) Ω (2.11)
with (•, •) X denoting the L 2 inner product on X. For a discussion on the spaces H 1 and H 1/2 the reader is referred to [START_REF] Evans | Partial Differential Equations[END_REF]. The Bubnov-Galerkin discretization of (2.11) is:

⎧ ⎪ ⎨ ⎪ ⎩ Find u h ∈ V g,h such that a h (u h , v h ) = ℓ h (v h ), ∀v h ∈ V 0,h where a h (u h , v h ) = (∇u h , ∇v h ) Ω , ℓ h (v h ) = (f, v h ) Ω (2.12)
By comparing (2.12) with (2.11), it is clear that assuming a h = a, ℓ h = ℓ, and integrating by parts shows that the consistency error is zero:

a h (u, v h ) -ℓ(v h ) = 0, ∀v h ∈ V 0,h .
We used an abuse of notation here: by u, the (unknown) representative (e.g. an interpolant or projection) u ȗ ∈ V h of u ∈ V is meant. Therefore, the Lax-Milgram lemma gives the well-posedness of the discrete problem. A boundary condition is said to be strongly enforced if it is treated as an essential boundary condition, i.e. imposed to the approximation space, as done so in this case.

Constraints

As stated in the introduction, all the interfaces are assumed to be completely open from the beginning in our model. Therefore, constraints are needed to suppress these openings (displacement jumps) before the crack actually appears there. In this section we review some techniques used for weakly enforcing constraint equations. It is common in the different discretization methods which use non-matching meshes (fictitious domain methods, immersed boundary methods, X-FEM, CutFEM, etc.) that special care needs to be taken to obtain a stable solution when constraints are enforced weakly. After giving the rationale of dealing with constraints in computational mechanics, this section provides an overview of the most commonly used constraint enforcement techniques and their use in unfitted discretizations. As most of these methods demand stabilization to achieve a non-oscillatory solution, we compare the main stabilization approaches.

Applications leading to constraints

We focus on constraints in continuum solid mechanics. Such constraints may be met in various contexts. Rubber-like materials are incompressible, a phenomenon existing in solid mechanics that requires the use of constraints [START_REF] Bathe | Finite Element Procedures[END_REF]. This example is characteristic of the bulk behaviour. Conversely, jump conditions relate the states of a quantity on two sides of an interface. Material inhomogeneity in solids causes jump in the strain (also called weak discontinuity). On the other hand, a strong discontinuity is associated with a jump in the displacement field, i.e. with the opening of cracks. When two bodies are in contact, the displacement jump at the contact point is zero, a constraint that should be enforced to avoid interpenetration. As a reflection of such physical situations in which constraints are required, constraints are naturally present in computational methods too. It is often convenient both from the mesh generation and from the parallel solution viewpoints to generate unrelated meshes on subdomains. The continuity across these subdomains that would be physically required, has to be maintained by constraints. However, enforcing such boundary or interface conditions in the function space is not possible in meshless and unfitted finite element methods. In those cases, the conditions are therefore weakly prescribed as constraints.

Methods to weakly enforce constraint equations

Non-conforming finite elements, e.g. the Crouzeix-Raviart element [START_REF] Crouzeix | Conforming and nonconforming finite element methods for solving the stationary Stokes equations I[END_REF] have long been used. An alternative weak formulation of (2.11) is

⎧ ⎪ ⎨ ⎪ ⎩ Find u ∈ V such that a(u, v) = ℓ(v), ∀v ∈ V, u| Γ = g where a(u, v) = (∇u, ∇v) Ω , ℓ(v) = (f, v) Ω (2.13)
In the discretization, one may choose to select a h ̸ = a or ℓ h ̸ = ℓ. Such a choice is made for instance when the integration is performed with a non-exact quadrature rule or when stabilizing terms are added to the weak form. In some cases, it renders the method inconsistent as in the penalty method. However, a wise choice of a h and ℓ h can restore the consistency as done in the Nitsche method. The lack of conformity or consistency is sometimes referred to as a "variational crime" [START_REF] Strang | Variational crimes in the finite element method[END_REF]. For further information on such issues, the reader is referred to Chapter 10 in [START_REF] Brenner | The Mathematical Theory of Finite Element Methods[END_REF].

A particular discretization method needs to construct a basis for V h . In the FEM, the goal is to choose the mesh such that it fits to the boundaries and the internal material interfaces. It that case, due to the Kronecker-delta property of the shape functions, the constraints (including the boundary conditions) can be strongly prescribed by simply modifying the nodal values. This was done in (2.12), which is the standard H1 0 -conforming discretization of the Poisson equation. On the other hand, in many practical cases, it is inconvenient to work with subspaces of H 1 0 because the Dirichlet boundary conditions need to be satisfied a priori. In FEM, if the mesh is not fitted to the geometry (called unfitted, non-matching or non-conforming 1 mesh), nodes no longer reside on the material interfaces or the boundaries. The comparison of a conforming and a non-conforming mesh is illustrated in Fig. 2.3. The high quality conforming mesh in Fig. 2.3a was obtained with the Matlab code distmesh [START_REF] Persson | A Simple Mesh Generator in MATLAB[END_REF]. To the right, the same domain is embedded to a simple non-matching background mesh.

One way to build in the constraints to the function approximation is (local) remeshing so as to create a fitted mesh. Another possibility is the modification of the basis functions so that collocation is made possible. However, neither of these methods are convenient, especially for evolving interfaces. The remaining choice is therefore to weakly enforce the constraints by modifying the variational formulation instead of imposing them to the finite element space. Although it provides great flexibility, the well-posedness of the discrete problem in this case cannot be derived from the well-posedness of the continuous problem. In what follows, we show how the three most common methods -penalty, Lagrange multiplier and Nitsche methods -handle the constraints. The example is the Poisson problem analysed before, but other constraints could have been regarded as well, e.g. the jump condition.

Penalty method

The penalty method was introduced in link with the FEM in [START_REF] Babuška | The Finite Element Method with Penalty[END_REF] and later used in non-matching discretization methods, e.g. in X-FEM [START_REF] Chessa | The extended finite element method (XFEM) for solidification problems[END_REF]. The penalty discretized formulation of (2.13) is the following

⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ Find u h ∈ V h ⊂ V such that a h (u h , v h ) = ℓ(v h ), ∀v h ∈ V h where a h (u h , v h ) = (∇u h , ∇v h ) Ω + τ ⟨u h , v h ⟩ Γ ℓ(v h ) = (f, v h ) Ω + τ ⟨g, v h ⟩ Γ (2.14)
Although simple to implement, one problem with the penalty method is that it is inconsistent. It is easy to show that the consistency error is

a h (u, v h ) -ℓ(v h ) = ⟨∇ n u, v h ⟩ Γ , (2.15)
which does not vanish for arbitrary v h ∈ V h . The other disadvantage of the penalty method is that a sufficiently large penalty parameter τ is needed to prescribe the constraint accurately enough. However, for large penalty parameters, the interfacial flux oscillates and the linear system arising from the discretization becomes illconditioned. It was shown in [START_REF] Ferté | 3D crack propagation with cohesive elements in the extended finite element method[END_REF] that the penalty approach to enforce zero opening in the cohesive zone is not a robust method.

Lagrange multiplier method

As the penalty method, the Lagrange multiplier method was also first applied in the FEM by Babuška in [START_REF] Babuška | The finite element method with Lagrangian multipliers[END_REF]. The Lagrange multiplier method is an instance of the abstract mixed formulation

⎧ ⎪ ⎨ ⎪ ⎩ Find (u, λ) ∈ V × Λ such that a(u, v) + b(λ, v) = ℓ(v), ∀v ∈ V b(u, µ) = m(µ), ∀µ ∈ Λ (2.16)
where a(•, •) and b(•, •) are continuous bilinear forms [START_REF] Boffi | Mixed Finite Element Methods and Applications[END_REF]. The formulation (2.16) is also called a saddle point problem which comes from the observation that if a(•, •) is symmetric, the solution to problem (2.16) is the saddle point of the functional

L(u, λ) = 1 2 a(u, u) + b(u, λ) -ℓ(u) -m(λ) (2.17)
Applying (2.16) to our model problem (2.10) yields

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ Find (u, λ) ∈ H 1 (Ω) × H -1/2 (Γ) such that a(u, v) + b(λ, v) = ℓ(v), ∀v ∈ H 1 (Ω) b(u, µ) = m(µ), ∀µ ∈ H -1/2 (Γ) where a(u, v) = (∇u, ∇v) Ω , b(u, µ) = ⟨u, µ⟩ Γ ℓ(v) = (f, v) Ω , m(µ) = ⟨g, µ⟩ Γ (2.18)
where the Lagrange multiplier λ was introduced to enforce the Dirichlet boundary condition as a constraint. Integration by parts allows us identifying the Lagrange multiplier as the normal derivative of the unknown field u:

λ = -∇ n u| Γ . (2.19)
The existence and uniqueness of the solution to (2.16) is guaranteed for any (f, g) ∈ V ′ × Λ ′ if the following two conditions are fulfilled [START_REF] Brezzi | On the existence, uniqueness and approximation of saddle-point problems arising from Lagrangian multipliers[END_REF]:

1. coercivity of a(•, •) in the kernel of b(•, •): ∃α > 0 such that a(v, v) ≥ α∥v∥ 2 V , ∀v ∈ ker b(v, λ) := {v ∈ V | b(v, λ) = 0, ∀λ ∈ Λ} (2.20) 2. inf-sup condition on b(•, •) (LBB condition): ∃β > 0 such that inf λ∈Λ\{0} sup u∈V \{0} b(λ, v) ∥v∥ V ∥λ∥ Λ ≥ β (2.21)
Moreover, the a priori stability estimate

∥u∥ V + ∥λ∥ Λ ≤ c(∥f ∥ V ′ + ∥g∥ Λ ′ ) (2.22)
holds. The individual fields u and λ are bounded by the data, given as

∥u∥ V ≤ 1 α ∥f ∥ V ′ + 2∥a∥ αβ ∥g∥ Λ ′ , (2.23a
)

∥λ∥ Λ ≤ 2∥a∥ αβ ∥f ∥ V ′ + 2∥a∥ 2 αβ 2 ∥g∥ Λ ′ , (2.23b)
which reflects the importance of the stability constants α and β. If these constants are not bounded away from zero, the errors in the data are significantly magnified and the a priori estimates (2.23a)-(2.23b) can only guarantee a poor approximation. The conforming discretization of (2.16) is

⎧ ⎪ ⎨ ⎪ ⎩ Find (u h , λ h ) ∈ V h × Λ h such that a(u h , v h ) + b(λ h , v h ) = ℓ(v h ), ∀v h ∈ V h b(u h , µ h ) = m(µ h ), ∀µ h ∈ Λ h (2.24)
with V h ⊂ V and Λ h ⊂ Λ.

Remark 1.

A distinction is made between hybrid and mixed finite elements [START_REF] Oden | An Introduction to the Mathematical Theory of Finite Elements[END_REF].

A hybrid finite element method uses independent approximations for the unknown within the elements and for its trace on the boundary. A mixed finite element method uses independent approximations for the unknown and for its derivatives within the elements. Since more widespread in the literature, we will use the term mixed for the rest of the work.

Due to conformity and the unchanged linear and bilinear forms (compare (2.16) with (2.24)), the discrete version of (2.18) is formally obtained by the writing the h subscript. To have a stable discrete method, the constants α h , β h > 0 must be bounded irrespective of the mesh size h. In practice, satisfying the discrete version of (2.20) is usually simple. On the other hand, proving that the discrete inf-sup condition (2.21) holds is often very difficult because it involves the interplay of two spaces, V h and Λ h . In many cases, the most convenient choices of these spaces result in an unstable approximation because lim h→0 β h = 0 upon mesh refinement. The inf-sup test [START_REF] Chapelle | The inf-sup test[END_REF], a numerical test, is a necessary condition for (2.21) to be satisfied. This numerical validation proved to be robust in practice. In particular, the inf-sup test has been used successfully to assess the stability of formulations for incompressible materials [START_REF] Chapelle | The inf-sup test[END_REF], a contact algorithm with the Lagrange multiplier method [START_REF] Bathe | The inf-sup condition and its evaluation for mixed finite element methods[END_REF], a shell element [START_REF] Bathe | An inf-sup test for shell finite elements[END_REF], a thermo-hydro-mechanics problem [START_REF] Na | Computational thermo-hydro-mechanics for multiphase freezing and thawing porous media in the finite deformation range[END_REF], and a four-field formulation of fluid flow in porous media [START_REF] Haga | On the causes of pressure oscillations in low-permeable and low-compressible porous media[END_REF], among others. We will present it in details for our problem in Section 3.6.2. Well-known stable approximation schemes of saddle point problems include the Taylor-Hood elements for the Stokes equation that use one degree higher piecewise continuous polynomials for the velocity field than for the pressure field. Stabilization techniques, discussed in Section 2.4.3, enable selecting more versatile spaces V h , Λ h .

Two methods which use Lagrange multipliers to enforce the inter-subdomain continuity of fields are the mortar method [START_REF] Belgacem | The Mortar finite element method with Lagrange multipliers[END_REF] and the FETI method [START_REF] Farhat | A method of finite element tearing and interconnecting and its parallel solution algorithm[END_REF]. The main difference between them is how the Lagrange multiplier field is approximated [START_REF] Lacour | Two different approaches for matching nonconforming grids: The Mortar Element method and the Feti Method[END_REF]. Traditionally, the discrete Lagrange multiplier space in the mortar finite element method consists of piecewise linear functions (P1), modified to constants at junctions (where more than two subdomains meet). The corresponding mathematical theory is mature for planar interfaces, while [START_REF] Puso | Mesh tying on curved interfaces in 3D[END_REF] demonstrated numerically how to do mortaring on curved interfaces. The authors came to the conclusion that the piecewise constant (P0) approximation is stable for closed interfaces (coinciding endpoints) only, whereas P1 is stable with the constant modification at the endpoints.

The penalty method can be derived from the Lagrange multiplier method by substituting λ = τ 2 (u -g) into (2.17) and taking the directional derivative. One advantage of the Lagrange multiplier method is that general nonlinear constraints can be prescribed. The Lagrange multiplier in interface problems often has the physical meaning of the interfacial flux, which is usually of physical interest. Obtaining it as part of the solution of a mixed method gives higher accuracy compared to if it was computed by post-processing the primary field. Although finding a stable pair (V h , Λ h ) is problem-dependent and challenging, once it is done it provides a high-performance element. Disadvantages of the Lagrange multiplier method include the already discussed problem of selecting a stable discretization scheme, the increased number of unknowns and that the matrix is indefinite which precludes the use of some fast iterative solvers and makes the preconditioning more complex.

Nitsche method

The Nitsche method, 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], is a variationally consistent modification of the penalty method. Its idea is to remove the inconsistency (2.15) by adding the so-called consistency term

⟨∇ n u, v h ⟩ Γ to a h . The symmetry term ⟨u -g, ∇ n v h ⟩ Γ
is also added in the original version of Nitsche's method. As u| Γ = g, it does not introduce further inconsistency. The penalty term (often referred to as stability term) is kept to provide the coercivity of a h , and is scaled with the reciprocal of the mesh size. So the discretized expression of Nitsche's formulation is

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ Find u h ∈ V h ⊂ V such that a h (u h , v h ) = ℓ(v h ), ∀v h ∈ V h where a h (u h , v h ) = (∇u h , ∇v h ) Ω -⟨∇ n u h , v h ⟩ Γ -⟨u h , ∇ n v h ⟩ Γ + τ h ⟨u h , v h ⟩ Γ ℓ(v h ) = (f, v h ) Ω -⟨g, ∇ n v h ⟩ Γ + τ h ⟨g, v h ⟩ Γ
(2.25) There exist also unsymmetric and weighted versions of the original Nitsche formulation.

Nitsche's method bears similarities with other discretizations. It can be obtained from the Barbosa-Hughes stabilization as shown in [START_REF] Stenberg | On some techniques for approximating boundary conditions in the finite element method[END_REF] for matching meshes and in [START_REF] Juntunen | On the connection between the stabilized Lagrange multiplier and Nitsche's methods[END_REF] for coupling non-matching meshes. It can also be derived from the Lagrange multiplier method with the substitution λ = -∇ n u + τ 2h (u -g) in (2.17).

Stabilizing constrained formulations

In the FEM, special elements were designed specifically in order to satisfy the discrete version of the inf-sup condition for certain PDEs. Here, we restrict our attention to mixed methods arising from the introduction of constraints with Lagrange multipliers, especially along interfaces when using non-matching meshes, as will be used in the rest of the thesis. A straightforward way to construct the discrete Lagrange multiplier space would be to create a mesh on the interface, consisting of nodes obtained by the intersection of the element edges with the interface, see Fig. 2.4. This is a convenient choice because the quadrature and the assembly can be readily done for each cut element. However, such a choice results in oscillations in the Lagrange multiplier field, as first demonstrated in [START_REF] Ji | On strategies for enforcing interfacial constraints and evaluating jump conditions with the extended finite element method[END_REF]. The reason of this instability is that the dimension of the Lagrange multiplier space is too large compared to the dimension of the displacement space. By properly selecting a subset of the nodes on the interfacial mesh, a stable approximation can be obtained [START_REF] Moës | Imposing Dirichlet boundary conditions in the extended finite element method[END_REF]. This method therefore decreases the dimension of Λ h with respect to V h . An improvement of that algorithm, based on a bulk mesh, was published in [START_REF] Béchet | A stable Lagrange multiplier space for stiff interface conditions within the extended finite element method[END_REF].

Γ ψ 1 2 3 4 Γ ψ 1 ψ 2 ψ 3 ψ 4 1 2 3 4
Figure 2.4: Interfacial mesh for the discretization of the Lagrange multiplier field. This naive approach is not stable [START_REF] Ji | On strategies for enforcing interfacial constraints and evaluating jump conditions with the extended finite element method[END_REF].

Another way to avoid a too rich approximation space for Λ h compared to V h is to keep the space Λ h intact and increase dim V h . This strategy is followed in methods referred to as bubble stabilization. Bubble stabilization of classical finite elements dates back to [START_REF] Arnold | A stable finite element for the Stokes equations[END_REF], in which the famous MINI element was invented for the Stokes equation. The idea was to add a so-called bubble function, a function taking zero value on the boundary of the element, to the finite element approximation such that u = u FEM + u bubble (2.26) This approximation can be interpreted as a subgrid scale model, where the bubble function u bubble represents the unresolved scale (subscale) [START_REF] Hughes | Multiscale phenomena: Green's functions, the Dirichletto-Neumann formulation, subgrid scale models, bubbles and the origins of stabilized methods[END_REF]. The bubble functions are often polynomials and are chosen such that the finite element discretization becomes stable. The variational formulation leads to two unknown fields, the standard finite element field u FEM and the bubble-enrichment u bubble . A single-field formulation can be obtained by the static condensation of the bubbles. In [START_REF] Mourad | A bubble-stabilized finite element method for Dirichlet constraints on embedded interfaces[END_REF], the authors investigated the Laplace equation with interfacial constraints treated by the Lagrange multiplier method on an unfitted mesh, and with stabilization ensured via bubble functions. The method was extended to the linear elasticity problem in [START_REF] Sanders | On methods for stabilizing constraints over enriched interfaces in elasticity[END_REF]. The stability parameter of the model is computed element-wise from the element-level bubble function. Hence, as soon as the bubble function is chosen, the stability can be assured by computing the stability parameter at the element level. The residual-free bubbles, developed in [START_REF] Dolbow | Residual-free bubbles for embedded Dirichlet problems[END_REF] on non-matching meshes, follows a different approach. Instead of defining the bubble function a priori, it is computed by assuming the coarse field u FEM and the Lagrange multiplier field given. Common in the Nitsche and the bubble stabilization methods is that the Lagrange multiplier, if necessary, is obtained through post-processing. Similarly to lowering the number of discrete Lagrange multipliers discussed above, the bubble stabilized method is also free of user-defined parameters.

The Lagrange multiplier method with carefully chosen discrete spaces and the bubble stabilization methods are a priori stable: no tunable parameter remains in the weak form. On the other hand, many stabilized methods contain one or more stabilization parameters. Stabilized methods are often constructed from the pure Lagrange multiplier formulation (2.24) [START_REF] Burman | Deriving Robust Unfitted Finite Element Methods from Augmented Lagrangian Formulations[END_REF]. In the model problem, the Lagrange multiplier was identified as the normal derivative (see (2.19)). Their difference at the discrete level can be used as a penalization, as done in the Barbosa-Hughes stabilization, leading to a consistent mixed formulation:

⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ Find (u h , λ h ) ∈ V h × Λ h such that ∀(v h , µ h ) ∈ V h × Λ h (∇u h , ∇v h ) Ω + ⟨λ h , v h ⟩ Γ + ⟨u h , µ h ⟩ Γ -γ ∑ e∈E h e ⟨λ + ∇ n u, µ + ∇ n v⟩ Γ = = (f, v) Ω + ⟨g, µ⟩ Γ (2.27)
where

V h ⊂ H 1 (Ω), Λ h ⊂ H -1/2 (Γ)
and h e is the element size on the mesh E on Γ. Yet another modification of the Lagrange multiplier method is the augmented Lagrange multiplier method in which the Lagrangian (2.17) is extended with the penalty term

1/2∥ √ τ 0 (u -g)∥ 2 (cf. (2.14)): ⎧ ⎪ ⎨ ⎪ ⎩ Find (u h , λ h ) ∈ V h × Λ h such that ∀(v h , µ h ) ∈ V h × Λ h (∇u h , ∇v h ) Ω + ⟨λ h , v h ⟩ Γ + ⟨u h , µ h ⟩ Γ + τ 0 ⟨u h , v h ⟩ Γ = = (f, v h ) Ω + ⟨g, µ h ⟩ Γ + τ 0 ⟨g, v h ⟩ Γ , (2.28)
where the extra terms compared to the pure Lagrange multiplier method (2.18) are the terms with coefficient τ 0 . Setting

λ h = -∇ n u h , µ h = -∇ n v h and τ 0 = τ /h
gives the Nitsche method (2.25). In recent versions of CutFEM, the coercivity of the bilinear form is maintained by the ghost penalty method [START_REF] Burman | Ghost penalty[END_REF] and the continuity is weakly enforced by Nitsche's method [START_REF] Nitsche | Über ein Variationsprinzip zur Lösung von Dirichlet-Problemen bei Verwendung von Teilräumen, die keinen Randbedingungen unterworfen sind[END_REF], although a stabilized Lagrange multiplier method was also proposed in [START_REF] Burman | Fictitious domain finite element methods using cut elements: I. A stabilized Lagrange multiplier method[END_REF]. When applied to unfitted meshes, this parameter strongly depends not only on the size of the cut element but also on the position of the cut [START_REF] De Prenter | A note on the stability parameter in Nitsche's method for unfitted boundary value problems[END_REF]. To achieve a robust method for a large range of values of τ , the ghost penalty method was introduced in [START_REF] Burman | Ghost penalty[END_REF], which consistently penalizes the jump of the normal derivative along some selected edges near the boundary or the interface. Thanks to it, the condition number scales with mesh refinement in the same rate as in the standard fitted FEM. As a summary, Tab.

H 1 , H -1/2 yes yes [15] bubble (8) in [60] H 1 (Ω Ω) † yes no 2 [60]
1 Courant's method was not yet called finite element method, but it was the first time when continuous piecewise linear functions, defined on a triangular mesh, were used. 2 after the auxiliary fields (bubble, Lagrange multiplier) have been eliminated at the element level † embedded domain formulation with Ω Ω ⊃ Ω Table 2.1: Comparison of some discretization methods for the BVP (2.10)

Enriched finite element methods

There exist numerous closely related discretization techniques relying on similar principles. To understand these methods, some of these techniques are described, for each of them giving the basic idea, the challenges and some applications (mainly related to cracks due to the objectives of the thesis). An emphasis is put on the connection among these methods.

Partition of unity method (PUM)

If the solution of a boundary value problem admits singularities (jump discontinuities, branch points, etc.), higher-order approximations do not necessarily have higher order of convergence. The non-smoothness of the solution can be caused by, among others, the rapidly changing coefficients in the PDE, non-convex domains and the boundary data. One solution is a local refinement around the singularity, which can be challenging for mesh-based discretization methods. In the extreme case, the convergence can be arbitrary slow even upon refinement [START_REF] Babuška | Can a finite element method perform arbitrarily badly?[END_REF]. Another approach to retain the optimal convergence order of an approximation is to include the type of the singularity to the function approximation. In general -not only for singularities -if information about the behaviour of the solution is known, it can be used to improve the approximation properties of discretization methods. Although the use of special functions to improve the accuracy in the presence of singularities was already considered in [START_REF] Fix | On the use of singular functions with finite element approximations[END_REF] in 1973, the use of global basis functions caused the loss of sparsity. To the knowledge of the author, the first contribution which incorporated special functions with local support into the approximation was [START_REF] Babuška | Special Finite Element Methods for a Class of Second Order Elliptic Problems with Rough Coefficients[END_REF] and was coined as special finite element method. The special function in that paper was the product of the finite element hat functions with non-polynomial functions. It was realized in [START_REF] Melenk | The partition of unity finite element method: Basic theory and applications[END_REF] and in [START_REF] Babuška | The partition of unity method[END_REF] that instead of the hat functions, any set of functions satisfying the partition of unity (PU) property suffices. Hence the method was given the name partition of unity method (PUM). These two papers laid the mathematical foundation of PUM. In the PUM, the domain Ω is covered by a family of sets {ω i } N i=1 on which the local approximation spaces V i (i = 1, . . . , N ) can capture the local behaviour of the solution [START_REF] Babuška | Generalized finite element methods -main ideas, results and perspective[END_REF]. The local spaces are then pasted together with the PU functions {ψ i } N i=1 such that the good local approximability is preserved in the global approximation. The precise definition of the partition of unity is due to [START_REF] Babuška | The partition of unity method[END_REF][START_REF] Melenk | The partition of unity finite element method: Basic theory and applications[END_REF]:

Definition 1. Let Ω ⊂ R n be an open set, {ω i } be an open cover of Ω satisfying the pointwise overlap condition ∃M ∈ N, ∀x ∈ Ω card {i | x ∈ ω i } ≤ M.
Let {ψ i } be a Lipschitz partition of unity subordinate to the cover {ω i } satisfying

supp ψ i ⊂ ω i , ∀i (2.29) ∑ i ψ i ≡ 1 on Ω (2.30) ∥ψ i ∥ L ∞ (R n ) ≤ C ∞ (2.31) ∥∇ψ i ∥ L ∞ (R n ) ≤ C G diam ω i , (2.32)
where C ∞ and C G are two constants. Then

{ψ i } is called an (M, C ∞ , C G ) partition of unity subordinate to the cover {ω i }. The partition of unity {ψ i } is said to be of degree m if {ψ i } ⊂ C m (R n ).
The covering sets {ω i } are called patches.

The PU spaces are constructed as

V PU = N ∑ i=1 ψ i V i = span{ψ i φ i j }, (2.33)
where φ i j constitute a basis for V i . In other words, the PU functions ψ i provide the locality, whereas φ i j ensure the local approximating capabilities. Common choices for PU functions include the finite element hat functions, Shepard functions and certain particle shape functions used in meshless methods. The error estimates for the PUM are problem-dependent, but the main achievement of the method can be concluded as follows: the global error is bounded by the local/patch-wise errors. It allows constructing finite element methods in which the convergence is not influenced by the regularity of the solution.

The PUM has the ability to construct high regularity finite element spaces, being an alternative to conforming elements and mixed formulations. The PUM can be used to include so-called handbook functions in the approximation. Handbook functions are solutions to handbook problems -simplified versions of the complex problem at hand, for which the analytical solution is known. In the meshless method community, the concept of PU is represented by the method of finite spheres [START_REF] De | The method of finite spheres[END_REF] and the h-p cloud method [START_REF] Duarte | An h-p adaptive method using clouds[END_REF].

Generalized finite element method (GFEM)

In the context of the finite element method (FEM)2 , the number of patches in the PUM equals the number of nodes and each patch consists of elements incident to a node. This combination of PUM and FEM is known under the name of generalized finite element method (GFEM) [START_REF] Strouboulis | The design and analysis of the Generalized Finite Element Method[END_REF] and the piecewise linear shape functions take the role of the PU basis functions. It is easy to show that the FEM with C 0 linear shape functions is a special case of GFEM (cf. (2.33)):

V FEM = span{ψ i } N i=1 = span ( {ψ i } N i=1 × {1} ) = V GFEM . (2.34)
In 1D, a GFEM with a suitably chosen local basis can be constructed such that it is equivalent to FEM based on continuous piecewise polynomial shape functions of degree k [START_REF] Babuška | Generalized finite element methods -main ideas, results and perspective[END_REF]. The GFEM allows the standard FEM to become largely mesh-independent and improves the approximation properties.

Although the local approximation and the "pasting" of the local approximants are separated in the PUM-based methods, one has to take care of the linear independence of the product shape functions ψ i φ i j . There are two main approaches to solve this. The first one uses specialized solvers which are not susceptible to linear dependency. In [START_REF] Strouboulis | The design and analysis of the Generalized Finite Element Method[END_REF] a specialized direct solver was selected that can handle singular matrices.

In the same paper, as an alternative, the singular system matrix was perturbed, and based on that an iterative solution method was proposed. The second class of approaches ensures that the PUM results in a linearly independent spanning set.

The GFEM was used in [START_REF] Sun | A new family of high regularity elements[END_REF] to create an H 2 -conforming approximation for the biharmonic equation. The GFEM allows constructing high-order approximations without the complexity of the implementation of p-FEM [START_REF] Duarte | Generalized finite element methods for three-dimensional structural mechanics problems[END_REF] (in GFEM, all the degrees of freedom are associated with the nodes). If the solution rapidly changes only in one direction, as in case of boundary layers, it is possible to create higher-order approximation only along that coordinate direction. This directional p-enrichment was used in [START_REF] Duarte | Mesh-independent p-orthotropic enrichment using the generalized finite element method[END_REF], allowing working with a regular mesh, therefore easing the mesh generation process. In some problems, analytical enrichments are not known a priori (e.g. the asymptotic displacement field at a crack tip for linear elastic fracture mechanics in 3D). In such cases, numerical enrichments can be constructed by solving local problems. The local problem is a BVP for which the boundary conditions come from the solution of the global problem. Once the local solution is known, it can be included in the global approximation with the partition of unity [START_REF] Duarte | Analysis and applications of a generalized finite element method with global-local enrichment functions[END_REF]. The advantage of this method is that the global mesh can be coarse and that the local problems can be solved in parallel as they can be assumed to be independent. This approach is especially useful for propagating cracks because it avoids the need to transfer the solution from one mesh to another, which is costly and inaccurate.

Extended finite element method (X-FEM)

Independently of GFEM, the extended finite element method (X-FEM) was developed in [START_REF] Moës | A finite element method for crack growth without remeshing[END_REF]. In its original version, it was designed to allow simulating brittle fracture in a mesh-independent way. The crack is modelled as a displacement jump across Γ with the Heaviside function H, and the following four functions give the asymptotic displacement field around the crack tip:

{F ℓ (r, θ)} = { √ r sin θ 2 , √ r cos θ 2 , √ r sin θ 2 sin θ, √ r cos θ 2 sin θ } (2.35)
The first of these functions has a branch cut, while the other three are continuous, and (r, θ) is the polar coordinate system attached to the crack tip. The X-FEM approximation of the displacement field is

u = ∑ i∈I ψ i u i + ∑ j∈J ψ j H(Γ)a j + ∑ k∈K ψ k ( 4 ∑ ℓ=1 F ℓ b k ℓ ) , (2.36)
where the first sum is the continuous part of the approximation, while the other two are responsible for the discontinuous parts. The general form of X-FEM with P number of enrichments f j is

u = ∑ i∈I ψ i u i + P ∑ j=1 ∑ k∈J j ψ ψj k f j a jk (2.37)
and we refer to it as standard X-FEM. Here, I is the set of all the nodes in the mesh, and J j ⊆ I. The shape functions used for the enrichments are often taken as the Γ H(Γ)

ψ i ψ i H(Γ) 1 2 3 4 5
ψ 3 ψ 4 h Figure 2
.5: Some of the basis functions in the X-FEM for a 1D configuration standard shape functions:

ψ ψj k = ψ k . It is evident from (2.37
) that the X-FEM falls back to the FEM if there are no enrichments. The enrichment procedure for a single discontinuity is illustrated in Fig. 2.5.

Remark 2.

The Heaviside-enriched X-FEM and the GFEM approximations are the same as shown with the example depicted in Fig. 2.5:

V XFEM = span{ψ i } 5 i=1 + span{(ψ 3 + ψ 4 )H(Γ)} = = {ψ 1 } × {1} + {ψ 2 } × {1} + {ψ 3 } × {1, H(Γ)}+ + {ψ 4 } × {1, H(Γ)} + {ψ 5 } × {1} = 5 ∑ i=1 ψ i V i = V GFEM (2.38)
Compared to the standard X-FEM, the approximation in the intrinsic X-FEM [START_REF] Fries | The intrinsic XFEM: a method for arbitrary discontinuities without additional unknowns[END_REF], yet another version of X-FEM, has the form of

u = ∑ i∈I ψ ⋆ i u i , (2.39)
showing that no additional unknowns are involved as in X-FEM (cf. Eq. (2.37)). Rather, the same enrichments that can be used in X-FEM are incorporated to the basis ψ ⋆ i with the moving least squares technique -hence the name intrinsic. Ill-conditioning in X-FEM applied to fracture mechanics may stem from the Heaviside enrichment and from the singular crack tip enrichments. The former one happens when a crack almost goes through a node. Consider the discontinuity shown in Fig. 2.5. If ϵ = 0, H(Γ) = 0 on the element containing Γ. As the constant 0 is already in span{ψ i } N i=1 (cf. Eq. (2.37)), linear dependency is introduced between the standard part and the Heaviside-enriched part of the approximation. When 0 ̸ = ϵ ≪ h, the condition number becomes very high. The situation is similar in 2D and 3D too. One way to avoid very small cuts is to modify the interface to go through the node [START_REF] Moës | Non-planar 3D crack growth by the extended finite element and level sets -Part I: Mechanical model[END_REF]. This method is applicable if the interface originally passes close to a node, otherwise the modification would significantly change the geometry. Moreover, in the presence of several interfaces in 2D and 3D, small cuts can occur due to junctions and it cannot be cured by this method. Instead of moving the interface to the node, the node can also be shifted to the interface. However, it precludes the use of Cartesian meshes. Junctions pose no problem if the volumetric criterion is applied [START_REF] Daux | Arbitrary branched and intersecting cracks with the extended finite element method[END_REF]. In that approach, a node is not enriched with the Heaviside function (or if already enriched, that degree of freedom is eliminated) if its support is small on one side of the interface with respect to the support on the other side of the interface. A more precise condition is the stiffness criterion introduced in [START_REF] Siavelis | Large sliding contact along branched discontinuities with X-FEM[END_REF]. The second class of methods consists of dealing with the ill-conditioning at the matrix level by applying specialised preconditioners as the Cholesky decomposition used in [START_REF] Béchet | Improved implementation and robustness study of the X-FEM for stress analysis around cracks[END_REF], or a simple diagonal scaling used in [START_REF] Lang | A simple and efficient preconditioning scheme for heaviside enriched XFEM[END_REF]. The third type of techniques for the high condition number is the stable GFEM (SGFEM) [START_REF] Babuška | Stable Generalized Finite Element Method (SGFEM)[END_REF] which modifies the enrichment functions in such a way that they become linearly independent from the finite element shape functions. This works in 1D but is not too effective in 2D and 3D. Yet another technique is the orthogonalization of (some of) the enrichment functions [START_REF] Sillem | The Orthonormalized Generalized Finite Element Method-OGFEM: Efficient and stable reduction of approximation errors through multiple orthonormalized enriched basis functions[END_REF]. The quasi-orthogonalization method in [START_REF] Agathos | Improving the conditioning of XFEM/GFEM for fracture mechanics problems through enrichment quasiorthogonalization[END_REF], applied for two and three-dimensional fracture mechanics problems, removes the linear dependency among the enrichment functions, but not the dependency among the PU functions and the enrichment functions. The latter dependency is provably removed when the so-called flat-top partition of unity functions are used [START_REF] Schweitzer | Generalizations of the Finite Element Method[END_REF]. However, this approach is less suitable for existing finite element codes where the standard finite element shape functions are the PU functions.

Ill-conditioning emerges from the crack tip enrichments as well. The stress singularity at the crack tip is characteristic to the PDE. If only the nodes of the element containing the crack tip are enriched, the behaviour of the singularity is localized to a smaller and smaller region as the mesh is finer and finer. This initial strategy was called topological enrichment in [START_REF] Béchet | Improved implementation and robustness study of the X-FEM for stress analysis around cracks[END_REF] and it results in non-optimal convergence. The optimal order of convergence is retained if the nodes in a disk of fixed radius are enriched. It was named geometrical enrichment [START_REF] Béchet | Improved implementation and robustness study of the X-FEM for stress analysis around cracks[END_REF]. Instead of associating the four crack tip enrichment functions (2.35) to vectorial degrees of freedom, the vectorial enrichment strategy [START_REF] Chevaugeon | Improved crack tip enrichment functions and integration for crack modeling using the extended finite element method[END_REF] associates three vectorial enrichment functions

K 1 = √ r cos θ 2 (κ -cos θ)(e 1 + e 2 ) K 2 = √ r sin θ 2 (κ + 2 + cos θ)e 1 + √ r cos θ 2 (κ -2 cos θ)e 2 K 3 = √ r sin θ 2 e 3 (2.40)
to scalar degrees of freedom, thereby decreasing the number of unknowns while still keeping the same level of accuracy. In Eq. (2.40) κ = 3 -4ν, where ν is the Poisson's ratio.

Tracking (possibly non-straight, 3D) evolving interfaces (e.g. cracks) becomes very difficult if the interfaces are explicitly described. The level set method is an interface capturing technique in which a surface is represented by the zero isocontour of a higher dimensional surface. It was first coupled with the X-FEM in [START_REF] Sukumar | Modeling holes and inclusions by level sets in the extended finite-element method[END_REF] to model voids and inclusions. Its principle is to discretize the level set function (often the signed distance function) on the finite element mesh. Describing a crack which does not completely cut the domain requires two level set functions [START_REF] Stolarska | Modelling crack growth by level sets in the extended finite element method[END_REF]: the normal level set is used to localize the crack discontinuity and the tangential level set gives the position of the crack tip.

The X-FEM has been used for various applications, especially in which weak or strong discontinuities are expected in the solution. It was invented for simulating 2D brittle fracture in [START_REF] Moës | A finite element method for crack growth without remeshing[END_REF]. The use of X-FEM for cracks was already covered in Section 2.2.3. The X-FEM was successfully applied to multiphase fluid and fluidstructure interaction as well. It is a useful tool for computational homogenization problems because there is no need to mesh a complicated microstructure in a unit cell. The first such study along this line was [START_REF] Moës | A computational approach to handle complex microstructure geometries[END_REF] which described the microstructure with level sets and used a periodic mesh to account for the periodicity. In [START_REF] Legrain | An X-FEM and level set computational approach for image-based modelling: Application to homogenization[END_REF] the microstructure was not generated but directly obtained by a level set-based image segmentation process from an existing 3D image. That level set was also used to create the enrichment function, allowing an automatic finite element analysis. Probably the earliest contribution to cope with frictional contact in X-FEM was [START_REF] Dolbow | An extended finite element method for modeling crack growth with frictional contact[END_REF]. Large sliding was considered in [START_REF] Nistor | An X-FEM approach for large sliding contact along discontinuities[END_REF], and [START_REF] Siavelis | Robust implementation of contact under friction and large sliding with the eXtended finite element method[END_REF] combined the previous two approaches using the stable discrete displacement-Lagrange multiplier pair of [START_REF] Béchet | A stable Lagrange multiplier space for stiff interface conditions within the extended finite element method[END_REF]. The model was further improved in [START_REF] Siavelis | Large sliding contact along branched discontinuities with X-FEM[END_REF] by allowing crack branching.

Since the early contributions to GFEM/X-FEM, several papers were published on their computer implementation (see e.g. [START_REF] Bordas | An extended finite element library[END_REF][START_REF] Sukumar | Modeling quasi-static crack growth with the extended finite element method Part I: Computer implementation[END_REF]). In the beginning, they were research codes mostly implemented from scratch. Indeed, the node-based programming approach of X-FEM makes it difficult to integrate into existing FEM software, which is element-based. Researchers then provided X-FEM plugins for proprietary finite element programs such as Abaqus [START_REF] Giner | An Abaqus implementation of the extended finite element method[END_REF]. Finally, as of version 6.9, an official X-FEM module is available in Abaqus, demonstrating the ubiquity, maturity and industrial strength of the method.

Cut finite element method (CutFEM)

Another way to represent a jump in a function is the method first proposed in [START_REF] Hansbo | A finite element method for the simulation of strong and weak discontinuities in solid mechanics[END_REF] and later called under several names: phantom node method [START_REF] Song | A method for dynamic crack and shear band propagation with phantom nodes[END_REF], virtual node method [START_REF] Molino | A virtual node algorithm for changing mesh topology during simulation[END_REF] and cut finite element method (CutFEM) [START_REF] Burman | CutFEM: Discretizing geometry and partial differential equations[END_REF] -all referring to the specific way it is implemented. For the rest of the thesis, we refer to it as CutFEM. The idea of the method is to introduce the discontinuity across an interface Γ by not enforcing the continuity of a function u between the two subdomains Ω + and Ω - separated by Γ. Therefore,

u ∈ V = V + ⊕ V -. (2.41)
In the finite dimensional setting, the finite element spaces V + h and V - h can be constructed either on matching or on non-matching meshes. The non-matching subdomain meshes can either be generated subdomain-wise or cut out from a common background mesh. CutFEM was generalized to multiple intersecting meshes in a recent paper [START_REF] Johansson | Multimesh finite element methods: Solving PDEs on multiple intersecting meshes[END_REF]. Their multimesh finite element method allows arbitrary (possibly moving) meshes to be superimposed on a background mesh. In this way, the individual body parts can be meshed independently.

Since CutFEM uses the standard finite element approximation on Ω + and Ω -, the patches of the PUM are taken as shown in Fig. 2.6a. Note that no patches overlap across Γ. This formulation would need a special construction of the shape functions whose supports are intersected by Γ, i.e. ψ ψ3 , ψ ψ4 , ψ ψ5 and ψ ψ6 in Fig. 2.6a. The explicit construction of these cut shape functions can be avoided by realizing that one requires the integral of these shape functions and not the shape functions themselves. Therefore, the cut-off can be taken into account during the integration by performing the quadrature only on part of a cut element. This observation allows an easy implementation of CutFEM to existing finite element codes. One can see in Fig. 2.6a that the shape functions ψ ψ4 and ψ ψ5 can be associated with the nodes 4 and 5, respectively, which are not part of the physical domain. This is why CutFEM is also called phantom node method. As opposed to the phantom node method, the floating node method [START_REF] Chen | A floating node method for the modelling of discontinuities in composites[END_REF] creates sub-elements from a cut element by local remeshing, see Fig. 2.6b. CutFEM is also affected by ill-conditioning if a cut element part is small compared to the whole element. Figure 2.7 shows such a configuration. In this figure, filled circles represent the original nodes part of each subdomain, while hollow circles denote the phantom nodes. Ill-conditioning originates from the stiffness matrix of the second submesh of Fig. 2.8. We can see that supp ψ 1 ′ ∩ Ω 2 is very small, resulting in an almost linearly dependent basis. In the CutFEM community this is solved by adding additional terms to the weak form, specifically tailored to avoid ill-conditioning in the presence of small cuts. We saw this approach in Section 2.4.3 for the CutFEM method with ghost penalty. Using CutFEM, the discontinuity along the interface Γ is provided by overlapping meshes as depicted in Fig. 2.8.

Both in CutFEM and in the Heaviside-enriched X-FEM, one has to compute integrals on cut elements. There are three main approaches to achieve this [START_REF] Sudhakar | Quadrature schemes for arbitrary convex/concave volumes and integration of weak form in enriched partition of unity methods[END_REF]. The earliest and most intuitive method was used in the first X-FEM paper [START_REF] Moës | A finite element method for crack growth without remeshing[END_REF]. It works by decomposing the cut element part into geometrical primitives (most often triangles in 2D and tetrahedra in 3D) on which quadrature rules are well-known. The technique gives accurate results because the quadrature triangles/tetrahedra are conforming. On the other hand, particularly in 3D and for the integration of high degree polynomials, it results in a lot of quadrature points, which makes it a computationally expensive strategy. Octree-based adaptive quadrature is widely used in the finite cell method (introduced in [START_REF] Parvizian | Finite cell method[END_REF], see [START_REF] Schillinger | The Finite Cell Method: A Review in the Context of Higher-Order Structural Analysis of CAD and Image-Based Geometric Models[END_REF] for a detailed review), which creates a non-conforming subdivision of integration cells near an interface. It is simple to implement but requires a lot of quadrature cells -and therefore many function evaluations -for sufficient accuracy. Moment fitting is a promising way to construct quadrature rules to integrate polynomials on arbitrary convex or concave polyhedra [START_REF] Sudhakar | Quadrature schemes for arbitrary convex/concave volumes and integration of weak form in enriched partition of unity methods[END_REF]. The moment fitting equation is a nonlinear least squares problem, being linear in the quadrature weights and nonlinear in the quadrature point locations.

1 Γ Ω 1 Ω 2 Ω 1 Ω 2 Γ 1 2 3 4
The common way to handle the moment fitting equations is to fix the quadrature points and then solve a linear least squares problem for the weights. So far we saw that the CutFEM approximation uses broken polynomial spaces implemented by taking the restriction of the continuous finite element shape functions on the physical part of the domain. The same idea can be followed to construct an approximation space on a surface by taking the restriction (trace) of the bulk shape functions on the surface. This method was introduced in [START_REF] Béchet | A stable Lagrange multiplier space for stiff interface conditions within the extended finite element method[END_REF] for X-FEM and independently in [START_REF] Olshanskii | A Finite Element Method for Elliptic Equations on Surfaces[END_REF] for surface PDEs. The latter formulation was later coined as TraceFEM [START_REF] Olshanskii | Trace Finite Element Methods for PDEs on Surfaces[END_REF].

For elements containing multiple elaborate discontinuity schemes in 2D or 3D, the use of branch-enrichment in X-FEM is particularly complex due to dependency issues among the basis functions. Conversely, CutFEM generalizes very easily for fully general discontinuity configurations as one just needs to create as many overlapping elements as the number of subdomains present in the element. Another useful property of CutFEM is that all the degrees of freedom have the meaning of function evaluation at a node, as in the standard FEM. In contrast, the coefficients a jk in X-FEM have no direct physical meaning. In CutFEM, all the nodes have the same number of degrees of freedom, making the bookkeeping in the implementation easier. Although the shape functions are different for the two methods, they span the same linear space as it was shown for a single discontinuity in [START_REF] Belytschko | A comment on the article "A finite element method for simulation of strong and weak discontinuities in solid mechanics[END_REF]. This is no longer the case when junctions are present [START_REF] Richardson | An XFEM method for modeling geometrically elaborate crack propagation in brittle materials[END_REF]. CutFEM concentrates more on the domain decomposition nature of the problem, as shown by (2.41), while in the X-FEM more versatile local features can be incorporated, not only discontinuities. Partially cracked elements are also easier to use with X-FEM, although the phantom node method (CutFEM) was successfully applied to that case too [START_REF] Rabczuk | A new crack tip element for the phantom-node method with arbitrary cohesive cracks[END_REF]. One can find publications about the robust implementation of CutFEM as well [START_REF] Burman | CutFEM: Discretizing geometry and partial differential equations[END_REF].

We remind the reader that CutFEM in this thesis is used to represent the kinematics of strong discontinuities and must be combined with constraints to represent continuity (in a weak sense) of the displacements when cracks did not propagate yet. In this chapter the strong and weak formulations of the linear elasticity equations are derived for multiple non-overlapping domains.

The chapter starts with the definition of commonly used terms and notations.

For an initially intact body (containing no cracks), the linear elasticity equation is written in the strong form. To enable potential cracks to appear along every interface, the displacement field must allow a jump across the interfaces. Such a discontinuity can be provided by considering the linear elasticity equations for each subdomain (representing a grain) independently. For this formulation to be equivalent to the one written for the whole domain (as if no cracks were present), it is necessary to impose the continuity of the displacement and the traction fields.

The differential equations on the subdomains, the continuity conditions and the boundary conditions will be weakly imposed using the Lagrange multiplier method, which turns the constrained minimization problem into a saddle point problem. It turns out that the Lagrange multiplier takes the role of the traction vector along 43 interfaces between subdomains. The traction vector will be used in constructing a crack propagation scheme, one of the major reasons why the Lagrange multiplier method was chosen in the thesis. Based on smoothness requirements, we identify the function spaces in which the displacement and the Lagrange multiplier solutions lie to build an a priori stable framework. The scope of this chapter is the construction of local finite element spaces for the two unknown fields: the displacement and the Lagrange multiplier, assessing the stability of the corresponding formulations, and creating robust algorithms.

The meshing of a complicated network of cracks is challenging even in 2D. Therefore, the mesh generation and the function approximation are separated as much as possible in this work. The method uses the cut finite element method (CutFEM) to construct the local finite element spaces for the displacement field on each grain, which are then patched together to form the global finite element space. The Lagrange multiplier space is constructed in a similar vein, with the exception that a reduction is introduced to decrease its dimension and therefore render the mixed method stable. One advantage of CutFEM is its embarrassingly parallel nature. Indeed, very little communication is needed between the subdomains, therefore the geometrical algorithms and the creation and assembly of the local finite element matrices and vectors can be done independently for each grain and interface. The methodology introduced so far is general, but for simplicity and computational efficiency, a global Cartesian mesh is used as the background mesh.

The Lagrange multiplier space reduction algorithm is based on the relative position of the interfaces with respect to the mesh. The stability of the resulting formulation is verified with the inf-sup test, and the patch test proves that the constant stress state can be reproduced by the designed method.

The ease of mesh generation comes at the price of mesh cutting requirements and more complicated quadrature. Cutting the elements with the grains and with the interfaces is a polygon clipping problem in computational geometry. Special care is needed to carry out these cutting operations robustly so as not to introduce kinematic inconsistency. Two detailed algorithms are given to tackle this problem. The first one modifies the geometry, while the second one leads to an optimization problem.

Definitions

Let us consider the domain Ω ⊂ R d (d = 2, 3) partitioned into M non-overlapping subdomains Ω i of polytope shapes, identified by the index set

I s = {1, . . . , M } such that Ω = ⋃ i∈Is Ω i , Ω i ⋂ Ω j = ∅, ∀i, j ∈ I s , i ̸ = j (3.1)
where the overline denotes the closure of a set. Without the loss of generality, it is assumed that M > 1. When we talk about the physical nature of the problem, we will use the word grain instead of subdomain. A straight segment between two subdomains is called an interface and is denoted by

Γ i .
Γ i N (i ∈ I N ⊆ I s ) and Γ i D (i ∈ I D ⊆ I s ) such that Γ i N = Γ N ⋂ Ω i , ⋃ i∈I N Γ i N = Γ N , (3.2)
Γ i D = Γ D ⋂ Ω i , ⋃ i∈I D Γ i D = Γ D . (3.3)
The subdomain boundaries are oriented counter-clockwise and the interface and boundary segments (Γ i and Γ i N , Γ i D ) are oriented according to their outward unit normal attached to Ω i+ . In other words, the orientation of the subdomain determines the normal and tangent unit vectors of its boundary, and an interface inherits these unit vectors from the first neighbouring subdomain it belongs to. These notations are visualized in Fig. 3.1 in a general situation. The intuition behind these definitions is the observation that the boundary conditions and the interface conditions can be handled within a unified framework which is advantageous when the constraint equations are weakly enforced, as will be explained in Section 3.3.

Ω 1 Ω 2 Ω 3 n 2 n 1 n 3 n 4 Γ 3 D Γ 1 D Γ 2 N Γ 4 Γ 3 Γ 1 Γ 2 Ω 1 = Ω 1+ = Ω 2+ = Ω 4- Ω 2 = Ω 3+ = Ω 4+ Ω 3 = Ω 1-= Ω 2-= Ω 3-

Governing equations

In this work, isotropic linear elastic grains are considered. Note that this assumption is not critical and that grain anisotropy could be considered. The equilibrium, constitutive and kinematic equations for these grains are therefore given as Eq. (3.4):

σ • ∇ = 0, x ∈ Ω σ = C : ε, x ∈ Ω ε = ∇ s u x ∈ Ω. (3.4)
To simplify the presentation, body forces are neglected. The Hooke tensor C can vary from grain to grain if the material residing in Ω is heterogeneous. The boundary conditions are written as

u = u D , x ∈ Γ D , σ • n = t N , x ∈ Γ N , (3.5)
where n is the outward unit normal to the boundary. Because of the (cohesive) granular structure of the assumed material, it is straightforward to perform the calculation grain-wise as done in [START_REF] Simone | A Generalized Finite Element Method for polycrystals with discontinuous grain boundaries[END_REF]. To this end, the displacement field is decomposed as

u = ∑ m∈Is χ m u m , (3.6)
where χ m is the indicator function for grain m and is defined as

χ m (x) := { 1 x ∈ Ω m 0 otherwise. (3.7)
The boundary value problem Eq. (3.4)-(3.5) is then written for each grain as:

σ m • ∇ = 0, x ∈ Ω m σ m = C m : ε m , x ∈ Ω m ε m = ∇ s u m , x ∈ Ω m (3.8) u i = u i D , x ∈ Γ i D , (3.9)
σ j • n j = t j N , x ∈ Γ j N , (3.10
)

∀m ∈ I s , ∀i ∈ I D , ∀j ∈ I N .
In addition, we need to provide the continuity of the primary variable (the displacement) and its derivative along the interfaces. Physically, this means prescribing along each interface the continuity of the displacement field

u i+ ⃓ ⃓ Γ i = u i-⃓ ⃓ Γ i (3.11)
and of the interfacial traction

t j+ ⃓ ⃓ Γ j = t j-⃓ ⃓ Γ j
(3.12)

σ j+ • n j ⃓ ⃓ Γ j = σ j-• n j ⃓ ⃓ Γ j (3.13)
Introducing the shorthand notation • j = (•) j+ -(•) j-, the above constraints take the form:

u i = 0, x ∈ Γ i , (3.14) σ j • n j = 0, x ∈ Γ j . ( 3

Mixed continuous formulation

Equations (3.8)-(3.15) are cast into the weak form for subsequent discretization. First, we formally derive the equations and then choose the appropriate function spaces based on mechanical reasoning. The strain energy of the deformable body and the work done by the surface tractions are the sum of the individual contributions on the grains, therefore the energy functional takes the form:

Π b = ∑ m∈Is 1 2 ∫ Ω m σ m (u m ) : ε m (u m ) dΩ - ∑ i∈I N ∫ Γ i N u i+ • t i N dΓ. (3.16)
Constraint equations are required to enforce the continuity equations (3.14)-(3.15). During the discretization, we will see that the finite element nodes do not coincide with the interfaces or the boundaries. Therefore, the continuity equations and the Dirichlet BCs are taken into account weakly. An additional advantage of the weak imposition of Dirichlet boundary conditions is that the system matrix does not need to be modified. The weak solution to (3.8)-(3.15) is found by minimizing the total potential energy (3.16) with respect to the weak constraints

Π c = ∑ i∈I int ∫ Γ i u i • λ i dΓ + ∑ i∈I D ∫ Γ i D ( u i+ -u i D ) • λ i dΓ, (3.17)
where λ i are the vector-valued Lagrange multiplier functions. The constrained optimization problem can be turned into a saddle-point formulation by constructing the Lagrangian Π as

Π = Π b + Π c (3.18)
and taking its variation with respect to the two fields, i.e. the displacement and Lagrange multiplier fields. The variational problem is then (cf. (2.16)):

Find u ∈ V and λ ∈ Λ such that a(u, v) + b(λ, v) = f (v), ∀v ∈ V (3.19a) b(u, µ) = g(µ), ∀µ ∈ Λ (3.19b)
As interfaces and Dirichlet boundaries are handled in the same way, with an abuse of notation, both will be referred to as interfaces. In this way, Γ := Γ i ∪ Γ D is introduced. One can unify the interface and Dirichlet boundary conditions by noticing that boundaries are only surrounded by subdomains from one side. Therefore, the displacement jump on an arbitrary internal interface Γ i or Dirichlet boundary Γ D is redefined as

u i = { u i+ -u i-on Γ i u i+ on Γ i D
The symmetric bilinear and linear forms introduced in

(3.19a)-(3.19b) are expressed by a(u, v) = ∑ m∈Is ∫ Ω m ε m (v m ) : C m : ε m (u m ) dΩ (3.20) b(λ, v) = ∑ i∈I int ∪I D ∫ Γ i ∪Γ i D v i • λ i dΓ (3.21) f (v) = ∑ i∈I N ∫ Γ i N v i+ • t i N dΓ (3.22) g(µ) = ∑ i∈I D ∫ Γ i D µ i • u i D dΓ (3.23)
It can be identified from the equilibrium equation (3.19a) that the Lagrange multiplier represents the traction vector on the interface:

σ i+ • n i+ = -λ i , σ i-• n i-= λ i , (3.24)
where n i+ = -n i-. The corresponding derivation can be found in Appendix A.

In the presence of cracks between the grains, the displacement field is discontinuous. At a junction of interfaces (e.g. when two subdomains have multiple common interfaces), the normal vector to the interfaces changes direction in a discontinuous fashion. Hence, the Lagrange multiplier field representing the tractions across interfaces must be discontinuous as well at such junction points according to (3.24). Therefore, the weak solution (u, λ) ∈ (V, Λ) to the equations (3.8)-(3.15) is sought in the broken Sobolev spaces

V = [ H 1 (Ω) ] d = × i∈Is [ H 1 (Ω i ) ] d (3.25) Λ = [ H -1/2 (Γ) ] d = × i∈I int ∪I D [ H -1/2 (Γ i ) ] d (3.26)
where × is the direct product and d is the spatial dimension.

Approximations

The equations of the bulk (Eq. (3.8)) are formulated on subdomains. This strategy fits well for CutFEM in which the discrete displacement and Lagrange multiplier spaces are created as broken polynomial spaces, i.e. the spaces of piecewise continuous functions of low-degree polynomials, defined on certain submeshes. To create approximation spaces in the bulk and on the surfaces, the continuous finite element basis functions defined on a background mesh are restricted to the bulk and to the interfaces to obtain approximations for the displacement and Lagrange multiplier fields, respectively. In what follows, this intuition is made precise. Generalizing (2.41) for multiple non-overlapping subdomains, the piecewise continuous space V h on Ω can be constructed by patching the subspaces together ( ⨁ denoting the direct sum):

V h = ⨁ i∈Is V i h (3.27)
and hence it holds that u h | Ω i = u i h . A similar decomposition was used in [START_REF] Simone | A Generalized Finite Element Method for polycrystals with discontinuous grain boundaries[END_REF] for polycrystals.

Remark 3. The displacement approximation can be constructed subdomain-wise (see Eq. (3.27)). Therefore, it naturally results in a scheme prone to domain decomposition, and hence to parallelism.

The basis for V i

h is constructed from a background mesh M as follows (the notations can be followed in Fig. 3.2a). The background mesh M is chosen to be large enough such that it contains ∪ i∈Is Ω i . For efficiency, M is chosen as a structured quadrilateral mesh. To approximate the displacement field on subdomain i, we introduce the approximation space constructed as

V i h := span { ψ ψi j | supp ( ψ ψi j ) ∩ Ω i ̸ = ∅, j ∈ M } , (3.28)
where ψ ψi j are the vectorial nodal shape functions cut off on the boundaries of Ω i :

ψ ψi j := { ψ i j on Ω i 0 on Ω \ Ω i . (3.29)
Since the support of the nodal finite element shape functions we use stretch to the neighbouring elements, we define the submesh M i ⊆ M, consisting of those elements e which are at least partially cut by Ω i :

M i := {e ∈ M | e ∩ Ω i ̸ = ∅}. (3.30)
Therefore, a basis for V i h can be constructed from the nodal basis functions corresponding to the mesh M i . It is then sufficient to consider

V i h := span{ψ ψi j | j ∈ M i }. (3.31)
As explained in Section 2.5.4, the cut shape functions ψ ψi j are replaced by the standard shape functions ψ i j and the integrals are evaluated on Ω i . Note that without constraint equations, u h is discontinuous between any two subdomains

Ω i ̸ = Ω j . M i M j M Ω i Ω j
(a) Submeshes for the bulk Similarly to (3.30), the submesh containing the interface Γ i is defined as (see Fig. 3.2b)

M i Γ M Γ i (b) Submesh for the surface
M i Γ := {e ∈ M | e ∩ Γ i ̸ = ∅}. (3.32)
The discrete Lagrange multiplier space, following [START_REF] Béchet | A stable Lagrange multiplier space for stiff interface conditions within the extended finite element method[END_REF], is assumed to be

Λ i h = span { ψ ψi j | j ∈ M Mi Γ } , (3.33)
where M Mi Γ is the same mesh as M i Γ but not every node of it defines an independent basis function. The Lagrange multiplier basis functions are created as the weighted trace of the nodal basis functions on the interface:

ψ ψi j = ∑ k∈M i Γ α jk ψ k | Γ i . (3.34)
The global approximation space is therefore

Λ h = ⨁ i∈I int ∪I D Λ i h . (3.35)
Our task is then to find M Mi Γ and α jk for all interfaces such that a stable formulation is obtained. An algorithm to accomplish this goal is detailed in Section 3.6.1. The construction of the discrete spaces V h in (3.27) and Λ h in (3.35) allows seeking the approximation on each subdomain and on each interface independently.

Discrete weak form

The saddle point problem (3.19a)- (3.19b) is now discretized. Using the spaces (3.31) and (3.33), the approximations are given by

u i h = ∑ k∈M i ψ i k U i k , v i h = ∑ j∈M i ψ i j V i j , i = 1, . . . , M λ i h = ∑ k∈M Mi Γ ψ ψi k L i k , µ i h = ∑ j∈M Mi Γ ψ ψi j M i j , i = 1, . . . , N (3.36)
where

U i k , V i k , L i k , M i k are
⎧ ⎪ ⎨ ⎪ ⎩ Find (u h , λ h ) ∈ V h × Λ h such that a(u h , v h ) + b(λ h , v h ) = f (v h ), ∀v h ∈ V h b(u h , µ h ) = g(µ h ), ∀µ h ∈ Λ h (3.37)
with the discretized symmetric bilinear and linear forms expressed as

a(u h , v h ) = ∑ i∈Is ∑ j,k∈M i V i j K i jk U i k , (3.38) b
(u h , µ h ) = ∑ i∈I int ∑ j∈M Mi Γ ∑ k∈M i± M i j B i jk (U i+ k -U i- k ) + ∑ i∈I D ∑ j∈M Mi Γ ∑ k∈M i+ M i j B i jk U i+ k , (3.39)
f (v h ) = ∑ i∈I N ∑ j∈M i+ V i+ j F i j , (3.40)
g(µ h ) = ∑ i∈I D ∑ j∈M Mi Γ M i j G i j , (3.41)
in which M i± are the two submeshes corresponding to the neighbouring subdomains of interface Γ i .

Remark 4. Note that V h ̸ ⊂ V and Λ h ̸ ⊂ Λ as the geometry is taken into account during the quadrature. As a result, our discretization is non-conforming, making the stability more difficult to achieve and prove.

The stiffness matrix on Ω i , with ε i ( ψ i j ) = ∇ s ψ i j , is given by

K i jk = ∫ Ω i ε i ( ψ i j ) : C ( E i , ν i ) : ε i ( ψ i k ) dΩ, (3.42)
the coupling matrix is written as

B i jk = ∫ Γ i ∪Γ i D ψ ψi j • ψ i+ k dΓ (3.43)
and the nodal force vectors are

F i j = ∫ Γ i N t i N • ψ i j dΓ (3.44) G i j = ∫ Γ i D u i D • ψ ψi j dΓ. (3.45)
Note, that we use ψ i+ k in (3.43) as it is the same as ψ i- k . The assembled matrices have the saddle-point structure

[ K B ⊤ B 0 ] [ u λ ] = [ f g ] , (3.46)
in which u contains the unknowns U i k and λ contains the unknowns L i k . The algebraic saddle point system (3.46) is uniquely solvable if the following conditions hold [START_REF] Benzi | Numerical solution of saddle point problems[END_REF] ker

B ⊤ = 0, (3.47) ker K ∩ ker B = 0 (3.48)
Equation (3.47) expresses that the Lagrange multipliers are linearly independent, while (3.48) is the discrete version of the condition (2.20). We emphasize that the second equation is a relaxed version of K being positive definite, which would not hold in our case. In fact, the null space of K is spanned by the rigid body modes of the individual grains.

Stable mixed formulation

Although a unique solution exists if the equations (3.47)-(3.48) hold, a robust numerical method must be stable too. Consider a discrete linear operator L mapping from the solution space to the data space (image of L). The continuous dependence of the solution on the data is given by the inverse mapping L -1 . To quantify stability, L -1 is measured in a suitable norm [START_REF] Arnold | Stability, Consistency, and Convergence of Numerical Discretizations[END_REF]. This definition applied to mixed formulations of the type (2.16) has the form (2.22). For a conforming or non-conforming discretization with discretization parameter h, (2.22) becomes

∥u h ∥ V + ∥λ h ∥ Λ ≤ c∥f h ∥ V ′ + ∥g h ∥ Λ ′ . (3.49)
It is relevant that stability is a stronger concept than uniqueness. If (3.49) holds uniqueness is granted but the reverse is not true. One must take a series of discretizations and (3.49) must be satisfied for c independent of h.

In the first part of this section, an algorithm is developed to create the Lagrange multiplier space. The second part discusses the stability of the resulting mixed method.

Construction of the Lagrange multiplier space

Some additional notions used in the creation of a stable Lagrange multiplier space are now defined. A cut edge is an element edge cut by an interface. A cut element is an element the interior of which is fully or partially cut by an interface. An intersection point is the intersection of an interface with an edge. An isolated node is a node of a cut element that does not lie on a cut edge. A junction is the intersection of interfaces. The connectivity graph is an undirected graph, the vertices of which are the cut edges, and in which two of such vertices are connected if the corresponding cut edges meet at a common node.

These concepts are represented in Fig. 3.3. The bottom of Fig. 3.3 shows the connectivity graph of the intersection points. The vertices of this graph are circled. We can see that vertex 1 of the graph matches a cut edge which does not connect to any other cut edges. Therefore, vertex 1 is isolated in the graph (not to be confused with the isolated nodes). This graph is used in this section to construct a stable interpolation for the Lagrange multiplier field.

In order to avoid interface meshing, which can be complicated in 3D, the Lagrange multiplier unknowns are defined at the nodes of the background finite element mesh as in [START_REF] Béchet | A stable Lagrange multiplier space for stiff interface conditions within the extended finite element method[END_REF] (see also [START_REF] Olshanskii | A Finite Element Method for Elliptic Equations on Surfaces[END_REF]). The corresponding question that needs to be addressed is how to distribute the Lagrange multiplier degrees of freedom to the nodes of the mesh M Mi Γ . If each node of M Mi Γ held an independent Lagrange multiplier, the dimension of Λ i h would be too high. Therefore, a reduction is performed among them by associating a certain Lagrange multiplier to a group of nodes. The identification of these groups is done not element-by-element but at the global level, i.e. regarding all nodes of M Mi Γ . The groups are determined according to topological information, i.e. the location of the nodes relative to the interface Γ i . Once these groups are formed, the interpolation within an element is built by taking the restriction of the bulk element basis functions on Γ i . An example is provided in Fig. 3.4, where the Lagrange multiplier on element e is interpolated as Triangle elements with linear interpolation were used in [START_REF] Béchet | A stable Lagrange multiplier space for stiff interface conditions within the extended finite element method[END_REF] with a single interface. As a consequence, the trace of these linear shape functions resulted in a linear interpolation for the corresponding Lagrange multiplier along the interface. The procedure was next extended to quadratic triangles in [START_REF] Ferté | Interface problems with quadratic X-FEM: design of a stable multiplier space and error analysis[END_REF]. In [START_REF] Jin | Crack Lip Contact Modeling Based on Lagrangian Multipliers with X-FEM[END_REF] the idea of [START_REF] Béchet | A stable Lagrange multiplier space for stiff interface conditions within the extended finite element method[END_REF] was used with an additional crack tip enrichment. They concluded that both the standard tip enrichments (2.35) and the vectorial tip enrichments (2.40) applied with the Heaviside enrichment pass the inf-sup test, provided that constant interpolation is used for the Lagrange multiplier field on the element containing the crack tip.

λ e = ψ 2 | Γ i L i 1 + (ψ 1 + ψ 3 + ψ 4 )| Γ i L i 2 .
The present work generalizes this type of approach by considering the following further improvements: (i) general multiple interface configurations, furthermore (ii) incorporated within bilinear quadrilateral elements forming a Cartesian mesh. Before outlining the algorithm for creating the shape functions for the Lagrange multipliers, let us look at the requirements to be fulfilled by the Lagrange multiplier space. The rationale behind it is that the best approximation property is achieved upon polynomial completeness. We, therefore, seek a discretization that can at least represent a uniform stress field Σ on Ω. This can be achieved if the material is homogeneous (all grains have the same elastic properties, E and ν) and linear elastic. Setting the Poisson ratio to 0, if a rectangular specimen is clamped on its left edge and pulled horizontally from the right, the displacement field in Cartesian coordinates is u = (αx + β)e x .

(3.50)

Since this analytical expression for u is continuous, u ij = 0, and (3.19b) is automatically fulfilled. After substitution of (3.50) into (3.37) and integration by parts, we can observe that

• as Σ is uniform, the interface traction and therefore the Lagrange multiplier shape functions must form a partition of unity. To be able to represent a uniform stress field with the interpolation, due to the property that ∑ 4 k=1 ψ k = 1 on an element, all nodes of M Mi Γ must contain a Lagrange multiplier.

• The integrand of ∫

Γ N v • t N dΓ (3.51)
is discontinuous even if the traction is uniform on Γ N because the test functions

v m ∈ [H 1 (Ω m )]
d are discontinuous across the interfaces. Therefore, if we use the standard Gauss-Legendre quadrature and want to integrate exactly, we must integrate on Γ N segment-by-segment as (3.40) and (3.44) suggest.

• As noted at the end of Section 3.3, the Lagrange multiplier field is discontinuous at junctions.

• Analogously to mortar methods, in which the junction points carry no Lagrange multipliers [START_REF] Puso | Mesh tying on curved interfaces in 3D[END_REF], we use a constant interpolation for the Lagrange multiplier on an element containing a junction.

• In order to avoid over-determination, Dirichlet boundaries are handled in the same way as internal interfaces.

First, a method for building a discrete Lagrange multiplier space is given in case of a single interface in Ω.

To avoid defining a too rich Lagrange multiplier space that causes the instability, some of the multipliers are linked together and are prescribed to have the same values. The nodes of the cut elements on the two sides of an interface are connected. A straightforward approach that consists of using all the cut edges to link the nodes on the two sides of an interface fails to satisfy the LBB condition in the general case. This wrong, naive strategy is illustrated in Fig. 3.5. In such a strategy, the number of tyings/constraints is too high with respect to what is strictly needed to enforce the continuity of the fields. Therefore, a reduction of the Lagrange multiplier space has to be performed similarly to [START_REF] Béchet | A stable Lagrange multiplier space for stiff interface conditions within the extended finite element method[END_REF] but tailored to quadrilateral elements. This is done by first selecting some cut edges. The selection of the set of these cut edges to reach a stable approximation is not unique hence we have some freedom in choosing them. The following four rules are applied in this selection procedure based on the connectivity graph illustrated in Fig. 3.3.

1. A cut edge corresponding to an isolated vertex of the connectivity graph is always selected.

2. A selected cut edge is not allowed to be connected to any other selected cut edge in the connectivity graph.

3. A non-selected cut edge must be connected to at least one selected cut edge in the connectivity graph.

4. A cut edge is never selected if the interface cutting it has a junction in an element containing the cut edge.

The third rule makes the discrete Lagrange multiplier space as large as possible but due to the second rule not too large to avoid dependencies. The fourth rule is responsible for not introducing independent Lagrange multipliers at a junction. Finally, the first rule guarantees that all nodes of the cut elements contain a Lagrange multiplier, thereby satisfying the partition of unity requirement. Setting the discrete Lagrange multiplier space is performed according to Alg. 3.1 (follow also in Fig. 3.6). The links (ellipses in Fig. 3.5) can also be represented as a graph, shown in Fig. 3.6.

L i 1 L i 3 L i 1 L i 2 , L i 3 L i 2 Γ i

Algorithm 3.1 Linking Lagrange multipliers

1: Nodes of an element which contains a junction are linked together to fulfil the partition of unity 2: Determine the set of selected edges from the cut edges not part of the already processed edges 3: Link the two nodes of a selected edge 4: The remaining node of a cut edge is linked the other node of the cut edge 5: An isolated node is linked to the opposite node of the cut element it belongs to With these chosen links, the Lagrange multiplier shape functions inside the elements are built as the trace of the nodal (bulk) displacement shape functions on the interface. They are • constant if there is a change in the normal vector in the underlying element • linear if the interface passes through opposite edges

• quadratic without the linear term if the interface passes through neighbouring edges (cf. the third plot in Fig. 3.6)

Remark 5. The Lagrange multiplier shape functions constructed in this way are not local to the considered element, their support can extend to the neighbouring elements.

Note that the use of the intersection points on the selected cut edges would not constitute a robust solution to define a basis for Λ i h . This is due to the fact that

• the shape functions need to be created globally and not at the element level since their support extends to neighbouring elements,

• for non-straight interfaces this process would be even more difficult

• the extension would become very complicated for 3D configurations in which surface interfaces have to be considered, requiring global 2D shape functions to be defined.

Remark 6.

The algorithm proposed in [START_REF] Hautefeuille | Robust imposition of Dirichlet boundary conditions on embedded surfaces[END_REF] forms piecewise continuous Lagrange multiplier shape functions, mentioning that the computation of the weights can be done element-wise. In our case, the weights are also calculated element-wise and the created shape functions are globally C 0 continuous except inside those elements where the normal vector to the interfaces changes in a discontinuous fashion. Moreover, they are piecewise C 1 continuous (see the third plot of Fig. 3.6).

For multiple interfaces, the procedure described above is performed on each of them, independently from each other. This is illustrated in Fig. 3.7 where the Lagrange multiplier field is constant on the element containing the junction. The proposed general algorithm Alg. 3.1 for the selection of Lagrange multipliers is now represented for different cases in Fig. 3.6, where the cut edges selected to define them are drawn in green, the other non-selected cut edges are depicted in blue, and the isolated nodes are connected in red. Figure 3.6 depicts this general case in which the shape function ψ 5 is used to interpolate both L 1 and L 2 . If we chose ψ 5 to interpolate only λ 1 for example, then ψ ψ2 would be constant on the element with nodes 1-2-5-4. In Fig. 3.7 the selection procedure is shown in presence of a triple junction. As mentioned before, the building of the discrete Lagrange multiplier space is done interface-wise, therefore the method is represented on one of the three interfaces (cf. first figure of Fig. 3.7). It is shown in the third figure of Fig. 3.7 that the Lagrange multiplier field is constant on the element containing the junction.

A remaining situation to tackle is when a subdomain is fully embedded inside an element. There is then no intersection point on the interfaces. As all the sides of this embedded subdomain (i.e. its interfaces) are inside an element, a constant Lagrange multiplier field is used on all these sides by linking all four nodes of the element together.

In conclusion, the process of mesh tying and stiffness assembly (in the linear elastic case) is the following. 2. Compute the stiffness matrices on each subdomain (possibly in parallel).

L 1 L 1 L 2 L 2 4 9 8 1 2 L 1 , 1 2 L 2 1 2 L 1 , 1 2 L 2 L 1 L 2 A B C D A B C D 0 1 ψ ψ1 ψ ψ2
3. Select some of the cut edges to define independent Lagrange multipliers.

4. Create the Lagrange multiplier space according to Alg. 3.1.

5. Calculate the coupling integrals in Eq. (3.43).

As the Lagrange multipliers are defined at the nodes of the background mesh, it is not immediately apparent what their values should be on an interface. Associating a geometrical point on an interface to a Lagrange multiplier is however important because physically relevant quantities (traction vector, damage, energy release rate) are defined on interfaces. As we will see in Chapter 4, the computations do not need to know about these geometrical points, but useful insight can be gained by plotting such a physical quantity on an interface. Let us use the variable i to index a vectorial Lagrange multiplier DOF. A reasonably chosen geometrical point P i corresponding to an index i is called location and is determined by Alg. 3.2.

Stability of the discretization scheme

We now assess the stability of the designed mixed formulation. Two kinds of tests are performed in this part: the patch test and the inf-sup test.

Algorithm 3.2 Determine the geometric location of a Lagrange multiplier 1: Find the interface Γ i the Lagrange multiplier indexed by i is associated to 2: Fetch all the nodes N i,k which hold the Lagrange multiplier indexed by i 3: For all N i,k find the closest points N * i,k on Γ i 4: Take the mean of the coordinates of N * i,k , the resulting point P i lies on Γ i

The patch test is used to evaluate the polynomial completeness of the finite element. For the elasticity equation of compressible solid, the Q1 quadrilateral elements satisfy the patch test, i.e. the linear displacement solution can be reproduced exactly (assuming no under-integration). We investigate here whether the Lagrange multipliers can ensure perfect tying across the interfaces, i.e. zero displacement jump can be achieved. This procedure is applied considering an assembly of ten subdomains, generated using a Voronoi tessellation with uniformly distributed random seeds on the unit square (Fig. 3.8a). Zero displacement is prescribed on the left side of the square and uniform normal traction is applied on the right side. The elastic properties of the grains were collectively set to E = 1 and ν = 0. As there is no contraction, the vertical displacement component must be zero, and we expect a linear variation of the horizontal displacement along x. Figure 3.9 shows that the included discontinuities are successfully deactivated by the Lagrange multipliers and that the numerical result matches with the exact solution (3.50). Hence, the patch test on a Cartesian mesh is passed by the formulation. As an additional verification test, the numerically determined strain energy is compared with the exact value, available in this simple case. Since the strain and stress fields are uniform, the exact strain energy is

U = 1 2 ∫ Ω σ : ε dΩ = 1 2 ΩΣ : E = 1 2 , (3.52)
while the numerical value obtained is

1 2 u ⊤ Ku = 0.499999999999997. (3.53)
The computation of the domain integrals on polygonal element parts is handled by triangulation (more details on quadrature can be found in Appendix B). It can be seen in Fig. 3.10 that only the cut elements are decomposed into triangles, standard Gauss-Legendre integration is performed on the uncut bilinear elements. We note that the subdivision of an element into triangles is done only for the sake of computing the integrals in the stiffness matrices (3.42), no additional degrees of freedom are introduced.

The normal and tangential tractions are directly calculated from the Lagrange multipliers knowing the interface normals. Note that the traction is also computed on the Dirichlet boundary Γ D (Fig. 3.11). It can be seen in these figures that, as expected, the traction vector is close to zero on interfaces with an orientation parallel to the direction of the loading. On Γ D , t t = 0 and t n = 1. We now focus on the inf-sup conditions. In case of mixed problems of the form (2.16), the finite dimensional version of the inequalities (2.20)-(2.21) are

Computing ∫ Γ λ • u dΓ = λ ⊤ B(u + -u -) (3.54)
0 < α ⋆ ≤ inf 0̸ =u h ∈Z h sup 0̸ =v h ∈Z h a(u h , v h ) ∥u h ∥ V ∥v h ∥ V , (3.55)
0 < β ⋆ ≤ inf 0̸ =λ h ∈Λ h sup 0̸ =v h ∈V h b(λ h , v h ) ∥v h ∥ V ∥λ h ∥ Λ (3.56)
with α ⋆ and β ⋆ being independent of the mesh size h [START_REF] Brezzi | On the existence, uniqueness and approximation of saddle-point problems arising from Lagrangian multipliers[END_REF]. The second condition (3.56) is often called the inf-sup condition. We will refer to it as LBB3 condition (which is another commonly used name) to distinguish it from the first one we will refer to as coercivity condition.

Remark 7. The continuity of a(•, •) and b(•, •) need not be checked, it automatically comes from the fact that we used the proper norms (see [START_REF] Boffi | Mixed Finite Element Methods and Applications[END_REF]).

Checking the fulfilment of the pair of conditions (3.55)-(3.56) is very difficult for most BVPs. Even if the stability was proved for the continuous problem, it would not imply the stability of the discrete problem due to the non-conforming discretization. The inf-sup test provides practical means to determine the LBB constant. The original version of the test was created for constraints prescribed on domains [START_REF] Chapelle | The inf-sup test[END_REF]. For interfacial constraints, the test was slightly modified in [START_REF] Bathe | The inf-sup condition and its evaluation for mixed finite element methods[END_REF]. Its idea is to express the LBB constant with the help of an eigenvalue problem. The square root of the smallest non-zero eigenvalue of

{ Find β ∈ R and 0 ̸ = (u h , λ h ) ∈ V h × Λ h such that ∀(v h , µ h ) ∈ V h × Λ h (u h , v h ) V ;Ω + b(λ h , v h ) + b(u h , µ h ) = -β(λ h , µ h ) Λ;Γ . (3.57)
is identified as the LBB constant:

β ⋆ = √ β min . In our problem, V i = [H 1 (Ω i )] 2 , ∀i ∈ I s and Λ i = [H -1/2 (Γ i )] 2 , ∀i ∈ I int ∪ I D .
As it is difficult to deal with fractional Sobolev spaces, they are replaced by mesh-dependent L 2 -norms as in [START_REF] Béchet | A stable Lagrange multiplier space for stiff interface conditions within the extended finite element method[END_REF]:

∥µ h ∥ 2 -1/2;Γ := ∑ e∈E h e ∥µ h ∥ 2 0,e , (3.58
)

∥v h ∥ 2 1/2;Γ := ∑ e∈E 1 h e ∥v h ∥ 2 0,e , (3.59)
where E is the one-dimensional mesh on Γ and h e is the size of element e in E. We note that in our Lagrange multiplier selection algorithm, no such situation arises that h e is arbitrarily close to zero, whatever the position of the interface and the mesh used. This property and the norm equivalence allows replacing the one-dimensional mesh segments h e with the background mesh size h. We also use norm equivalence between the H 1 -norm ∥ • ∥ 1;Ω and the energy norm ∥ • ∥ E;Ω := √ a(•, •) to replace (u h , v h ) 1;Ω with a(u h , v h ), therefore being able to use the already available stiffness matrix. Using these norms, the matrix representation of the eigenvalue problem (3.57) reads:

[ K B ⊤ B 0 ] [ u λ ] = β [ 0 0 0 hG ] [ u λ ] , (3.60)
where the stiffness matrix K and the coupling matrix B are already defined in (3.42) and in (3.43). The Gram matrix G = diag G i kℓ is block diagonal, and the block matrix corresponding to interface i is formed as

G i kℓ = (ψ ψi k , ψ ψi ℓ ) Γ i .
The assembly of G i kℓ is similar to that of B i kℓ . It is a common practice to solve (3.60) either for u or for λ. As K is singular, we chose to eliminate λ from (3.60) to obtain

1 h ( B ⊤ G -1 B ) u = βKu. (3.61)
A series of generalized linear eigenvalue problems (3.61) are solved on increasingly refined meshes. Then the inf-sup constant is approximated as

β ⋆ = min h i √ β h;min . (3.62)
Computing β h;min on many meshes increases the reliability of the inf-sup test. However, obtaining all the eigenvalues for matrices of size K is very demanding. This is even so because Matlab's eig function can compute all the eigenvalues of a generalized eigenvalue problem for full matrices only. Iterative methods are able to return only a selected few eigenvalues. The strategy used here is therefore the following. Assume that our discretization is stable. Then the inf-sup constant should not decrease significantly. Therefore, we determine β H,min out of all the eigenvalues for a coarse mesh with mesh size H. In the model problem depicted in Fig. 3.12a, H was selected to be 1/32, and the full spectrum could quickly be obtained for that mesh size on a laptop. For meshes with h < H, we searched N eigenvalues in the vicinity of β H,min , keeping the matrices in sparse storage format. If we had only looked for one eigenvalue, the iterative method might have found a zero eigenvalue, corresponding to a rigid body mode of K.

Z h = ker b(v h , µ h ) := {v h ∈ V h | b(v h , µ h ) = 0, ∀µ h ∈ Λ h }. (3.63)
The smallest in modulus eigenvalue of the generalized eigenvalue problem

{ Find α ∈ R and 0 ̸ = u h ∈ Z h such that ∀v h ∈ Z h a(u h , v h ) = α(u h , v h ) V ;Ω (3.64)
gives the Brezzi coercivity constant: α ⋆ = |α min |. In most problems, Z h is not known explicitly. An alternative eigenvalue problem is

{ Find α ∈ R and 0 ̸ = (u h , λ h ) ∈ V h × Λ h such that ∀(v h , µ h ) ∈ V h × Λ h a(u h , v h ) + b(λ h , v h ) + b(u h , µ h ) = α(u h , v h ) V ;Ω . (3.65)
It was proved in [6] that the smallest in modulus eigenvalue of (3.65) and the smallest in modulus eigenvalue of (3.64) agree. This time we cannot replace (u h , v h ) V ;Ω with a(u h , v h ) as we did before because all the eigenvalues would be infinity or 1. Therefore, the new term is directly evaluated as follows. By the definition of the broken Sobolev spaces and given that in our problem

V i = [H 1 (Ω i )] 2 for i ∈ I s , we have (u h , u h ) 1;Ω = ∥u h ∥ 2 1;Ω = ∑ i∈Is ∥u h ∥ 2 1;Ω i = = ∑ i∈Is ( ∥u h ∥ 2 0;Ω i + |u h | 2 1;Ω i ) = u ⊤ Su.
(3.66)

The subdomain matrices are formed as

S i kℓ = (ψ i k , ψ i ℓ ) Ω i + (∇ψ i k , ∇ψ i ℓ ) Ω i , i ∈ I s , (3.67)
where the global matrix is S = diag S i kℓ and u is the global vector of the displacement unknowns. The assembly of the subdomain matrices S i kℓ and the global matrix S is done exactly as for the stiffness matrices. However, more Gauss points were used because S contains the shape functions as well, not only their derivatives as K. The eigenvalue problem (3.65) with these matrices is

[ K B ⊤ B 0 ] [ u λ ] = α [ S 0 0 0 ] [ u λ ] (3.68)
The calculation of |α min | is performed with the same method as used for Eq. (3.61). Figure 3.12b shows that for the investigated meshes the Brezzi coercivity constant is also bounded away from zero. Although numerical evidence is not equivalent to a proof, the patch test and the careful examination of the inf-sup constants give us confidence about the stability of the discretization scheme.

Robustness

Mesh-independent methods, such as X-FEM, alleviate the mesh generation phase at the expense of some additional care that must be taken. Ill-conditioning coming from small cuts is an issue with X-FEM-like methods as was discussed in Section 2.5.3. In this section two strategies are proposed, both having their advantages and drawbacks.

Interface adaptation

Care must be taken when an interface passes through a node or even coincides with an edge. Such situations have to be considered if the interface is closer to a node than a given tolerance. When an interface goes through a node (i.e. the intersection point coincides with a node), one should decide which cut edge to use for building the tyings. Although the intersection point could be associated with the node it coincides with, that special case would make the implementation more difficult. Therefore, we always associate an intersection point to an edge. One can always select a cut edge containing that node and independent of the other selected edges, showing that the special node always holds an independent Lagrange multiplier. However, to handle this in a general setting, we rather choose to move the intersection point along one of the four edges to which the node belongs. This is denoted in the sequel as geometry adaptation. The corresponding algorithm is given in Alg. 3.3. Steps 3 to 6 of Alg. 3.3 are easy to understand but steps 1 and 2 require some explanation. Two strategies can be followed. The first one is step 3 of Alg. 3.3. The second strategy consists of steps 1-3 of Alg. 3.3. Assuming an interface that almost goes through a node (dashed line in Fig. 3.13), two connecting edges to the centre node are cut. According to the first strategy, the two intersection points are moved away from the node (Fig. 3.13a, circles). The second strategy first moves the intersection points to the centre node, deletes one of them, and finally shifts the remaining one along an edge (Fig. 3.13b, middle circle). The newly created interface then cuts the other edge (Fig. 3.13b, square). Both strategies work well in this case.

Algorithm 3.3 Geometry adaptation

1: Intersection points close to a node are moved to that particular node 2: Intersection points being at a node, stemming from the same interface, are deleted 3: Intersection points sitting at a node are moved away from that node on an edge 4: Construct new interfaces in place of those, whose intersection points were moved 5: Modify the subdomains according to the newly included interfaces 6: Find the new intersection points which emerge as the intersection of the new interfaces with the mesh Now let us consider the situation with an interface going exactly through a node (Fig. 3.14). For the depicted configuration, four coinciding intersection points appear on the edges incident to that node. According to the first strategy, these four intersection points are moved along those edges, away from the node. The circles in Fig. 3.14a show the displaced intersection points. The question is how to reconstruct the modified interface. There is no unique choice as it was in the previous case (cf. given in Fig. 3.14a. One can notice that this method results in a very jagged interface configuration. Conversely, the second strategy works equally well as in the first case (see Fig. 3.14b) because the number of intersection points is reduced to one, tracing back this configuration to the previously discussed setup. This second strategy is therefore selected.

Optimal background mesh generation

The interface adaptation method of the previous section has some drawbacks. Modifying the interfaces based on the current background mesh introduces mesh-dependency to the otherwise mesh-agnostic CutFEM. The other downside of this strategy is that creating new, shorter interfaces without mesh refinement results in a worse approximation of the Lagrange multipliers, which in turn deteriorates the accuracy of the traction vector. Finally, it is not so simple to implement, especially when considering future 3D extensions.

These deficiencies gave the idea to search an alternative, possibly simpler method.

Modifying the interfaces served two purposes. First, it avoids the special cases in the mesh connectivity, allowing every possible interface-mesh intersection being carried out with a common computational geometry algorithm. Second, it allows suppressing the extreme ill-conditioning caused by the small cut areas. The aim is therefore to create a background mesh for which small cuts do not occur. This can be achieved by keeping the geometry fixed and moving the mesh optimally. This suggests the existence of an optimization problem to determine the best position of the mesh. Solvers evaluate the objective function for such an optimization problem at certain points as part of an optimization algorithm. If the objective function is too costly to evaluate, this new technique becomes unfeasible in practical situations with many subdomains and a fine mesh. Therefore, instead of cutting the elements of the mesh with the subdomains, an alternative measure is sought which characterizes the cut part. In 2D, tiny cut element parts are present in the following two cases:

(i) an interface is close to a node (Fig. 3.15a),

(ii) an interface endpoint is close to an edge (Fig. 3.15b)

The goal is therefore to maximize the distance among interfaces and nodes, and among interface endpoints and element edges, which leads to an optimization problem. This heuristic approach is much cheaper than if an element-subdomain cutting was executed. Such a measure based on distances clearly cannot provide good results if a polygonal subdomain is very elongated, but in such cases, a finer or graded mesh would anyway be required. As a single small element part causes a high condition number for the global system matrix, it is not sufficient to maximize the sum of the square-distances. Our aim is to maximize the minimum distance between each of the geometrical entities involved (i.e. nodes, edges, interfaces). The general optimization problem is therefore formulated as: find the location x of the nodes of the mesh such that

max x ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ min i∈I j∈N d(I i , n j (x)) min i∈I j∈E k=1,2 d(P ik , e j (x)) , such that Ω Ω ⊇ Ω, (3.69)
holds, where P ik are the two endpoints of interface I i indexed by I, e j are edges indexed by E, while n j are nodes of the mesh indexed by N . The constraint expresses that the domain of the background mesh must cover the domain of the geometry, containing the interfaces. Task (3.69) can be simplified. Thinking of the optimization problem as moving the background mesh in horizontal (x) and vertical (y) directions, and taking into account that we deal with a uniform Cartesian mesh, it is easy to see that the translation of the mesh must be in the range [-h x , h x ] × [-h y , h y ], where h x and h y are the mesh sizes in x and y directions. These intervals give lower and upper bound constraints to the optimization problem.

As the distance calculation needs to know only about the interfaces and the element nodes and edges, the algorithm takes the interfaces and the mesh as inputs. The background mesh stretches on a rectangular domain determined by the axis-aligned bounding box (AABB) of the set of interfaces. Let this box be given by [a, b] × [c, d], so that its lower-left corner is (a, c). As the mesh can translate by h x and h y distances at most, one layer of elements is placed along the left, right, bottom and top sides of the mesh defined on the AABB of the interfaces. The lower-left corner node is called reference point. The location of all the nodes are determined by this reference point and the uniform mesh size h x , h y . Hence, the background mesh at any iteration of the solver is completely determined by the number of elements in x and y directions, being user inputs. The next step of the algorithm is the distance computation. It is not necessary to compute the distance between every pair of entities of the mesh and the interfaces; most of these distances will be greater than the minimum distance. Therefore, it is enough to consider the nodes and edges of the neighbouring elements to a given interface. As we do not wish to cut the elements with the interfaces, the candidate elements to an interface are taken as the elements covering the AABB of the interface plus an additional layer of elements around them -the same strategy we already used in the initial background mesh construction. As only the position of the candidate elements change when the mesh is moved (this was the goal), the ID of their edges and nodes are stored. This offline part makes it possible to skip the identification of candidate edges and nodes whenever the optimization code invokes the distance computation function.

To solve the maximin problem (3.69), the fminimax solver of Matlab's Optimization Toolbox was used. As we need to solve a maximin problem, the negative of the objective function was passed to the solver. There are only two scalar degrees of freedom (translation modes), so the objective function has two variables: x and y. When finding the optimum, the lower and upper bound constraints must be met: x ∈ [a-2h x , a], y ∈ [c -2h y , c]. The initial guess was set to (x 0 , y 0 ) = (a -h x /2, c -h y /2). The algorithm is assessed by two examples. We are interested in how well the proposed heuristic strategy works. The condition number of the block matrix in Eq. (3.46) is estimated with Matlab's condest function both for the initial and for the optimal mesh. To eliminate the ill-conditioning of the saddle point problem coming from the different scales present in matrices K and B, we choose the material properties such that the entries of K are of the same order of magnitude as the entries We also measure the time required for generating the optimal mesh. In the first test, the unit square domain is considered with an interface separating it into two subdomains. The coordinates of the interface endpoints were randomly chosen. Assume that the user wants a 2 × 2 mesh on the domain on average. The reference point of the initial mesh is at (-0.25, -0.25), shown in Fig. 3.16a. We can see that the interface endpoint on the top side of the domain is close to an edge. The cut element parts are of course not small in this case; this example merely wants to visualize the working of the algorithm. The output of the optimization procedure is the optimal mesh, shown in Fig. 3.16b. The optimal mesh in this simple test example could have been found by "hand". Indeed, the initial mesh should not be moved vertically, it is already in perfect position along this direction. On the other hand, it must be translated horizontally such that its distance from the closest two vertical mesh lines becomes equal. It demonstrates that the algorithm acts reasonably, at least for one interface.

To study the effect of mesh refinement on the conditioning and on the elapsed time, equispaced meshes with 2 i × 2 i number of elements were used for i = 1, . . . , 6. The results are collected in Tab. 3.1. Analysing the condition numbers, the following consequences are apparent:

• the optimal mesh always results in a smaller condition number than the initial mesh

• the optimal mesh is consistent in the sense that the condition number monotonically increases with mesh refinement as opposed to the initial mesh in which the condition number depends on how the mesh is positioned with respect to the geometry

• for an unfortunate cutting situation, even a coarse initial mesh can result in huge condition numbers By measuring the elapsed time, it turned out that the cost of generating an optimal initial mesh is negligible compared to creating the finite element spaces V i h and Λ i h . In case of multiple interfaces, moving the mesh away from an interface can bring it closer to another interface, making the mesh positioning dependent on the relative location of the interfaces. It is therefore expected that the algorithm performs worse than in the single interface case. To assess the capabilities of our algorithm in such a situation, a more realistic example is taken. A circle of diameter 150 mm is tessellated with a bounded Voronoi diagram consisting of 100 seeds (see Fig. 3.17). The tessellations were produced by the VoronoiLimit function [START_REF] Sievers | VoronoiLimit[END_REF].

Mesh

Several mesh sizes are considered to see how the optimization strategy performs with respect to mesh refinement. Table 3.2 reports the obtained data. The previous observations for the condition numbers are valid here as well, the situation being more severe in this example. The computational cost is mild compared to other operations (mesh-subdomain intersection, assembly, solving). Based on these two examples, the following conclusions can be drawn.

• The method is purely geometric, does not need any information about the particular PDE. Therefore, it can be applied in fictitious domain methods, X-FEM, CutFEM, etc.

• This method alone is not sufficient to handle the ill-conditioning coming from the small cuts.

• Compared to the interface adaptation method, the mesh optimization technique is simpler and naturally extends to 3D, as it is based on simple distance calculations does not modify the geometry and hence * it retains the mesh independence of the discretization * the Lagrange multiplier approximation is not deteriorated cannot reduce the condition number as effectively because the relative position of the mesh and the geometry is taken into account globally and not locally

The current algorithm could be improved by the following future developments.

• Instead of computing the condition number of the saddle point matrix, compute the effective condition number of the stiffness matrix. This allows completely getting rid of the possible ill-conditioning created by the saddle point structure.

• The squared distance between a point and a line gives a quadratic function.

The derivative of this objective function is computable analytically, therefore the Jacobian matrix can be given by the user instead of relying on the finite difference computation by the minimax solver.

• The distance computations are independent and hence easy to parallelize.

• An assessment with different initial mesh positions could be used to try finding a global optimum.

However, these enhancements are not yet implemented; this method should be considered as a proof-of-concept in its current state. The designed optimal background mesh strategy is purely geometric (independent of the PDE), not complicated to implement (only requires distance calculation and a black-box minimax solver) and the core idea does not depend on the spatial dimensions (the objective function differs in 1D, 2D and 3D due to the different cause of small cut parts, but the general idea remains the same).

Remark 8. Note that the randomly generated Voronoi cells may have very short edges. This is a difficulty for mesh generators because even if they manage to produce a mesh, it can be of low quality or excessively refined. In [START_REF] Abdelaziz | Grain based modelling of rocks using the combined finite-discrete element method[END_REF] the edges with length under a specified threshold were removed by a two-step algorithm. In our method, small edges do not make the simulation fail. It is true that the Lagrange multiplier approximation is poor along very short edges but that can be remedied by a local non-conforming mesh refinement. We did not experience any numerical issues in the presence of short Voronoi edges.

Brazilian test configuration

The flexibility of the discretization is demonstrated here on the numerical solution for a geometrical configuration matching a classical experimental test. The Brazilian test, proposed in [START_REF] Carneiro | A new method to determine the tensile strength of concrete[END_REF] is an indirect tension test to determine the tensile strength of brittle materials. A cylinder is compressed along its diameter. Assuming ideal conditions (no friction between the loading device and the cylinder, concentrated force applied, homogenous material), an analytical formula exists for the tensile stress σ in the centre of the specimen [100]:

σ = 2P πDt , (3.70)
where P is the applied load, D is the diameter of the cylinder and t is its thickness. In the numerical model D = 5 mm, E = 40 MPa, ν = 0.36, f t = 6 MPa were chosen and a plane strain assumption was considered. The circle is approximated by a regular polygon of 200 sides. The discretized circle is clamped at the bottom along one polygon side and is subjected to distributed force in the vertical direction at the top. Both boundary conditions are prescribed along one polygon side of the discretized circle. These short boundary segments are marked in red in Fig. 3.18a. The cracking is assumed to start along the diameter of the circle along which the loading is applied (blue segment in Fig. 3.18a). Therefore, we create one interface separating the circle into two half-circles, shown in blue in Fig. 3.18a. The simulation was performed on a series of finer and finer meshes made up of the same number in horizontal and vertical directions. The exact critical force for crack initiation obtained from (3.70) as P crit exact = (1/2)πtDf t is compared with the numerically computed one P crit num , determined by the maximum normal traction on the interface. Equivalently, this boils down to estimating computationally the load at which a given stress is obtained at the centre of the specimen. The relative error is shown in Fig. 3.18b for different mesh densities. One can see that a sufficiently fine mesh is necessary to achieve a good result and that the relative error 10 -3 seems to be a plateau. This behaviour could be explained by the fact that the exact solution corresponds to a different BVP than the one considered in this model.

It is also interesting to study the distribution of the normal traction on the vertical diameter of the circle with mesh refinement. Assuming point loading, the normal traction under both plane stress and plane strain conditions is given by the following formula [START_REF] Hung | Theoretical analysis and digital photoelastic measurement of circular disks subjected to partially distributed compressions[END_REF]:

t n = P πRt ρ 8 + 4ρ 4 -4ρ 2 -1 + 2(-2ρ 6 + ρ 4 + 1) (ρ 2 -1) 4 , (3.71)
which is very close to a constant function. This constant function for the critical load P crit is the tensile strength f t . The non-dimensional variable ρ denotes r/R with r being the radial polar coordinate and R = D/2. The numerically computed traction profiles are shown in Fig. 3. [START_REF] Bathe | An inf-sup test for shell finite elements[END_REF]. One can see in the zoomed regions that for finer and finer meshes the numerical solution approaches the analytical one. A perfect match is not possible because the circle is approximated with a polygon and the loading is not prescribed at a single point. This logically results in large compressive stresses close to the perimeter of the circle. This example shows that polygonal domains other than squares, seen in the previous sections, can simply be embedded into a background mesh, and lead to a proper stress distribution in the simulation. We note that the current test case serves a verification purpose and is not a validation with experimental results. In Section 5.2, the Brazilian test with multiple grains will be studied, following the complete softening of the material. This chapter starts with the recognition that both from the kinematic and from the material behaviour points of view, the introduction of a cohesive zone is required. Then an unconventional damage-driven fracture solution procedure, instead of the more common load-controlled computations, is developed.

Diffuse microcracking in cohesive grain-based materials

For the more accurate enforcement of the contact conditions and to avoid the ill-conditioning coming from an intrinsic cohesive zone formulation with a high initial elastic stiffness, an extrinsic traction-separation relation is assumed. Following a recent work, the completely failed cohesive elements are not removed from the mesh, rather a modified Lagrangian formulation is applied in which only finite terms appear in the compliance matrix upon complete failure. This approach greatly increases the robustness of the method.

Based on available experimental evidence, the Mohr-Coulomb failure criterion is used in this work with a cut-off in the tensile regime and a cap in the compression region.

The extrinsic cohesive zone model is reparametrized with an internal variable -the interfacial damage -which (i) separates the crack propagation from the determination of the mechanical fields (displacement and Lagrange multiplier), (ii) ensures the irreversibility of the crack growth, (iii) allows an explicit treatment of both crack growth and contact enforcement in the cohesive zone and (iv) provides an energetic criterion for both crack initiation and crack growth. The damage variable is interpolated with the same basis functions as the Lagrange multiplier which makes the crack propagation easy: there is no need for tracking the crack tip.

The modified formulation includes some additional terms in the weak form. The experimentally observed increase of the mode II fracture energy with compression on interfaces is taken into account. However, it implies that the damage parameter no longer corresponds to the level of degradation. Therefore, a new, physically meaningful quantity related to energy dissipation is defined, through which the damage field is updated.

Finally, the solution scheme of the crack propagation problem is summarized.

Crack propagation along interfaces

To propagate cracks in the proposed discretization scheme, it seems natural to deactivate some constraints (Lagrange multipliers) to let the interface open based on some energy release principle. However, several aspects need to be tackled to reach this objective. When simply fully deactivating a tying constraint, the fictitious node method indeed creates a discontinuity along the whole element at once. Consequently, the crack propagates along the whole element length. Moreover, if a Lagrange multiplier is set to vanish, it affects the entire support of its basis function, meaning that the crack will (partially) open on the neighbouring cut element as well. A remedy for this would consist of using piecewise-constant Lagrange multiplier approximations on the elements. A crack could then be opened on each element separately by deactivating the corresponding Lagrange multipliers without affecting the neighbouring element. However, this would deteriorate the accuracy of the Lagrange multiplier approximation; and the piecewise constant Lagrange multiplier approximation is not stable when used with a linear displacement interpolation [START_REF] Puso | Mesh tying on curved interfaces in 3D[END_REF].

Conversely, one could open a crack on a whole interface, but this would yield an imprecise propagation for large grains. Opening a crack at a junction is also another delicate issue in such cases as it results in a kinematic inconsistency when two crack segments meet, while shifting the junctions towards element edges would introduce significant geometrical errors and mesh dependency. An equally important aspect is the length scale. In quasi-brittle materials, the nonlinear zone around the crack tip is not negligible. This motivates us to use a method which does not break an interface at once but provides the possibility of the progressive opening of a crack. The cohesive model allows convenient implementation, smoother crack opening and a better estimation of the released energy from a traction-separation relation. This enables tracing the softening response in a quasi-brittle failure context. To avoid the costly incrementaliterative procedure most often used in cohesive approaches, an explicit solution procedure based on a dissipation-driven scheme will be used to circumvent the need to solve iteratively nonlinear problems.

Quasi-static crack simulations can be performed by fixing the load (traction or displacement loading) and propagating the cracks by a distance stemming from equilibrium solving. Alternatively, it can be implemented by propagating the cracks by a certain distance and deducing the required load for this propagation to have equilibrium [START_REF] Ferté | 3D crack propagation with cohesive elements in the extended finite element method[END_REF]. The latter approach enables the use of an explicit algorithm, i.e. to set the crack opening increments at the beginning of a step and to deduce the load factor from this information using equilibrium. This principle will be used in the sequel.

Mixed method for cohesive cracks

With at hand a mixed discretization involving displacement and Lagrange multiplier unknowns, an extrinsic CZM can be used. The use of a mixed formulation involving displacements and tractions will also prove useful to enforce contact conditions precisely. To this end, the two-field formulation of [START_REF] Cazes | A two-field modified Lagrangian formulation for robust simulations of extrinsic cohesive zone models[END_REF] is used within the CutFEM setting instead of mesh-conforming cohesive elements.

To simplify the notations, we omit the interface and subdomain indices. The weak form (3.19a)- (3.19b) then reads:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ Find (u, λ) ∈ V × Λ such that ∀(v, µ) ∈ V × Λ ∫ Ω ε(v) : C : ε(u) dΩ + ∫ Γ v • λ dΓ = ∫ Γ N v • t N dΓ, ∫ Γ µ • u dΓ = ∫ Γ D µ • u D dΓ. (4.1)
The cohesive traction λ( u ) (which is a purely softening function for an extrinsic formulation) is introduced in the weak formulation in its compliance form (see [START_REF] Cazes | A two-field modified Lagrangian formulation for robust simulations of extrinsic cohesive zone models[END_REF]) as follows:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ Find (u, λ) ∈ V × Λ such that ∀(v, µ) ∈ V × Λ ∫ Ω ε(v) : C : ε(u) dΩ + ∫ Γ v • λ dΓ = ∫ Γ N v • t N dΓ, ∫ Γ µ • ( u -u (λ)) dΓ = ∫ Γ D µ • u D dΓ. (4.2)
The second equation of (4.2) compares the displacement jump calculated from the displacement field with the displacement jump obtained from the cohesive model based on the cohesive traction. Until the activation of a cohesive zone, Eq. (4.1) and Eq. (4.2) are the same. Using this formulation, an incremental-iterative solution procedure would be required to solve this nonlinear problem. Since the opening is computed from the Lagrange multiplier in the compliance form as u = Rλ, some terms in the compliance matrix R tend to infinity when λ tends to 0 (i.e. failure of the cohesive zone). This was identified in [START_REF] Cazes | A two-field modified Lagrangian formulation for robust simulations of extrinsic cohesive zone models[END_REF], which motivated the use of a modified Lagrange multiplier method. An effective cohesive traction ζ is therefore introduced as a Lagrange multiplier defined by

ζ = λ + k • u , (4.3)
where

k = [ k n 0 0 k t ] (4.4)
with the parameters k n and k t being positive. This allows evaluating the compliance matrix even when λ = 0. In Section 4.7, lower bounds will be determined for these parameters. From (4.2), using (4.3), one obtains a new weak form expressed as:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ Find (u, ζ) ∈ V × Λ such that ∀(v, η) ∈ V × Λ ∫ Ω ε(v) : C : ε(u) dΩ + ∫ Γ v • (ζ -k • u ) dΓ = ∫ Γ N v • t N dΓ - ∫ Γ D v • k • u D dΓ, ∫ Γ η • ( u -u (ζ)) dΓ = ∫ Γ D µ • u D dΓ.
(4.5) The opening of cracks can cause other cracks to close due to stress redistribution. A contact condition must, therefore, be provided to avoid interpenetration of the crack faces. During contact the conditions u n = 0 and t n ≤ 0 hold. Therefore, according to Eq. ( 4.3), ζ n ≤ 0. To enforce this in Eq. (4.5), the interface stiffness is only decreased when the cohesive zone is in loading state:

∫ Ω ε(v) : C : ε(u) dΩ + ∫ Γ v • (ζ -k • u ) dΓ = ∫ Γ N v • t N dΓ - ∫ Γ D v • k • u D dΓ, ∫ Γ η • ( u -χ(ζ n > 0) u (ζ)) dΓ = ∫ Γ D η • u D dΓ, (4.6)
where χ is the indicator function acting only in normal direction to activate contact upon closure.

In Section 4.4, the mixed formulation (4.6) is recast with the introduction of interfacial damage. It will provide a damage propagation criterion and will allow for a segregated solution procedure, i.e. the system solving for (u, ζ) is separated from the damage update.

Failure criterion

Mode I crack propagation was mainly considered in previous efforts [START_REF] Lé | Coupling damage and cohesive zone models with the Thick Level Set approach to fracture[END_REF]. Experiments are however often carried out with triaxial tests for rock and geomaterials. The modelling of such tests involves compressive-shearing loading on interfaces, and a pure mode I crack propagation model would not be sufficient for such cases. Furthermore, even if uni-axial tension is applied to heterogeneous materials, the slanted interfaces undergo shearing as well. Assuming only mode I propagation is therefore not justified for our purpose.

A model uncoupling failure along the normal and tangential directions is to be avoided as residual normal (tangential) tractions should not be present after the interface has failed in tangential (normal) direction.

A single Mohr-Coulomb criterion would predict an overestimated uni-axial tensile strength. Such a criterion is therefore combined with a Rankine criterion to construct the Mohr-Coulomb criterion with tension cut-off (see e.g. [START_REF] Ottosen | The Mechanics of Constitutive Modeling[END_REF]). Expressed in the space of normal and tangential cohesive tractions, t n and t t , the initial damage criterion consists of the equations

t n -f t = 0, tan φ t n + |t t | -c = 0, (4.7)
giving the tension cut-off and the Mohr-Coulomb parts of the criterion, where f t is the tensile strength, c is the cohesion and φ is the angle of friction. The criterion is depicted in Fig. 4.1a.

t t t n f t c φ (a) tension cut-off t t t n f t c φ ψ -f c (b)

tension cut-off and compression cap

t n -f t = 0, tan φ t n + |t t | -c = 0, (t n + f c ) tan ψ -|t t | = 0, (4.8)
where f c is the compressive strength. In [START_REF] Schellekens | Computational Strategies for Composite Structures[END_REF], the compression cap was an ellipsoid. A compression cap will be shown necessary so that the parameter k to be selected for defining the augmented Lagrangian remains finite for all traction component combinations (see Section 4.7).

The failure criterion gives what combination of the tractions causes damage initiation and therefore the start of softening.

p 1 = f t -f c t t t n q 1 q 2 q 3 q 4 p 2 q 4 q 3 t t u t t t G t G t G t q 1 t t t n p 2 G n t t q 2 G t u t u t u t u n -p 34

Damage-based cohesive zone formulation

In what follows, we consider an extrinsic cohesive model (i.e. without any elastic behaviour before softening), suitable for initially perfect interfaces. Following [START_REF] Lé | Coupling damage and cohesive zone models with the Thick Level Set approach to fracture[END_REF], the CZM is reformulated with an interfacial damage variable d. This allows formulating evolution equations, which simplifies the numerical implementation of the CZM. Using only one damage variable both for the normal and tangential directions ensures that whenever the interface is fully failed, both the normal and the tangential cohesive tractions vanish. The free energy functional is defined as

ϕ( u , d) = 1 2 ( 1 d - 1
) u • k • u (4.9)
and k is given in Eq. (4.4). Another internal variable y is defined as the dual of d, so that evolution laws can be written. The state laws, i.e. the expressions of the cohesive traction t and of the cohesive energy release rate y, respectively, are derived from the potential (4.9) according to:

t = ∂ϕ ∂ u = ( 1 d - 1
) k • u , (4.10)
y = - ∂ϕ ∂d = 1 2d 2 u • k • u . (4.11)
Equation (4.10) justifies the choice of the Helmholtz free energy expression: for a vanishing damage, the interface stiffness is infinite, i.e. an extrinsic CZM is obtained. We do not define the equivalent displacement jump and the equivalent traction as it is usually done for mode-mixity. The reason is that the failure criterion discussed in Section 4.3 mostly resides in the compression regime where u = u t e t and so no equivalent displacement jump is required.

Remark 9. It is noted that using the damage to describe the traction-separation relation was already used in [START_REF] Cuvilliez | A finite element approach coupling a continuous gradient damage model and a cohesive zone model within the framework of quasi-brittle failure[END_REF], but in another context. They used the damage to track the degradation of the bulk and when a certain damage value was reached, they switched to a cohesive zone formulation. In our case, the damage d exists on the interfaces and parametrizes the cohesive law in the whole softening regime (i.e. from intact interface (d = 0) to the completely failed cohesive zone (d = 1)).

The Lagrange multiplier λ in the pure Lagrange multiplier formulation (4.2) was identified as the traction vector: t = λ. The augmented Lagrange multiplier ζ in Eq. (4.3) does not carry a specific meaning. On the other hand, using equations (4.3) and (4.10), one recognizes the modified cohesive law

ζ = 1 d k • u . (4.12)
Determining the fracture energy -discussed in Section 4.5.2 -requires the components of the traction vector. However, only u and ζ are available after solving (4.6). The previous equations allow us to determine the traction vector as

t = (1 -d)ζ. (4.13)
In this context, the Lagrange multiplier ζ can be thought of as the effective traction vector. With the introduction of the damage variable, the damage indicator χ from the second equation of (4.6) can be computed as

χ := { 1, ζ n > 0 or ζ n ≤ 0 and d = 0 0, otherwise (4.14)
Written in terms of the normal and tangential components, and taking into account contact, Eq. (4.13) is computed as

t n = (1 -χd)ζ n , t t = (1 -d)ζ t . (4.15)
Substituting u (ζ) from Eq. (4.12) into Eq. (4.6) defines the final set of equations that will be used in the discretization procedure:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ Find (u, ζ) ∈ V × Λ such that ∀(v, η) ∈ V × Λ ∫ Ω ε(v) : C : ε(u) dΩ + ∫ Γ v • (ζ -k • u ) dΓ = ∫ Γ N v • t N dΓ - ∫ Γ D v • k • u D dΓ, ∫ Γ η • ( u -χdk -1 • ζ ) dΓ = ∫ Γ D η • u D dΓ.
(4.16) The second term of the integrand of the second equation of Eq. (4.16) written in the basis of the coordinate system local to the interface should be interpreted as

[ k -1 n 0 0 k -1 t ] [ d 0 0 d ] [ χζ n ζ t ] ,
in which the contact condition makes the problem nonlinear. To solve the global system in a linear fashion, at a given damage step s, ζ (s-1) is used as an approximation in (4.16).

The discretization of equations (4.16) results a similar matrix as in (3.46). The CutFEM discretization of the displacement and the Lagrange multiplier fields does not require a mesh on the interface. To keep the robustness, we use the same basis for the damage interpolation as for the Lagrange multiplier interpolation:

d i = ∑ j∈M Mi Γ D i j ψ ψi j , (4.17)
where ψ ψi j is a component of ψ ψi j . The new discretized formulation then reads:

a(u h , v h ) = ∑ i∈Is ∑ j,k∈M i V i j K i jk U i k + ∑ i∈I int ∑ j,k∈M i± ( V i+ j -V i- j ) P i jk ( U i+ k -U i- k ) + + ∑ i∈I D ∑ j,k∈M i+ V i+ j P i jk U i+ k , (4.18) b
(u h , η h ) = ∑ i∈I int ∑ j∈M Mi Γ ∑ k∈M i± E i j B i jk (U i+ k -U i- k ) + ∑ i∈I D ∑ j∈M Mi Γ ∑ k∈M i+ E i j B i jk U i+ k , (4.19) c(ζ h , η h ) = ∑ i∈I int ∑ j,k∈M Mi Γ E i j C i jk Z i k , (4.20)
f (v h ) = ∑ i∈I N ∑ j∈M i+ V i+ j F i j + ∑ i∈I D ∑ j∈M i+ V i+ j H i j , (4.21)
g(η h ) = ∑ i∈I D ∑ j∈M Mi Γ E i j G i j , (4.22)
where the new matrices and vector are given by

P i jk = - ∫ Γ i ∪Γ i D ψ i j • k -1 • ψ i k dΓ, (4.23)
C i jk = - ∫ Γ i ∪Γ i D d i ψ ψi j • k -1 • ψ ψi k dΓ, (4.24)
H i j = - ∫ Γ i D u i D • k • ψ i j dΓ. (4.25)
The resulting algebraic system reads

[ K + L B ⊤ B C ] [ u z ] = [ f g ] , (4.26)
where the stiffness matrix K, the coupling matrix B and the vector g are the same as in Section 3. Remark 10. The Lagrange multiplier field ζ is approximated in the same way as the original Lagrange multiplier λ. This is made possible because both λ and u (ζ) are defined in the same basis (as opposed to u , which is given in the nodal basis).

Remark 11. The matrix L couples displacement degrees of freedom from neighbouring subdomains. It implies that the block-diagonal structure of K is not preserved in the system matrix.

The matrix of k in the {e n , e t } basis attached to the interface is diagonal:

k = [ k n 0 0 k t ] . (4.28)
Note that k n and k t do not represent by any means the elasticity of the CZM, as they are only used to deal with completely open cracks. Equations (4.16) are discretized in the standard basis {e x , e y } to obtain (4.26). However, the cohesive model is formulated in the local basis {n, t}, aligned with the interface. One must, therefore, express k and k -1 in the global frame. As the quadratic form u

• k • u is invariant, one has u ⊤ (x,y) k (x,y) u (x,y) = u ⊤ (n,t) k (n,t) u (n,t) . (4.29)
The change-of-basis matrix Q is constant for a given straight interface and is given by

[ n t ] = [ n x n y -n y n x ] [ e x e y ] = Q [ e x e y ] . (4.30)
Equations (4.29)-(4.30) yield

k (x,y) = Q ⊤ k (n,t) Q. (4.31)
and similarly in the matrices C and P:

k -1 (x,y) = Q ⊤ k -1 (n,t) Q. (4.32)
At the start of the simulation, assuming no initial damage (C = 0) and no penalty term (L = 0), the stable u-λ approximation described in Section 3.6 ensures an accurate calculation of the interfacial traction field. Note that only the C matrix contains the damage variable. Therefore, only this matrix needs to be updated during crack propagation, with all other matrices to be computed only once. Furthermore, for efficiency reasons, the integral of the shape functions in (4.24) can be precomputed so that only the multiplication by the nodal damage values is required during the update part. More specifically, using the damage interpolation (4.17), matrix (4.24) can be written as

C i jk = ∑ ℓ∈M Mi Γ C Ci jkℓ D i ℓ , (4.33)
where One can see that both the local non-conforming mesh construction and the local assembly are completely independent for each subdomain and interface, allowing straightforward parallelization. Create the background mesh M i for the subdomain Ω i 8:

C Ci jkℓ = - ∫ Γ i ∪Γ i D ψ ψi j • k -1 • ψ ψi k ψ ψi ℓ dΓ (4.
Find and store

∅ ̸ = Ω i e = Ω i ∩ e 9:
for e ∈ M i do 10:

K i e ← (ψ, ψ) Γ i e (3.42) 11:
K i ← element-to-local(K i e ) 12:
end for

K ← local-to-global(K i ) 14: end for 15: ▷ Construct the matrices defined on the interfaces 16: for i ∈ I int ∪ I D do 17:

Create the background mesh M i Γ for the interface Γ i 18:

Find and store

∅ ̸ = Γ i e = Γ i ∩ e 19:
for e ∈ M i Γ do 20:

T i e ← (ψ, ψ) Γ i e 21:

P i e ← R i P • T i e ▷ R i P : reduction operator to get (4.23)

22: end for 26:

B i e ← R i B • T i e ▷ R
f i e ← R i f • T i e ▷ R
g i e ← R i g • T i e ▷ R
P ← local-to-global(P i ) 27: B ← local-to-global(B i ) 28: f ← local-to-global(f i ) 29:
g ← local-to-global(g i ) 30: end for 31: return K, P, B, f , g

Computing the energetic quantities

Given the damage state, the discrete displacement and Lagrange multiplier values are obtained by solving the generalized saddle point problem (4.26). The new damage state at each point is computed from the evolution laws given by the classical Kuhn-Tucker relations

d ḋ ≥ 0, (4.35) y -y c (d) ≤ 0, (4.36)
(y -y c (d))d ḋ = 0, (4.37)
prescribing damage growth as an irreversible process, and expressing that damage does not increase when the energy release rate remains under a critical value characteristic of the interface. Damage propagation occurs when the critical energy release rate is reached. In (4.36)-(4.37), y c (d) is an increasing function of d, so that an increasing energy is required to further increase the damage on an interface. Section 4.6 deals with the discretization of (4.35)-(4.37). Depending on the normal and tangential traction components at a given location, four characteristic loading modes can be distinguished. The notations can be followed in Fig. 4.2. The remainder of this section is separated into three parts. First, the energy release rate is determined, then assumptions on the mixed mode fracture energy are introduced and finally, the critical energy release rate is computed.

Energy release rate

As k is diagonal according to Eq. (4.4), the energy release rate in Eq. (4.11) can be rewritten as

y = 1 2d 2 k ( u n 2 + u t 2 ) , (4.38)
where the assumption k = k n = k t was made. Note that the displacements are available in the mesh nodes, therefore u n and u t are not directly available on the interface. Taking their restriction on the interface results in an oscillatory approximate displacement field. Therefore, we want to express u n and u t with the help of the Lagrange multipliers because they are free of oscillations. This is doable because the Lagrange multipliers are in connection with the tractions, which are in turn known from the cohesive law. Since the softening in normal and tangential directions is separated according to (4.10), we have

t n = ( 1 d -1 ) k u n , t t = ( 1 d -1 ) k u t . (4.39)
Based on this, the energy release rate can be decomposed to normal and tangential parts as

y = y n + y t , y n = 1 2k ( max{0, t n } 1 -d ) 2 , y t = 1 2k ( t t 1 -d ) 2 , (4.40)
where we postulated that no energy is released in pure compression. Although y is not defined for d = 1 above, it does not pose any problem because a completely failed cohesive zone cannot release energy anyway. In the implementation, this issue can be avoided by noticing that d can be eliminated from y by using (4.13):

y = y n + y t , y n = 1 2k (max{0, ζ n }) 2 , y t = 1 2k ζ 2 t . (4.41)

Decomposing the mixed mode fracture energy

We make the hypothesis that the mixed mode fracture energy G c can be additively decomposed into a normal and a tangential part (G n and G t , respectively):

G c = G n + G t , (4.42)
where G n and G t are determined from the failure criterion.

For geomechanical applications, the mixed mode fracture energy under compressive stress states increases with compression. Therefore, a model inspired by [START_REF] Lourenço | Computational strategies for masonry structures[END_REF][START_REF] Van Zijl | Computational Modelling of Masonry Creep and Shrinkage[END_REF] is applied, imposing a linear increase of the fracture energy as a function of the compressive stress on the interface. The expressions of G n and G t , respecting (4.42), are collected in Tab. 4.1 for the four possible loading modes. For a general loading mode, the point (p, q) on the failure surface gives the strength values. Here, G I and G II are the pure mode I and mode II fracture energies, and a > 0 is a coefficient which sets the growth of the mixed mode fracture energy with increasing compressive confinement on the interface. As a simplification, it is assumed that once the damage initiates at a location, the fracture energy for that location does not change with t n if t n further evolves in the course of degradation. loading mode (see Fig. 4.2)

G n G t 1 ( 1 -q c ) G I q c G I 2 ( 1 -q c ) G I G II -c-q c-qc ( G II -qc c G I ) 3 0 G II -at n 4 0 G II -ap 34

Critical energy release rate

Given a cohesive relation expressed in traction-separation form t( u ), the critical energy release rate y c (d) is determined such that the potential-based cohesive model is equivalent to the physical cohesive relation. The linear extrinsic cohesive relation can be given by two parameters, the strength and the fracture energy (see Fig. 4.2):

t n = p - p 2 2G n u n , (4.43)
t t = q - q 2 2G t u t . (4.44)
Equating them with the tractions derived from the free energy in Eq. (4.39), the cohesive openings u n and u t can be isolated. They are then substituted to Eq. (4.38), and using the fact that y = y c during damage propagation, the critical energy release rate is obtained:

y c (d) = ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ k 2 [ ( A n (1 -d)k + B n d ) 2 + ( A t (1 -d)k + B t d ) 2
] tension and shear k 2

( f t (1 -d)k + f 2 t /(2G I )d ) 2 pure tension k 2 ( c (1 -d)k + c 2 /(2G II )d ) 2 pure shear k 2 ( A t (1 -d)k + B t d ) 2
compression and shear (4.45) where A n = p, B n = p 2 /(2G n ), A t = q, B t = q 2 /(2G t ). The restrictions on the value of k can be found in Section 4.7. Table 4.2 summarizes the computation of the energetic quantities for the four loading modes. quantity loading mode (see Fig. 4.2)

1 2 3 4 G n ( 1 - q c ) G I ( 1 - q c ) G I 0 0 G t q c G I G II - c -q c -q c ( G II - q c c G I ) G II -at n G II -ap 34 y n ζ 2 n 2k ζ 2 n 2k 0 0 y t ζ 2 t 2k ζ 2 t 2k ζ 2 t 2k ζ 2 t 2k y c,n k 2 ( A n (1 -d)k + B n d ) 2 with A n = p, B n = p 2 2G n y t,n k 2 ( A t (1 -d)k + B t d ) 2 with A t = q, B t = q 2 2G t

Table 4.2: Energetic quantities for the different loading modes

It is noted that the proposed damage formulation automatically handles damage initiation except in the tension+shear case. To realize this, note that when the damage grows from 0, y(d = 0) = y c (d = 0). Using Eq. (4.40) and Eq. (4.45), this yields

⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ t 2 n + t 2 t = p 2 + q 2 tension
and shear t n = f t pure tension t t = c pure shear t t = q compression and shear (4.46) Except in the tension+shear region, Eq (4.46) is the stress criterion for crack initiation.

Damage update algorithm

The crack propagation is governed by the evolution equations (4.35)-(4.37). Upon crack initiation or crack extension, certain amount of energy dissipates, governed by the damage variable d. To explain the shortcoming of using d for energy dissipation, let us introduce the so-far dissipated energy at any location:

G : [0, 1] → [0, G c ] defined as G(d) = d ∫ 0 y c (d d) dd d. (4.47)
Upon complete failure, G(1) = G c . The dissipated energy increment is

∆G = d+∆d ∫ d y c (d d) dd d. (4.48)
Substituting y c from (4.45) into (4.48), we realize that ∆G depends on the strength and fracture energy parameters p, q, G n , G t . It means that the energy dissipation is in general not the same for two locations i and j even if d i = d j and ∆d i = ∆d j . In case of a pure Rankine failure criterion (without Mohr-Coulomb failure in compressive states), and only one interface, as considered in [START_REF] Lé | Coupling damage and cohesive zone models with the Thick Level Set approach to fracture[END_REF], the strength and the fracture energy is spatially uniform for all locations and therefore the dissipated energy at a given step is proportional to the damage increment. This is no longer the case for multiple interfaces having different material parameters or for non-uniform fracture energy on a single interface. This would cause that with a pure control by the damage d, the majority of the energy dissipation may not be attributed to the favourable locations.

To cure this behaviour, a physically meaningful quantity, the so-far dissipated energy G introduced above (rather than the damage), is chosen to describe the degradation of the cohesive interfaces, in order to drive the computation. We wish to express the evolution equations with G. Since G ′ (d) = y c from (4.47) and y c > 0, Eq. (4.35) is replaced by

G Ġ ≥ 0. (4.49)
By parametrizing y c with G, the other two evolution equations (4.36) and (4.37) are

y -y c (G) ≤ 0, (4.50)
(y -y c (G))G Ġ = 0. (4.51)
As the function G is invertible,

∂f ∂G = -y ′ c (G) = -y ′ c (d(G)) = -y ′ c (d)d ′ (G) = -y ′ c (d) 1 G ′ (d) = - y ′ c (d) y c (d) . (4.52)
Note that we do not explicitly need G -1 . Combining these results, the discrete versions of evolution equations at step n read

∆G (n) i ≥ 0, (4.53)
f (n) i (µ, G i ) ≈ f (n-1) i + α (n-1) i ∆µ (n) + β (n-1) i ∆G (n) i ≤ 0, (4.54)
f (n) i ∆G (n) i = 0 (4.55) in which α (n-1) i = 2µy i ( d (n-1) i , t (n) i ) , β (n-1) i = - y ′ c ( d (n-1) i ) y c ( d (n-1) i ). (4.56)
Step numbers are in superscript and between parentheses, while location indices are in subscripts. In this section, the index i refers to a location on the set of all the interfaces where the Lagrange multipliers are defined (i.e. on the internal interfaces and on the Dirichlet boundaries). In other words, we do not make a distinction on which interface a location lies on. This allows the damage field to evolve automatically, no topological information is needed. The set of all locations is denoted by I. Some of these locations have completely failed by step n.

The not yet failed ones are indexed with J ⊆ I.

Computing ∆G (n) i
is performed similarly, but with relevant modifications, as it was done in [START_REF] Lé | Coupling damage and cohesive zone models with the Thick Level Set approach to fracture[END_REF] for ∆d (n) i . Rather than prescribing load increments, the quasi-static problem is solved stepwise by prescribing fractions of G c increments4 , and deducing the load level µ multiplying a given reference load F ref ,

F = µF ref (4.57)
such that equilibrium is obtained for the prescribed dissipation state. The step indices are omitted in the sequel for clarity, unless otherwise stated. From Eq. (4.41) we see that y ∝ ζ 2 and the linearity of Eq. (4.16) implies that ζ ∝ F ∝ µ, where the last proportionality comes from the definition of the load factor in Eq. (4.57). It shows that the energy release rate scales quadratically with the load factor. The energy release rate obtained for the reference load can then be scaled in such a way that the evolution criterion (4.50) is satisfied for all discrete Lagrange multipliers:

f i (µ, G i ) = µ 2 y i -y c (G i ) ≤ 0. (4.58)
Based on this information, the value of the load factor can then be identified as corresponding to the most critical location:

µ = min i∈J √ y c (G i ) y i . (4.59)
Here, we note that y c must be computed with the true fracture energy, not with the one which comes from the arbitrarily chosen reference load. This means that for the calculation of G c in (4.42), with G n and G t given in Tab. 4.2, the true normal traction component is taken into account to evaluate the mixed mode fracture energy. This is obtained by finding the intersection of the initial damage surface with the half-line stemming from the origin and given by the loading direction.

Although (4.58) is fulfilled with the chosen µ, it does not optimize the new damage front position, as one could wonder whether the damage should grow only at one position (i.e. the most critical one), or at other positions too. An iteration between the old and the new damage states could be used to answer this question [START_REF] Bernard | Damage growth modeling using the Thick Level Set (TLS) approach: Efficient discretization for quasi-static loadings[END_REF] however, a simpler scheme is used here. At this point, α i are known and β i are determined from the cohesive model. Instead of the inequality (4.54), a prediction is used based on

pred f i (µ, G i ) = f i + α i pred ∆µ + β i ∆G i = 0. (4.60)
To avoid the use of an iterative scheme, an explicit approach is used and we solve for pred ∆µ based on a chosen maximum energy dissipation increment ∆G max . First, the predicted load factor increment is expressed as

pred ∆µ = min i∈J -f i -β i ∆G max α i (4.61)
and the corresponding dissipation increments are computed using

∆G i = max { 0, -f i -α i pred ∆µ β i } . (4.62)
Note that Eq. (4.62) filters out the possibly negative increments, thereby obeying Eq. (4.53).

The remaining question is how to choose the step size ∆G max . Since the fracture energy, in general, is different from location to location, dissipating a fixed amount of energy would not be reasonable. Therefore, ∆G max is chosen to be a fraction of the fracture energy corresponding to the most critical location (i.e. where the load factor is determined for in Eq. (4.59)) at a given step:

∆G max = ∆G Gmax G c,ℓ , 0 < ∆G Gmax ≤ 1, ℓ is the most critical location (4.63)
We call the constant of proportionality ∆G Gmax maximum allowed relative energy dissipation increment and it is a parameter of the solution procedure. As opposed to ∆G Gmax , the step size ∆G max is not a global parameter. Generally, it changes from step to step because the most critical location also changes as the body continues to deform and degrade. Given the explicit nature of the approach, a jagged loaddisplacement response is expected, with larger oscillations for a larger ∆G Gmax . This choice for ∆G max may cause another location i to dissipate more than what is available:

G (n) i = G (n-1) i + ∆G (n) i > G c,i . In that case, we set G (n) i = G c,i
. This inaccuracy is due to the explicit nature of our method; decreasing ∆G Gmax in (4.63) results in a better resolution of the true dissipation.

One more ingredient is required for the algorithm to work properly. At any location, as soon as the damage initiates, the strength and the fracture energy are assumed to be constant for that location during the subsequent softening (cf. Section 4.5.2). One can imagine a situation in which later in the course of softening, due to stress redistribution, the local loading state at a location significantly changes. For instance, a location on an interface was under compression, when G c was fixed. Later in the simulation, if that location experiences tension for which the corresponding fracture energy is much smaller, fixing G c to the compressive value would result in a major inaccuracy. To avoid these cases, a tolerance is set such that energy starts being dissipated only if this dissipation is greater than that tolerance. Of course, the tolerance b must be relative, and is therefore set to be a fraction of ∆G Gmax :

b = 1 f 1 ∆G Gmax , f 1 > 0. (4.64)
By introducing the relative energy dissipation increment, defined as

∆G Gi = ∆G i G c,i , (4.65)
the dissipation increment predicted by (4.62) is not accepted if ∆G Gi < b, and G i = 0. (4.66)

As the weak form (4.16) is coupled with d and not with G, we determine ∆d

(n) i
from the following equation:

∆G (n) i = d (n-1) i +∆d (n) i ∫ d (n-1) i y c (d d) dd d. (4.67)
Solving this equation for ∆d

(n) i
in the tension+shear case results in a very long expression we obtained by version 16 of Maple. We only show the result when dissipation occurs only in one direction (normal or tangential). Introducing the dimensionless quantities

K n = k 2G n p 2 , K t = k 2G t q 2 , (4.68)
and K denoting one of them, the damage increment is given by ∆d

(n) i = ( K (n) i + d (n-1) i -K (n) i d (n-1) i ) 2 K (n) i ( K (n) i -1 + G c,i ∆G (n) i ) - ( K (n) i -1 ) 2 d (n-1) i . (4.69)
Finally, the damage field is updated to obtain the new damage state at every location: Compute y

d (n) i = d (n-1) i + ∆d (n) i . ( 4
(n) i , y (n) c,i , y ′ (n)
c,i from Tab. 4.2 13: end for 14: Find locations j ∈ J ⊆ I which are prone to further damage, i.e. d The positivity of y ′ c (d) is required on all parts of the failure criterion. The maximum values are to be found so that k is independent of the stress state. If k depended on the local stress values, a complicated and very costly algorithm would have to be used which results in a non-symmetric matrix and all global matrices would need to be assembled at each step. This would clearly destroy all the advantages of the presented method.

Store p i , q i , G n,i , G t,i , ∆G (n)
In the compression region of the failure criterion, along the Mohr-Coulomb line (label 3 in Fig. 4.2), only B t is active and is given by

B t (q) = q 2 2(G II + ap(q)) = q 2 2 ( G II + a q-c tan φ ). (4.72)
Asymptotically, B t = O(q) and so it is unbounded. That is why the compression cap was introduced to the model, which restricts the strength q to the interval [c, q 34 ].

It is emphasized that the use of a compression cap is not physically motivated and not even phenomenological, but rather numerical that allows choosing a finite value for the k parameter. Function B t takes its global maximum either at the end points (0, c) and (-p 34 , ±q 34 ) or where its derivative vanishes. Therefore,

B t,max = max { c 2 2G II , q 2 34 2(G II + ap 34 ) , 2 tan φ(ac -G II tan φ) a 2 } . (4.73)
It can be seen in Fig. 4.2 that

p 34 = ⃓ ⃓ ⃓ ⃓ c -f c tan ψ tan φ + tan ψ ⃓ ⃓ ⃓ ⃓ , q 34 = c + p 34 tan φ (4.74)
In the tension region of the failure criterion, the safe choice For optimal conditioning, k must be selected such that the conditioning of (4.26) is good, in which matrices L and C contain k. For the (1, 1) block be well-conditioned, entries of similar magnitude must be present in K and in L. In one-dimensional finite elements, k ≈ Eh is a good choice, where h is the characteristic mesh size.

max{B n,max , B t,max } = c 2 2G I (4.

Outline of the global solution scheme

Combining the discretization proposed in Chapter 3 with the present damage update procedure, a crack propagation problem for the heterogeneous grain assembly is solved according to Alg. 4.3. if n = 0 then 7:

Precompute and store C C (4.34)

8:
end if

9:

Create C (n) from C C using (4.33) 10:

Solve system (4.26) with C (n)

11:

Save variables of interest, if desired 12:

Update the damage field (Section 4.6, Alg. 4.2) 13: end for Finally, we note that when computing the energetic quantities, we do not need to know which interface a certain Lagrange multiplier corresponds to. In other words, the update procedure is automatic, no geometrical or topological information is required.

In this chapter different grain assemblies of increasing complexity are tested in configurations related to real tests. First, a classical three-point bending test is performed which validates the model for pure mode I. As the crack is expected to propagate vertically along the middle of the specimen, only one interface is considered there. The next example consists of a configuration inspired by a Brazilian test, widely used to obtain the tensile strength of quasi-brittle materials, and is used in the literature for rocks as well. The influence of the numerical and model parameters on the crack path obtained and on the load-displacement curve is thoroughly analysed in this example. To show the ability of the methodology to be used in large parametric studies, the effect of the heterogeneity is scrutinized. The last computational test is a uniaxial compression test which demonstrates that the model is valid in compression as well. Finally, based on these three examples, global conclusions are drawn. 101

Three-point bending test

The proposed framework is first validated on a quite simple configuration, the threepoint bending test, which is a classical mode I problem for the fracture of quasi-brittle materials. Although it consists of a single vertical interface, it incorporates two main components developed so-far: a non-matching mesh and the damage formulation with explicit solution procedure of the CZM. The three-point bending of an unnotched rectangular homogeneous concrete body is considered as sketched in Fig. 5.1. The body has length L = 1088 mm, height c = 500 mm, depth 40 mm and is supported in its bottom left and bottom right corners over a width v/2 = L/72. It is loaded in the middle of the top edge over a width v = L/36. The geometrical and material data of this test are taken from [START_REF] Hoover | Cohesive crack, size effect, crack band and work-of-fracture models compared to comprehensive concrete fracture tests[END_REF]. The concrete is considered to be isotropic and linear elastic with Young's modulus E = 37 GPa and Poisson's ratio ν = 0.21. The linear cohesive model depicted in Figs. 2.2b and 4.2 is not a suitable model for characterizing the softening of concrete. Using a bilinear softening law is appropriate to approximate the load-displacement response measured in laboratories [START_REF] Bažant | Concrete fracture models: testing and practice[END_REF]. The slope of the first part of the model, shown in Fig. 5.2, influences the peak load, while the second part characterizes the pull-out process of the grains [START_REF] Wittmann | Fracture energy and strain softening of concrete as determined by means of compact tension specimens[END_REF]. The cohesive properties can be found in Tab. 5.1. Two fracture energy quantities are defined in [START_REF] Bažant | Concrete fracture models: testing and practice[END_REF]: G F is the total fracture energy, i.e. G F = 1/2(f t u n,k + t n,k u n,c ), while G f = 1/2f t u n,1 corresponds to the area under the initial slope. Table 5.1: Cohesive properties of the three-point bending specimen from [START_REF] Hoover | Cohesive crack, size effect, crack band and work-of-fracture models compared to comprehensive concrete fracture tests[END_REF]. Notations can be followed in Fig. 5.2.

f t t n,k t n,2 u n,c u n t n u n,1 u n,k
f t (MPa) t n,k (MPa) u n,c (µm) u n,1 (µm) G f (N m -1 ) G F (N m -
The numerical parameter k, involved in the definition of the effective cohesive traction vector ζ, appears in the computation of the energy release rate and of the critical energy release rate (Tab. 4.2), as well as in the weak form (4.16). As the fracture energy is uniform in this test setup, the minimum value of k is determined as

k min = f 2 t 2G f ≈ 155 N mm -3 (5.1)
according to (4.71). In this problem k = 2000 N mm -3 was chosen, which is large enough to avoid any numerical artefacts. We observed that k could take values from a large interval without causing numerical difficulties. The simulation is performed in plane strain conditions with different maximum damage increments ∆d max . This is possible because the single interface is in mode I, and therefore a computation based on ∆G Gmax would lead to the same solution. A snapshot of the deformed state with a magnification factor of 100 is shown in Fig. 5.3. The zero vertical displacement BCs are taken into account by modifying the assembled matrix. For that, a conforming mesh is required in our current implementation, that is why the element edges almost coincide with the crack interface.

Brazilian test inspired simulations

The Brazilian test, also known as indirect tension test, is a widely used laboratory test to indirectly evaluate the tensile strength of quasi-brittle materials under the assumptions of uniform homogeneous material properties [START_REF] Li | The Brazilian Disc Test for Rock Mechanics Applications: Review and New Insights[END_REF]. An elastic simulation of it was used as a verification of our discretization scheme in Section 3.8, where a homogeneous body was considered. Assuming homogeneous material properties, the centre of the specimen is subjected to a mode I stress state in the horizontal direction if loaded by vertical concentrated forces. We determined the critical stress for the onset of damage in Section 3.8.

Single interface

After having computed the critical load for damage initiation in Section 3.8, we now follow the damage evolution. Similarly to the three-point bending test, the Brazilian test with a homogeneous bulk and a single interface is a mode I problem from the point of view of crack propagation. Therefore, the original damage update algorithm [START_REF] Lé | Coupling damage and cohesive zone models with the Thick Level Set approach to fracture[END_REF] is used here, prescribing the maximum allowed damage increment ∆d max at any step. The normal traction reaches its maximum in the centre of the interface, therefore the largest energy dissipation should occur there.

Grain assemblies

Using the Brazilian test for complex heterogeneous rocks is however made more difficult to interpret as the heterogeneity introduces deviations with respect to a pure mode I tensile stress state in the centre of the body, and the resulting cracking may be more jagged and distributed. In this section, we follow the softening response of the sample until complete failure.

From a computational point of view, the Brazilian test on grain-based materials is a good example to test the robustness of our discretization method and the damage propagation criterion. It involves interfaces that are slanted with respect to the main loading direction and that are therefore locally subjected to mixed mode loading. Moreover, the contact formulation is important as well in this second problem. A simplified 2D model of this test configuration is depicted in Fig. 5.6. As in the homogeneous case in Section 3.8, the circle is approximated as a regular polygon. The number of sides was set to 44 and the boundary conditions were explicitly given along the top and bottom interfaces. This choice allows a proper representation of the geometry, while at the same time not causing numerical problems in the presence of excessively concentrated forces at the top. The zero Dirichlet boundary condition is also prescribed along one segment at the bottom. We observed that a too small Dirichlet boundary fails to properly fix the body and the disk rotates as the load is applied. The application of the loading in the described manner corresponds to a central angle of about 8.18°, which is close to 10°, recommended by the standard [100]. We consider a disk of diameter D = 150 mm. The simulations are performed under plane strain conditions. All grain configurations in this chapter were produced by the Matlab function VoronoiLimit [START_REF] Sievers | VoronoiLimit[END_REF]. To make the example representative, samples with 500 grains were generated. Figure 5.7 shows that it provides more than ten grains along the diameter, as recommended by the standard [100]. The circle, or more precisely its polygonal approximation, was embedded into a background mesh. In the whole Section 5.2, the same background mesh with different element sizes will be used. The one with 50 × 50 elements is depicted in Fig. 5.7. Note that this overly coarse mesh is used only for the purpose of illustration, and not in actual computations. As opposed to the three-point bending test in which a local mesh refinement would have contributed to a major speed-up, this test problem would not benefit significantly from a locally refined mesh. We also note that the large number of elements in the background mesh which do not belong to any grain does not increase the computational time because degrees of freedom are associated only with the nodes of the cut elements. We are interested in the evolution of the microcracking patterns, in the final crack pattern and in the load-displacement response of the samples. Three types of plots will be used in the remaining of this chapter. For the load-displacement plots, the vertical displacement is monitored at the point (0.5D, 0.99D) in the coordinate system seen in Fig. 5.7. The second plot type will show the local relative energy dissipation G/G c (fraction of the fracture energy dissipated so far) on the interfaces, where G is the so-far dissipated energy given in Eq. (4.47) and G c is the mixed mode fracture energy expressed by Eq. (4.42). These plots will come in three forms. The coloured version maps data to colours with the parula colormap, see Fig. 5.8a. It is useful to display all the interfaces, therefore visualizing which possible crack paths are actually realized. The grayscaled version, shown in Fig. 5.8b, eases the comprehension of the progressive development of cracking but does not show the intact interfaces. The redscale colormap, see Fig. 5.8c, will be used to highlight the crack paths for elastically heterogeneous microstructures. Finally, the third plots will show the traction-free interfaces, i.e. those positions of the fully failed interfaces that have d = 1 at the end of the simulations. Such completely failed locations will be drawn in black according to the grayscale colormap.

Elastically homogeneous sample

The rock sample is first considered elastically homogeneous with Young's modulus E = 40 GPa and Poisson's ratio ν = 0.36, corresponding to the elastic properties of biotite [START_REF] Massart | Computational modelling of crackinduced permeability evolution in granite with dilatant cracks[END_REF]. In Section 5.2.3, the simulations will be performed on the grain configuration seen in Fig. 5.7.

A parametric study was first carried out to investigate the effect of some of the parameters of the model. In all cases, only one parameter was changed at a time. To make the numerical study consistent, a so-called master configuration was identified, with material parameter values inspired by the literature [START_REF] Vasconcelos | Experimental characterization of the tensile behaviour of granites[END_REF][START_REF] Vasconcelos | Ultrasonic evaluation of the physical and mechanical properties of granites[END_REF][START_REF] Vasconcelos | Compressive Behavior of Granite: Experimental Approach[END_REF], and with prescribed values of the solution procedure parameters. The values are contained in Tabs. 5.2 and 5.3, respectively. In Tab. 5.2, K K = k/k min and nElem denotes the number of elements in the background mesh. The cohesive properties of the interfaces are assumed to be all identical and are given in Tab. 5.3. The first five parameters characterize the Mohr-Coulomb criterion with tension cut-off and compression cap (cf. Fig. 4.1b), the next two parameters denote the pure mode I and mode II fractures energies, while the last one governs the growth of the fracture energy with increasing compressive stress. Although the pure mode I fracture energy is larger than the one measured in [START_REF] Vasconcelos | Experimental characterization of the tensile behaviour of granites[END_REF] for granite, the fracture energy given in our model is characteristic to the meso-scale, i.e the individual interfaces. Experiments, on the other hand, measure the fracture energy at the sample scale, which therefore corresponds to the macro-scale. Similarly to the fracture energies, the strengths f t , c and f c may differ from what is measured in laboratories. If quantitative results are desired, the meso-scale parameters must be obtained by calibration.

K K (-) ∆G Gmax (-) nElem (-)
f t (MPa) c (MPa) f c (MPa) φ (°) ψ (°) G I (N mm -1 ) G II (N mm -1 )

Effect of the material parameters

The Brazilian test is expected to characterize a globally mode I failure. Therefore, we were first interested in the effect of the interfacial tensile strength. As we observed only minor differences in the crack path, only the global load-displacement response is shown in Fig. 5.9. It shows that the peak load increases with increasing f t , implying that interfaces indeed fail predominantly in mode I. One can also see that the peak loads are proportional to the prescribed interfacial f t values, suggesting that using the macroscopic tensile strength in the microscopic model is probably acceptable.

The snapback is automatically handled by our damage framework because it works by imposing damage increments and finding the corresponding load that satisfies the equilibrium equations at that step. The jaggedness of the response curve is due to the solution procedure. In classical load-controlled schemes, the load-displacement response is smooth because the load is prescribed and the degradation is the end result of the equilibrium state. In our case, finite energy dissipation steps are chosen, for which the necessary load is deduced. A smooth response therefore cannot be expected. Both the load-controlled and the damage-controlled solution procedures have their advantages and drawbacks. ). This can be explained as follows. A small a parameter only marginally increases the mixed mode fracture energy due to compression, therefore interfaces which happen to be under compressive loading locally are prone to fail. On the contrary, large values for a increase the fracture energy in compression so much that energy dissipation becomes possible only in tension. Since the central vertical interfaces experience the highest tensile stress levels, those interfaces become favourable for cracking. Close to the application of the load, the cracking is more distributed in case of small a, which is logical because failure in shear is more favourable for small values of a than for large values. These experiments suggest that the major crack path remains similar for a wide range of a. For K K = 50, the simulation was stopped because due to the large k, as observed for the three-point bending as well, only a few locations dissipated energy at a step, leading to very slow propagation. The crack path in Fig. 5.12d and Fig. 5.13d is, of course, unrealistic and substantially differs from the ones obtained with smaller K K. Our recommendation is therefore to use K K close to 1. The effect of the step size ∆G Gmax on the crack path can be seen in Fig. 5.14 (see also 5.15 for the fully open cracks). The resulting fracture patterns are similar, the difference being particularly small for the two smallest step sizes. We are also interested in the load-displacement curves, which are plotted in the same figure 5.16. Similar conclusions can be drawn as for the crack path: the largest step size is only suitable for capturing the qualitative response, a quantitatively matching solution requires smaller steps, as a too large step size underestimates the peak load. Next, the effect of the mesh density on the accuracy of the solution is studied. The background meshes consist of 100, 150, 250 and 350 number of elements in both directions. Figure 5.17 shows the failure patterns on these different meshes. Even the coarsest mesh predicts the main crack path reasonably well. Finally, the progressive energy dissipation evolution with the smallest step size ∆G Gmax = 0.001 can be followed in Fig. 5.18. One can see in Fig. 5.18a that interfacial failure is preceded by diffuse microcracking along the whole vertical diameter of the sample. Some interfaces experience higher tensile stresses than others as can be observed in Fig. 5.18b. These interfaces are typically parallel to the loading direction and are located near the place where the external load is applied. Figure 5.18c suggests that the macro-crack propagates towards the centre. Once the centre is reached, the crack continues downwards as illustrated in Fig. 5.18d. As the bottom Dirichlet BC hinders the complete vertical separation of the disk, a secondary crack appears on the boundary and merges with the existing crack. Based on the parametric studies above, we can state that crack evolution is more complex to interpret than it was in the fully homogeneous, single interface case in Section 5.2.1.

Heterogeneity

We now investigate how the effect of the heterogeneity influences the preferential crack paths. For the sake of illustration, this example will be simulated with properties of the phases similar to the Lac du Bonnet granite the material properties of which are taken from reference [START_REF] Massart | Computational modelling of crackinduced permeability evolution in granite with dilatant cracks[END_REF]. The three main mineral species of this type of granite are biotite, quartz and feldspar. The assumed Young moduli and Poisson's ratios are E biotite = 40 GPa, E quartz = 100 GPa, E feldspar = 80 GPa, ν biotite = 0.36, ν quartz = 0.07 and ν feldspar = 0.32. Volume fractions of 0.07, 0.2, 0.73 are assumed for each phase. Ten grain realizations were generated, two of which are shown in Fig. 5.20. The constituents were associated with grains such that the same minerals form clusters. For the remainder of this chapter, the colour codes in Fig. 5.19 will represent the three constituents of the Lac du Bonnet granite. With the elastic heterogeneity itself, no significant difference in the solution was found with respect to the elastically homogeneous sample in Section 5.2.3. Therefore, higher tensile strengths and fracture energy values were set for those interfaces that are located between two quartz grains. The combined effect of these two kinds of heterogeneities is expected to influence the crack path. Indeed, cracks will develop round the hard phase, i.e. the quartz grains. In the simulations the tensile strength and the cohesion of these strengthened interfaces were set to 1.25 times the values given in the master configuration in Tab. 5.3, while the mode I and mode II fracture energies were multiplied by 10. 5.20b, the following observations can be made. As expected, the main crack bypasses the hard phase and tends to go along the boundary of the quartz clusters. Two crack propagation mechanisms can be identified for the two grain assemblies. In the first realization in Fig. 5.20a, the quartz phases form a large cluster just below the centre of the disk. This causes the crack to stop (see Fig. 5.21b). Instead of breaking through the quartz group, it is more favourable to have new crack initiations at the bottom, which is shown in Fig. 5.21c. Then this new crack propagates upwards, finally leading to the failure of the quartz cluster (5.21d). In the second realization in Fig. 5.20b, the small quartz cluster at the top and the large quartz group at the bottom stand in the way of a vertical cracking, observed in homogeneous configurations. As Fig. 5.23a shows, the crack does not initiate at the top as in the elastically homogeneous case (cf. Fig. 5.18b), but slightly below. As opposed to the previous configuration when the main crack stopped when reaching the large quartz island, it now takes a sharp turn at about 45°to avoid the hard phase. Also note that compared to the elastically homogeneous example, the unloading is more pronounced when heterogeneity is included in the model. Due to the heterogeneity, we expect a larger influence of the parameter a, governing the mixed mode compressive fracture energy, than what was experienced for the elastically homogeneous disk. In both configurations, the crack paths in Fig. 5.25 demonstrate that when a is small, the cracks can easily bypass the quartz grains. On the other hand, for large values of a, the crack tends to remain vertical. This is logical because with increasing a the fracture energy increases in compression and shear, making crack propagation possible only in tension. As the central vertical interfaces are dominantly in tension, the cracking algorithm chooses those interfaces. This is only qualitative reasoning because the strong heterogeneity changes the stress state and the central vertical interfaces may not remain is tension. To show the versatility of our crack propagation framework, the remaining eight grain realizations are briefly discussed. Figure 5.26 shows the final crack paths for the master parameters. We can observe that the cracks tend to propagate in the softer minerals and that the quartz cluster boundaries "attract" the cracks. In places where large quartz clusters stand in the way of the propagating macrocrack, the macrocrack either bypasses or if it cannot, new macrocracks initiate.

Uniaxial compression

A uniaxial compression test is now modelled in this section. The size of the rectangular specimen is 108 mm×54 mm. To avoid localized failure in compression and shear near the top and bottom edges, platens are attached, modelled as rigid grains (i.e. no relative displacement is allowed between the specimen and the platen). The interfaces of the grains touching the platens are assumed to be free of damage. A distributed force is applied along the top edge, and the bottom edge is fixed. We note that such boundary conditions do not impose locally a uniaxial stress state, but they are meant to represent at best the real boundary conditions in a real test. Nevertheless, the configuration is suitable to test whether our model works in compression-dominated problems. The elastic, cohesive and numerical parameters are the same as used in the Brazilian test in the previous section. This time 1000 grains were generated for the specimens. The mechanical model can be seen in Fig. 5.27a, while the distribution of the constituents for the heterogeneous version of the test is shown in Fig. 5.27b. The load-displacement response is reported in terms of distributed force on the top edge versus vertical displacement of the loading platen. All the uniaxial compression test simulations were performed on a 100×200 mesh and under plane strain assumptions. , where the assumption is made that the lateral principal strains are equal to the unique lateral strain computed in a plane strain simulation. In a perfectly uniaxial stress state, ε 1 and ε 2 could be determined by keeping track of the axial and lateral displacements, respectively. Due to our boundary conditions that do not allow rollers on the top and bottom boundaries, the stress state is neither uniaxial nor uniform, and the body deforms with barreling. With these boundary conditions, the elastic volume change is expected to be smaller than in the case when the lateral displacement on the horizontal boundaries is made possible (e.g. with Teflon sheets). More than that, the current boundary conditions do not allow parallel vertical crack openings, hence the volume of the opening cracks is also smaller than for the "proper" boundary conditions.

Elastically homogeneous sample

In the uniaxial compression test, the whole body is loaded with a dominant compressive loading due to the sample geometry. This is in contrast to the Brazilian test in which the body is locally subjected to mode I loading at the centre of the specimen, which triggered mode I cracking. On the other hand, the uniaxial stress state is not ensured, so the cracking is assumed to be less distributed than what a pure mode loading would give. The loading and the boundary conditions suggest that even in the case when the whole domain consists of only one constituent, the (compressive) material parameters have a more significant effect on the cracking pattern than what we experienced in the Brazilian test. Since the various material and solution procedure parameters were thoroughly tested for the Brazilian test, here we restrict ourselves to the study the effect of those parameters that are presumed to play a significant role.

Different combinations of the tensile strength f t and the cohesion c were triedthe crack paths being reported in Fig. 5.28. In the first setup, f t was set to 3 MPa and c was given the values 26, 40 and 45 MPa. These values are within the interval of the measured (macroscopic) data [START_REF] Massart | Computational modelling of crackinduced permeability evolution in granite with dilatant cracks[END_REF]. The major crack initiates either from the bottom left or from the bottom right corner of the specimen, which is probably due to the prescribed zero horizontal displacement on the bottom edge. The figures show that the orientation of the crack paths slightly changed with increasing c, but not substantially between f t = 3 MPa and f t = 6 MPa with a given c. As c increases, compressive-shearing failure of slanted interfaces is less likely, and the cracks tend to grow on interfaces under tension. This explains why the crack path deflects towards the vertical loading direction. Since the platen is modelled as a very stiff grain, the crack could not break it, resulting in the onset of an auxiliary crack, growing parallel to the platen. The load-displacement response for the investigated combinations of f t and c are shown in Fig. 5.29.

The peak stress for the investigated parameter values is provided in Tab. 5.4. They indicate that the macroscopic compressive strength is proportional to the interfacial tensile strength and cohesion values, fed to the microstructure. Moreover, these computationally determined compressive strengths are in the range 30 MPa to 160 MPa, measured for the Lac du Bonnet granite [START_REF] Vasconcelos | Compressive Behavior of Granite: Experimental Approach[END_REF].

We can conclude that the uniaxial test, or at least the way it has been modelled in this section, makes the results difficult to interpret, let alone predicting the crack path before executing the simulations.

Heterogeneity

As in the case of the Brazilian test, we also carry out simulations on heterogeneous samples. First, the microstructure in Fig. 5.27b is considered. Several computations were done with different parameters; the one with f t = 6 MPa and c = 40 MPa gave an interesting result. To test another realization, the one in Fig. 5.32 has been generated. With the same parameter values as before, the response of the body to the loading is visualized in Fig. 5.33. The large snapback between "times" A and B corresponds to the abrupt crack arrest in Fig. 5.34b. The complete diagonal crack could not develop because of the hard phase. Hence, cracks grow in the top left corner (Fig. 5.34c). This crack growth becomes more intensive in Fig. 5.34d, and finally the quartz cluster at the bottom right breaks. Similarly to the heterogeneous Brazilian test, the reinforced interfaces are also not prone to break as it can be followed in Fig. 5.35 for the two investigated grain realizations. This implies different damage initiation points, depending mainly on the distribution of the phases but also on the grain configuration.

Concluding remarks

Based on the computational experiments for the three test problems, the following consequences can be drawn.

• Accurate crack paths are obtained on relatively coarse meshes.

• A large step size was sufficient for all the test cases to predict the main crack path and to obtain a qualitatively correct load-displacement curve. Decreasing the step size proportionally increased the computational time, but resulted in a more accurate and smoother load-displacement response.

• Although k is a numerical parameter, it influences which locations dissipate energy at a given step. The experiments show that this parameter should be close to k min (i.e. K K ≈ 1) to get correct results.

• The elastic heterogeneity alone had a minor effect, therefore we tried different cohesive properties on some selected interfaces. These two types of heterogeneities produced physically meaningful crack paths, indicating that a real microstructure also owes distinct local fracture energies, which are complicated to measure.

• The boundary conditions play a fundamental role. As a crack approaches a Dirichlet boundary, it would require more and more force to propagate further. Therefore, other cracks tend to initiate at other parts of the body. The inhomogeneous BCs, especially if they act on small surfaces, cause stress concentration, which requires many elements to be evaluated accurately.

• The oscillations have two causes. One is physical and is attributed to the stress redistribution in the body (occurs e.g. in case of heterogeneities). The other kind of oscillation has a numerical origin and is particularly relevant in configurations with more than one interface. It is caused by the fact that the most critical location "jumps" from one place to another between two steps.

The proposed computational tool is useful to carry out numerical investigations on microstructures in an automated way, without any intervention from the user. One can freely vary the distribution of the mineral species, the grain size distribution, etc. The explicit damage stepping can handle the complete softening behaviour of the brittle Brazilian test, without any complicated path-following techniques. This page is intentionally left blank

Summary and global comments

The principal aim of the thesis was to enable efficient and flexible microcrack simulations in grain-based materials. As a first step towards this goal, simplifying assumptions were made on the possible crack paths and on the constitutive model of the grains. Based on these assumptions, a two-field variational formulation was developed to represent the displacement field on each grain, and the traction field on each interface and Dirichlet boundary segment. The Lagrange multipliers provided the (weak) continuity of the displacement field on the whole assembly of grains. This construction allows independent local calculations, prone to parallelization. The discretization of this mixed formulation is achieved with CutFEM, which makes the meshing trivial. All the unknowns are defined at the nodes of the rectangular elements. The displacement space is spanned by the bilinear shape functions of rectangular elements. The Lagrange multiplier shape functions are constructed as a restriction of the nodal shape functions on the interfaces, with some connections among them. Establishing the reduction algorithm for the Lagrange multipliers on Q1 elements is a novelty of our work. Interfaces are allowed to cross the elements arbitrarily. Numerical evidence was given of the stability of the resulting discrete mixed method.

The gradual degradation due to microcracking was then modelled by an extrinsic cohesive law, derived from a potential. The potential was defined as a function of the cohesive opening, from which the traction-separation relation can be obtained, and an internal variable, which plays a role in formulating the evolution equations. To avoid the problem of infinite terms in the compliance matrix upon the failure of a cohesive zone, an augmented Lagrangian method was used, modifying the original pure Lagrangian formulation. Self-contact between the (partially) open crack lips was enforced explicitly. It turned out that by discretizing this augmented Lagrangian method, most of the matrices can be precomputed, therefore the additional matrix assembly at each step is very cheap. As the cohesive parameters are not uniform along the interfaces, the same damage increment at two different locations is not associated with the same amount of degradation in general. To take this into account, a physically meaningful quantity, the energy dissipation, was introduced to drive the energy release. The damage increment, present in the weak form, was then subsequently computed from the energy dissipation increment. Algorithms shed light on the details of the overall solution process.

To harness the capabilities of the developed framework, three examples were considered. The three-point bending test verified the numerical algorithms for mode I cracking. A thorough analysis of the Brazilian test revealed the influence of the material and solution procedure parameters. The number of grains in the Brazilian disk allowed a meaningful representation of realistic microstructures. Our method enables the user to easily perform statistical analysis by generating random grain packings and running computations on them. The third test was the uniaxial compression test, where even more grains were taken. Qualitatively correct results were obtained both for the Brazilian test and for the uniaxial compression test, using crude meshes and large step sizes.

We consider the discretization method and the energy-based solution procedure for the propagation of multiple cracks as the main strength of the thesis. The discretization itself, without including the damage, is worth highlighting due to its attractive properties: (i) no meshing difficulties, (ii) accurate computation of the traction vector, (iii) block diagonal subdomain stiffness matrices K i , and (iv) known null-space of K i , ∀i ∈ I s . The mesh-independence has the additional merit that the same background mesh can be used to run simulations on different grain distributions, which becomes useful for RVE computations or for parametric investigations as illustrated in Chapter 5. The quasi-static stepping based on the energy dissipation has the advantage over the load-controlled solution procedures that a complex load-displacement response can be tracked without sophisticated path-following techniques. It further increases the automated nature of the framework since the user does not have to experiment with the proper choice of the load-controlling parameters so that the nonlinear iteration procedure does not fail. Its combination with an extrinsic cohesive zone formulation, parametrized by the energy dissipation, results in the automatic tracking of microcrack initiation and propagation. Another nice property is that the damage update is separated from the solving of the mechanical system, making a segregated solution procedure possible. Although the explicit update algorithm produces jagged load-displacement curves, this is mostly aesthetic as an envelope of the graph is easily recognizable.

The developed framework in its current stage has some deficiencies. Although attempts were made in Section 3.7 to achieve matrices with smaller condition number, a more efficient strategy is needed to be able to use iterative solvers. Some prospects are recommended in item N.1 of the next section. The other main problem is that with the introduction of the augmented Lagrangian method, the penalty term destroys some good characteristics of the original saddle point system: (i) the choice of the penalty-like parameter k is not dictated by the conditioning only, but its minimum value is determined by the cohesive material constants and the failure criterion; and (ii) the subdomain stiffness matrices K i are complemented by interfacial terms, increasing the bandwidth and making the preconditioning more difficult.

Future work

Now, that the initial ideas were realized in a mature computational framework, several further improvements are worth investigating in the future. These ideas are classified as follows.

P. Additional physical phenomena

P.1 Frictional contact. Friction between surfaces of fractured rock may have a significant effect and could be included in the framework.

P.2

Coupling with fluid flow in the microcracks for permeability computations.

In the introductory section 1.1, we saw that microcracking-induced permeability change is important in practice. We plan to couple our damage model with low-velocity fluid flow in the channels formed by open cracks.

Existing attempts, such as the ones in [START_REF] Massart | Computational modelling of crackinduced permeability evolution in granite with dilatant cracks[END_REF] and [START_REF] Van Den Eijnden | A FE 2 modelling approach to hydromechanical coupling in cracking-induced localization problems[END_REF], can benefit from the non-conforming discretization.

P.3

Intrinsic cohesive zone model. Rock grains investigated in this thesis are closely packed, therefore it was reasonable to apply extrinsic cohesive zone.

Our discretization was successfully tried on rubble masonry structures with non-Voronoi grains, but the resulting crack path was not satisfactory. This can be explained by the finite width and compliance of the mortar joints around the stones, requiring the use of an intrinsic cohesive zone model.

P.4 Different angles of friction among the phases.

For the Brazilian test and the compression test, the interfaces within a quartz cluster were set to higher strength and fracture energy values than for the other interfaces. It is worth checking what happens if the angle of friction is also changed on the interfaces. The different angles of friction on the interfaces might be obtained experimentally with a direct shear test machine.

P.5 Determine the dilatancy.

Due to its practical relevance, the characterization of dilatancy is on high priority on our list of future investigations. We noted in the discription of the uniaxial compression test in Section 5.3 that computing the dilatancy requires "roller-type" boundary conditions. See also N.2.

M. Modelling questions

M.1 3D modelling. Three spatial dimensions allow richer kinematics, essential to capture the permeability increase of fractured geomaterials before the occurrence of the full loss of the load-bearing capacity. The construction of the discrete Lagrange multiplier space remains the same in 3D, i.e. the cut edges play a role and not the intersected faces. The failure criterion and the damage update scheme also remain almost identical. However, it will be more difficult to cut the possibly concave polyhedral grains with the elements. The development to perform is therefore essentially of geometrical nature.

M.2 Initial microcracks. Microcracks in rocks are results of either natural processes (tectonic movements, etc.) or man-made activities (excavation, boreholes, etc.) [START_REF] Hamdi | Characterizing the influence of stressinduced microcracks on the laboratory strength and fracture development in brittle rocks using a finite-discrete element method-micro discrete fracture network FDEM-µDFN approach[END_REF]. They act as stress-concentrators, significantly lowering the strength of rocks. Higher microcrack density causes the increase of permeability. The incorporation of pre-existing microcracks into our model can be easily done by prescribing non-zero damage values for some locations on some interfaces. Completely open cohesive cracks can be easily incorporated by setting d = 1 at the desired locations (e.g. along the interfaces lying between feldspar grains). This is an advantage of our dissipation-driven solution method over the classical load-controlled schemes in which it is not straighforward how to place the initial microcracks.

M.3 Effect of grain concavity. Our method works without further development for concave grains too. It would be interesting to characterize the convexity of grains and investigate how macroscopic dilatancy is influenced by the uplifting effect of the shearing of convex-concave grain connectivities.

M.4 Quantification of the microstructure. Use statistical methods, such as correlation functions, to quantify the clustering of the mineral species.

M.5 Homogenization. Computational homogenization is typically an application our tool would excel at. With a single background mesh, a set of grain realizations for statistical studies can be generated to obtain macroscopic responses.

M.6

Reduced-order modelling. The calculations could be made much fasterthough sacrificing accuracy -by decreasing the number of unknowns in the bulk and on the interfaces. Fewer degrees of freedom for the bulk can be achieved by allowing fewer deformation modes for the elastic grains. The number of interfacial unknowns (Lagrange multipliers) might be lowered to three. These three scalar Lagrange multipliers would weakly enforce the three interface opening modes: opening in the normal direction, opening in the tangential direction and rotating around a point on the interface. Finally, one could try reducing all the unknowns onto the interfaces by static condensation. Reworking the current discretization framework to be compatible with model reduction is probably not easy and therefore it has low priority.

N. Numerical tasks

N. [START_REF] Abdelaziz | Grain based modelling of rocks using the combined finite-discrete element method[END_REF] Preconditioning. An alternative to the attempts we tried to decrease the condition number is adding extra terms to the weak form which efficiently improves the conditioning, such as the ghost penalty method [START_REF] Burman | Ghost penalty[END_REF]. Generalpurpose (black-box) preconditioners are often not powerful enough for saddle point matrices [START_REF] Benzi | Numerical solution of saddle point problems[END_REF]. With a deeper understanding of the continuous operators, one increases the chance to construct effective preconditioners [START_REF] Hiptmair | Operator Preconditioning[END_REF].

N.2 Support other types of boundary conditions. The use of periodic boundary conditions in computation homogenization allows one generating a smaller RVE compared to the case when displacement or traction boundary conditions are prescribed. The literature already contains solutions to weakly enforce periodicity [START_REF] Svenning | Computational homogenization of microfractured continua using weakly periodic boundary conditions[END_REF]. Periodicity conditions can be implemented similarly to jump conditions, which are already contained in our method.

Another improvement is to be able to prescribe only one component of the displacement vector. This corresponds to the mechanical model of roller support, used for example in the three-point bending test and the uniaxial compression test. A recent work discussing it in the non-conforming case is [START_REF] Lu | Weak impositions of Dirichlet boundary conditions in solid mechanics: A critique of current approaches and extension to partially prescribed boundaries[END_REF].

N.3 Higher-order Lagrange elements. In [START_REF] Ferté | Interface problems with quadratic X-FEM: design of a stable multiplier space and error analysis[END_REF] a stable Lagrange multiplier space was designed for conventional P2 triangles in X-FEM. One could try using higher-order elements with all degrees of freedom associated with the nodes, as it was done in GFEM [START_REF] Duarte | Generalized finite element methods for three-dimensional structural mechanics problems[END_REF]. A nice generalization of the current discretization would be a Q k -Q k ⋆ interpolation, i.e. k-th order polynomials for the displacements and reduced k-th order polynomials for the Lagrange multipliers -all degrees of freedom placed at the nodes.

I. Implementation

I.1 Distributed computing capabilities. Most operations are local, therefore little communication would be necessary among the distributed computing nodes.

I.2

Local (non-conforming) mesh refinement. Currently, only a uniform mesh is supported. In certain examples (e.g. the three-point bending test), a local refinement would significantly reduce the memory and CPU resource requirements. The refinement could either be linked to an a posteriori error estimate or simply be purely geometric. For smooth boundary data (e.g. L 2 -regularity), we expect a smooth displacement solution inside the grains and jumps across the interfaces. Therefore, refining the mesh only in the neighbourhood of the interfaces seems reasonable. It has the advantage that the mesh can be fixed in the beginning, so the matrix sizes do not change and therefore the matrices can be precomputed, as it is done in the current implementation. Note that a non-conforming quadtree-based refinement is sufficient in our framework.

I.3

Performant, non-proprietary software. Using Matlab was convenient during the prototyping phase, but for reproducibility reason and to promote open science, free tools are welcome. The C++ language is free, compiles to fast executables, and is widely used in scientific computing with a large number of existing libraries.

I.4

Be able to track the loading mode. To obtain more insight on why a given crack path develops, it is helpful if we know which part of the failure criterion is reached at a given location.

Some of these issues are part of ongoing work.

B Quadrature

The use of formulation (4.16) requires computing domain (area) integrals on cut or uncut elements and interface (line) integrals. The stiffness matrix on uncut elements can be exactly integrated by the 2 × 2 Gauss-Legendre quadrature. A vast literature is available about quadrature on complicated domains. The most relevant ones were discussed in Section 2.5.4. The most straightforward method is to decompose the cut element into simpler shapes for which the integration is easier. A constrained Delaunay triangulation, visualized in Fig. 3.10, is used here for this purpose, together with a three-point Gauss-Legendre quadrature on each triangle. The reason we applied subtriangulation is due to its simplicity and accuracy. Low order elements (Q1 in our case) require few quadrature points, especially in 2D where the number of subtriangles is low. Hence, not much would be gained by using advanced techniques. The triangulation can be carried out either on the physical domain or on the reference domain. If the interfaces bounding the cut element part are mapped to the reference square, they become curved due to the bilinear mapping. In [START_REF] Sabada | A stable X-FEM in cohesive transition from closed to open crack[END_REF], an improved integration is used to avoid this problem for bilinear quadrangular elements. Therefore, the cut elements are partitioned into triangles on the physical domain, these triangles being then mapped to a reference domain (see Fig. B.1). Although the inverse bilinear mapping χ -1 2 exists in closed form [START_REF] Hua | An inverse transformation for quadrilateral isoparametric elements: Analysis and application[END_REF], the affine mapping between the reference square and the physical element allows a simpler inversion for the Cartesian mesh we use in this work. The integrand in the interface integrals is the product of the shape functions, restricted on the interface. When a bilinear shape function is evaluated on a straight segment, a complete linear or an incomplete quadratic (with a missing linear term) polynomial is obtained. The integrand is therefore described by four parameters. A two-point 1D Gauss-Legendre quadrature is sufficient to calculate the integral exactly. On the other hand, the full integration of the Lagrange multipliers can produce oscillations, as shown in [START_REF] Schellekens | On the numerical integration of interface elements[END_REF]. As noted in [START_REF] Nguyen | An open source program to generate zero-thickness cohesive interface elements[END_REF], the fictitious node method has the same interface matrices as that of interface elements, which motivates using the twopoint Gauss-Lobatto integration here as well. Note that the two-point Gauss-Lobatto integration points in 1D are located at the interval endpoints. Hence, a quadratic polynomial without the linear term cannot be integrated exactly. This means that the line integrals are not computed exactly when the interface goes through the neighbouring edges of a quadrilateral element. For the precomputed matrix (4.34), the product of the three shape functions requires three Gauss-Legendre quadrature points for exact integration. Title : Mesh-independent modelling of diffuse cracking in cohesive grain-based materials Keywords : mesh-independent, damage, cohesive zone, CutFEM, Lagrange multiplier, cracking Abstract : In this thesis a flexible and general stable displacement-Lagrange multiplier mixed formulation is developed to model distributed cracking in cohesive grain-based materials in the framework of the cut finite element method. The displacement field is discretized on each grain separately, and the continuity of the displacement and traction fields across the interfaces between grains is enforced by Lagrange multipliers. The design of the discrete Lagrange multiplier space is detailed for bilinear quadrangular elements with the potential presence of multiple interfaces/discontinuities within an element. We give numerical evidence that the designed Lagrange multiplier space is stable and provide examples demonstrating the robustness of the method. Relying on the stable discretization, a cohesive zone formulation equipped with a damage constitutive model expressed in terms of the traction is used to model the propagation of multiple cracks at the interfaces between grains. To prevent the crack faces from self-penetrating during unloading, a contact condition is enforced. The solutions for the mechanical fields and the damage field are separately obtained and an explicit damage update algorithm allows using a non-iterative approach. The damage formulation couples the normal and tangential failure modes, accounts for different tension and compression behaviours and takes into account a compression-dependent fracture energy in mixed mode. The framework is applied to complex 2D problems inspired by indirect tension tests and compression tests on heterogeneous rock-like materials.

Figure 2 . 1 :

 21 Figure 2.1: Cohesive zone

Figure 2 . 2 :

 22 Figure 2.2: Two types of cohesive zone models

Figure 2 . 3 :

 23 Figure 2.3: Conforming and non-conforming meshes for a square domain with a hole

Figure 2 . 6 :

 26 Figure 2.6: Patches and some of the basis functions in the phantom node and floating node methods

Figure 2 . 7 : 2 Figure 2 . 8 :

 27228 Figure 2.7: An interface passing very close to a node

Figure 3 . 1 :

 31 Figure 3.1: Notations used for the geometrical description of the domain of interest

. 15 )

 15 Equations (3.4)-(3.5) are equivalent to Equations (3.8)-(3.15).

Figure 3 . 2 :

 32 Figure 3.2: Notations for the local meshes

  vectorial degrees of freedom. By substituting (3.36) in(3.19a)-(3.19b), the discrete saddle-point system reads:

Figure 3 . 3 :Figure 3 . 4 :

 3334 Figure 3.3: Notions used for the stable Lagrange multiplier space

Figure 3 . 5 :

 35 Figure 3.5: Naive strategy to equate Lagrange multipliers. The node links are denoted by an elliptical shape.

Figure 3 . 6 :Figure 3 . 7 :

 3637 Figure 3.6: Constructing the Lagrange multiplier basis functions for a single interface

Figure 3 . 8 :

 38 Figure 3.8: Geometry and mesh for the simple tension test. The interface endpoints are marked with filled blue circles.

Figure 3 . 9 :Figure 3 . 10 :

 39310 Figure 3.9: Displacement field for the tension test

Figure 3 . 11 :

 311 Figure 3.11: Traction field for the tension test

Fig. 3 .

 3 13a). One way of connecting the intersection points to form new interfaces is Γ (a) First strategy Γ (b) Second strategy

Figure 3 . 13 :

 313 Figure 3.13: The original interface (dashed line) almost goes through a node. Two strategies to move an intersection point away from a node.

Figure 3 . 14 :

 314 Figure 3.14: The original interface (dashed line) exactly goes through a node. Two strategies to move an intersection point away from a node.

Figure 3 . 15 :

 315 Figure3.15: 2D geometry-mesh configurations leading to small cuts. Zoom on a triple junction with three subdomains. The subdomain filled in with green tilted lines cuts a small part out of the element marked in yellow.

Figure 3 . 16 :

 316 Figure 3.16: Initial and optimal background meshes for a single interface

Figure 3 . 17 :

 317 Figure 3.17: Background mesh for multiple interfaces

Figure 3 . 18 :

 318 Figure 3.18: Brazilian test

Figure 3 . 19 :

 319 Figure 3.19: Normal traction distribution in the Brazilian test

Page 4 . 1

 41 Crack propagation along interfaces . . . . . . . . . . . . . . . . 79 4.2 Mixed method for cohesive cracks . . . . . . . . . . . . . . . . . 80 4.3 Failure criterion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82 4.4 Damage-based cohesive zone formulation . . . . . . . . . . . . 85 4.5 Computing the energetic quantities . . . . . . . . . . . . . . . . 90 4.5.1 Energy release rate . . . . . . . . . . . . . . . . . . . . . . . . 90 4.5.2 Decomposing the mixed mode fracture energy . . . . . . . . . 91 4.5.3 Critical energy release rate . . . . . . . . . . . . . . . . . . . . 91 4.6 Damage update algorithm . . . . . . . . . . . . . . . . . . . . . . 93 4.7 Choosing the parameter k . . . . . . . . . . . . . . . . . . . . . . 98 4.8 Outline of the global solution scheme . . . . . . . . . . . . . . . 99

Figure 4 . 1 :

 41 Figure 4.1: Mohr-Coulomb damage criterion with tension cut-off and compression cap in the space of interfacial traction components

Figure 4 .

 4 2 shows the Mohr-Coulomb criterion for positive shear stresses with cut-off in the tension and cap in the compression regime. One can differentiate four loading modes, depending on which part of the criterion becomes active when the softening begins. Label 1 marks the tension cut-off part, label 2 the part of the Mohr-Coulomb line being in the tension region, label 3 the Mohr-Coulomb line in the compressive regime, while label 4 indicates the compression cap. The corresponding traction-separation relations are also indicated in the figure.

  c

Figure 4 . 2 :

 42 Figure 4.2: Failure criterion

  5 (cf. Eq. (3.42), (3.43), (3.45)). The matrix L has a block structure

34 )

 34 can be precomputed and stored. Equation(4.34) shows that computing C requires more quadrature points than what is needed for the computation of B in (3.43) and P in (4.23). Algorithm 4.1 gives the procedure to precompute the global matrices that do not change with propagating damage (• denotes the Hadamard product).

Algorithm 4 . 1

 41 Precomputing global matrices and vectors 1: ▷ Preprocessing 2: Input: grain configuration, elastic properties of grains, cohesive properties of interfaces, boundary conditions 3: Generate a uniform Cartesian mesh covering all the grains 4: Adapt the geometry (Section 3.7.1) or optimize the mesh (Section 3.7.2) 5: ▷ Construct the subdomain stiffness matrices 6: for i ∈ I s do 7:

  i g : reduction operator to get (3.45) 25:

15 : 17 :

 1517 Compute the load factor µ (n) from (4.59) 16: Compute f (i) i = µ (n) 2 y Calculate the maximum energy dissipation ∆G (n) max from (4.63) 18: Compute the load factor increment pred ∆µ (n) from (4.61) with j ∈ J 19: Compute the energy dissipation increments ∆G

Algorithm 4 . 3 3 : 4 :

 4334 Quasi-static simulation 1: Inputs: parameters required by the algorithms called by this algorithm 2: Determine k min by (4.76) in Section 4.7 Get global matrices K, P, B and vectors f , g by Alg. 4.1 Initialize the damage field 5: for n from 0 to N step do ▷ damage stepping loop 6:

Figure 5 . 1 :

 51 Figure 5.1: Sketch of the three-point bending problem

Figure 5 . 2 :

 52 Figure 5.2: Bilinear softening function for pure mode I loading

Figure 5 . 3 :Figure 5 . 4 :

 5354 Figure 5.3: Crack opening in the three-point bending test (magnification: 100)

Figure 5 . 5

 55 shows the damage evolution for different step sizes. The damage initiates at the centre of the disk and symmetrically propagates outwards -captured even by the crudest steps size. The results were obtained for the configuration shown in Fig.3.18a, with a mesh of 80×80 elements. Plane strain conditions were assumed and the material properties were set to E = 40 GPa, ν = 0.36, f t = 6 MPa, G I = 0.37 N mm -1 . ∆d max = 0.001

Figure 5 . 5 :

 55 Figure 5.5: Damage front evolution for the Brazilian test in Fig. 3.18a for various damage step sizes. The distance along the vertical diameter of the disk (in mm) is measured on the vertical axis.

F 8 Figure 5 . 6 :

 856 Figure 5.6: Sketch of the Brazilian test problem

Figure 5 . 7 :

 57 Figure 5.7: Brazilian disk embedded to a background mesh

Figure 5 . 8 :

 58 Figure 5.8: Colormaps to represent the relative energy dissipation on the interfaces

Figure 5 . 9 :

 59 Figure 5.9: Load-displacement response for various tensile strengths f t

1 Figure 5 . 10 :

 1510 Figure 5.10: Relative dissipated energy on the interfaces for different parameters a

(a) a = 1 × 1 Figure 5 . 11 : 5 . 2 . 3 . 2 Figure 5 .

 1151152325 Figure 5.11: Traction-free cracks for different parameters a

  (a) K K = 1.01 (b) K K = 5 (c) K K = 10 (d) K K = 50

Figure 5 . 12 :

 512 Figure 5.12: Relative dissipated energy on the interfaces for different parameters K K

Figure 5 . 13 :

 513 Figure 5.13: Traction-free cracks for different parameters K K

3 Figure 5 . 14 : 3 Figure 5 . 15 :Figure 5 . 16 :

 35143515516 Figure 5.14: Relative dissipated energy on the interfaces for different parameters ∆G Gmax

Figure 5 . 17 :

 517 Figure 5.17: Relative dissipated energy on the interfaces for different mesh sizes

Figure 5 . 18 :

 518 Figure 5.18: Progressive crack propagation in the configuration corresponding to Fig. 5.14c

Figure 5 . 19 :

 519 Figure 5.19: Colours representing the mineral species in the Lac du Bonnet granite

2 Figure 5 . 20 :

 2520 Figure 5.20: Distribution of the mineral species in the Lac du Bonnet granite used for the Brazilian test. The interfaces among the quartz grains, marked in white, are made stronger.

Figure 5 . 21 :Figure 5 . 22 :

 521522 Figure 5.21: Progressive crack propagation in the configuration corresponding to Fig. 5.20a

Figure 5 . 23 :Figure 5 . 24 :

 523524 Figure 5.23: Progressive crack propagation in the configuration corresponding to Fig. 5.20b

Figure 5 . 25 :

 525 Figure 5.25: Dependence of the crack path on the parameter a

Figure 5 . 26 :

 526 Figure 5.26: Failure patterns in further grain realizations for the Brazilian disk

F

  Sketch of the uniaxial compression test (b) Distribution of mineral species inspired by the Lac du Bonnet granite used for the uniaxial compression test. The interfaces among the quartz grains, marked in white, are made stronger.

Figure 5 . 27 :

 527 Figure 5.27: Uniaxial compression test

  (a) f t = 3 MPa, c = 26 MPa (b) f t = 3 MPa, c = 40 MPa (c) f t = 3 MPa, c = 45 MPa (d) f t = 6 MPa, c = 26 MPa (e) f t = 6 MPa, c = 40 MPa (f) f t = 6 MPa, c = 45 MPa

Figure 5 . 28 :

 528 Figure 5.28: Relative dissipated energy on the interfaces for different parameters f t and c

Figure 5 .

 5 Figure 5.29: Load-displacement response for the investigated elastically homogeneous uniaxial compression test

Figure 5 .

 5 30 shows the microcracking. The damaging process starts with microcracks initiation all over the domain. A completely open crack appears near the centre. (Fig. 5.30a). The crack then tends towards the bottom right corner (Fig. 5.30b). Concurrently, a second main crack appears near the left edge of the rectangle, probably due to the quartz cluster nearby. According to Figs. 5.30c and 5.30d, crack branching near the left side is more favourable than further crack growth along the diagonal.

Figure 5 . 30 :

 530 Figure 5.30: Progressive crack propagation in the configuration corresponding to Fig. 5.27b

Figure 5 . 31 :

 531 Figure 5.31: Relative dissipated energy on the interfaces for different parameters f c for the configuration in Fig. 5.27b

Figure 5 . 32 :

 532 Figure 5.32: Another grain realization for the uniaxial compression test

Figure 5 . 33 :Figure 5 . 34 :

 533534 Figure 5.33: Load-displacement response for the configuration in Fig. 5.20b

Figure 5 . 35 :

 535 Figure 5.35: Failure patterns due to heterogeneity in the uniaxial compression test

Figure B. 1 :

 1 Figure B.1: Improved integration on polygons

Titre:

  Modélisation de la fissuration dans les matériaux cohésifs à grains par discrétisations non conformes Mots clés : maillage-indépendante, endommagement, zone cohésive, CutFEM, fissuration, multiplicateur de Lagrange Résumé : Dans cette thèse, une formulation stable mixte de déplacement-multiplicateur de Lagrange est développée pour modéliser la fissuration dans les matériaux cohésifs à grains dans le cadre de la méthode des éléments finis etendus (CutFEM). Le champ de déplacement est discrétisé sur chaque grain individuellement, et la continuité des champs de déplacement et de traction aux interfaces entre grains est assurée par des multiplicateurs de Lagrange. La construction de l'espace discret des multiplicateurs de Lagrange est détaillée pour les éléments quadrangulaires bilinéaires avec la présence d'interfaces multiples dans un élément. Des preuves numériques sont données que cet espace de multiplicateurs de Lagrange est stable, et des exemples démontrant la robustesse de la méthode sont fournis. Avec cette discrétisation stable, une formulation de zone cohésive permet de modéliser la propagation de fissures multiples aux interfaces entre grains. Pour éviter des interpénétrations aux faces des fissures pendant le déchargement, une condition de contact est imposée. Les solutions pour les champs mécaniques et le champ d'endommagement sont obtenues séparément et un algorithme explicite permet d'utiliser une approche non itérative. La formulation de l'endommagement associe les modes de rupture normal et tangentiel, tient compte de différents comportements de tension et de compression et prend en compte une énergie de rupture dépendante de la compression en mode mixte. La méthode est appliquée à des problèmes 2D complexes inspirés par des tests de tension indirecte et des tests de compression sur des matériaux hétérogènes ressemblant à de la roche.

  2.1 classifies the reviewed discretization techniques for the model problem (2.10).

	formulation	equation	conformity consistent mixed reference
	unconstrained	(2.12)	H 1 0	yes	no	[53] 1
	penalty	(2.14)	H 1	no	no	[8]
	Nitsche	(2.25)	H 1	yes	no	[138]
	Lagrange multiplier	(2.18)	H 1 , H -1/2	yes	yes	[9]
	augmented Lag. mult. (2.28)	H 1 , H -1/2	yes	yes	[82]
	Barbosa-Hughes	(2.27)

67 3.8 Brazilian test configuration

  . . . . . . . . . . . . . . . . . . . . . 74

3.1 Definitions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 3.2 Governing equations . . . . . . . . . . . . . . . . . . . . . . . . . 46 3.3 Mixed continuous formulation . . . . . . . . . . . . . . . . . . . 48 3.4 Approximations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 3.5 Discrete weak form . . . . . . . . . . . . . . . . . . . . . . . . . . 52 3.6 Stable mixed formulation . . . . . . . . . . . . . . . . . . . . . . 54 3.6.1 Construction of the Lagrange multiplier space . . . . . . . . . 54 3.6.2 Stability of the discretization scheme . . . . . . . . . . . . . . 59 3.7 Robustness . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66 3.7.1 Interface adaptation . . . . . . . . . . . . . . . . . . . . . . . 66 3.7.2 Optimal background mesh generation . . . . . . . . . . . . . .

  Note that two subdomains can have multiple common interfaces. The interfaces are indexed by I int = {1, . . . , N }. The neighbouring subdomains to Γ i are denoted by Ω i+ and Ω i-, or Ω i± when dealing with them as a collection. Parts of the boundary of Ω where non-zero Neumann boundary conditions (BC) are prescribed are denoted by Γ N , whereas Γ D denotes parts where zero or non-zero Dirichlet BCs are given. It is assumed that Γ N ∩ Γ D = ∅. These boundary parts are replaced by straight segments

Table 4 .

 4 1: Normal and tangential component of the mixed mode fracture energy depending on the loading mode

  , q i , G n,i , G t,i already fixed Use p i , q i , G n,i , G t,i Determine p i , q i from tDetermine G n,i and G t,i from Tab. 4.2

	Algorithm 4.2 Damage update
	1: Inputs: ζ	(n) i , d (n-1) i	∀i ∈ I, numerical parameters: ∆G Gmax (4.63), f 1 (4.64)
	2: for i ∈ I do 3: Determine ζ n,i and ζ (n) t,i from ζ (n)	(n) i	using (4.30)
	4: 5:	Calculate the damage indicators χ Determine the cohesive tractions t (n) (n) i n,i and t according to (4.14) (n) t,i from (4.15)
	6: ▷ p i 7: if G n,i exists then
	8: 9:	else		(n) n,i , t (n) t,i and Fig. 4.2
	10:
	11:	end if
	12:
				.70)
		Algorithm 4.2 summarizes the update of the damage field.

  As mentioned in Section 4.5, the critical energy release rate must be an increasing function of the damage. Taking the derivative of y c (d) in Eq. (4.45), one can easily derive the sufficient condition

		4.7 Choosing the parameter k
			k > k min = max	⎧ max{B n,max , B t,max } tension and shear ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ B t,max compression and shear ⎪ f 2 t /(2G I ) pure tension c 2 /(2G II ) pure shear	(4.71)
				i
	26: 27: 28:	end if Compute ∆d Determine d (n) (n) i i	from (4.67) from (4.70)
	29: end for 30: return d	(n) i	∀i ∈ I

Table 5 .

 5 2: Master configuration for the Brazilian test

	1.01	0.1	150 × 150

Table 5 .

 5 3: Cohesive properties of the Brazilian test specimen

	a (mm)

Table 5 .

 5 4: Peak stress values obtained from the uniaxial compression test (all units are in MPa)

This page is intentionally left blank

The mesh-conformity is not to be confused with a conforming discretization. It is possible to define a non-conforming discretization on a conforming mesh, as it is done in the discontinuous Galerkin method.

Although the widely accepted definition of the finite element by Ciarlet[START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF] is general, under standard FEM, or sometimes just FEM, we mean mesh-conforming finite element discretization based on low order Lagrangian finite elements in this thesis.

the acronym of Ladyzhenskaya-Babuška-Brezzi

We will show soon how damage increments can be obtained from dissipated energy increments.

Acknowledgements

This page is intentionally left blank Bibliography

To my parents

Contents

Appendices

A Physical meaning of the Lagrange multiplier

Integration by parts on Eq. (3.19a) yields

where S m is the boundary of Ω m . From the above domain integral we gain the equilibrium equations (first equation of Eq. (3.8)). Each subdomain boundary S i can be decomposed into a disjoint set of segments, where one of interface, Neumann or Dirichlet conditions are prescribed. Regarding the Neumann BCs, from Eq. (A.1) we get

which yields Eq. (3.10) for arbitrary v i+ . For the Dirichlet BCs, we obtain

The remaining terms in Eq. (A.1) represent the interface integrals:

as stated in (3.24). Since n i+ = -n i-, (A.6) provides the traction continuity across the interfaces, Eq. (3.15). Finally, (3.19b) reads

which directly gives the displacement continuity across the interfaces (cf. Eq. (3.14)) and the Dirichlet boundary condition (cf. Eq. (3.9)).