Véronique Lazarus

Omar Grazie

Keywords: fracture, variational methods, laminate composites, matrix-cracking, transverse cracks, damage models, finite elements, phase-field models. 7 rupture, approches variationnelles, composites stratifiés, fissuration matricielle, fissures transverses, élements finis, méthode de champ de phase composites: fibre-breakage . . . . . . .

I want to seize this opportunity to thank the people who've made an important impact by encouraging and supporting me in this journey.

First and foremost, I want to thank my advisors, Corrado Maurini and Julien Berthe, for believing in me and this work. Your constant support played a key role in this journey and was fundamental for my growth both as a scholar and a person. It was really a pleasure to work with you, both on a human and professional level. This journey has been challenging at times, but definitely worth it and you

List of Figures

List of Tables

4.

1 Key properties of the gradient damage models considered in this work as a function of the material parameters E 0 , w 1 and η appearing in (4.11): the elastic limit σ e , the material strength σ c , the fracture toughness G c , and the width of the damage band D loc .

Thesis outline

The use of composite materials has skyrocketed in the past two decades thanks to their small weight and high stiffness. A few popular fields of application are: civil engineering, aviation, space, sports, bio-medicine and so on. We study in this thesis the fracture of cross-ply unidirectional (UD) laminates subjected to uni-axial tension, a classical topic in the field. We tackle the problem within the variational approach to fracture initiated by [START_REF] Griffith | The phenomenon of rupture and flow in solids[END_REF] and revisited in the framework of the mathematical theory of the calculus of variations by [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF] and Bourdin et al. [2000b]. We focus on the problem of the determination of the crack spacings in the networks of transverse cracks appearing in the core layer of symmetric laminates.

The content of this thesis is organised in the four chapters outlined below.

Chapter 1

State of the art of matrix cracking in unidirectional (UD) laminated composites: in this chapter we introduce the basic notions of a composite material, following the typical damage mechanisms occurring in laminated composites and finally a state of the art on matrix-cracking in UD cross-ply fibre-reinforced composite laminates from an experimental viewpoint mainly. A brief introduction on numerical approaches, explaining the main differences between cohesive-zone models, continuum damage models and variational approaches to fracture allowing the reader to understand the interest we take in applying a variational approach.

Chapter 2

Experimental procedure and results on carbon fiber reinforced UD cross-ply laminates through passive infrared thermography: a first analysis of the physical cracking patterns: in this chapter we introduce the experimental results of a campaign conducted prior to this thesis. A new experimental protocol, relying on infrared thermography, is presented as well as the results that will be taken into account for comparison in the following chapters.

Chapter 3

Analytical and numerical approaches implemented to study transverse cracking patterns in UD cross-ply laminates: in this chapter we introduce the analytical model we developed to study transverse cracks.

We introduce a sharp interface model à la Griffith focused on a cell-based approach. We discuss the results obtained for the crack spacings, the competition between failure modes with uniform and non-uniform spacings, and the influence of the material and the geometric parameters on the resulting crack pattern.

Chapter 4

Phase-field approach: application of a local minimisation approach to a composite bar: in this chapter we introduce in depth the phase field approach and show the classical example of a one-dimensional bar under tension modelled through gradient damage model. We pursue the analysis on the composite bar, discussing the multiple issues related to the problem. Finally, we discuss our numerical results in the light of the sharp interface model.

1

State of the art of matrix cracking in unidirectional (UD) laminated composites

The purpose of this chapter is to introduce laminated composites with their failure mechanisms and to contextualize the research carried out in this thesis among the existing works present in the state of the art.

Introduction and general overview of composite materials

Introduction and general overview of composite materials

This section aims at giving a few key and general notions on composite materials to help better introduce the topic of this thesis.

What are composite materials ?

By general definition a composite material is the composition of two or more materials exhibiting different material and mechanical properties, which, taken as a unit reveal higher material and mechanical properties than if they were considered separately [START_REF] Gay | Matériaux Composites -6ème édition révisée[END_REF], [START_REF] Berthelot | Matériaux Composites-Comportement Mécanique et Analyse des Structures[END_REF]). The type of interaction between these two (or multiple) materials in the most generic case can be imagined as a discontinuous phase (or multiple discontinuous phases) distributed in a continuous phase. The continuous phase is called the matrix, while the discontinuous phase is known as the reinforcement due to its classically higher mechanical properties. Depending on the nature of the matrix one can classify them into three main categories:

• Organic matrix (epoxy resins, polyurethane resins, polyamide resins, phenolic resins.),

• Metallic matrix (aluminium alloy, magnesium alloy, titanium alloy)

• Mineral matrix (ceramics)

The discontinuous phase can either be a fibre reinforcement or a particle reinforcement. The same classification, as for the matrices, can be applied to the fibres resulting into multiple different fibrematrix combinations. From a macroscopic point of view, the nature-describing architecture can be classified in three categories:

• Monolithic structures • Laminated structures • Sandwich structures However, reinforced plastics and volume composites also exist. These are only few generalities, the reader will find a more exhaustive and detailed introduction to composite materials in [START_REF] Berthelot | Matériaux Composites-Comportement Mécanique et Analyse des Structures[END_REF] and [START_REF] Gay | Matériaux Composites -6ème édition révisée[END_REF]. In this thesis we consider and discuss high performance composite materials for aeronautical approaches, i.e. carbon-fibre reinforced epoxy-resin laminated structures.

Unidirectional plies and laminates

A laminate is essentially the result of a stacking sequence of layers called laminae. These laminae consist of fibres embedded in a matrix, as shown in Figure 1.1c.

Different kinds of resins constituting the matrix exist and their role is to transfer mechanical stresses to the fibres shielding them from the external environment. In order to optimise their purpose, the resin-fibre combinations have to be compatible. The resins classically used in our context can be separated in two main categories:

• Thermo-setting resin: polyester resin, vinyl ester resin, epoxy resin A taxonomical classification can be seen in Figure 1.2 which sum up all the different discussed categories of composite materials based on reinforcements and matrices. In this context we will focus on epoxy resin composites mainly, which are highlighted in red in the charts.

Epoxy resins are the most popular after polyester resins thanks to their high mechanical performances, in fact matrices with epoxy resins are classified as high performing. However, due to their high cost, they only represent a small percentage of the market. As listed in [START_REF] Berthelot | Matériaux Composites-Comportement Mécanique et Analyse des Structures[END_REF], some of the advantages are: good mechanical performances in tensile, compressive, flexion, shock and creep tests; good temperature resistance (up to 190 • C); excellent chemical resistance; low moulding shrinkage (between 0, 5 -1%); excellent adhesion to metallic materials. Meanwhile some of the main inconveniences are: long curing time; high cost; precautions for use and sensitivity to cracking.

The embedded fibers can be oriented in different directions. When they follow the same direction within the ply they are called unidirectional (UD). The layers can be oriented in different directions as well, creating different stacking sequences. A key concept in composites is that the material anisotropy can be managed through the orientation of the fibres.

Main stacking sequences

Laminates are generally designated by three criteria: the number of layers, the angles of each layer with respect to a plane of reference and finally any symmetry related to the stratification. The most studied stacking sequences, to better understand the behaviour of the laminate subjected to specific loadings, are: UD-ply laminates, angle-ply laminates, cross-ply laminates and quasi-isotropic laminates.

UD-plies

UD-ply laminates are widely used for characterisation of mechanical properties in order to obtain for example the Young Modulus in the different material directions.

Angle-plies

When considering angle plies, each ply is typically constituted of the same material and has the same thickness. Their main characteristic is that the sequence order alternates the fiber's orientation of ±θ with respect to the plane of reference.

Cross-plies

Cross plies are as well typically constituted of the same material and each ply has the same thickness. The particular feature to them is that the sequence uses only 0 • and 90 • oriented fibers.

Quasi-isotropic laminates

Quasi-isotropic laminates correspond to stacking configurations of the type: [-45 m /0 m /45 p /90 q ] s . These configurations meet the various design rules of the aeronautical industry (angular distance between the different plies, presence of fibres in different directions, etc.) and are therefore widely used for the production of laminated structures.

Main damage and failure mechanisms in cross-ply laminates

Damage in composite materials occurs at different length-scales from microscopic fibre-matrix debonding to macroscopic matrix-cracking [START_REF] Ogin | 1 -composite materials: constituents, architecture, and generic damage[END_REF]). In this context we will focus on the macroscopic damage modelling viewpoint of laminated composites. Depending on the loading tests, laminates generally exhibit three main failure or damage mechanisms. These mechanisms are: intralaminar failure, interlaminar failure and finally translaminar failure. Each failure mechanism will be

Intralaminar failure in laminated composites: matrix-cracking

Intralaminar damage forms in the matrix and at the interface between the fibers and the matrix. It is one of the most common damage mechanisms to occur in laminated composites and tends to appear mostly under tensile loadings [START_REF] Okabe | Numerical modeling of progressive damage in fiber reinforced plastic cross-ply laminates[END_REF], [START_REF] Maimí | Numerical modeling of matrix cracking and intralaminar failure in advanced composite materials[END_REF]). This damage mechanism observable at the ply scale is due to microcracks that coalesce and grow (fiber-matrix debonding) to finally form a macroscopic transverse or matrix cracks. This denomination is due to its characteristic shape and orientation. It is typically considered to be the first damage mechanism at the ply scale to occur during low-velocity impacts [START_REF] Vinogradov | Probabilistic energy based model for prediction of transverse cracking in cross-ply laminates[END_REF], [START_REF] Trousset | Prévision des dommages d'impact basse vitesse et basse énergie dans les composites à matrice organique stratifiés[END_REF], [START_REF] Zhang | Predicting low-velocity impact damage in composites by a quasistatic load model with cohesive interface elements[END_REF], [START_REF] Choi | A model for predicting damage in graphite/epoxy laminated composites resulting from low-velocity point impact[END_REF]). These cracks generally propagate through the thickness of the ply and not through the thickness of the laminate (unless induced by fatigue [START_REF] Ogin | 1 -composite materials: constituents, architecture, and generic damage[END_REF]) as the outer layers are usually tougher and tend to contain the damage or at least prevent the fibers from breaking up to a certain extent. These cracks are usually non-catastrophic as the laminate can continue to endure larger loads long before catastrophic fibre breakage, however, they can induce micro-damage mechanisms and delamination between the plies. Delaminations are associated to a loss of flexural and compressive strength which can provoke dangerous failure mechanisms. To sum up, here below are listed some of the structural damages caused by matrix-cracking reported from [START_REF] Manders | Statistical analysis of multiple fracture in [0/90/0] glass fiber/epoxy resin laminates[END_REF]:

• reduction of the effective stiffness of the 90 • ply,

• allowing of fluid weepage,

• allowing moisture and corrosive fluids to penetrate,

• may lead to delamination between plies,

• reduction of the effective structure's toughness.

Interlaminar failure in laminated composites: delamination

Interlaminar failure is most commonly referred to as delamination, which is essentially the separation between plies in the laminate. This type of damage is considered to be the most critical since it brings the structure to catastrophic failure if the delamination propagates. These cracks originate 1.3. State of the art on prediction of matrix-cracking patterns in UD cross-ply laminates 29 at the interface between plies that have different orientations since there is a difference of stiffness [START_REF] Choi | A new approach toward understanding damage mechanisms and mechanics of laminated composites due to low-velocity impact: Part i-experiments[END_REF]) show that initiation and propagation of interlaminar fractures occurs only in the presence of matrix-cracking. After propagation starts, delaminations are most likely to continue along the interface since the encountered fibres in the transverse plies have greater toughness they will not allow cracks to propagate perpendicularly. However, delaminations can propagate in two different ways which are typically called cohesive delamination and adhesive delamination. The first one consists in the propagation of the crack within the matrix between plies, whereas the latter consists in the propagation of the crack between the fibers and the matrix. The propagation usually includes both ways in alternation, creating irregular profiles. The stacking sequence is one of the first key elements influencing the shape of the delaminated area.

Translaminar failure in laminated composites: fibre-breakage

Translaminar failure is usually the last damage mechanism to occur after matrix-cracking and delamination. It is typically associated with tensile or compressive tests leading to failure of the fibres in the laminate. The fracture toughness depends on three different failure mechanisms, as stated in [START_REF] Oterkus | 4 -fracture modes, damage tolerance and failure mitigation in marine composites[END_REF], fibre/matrix debonding, fibre failure and fibre pullout, however it usually has larger magnitude values which is why is often bypassed in modeling. This type of damage leads to a total rupture of the structure and therefore is the most serious type of damage. Finally, we can see in Figure 1.4 a schematic view of three main damage mechanisms in a [0 • /90

State of the art on prediction of matrix-cracking patterns in UD crossply laminates

In this section a short overview of matrix-cracking in composites, mainly fibre-reinforced cross-ply laminates, is presented. This section focuses on regrouping, following a chronological order, the main studies that have been published in this context to predict and study transverse cracks in UD laminates and how these works motivated and fit in the line of work of this thesis.

Chapter 1. State of the art of matrix cracking in unidirectional (UD) laminated composites

Overview of matrix-cracking in composite fibre reinforced cross-ply laminates

Matrix-cracking in fibre-reinforced laminates has been investigated for many years from the late '70s to today, experimentally, analytically and numerically. Starting with Parvizi and Bailey [1978], Parvizi et al. [1978] who conducted a series of experiments on [0 • n /90 • m ] sym glass reinforced polyester resin cross-ply laminates under uniaxial loading in the direction parallel to the fibres. They investigated the effect of matrix cracks on the specimens by measuring the transverse tensile strength, knowing that in general no delamination was observed for relatively thin plies (with respect to the specimen length), however for thicker specimens, where the number of 90 • layers were increased, both small debondings and oblique cracks were observed, see Figures 1.5-1.6. [1978] Their results show a remarkable even distribution of the cracks in the 90 • layers and an influence on crack spacing dependent on the thickness of the 90 • layers. In fact, transverse crack spacing was found decreasing while applied stress was being increased and ply thickness decreased. However, crack spacing values were underestimated for given applied stresses, which was assessed in Parvizi and Bailey [1978], Parvizi et al. [1978], who reported more detailed studies on multiple transverse cracks in unidirectional glass fibre reinforced epoxy cross-ply composites aiming at explain the discrepancy in the previous works. Transverse crack spacings were measured through the use of a travelling microscope. Both studies used a modified shear lag analysis which is a technique first suggested by Parvizi and Bailey [1978] and is essentially based on the transfer of longitudinal stress between the 0 • and 90 • layers. More specifically this technique considers the axial load to be transferred by shear of an interlaminar 'resin-rich' transverse layer as emphasized in [START_REF] Manders | Statistical analysis of multiple fracture in [0/90/0] glass fiber/epoxy resin laminates[END_REF], [START_REF] Silberschmidt | Matrix cracking in cross-ply laminates: Effect of randomness[END_REF]. From these studies a dependence on the thickness of the damageable layers in the laminate was finally proved. In fact Parvizi et al. [1978], Parvizi and Bailey [1978] show that the critical (nucleation) strain at which multiple transverse cracking begins is higher for thinner 90 • plies thanks to an energetic principle defined in [START_REF] Aveston | Proceedings of national physical laboratory conference "the properties of fibre composites[END_REF] stating that the energy required by the crack surface decreases linearly with the ply thickness, while the energy release decreases more rapidly. In Figure 1.7, it can be seen that the strain initiating micro-cracking increases as the thickness of the 90 • layer decreases. The data shows that for 90 • plies with thickness between 0.1mm and 0.4mm partial micro-cracks initiate at the free-edges and then propagates through the width, see Parvizi et al. [1978]. The influence of ply thickness on matrix-cracking is also discussed in [START_REF] Flaggs | Experimental determination of the in situ transverse lamina strength in graphite/epoxy laminates[END_REF], where results show that the transverse strain failure increases as the ply thickness decreases. [START_REF] Manders | Statistical analysis of multiple fracture in [0/90/0] glass fiber/epoxy resin laminates[END_REF] investigated as well on glass fibre epoxy resin cross-ply laminates under tensile loading in the direction of the fibres and measured crack positions by travelling microscope. Their study shows how the 90 • ply has a variable strength, this implies that the critical strain at which transverse cracks appear and the crack spacings also depend on the specimen's length. They observe that crack distribution is not consistent with the assumption of Parvizi et al. [1978], Parvizi and Bailey [1978] of uniform stress in the 90 • ply. This is one of the first studies to point out the importance of the crack spacings. A statistical approach based on the Weibull distribution of the stress in the 90 • layer is applied. Although [START_REF] Manders | Statistical analysis of multiple fracture in [0/90/0] glass fiber/epoxy resin laminates[END_REF]'s model correct the Garret-Bailey-Parvizi approach, which is based on the assumption of uniform strength in the 90 • plies implying that new cracks will FIGURE 1.7: "The strain to initiate micro-cracking in E glass/Shell Epikote epoxy [0/90/0] laminates as a function of the total thickness of the 90 • plies. The 0 • plies each have a constant thickness of 9.5mm". Figure and description from [START_REF] Nairn | Matrix microcracking in composites[END_REF] whom reproduced the data from Parvizi et al. [1978].

always form midway between two cracks with the largest separating space, they still assume that the strengths of all the constituent volumes are independently and identically distributed. This method only works in certain conditions, for instance only for widely spaced cracks. In fact, "When the crack spacing is significantly larger than the unstressed length [...] the probabilistic model is appropriate, and when it is of similar magnitude the Garret-Bailey-Parvizi model is more appropriate" [START_REF] Manders | Statistical analysis of multiple fracture in [0/90/0] glass fiber/epoxy resin laminates[END_REF]).

Most experimental quasi-static tensile tests on cross-ply laminates show that transverse cracks appear instantaneously from edge to edge [START_REF] Berthelot | Transverse cracking of cross-ply laminates: Part 1[END_REF]). All of the studies mentioned above have based their approach on evaluating stress distributions in order to study the progression of transverse cracking. [START_REF] Berthelot | Transverse cracking of cross-ply laminates: Part 1[END_REF] studied as well the progression of transverse cracking through a stress based approach, more specifically through a parabolic shear-lag analysis which led to a good approximation of the strain-stress distributions that they've obtained through FEM (Finite Element Analysis). Their experimental and analytical results were in good agreement with the ones by Garret-Bailey-Parvizi.

In general, crack spacing is assumed to be mostly uniform and the focus, over the years, was mainly shifted on stress distributions and crack density rates. Different cracking patterns have been investigated in [START_REF] Silberschmidt | Matrix cracking in cross-ply laminates: Effect of randomness[END_REF], where these cracks are described to be "randomly" distributed. They investigated, a very similar case to the one that will be treated in this work, crack spacings on carbon-fibre epoxy resin [0 n /90 m /0 n ] cross-ply laminates subjected to uniaxial tension. The article emphasizes that, although a variety of analytical methods to describe the experimental process of matrix cracking exist, most of them share the assumption of uniform longitudinal distribution of transverse cracks and that the homogenized schemes used for simplification are reduced to a unit cell (representative unit volume) approach. They argue that reducing the laminate to a cell based approach implies basing the work on the concept of uniformity, or as stated "periodicity" of crack distribution. The concept of periodicity is mentioned in their work as the phenomenon where all cracks are evenly distributed. In this context, two main loading cases were treated, uniaxial (tensile) loading tests and fatigue (cycles). [START_REF] Silberschmidt | Matrix cracking in cross-ply laminates: Effect of randomness[END_REF] tested cross-ply carbon fibre-epoxy T300/914C laminates, and emphasizes four different cases of crack distribution, three of which are periodic with different crack densities and one described to have a "random" distribution, hence aperiodic which can be seen in Figure 1.8. X-ray radiography results were approximated by a three-parameter Weibull approximation. Despite the use of probability density and cumulative distribution functions, the number of statistical distributions of cracks is infinite and the values of the three parameters of the distribution are obtained through analysis of experimental data.

One of the main key arguments tackled in their work, and all of the works previously mentioned, are based on the description analysis of the influence that crack spacing has on stress distribution in the 0 • and 90 • layers. The reason for the interest can be easily explained: when considering uniaxial longitudinal loading where ply thickness is considered constant, the localisation of cracks that manifest in matrix-cracks are linked to axial variations on the material properties. Axial variations can either be introduced through flaws, which are managed through statistics for micro-cracks, or Chapter 1. State of the art of matrix cracking in unidirectional (UD) laminated composites through strength distributions (e.g. shear lag).

In conclusion, numerous studies have dealt with microcracking coalescing into transverse cracks in cross-ply laminates. In the present literature to date, little if any studies have shown an energetic or variational approach allowing to study and predict the cracking pattern or regime of crack propagation marking the passage between periodic cracking and aperiodic cracking. We have based and developed our approach on an interesting work to this approach presented in [START_REF] Bourdin | The Variational Approach to Fracture[END_REF]. Furthermore, we've seen that cell-based approaches are considered to be unfit to predict aperiodic cracking patterns. In this thesis, we aim at studying, predicting and analyzing through different variational approaches, based on a one-cell approach, without the use of any statistical distribution of the stress, the different cracking patterns that can be obtained in both isotropic and orthotropic settings for different material parameters.

A short introduction to numerical approaches for the simulation of damage mechanisms at different length scales

From a numerical viewpoint we can divide in three main categories numerical methods. On the one hand we have discrete models, or sharp models that seek to describe cracks as surface discontinuities (eg. Cohesive Zone Models). On the other hand, we have continuum damage models. These models seek to describe damage through continuum mechanics and the homogenization of the micro structure properties. In the middle, smeared approaches have been developed (e.g. variational approaches). These describe damage as a scalar variable assuming either a zero value or unitary, depending on the damage state. In this approach sharp discontinuities are smeared across bands of non-vanishing width.

Damage mechanisms can intervene at different length scales, as the material itself can be considered as multi-scale [START_REF] Ladevèze | Multiscale modelling of damage and fracture processes in composite materials, chap. multiscale computational damage modelling of laminated composites[END_REF]). For this reason different damage describing models for different scales have been developed. These can be divided into three main categories:

• Micro models: the common dimensions are approximately around ≈ 10µm, which correspond to the fibre scale. At this scale three main damage mechanisms can be observe.

1. Diffused damage in the ply: this is manifested as matrix-fibre decohesion [START_REF] Li | Modeling interfacial debonding and matrix cracking in fiber reinforced composites by the extended voronoi cell fem[END_REF], [START_REF] Segurado | A computational micromechanics study of the effect of interface decohesion on the mechanical behavior of composites[END_REF], [START_REF] Caporale | Micromechanical analysis of interfacial debonding in unidirectional fiber-reinforced composites[END_REF]). Also the appearance of cracks in the matrix (micro-cracking, [START_REF] Nairn | A critical evaluation of theories for predicting microcracking in composite laminates[END_REF], [START_REF] Nairn | Matrix microcracking in composites[END_REF], [START_REF] Boniface | Transverse ply cracking in cross-ply cfrp laminates-initiation or propagation controlled[END_REF]) is also studied at this scale.

2. Diffused damage at the interface or micro-delamination: it shows through the appearance of micro-cracks at the resin interface connecting plies.

3. Fibre breakage: it consists in the rupture of fibres.

• Meso models: The common dimensions are approximately around (≈ 0.1mm), which corresponds to the ply scale. The main hypothesis implemented in these models is that a laminate can be modeled from a stacking sequence of two main constituents generally corresponding to the interface and the ply. The main damage mechanism is transverse cracking. These transverse cracks are modelled mainly through two methods: Continuum Damage Models and Discrete Damage Models. At this scale, delamination and final fibre breakage also occur if the specimen is subjected to critical loadings.

1.4. A short introduction to numerical approaches for the simulation of damage mechanisms at different length scales 35

• Macro models: The common dimensions are approximately around (≈ cm) which corresponds to the scale of the laminate. At this length scale problems are usually modeled bi-dimensionally following the Classical Laminated Plate Theory (CLPT). These models are typically based on a two step process which involves calculating the displacement/strain and stress fields and the use of damage criteria.

Continuous models aim at describing damage through classical continuum mechanics approaches and principles, based mainly on the homogenisation of the damaged ply behaviour (usually) at the micro scale. A popular approach, as discussed in Abisset [2012], defined as micro-hybrid, is the one developed by [START_REF] Ladeveze | A computational damage micromodel of laminated composites[END_REF] at the LMT Cachan allowing to describe damage mechanisms up to millimeter-precision see [START_REF] Violeau | Micromodel-based simulations for laminated composites[END_REF]. Discrete damage models still treat the ply behaviour in a continuous context, however the cracks are introduced through a discrete approach.

In this sequel a short introduction to cohesive zone models and continuum damage models will be given as well as for to the variational approach to fracture applied to composites. For a more in-depth introduction to the Variational approach to fracture, see Chapter 4.

Cohesive zone models (CZM)

The concept of cohesive zone model was first introduced by [START_REF] Barenblatt | The mathematical theory of equilibrium cracks in brittle fracture[END_REF] in order to cope with stress singularities around the crack tip where, by definition in linear elastic fracture mechanics (LEFM), stresses tend to infinity. This method consists in avoiding a zone of vanishing thickness ahead of the crack tip. The concept is based on considering surfaces to be held together by a cohesive traction. When the material fails, these cohesive surfaces separate and the process is described by a cohesive law.

Barenblatt's model acts at an atomic level, where these crack surfaces are considered to be held by atomic bonding forces. The advantage, as stated in [START_REF] Sun | Chapter 9 -cohesive zone model[END_REF], is that it has an intrinsic fracture energy dissipation mechanism compared to LEFM, however the direct disadvantage is the problem in linking different scales, from the atomic to the continuum.

Cohesive zone models act in a phenomenological framework that aims at modelling crack nucleation and its propagation. The main advantages that arise compared to LEFM are the non-necessity of a pre-crack and, as mentioned previously, the avoidance of stress singularities. More precisely, this is achieved through the introduction of Traction-Separation-Laws (TLS). In CZM the creation of a crack implies the creation of new unconnected surfaces. Let's suppose that these surfaces are called Γ + and Γ -and we call δ the distance between them, see Figure 1.9. The zone around the crack tip is called process or cohesive zone.

Generally, the cohesive traction and the distance between the two surfaces (separation) are the key arguments implemented in a dimensionless function, such that this function can be written as:

σ = σ c f (δ/δ c ) (1.1)
A general cohesive law is shown in Figure 1.10 Where σ is the cohesive traction, σ c is the peak cohesive traction, δ is the separation or displacement between Γ + and Γ -and δ c is the critical or characteristic separation value. This of course depends on the failure mechanism (fragile, ductile...), the opening mode and the material in question. However, [START_REF] Barenblatt | The mathematical theory of equilibrium cracks in brittle fracture[END_REF] had originally formulated a function of the total distance (end to end) of the crack in the cohesive zone, which complexified the model leading to its stand-by for several years until Finite Element Methods (FEM) allowed to implement and solve these partial differential systems, as emphasized in [START_REF] Nordmann | Cohesive zone models-theory, numerics and usage in high-temperature applications to describe cracking and delamination[END_REF].

There is a vast multitude of TSLs used to implement different CZMs, to mention a few:

• The Dugdale model, introduced by [START_REF] Dugdale | Yielding of steel sheets containing slits[END_REF] primarily to describe ductile failure

• Linear Softening Model, introduced to describe quasi-brittle failure (ceramics, concrete)

• Trapezoidal Model: this an extension of both the Dugdale model and the linear softening model. Depending on the chosen parameters, it can either represent one or the other method.

• Exponential Model: it's based on the universal binding energy curve, introduced by Rose et al.

[1981], used to simulate crack growth in ductile metals.

To conclude, cohesive zone models are often implemented to study delamination as the interface can be modeled as a cohesive zone, therefore they are appropriate and efficient in this line of work. The main drawback of this technique, however, is the necessity to know the crack path in advance [START_REF] Nordmann | Cohesive zone models-theory, numerics and usage in high-temperature applications to describe cracking and delamination[END_REF]).

1.4. A short introduction to numerical approaches for the simulation of damage mechanisms at different length scales 37

Continuum damage models

The main idea of continuum damage models relies on the homogenization at the chosen scale of the REV. The choice of the REV will determine the nature of the model. For instance, the REV can be modelled at the fibre scale, where the elementary volume consists of the fibre and the matrix mainly, [START_REF] Costanzo | Micromechanics and homogenization of inelastic composite materials with growing cracks[END_REF], [START_REF] Murari | Micromechanics based ply level material degradation model for unidirectional composites[END_REF]. The scale being focused only on these two constituents, the main drawback is the inability to foresee and model damage mechanisms related to the fibre-matrix interaction and the ones related to a higher scale. The next step would be a full mesoscale approach. At this level, the ply should represent one or more unidirectional individual plies modelled through a homogenized volume, while the interface should represent the resin rich zone found between individual plies stacked in different directions. Many models at this scale can be found in the following works [START_REF] Williams | Application of a damage mechanics model for predicting the impact response of composite materials[END_REF], [START_REF] Ladevèze | On a damage mesomodel for laminates: micro-meso relationships, possibilities and limits[END_REF], [START_REF] Kashtalyan | Analysis of composite laminates with intra-and interlaminar damage[END_REF], [START_REF] Maimí | A continuum damage model for composite laminates: Part ii -computational implementation and validation[END_REF], [START_REF] Laurin | A multiscale progressive failure approach for composite laminates based on thermodynamical viscoelastic and damage models[END_REF], [START_REF] Van Der Meer | A level set model for delamination -modeling crack growth without cohesive zone or stress singularity[END_REF]. Depending on the damage mechanism, different approaches are implemented. For instance, delamination is modelled by placing cohesive interfaces between the neighbouring individual layers to connect them, see [START_REF] Shi | Modelling damage evolution in composite laminates subjected to low velocity impact[END_REF], [START_REF] Li | Modeling interfacial debonding and matrix cracking in fiber reinforced composites by the extended voronoi cell fem[END_REF]. Still the main drawback to the plybased models is cost-related, due to the fact that each layer of the laminate is numerically modelled and simulated, the computational cost becomes quickly expensive and the application to a real-sized structure quite difficult. The second drawback, as emphasized in [START_REF] Forghani | An overview of continuum damage models used to simulate intralaminar failure mechanisms in advanced composite materials[END_REF], is linked to the independent behaviour of the layers, since the mesoscale model assumes that the behaviour of a ply is independent of its neighbours, while individual plies behave differently when stacked in multidirectional sequences, see [START_REF] Camanho | Prediction of in situ strengths and matrix cracking in composites under transverse tension and in-plane shear[END_REF]. The third drawback is related to the definition of the damage evolution law which is at this scale empirical and/or phenomenological. This is mainly due to the fact that transverse cracks rely on the coalescence and growth of damage at the microscale. This could be overcome by bridging micro and mesoscales, as it was done in [START_REF] Ladevèze | On a damage mesomodel for laminates: micro-meso relationships, possibilities and limits[END_REF], [START_REF] Ladevèze | Towards a bridge between the micro-and mesomechanics of delamination for laminated composites[END_REF], [START_REF] Ladeveze | A computational damage micromodel of laminated composites[END_REF]. They have developed a micro-hybrid model that allows through homogenisation at the microscale to generalise the mechanical response to the mesoscale. More precisely, the model is based on a micro-mechanical approach that is coupled with numerical approaches to generalise the model and treat structural problems. This model is defined to be hybrid since the cracks are introduced discreetly and introduced through cracking surfaces. However, continuum damage mechanics is applied to describe progressive degradation mechanics such as fibre-matrix debonding. Furthermore, the evolution of transverse micro-cracks is described through a simplified homogenized law. In fact, the main idea of the model is to homogenise the strain and stresses at a microscale on a representative elementary volume (REV) in order to reduce the numerical cost when implemented at mesoscales. The chosen REV can either be at the fibre scale or at the transverse crack scale. The fibre scale consists in a REV constituted of mainly the fibre and the matrix. The main drawback to this scale, chosen for homogenisation, is that it is limited to the prediction of fibre-matrix damage types. At the transverse crack scale, suggested by [START_REF] Ladevèze | On a damage mesomodel for laminates: micro-meso relationships, possibilities and limits[END_REF], [START_REF] Ladeveze | A computational damage micromodel of laminated composites[END_REF], a connection between micro and mesoscale has been developed. On a side note, the numerical model assumes that the micro-cracking patterns are periodic since as stated in the publication the assumption is "consistent with most practical situations" [START_REF] Ladevèze | Towards a bridge between the micro-and mesomechanics of delamination for laminated composites[END_REF]). The benefit of this model is the precise description of the damage evolution, however the main drawback is the limitation to simulate to total thicknesses of just a few millimetres.

Finally, macro models are the last approach, also called sub-laminate approaches. These are the only ones that can be applied to simulate large-scale structures. On the downside, the damage mechanism cannot be predicted in detail, down to the micro or mesoscale. The choice of the model is made Chapter 1. State of the art of matrix cracking in unidirectional (UD) laminated composites depending on the problem one seeks to solve.

To conclude, continuum damage models have found success in the modelling of damage in composites, specifically in the case of intralaminar damage involving fibre breakage and matrix-cracking at the ply scale. Mesoscale damage behaviour laws are now widespread in finite element modelling tools. In this context, numerical simulations are performed with at least one element in the ply thickness, which leads to the main limitation of such an approach related to the computational cost to simulate a full scale structure. Furthermore, stress-softening material behaviour lead to numerical issues such as the choice for mesh orientation, see [START_REF] Forghani | An overview of continuum damage models used to simulate intralaminar failure mechanisms in advanced composite materials[END_REF] for a complete review.

Variational approaches to fracture

Phase field models belong to a family of variational approaches and are often referred to as gradient damage models or smeared crack models. These models are particularly powerful as they are able to predict crack paths and branching without needing to know the cracking patterns a priori, an innovating feature since Griffith. They rely on energy minimization to describe the damage field, which is described though a scalar function that allows to modulate the material stiffness, smeared across a band of a specific width, determined by an internal length variable, as opposed to continuum models where "cracks are viewed as zones where a loss of stiffness continuously occurs along localized bands of non-zero thickness" [START_REF] Marigo | An overview of the modelling of fracture by gradient damage models[END_REF]. They differ from local damage models, as in fact in local damage models the material failure is managed at a material point through a critical stress that allows to define the elastic critical threshold and damage criterion. Predicting stress-hardening problems with local models does not raise any issues as the evolution problems are well-posed and the response unique, however the issue is related to stress-softening problems as they appear to be ill-posed, [START_REF] Marigo | An overview of the modelling of fracture by gradient damage models[END_REF]. These issues are overcome by adding a gradient term of the damage variable, thus a non-local term, to the definition of the energy functional. An in-depth introduction to these models can be found in Chapter 4.

Phase field models and their application to composite materials

In this section we will discuss the application of phase field models, which are variational approaches to fracture, in the case of composite applications and what are the main issues related to the accurate simulation of damage mechanisms such as matrix-cracking and delamination. An in-depth introduction to variational approaches and the phase field model is done in Chapter 4.

The application of phase-field models in brittle fracture is a delicate and still open question when it comes to anisotropic materials, such as composites. The application of these models has seen many fruitful results in the context of quasi-static, isotropic and homogeneous cases, even though applications to dynamic and anisotropic media is possible. We will not be concerning ourselves with the dynamic approach to phase field, as it goes beyond the scope of this study. However, a modelling of an orthotropic case will be discussed further on and therefore a short overview of a few studies in the field will be mentioned in this section.

In the quasi-static setting, the classical theory for predicting crack propagation includes Griffith's propagation criterion and a crack path criteria such as the maximum energy release rate, the principal of local symmetry and the maximum circumferential stress. These different criteria, predict slightly different results despite the fact that the material has isotropic surface energy [START_REF] Sumi | Mathematical and Computational Analyses of Cracking Formation[END_REF]), nevertheless, as stated in [START_REF] Li | Crack kinking in a variational phase-field model of brittle fracture with strongly anisotropic surface energy[END_REF], the relative deviations of the experimental evidences are nearly indiscernible.

In the scope of applying phase field models to composites, several studies have been published, to cite a few remarkable ones [START_REF] Quintanas-Corominas | A phase field approach to simulate intralaminar and translaminar fracture in long fiber composite materials[END_REF] proposed as well a new phase field model suitable for interlaminar and intralaminar failure in long fibre-reinforced composites. The model is capable of reproducing complex failure mechanisms within these laminates. This is achieved though the implementation of Gibb's free energy form release rates. In 2018 [START_REF] Msekh | Fracture properties prediction of clay/epoxy nanocomposites with interphase zones using a phase field model[END_REF] proposed a phase field model to predict damage mechanisms in clay/epoxy nanocomposites. [START_REF] Bleyer | Phase-field modeling of anisotropic brittle fracture including several damage mechanisms[END_REF] proposed a model that aims at extending the phase field modeling to orthotropic elastic media to represent longitudinal and transverse damage (LTD) for failure on UD laminates, they use two phase field variables (different internal lengths) to characterize longitudinal damage and transverse damage. [START_REF] Zhang | A double-phase field model for multiple failures in composites[END_REF] suggest as well an approach based on a double phase field variable model, each one corresponding to the type of damage. They introduce a new kind of strain energy which includes four effective strain variables that have invariance with respect to the fibre orientation, as well as a new degradation tensor.

From a more general perspective, based on the fracture properties, materials can be categorized either as weakly anisotropic or strongly anisotropic. Weakly anisotropic materials are characterized by a convex surface energy function dependent on the normal to the crack, therefore cracks are orientation-dependent. Whereas, strongly anisotropic materials have a non-convex surface energy function dependent on the angle of the crack and fracture toughness is a non-convex function of the normal as well, see [START_REF] Li | Crack kinking in a variational phase-field model of brittle fracture with strongly anisotropic surface energy[END_REF]. Microstructure and manufacturing influence a great deal the surface energy and are mostly the reason why many materials have anisotropic surface energy. In the case of a material with isotropic elasticity and anisotropic surface energy the principle of local symmetry and stress based criteria lose pertinence, as stated in [START_REF] Li | Crack kinking in a variational phase-field model of brittle fracture with strongly anisotropic surface energy[END_REF], [START_REF] Nguyen | Phase field modelling of anisotropic crack propagation[END_REF]. [START_REF] Li | Crack kinking in a variational phase-field model of brittle fracture with strongly anisotropic surface energy[END_REF] propose a variational regularized approach for both strongly and weakly anisotropic materials. In both cases a term is added in the surface energy, having an energetically regularizing effect through the penalisation of the crack's curvature. They provide analytical formulas for the energy as a function of the crack's angle and an open source finite element solver for high order phase field problems. Their result shows how cracking patterns can be strongly affected by anisotropy in two main ways:

1. Anisotropy can stabilize straight cracks 2. Anisotropy can induce zig-zag cracking patterns.

It is worth mentioning that for weakly anisotropic materials, a second order tensor symmetric positive-definite is introduced in the surface energy and Γ-convergence was proved in [START_REF] Focardi | On the variational approximation of free-discontinuity problems in the vectorial case[END_REF]. Meanwhile in [START_REF] Li | Phase-field modeling and simulation of fracture in brittle materials with strongly anisotropic surface energy[END_REF], a model for both strongly and weakly anisotropic surface is given and in the case of strongly anisotropic surface energies a fourth order symmetric positive definite tensor is introduced. They emphasized that, in this case, the sharp interface model (Griffith's approach) is not well-posed due to the lack of semi-continuity, Γ-convergence can be proven only on the convexified variational approach.

2

Experimental procedure and results on carbon fiber reinforced UD cross-ply laminates through passive infrared thermography: a first analysis of the physical cracking patterns

The scope of this chapter is to dive into the core thematic addressed in this research. An experimental campaign prior to this work has been carried at the ONERA centre in Lille, where a new and innovative protocol has been implemented. A particular material is considered and will be at the core of the next chapters as well. An analysis of the data is discussed, with particular focus on the distribution of cracks and the parameters that play a role in it. In conclusion, this chapter contextualizes and motivates the physical and mathematical problem assessed further on.

Introduction

In this chapter, a particular experimental procedure, introduced by Berthe and Ragonet at ONERA in the research unit CRD of the DMAS department, is discussed. This procedure is applied to the T700GC/M21 carbon fiber-reinforced UD cross-ply laminates. These experiments were conducted prior to the works of this thesis and they constitute the main data set used for comparisons with the numerical simulations. A brief introduction of the material can be found in section 2.2, where a table sums up the mechanical and material properties. The context of the study is briefly introduced in section 2.3 while the experimental procedure and results of passive infrared thermography conducted by [START_REF] Berthe | Passive infrared thermography measurement of transverse cracking evolution in cross-ply laminates[END_REF] are reported in section 2.4. Finally, in section 2.5, we show the post-treatment analysis of the experimental data and a conclusion, in section 2.6, motivating the analytical and numerical studies of the next chapters.

Material properties of T700GC/M21

In this section we introduce the mechanical and material properties of the studied material, which is the T700GC/M21 carbon/epoxy. The laminate in question is a symmetrical [0 n /90 m ] s classical UD cross-ply, where the resistance is improved thanks to the inclusion of thermoplastic nodules in the epoxy resin. The material properties are detailed in table 2.1. This table collects all material properties of the tested specimens. These values are to be considered for both the experimental tests and the numerical implementations that will be discussed in later chapters.

-I G IC J/m 2 350 (0 • /0 • ) Mode-II G I IC J/m 2 1250 Tenacity*** Mode-I G IC J/m 2 300 (0 • /45 • ) Mode-II G I IC J/m 2 1600
Chapter 2. A first analysis of the physical cracking patterns

Context of the study and motivation for the proposed experimental protocol

The most commonly studied stacking sequence to characterize transverse cracking, as seen in Chapter 1, is the [0 m /90 n ] s cross-ply laminate. Different studies have been conducted over the years, and they all show that transverse cracking in [0 m /90 n ] s is instantaneous from side to side in the 90 • layers given a sufficiently high ply thickness , see [START_REF] Berthelot | Transverse cracking of cross-ply laminates: Part 1[END_REF], Parvizi and Bailey [1978], Parvizi et al. [1978]. Multiple different experimental techniques over the years have been applied, consisting in the use of: visible imaging sensor (see Parvizi et al. [1978], Parvizi and Bailey [1978]), X-ray radiography [START_REF] Crossman | Initiation and growth of transverse cracks and edge delamination in composite laminates part 2. experimental correlation[END_REF], [START_REF] Groves | An experimental and analytical treatment of matrix cracking in cross-ply laminates[END_REF]), optic microscope (see [START_REF] Varna | Damage in composite laminates with off-axis plies[END_REF], [START_REF] Huchette | A multiscale damage model for composite laminate based on numerical and experimental complementary tests[END_REF], [START_REF] Ogi | Empirical models for matrix cracking in a cfrp cross-ply laminate under static-and cyclic-fatigue loadings[END_REF]), digital image correlation (see [START_REF] Farge | Damage characterization of a cross-ply carbon fiber/epoxy laminate by an optical measurement of the displacement field[END_REF] , and finally through the use of high resolution cameras (see [START_REF] Sebaey | Measurement of the in situ transverse tensile strength of composite plies by means of the real time monitoring of microcracking[END_REF]).

These techniques require static loading conditions in order to be applied, therefore dynamic impacts or intermediate speed loadings are not experimental settings where damage mechanisms can be captured or monitored to characterize the damage behaviours. A valid alternative would be the application of high resolution rapid cameras (20MPix), however these kind of devices are not currently available in the market as they have not been developed. Hence the necessity for a new experimental protocol. In this context, a new experimental protocol has been developed at ONERA, introduced in [START_REF] Berthe | Passive infrared thermography measurement of transverse cracking evolution in cross-ply laminates[END_REF]) in which infrared thermography is used to capture intermediate speed loading tests causing transverse cracks in the core layers. The main concept of this method is to exploit infrared radiations emitted by the object and link them to the body temperature. The thermal field of the specimen can be obtained if the emissivity of the body is known and based on the calibration of the infrared camera.

Passive infrared thermography results and their analysis based on the influence of thickness and loading rates

Through the use of passive infrared thermography, in a first moment [START_REF] Berthe | Passive infrared thermography measurement of transverse cracking evolution in cross-ply laminates[END_REF] verified the hypothesis of instantaneous crack propagation from edge to edge in the specimens, the reader can find further details on the experimental procedure conducted on the specimens in [START_REF] Berthe | Passive infrared thermography measurement of transverse cracking evolution in cross-ply laminates[END_REF]. In a second moment, the captured images from the infrared cameras were analyzed through the use of an algorithm developed by [START_REF] Coussa | Rate dependency and ply thickness influence on transverse cracking evolution in cross-ply laminates[END_REF]. This numerical post-processing technique has been implemented to reproduce the evolution of matrix cracking patterns as well as showing the localisation of each crack on the specimens as a function of the applied loading, as it can be seen in Figure 2.1. The tested specimens are subject to uniaxial tension where different constant speeds are applied, four being the main ones applied in these tests: 5mm/min, 50mm/min, 500mm/min and 0.1m/s. Three main stacking sequences were tested, respectively: [0/90 1 ] s , [0/90 2 ] s and [0/90 3 ] s , each at the four different speeds previously mentioned. For each stacking sequence, at each testing speed, three specimens were tested. Following this procedure, an analysis on the crack density rates is done and on the influence of both, the thickness of each ply and, respectively, the influence of the applied speed loading condition. The crack density curves associated to these tests can be seen in Figures 2.2 and 2.3.

The influence of the thickness shows, that for the thinner configurations ([0/90/0]) at the tested rates, a limitation in the experimental set up appears as a greater number of pixels is needed. The spatial resolution is too low to capture crack occurrences in a consistent way. [START_REF] Coussa | Rate dependency and ply thickness influence on transverse cracking evolution in cross-ply laminates[END_REF] suggest to reduce the specimens free-length to overcome the spatial resolution issue. The curves in Figure 2.3 show the crack density curves corresponding to different loading rates and different total thicknesses. As it can be noticed, for the thinner stacking sequence the crack density curves show a different behaviour than for the thicker configurations. In fact, in Figure 2.2 we can notice that the crack density reaches a saturation for greater loadings, also transverse cracks appear in the specimen for greater loadings than for thicker configurations, knowing that the total free length of the specimen is kept constant. The trend of the crack density curves seems to follow a more "linear" evolution for the [0/90 2 ] s and [0/90 3 ] s than for the [0/90 1 ] s specimen, this can be easily seen in Figure 2.3. The critical stress value for crack nucleation, as it can be seen in Figure 2.3, is "delayed" as the specimen thickness increases. These trends are also reported in Parvizi and Bailey [1978], [START_REF] Nairn | Matrix microcracking in composites[END_REF], [START_REF] Berthelot | Transverse cracking and delamination in cross-ply glass-fiber and carbon-fiber reinforced plastic laminates: Static and fatigue loading[END_REF], as emphasized in [START_REF] Coussa | Rate dependency and ply thickness influence on transverse cracking evolution in cross-ply laminates[END_REF] and hence assure coherence in the obtained experimental results. Furthermore, the influence of the loading rate appears to be minimal on these curves for the thicker plies and is more remarkable for the thinner plies. However the loading rates .

do play an important role on the cracking patterns, even on the thicker configurations and this is observable in Figure 2.4, where it is visible that the cracking pattern evolves as oblique cracks appear for higher loading rates.

The carried study on the three stacking sequences proves that both the thickness and the loading rates influence the crack density rate and the cracking pattern. However these phenomena are not fully understood. On the basis of these results, the work of this thesis focuses on investigating the influence that the ply thickness, and the total specimen thickness, play on the cracking patterns. The influence of the lamina's thickness on the critical stress value that enables nucleation decreases as the the laminate increases in thickness, see Figure 2.2 where in blue, green and red the following stacking sequences [0/90] s , [0/90 2 ] s and [0/90 3 ] s are respectively reported. In fact, for the [0/90] s specimen nucleation starts for a stress of approximately 300 MPa, while for the thicker specimens the critical stress is around 400 MPa. In the ordinates, the crack density sees an increase proportional to the applied stress and exceptional high values are reported for the thicker specimen.

As reported in Parvizi et al. [1978], depending on the thickness of the 90 • layers, the damage criterion influencing nucleation can either be described or controlled through a stress or an energetic criterion. The limit value between the criteria can be established based on the material parameters. In Figure 2.5, the image illustrates an example of such criterion proposed by Parvizi et al. [1978] for a cross-ply glass fiber-reinforced [0 • /90 • ] s composite, where the thickness of the 90 • plies evolves up to 4mm. In the experiments presented in this Figure, the transverse cracking behaviour changed as the inner-ply thickness decreased below about 0.4mm. Based on Parvizi et al. [1978]'s results, their conclusion affirms that for low-ply thicknesses an energetic criterion would be more suitable than a strength-based criterion, since they observe that energetics control crack growth rather than the mechanism.

Post-Processing analysis of matrix-cracking results in the tested specimens

In this section, a new analysis of the available experimental data is proposed in order to provide a complementary view to the existing results in the literature. In all of the literature presented in 2.5. Post-Processing analysis of matrix-cracking results in the tested specimens 49 Chapter 1, although experiments were mostly conducted on glass reinforced polyester resin crossply laminates, the cracking patterns observed are uniform, periodic. Based on the presented literature, we take an interest in observing the mechanism in which cracks appear and are distributed in the specimen, this could either happen exactly midway, which is what is mostly observed [see e.g. Parvizi and Bailey, 1978]), or, in a longer segment and a shorter segment. The latter, is what is mostly observed in our experimental procedures. This type of analysis of the statistical distribution of segments' lengths before and after cracking is quite rare in the literature. In fact, it requires a fine monitoring of the appearance of the different cracks and of the measuring of the generated segments in the 90 • ply. Moreover, in the case of carbon fibre composites, the analysis of transverse cracking requires the test to be interrupted in order to take micrographs or X-rays for example. In this context, it is very difficult, if not impossible, to do this perfectly for each crack. In this context, the new experimental protocol previously presented also makes it possible to propose an analysis of the statistical distribution of the segments formed. This new data analysis has been proposed in this thesis to enrich the modelling that will be carried out later.

In the sequel, the analysis is mostly focuses on low loading rate tests (5 mm/min, 50 mm/min and 500 mm/min) and for the thicker configuration ([0/90 3 ] s ). Other speeds were tested, as shown in Figures 2.3-2.4, however these are the main ones we focus on this work, since our entire study will be based on the development of a quasi-static approach. First, how the length of each segment is calculated from the data is explained. Following, crack distribution for the three tested speeds is shown. An analysis on how a segments breaks over time and distributes over the free length is made as well as the percentage long and short segments constitute. From the collected data, the pixel position for each appearing crack is known. From these values we calculate the normalized length of each segment, that here henceforth will be named ℓ i , so that: ℓ i = (p fp i )/p tot . Where p f is the pixel position at which the current crack appears, while p i indicates the pixel position of the previous appearing crack. If no pre-existing crack exists after the first apparition, then the extremity values of the whole effective zone are considered. Finally, p tot indicates the total number of pixels used to capture the corresponding total length of the specimen. For instance, let's consider the case where the effective area analyzed from the camera is within the range x ∈ [25, 306], as shown in Figure 2.6. The total length of the analyzed zone is: p tot := 306 -25 = 281. A new crack appears at x = 284, the lengths, normalized with respect to the total length of the specimen, of the two generated segments are: ℓ 1 := (284 -25)/281 = 0.9217 and For each appearing crack, it appears that a longer segment and a shorter segment are always generated. This creates a cracking pattern suggesting an aperiodic distribution of the cracks could be the final appearing network in the specimen. This pattern, obtained from the experiments conducted by [START_REF] Berthe | Passive infrared thermography measurement of transverse cracking evolution in cross-ply laminates[END_REF], has been reconstructed and is shown in Figure 2.7. For this specific thickness configuration, at none of the tested speeds, a periodic network of cracks is obtained.

The evolution of the greatest segment for each appearing crack is shown in the histogram in Figure 2.8. These three figures show the length of the greatest segment as a function of the number of cracks for three different upper holder speeds, respectively 5 mm/min, 50 mm/min and 500 mm/min. The first observation that can be made, is that, for slower tests (closer to quasi-static testing) such as the 5 mm/min and the 50 mm/min the tendency is to have a whole segment that will break in a clear longer and shorter one. In fact the corresponding values in this case, for one appearing crack n = 1, are: ℓ 1 = 0.92, ℓ 2 = 0.08 for a speed test of 5mm/min and ℓ 1 = 0.96, ℓ 2 = 0.04 for a speed of 50mm/min. Whereas for the higher speed tests the segments tend to break midway or create two segments of similar size. For each specimen, three tests were conducted and the same behaviour was observed in each graph, the value of each segment changes as the position of each crack is "random" and is influenced by multiple factors, however the trend of these histograms remains the same.

In these histograms we report the longest segment at each appearing crack, and not the evolution of the length of the segment being broken into multiple pieces. We can observe that the longest segment can remain integer for two or multiple appearing cracks, meaning that a shorter segment does break in its place. Therefore, the longest segment does not always separate into a long and short piece first, however it does occur eventually. We want to better understand how the segment is divided in two pieces. For that purpose, the ratio between the longest generated segment after breakage and the full length of the segment prior to breakage is calculated. The percentage thus obtained is plotted in the histograms shown in Figure 2.9. These figures show the percentage for which the largest segment will break first, and as it can be seen, as the upper holder speed increases, the probability of breaking into a long and a short segment in the range between 90 -100% goes to zero. As the upper holder speed increases we can observe a decrease in all ranges from 60 -100%, this means that the cracks will tend to be more evenly spaced. This is also influenced by the density of the network as well, in fact, as the number of cracks increases, for all cases, the segments will be more evenly spaced. We can expect this result also by looking at Figures 2.8, in fact we see that as the number of cracks increases and the density reaches a saturation "peak", we can observe that the greatest segment remains the same.

Conclusion

The crack density rate is clearly influenced by the thickness of the entire laminate, therefore by the number of plies in the core. We can closely verify this statement by looking at Figure 2.3. As the layers increase, it appears that the loading rate has a slight influence on the crack density rate, see Figure 2.3c, and the cracking pattern seems to be influenced as well, see Figures 2.7, 2.8. Many are the cited studies in Chapter 1 that have assumed a periodic pattern to be granted or to be a standard occurrence. Few studies adventured in the description of aperiodic distribution. No real analysis on aperiodic patterns has been published based on the study of an elementary cell without the use of any statistical tools. We analyzed the experimental data by looking at: the crack density rate, the reconstructed cracking pattern, the evolution of the cracked segment's length and the percentage for which the longest segment tends to break first as the applied speed increases. Based on the experimental results, it is difficult to generalize for which parameters a certain behaviour occurs.

From these results we can conclude that both the thickness and the loading rate do influence both the cracking pattern and the crack density rate. However it is difficult to conclude a criterion for which the pattern will be periodic or aperiodic based only on this analysis.

The work presented in the next chapters is based on both an analytical and numerical bi-dimensional approach that identifies how material and geometrical parameters affect the cracking pattern and their distribution in a quasi-static setting by simply looking at the convexity of the energy functional describing the system. Based on a theorem from [START_REF] Puglisi | Mechanics of a discrete chain with bi-stable elements[END_REF], we will prove how the segmentation occurring in a long and short piece is stable and thus energetically favoured for a certain set of parameters and therefore can be implemented for the experimental set of parameters for comparison. 3

Analytical and numerical approaches implemented to study transverse cracking patterns in UD cross-ply laminates

We study the problem of transverse cracking in cross-ply laminates through a variational approach to fracture. In this chapter, we propose a new approach to the study of this issue, bypassing any statistical criteria of defect distribution. Instead, we prove that crack distribution depends on the concavity of the elastic compliance function of the elementary cell between two cracks, as a function of the cell length. We will introduce a Griffith like functional of the total energy of the modelled composite bar, and formulate the optimisation problem on a cell-based approach. Solutions of the optimisation problem are computed numerically, through an algorithm that we've implemented ourselves, and appear in bifurcation diagrams where their stability is also studied. Furthermore phase diagrams to show the evolution of the energy with respect to the material properties have also been computed. Finally, a reduced one-dimensional shear-lag approach is also presented.

Introduction

As introduced in the bibliographical study, we know that cross-ply laminates under uni-axial tension exhibit primarily transverse cracking patterns that have been exhaustively investigated experimentally in reviews such as [START_REF] Nairn | Matrix microcracking in composites[END_REF] and [START_REF] Berthelot | Transverse cracking and delamination in cross-ply glass-fiber and carbon-fiber reinforced plastic laminates: Static and fatigue loading[END_REF]. Therefore, the appearance of multiple transverse cracks in the brittle core is the predominant failure mechanism, and it is the object of the present work.

We remind the reader that the available experimental results [see e.g. Parvizi and Bailey, 1978[START_REF] Flaggs | Experimental determination of the in situ transverse lamina strength in graphite/epoxy laminates[END_REF][START_REF] Peters | Statistical analysis of multiple fracture in [0/90/0] glass fiber/epoxy resin[END_REF][START_REF] Berthe | Passive infrared thermography measurement of transverse cracking evolution in cross-ply laminates[END_REF] indicate the following fundamental features: (i) cracks appear for an in situ stress in the core larger than the nominal strength of the standing alone lamina; (ii) the average crack spacing is decreasing with the applied loading; (iii) the in situ strength increases when the relative core thickness decreases. Energetic fracture mechanics models predict reasonably well the in situ strength and the evolution of the crack spacing with the loading [START_REF] Nairn | Matrix microcracking in composites[END_REF], [START_REF] Berthelot | Transverse cracking and delamination in cross-ply glass-fiber and carbon-fiber reinforced plastic laminates: Static and fatigue loading[END_REF].

The available models range from basic shear-lag models Parvizi et al. [1978] to full-field finite element numerical computations, passing through sophisticated variational approximations [START_REF] Hashin | Analysis of cracked laminates: A variational approach[END_REF], [START_REF] Berthelot | Transverse cracking of cross-ply laminates: Part 1[END_REF]. The energetic models are in agreement with experiments until reaching a crack saturation regime that can be explained only by introducing additional length scales [START_REF] Dvorak | Analysis of first ply failure in composite laminates[END_REF], [START_REF] García | Transverse crack onset and growth in cross-ply [0/90]s laminates under tension. Application of a coupled stress and energy criterion[END_REF], [START_REF] Leguillon | Strength or toughness? A criterion for crack onset at a notch[END_REF].

The existing works assume the crack distribution in the core to be either uniform, or determined by statistical criteria of defect distribution [START_REF] Manders | Statistical analysis of multiple fracture in [0/90/0] glass fiber/epoxy resin laminates[END_REF], [START_REF] Berthelot | Statistical analysis of the progression of transverse cracking and delamination in cross-ply laminates[END_REF], [START_REF] Javaland | Numerical modeling of diffuse transverse cracks and induced delamination using cohesive elements[END_REF], [START_REF] Okabe | Numerical modeling of progressive damage in fiber reinforced plastic cross-ply laminates[END_REF], [START_REF] Silberschmidt | Matrix cracking in cross-ply laminates: Effect of randomness[END_REF], [START_REF] Vinogradov | Probabilistic energy based model for prediction of transverse cracking in cross-ply laminates[END_REF], [START_REF] Maimí | Numerical modeling of matrix cracking and intralaminar failure in advanced composite materials[END_REF].

In the present work, we tackle the optimal crack spacing problem within the variational approach to fracture Francfort and [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF] and a Griffith crack model, leveraging the ideas presented by [START_REF] Bourdin | The Variational Approach to Fracture[END_REF]. Considering a layered composite bar with a set of n transverse cracks in the core, we show that whether their spacing is uniform or not depends on the concavity of the elastic compliance function of the elementary cell between two cracks, S(ℓ), as a function the cell length ℓ. We accurately compute the first and second derivatives of the elementary elastic compliance by the finite element method and domain-derivative techniques. Our results show that its second derivative change in sign at a critical length ℓ c . Assuming local energy minimality as a stability criterion, we show that the solution with uniform crack spacing is stable only when the cell length is smaller than ℓ c . When several cracks appear simultaneously, their spacing is either homogeneous or built on two different cell lengths, a short and a long one. In analogy with the behavior of a chain of bi-stable springs studied by [START_REF] Puglisi | Mechanics of a discrete chain with bi-stable elements[END_REF], we show that only solutions with a single long cell are energetically stable. This chapter is organized as follows. Section 3.2 formulates the general problem of optimal crack spacing and illustrates the possible bifurcation from uniformly to non-uniformly crack spacing through direct numerical computations. Section 3.3 presents the modular approach, the related optimality conditions, and the numerical method for the calculation of the compliance of the basic cell and its derivative by domain derivate techniques. Section 3.4 discusses the basic properties of the solution of the crack spacing problem and Section 3.5 reports on the influence of the material and numerical parameters. Section 3.6 outlines standard results for calculating the optimal number of Chapter 3. Analytical and numerical approaches implemented to study transverse cracking patterns in UD cross-ply laminates cracks. Section 3.7 introduces a reduced shear-lag model of the problem and shows the limitations of this approach. Section 3.8 closes the work.

3.2 2-Dimensional approach: Static variational problem for a given number of transverse cracks in the core

Problem formulation

We consider a layered composite bar of length L in Figure 3.1 composed of two identical outer layers of thickness H o and a core layer of thickness H c . We study transverse cracking in the core under the following assumptions:

• (i) the outer layers are linear elastic and unbreakable, with an elastic stiffness tensor A o ;

• (ii) the core is perfectly brittle, being linear elastic up to fracture, with an elastic stiffness tensor A c and a fracture toughness G c .

• (iii) cracks are possible only in the core (c) and are supposed to be in the form of a finite number of n -1 straight cracks orthogonal to the axis of the bar and spanning the whole thickness of the core (see Figure 3.1);

• (iv) we assume a 2d plane-stress model in the xy plane and neglect geometric non-linearities.

Oblique cracks and delamination, usually appearing after transverse cracks in the core, are not considered here. The 2d domain

Ω ≡ {(x, y) ∈ [0, L] × [-H c -H o , H c + H o ]}
is the reference configuration of the solid and

Γ(ℓ) ≡ ∪ n-1 i=1 Γ i (ℓ) with Γ i (ℓ) ≡ x = xi := i ∑ j=1 ℓ j , y ∈ [-H c , H c ] (3.1)
is the crack set according to the assumption (iii) above, where xi are the crack axial positions and ℓ := {ℓ i } n i=1 the corresponding crack spacings, see Figure 3.1. We study the fracture problem in the framework of the variational approach to fracture Francfort and [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF], [START_REF] Bourdin | The Variational Approach to Fracture[END_REF], where the cracked status of the system is determined by a minimality condition on a global energy functional, sum of the elastic energy of the cracked body and the surface energy due to the cracks.

3.2. 2-Dimensional approach: Static variational problem for a given number of transverse cracks in the core 59

Within the Griffith fracture model, the surface energy required to create the cracks is proportional to the crack surface. For our case, the surface energy is 2 (n -1)G c H c . Denoting by u : x = (x, y) ∈ Ω → u(x) ∈ R 2 the displacement field, the total energy functional to be minimized writes as

F (u, ℓ) := P tot (u, ℓ) + 2 (n -1) G c H c , P tot (u, ℓ) := Ω\Γ(ℓ) 1 2 A(x)ε(u(x)) • ε(u(x)) dx, (3.2)
where the first term is the elastic energy of the cracked solid, ε(u) := sym(∇u) being the linearized strain tensor. We consider the general case of orthotropic constituent materials under plane-stress condition, for which, using the Voigt notation:

A m∈{o,c} = 1 1 -ν xy m ν yx m   E x m E x m ν yx m 0 E y m ν xy m E y m 0 0 0 (1 -ν xy m ν yx m )G xy m   , ε(u) =     ∂u x ∂x ∂u y ∂y ∂u x ∂y + ∂u y ∂x     . (3.3)
For isotropic materials with Young modulus E m and Poisson ratio ν,

E x m = E y m = E m , ν yx m = ν xy m = ν m , and
G xy m = E m /(2(1 + ν m )).
We assume that the bar is loaded in tension by a hard device imposing a vanishing axial displacement on the left-end side and the axial displacement ū on the right-end side. The space of admissible displacement field C ū must respect these Dirichlet boundary conditions and should be regular except on the cracks, where it can jump, namely

C ū(ℓ) = {u ∈ H 1 (Ω \ Γ) : u 1 (0, x 2 ) = 0, u 1 (L, x 2 ) = ū, ∀x 2 ∈ [-H c -H o , H c + H o ]}, (3.4)
where H 1 (Ω, R 2 ) is the usual Sobolev space of functions with square integrable first derivatives defined on Ω and with values in R 2 .

We apply here the variational approach to fracture [START_REF] Bourdin | The Variational Approach to Fracture[END_REF], to determine the crack spacings ℓ's for a given loading and number of cracks n -1. We will discuss only briefly in Section 3.6 how to determine the optimal number of cracks as a function of the loading.

The static variational problem

We denote as the static problem the determination of the cracked state of the solid for a fixed loading (here the imposed displacement ū). In the static problem, we ignore any history effect and irreversibility condition on the evolution of the cracks.

Adopting the variational approach to fracture and the hypotheses above, this problem consists in solving the following minimization problem for the total energy functional (3.2):

min u,ℓ,n {F (u, ℓ), with u ∈ C ū(ℓ), ℓ ∈ R n , n ∑ i=1 ℓ i = L, n ∈ N}.
(3.5)

Differently from previous works [START_REF] Nairn | Matrix microcracking in composites[END_REF], [START_REF] Berthelot | Transverse cracking and delamination in cross-ply glass-fiber and carbon-fiber reinforced plastic laminates: Static and fatigue loading[END_REF], we do not assume any hypothesis on the crack spacings, neither we introduce statistical distribution of defects that can influence the fracture toughness or the energy release associated with a crack. Thanks to the hypothesis of linear elastic behavior in Ω \ Γ, we can uniquely determine the displacement field and the corresponding elastic energy as a function of the crack spacings ℓ, as follows:

u * (ℓ) = argmin{P tot (u, ℓ), u ∈ C ū(ℓ)}, P tot (ℓ) = P tot (u * (ℓ), ℓ).
(3.6)

Chapter 3. Analytical and numerical approaches implemented to study transverse cracking patterns in UD cross-ply laminates

Hence, the variational problem (3.5) can be reduced to the following constrained optimization problem in R n :

min ℓ,n F(ℓ), ℓ ∈ R n , n ∑ i=1 ℓ i = L, n ∈ N , with F(ℓ) := P tot (ℓ) + 2 (n -1) G c H c . (3.7)
As far as the number of cracks n -1 is given, the surface energy is fixed, being independent of the crack spacing. Hence, the key problem is to find the optimal spacing ℓ to minimize the elastic energy P tot (ℓ).

An example of solutions with non-uniform spacing: a bar with a single crack

Let us consider the simple example of a bar of length L with a single crack in the core layer, splitting it in two segments of length ℓ 1 = ℓ and ℓ 2 = L -ℓ, as shown in Figure 3.2. For this case, the variational problem (3.7) reduces to min 0<ℓ<L

P L (ℓ) = P tot (ℓ, ℓ -L). (3.8)
For a given bar length L and crack spacing ℓ, we compute the elastic energy P L (ℓ) by solving the linear elastic problem (3.6). We use a standard finite element approach with an unstructured uniform mesh and the geometry shown in Figure 3 • For long bars (L > L (c) ≃ 7H c ), the elastic energy has two minima at ℓ * 1 , ℓ * 2 = L -ℓ * 1 , corresponding to a crack dividing the bar in a short and a longer segment.

This simple example shows that the problem of the optimal crack spacing is non-trivial. Solutions with uniform spacing can be unstable for sufficiently long bars. The commonly accepted assumption [see e.g. [START_REF] Nairn | Matrix microcracking in composites[END_REF][START_REF] Berthelot | Transverse cracking and delamination in cross-ply glass-fiber and carbon-fiber reinforced plastic laminates: Static and fatigue loading[END_REF]] that the optimal crack position is at the center of the bar is not always in agreement with the energy minimality criterion: for sufficiently long bars the optimal crack spacing can be not uniform, even in perfect systems. This observation calls for a more general analysis to better understand the case with several cracks and to investigate the dependence on the material and geometric parameters.

Modular energetic approach for the case with n cracks

When considering a bar with several cracks, calculating the function P tot (ℓ) requires solving the linear elasticity problem (3.6) for each combination of spacing ℓ. A direct numerical strategy is not practical for more than one crack. More importantly, it does not allow us to unveil the key property of the solution. We develop here a modular approach, where the elastic energy is expressed as a function of the elastic compliance of an elementary cell.

The elastic energy as a function of the crack spacing

We regard the cracked beam as modular system composed of the boundary and bulk modules sketched in Figure 3.1. This decomposition requires the axial displacement to be constant throughout the outer layers in the cracked cross-sections, a mild assumption that is closely verified in practice. Let us denote by u i the displacement for the cross section x = xi of the crack Γ i and by ∆u i = u i+1u i the elongation imposed to the i-th cell of length ℓ i . Exploiting the symmetries, both modules can be regarded as a special cases of the elementary cells in Figure 3.1. By the linearity of the problem, the elastic energy of the elementary cell of length ℓ is quadratic with respect to the imposed elongation ∆u and can be written in the form

P cell (ℓ) = 1 2 K(ℓ) ∆u 2 , (3.9)
where the scalar valued function K(ℓ) is the global stiffness of the elementary cell. The function K(ℓ) can be calculated by solving an elementary cell problem, as it will be discussed in the next section.

Because of the symmetries, the potential energy of the boundary and bulk cells in Figure 3.1 are then given by

P 1 (ℓ) = 2 K(ℓ) ∆u 2 1 2 , P n (ℓ) = 2 K(ℓ) ∆u 2 n 2 , P i (ℓ) = 4 K(ℓ/2) (∆u i /2) 2 2 for i = 2, . . . , n -1.
For the boundary cell it is sufficient to multiply by two the energy of the elementary cell, whilst for the energy of a bulk modulus of length ℓ i and applied elongation ∆u i is four times the energy of elementary cell of length ℓ i /2 with an applied elongation ∆u = ∆u i /2. Hence, the total energy of the composite bar as a function of the segment lengths and elongations writes as:

P tot (ℓ; ∆u) = K(ℓ 1 ) ∆u 2 1 + K(ℓ n ) ∆u 2 n + n-1 ∑ i=2 1 2 K(ℓ i /2) ∆u 2 i . (3.10)
The optimal elongation of the segments to approach the solution of the problem (3.6) can be determined by minimizing P tot (ℓ; ∆u) under the constraint ∑ n i=1 ∆u i = ū. This gives the following Chapter 3. Analytical and numerical approaches implemented to study transverse cracking patterns in UD cross-ply laminates first-order optimality conditions, [see e.g. [START_REF] Luenberger | Linear and Nonlinear Programming[END_REF]:

∆u 1 = F 2 K(ℓ 1 ) , ∆u i = F K(ℓ i /2) for i = 2, . . . , n -1, ∆u n = F 2 K(ℓ n ) , n ∑ i=1 ∆u i = ū,
where F is the Lagrange multiplier representing the total reaction force acting on the bar to impose the end-displacement. After elementary algebraic manipulations, one finds the following final expression for the total elastic energy as a function of the crack spacing:

P tot (ℓ) = K tot (ℓ) ū2 2 = 1 S tot (ℓ) ū2 2 , (3.11)
where

S tot (ℓ) = S(ℓ 1 ) 2 + n-1 ∑ i=2 S(ℓ i /2) + S(ℓ n ) 2 , K tot (ℓ) = 1 S tot (ℓ) , (3.12)
are the equivalent total compliance and stiffness of the composite bar and

S(ℓ) = 1 K(ℓ) = ∆u 2 2P cell (ℓ) (3.13)
is the compliance of the elementary cell, the fundamental brick of our approach.

The optimal crack spacing problem

The optimal crack spacing is given by the solution of the problem (3.7), which reduces to the minimization of the elastic energy P tot (ℓ) under the constraint ∑ n i=1 ℓ i = L. Using the expression (3.11), this can be recast as the following maximization problem for the total elastic compliance (3.12):

max ℓ S tot (ℓ), ℓ ∈ R n , n ∑ i=1 ℓ i = L (3.14)
Eliminating the constraint by setting ℓ n = L -∑ n-1 i=1 ℓ i , gives an unconstrained maximisation problem in R n-1 . A solution {ℓ * i } n-1 i=1 must verify the following first order optimality conditions:

∂S tot (ℓ * 1 , ℓ * n-1 , L -∑ n-1 i=1 ℓ i ) ∂ℓ i = 0 n-1 i=1 ⇔ S ′ (ℓ * 1 ) = S ′ (L -∑ n-1 i=1 ℓ * i ), S ′ (ℓ * i /2) = S ′ (L -∑ n-1 i=1 ℓ * i ) for i = 2, . . . , n -1.
(3.15) Second-order optimality conditions are given by the sign of the (n -1) × (n -1) Hessian matrix H of second derivatives, which is calculated to get:

H = 1 2     H 1 + H n H n . . . H n H n H 2 + H n . . . H n . . . . . . . . . . . . H n H n . . . H n-1 + H n     , with      H 1 = S ′′ (ℓ * 1 ) H i = S ′′ (ℓ * i /2)/2, i = 2, . . . , n -1 H n = S ′′ (L -∑ n-1 1 ℓ * i ) (3.
16) The sign of the Hessian matrix determines the stability of the solution: if the Hessian matrix evaluated at a given solution is negative definite, the solution is stable, being a maximizer of the total x FIGURE 3.3: Composite bar alternative with displacement imposed only on the outer layer. In this case, differently from the case in Figure 3.1 there is no need for a distinction between boundary modules and bulk modules. All the modules are of the bulk type of Figure 3.1, see also Remark 1.

compliance, hence a minimizer of the total potential energy. Vice-versa, if it is not negative semidefinite, the solution is unstable. The Hessian matrix happens to be in the same form as the one of the chain of bi-stable springs studied by [START_REF] Puglisi | Mechanics of a discrete chain with bi-stable elements[END_REF]. We apply here the results of their stability analysis:

(S1) If H i < 0 for all i = 1, . . . , n, then the matrix is negative definite and the solution is stable.

(S2) If there are two or more than H i > 0, the matrix has at least one positive eigenvalue and the solution is unstable.

(S3) If only one of the H i is positive and n -1 are negative, the solution is stable if

r := n ∑ i=1 1 H i < 0 (3.17)
and unstable if the same quantity is r > 0.

A trivial solution of the problem (3.15) is the solution with uniform spacings, for which

ℓ 1 = ℓ n = ℓ * , ℓ i = 2ℓ * , with ℓ * = L 2 (n -1) . (3.18)
This solution exists independently of the property of the function S ′ (ℓ * ). ). To assess whether this is possible, we need first to compute the function S ′ (ℓ * ) and S ′′ (ℓ * ) and determine their qualitative properties.

Remark 1 A simpler variant of the problem is obtained by considering the case where only the displacement of the upper layer is imposed at the two ends, as shown in Figure 3.3. In that case, there are no special boundary modules and the total compliance (3.12) reduces to:

S tot (ℓ) = n ∑ i=1 S(ℓ i /2).
(3.19)

Chapter 3. Analytical and numerical approaches implemented to study transverse cracking patterns in UD cross-ply laminates

Elementary cell problem: computing the elastic compliance of the elementary cell and its derivatives

The modular approach defines the total elastic energy of a bar with several cracks as a function of the elastic energy of the elementary cell in Figure 3.1. The solution of the optimization problem (3.14) requires the accurate computation of the compliance (3.13) of the elementary cell S(ℓ) as a function of its length and of its first and second derivatives S ′ (ℓ), S ′′ (ℓ). We compute these functions through a domain-derivative approach and by solving the elementary cell problem through a classical finite element method.

To accurately evaluate the derivatives of S(ℓ) with respect to the domain length, we reformulate the elastic cell problem for a cell of length ℓ on a ℓ-independent domain Ω ≡ Ω c ∪ Ω o of length l by performing the following change of coordinates

x = (x, y) ∈ Ω ℓ → x = ( x, y) ∈ Ω with x = ℓ ℓ x, y = y, (3.20) for which ∂(•) ∂x = ℓ ℓ ∂(•) ∂ x , ∂(•) ∂y = ∂(•) ∂ y .
This change of coordinate defines a new energy density

W ℓ (u) = 1 2 ℓ ℓ A(x)ε(u(x)) • ε(u(x))
on the fixed domain Ω such that:

P ℓ (u) := Ω ℓ 1 2 A(x)ε(u(x)) • ε(u(x)) = Ω W ℓ (u) d x,
where

ε(u) = l ℓ ∂u x ∂ x , ∂u y ∂ y , ∂u x ∂ y + l ℓ ∂u y ∂ x T .
The solution of the following minimization problem on the ℓ-independent domain gives the elastic energy of the cell:

u ℓ := Argmin u∈ C P ℓ (u), P cell (ℓ) = P ℓ (u ℓ ), (3.21)
where the space of admissible displacements includes the Dirichlet boundary conditions

C = {u ∈ H 1 ( Ω) : u 1 (0, x 2 ∈ [0, H c + H o ]) = u 2 ( x 1 ∈ [0, ℓ], 0) = 0, u 1 ( ℓ, x 2 ∈ [H c , H c + H o ]) = ∆u}. (3.22)
The first order optimality condition for (3.21) gives the linear system to solve, a standard linear elastic problem on the unit cell:

DP ℓ (u ℓ )[v] = 0, ∀v ∈ C 0 . (3.23)
Here and henceforth, we denote by

DF (u)[v] := d dh F (u + h v) h=0 , D 2 F (u)[w, v] := d dh DF (u + h w)[v] h=0 (3.24)
3.3. Modular energetic approach for the case with n cracks 65 the first and second directional derivative of the functional F with respect to u in the directions v and w and by C 0 the space of admissible variations associated to C.

With the change of coordinates (3.20), the cell length ℓ appears directly as a parameter in the energy. Hence, the domain derivatives can be straightforwardly computed as

P ′ cell (ℓ) := dP cell (ℓ) dℓ = Ṗℓ (u ℓ ) + : 0 DP ℓ (u ℓ )[ uℓ ] = Ṗℓ (u ℓ ),
(3.25a)

P ′′ cell (ℓ) := d 2 P cell (ℓ) dℓ 2 = Pℓ (u ℓ ) + D Ṗℓ (u ℓ )[ uℓ ] (3.25b)
where

Ṗℓ (u) := Ω ∂ W ℓ ∂ℓ (u) d x, Pℓ (u) := Ω ∂ 2 W ℓ ∂ℓ 2 (u) d x.
and uℓ is derivative of the solution with respect to ℓ. The computation of the first derivative does not require the knowledge of uℓ because the second term in (3.25a) is equal to zero thanks to the stationarity condition (3.23). Vice-versa, the evaluation of the second derivative requires the sensitivity of the solution with respect to the cell length uℓ . Taking the derivative with respect to ℓ of the weak form of the equilibrium condition (3.23) gives the following linear tangent problem to solve for uℓ :

uℓ ∈ C 0 : D 2 P ℓ (u ℓ )[ uℓ , v] + D Ṗℓ (u ℓ )[v] = 0, ∀v ∈ C 0 .
(3.26)

Using these results, we evaluate S(ℓ), S ′ (ℓ), S ′′ (ℓ) for given material parameters and geometry with a standard finite-element technique as follows:

1. Mesh the fictitious domain with a unstructured triangular mesh, suitably refined around the singularity. A typical mesh is shown in Figure 3.5.

2. Solve the linear elastic problem (3.23) on the fictitious domain to get u ℓ for ∆u = 1 and evaluate S(ℓ) = 2/P cell (ℓ) and S ′ (ℓ) = -2 P ′ cell (ℓ)/P 2 cell (ℓ) with P cell (ℓ) = P ℓ (u ℓ ) and P ′ cell (ℓ) = Ṗℓ (u ℓ ).

3. Solve the linear tangent problem (3.26) on the fictitious domain to get uℓ and evaluate S ′′ (ℓ) = -2P ′ cell (ℓ) 2 /P 2 cell (ℓ) + 2P ′′ cell (ℓ)/P cell (ℓ) with P ′′ cell (ℓ) given by (3.25b). In practice the linear problems (3.23) and (3.26) are solved using the finite element library FEniCSx and the derivatives implied in their expressions are computed automatically exploiting the automatic symbolic differential capabilities of the UFL package.

Figure 3.5 reports the graph of P ′ cell (ℓ), S ′ (ℓ) and S ′′ (ℓ) as a function of the cell length ℓ for a bar made of isotropic materials with E o /E c = 1, H o /H c = 0.5, ν = 0.3, as for the example in Figure 3.2. These curves are the result of the interpolation of 100 data points corresponding to finite element simulations for uniformly spaced cell lengths ℓ ∈ [0.1, 10]. In our numerical experiments, we observed that the qualitative behavior of these functions do not depend on the specific value of the thicknesses and Young modulus ratio (see the following Section 3.5). Their key properties are: (P1) For ℓ → 0 one can neglect the influence of the core layer, getting the following limit elastic energy and compliance:

P 0 (ℓ) = E x o H o ℓ ∆u 2 2 , S 0 (ℓ) = ℓ (E x o H o ) . ( 3
P ∞ (ℓ) = E x c H c + E x o H o ℓ ∆u 2 2 , S ∞ (ℓ) = ℓ E x c H c + E x o H o (3.

28)

(P3) There exist a finite critical cell length ℓ c such that the derivative of the compliance, S ′ (ℓ), is monotonically decreasing for ℓ ∈ (0, ℓ c ), where S ′ (ℓ) < 0, and monotonically increasing for ℓ > ℓ c , where S ′ (ℓ) > 0. The cell length for which the compliance attains its minimal value is the solution of S ′′ (ℓ c ) = 0.

In our numerical experiments, we observed that these key properties stay true for a wide range of material and geometric parameters. We will discuss in Section 3.5 the dependence of ℓ c on the stiffness ratio and the thickness ratio.

Solutions for the crack spacing problem and their stability

The case of one crack dividing the bar in two segments

We revisit here at the light of our analysis the results anticipated in Section 3.2.3 for a bar divided in two sub-segments by a single crack. This is the simplest case, but also the most important. The previous works on composite cracking [see e.g. [START_REF] Nairn | Matrix microcracking in composites[END_REF][START_REF] Berthelot | Transverse cracking and delamination in cross-ply glass-fiber and carbon-fiber reinforced plastic laminates: Static and fatigue loading[END_REF] commonly assume that no more than one crack can appear simultaneously. In absence of imperfections, their n-crack solution is the result of an evolution process with a cascade of cracking events breaking the longest segment in two equal parts. We study the competition between the homogeneous solution with equal segment lengths ℓ 1 = ℓ 2 = L/2 and solutions with ℓ 1 ̸ = ℓ 2 .

For n = 2, the system (3.15) of first order optimality conditions for the crack spacing problem (3.7) reduces to the condition S ′ (ℓ 1 ) = S ′ (ℓ 2 ). When the properties (P1-P3) of Section 3.3.3 are verified, there can be at most two distinct cell lengths ℓ I and ℓ I I > ℓ I such that S ′ (ℓ I ) = S ′ (ℓ I I ). For each ℓ I I > ℓ c there exists a unique ℓ I (ℓ I I ) < ℓ c solution of S ′ (ℓ I ) = S ′ (ℓ I I ). Figure 3.6 illustrates the problem graphically (left) and reports the corresponding function ℓ I (ℓ I I ) (right). To compute it, we construct an interpolating function from the data of S ′ (ℓ) coming from the finite element computations of Section 3.3.3, we solve the nonlinear equation numerically for a set of given ℓ I I , and tabulate the interpolating function for ℓ I (ℓ I I ) once for all. When looking for ℓ I given ℓ I I , the numerical root-finding procedure converges robustly even in presence of the small scale oscillations of the finite element results, as those visible in the plot of Figure 3.6-left. Vice-versa, the complementary approach of finding the long segment given the short one leads to hill-conditioned problems and must be avoided.

As clear from the graphical illustration of Figure 3.6-left and the properties (P1-P3) of Section 3.3.3, the function ℓ I (ℓ I I ) monotonically decreases from ℓ c to the asymptotic value ℓ ∞ := ℓ I (∞).

To solve the crack spacing problem, it remains to determine the segment lengths ℓ I (L) and ℓ I I (L) as a function of the total bar length L such that L = ℓ I (L) + ℓ I I (L). Given ℓ I (ℓ I I ), this is easily obtained by solving numerically a root finding problem, which has a unique solution for ℓ I I > ℓ c . Finally, we get the diagram in Figure 3.7-left, showing a bifurcation from the fundamental homogeneous solution ℓ I = ℓ I I to the solution with ℓ I I > ℓ I . Stable solutions are in blue and unstable solutions in pink, where the stability is tested with the second-order condition of Section 3.3.2. For n = 2 the latter simplifies to S ′′ (ℓ I I ) < -S ′′ (ℓ I ). We can summarize the results for a single crack as follows:

• The homogeneous solution with equal segment lengths ℓ 1 = ℓ 2 = L/2 is stable for L < 2ℓ c and unstable for L > 2ℓ c . The corresponding total compliance of the bar is • For L > 2ℓ c there are two energetically equivalent stable solutions with ℓ 1 = ℓ I and ℓ 2 = ℓ I I or ℓ 1 = ℓ I I and ℓ 2 = ℓ I , where ℓ I is monotonically decreasing with L from ℓ c to ℓ ∞ . The total compliance of the bar for this solution is

S tot (L) = S(L/2). ( 3
S tot (L) = S(ℓ I (L)) + S(ℓ I I (L)) 2 .
(3.30)

For long bars (L ≫ ℓ c ), ℓ I → ℓ ∞ and ℓ I I → L -ℓ ∞ constitutes a useful asymptotic approximation of the solution.

Figure 3.7-right compares the total compliance of the two solutions, which is inversely proportional to the corresponding elastic energy: their difference is immaterial in absolute value for the considered numerical example. This is related to the fact that the change of concavity of S(ℓ) is quantitatively very small and implies a small margin of stability of the solutions. This point will be discussed further in the Section 3.5 dedicated to the study of the influence of the geometric and material parameters and in the conclusions.

The results obtained with the modular approach are in agreement with the preliminary analysis of Section 3.3.1 obtained via direct numerical simulations. In turn, we can generalize the modular approach to treat the case with several cracks without further numerical burden.

General results for n-1 cracks dividing the bar in n segments

When the properties (P1-P3) in Section 3.3.3 are verified, there can be at most two distinct cell lengths ℓ I and ℓ I I > ℓ I such that S ′ (ℓ I ) = S ′ (ℓ I I ). Moreover, this is possible only for ℓ I I > ℓ c , and S ′′ (ℓ I ) < 0, S ′ (ℓ I I ) > 0. This implies three fundamental consequences:

1. There can be at most two different elemental segment lengths ℓ I and ℓ I I co-existing in a solution of the optimality conditions (3.15).

2. Because of the stability criterion (S2) in Section 3.3.2, stable solutions can have at most one long elemental cell length ℓ I I .

Chapter 3. Analytical and numerical approaches implemented to study transverse cracking patterns in UD cross-ply laminates 3. Applying the stability condition (S3) in Section 3.3.2, solutions with n -1 short segments and one long segment are stable if

r := 1 S ′′ (ℓ * 1 ) + 1 S ′′ (ℓ * n ) + n-1 ∑ i=2 2 S ′′ (ℓ * i /2) < 0. (3.31)
We conclude that the possible solutions of the problem of Section 3.3.2 for a bar of length L divided in n segments by n -1 cracks can be classified as sketched in Figure 3.8:

(H) Homogeneous solutions with uniform cell lengths as in (3.18). These solutions are stable for

L < L c = 2 (n -1) ℓ c (3.32)
and the corresponding total compliance is

S (H) tot (ℓ) = (n -1) S L 2 (n -1) . ( 3

Effect of the material and geometric parameters

The stability of the solution with homogeneous crack spacing depends on the concavity property of the derivative of the compliance function of the basic cell S(ℓ). As shown in Figure 3.5, S(ℓ) is asymptotically linear for ℓ → 0 and ℓ → ∞. The change of concavity in the transition zone around ℓ = ℓ c is not clearly visible in the plot of S(ℓ). It is detected only by the change of sign of the slope of S ′ (ℓ) and the change of sign of S ′′ (ℓ) in the diagrams of Figure 3.5. The ratio S ′ (ℓ c )/S ′ (∞) between the minimal value of S ′ , attained at the critical length ℓ c , and its value for ℓ → ∞ is a possible quantitative signature of this change of concavity. The previous sections reported the results for the specific example of a laminated bar made of isotropic materials with a stiffness ratio E o /E c = 1 and a thickness ratio H o /H c = 0.5. In this case, the change of concavity is real, but quantitatively very small, as shown by the values of S ′ (ℓ c ) and S ′ (∞) in Figure 3.5. This implies that the margin of stability of the solutions is small, as discussed also in the comments to Figure 3.7. We investigate here how this change of concavity varies with the material and geometrical parameters. Figure 3.9 illustrates the results obtained when varying the thickness (left column) and stiffness ratios (right column) for the case of isotropic materials with ν = 0.3. They have been obtained with the numerical method described in Section 3.3.3. Each curve is the outcome of an interpolation of 100 finite element computations for different values of the bar length, keeping constant the other parameters. We managed to reliably investigate numerically only stiffness and thickness ratios of the order of 0.1 -10, the numerical problem becoming ill-conditioned otherwise. The small scale oscillations of the curves of S ′ and S ′′ for H o /H c = 0.1 and E o /E x = 0.1 are a marker of such numerical issues.

We observe the following facts:

• The ratio S ′ (ℓ c )/S ′ (∞) is larger for small stiffness ratios and small thickness ratios, i.e. soft and thin outer layers. For H o /H c = 0.1 and E o /E x = 0.1 this ratio is of the order of the unity.

• The value of ℓ c for which the minimum of S ′ is attained varies almost linearly with the thickness of the outer layer, while its dependence on the stiffness ratio is very weak.

Chapter 3. Analytical and numerical approaches implemented to study transverse cracking patterns in UD cross-ply laminates Figure 3.10 reports the curve of S ′ (ℓ) for the case of [0/90/0] composite laminates made of unidirectional fiber reinforced materials. We consider typical carbon (left) and glass (right) fiber composites with different thickness ratios, see the caption for the specific material parameters. The change of slope is either absent or barely visible for all the cases. As in the isotropic case, it is more evident for small thickness ratios. This can be explained by the high stiffness ratio between 0 and 90 laminae which is of the order 4 and 10 for the glass and carbon composites, respectively. Figure 3.11 illustrates the role of the shear stiffness for the case of glass fibers. The ratio S ′ (ℓ c )/S ′ (∞) increases with increasing shear stiffness. The small values of the shear stiffness of fiber-reinforced composites is a possible explanation of the difference with respect to the isotropic case of Figure 3.9.

Optimal number of cracks

Our work focuses on the study of the crack spacing for a given number of cracks. We quickly sketch here how to determine the optimal number of cracks as a function of the loading. For the sake of conciseness, we limit the discussion to the case of uniform crack spacing and the simplified boundary conditions in Figure 3.3-bottom. For this case, the total Griffith energy of the bar of total length L with n cracks writes as:

F n = 1 2 n S L 2n ū2 2 + G c (n -1) (3.37)
Figure 3.12-left reports the plots of F n for n = 1, . . . , 10 and L = 10. According to the global energy minimization principle, the cracked state is the one with minimal energy. The critical load ūn to pass from n -1 to n cracks is simply obtained by imposing F n = F n+1 . The grid lines and the black points in Figure 3.12 mark these critical loadings. Physically, the number of segments n is a discrete variable taking integer values. In a smeared model, one can also regard n as a real number, approximate indicator of the damage level in the core. In this setting, it is possible to get the optimal n by solving the equation dF n /dn = 0. This gives the red dashed curve in Figure 3.12-left, which is an excellent approximation of the lower envelope of the F n curves, and the plot in Figure 3.12-right for the crack density n/L as a function of the average stress in the bar σ average = ū/ (2 n S(L/2n)) . E c H c . The grid lines and the black dots denote the critical loadings for which F n = F n + 1. The dashed red curve is F n * with n * such that dF n /dn = 0, which gives an approximation of the lower envelope of the F n 's. Right: crack density as a function of the average stress in the bar corresponding to the red lower envelope in the left plot.

This kind of results have already been obtained in the literature. Previous works have shown them to be in agreement with experimental findings even with uniform crack spacing hypothesis [see e.g. [START_REF] Nairn | Matrix microcracking in composites[END_REF][START_REF] Berthelot | Transverse cracking and delamination in cross-ply glass-fiber and carbon-fiber reinforced plastic laminates: Static and fatigue loading[END_REF]. Figure 3.13 further provides a state diagram showing the optimal number of segments (hence cracks) as a function of the loading and the bar length.

The energy gap between solutions with uniform and non-uniform spacings is generally low for the numerical case tested in this work. Hence, considering non-uniform crack spacings does not significantly influence the results for the optimal number of cracks. For this reason, we do not report results for non-uniform spacings. The effect of the crack irreversibility is a further point that we do not address here in details. During a quasi-static evolution, not all the cracked states of Figures 3.12 are accessible, because previously created cracks cannot heal. Similar analysis are reported in León [START_REF] León Baldelli | Fracture and debonding of a thin film on a stiff substrate: Analytical and numerical solutions of a one-dimensional variational model[END_REF] for the case of a one-dimensional model of a bar-substrate system. The discussion is not duplicated here.

One-dimensional shear-lag model

In the previous sections we have calculated the elastic compliance of the single module of the system by a full-field finite element calculation. Several approximate solutions for the problem of a twolayered beam are available in the literature, see [START_REF] Nairn | A critical evaluation of theories for predicting microcracking in composite laminates[END_REF][START_REF] Hashin | Analysis of cracked laminates: A variational approach[END_REF][START_REF] Berthelot | Transverse cracking and delamination in cross-ply glass-fiber and carbon-fiber reinforced plastic laminates: Static and fatigue loading[END_REF]. The different theories are based on simplified hypothesis that are valid in a specific regime only. Here, we review only the simplest shear-lag model approach, which is also at the basis of simplified phase-field theories [START_REF] Alessi | Phase-field modelling of failure in hybrid laminates[END_REF]. We will show that the simple shear-lag approach is not capable to predict the change of concavity of the elastic compliance function, which in terms translates in the inability to predict non-uniform crack spacing solutions.

A simplified shear-lag approach

We develop a one dimensional model for the problem in Figure 3.1, where each layer is regarded as one-dimensional beam. The axial displacements of the two beams are elastically coupled together through distributed springs of stiffness k = G xy /H c , where G xy is the shear modulus of the bar, see Figure 3.14. The horizontal axis of the laminate undergoes a uni-axial displacement condition. We distinguish two axial displacements depending on the layer: u c (x) for the damageable core, and u o (x) in the unbreakable upper layer. The distributed springs are a simplified model of the coupling between the two layer due to shear. This kind of model is referred in the literature as a shear-lag approach. Following the same reasoning of the previous sections, we develop a variational approach à la Griffith and will explicit the analytical solution for the displacement field.

Elastic compliance with the shear-lag model

We apply the same assumptions of section 3.2 and define, within the Griffith fracture model the elastic energy functional of a cell to be minimized. Denoting by u c : x ∈ [0, ℓ] → R and u o : x ∈ [0, ℓ] → R the displacement fields, the total energy functional to be minimized writes as

P 1D cell (u c , u o ) := ℓ 0 1 2 E c H c u ′2 c (x)dx + ℓ 0 1 2 E o H o u ′2 o (x)dx + ℓ 0 1 2 G xy H c (u c (x) -u o (x)) 2 dx (3.38)
Chapter 3. Analytical and numerical approaches implemented to study transverse cracking patterns in UD cross-ply laminates

We introduce the following non-dimensional variables to obtain a non-dimensional version of equation (3.38):

x = ℓ x, u c (x) = ū ûc (x/ℓ), u o (x) = ū ûo (xℓ),
The elastic energy writes as

P 1D cell ( ûc , ûo ) := E c H c ℓ 1 0 û′2 c ( x) 2 + E o H o E c H c û′2 o ( x) 2 d x + G xy ℓ 2 E c H 2 c ( ûc (x) -ûo (x)) 2 2 dx (3.39)
Defining the dimensionless energy

P1D cell ( ûc , ûo ) := ℓ E c H c ū2 P 1D tot (u c , u o )
and the dimensionless parameters

ρ = E o H o E c H c , κ = G xy ℓ 2 E c H 2 c ,
representing the relative stiffness of the outer layer and of the interfacial springs, respectively, we finally get

P1D cell ( ûc , ûo ) := 1 0 û′2 c ( x) 2 + ρ û′2 o ( x) 2 d x + κ ( ûc (x) -ûo (x)) 2 2 dx (3.40)
Defining the space of admissible displacement field C ū which must respect the Dirichlet boundary conditions (3.41) the solution of the following minimization problem gives the elastic energy of the bar:

C ū = { ûc ( x), ûo ( x) ∈ H 1 ([0, 1]) : ûc (0) = 0, ûo (0) = 0, ûo (1) = ū},
ûc ( x), ûo ( x) := Argmin ûc ( x), ûo ( x)∈C ū P1D cell ( ûc , ûo ) (3.42)
The first order optimality condition for (3.42) gives the linear system to solve:

D P1D cell ( ûc , ûo )[v c , v o ] = 0, ∀(v c , v o ) ∈ C 0 ū. (3.43) where C 0 ū = { ûc ( x), ûo ( x) ∈ H 1 ([0, 1]) : ûc (0) = 0, ûo (0) = 0, ûo (1) = 0}. (3.44)
We obtain the following weak formulation from Equation (3.43):

1 0

[ û′ c ( x)v ′ c + ρ û′ o ( x)v ′ o + κ( ûc ( x) -ûo ( x))(v c -v o )] d x = 0 (3.45)
Integrating by parts gives:

1 0 [-ûc "( x)v c -ρ ûo "( x) v o ] d x + [ û′ c ( x)v c ] 1 0 + [ρ û′ o ( x)v o ] 1 0 + κ 1 0 [( ûc ( x) -ûo ( x))(v c -v o )] d x = 0.
(3.46)

The terms in brackets cancel out thanks to the Dirichlet boundary conditions except for the term in û′ c (1) which gives the Neumann boundary condition û′ c (1) = 0. Applying standard arguments of the calculus of variations, the strong formulation of Equation (3.46) writes:

ûc "( x) -κ( ûc ( x) -ûo ( x)) = 0, ρ ûo "( x) + κ( ûc ( x) -ûo ( x)) = 0 (3.47)
with the boundary conditions

ûo (0) = ûc (0) = ûo (1) = 0, û′ c (1) = 0,
Solving this system of linear equations, and replacing the solution in the (3.42) one gets the potential energy of the system. We proceed by calculating the second derivative with respect to ℓ of the corresponding elastic compliance S(ℓ) and we plot it as a function of the parameters in Figure 3.15. In these phase diagrams, we do not observe any change of sign of S ′′ (ℓ) with the material and geometrical parameters. This is due to the limitation of the reduced 1D approach that we overcome through the bi-dimensional approach.

Conclusion

We have studied transverse fracture in the core layer of symmetric laminate composites using a sharp interface Griffith model and a variational approach on both a reduced one-dimensional model and a more in-depth bi-dimensional model. We have shown that, even in a perfect system, the solution with uniform crack spacings, commonly accepted in the literature, is not always optimal. For sufficiently long bars, solutions with uniform spacings become unstable and stable solutions with a non-uniform crack spacings appear. With a modular model, we have related the bifurcation phenomenon to the concavity properties of the compliance S(ℓ) of the basic cell of the system as a function of the cell length ℓ, applying the ideas of [START_REF] Bourdin | The Variational Approach to Fracture[END_REF]. The change of concavity of the elementary cell compliance is a necessary condition for the existence of solutions with non-uniform crack spacings. We computed S(ℓ) and its first and second derivatives numerically by finite element method and an accurate domain derivative approach. Our results unveil a change of concavity of S ′ (ℓ), related to the instability of the solutions with uniform crack spacings. For the case of bar with n cracks, we have shown that two distinct cell lengths are possible, short and long, and that only solutions with a single long cell are stable.

Our results unveil a novel phenomenon in perfect symmetric laminated composites and describe its main qualitatively properties: the bifurcation from solutions with uniform crack spacings to non uniformly spaced solutions. However, the stability margin of the different solutions turns out to be quantitatively small. For usual material and geometric parameters, the dependency of the total energy of the laminate on the crack spacings is very weak, until the elementary cell length reaches a critical length of the order of the thickness. This implies that solutions with the same crack number and different crack spacings are almost equivalent from the energetic point of view. In particular, the change of concavity of the elemental compliance S is barely visible or absent for classical [0/90/0] laminated composites with unidirectional laminae made of fiber reinforced materials with large stiffness ratio and a relatively low shear stiffness. In this context random or deterministic imperfections will determine the crack spacings in real system, as supposed by previous models [see e.g. Silberschmidt, 2005, Vinogradov and[START_REF] Vinogradov | Probabilistic energy based model for prediction of transverse cracking in cross-ply laminates[END_REF]]. Yet, the results of our dedicated parametric analysis indicates that the change of concavity of S, related to energetic gap between solution with uniform and non-uniform crack spacings, becomes of the order of unit for laminated composites with thin and soft outer layers. We hope that our work can instigate further experimental works to understand whether the competition between uniform and non-uniform crack spacings can be deterministically linked to the geometrical and material properties in this specific regime, as predicted by our analysis. From the theoretical and numerical perspectives, multiple cracks in layered composites have been recently studied with cohesive Javaland et al. Phase-field approach to fracture of composite laminates

This chapter reports numerical simulations of the fracture mechanics problem of a sandwich beam using the phase-field approach. We focus on the crack spacing problem, trying to compare the numerical results obtained through the direct numerical phase-field simulation to the results of the Griffith approach introduced in Chapter 3. We start by introducing the main points of the variational approach to fracture starting from the Griffith approach and its regularized formulation leading to the gradient damage models applied nowadays as phase-field models of fracture. We give an example of its application to a classical problem, a bar under uniaxial tension, and the steps followed for the numerical implementation. We point out some issues of the classical phase field model when simulating a composite bar and how early non-physical delamination occurs. We discuss a possible strategy to reduce the delamination though the introduction of an anisotropic surface energy. Hence, we report the results of numerical simulations focusing on the crack spacing problem. We observe a transition from solution with uniform to non-uniform spacing in qualitative agreement with the results of the sharp interface approach of Chapter 3. Finally, we perform full-scale numerical simulations corresponding to experimental setting of Chapter 2 and discuss the influence of the stiffness ratio on the crack evolution.

4.1

The variational approach to fracture: from Griffith's theory to phasefield fracture models 4.1.1 The Griffith variational approach to fracture for a single crack on a given path [START_REF] Griffith | The phenomenon of rupture and flow in solids[END_REF] introduced a first approach to describe brittle fracture from an energetic viewpoint.

The theory was developed in linear elasticity considering an isotropic homogeneous material, and an energy based criterion for crack propagation was introduced for the first time. The main assumption of this model relies on the following concept: the energy spent to create a crack is proportional to the crack surface. The energy required to create a unit surface crack is called the material toughness G c , a new material constant. The propagation criterion translates then in to an energetic approach stating that a crack, present in an elastic isotropic homogeneous material, can propagate only if the released elastic energy can balance the surface energy required to create it.

To translate these ideas into a mathematical formulation, we start by defining the elastic and surface energies of a body in a bi-dimensional space. We start considering the case of a two dimensional body Ω ⊂ R 2 , where the potential crack path Γ is given in advance, and a displacement ū(t) is applied on a part of the boundary, denoted by ∂ u Ω, while the rest of the domain is free from any volume of boundary forces, see Figure 4.1. The position of the crack tip is tracked through the scalar variable l, the curvilinear abscissa along the crack path Γ, and we denote by l 0 the crack length before propagation. Griffith's original approach supposes the a priori knowledge of the crack path. In this case, describing crack propagation translates into a "time"-evolution problem, of a single crack. Therefore we denote variable l to be time dependant s.t. l → l(t). Here and henceforth, we denote by (•) ′ the derivative of a function with respect to its argument and by a dot (

•) the derivative with respect to t. Griffith's theory may be summed up into three main points:

1. Crack propagation criterion. An existing crack will propagate when the energy release rate G(l) := -P ′ (l) exceeds the critical threshold G c , commonly known as the material toughness:

P ′ t (l(t)) + G c ≥ 0.
2. Irreversibility condition. l(t) ≥ 0, Chapter 4. Phase-field approach to fracture of composite laminates which translates the physical condition of damage irreversibility.

3. Energy balance condition. The variation of total energy, sum of the elastic energy and the dissipated energy, is equal to the work of the external load. This implies that

P ′ t (l(t)) + G c l(t) = 0.
These three principles are the basis of the variational approach to fracture in the context of quasistatic rate-independent evolution problems.

From a numerical perspective, the problem must be formulated in the time-discrete settings. Considering the discretisation of the time (loading) variable in N steps {t i } N i=1 and denoting by l i = l(t i ), the three-item principle above become

l i -l i-1 ≥ 0, P ′ t i (l i ) + G c ≥ 0, P ′ t i (l i ) + G c (l i -l i-1 ) = 0. (4.1)
The set of inequalities (4.1) are the first order optimality conditions for the following optimization problem including the unilateral constraint due to the irreversibility condition l ≥ l i-1 :

l i = argmin l≥l i-1 P(l). (4.2)
Hence, in the time-discrete setting, the crack propagation problem is formulated as a simple energy minimisation principle. This principle will be the basis for the numerical developments.

The variational approach to fracture: Griffith revisited

The variational approach of brittle fracture was introduced by Francfort and [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF], extending to the case with arbitrary crack patterns the Griffith approach of section 4.1.1. It seeks to model failure through a discrete approach where cracks are modeled as surface discontinuities of the displacement field that dissipate an energy proportional to their surface.

Let us consider the same notation defined in Figure 4.1. The crack is designated through the variable Γ. The associated surface energy can be written as

S(Γ) = Γ G c dH N-1 , (4.3)
where dH N-1 is the N -1 Hausdorff measure, which is a surface measure [C. [START_REF] Evans | Measure theory and fine properties of functions[END_REF], N being the dimension of the space in which the domain exists. We assume that the body is linear elastic before cracking and we neglect geometrical non-linearities. We denote by u :

x ∈ Ω → u(x) ∈ R 2 the displacement field, ε(u) = 1 2 (∇u + ∇ T u)
the linearized strain, a second order tensor, and by A 0 the fourth order elastic stiffness tensor. Assuming that the displacement ū is prescribed on a part of the boundary and the absence of body and surface loads, the potential energy of the system is given by:

P (u, Γ) := Ω\Γ 1 2 A 0 (x)ε(u(x)) • ε(u(x)) dx, (4.4)
The total energy according to the Griffith's viewpoint can then be seen as the contribution of an elastic and a dissipation term. The dissipation term is thus proportional to the crack surface and the total energy writes:

E (u, Γ) := P (u, Γ) + S(Γ) (4.5)
We define the admissible displacements field as the sufficiently regular fields u ∈ H 1 (Ω\Γ) that satisfy the kinematic boundary condition such that u = ū on ∂ u Ω, where H 1 (Ω \ Γ) denotes standard Sobolev space of function with square integrable first derivatives on (Ω \ Γ). We denote by C t (Γ) the corresponding affine space:

C t (Γ) ≡ {u ∈ H 1 (Ω \ Γ), u = ū(t) on ∂ u Ω} (4.6)
Given an initial state (u 0 , Γ), the imposed displacement ū on the boundary ∂ u Ω, the Griffith revisited variational quasi-static rate-independent evolution re-frames the three pillars in the following statements, where the pair (u t , Γ t ) must verify at every time t:

(ir) Irreversibly: ∀τ ∈ [0, t) : Γ τ ⊆ Γ t (gst) Unilateral global stability: E(u t , Γ t ) ≤ E( û, Γ), ∀( û, Γ) such that Γ ⊇ Γ t , û ∈ C t ( Γ) (eb) Energy balance: E(u t , Γ t ) = E(u 0 , Γ 0 ) + t 0 ∂Ω\Γ(τ) [A 0 ε(u τ )ν] • uτ dH N-1 dτ,
where ν is the outer normal to Ω.

As opposed to the classical Griffith approach, the crack path is a genuine unknown of the problem.

In the time-discrete variational formulation, the cracked state of the body (u i , Γ i ) at time t i is found by solving the unilateral minimization problem for each time iteration:

(u i , Γ i ) = Argmin u∈C t (Γ i ), Γ⊇Γ i-1 E(u, Γ). (4.7)
under the irreversibility condition with respect to the crack set Γ i-1 at the previous time step.

The problem (4.7) is extremely complex to solve directly because the displacements jump and therefore the crack sets are part of the unknowns. The original mathematical theory associated to the variational approach to fracture can be found in Ennio De Giorgi's formulation of free discontinuity problems [De Giorgi [START_REF] Giorgi | Un nuovo tipo di funzionale del calcolo delle variazioni[END_REF]. Similar problems arise of in completely unrelated contexts, such as for signal reconstruction, curvature problems, and what is most known as the Mumford and Shah functional for image reconstruction [START_REF] Mumford | Optimal approximation by piecewise smooth function and associated variational problems[END_REF]].

Regularized formulation

In the context of image segmentation, [START_REF] Ambrosio | On the approximation of functionals depending on jumps by elliptic functionals via Gamma-convergence[END_REF] proposed a regularized version of the free discontinuity problem for the Mumford and Shah functional. The same approach has been adapted to the fracture mechanics problem (4.7) by Bourdin et al. [2000a], proposing an original solution strategy for the fracture mechanics problem. In this approach sharp discontinuities are smeared Chapter 4. Phase-field approach to fracture of composite laminates across bands of non-vanishing width [START_REF] Braides | Approximation of Free-Discontinuity Problems[END_REF] introducing a smooth scalar field α representing the damage of the material. The regularized energy functional writes as:

E η (u, α) = Ω (1 -α) 2 + κ η 2 A 0 ε(u) • ε(u) dx Elastic Energy + G c Ω α 2 η + η|∇α| 2 dx Dissipated Energy , (4.8)
where κ η = o(η) is a residual stiffness term. In this context, the analogue of the quasi-static problem (4.7) writes as

(u i , α i ) = Argmin u∈C t , α∈D(α i-1 ) E (u, α), ( 4.9)
where the spaces C t and D(β) are defined such that the damage field respects the irreversibility condition and vanishes where the displacement in imposed:

C t ≡ {u ∈ H 1 (Ω), u = ū(t) on ∂ u Ω}, (4.10a) D(β) ≡ {α ∈ H 1 (Ω), α ≥ β, α = 0 on ∂ u Ω}. (4.10b)
The functional (4.8) introduces an additional length parameter η. Gamma convergence results show that, for η → 0, the global minimizers of the regularised energy (4.11) tend toward the global minimizers of the Griffith energy (4.5) [Chambolle, 2004, Ambrosio and[START_REF] Ambrosio | Approximation of quasi-static evolution of brittle fractures[END_REF]].

Gradient damage models

The functional (4.8) can be written in a more general form as follows [START_REF] Benallal | Bifurcation and stability issues in gradient theories with softening[END_REF], Marigo and Pham, 2010, Pham and Marigo, 2010[START_REF] Marigo | An overview of the modelling of fracture by gradient damage models[END_REF]]

E (u, α) = Ω 1 2 a(α)A 0 ε(u) • ε(u) dx Elastic Energy + Ω w 1 w(α) + η 2 |∇α| 2 dx Dissipated Energy , (4.11)
where a(α) is a stiffness modulation function through the damage variable α and the function w(α) is a local dissipation term. Damage models assume that a(α) monotonically decreases, reducing the elastic stiffness, and the energy dissipation function w(α) increases monotonically.

In the variational formulation, a rate-independent evolution of the displacement and damage fields (u t , α t ) pair must respect the following conditions [START_REF] Pham | Approche variationnelle de l'endommagement : Ii. les modèles à gradient[END_REF][START_REF] Pham | The issues of the uniqueness and the stability of the homogeneous response in uniaxial tests with gradient damage models[END_REF][START_REF] Marigo | An overview of the modelling of fracture by gradient damage models[END_REF]:

• (ir) Irreversibility: the damage variable α t is a non-decreasing function of time: αt (x) ≥ 0.

(4.12)

• (uls) Unilateral Local Stability: ∀(v, β) ∈ C t × D(α t ) ∃ h ≥ 0 : ∀h ≤ h, E t (u t + h(v -u t ) + h(β -α t )) ≥ E t (u t , α t ). (4.13)
• (eb) Energy Balance:

E t (u t , α t ) -E 0 (u 0 , α 0 ) = t 0 ∂Ω (σ τ ν) • uτ ds dτ (4.14)
where we introduced the stress tensor

σ t = a(α t )A 0 ε(u t ).
Developing at the first order the unilateral local stability condition (4.13), one finds that the solution should respect the mechanical equilibrium and damage criterion:

• (eq) Equilibrium:

divσ t = 0 in Ω σ t n = 0 on ∂Ω\∂ u Ω (4.15)
• (dc) Damage Criterion defined by the Kuhn-Tucker conditions:

αt ≥ 0, g ε (u t ) , α t , ∇ 2 α t ≥ 0, g ε (u t ) , α, ∇ 2 α t αt = 0 in Ω, (4.16a) αt ≥ 0, ∂α t ∂n ≥ 0, ∂α t ∂n αt = 0 on ∂Ω\∂ u Ω (4.16b) where g ε, α, ∇ 2 α = a ′ (α) 2 A 0 ε • ε + w ′ (α) -2 η 2 ∇ 2 α w 1 .
Based on the behaviour the materials exhibit at complete failure, the damage models can be divided in two categories [START_REF] Marigo | An overview of the modelling of fracture by gradient damage models[END_REF][START_REF] Pham | The issues of the uniqueness and the stability of the homogeneous response in uniaxial tests with gradient damage models[END_REF]]:

• Strongly brittle materials: the dissipated energy is finite when α evolves from 0 to 1, hence w(1) < ∞.

• Weakly brittle materials: the dissipated energy at complete failure is infinite.

To understand the role of the functions a(α) and w(α), one can study first the homogeneous response, where α and ε are assumed to be constant in space. The domains of admissible strains and stresses corresponding to the damage criterion (4.16b) are:

A ε (α) = ε : A 0 ε • ε ≤ - 2w 1 w ′ (α) a ′ (α) , A σ (α) = σ : A -1 0 σ • σ ≤ 2w 1 w ′ (α) s ′ (α) (4.17)
Based on the growth properties of these domains the material is said to be:

• strain hardening (resp. softening) if A ε increases (resp. decreases) as α increases (resp. decreases),

• stress hardening (resp. softening) if A σ increases (resp. decreases) as α increases (resp. decreases).

Strongly brittle materials with stress-softening are required to obtain an energetic equivalence between the Griffith brittle fracture and a gradient damage model [START_REF] Marigo | An overview of the modelling of fracture by gradient damage models[END_REF]. An example of the homogeneous and localized solutions is given in the following section on a one-dimensional bar subject to uni-axial tension.

One-dimensional bar modelled through a gradient damage model

The case of a one-dimensional bar under tension is the most commonly studied example used to disclose the properties of the damage gradient models. We consider a bar of total length L with imposed end-displacements: u(0), u(L) = t L, as shown in Figure 4.2.

In this case, the functional (4.11) simplifies to

E (u, α) = L 0 a(α)E 0 u ′ (x) 2 2 + w 1 w(α) + η 2 α ′ (x) 2 dx, (4.18)
where E 0 denotes the Young modulus of the undamaged material, the stress being σ = a(α)E 0 u ′ (x).

Introducing the compliance function s = 1/a(α) the damage criterion writes as:

g ε, α, α ′′ := - s ′ (α) 2E 0 σ 2 + w ′ (α) -2 η 2 α ′′ w 1 ≥ 0. (4.19)
Two classical choices for the constitutive functions w(α), a(α) are

• (AT1): w(α) = α, a(α) = (1 -α) 2 , • (AT2): w(α) = α 2 , a(α) = (1 -α) 2 .
The basic properties of the damage models can be unveiled by looking at two fundamental damaging modes:

• homogeneous solutions, with damage uniform in space, and

• localised solutions, with a single damage localisation reaching the fully damaged state α = 1 on a point of the domain.

For homogeneous solutions, the damage criterion (4.16b) writes as:

σ 2 ≤ 2 E 0 w 1 w ′ (α) s ′ (α)
Hence, the maximum allowable stress as a function of the damage level is

σ c (α) := 2 E 0 w 1 w ′ (α) s ′ (α) (4.20)
The elastic limit for an undamaged bar (α = 0) is

σ e = w 1 E 0 2w ′ (0) s ′ (0) . (4.21)
Moreover, it is useful to define the critical stress σ M for which the response pass from the stresshardening (σ ′ c (α) > 0) to the stress-softening (σ ′ c (α) < 0) regimes. This happens either at the elastic limit, where the response may be not smooth, or at a critical damage level α M such that σ ′ c (α M ) = 0. Figure 4.3, taken from [START_REF] Pham | The issues of the uniqueness and the stability of the homogeneous response in uniaxial tests with gradient damage models[END_REF], shows the homogeneous stress-strain responses of the AT1 and AT2 models and their corresponding stability diagrams [START_REF] Pham | The issues of the uniqueness and the stability of the homogeneous response in uniaxial tests with gradient damage models[END_REF]. The stability of the solution is assessed by studying the sign of the second derivative of the energy functional in a suitable functional space accounting for the damage irreversibility [START_REF] Pham | The issues of the uniqueness and the stability of the homogeneous response in uniaxial tests with gradient damage models[END_REF][START_REF] Marigo | An overview of the modelling of fracture by gradient damage models[END_REF].

For the AT1 model

σ AT 1 e = σ AT 1 M = w 1 E 0 > 0.
For the AT2 model, since w ′ (0) = 0, there is no elastic phase and the damage is non-null for any loadings t > 0; the critical stress for the passage to the stress softening regime is

σ AT 2 M = 3 8 3 2 .
From the stability diagrams it appears that the length scale η introduces a size effect. In fact the stability limit depends on the ratio between the bar length and the internal length L/η. For bars short enough (L ≪ η), the solution in the stress-softening regime can be stable. However, for sufficiently long bars (L ≫ η) the solution is always unstable in the stress softening regime.

The stability loss of the homogeneous solution will corresponds to the appearance of a localised solution. In the fully localised solution, reaching the damage level α = 1, the stress vanishes: σ = 0. We can understand this solution as the regularized approximation of a crack within the damage model. The corresponding damage profile can be readily calculated by looking for solution satisfying the damage criterion (4.19) as an equality at least on a part of the domain and integrating the corresponding differential equation:

- s ′ (α(x)) 2E 0 w 1 σ 2 + w ′ (α(x)) -η 2 α ′′ (x) = 0 (4.22)
where the stress σ is constant in space because of the equilibrium condition.

For the case of the AT1 model, considering that the damage evolves only in one segment of the bar, the damage profile and the length of the segment with non-vanishing damage are found by solving the non-linear boundary value problem (4.22) with the boundary condition α = α ′ = 0. This solution is shown in Figure 4.4 which is made of two symmetric parabolic profiles.

In general, for this crack-like solution, the length of the damaged zone D is calculated to be

D loc = 2η 1 0 1 w(α) dα. ( 4

.23)

A further fundamental quantity to identify the parameters of the damage model is the energy dissipated in a crack-like solution. We assimilate this energy to the fracture toughness. By integrating the dissipated energy in (4.18) using the optimal damage profile solution of (4.23) one finds that the 4.1. The variational approach to fracture 89 dissipated energy in a fully localised solution is [START_REF] Marigo | An overview of the modelling of fracture by gradient damage models[END_REF]:

G c = 4 1 0 w(α)dα η. (4.24)
We will identify the material strength as the critical stress for the stability loss of the homogeneous solution and denote it by σ c . Assuming that the internal length η is small enough, we will consider that it coincide with the stress σ M for the passage to the stress-softening regime (this is not true for short bars as shown in the stability diagrams of Figure 4.3).

To conclude, two are the main material parameters that can describe brittle fracture: the fracture toughness G c and its critical stress σ c . We resume in the Table 4.1 resumes the expressions of the key properties of the damage models used in this work on the material parameters introduced in the total energy functional (4.11). We add to the classical AT 1 and AT 2 models a further parametric model introduced in [START_REF] Alessi | Gradient damage models coupled with plasticity and nucleation of cohesive cracks[END_REF] that we denote as AT m . This model is defined by the constitutive functions:

w(α) = 1 -(1 -α) 2 , a(α) = 1 -w(α) 1 + (m -1) w(α) (4.25)
Figure 4.5 reports its homogeneous response that it is characterised by a linear softening behaviour. This latter model is interesting because it allows for more flexibility through the parameter m. In the case of the AT 1 and AT 2 models, the internal length η is determined directly from the values of the so-called material characteristic length L c = G c E 0 σ c 2 [START_REF] Tanné | Crack nucleation in variational phase-field models of brittle fracture[END_REF]. For the AT m model we get

G c E 0 σ c 2 = mηπ 2 .
For example, this can allow the user to vary η by keeping σ c and G c fixed, by adjusting m [START_REF] Wu | A unified phase-field theory for the mechanics of damage and quasi-brittle failure[END_REF]. This flexibility has certain limits, which are not discussed here. In particular, one should assure m > 1 to keep the strain-hardening property of the homogeneous response.

w(α) a(α) c w σ e √ w 1 E 0 σ c √ w 1 E 0 G c w 1 η G c E 0 σ 2 c D loc η AT 1 α (1 -α) 2 2 3 1 1 8 3 8 3 η 4 AT 2 α 2 (1 -α) 2 1 2 0 3 8 3 2 2 256 27 η ∞ AT m 1 -(1 -α) 2 1-w(α) 1+(m-1) w(α) 1 2 2 m 2 m π mηπ 2 π TABLE 4
.1: Key properties of the gradient damage models considered in this work as a function of the material parameters E 0 , w 1 and η appearing in (4.11): the elastic limit σ e , the material strength σ c , the fracture toughness G c , and the width of the damage band D loc .

Numerical implementation of the gradient damage model

Alternate minimization

The numerical implementation of the regularized approach consists in looking for solutions for the problem in Equation (4.9). The method used in this work relies on the iterative algorithm introduced in Bourdin et al. [2000a]. The main idea of this numerical approach is based on the alternate minimization of the energy functional. Alternate minimization consists in splitting the minimization in problem in different sub-problems with respect to each argument that needs to be minimized. Problem (4.9) is solved in two steps:

1. Elastic problem: the damage field is fixed and the functional (4.11) is minimised with respect to the displacement u.

2. Damage problem: the displacement field is fixed and the functional (4.11) is minimised with respect to the damage variable α under the irreversibility constraint α ≥ α i-1 .

At each time steps, the two problems above are solved iteratively until a convergence criterion is met. Algorithm 1 summarizes the alternate minimization algorithm.

Algorithm 1: Alternate Minimization (u old , α old ) = (u i-1 , α i-1 ) (Initialize) 1. u new = argmin u E η (u, α old ) (linear solver) 2. α new = argmin α≥α i-1 E η (u new , α) (bound constrained linear solver) 3. if ∥α new -α old ∥ ∞ > tol repeat 1-3 (convergence test) 4. (u i , α i ) = (u new , α new ) (update)
Chapter 4. Phase-field approach to fracture of composite laminates

U t U e U G Energies G c W

Setting of the gradient damage problem

The elastic energy density of the undamaged materials composing the core (c) and the outer layer (o) are defined by

W 0 c (ε) = 1 2 A c ε • ε, W 0 o (ε) = 1 2 A o ε • ε.
(u i , α i ) = Argmin u∈C t , α∈D(α i-1 ) E (u, α), (4.27)
Assuming that the damage is non null only in the core layer, the total energy functional for the composite bar can be written as

E (u, α) = Ω c a(α)W 0 c (ε(u))dx + Ω o W 0 o (ε(u))dx + Ω c w 1 w(α) + η 2 |∇α| 2 dx
where ε 0 is an inelastic strain field (null in the original problem). The functional should be minimised imposing the following Dirichlet conditions:

• A symmetry condition on the bottom, u y (x, 0) = 0

• Null axial displacement is null at left u x (0, y) = 0

• Imposed axial displacement at right u x (L, y) = tL

• Null damage at the left and right ends and on the whole top layer.

To simplify the implementation and the visualization, instead of imposing the end-displacement, we impose the loading through an inelastic strain ε 0 = -te 1 ⊗ e 1 , by minimizing the modified functional

E ( ũ, α) = Ω c a(α)W 0 c (ε( ũ) -ε 0 )dx + Ω o W 0 o (ε( ũ) -ε 0 )dx + Ω c w 1 w(α) + η 2 |∇α| 2 dx
under the homogeneous version of the Dirichlet boundary conditions on ũ, ũx (L, y) = 0. The equivalence of the two loading devices is given by the change of variable

ũ = u -t x 1 ⇒ ε(u) = ε( ũ) + t e 1 ⊗ e 1
For visualization purposes, we will not include the constant term W 0 o (ε 0 ) in the plots reporting the numerical value of the elastic energy.

Key parameters

We use non-dimensional parameters, taking as a scaling length the thickness of the core layer H c and as scaling stress its Young modulus E c . The main parameters affecting the simulations are

L c = G c E c σ 2 c .
In practice, we set in our code H c = 1, E c = 1 and vary the value of H o , E o , L. In the AT1 and AT2 models, the length-scale L c is controlled by the internal length η, through the relations in Table 4.1.

Depending on the situations, we set either G c = 1 or σ c = 1 (this can be done without loss of generality, corresponding to a suitable re-scaling of the displacement field [START_REF] Tanné | Crack nucleation in variational phase-field models of brittle fracture[END_REF]).

We start by introducing a simple AT1 approach to solve problem (4.27) applied to the geometry of Figure 4.7. The issues encountered with this approach are discussed along with a solution developed in this work in the next section. An isotropic approach is discussed at first and an orthotropic application is tackled in a second moment. Results obtained through local minimization techniques are compared to the results of global minimization in section 4.6. A stability analysis of the results is also discussed. Finally, we give a conclusion on this approach and its comparison with global minimization.

Numerical simulations on a composite bar through the AT1 model

In this section, we apply the AT1 model to the composite bar in Figure 4.7, solving problem in Equation 4.27. We discuss here few examples to illustrate how different thickness ratios and internal lengths influence the damage field α and the appearance of transverse cracks. Two thickness ratios are studied and shown as the main examples showing the damage field evolution in the composite bar.

We consider the thickness ratios H o /H c = [0.25, 0.5] and different internal lengths. The simulations are for G c = 1 and with E c = E o = 1 so that G c is fixed and σ c evolves. We compute also simulations where σ c is fixed and G c evolves through η, in order to compare the evolution of the damage field between the two possible configurations. 10 reports different snapshots of the damage field and the evolution of the elastic and dissipated energy as a function of the loading.

For the thicker stacking sequence, and for a relative large internal length η = 0.5, we observe at the critical stress σ c = 3 G c E c /8η, the development of one distinct transverse crack, see Figure 4.8, on the left for G c fixed at 1 and on the right σ c is fixed at 1. For larger loadings, the damage patterns are less obvious, as shown in the figure with the snapshot at the last applied loading (σ c = 2.5). The energy curves, in both cases, in Figure 4.8 show the evolution of the dissipated and elastic energy, the jump occurs for a loading step of approximately σ c = 0.86 (left), this corresponds to the first nucleation and at σ c = 1. on the right figure since this parameter is fixed. As the loadings increase, the energy curves do not appear to jump abruptly, rather they increase (dissipated energy) and decrease (elastic energy) smoothly. This corresponds to the stable evolution of diffused damage. Figure 4.8 shows the damage field Different internal lengths have been simulated showing changes occurring in the cracking pattern as well as the start of delamination at the interface between the outer and inner layers. In Figure 4.9 the same configuration is tested for an internal length η = 0.1. When keeping G c fixed to 1, the internal length modulates the critical stress σ c . The image shows a different behaviour from the one visible with η = 0.5. Transverse cracks form at the extremities of the bar and quickly propagate at the interface. With this smaller internal length, we observe three main changes. Firstly, the position at which cracks appear seems to be influenced by the internal length value, however they maintain Figure 4.10 shows the case of a composite bar with thickness ratios: H o /H c = 0.5. Increasing the thickness ratios leads to the same issue, as an early delamination occurs at the interface between the two differently oriented layers. The first transverse cracks appear symmetrically one to the other on each extremity. As the crack approaches the interface, it develops following its horizontal direction creating a delamination type of behaviour. As the thickness ratio H c /H o increases, the crack continues to develop along this line preventing other types of damage to occur. Increasing the thickness of the core, the total number of transverse cracks is reduced and the delamination occurs for smaller loadings.

The delamination can be occasionally observed in physical experiments, but only at larger loadings, see Chapter 1 and Chapter 2. In our simulation, delamination is a parasite phenomena. It occurs for loading levels that are too small and it prevents us from correctly simulating the transverse cracking phenomenon.

Penalizing the delamination through an anisotropic toughness

The previous section illustrates the negative role of the parasite effect of delamination on the simulation of transverse cracking. The observed issue is the inability to withstand multiple transverse

AT1 E o = 1 E c = 1 ν = 0.3 Hc = 1 H o = 0.5, ψ = 1, η = 0.1 FIGURE 4.10: H o = 0.5, H c = 1, E o = E c = 1.0, η = 0.25
cracks before physical delamination, instead the transverse cracks nucleate and rapidly delaminate at the interface. The basic issue is probably related to the difficulty to accurately tune the phase-field nucleation criterion under multi-axial loadings. This issue has been discussed in a recent paper by [START_REF] Lorenzis | Nucleation under multi-axial loading in variational phase-field models of brittle fracture[END_REF]. In our situation, the delamination crack develops under a predominant shear loading. The classical phase-field model does not allow the user to independently set the critical shear-stress τ c and normal stress σ c for crack propagation. The τ c /σ c is fixed by the isotropy of the damage model and is too small with respect to realistic values. A parallel work is ongoing in our group to tackle this issue but a satisfactory solution is not yet available.

Here we attempt to develop an alternative method to avoid delamination, consisting in introducing an anisotropic fracture toughness, to penalise axial cracks corresponding to delamination.

Anisotropic surface energy

The energy of a sharp-interface Griffith model with an anisotropic surface energy writes as

E (u, Γ) = E d (u, Γ) + E s (Γ) = Ω\Γ W(u)dΩ + G 0 Γ γ(n)dH n-1 (4.28)
where γ(n) denotes the unit normal to the crack Γ and γ a dimensionless function giving the dependence of the fracture toughness on the orientation of the crack. In 2d, the fracture toughness for a crack propagating in the direction determined by the angle θ, see Figure 4.11, is G c (θ) = G 0 γ(n(θ)).

A simple way to introduce a weakly anisotropic surface energy [START_REF] Li | Crack kinking in a variational phase-field model of brittle fracture with strongly anisotropic surface energy[END_REF] in the phase field model is to introduce a gradient damage model in the form

E η (u, α) = Ω a(α)W(u)dΩ + w 1 Ω w(α) + η 2 ∂α ∂x 2 + (1 + ψ) ∂α ∂y 2 dΩ, (4.29)
where ψ introduces an anisotropic toughness in the gradient regularization. The effect of this term is to let the effective internal length depend on the direction:

• The effective internal length is η for a damage localisation with ∂α ∂y = 0, i.e. a transverse crack in our beam geometry.

• The effective internal length is η √ ψ for a damage localisation with ∂α ∂x = 0, i.e. an axial crack in our beam geometry.

For a given specific fracture energy w 1 this introduces in turns a dependence of the fracture toughness on the direction, which can be immediately deduced from the formulas in Table 4.1. Using the fact that for one-dimensional damage localizations G c ∝ w 1 η, we can deduce that the ratio between the axial and transverse fracture toughness introduced is exactly √ ψ. For the AT 1 models, the fracture energy for axial and transverse cracks is

G c (axial crack) = 3 8 w 1 η ψ, G c (transverse crack) = 3 8 w 1 η.
Similar results are obtained for the other models. Hence, increasing ψ it is possible to penalise axial cracks by increasing the fracture toughness in the axial direction.

Chapter 4. Phase-field approach to fracture of composite laminates

Influence of the anisotropic surface energy

Figures 4.12-4.13 report the numerical results obtained for different values of the anisotropic parameter ψ in the surface energy (4.29), reproducing the same test case as in Figure 4.10. The simulations are for the following parameters: total bar length L = 10, thickness ratio: H o /H c = 0.5, internal length η = 0.25, and the mesh size h is adjusted so that h = η/10. The main observations are the following:

• For relatively small values of ψ (ψ = {1, 5}), the critical load for the first crack nucleation and the corresponding initial crack pattern do not vary notably. However, the delamination phenomenon occurs similarly as in the isotropic case when ψ is too small.

• For ψ = 5 we observe a complex and unexpected array of cracks for large loadings (t = 5).

• Large values of the anisotropy prevents delamination to occur and we observe a third crack appearing between the first two at larger loading. However, this comes a phase of diffuse damage that significantly affects the solution, delaying the nucleation load, as shown in Figure 4.13d-4.13f.

The anisotropy penalizes variations of the damage in the y directions and adds a term to the surface energy coming from the crack tip that is not accounted for in the simplified one-dimensional analysis of Section 4.1.5. We speculate that the diffuse damage state and the delay of the initial crack initiation loading is related to the additional surface energy present at the tip of vertical cracks, where the damage field must go from 1 (crack tip) to 0 (undamaged upper layer).

Looking at the numerical results, we can deduce an intermediate value of the anisotropic parameter ψ ≃ 10 to be preferable to perform numerical simulations. The introduction of an anisotropic toughness allows to penalize damage along the y axis and therefore delay delamination at the interface. However, it can significantly perturb the damage evolution in a way that is only partially controlled and understood. Hence, we conclude the anisotropy cannot be regarded as a robust solution to delay the onset of the crack nucleation. For this reason, we will continue to use the isotropic model in the rest of this Chapter.

AT1 E o = 1 E c = 1 ν = 0.3 Hc = 1 H o = 0.5, ψ = 1, η = 0.1 E o = 1 E c = 1 ν = 0.3 Hc = 1 H o = 0.5, ψ = 5, η = 0.1 (A) ψ = 1 AT1 E o = 1 E c = 1 ν = 0.3 Hc = 1 H o = 0.5, ψ = 1, η = 0.1 E o = 1 E c = 1 ν = 0.3 Hc = 1 H o = 0.5, ψ = 5, η = 0.1 (B) ψ = 5 E o = 1 E c = 1 ν = 0.3 Hc = 1 H o = 0.5, ψ = 10, η = 0.1 E o = 1 E c = 1 ν = 0.3 Hc = 1 H o = 0.5, ψ = 20, η = 0.1 (C) ψ = 10 E o = 1 E c = 1 ν = 0.3 Hc = 1 H o = 0.5, ψ = 50, η = 0.1 E o = 1 E c = 1 ν = 0.3 Hc = 1 H o = 0.5, ψ = 100, η = 0.1 (D) ψ = 20 E o = 1 E c = 1 ν = 0.3 Hc = 1 H o = 0.5, ψ = 50, η = 0.1 E o = 1 E c = 1 ν = 0.3 Hc = 1 H o = 0.5, ψ = 100, η = 0.1 (E) ψ = 50 (F) ψ = 100

Analysis of crack-spacings obtained through phase-field simulations

In this section, we discuss the crack morphology and spacing at the first nucleation event. The main goal is to compare the results of the phase field approach with the analysis of Chapter 3, based on the sharp interface Griffith model and global energy minimality arguments. We perform numerical simulations for isotropic layers with thickness ratio H o /H c = 0.5, stiffness ratio E o /E c = 1, and ν c = ν o = 0.3, as for the case in Figure 4.20, setting E c = 1 and H c = 1 to define the scale of the stress and the lengths. In the phase field model, we perform a set of simulations varying the internal length η and the total length of the bar L. We adopt the AT 1 model and set σ c . When varying η the toughness varies with the relation of Table 4.1:

G c = 8 3 σ 2 c E η
During the simulation, we increment the loading quasi-statically and we stop the simulation at the first crack nucleation event. We report the result obtained with the isotropic surface energy (ψ = 0) and for the non-vanishing anisotropic parameter (ψ = {5, 10}).

Results with the anisotropic surface energy

Figures 4.17-4.18 report the results obtained for the anisotropic surface energy with ψ = 5 and ψ = 10, respectively. We consider two internal lengths, η = H c /2 and η = H c /16. For η = H c /2 there is a phase with homogeneous damage that largely delays the nucleation of the cracks. The apparent delamination is still present for ψ = 5 and not relevant for ψ = 10. The bifurcation from solutions with uniform spacing to solution with non-uniform spacing is robustly detected also in the case with anisotropic surface energy. However, as anticipated in a previous section, the effect of the anisotropic parameter ψ is important and not well-controlled. For this reason we classify the results as unreliable and we will base the discussion of the rest of this chapter on the results of the isotropic model.

Crack spacing: comparisons with the sharp interface model

We analyze quantitatively the crack spacings obtained in the phase-field simulations and compare them to the results obtained with the sharp interface model in Chapter 3. For each phase-field simulation of Figure 4.14, we measure the crack spacing observed at the first nucleation event. All the simulations are characterised by either two segments of equal length or a short segment of length ℓ I and a long segment of length ℓ I I . The latter case can be obtained with one or two cracks. When two cracks are present, they are symmetrically placed. As sketched in Figure 4.15, the case with two symmetric cracks dividing the bar in three segments, can be regarded as the symmetrical configuration of the case with two cracks and two different segment lengths. For this reason in the case with three segments, we note as the length ℓ I I half of the length of the long central segment, see Figure 4.15.

In Figure 4.19, we report the segment lengths measured for each total bar length, and this for the several possible internal lengths η tests. For each set of simulations for a given internal length η, we observe a qualitative agreement with the results of Chapter 3, Figure 4.20: for short bars the crack is placed in the middle of the bar dividing it in two segments of equal length L/2. Starting from a critical length L c , the crack originates two different segments of length ℓ I and ℓ I I .

The point of bifurcation, in this diagram, indicates the splitting of the cell in two different pieces of unequal length. We can observe how the internal length, η, influences the bifurcation point L c . We test four η values and for each obtain a different bifurcation length L c . For η = H c /2 the bifurcation occurs at L c ≈ 3.5. For the case of η = H c /4 the bifurcation is delayed to L c ≃ 2.5, and to L c ≃ 2.0 and L c ≃ 1.5 for η = H c /8 and η = H c /16, respectively. We do not observe a convergent behaviour for η → 0. This is also because for small η the delamination becomes important and the cracking pattern changes in nature.

Our results indicate an important size effect on the observed crack spacings, depending on the value of the internal length η. Obviously, this effect is absent in the results of the sharp interface model of Chapter 3, Figure 4.20, where the length-scale η is not included in the model.

The results of Chapter 3 are in quantitative agreement to those obtained with the phase-field model only for η = H c /4, as shown in Figure 4.20. Note that η = H c /4 is the case in which the regularised crack pattern produced by the phase-field model is in better agreement to the purely transverse cracks hypothesis at the basis of the analysis of Chapter 3. For η = H c /2, the internal length is not small enough to resolve cracks in the central layer of thickness H c . For smaller values of η, the effects of the cracks predicted by the phase-field approach are not anymore purely transversal, with the presence of delamination-like cracks, see Figure 4.14

Full-scale simulations for the case of the T700GC/M21 laminate

In this section, we report the results of the full scale numerical simulations for the experimental laminate presented in Chapter 2 and specifically for the symmetric laminate configuration [0/90 2 ] s .

The laminate is made of transversely isotropic laminae, with the material properties detailed in Table 2.1. For the sake of simplicity, we consider the layer as isotropic, with the material properties associated depending on their orientation with respect to the axial direction of the bar: 0 • for the outer layer and 90 • for the central layer. Assuming the symmetry with respect to the central beam plane, we model only half of the overall thickness.

For the [0/90 2 ] s configuration, the stiffness ratio is E o /E c ≃ 15, the thickness of the outer layer is H o = 0.27 mm and the thickness of the inner layer, 2H c , is 4 × 0.27 mm, corresponding to a thickness ratio of H o /H c = 0.5. We assume that the outer layer is unbreakable, setting α = 0 on it. For the brittle central layer, the nominal toughness and strength for cracks are taken from the properties in the 90-degree direction of the composite: The analysis of Chapter 3 shows an important influence of the stiffness ratio on the energetic properties of the system and potentially on the crack evolution. To discuss, the influence of stiffness ratios in the evolutions obtained with the phase-field approach, we report the numerical results for the case E o /E c = 15, see Figure 4.21, corresponding to the experimental value, and the case E o /E c = 1, see Figure 4.22.

The following observations apply:

• For E o /E c = 15, several cracks nucleate at once at the first nucleation events (dimensionless loading t = 1). Their distribution is not regular at the beginning. The successive nucleation events 'fill the gaps' to tend to a final crack spacing with regular spacing.

• For E o /E c = 1, there are a long series of nucleation events corresponding to two additional cracks at the extremities. There is a progressive and regular propagation of the cracked zone from the boundaries to the interior. After that, a crack reaches the central zone, additional cracks 'fill the gaps' and reduce the crack spacings.

• In both situations, the crack nucleation events stop when delamination starts, introducing a minimal crack spacing.

Our results confirm that the nucleation and propagation of cracks are highly influenced from the stiffness ratio. The results for E o /E c = 15 are qualitatively in agreement with the experimental findings obtained through infrared thermography. However, the position of the cracks predicted through the phase field is not the one obtained through experimental procedures. In the high-stiffness ratio regime the exact crack position is not deterministic, being mainly determined by the imperfections. The presence of heterogeneities and/or imperfections was not taken into account in our simulations. It is very difficult to conclude if the nucleation position could ever be predicted for an experimental setting, especially without the use of any statistical tool. However, we find acceptable results based on the crack density rate, which is given by the evolution of the surface energy in the phase-field simulation.

Conclusion of this chapter

In this chapter, we have studied transverse cracking of a composite cross-ply laminate through a local minimization approach, commonly known as the phase field approach. We've introduced the approach, and have chosen the application of the classical AT 1 model. The simulations have been carried on a bi-dimensional bar, introduced and studied in Chapter 3. The main simulations introduced in the previous chapter are tested with the phase field approach and ultimately compared. The results found through global minimisation unveil the bifurcation phenomenon from solutions with uniform crack spacings to non uniformly spaced solutions. A qualitatively similar effect is found through local minimisation. Phase field approaches introduce another scale, the internal length η.

Through testing different η values, we observed how this scale influences the originating point of nucleation, therefore the cracking patterns resulting in either equally spaced or unequally spaced segments.

We compare the results obtained with the phase field with the ones obtained through global minimization and see further proof of a dependence on the internal length on solutions with uniform crack spacing to non uniformly spaced solutions. These results are visible through plots representing the segments' length for each evolving bar length L, showing how global minimisation and local minimisation converge to the same results only for a certain internal length value. Finally, we show the obtained results on a [0/90 2 ] s laminate where the material parameters shown in Chapter 2 are taken into account. The damage field obtained is qualitative and reasonable compared against the experimental data we possess, cracks are unequally distributed and propagate instantaneously through the thickness. These simulations show reasonable results, despite the limits of the model.

We observe that phase-field model produces delamination-like cracks which are not expected in the experimental results at this level of loading. With the aim of overcoming this numerical delamination, we introduced a modified approach of the AT 1 model. A penalization term is added through an anisotropic toughness constant and derivative of the damage field along the ordinates in the dissipated energy term. Multiple cases have been tested and discussed leaving the delamination issue unsolved. We were able to overcome delamination only for a few steps of the time-evolution problem. Nonetheless, the results of the anisotropic model are unreliable, leading to a somewhat correct crack nucleation threshold and crack density rate, related to the presence of a phase with homogeneous damage that we are not able to explain. We conclude that the introduction of an anisotropic toughness does not solve the anticipated delamination issue, rather introduces more complexity to the problem. Hence, we based our analysis on the basic isotropic AT 1 model.

Conclusions and perspectives

Conclusions

We studied the phenomenon of transverse fracture in laminated, fibre-reinforced, cross-ply composites in the framework of variational fracture mechanics. The reference system is a bi-dimensional brittle elastic domain consisting in a symmetrical stacking sequence of two layers constituted of different geometrical and material parameters.

Chapter 2 introduced the innovative experimental protocol of [START_REF] Berthe | Passive infrared thermography measurement of transverse cracking evolution in cross-ply laminates[END_REF], where infrared thermography is applied to capture the damaging process of the specimen. We presented the main experimental results of carbon-fibre T700GC/M21 material and a post-processing analysis of the cracking pattern. The experiments show the progressive appearance of a series of purely transverse cracks in the core of the sandwich structure. The infrared thermography data allowed us for a complete characterization of the crack density and spacing as a function of the loading. The post-processing of the experimental data revealed that the crack spacing is aperiodic and dependent on the thickness of the laminate and the loading rate.

In the rest of the thesis, we analyzed the problem of the optimal crack spacing from the theoretical and numerical perspectives. The available fracture models of laminated composites assume that the crack spacing in nominally perfect systems is uniform and that the aperiodicity of the cracking pattern is due to imperfections only. They are based on stress-based or energetic-based crack nucleation, or a combination of both. However, they always assume that when a new crack appears, in the absence of imperfections, the longest segment of the core layer breaks into two pieces of equal length.

Chapter 3 studied the crack spacing problem with an energetic approach based on a global energy minimisation criterion and a sharp interface Griffith model, applying the ideas of [START_REF] Bourdin | The Variational Approach to Fracture[END_REF]. Considering a bar with n transverse cracks in core, their optimal spacings for a given loading has been determined by the minimisation of the total potential energy. Within a modular approach, the properties of the minimizers has been related to the property of the elastic compliance S of the elementary cell of the system, showing that solutions with uniform crack spacings are stable if the second derivative with respect to the cell length, S ′′ (ℓ), is negative. We evaluated the derivatives of the elastic compliance by two-dimensional finite element computations and an accurate domain derivatives approach. We found that S ′′ (ℓ) is negative for short segments and positive for long segments. Our results unveil that even in perfect systems, solutions with uniform spacings are unstable for sufficiently long bars and bifurcate toward solutions with non-uniform spacings. We showed that, when considering the static problem for n cracks appearing simultaneously, the stable solution must be composed of a single long segment and n -2 short segments of equal length, modulo boundary effects. We found that the change of concavity of the elemental compliance S is barely visible or absent for classical cross-ply fiber-reinforced UD composite laminates, leading to particularly small stability margin of the solutions. This translates to cracking patterns with the same number of cracks having almost equivalent energies. In this regime we expect a notable dependency on imperfections. Vice versa, our parametric analysis indicates that the stability margin is larger for bars with thin and soft outer layers, whose crack spacings are expected to show a minor less imperfection-sensitiveness. We have shown that the simple shear lag models cannot predict the instability of solutions with uniform spacings.

In Chapter 4 we pursue our analysis with a variational phase-field fracture model of the composite bar. The quasi-static evolution is found by numerical (local) minimisation of the energy functional of a gradient damage model, with internal length, under an irreversibility constraint on the damage variable. We solved the quasi-static evolution numerically with an in-house finite element code. We showed that standard phase-field models predict delamination-like cracks at the interface at loading levels that are too small with respect to the experimental observations. We attempted to fix this issue by introducing an anisotropic fracture energy. The results of this approach being of difficult interpretation, we decided to pursue the analysis with the classical isotropic model. We focused on the analysis of cracking patterns at the first nucleation event, when the effect of the delamination is negligible. Our numerical results with the phase-field approach are in qualitative agreement with the results of the sharp interface model of Chapter 3. For sufficiently long bars, they show a bifurcation from solutions with uniform spacings to solutions with non-uniform spacings composed by a single long segment. We performed full-size numerical simulations for the experimental system of Chapter 2, identifying the internal length of the damage model from the available material data for the composite critical toughness and critical stress values. The results confirmed that the crack spacings appear to be deterministic for bars with relative thin and soft outer layers, while it is of random nature in the case of stiff outer layers. We found that the internal length of the phase-field model, η, plays an unexpected crucial role on crack spacings, therefore the cracking patterns and has an equally important influence on the critical bar length for the bifurcation from uniform to non-uniform crack spacings. We are not able to furnish a clear explanation for this phenomenon.

Perspectives

The present work opens for several perspectives from the theoretical, numerical, and experimental viewpoints.

Experimentally, our results on one hand confirm that the crack spacings on the classical [0/90 n ] s composite laminate are most probably of imperfection-sensitive random nature. On the other hand, they suggest the presence of interesting regimes of deterministic and non-uniform crack spacings for the case of laminates with softer and thin outer layers. They naturally call for a wider experimental campaign to test this unexplored regime predicted by theory. For the theoretical and numerical point of view, we have identified several shortcomings of the classical phase-field models in solving the fracture mechanics problem in composite structures. They are essentially related to known issues in correctly predicting the crack nucleation threshold under multi-axial loading conditions. The solutions available in the literature, including energy splits of different nature or non variational approaches, present several drawbacks. Further works on this topic are necessary to be able to correctly predict composite fracture with phase-field approaches with full-field simulations. A possible attractive alternative is to apply specific structural models [START_REF] Alessi | Phase-field modelling of failure in hybrid laminates[END_REF] introducing ad hoc hypothesis on the cracking modes as the one proposed in [START_REF] Alessi | Phase-field modelling of failure in hybrid laminates[END_REF], [START_REF] Bouteiller | Simulation de la rupture des composites stratifiés dans le cadre de modèles de plaques multicouches[END_REF]. In particular, it would be specifically appealing to pursue the analysis of crack spacings and its dependence on the material and geometric properties by using the advanced laminate plate model of [START_REF] Bouteiller | Simulation de la rupture des composites stratifiés dans le cadre de modèles de plaques multicouches[END_REF], which is able to keep fine properties of the three-dimensional elastic solution. We expect that, differently from standard shear-lag-like approaches [START_REF] Alessi | Phase-field modelling of failure in hybrid laminates[END_REF], such a model should be able to correctly predict the presence of non-uniform crack spacings. This would be an ideal framework also for the numerical investigation of the competition between transverse fracture and delamination. Including the effect of the delamination in the sharp interface approach of Chapter 3 is a further straightforward interesting extension of the present work [see e.g. León [START_REF] León Baldelli | Fracture and debonding of a thin film on a stiff substrate: Analytical and numerical solutions of a one-dimensional variational model[END_REF][START_REF] León Baldelli | A variational model for fracture and debonding of thin films under in-plane loadings[END_REF].

A subtler unsolved problem is the understanding of the influence of the internal length η of the phase-field model on the crack spacings. We do not have any convincing explanation of our observations. They can be related to the subtle differences between the sharp interface approach of Chapter 3 based on a global energy minimisation approach and the phase-field model based on local energy minimisation introducing the scale effect related to the internal length.

Our comparisons between the experimental and numerical results is preliminary and incomplete, being almost of qualitative nature. Further efforts on both sides are surely required to get quantitatively accurate numerical predictions of the experimental findings. This can include the introduction of controlled imperfections in our variational theoretical and numerical models to assess their influence on the results. Figures A.3 and A.4 show the analytical derivatives of the potential energy with respect to the external radius. As it can be seen in these graphs, the analytical and numerical results overlap perfectly, which allows us to say that the numerical implementation of the Gθ method at both order 1 and 2 is a reliable tool to calculate these shape derivatives.

1. 1

 1 Example of fibers, matrix and laminae . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Taxonomies reproduced from Ibrahim et al. [2015]. . . . . . . . . . . . . . . . . . . . . . 1.3 Example of symmetrical cross-ply and angle-ply laminates . . . . . . . . . . . . . . . . 1.4 Schematic view of matrix cracking, delamination and fibre breakage in a [0 • /90 • ] UD laminate . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.5 "Transverse cracking in specimens with transverse ply thickness of (a) 0.42mm, (b) 1.2mm, (c) 2mm, and (d) 4 mm. All specimens have been loaded to 1.4% strain." Figure and description from Parvizi and Bailey [1978] . . . . . . . . . . . . . . . . . . . . . . . 1.6 "(a) Typical transverse crack for a specimen with transverse ply thickness of 1.2mm; there is no evidence of debonding between the plies. (b) Transverse cracking at a transverse ply thickness of 4mm; small amount of debonding between the plies and 45 • oblique cracks are observable." Figure and description from Parvizi and Bailey [1978] 1.7 "The strain to initiate micro-cracking in E glass/Shell Epikote epoxy [0/90/0] laminates as a function of the total thickness of the 90 • plies. The 0 • plies each have a constant thickness of 9.5mm". Figure and description from Nairn [2000] whom reproduced the data from Parvizi et al. [1978]. . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.8 "Four cases of modelled crack spacings: minimum (a), maximum (b), mean (c), and random (d)." Figure and description from Silberschmidt [2005]. The studied specimen has dimensions of 0.75mm × 10mm and total height of 6 × 0.125mm ply thickness. . . . 1.9 Example of modelling of cohesive zone ahead of crack tip. . . . . . . . . . . . . . . . . 1.10 Example of a general cohesive law showing the relationship between the cohesive traction σ and the displacement δ. Image from Sun and Jin [2012] . . . . . . . . . . . . . . 2.1 Spatial distribution of the thermal events for following loading rates (a): 5mm/min, (b): 50mm/min, (c): 500mm/min and (d): 0.1m/s, in a [0/90] s laminate. Images from Coussa et al. [2022]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 "Thickness effect on the crack density captured with the IR protocol and based on tensile tests performed at 5 mm/min.", Image and caption from Coussa et al. [2022] . . 2.3 Evolution of the crack density measured with respect to the applied stress, for three different stacking sequences: (a) [0/90] s , (b) [0/90 2 ] s , (c) [0/90 3 ] s . The stacking sequences were tested at 4 different upper holder speeds: 0.1 m/s, 500 mm/min, 50 mm/min and 5 mm/min. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 Transverse cracking patterns on the [0/90 3 ] s laminate prior the final failure of the specimen subject to various loading rates: 500mm/min, 0.1m/s, 1.4m/s. Figure from Coussa et al. [2022] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5 "Values of the variation of the strain ϵ tu , as a function of ply thickness 2d and experimental values of ϵ tu for various ply thicknesses in glass fibre/epoxy 90 • cross-ply laminates. The horizontal line depicts the limiting values of ϵ tu for large inner-ply thicknesses." Figure and caption from Parvizi et al. [1978]. . . . . . . . . . . . . . . . . 2.6 Schematic view of a specimen divided in two segments of length ℓ 1 and ℓ 2 defined by the pixel position of the corresponding appearing crack. . . . . . . . . . . . . . . . . . 2.7 Reconstruction of crack distribution in a [0/90 3 ] s laminate subject to three main different upper hold speeds. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.8 Length of the longest segment as the number of cracks increases in a [0/90 3 ] s specimen subjected to (a): 5mm/min, (b): 50mm/min and (c): 500mm/min upper hold speed test. 2.9 Histograms showing the percentage for which the largest segment breaks first, in a[0/90 3 ] s laminate. The percentage is calculated through the ratio of the normalized longest new broken segment over the length of its original segment. . . . . . . . . . . 3.1 Top: Sandwich bar with n -1 transverse cracks in the core layer. Bottom: boundary and bulk modulus and elementary cell of the modular system . . . . . . . . . . . . . . 3.2 Potential energy of the composite bar as a function of the crack position for different bar lengths. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Composite bar alternative with displacement imposed only on the outer layer. In this case, differently from the case in Figure 3.1 there is no need for a distinction between boundary modules and bulk modules. All the modules are of the bulk type of Figure 3.1, see also Remark 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Typical finite element mesh in the deformed configuration used for the computation of the basic cell compliance S(ℓ) for H c = 1, H o = 0.5, L = 4 (outer layer in cyan and core layer in yellow). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 First and second derivative of the compliance S(ℓ) for a composite bar made of isotropic materials with E o /E c = 1, H o /H c = 0.5, ν = 0.3. The numerical results are for obtained for a domain with H c = 1. Here and henceforth, all the lengths should be regarded as non-dimensional with respect to the scaling length H c , which is not reported explicitly in the plot labels. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6 Length of short segment ℓ I as a function of the length of the long segment ℓ I I . Left: graphical illustration of the problem of finding ℓ I < ℓ c such that S ′ (ℓ I ) = S ′ (ℓ I I ) given ℓ I I > ℓ c . Right: resulting function ℓ I (ℓ I I ), showing also the trivial solution ℓ I = ℓ I I , in blue. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.7 Solution for a bar divided in a short and long segment by a single crack (n = 2). Left: lengths of the short (ℓ I ) and long segment (ℓ I I ) as a function of the total bar length L with stable solution marked in blue and unstable solutions in pink. Right: Total elastic compliance of the homogeneous solution and the solution with two different segment lengths. Results for E o /E c = 1, H o /H c = 0.5. . . . . . . . . . . . . . . . . . . . . . . . . 3.8 Bar divided in n > 2 segments by n -1 cracks. Left: Different kinds of stable solutions for the case n = 5: homogeneous (H); long segment at the boundary (A); long segment not at the boundary (B). Right: Length of the short (ℓ I ) and long segment ℓ I I for the bar divided in two segments with stable solution in blue and unstable solutions in pink: possible solutions for different values of n ≥ 2 are reported. . . . . . . . . . . . . . . . 3.9 Dependence of derivatives of the compliance of the basic module on the material and geometric parameters. The numerical results for the first (S ′ (L), top row) and second (S ′′ (L), bottom row) derivatives for different thickness (left column) and stiffness ratios (right column) are obtained with the method of Section 3.3.3. Results for isotropic materials with ν = 0.3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.10 Basic module compliance derivative S ′ (L) for 0/90/0 composite laminate made of typical Glass Fiber laminae (E x = 43.5 GPa, E y = 11.5 GPa, G xy = 3.45 GPa, ν xy = 0.25) or typical carbon fiber laminae (E x = 132 GPa, E y = 10.8 GPa, G xy = 5.65 GPa, ν xy = 0.24). The results are for different thickness ratios. . . . . . . . . . . . . . . . . . 3.11 Influence of the shear stiffness G xy on the derivative of the compliance of the basic module. The results are a 0/90/0 composite laminate made of an ideal material with E x = 43.5 GPa, E y = 11.5 Gpa, G xy = 3.45 GPa, ν xy = 0.25 and varying G xy . . . . . . . 3.12 Left: Total energy F n (3.37) for different number of cracks, where the displacement is scaled with u 0 = G c E c H c . The grid lines and the black dots denote the critical loadings for which F n = F n + 1. The dashed red curve is F n * with n * such that dF n /dn = 0, which gives an approximation of the lower envelope of the F n 's. Right: crack density as a function of the average stress in the bar corresponding to the red lower envelope in the left plot. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.13 Optimal number of segment n cracks, as a function of the loading and the bar length, where n segments correspond to n -1 cracks. . . . . . . . . . . . . . . . . . . . . . . . . 3.14 1-Dimensional model of the composite bar presented in Figure 3.1. . . . . . . . . . . . 3.15 Second derivatives of the elastic compliance S"(L) obtained with the simplified 1D shear-lag approach. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1 Brittle elastic two-dimensional domain with a given crack path Γ. . . . . . . . . . . . . 4.2 One-dimensional example of a stress-softening bar subject to imposed displacement. . 4.3 Stress-strain solutions for AT1 and AT2 models (a-b). Stability diagrams for AT1 and AT2 models (c-d). Figure from Pham et al. [2011]. . . . . . . . . . . . . . . . . . . . . . 4.4 Damage with localised solution for AT1 model. Figure from Pham et al. [2011]. . . . . 4.5 Homogeneous response of the AT m model (4.25). . . . . . . . . . . . . . . . . . . . . . . 4.6 Results of a simulation on a bar under imposed displacement. Evolution of the elastic and dissipated energy of the traction test using AT1 model. Material, and geometrical parameters used: E 0 = 1., G c = 1., η = .1, bar's length L = 1., width W = .1 . . . . . . 4.7 Left: geometry and boundary conditions for the sandwich beam finite element simulations, taking into account the symmetry with respect to the axis of the beam to reduce it to a two-layer geometry. Right: typical finite element mesh for the two-layer geometry obtained after accounting for the symmetry. . . . . . . . . . . . . . . . . . . . . . . . . . 4.8 (Left): G c Fixed at 1. (Right): σ c fixed at 1. Simulations with the following material and geometric parameters: H o /H c = 0.25, internal length η = 0.5, Young moduli E c = 1.0 and E o = 1.0. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.9 H o /H c = 0.25, internal length η = 0.1, Young moduli E c = 1.0 and E o = 1.0 . . . . . . 4.10 H o = 0.5, H c = 1, E o = E c = 1.0, η = 0.25 . . . . . . . . . . . . . . . . . . . . . . . . . . 4.11 Schematic view of a crack in a bi-dimensional domain propagating along θ with n indicating the normal to the crack path. . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.12 Snapshots of the damage at the first appearance (top image) of cracks and at t = 5 (bottom image) for different values of the anisotropic parameter ψ in the surface energy (4.29). Results for L = 10, H o = 0.5, H c = 1, η = 0.25, G c = 1. . . . . . . . . . . . 4.13 Elastic and dissipated energies during the quasi-static evolution for the different values of the anisotropic coefficient ψ in the surface energy (4.29). Results for L = 10, H o = 0.5, H c = 1, η = 0.25, G c = 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.14 Damage field at the first nucleation as a function of the internal length η and the bar length L. Results obtained with the AT 1 model with isotropic surface energy (ψ = 0) and for E c /E o = 1, H o /H c = 0.5, ν c = ν o = 0.3, σ c = 1. . . . . . . . . . . . . . . . . . . 4.15 The bar with two symmetric cracks can be regarded as the symmetrical version of the case of one crack breaking the bar in a short and a long segment. . . . . . . . . . . . . . 4.16 Influence of the internal length (top-left: η = 0.5, top-right η = 0.25, bottom-left η = 0.125, bottom-right: η = 0.0625) on the evolution for a short bar (L = 1). Results for the AT 1 model with, L = 1, E o /E c = 1, H o /H c = 0.5, ν = 0.3, ψ = 0, σ c /E o = 1, H c = 1. 4.17 Phase field model with anisotropic surface energy (ψ = 5). Damage fields at the first nucleation as a function of the internal length η and the bar length L. Results obtained with the AT 1 model for E c /E o = 1, H o /H c = 0.5, ν c = ν o = 0.3, σ c = 1. . . . . . . . . . 4.18 Phase field model with anisotropic surface energy (ψ = 10). Damage fields at the first nucleation as a function of the internal length η and the bar length L. Results obtained with the AT 1 model for E c /E o = 1, H o /H c = 0.5, ν c = ν o = 0.3, σ c = 1. . . . . . . . . . 4.19 Segments' lengths obtained with phase field AT1 results for the isotropic model of Figure 4.14. Bar lengths varying from L = 1 to L = 15, the thickness ratio is H o /H c = 0.5 and the stiffness ratio E c /E o = 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.20 Comparison between results for sharp-interface model (solid red and blue lines) against phase field AT1 model for η : H c /4 (dotted, black line). Bifurcations found for L c = 2.5 in both cases, for tested bars evolving up to L = 15. Results for stiffness ratio E o /E c = 1 and thickness ratio H o /H c = 0.5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.21 Full scale numerical simulation with the phase-field approach for a [0/90 2 ] s composite bar made of T700GC/M21 laminae. Top: Evolution of the elastic and surface energies, where the discontinuity of the surface energy corresponds to nucleation events. Bottom: snapshots of the phase-field immediately after the different nucleation events. Material and geometric parameters (isotropic AT 1 model): E o /E c = 15, H o /H c = 0.5, ν c = ν o = 0.3, η/H c = 0.32, σ c /E c = 1, L/H c = 100. For visualization purposes, the elastic energy is cleared of a trivial term corresponding to E o H o t 2 /2, corresponding to the stiffness of the outer unbreakable layer. . . . . . . . . . . . . . . . . . . . . . . . . . 4.22 Full scale numerical simulation with the phase-field approach for a [0/90 2 ] s composite bar with stiffness ratio E o /E c = 1. All the other parameters as in Figure 4.21 . . . . . . A.1 Schematic drawing of one fourth of the total geometry, reproducing a cylinder subject to both internal and external pressure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.2 Plotting of the potential energy: comparison between analytical (solid line) solution and numerical (red circles) solution. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . A.3 Plotting of the first energy derivative with respect to the external radius: comparison between analytical (solid line) solution and numerical (red circles) solution. . . . . . . A.4 Plotting of the second energy derivative with respect to the external radius: comparison between analytical (solid line) solution and numerical (red circles) solution. . . . . 126

  are composite materials ? . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1.2 Unidirectional plies and laminates . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1.3 Main stacking sequences . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Main damage and failure mechanisms in cross-ply laminates . . . . . . . . . . . . . . 1.2.1 Intralaminar failure in laminated composites: matrix-cracking . . . . . . . . . . 1.2.2 Interlaminar failure in laminated composites: delamination . . . . . . . . . . . 1.2.3 Translaminar failure in laminated composites: fibre-breakage . . . . . . . . . . 1.3 State of the art on prediction of matrix-cracking patterns in UD cross-ply laminates 1.3.1 Overview of matrix-cracking in composite fibre reinforced cross-ply laminates 1.4 A short introduction to numerical approaches for the simulation of damage mechanisms at different length scales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4.1 Cohesive zone models (CZM) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4.2 Continuum damage models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4.3 Variational approaches to fracture . . . . . . . . . . . . . . . . . . . . . . . . . . 1.5 Phase field models and their application to composite materials . . . . . . . . . . . .

  FIGURE 1.1: Example of fibers, matrix and laminae

1. 2 .

 2 Main damage and failure mechanisms in cross-ply laminates 27 (A) A taxonomy of reinforcement materials in composites. (B) A taxonomy of matrices used in composites.

FIGURE 1 . 2 :

 12 FIGURE 1.2: Taxonomies reproduced from Ibrahim et al. [2015].

Chapter 1 .

 1 FIGURE 1.3: Example of symmetrical cross-ply and angle-ply laminates

FIGURE 1

 1 FIGURE 1.5: "Transverse cracking in specimens with transverse ply thickness of (a) 0.42mm, (b) 1.2mm, (c) 2mm, and (d) 4 mm. All specimens have been loaded to 1.4% strain." Figure and description from Parvizi and Bailey[1978]

FIGURE 1

 1 FIGURE 1.8: "Four cases of modelled crack spacings: minimum (a), maximum (b), mean (c), and random (d)." Figure and description from Silberschmidt [2005]. The studied specimen has dimensions of 0.75mm × 10mm and total height of 6 × 0.125mm ply thickness.

FIGURE 1 . 9 :

 19 FIGURE 1.9: Example of modelling of cohesive zone ahead of crack tip.

FIGURE 1 .

 1 FIGURE 1.10: Example of a general cohesive law showing the relationship between the cohesive traction σ and the displacement δ. Image from Sun and Jin [2012]

  FIGURE 2.1: Spatial distribution of the thermal events for following loading rates (a): 5mm/min, (b): 50mm/min, (c): 500mm/min and (d): 0.1m/s, in a [0/90] s laminate. Images from Coussa et al. [2022].

  FIGURE 2.2: "Thickness effect on the crack density captured with the IR protocol and based on tensile tests performed at 5 mm/min.", Image and caption from Coussa et al.[2022]

2. 4 .

 4 FIGURE 2.3: Evolution of the crack density measured with respect to the applied stress, for three different stacking sequences: (a) [0/90] s , (b) [0/90 2 ] s , (c) [0/90 3 ] s . The stacking sequences were tested at 4 different upper holder speeds: 0.1 m/s, 500 mm/min, 50 mm/min and 5 mm/min.

FIGURE 2 . 4 :

 24 FIGURE 2.4: Transverse cracking patterns on the [0/90 3 ] s laminate prior the final failure of the specimen subject to various loading rates: 500mm/min, 0.1m/s, 1.4m/s. Figure from Coussa et al. [2022]

  FIGURE 2.6: Schematic view of a specimen divided in two segments of length ℓ 1 and ℓ 2 defined by the pixel position of the corresponding appearing crack.

  ℓ 2 := (306 -284)/281 = 0.0782. (A) 5 mm/min (B) 50 mm/min (C) 500 mm/min

FIGURE 2

 2 FIGURE 2.7: Reconstruction of crack distribution in a [0/90 3 ] s laminate subject to three main different upper hold speeds.

  FIGURE 2.9: Histograms showing the percentage for which the largest segment breaks first, in a [0/90 3 ] s laminate. The percentage is calculated through the ratio of the normalized longest new broken segment over the length of its original segment.

Chapter 3 .

 3 Analytical and numerical approaches implemented to study transverse cracking patterns in UD cross-ply laminates formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58 3.2.2 The static variational problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 3.2.3 An example of solutions with non-uniform spacing: a bar with a single crack . 60 3.3 Modular energetic approach for the case with n cracks . . . . . . . . . . . . . . . . . . 60 3.3.1 The elastic energy as a function of the crack spacing . . . . . . . . . . . . . . . . 61 3.3.2 The optimal crack spacing problem . . . . . . . . . . . . . . . . . . . . . . . . . 62 3.3.3 Elementary cell problem: computing the elastic compliance of the elementary cell and its derivatives . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 3.4 Solutions for the crack spacing problem and their stability . . . . . . . . . . . . . . . 66 3.4.1 The case of one crack dividing the bar in two segments . . . . . . . . . . . . . . 66 3.4.2 General results for n-1 cracks dividing the bar in n segments . . . . . . . . . . simplified shear-lag approach . . . . . . . . . . . . . . . . . . . . . . . . . . . 74 3.7.2 Elastic compliance with the shear-lag model . . . . . . . . . . . . . . . . . .

FIGURE 3

 3 FIGURE 3.1: Top: Sandwich bar with n -1 transverse cracks in the core layer. Bottom: boundary and bulk modulus and elementary cell of the modular system

  .2 at the bottom. The Figure in the left reports P L (ℓ) as a function of ℓ for different bar lengths L for the case of two isotropic layers with Young moduli ratio E o /E c = 5, Poisson's coefficients ν c = ν o = 0.3, and thickness ratio H o /H c = 1. We observe two different regimes: • For short bars (L ≤ L (c) ≃ 7H c ) the elastic energy has a single minimum at ℓ * = L/2, corresponding to a crack dividing the bar in two segments of equal length.

  FIGURE 3.2: Potential energy of the composite bar as a function of the crack position for different bar lengths.

  FIGURE 3.4: Typical finite element mesh in the deformed configuration used for the computation of the basic cell compliance S(ℓ) for H c = 1, H o = 0.5, L = 4 (outer layer in cyan and core layer in yellow).

FIGURE 3 . 5 :FIGURE 3

 353 FIGURE 3.5: First and second derivative of the compliance S(ℓ) for a composite bar made of isotropic materials with E o /E c = 1, H o /H c = 0.5, ν = 0.3. The numerical results are for obtained for a domain with H c = 1. Here and henceforth, all the lengths should be regarded as non-dimensional with respect to the scaling length H c , which is not reported explicitly in the plot labels.

  FIGURE 3.7: Solution for a bar divided in a short and long segment by a single crack (n = 2). Left: lengths of the short (ℓ I ) and long segment (ℓ I I ) as a function of the total bar length L with stable solution marked in blue and unstable solutions in pink. Right: Total elastic compliance of the homogeneous solution and the solution with two different segment lengths. Results for E o /E c = 1, H o /H c = 0.5.

Remark 2

 2 FIGURE 3.8: Bar divided in n > 2 segments by n -1 cracks. Left: Different kinds of stable solutions for the case n = 5: homogeneous (H); long segment at the boundary (A); long segment not at the boundary (B). Right: Length of the short (ℓ I ) and long segment ℓ I I for the bar divided in two segments with stable solution in blue and unstable solutions in pink: possible solutions for different values of n ≥ 2 are reported.

FIGURE 3 . 9 :

 39 FIGURE 3.9: Dependence of derivatives of the compliance of the basic module on the material and geometric parameters. The numerical results for the first (S ′ (L), top row) and second (S ′′ (L), bottom row) derivatives for different thickness (left column) and stiffness ratios (right column) are obtained with the method of Section 3.3.3. Results for isotropic materials with ν = 0.3.

  FIGURE 3.10: Basic module compliance derivative S ′ (L) for 0/90/0 composite laminate made of typical Glass Fiber laminae (E x = 43.5 GPa, E y = 11.5 GPa, G xy = 3.45 GPa, ν xy = 0.25) or typical carbon fiber laminae (E x = 132 GPa, E y = 10.8 GPa, G xy = 5.65GPa, ν xy = 0.24). The results are for different thickness ratios.

FIGURE 3

 3 FIGURE 3.11: Influence of the shear stiffness G xy on the derivative of the compliance of the basic module. The results are a 0/90/0 composite laminate made of an ideal material with E x = 43.5 GPa, E y = 11.5 Gpa, G xy = 3.45 GPa, ν xy = 0.25 and varying G xy .

Chapter 3 .

 3 FIGURE 3.12: Left: Total energy F n (3.37) for different number of cracks, where the displacement is scaled with u 0 = G c E c H c . The grid lines and the black dots denote the critical loadings for which F n = F n + 1. The dashed red curve is F n * with n * such that dF n /dn = 0, which gives an approximation of the lower envelope of the F n 's. Right: crack density as a function of the average stress in the bar corresponding to the red lower envelope in the left plot.

  FIGURE 3.13: Optimal number of segment n cracks, as a function of the loading and the bar length, where n segments correspond to n -1 cracks.

FIGURE 3 .

 3 FIGURE 3.15: Second derivatives of the elastic compliance S"(L) obtained with the simplified 1D shear-lag approach.

  [2012] and phase-field fracture models Quintanas-Corominas et al. [2019], Bleyer and Alessi [2018], León Baldelli et al. [2014]. Investigating the crack spacing issue in the context of phase-field fracture models is the objective of the next chapter.

4

 4

FIGURE 4 .

 4 FIGURE 4.1: Brittle elastic two-dimensional domain with a given crack path Γ.

FIGURE 4 .

 4 FIGURE 4.2: One-dimensional example of a stress-softening bar subject to imposed displacement.

FIGURE 4 .

 4 FIGURE 4.3: Stress-strain solutions for AT1 and AT2 models (a-b). Stability diagrams for AT1 and AT2 models (c-d). Figure from Pham et al. [2011].

FIGURE 4 . 4 :

 44 FIGURE 4.4: Damage with localised solution for AT1 model. Figure from Pham et al. [2011].

FIGURE 4 . 5 :

 45 FIGURE 4.5: Homogeneous response of the AT m model (4.25).

  model

FIGURE 4 .

 4 FIGURE 4.6: Results of a simulation on a bar under imposed displacement. Evolution of the elastic and dissipated energy of the traction test using AT1 model. Material, and geometrical parameters used: E 0 = 1., G c = 1., η = .1, bar's length L = 1., width W = .1

FIGURE 4 .

 4 FIGURE 4.7: Left: geometry and boundary conditions for the sandwich beam finite element simulations, taking into account the symmetry with respect to the axis of the beam to reduce it to a two-layer geometry. Right: typical finite element mesh for the two-layer geometry obtained after accounting for the symmetry.

  m /(2(1 + ν m )).We solve the problem using the time discrete evolution principle based on the local minimisation irreversibility constraint on the damage:

•

  the relative thickness H o /H c • the relative Young modulus E o /E c Chapter 4. Phase-field approach to fracture of composite laminates • the dimensionless length of the bar L/H c • the length-scale

FIGURE 4 .

 4 FIGURE 4.8: (Left): G c Fixed at 1. (Right): σ c fixed at 1. Simulations with the following material and geometric parameters: H o /H c = 0.25, internal length η = 0.5, Young moduli E c = 1.0 and E o = 1.0.

1

  E o = 1 E c = 1 ν = 0.3 Hc = 1 H o = 0.5, ψ = 1, η = 0

FIGURE 4 .

 4 FIGURE 4.9: H o /H c = 0.25, internal length η = 0.1, Young moduli E c = 1.0 and E o = 1.0

FIGURE 4 .

 4 FIGURE 4.11: Schematic view of a crack in a bi-dimensional domain propagating along θ with n indicating the normal to the crack path.

FIGURE 4 .FIGURE 4 .

 44 FIGURE 4.12: Snapshots of the damage at the first appearance (top image) of cracks and at t = 5 (bottom image) for different values of the anisotropic parameter ψ in the surface energy (4.29). Results for L = 10, H o = 0.5, H c = 1, η = 0.25, G c = 1.

FIGURE 4 .

 4 FIGURE 4.14: Damage field at the first nucleation as a function of the internal length η and the bar length L. Results obtained with the AT 1 model with isotropic surface energy (ψ = 0) and for E c /E o = 1, H o /H c = 0.5, ν c = ν o = 0.3, σ c = 1.

FIGURE 4 .

 4 FIGURE 4.16: Influence of the internal length (top-left: η = 0.5, top-right η = 0.25, bottom-left η = 0.125, bottom-right: η = 0.0625) on the evolution for a short bar (L = 1). Results for the AT 1 model with, L = 1, E o /E c = 1, H o /H c = 0.5, ν = 0.3, ψ = 0, σ c /E o = 1, H c = 1.

FIGURE 4 .

 4 FIGURE 4.17: Phase field model with anisotropic surface energy (ψ = 5). Damage fields at the first nucleation as a function of the internal length η and the bar length L. Results obtained with the AT 1 model for E c /E o = 1, H o /H c = 0.5, ν c = ν o = 0.3, σ c = 1.

FIGURE 4 .

 4 FIGURE 4.18: Phase field model with anisotropic surface energy (ψ = 10). Damage fields at the first nucleation as a function of the internal length η and the bar length L. Results obtained with the AT 1 model for E c /E o = 1, H o /H c = 0.5, ν c = ν o = 0.3, σ c = 1.

FIGURE 4 .

 4 FIGURE 4.19: Segments' lengths obtained with phase field AT1 results for the isotropic model of Figure 4.14. Bar lengths varying from L = 1 to L = 15, the thickness ratio is H o /H c = 0.5 and the stiffness ratio E c /E o = 1.

FIGURE 4 .

 4 FIGURE 4.20: Comparison between results for sharp-interface model (solid red and blue lines) against phase field AT1 model for η : H c /4 (dotted, black line). Bifurcations found for L c = 2.5 in both cases, for tested bars evolving up to L = 15. Results for stiffness ratio E o /E c = 1 and thickness ratio H o /H c = 0.5.

  E o = 8.5 GPa σ c = 80 MPa, G c = 350 J/m 2 . (4.30)For the AT 1 model we find:

FIGURE 4 .

 4 FIGURE 4.21: Full scale numerical simulation with the phase-field approach for a [0/90 2 ] s composite bar made of T700GC/M21 laminae. Top: Evolution of the elastic and surface energies, where the discontinuity of the surface energy corresponds to nucleation events. Bottom: snapshots of the phase-field immediately after the different nucleation events. Material and geometric parameters (isotropic AT 1 model):E o /E c = 15, H o /H c = 0.5, ν c = ν o = 0.3, η/H c = 0.32, σ c /E c = 1, L/H c = 100.For visualization purposes, the elastic energy is cleared of a trivial term corresponding to E o H o t 2 /2, corresponding to the stiffness of the outer unbreakable layer.

110Chapter 4 .

 4 Phase-field approach to fracture of composite laminates

FIGURE 4 .

 4 FIGURE 4.22: Full scale numerical simulation with the phase-field approach for a [0/90 2 ] s composite bar with stiffness ratio E o /E c = 1. All the other parameters as in Figure 4.21

  FIGURE A.2: Plotting of the potential energy: comparison between analytical (solid line) solution and numerical (red circles) solution.

FIGURE A. 4 :

 4 FIGURE A.4: Plotting of the second energy derivative with respect to the external radius: comparison between analytical (solid line) solution and numerical (red circles) solution.

TABLE 2

 2

	Properties		Quantity Unity Value
	Traction**	Module	E 11 /E 22	GPa	130/8.5
		Maximum Stress σ 11 /σ 22	MPa 2126/80
	Compression** Module	E 11	GPa	114
		Maximum Stress σ 11	MPa -1409
	Shear**	Module	G 12	GPa	4.4
		Maximum Stress σ 12	MPa 89
		Poisson ratio**	ν 12	/	0.32
	Tenacity***	Mode

.1: Properties of UD T700GC/M21 carbon/epoxy laminate

TABLE 2

 2

.2: Table from Hautier [2010], ** "Properties identified during the program AMERICO (values are means of results of four laboratories (EADS, LGMT, ONERA, ENSICA) , *** Values from Prombut [2007]

  However, the Hessian matrix (3.16) is negative definite only if S ′′ (ℓ * ) ≤ 0. The solution with uniform spacing is stable only if S ′′ (ℓ * ) ≤ 0.Solutions with non-uniform spacings for (3.15) require the existence of cell lengths ℓ * I > 0 and ℓ *

II > ℓ * I such that S ′ (ℓ * I ) = S ′ (ℓ * I I

(A) 5 mm/min (B) 50 mm/min (C) 500 mm/min

(A) (B) (C)

Acknowledgements

Chapter 2. A first analysis of the physical cracking patterns

Finite element discretisation

The energy functional (4.11) is discretised using a standard finite element approach. We use a mesh made of triangular unstructured cells with a classical piece-wise linear interpolation of the displacement and the damage field (P 1 elements), having the nodal values as degrees of freedom. We use an unstructured mesh with uniform mesh size h. The mesh size is set to be sufficiently smaller with respect to the internal length η to be at least half its length, i.e. h = η/2 and no smaller than h = η/10.

Computational tools

The numerical results presented here are obtained through a finite element code developed in the group, damage-dolfinx. The code based on the finite element library FEniCSx and it is written in python. More in detail we use the following open-source tools:

• gmsh for mesh generation (python interface)

• FEniCSx for the finite element discretisation and data management

• PETSc for the linear (direct/iterative) and the nonlinear solvers. The bound-constrained minimisation problem is solved using the dedicated solvers SNESVI or TAO available in PETSc.

• Paraview and pyvista for visualization.

The code works in parallel (MPI) and runs on supercomputers. The computations reported here are of limited size and run on laptop in times spanning from minutes to 1h.

The bar traction test

We show the numerical solution to the minimization problem in Figure 4.6 for the AT1 model of the one-dimensional bar subjected to imposed displacement. We observe an initial elastic phase followed by a localisation in the form of a crack. The critical load for the crack nucleation predicted by the stability analysis of the damage model is

.6 reports also the critical load obtained when taking as evolution law the global minimisation of the Griffith energy. In the latter case, the nucleation is obtained when the total energy of cracked solution, G c W, becomes smaller than the total energy of the elastic solution E 0 U 2 t L/2, which gives the critical load:

We can see from this figure that the elastic limit predicted from the Griffith approach and the elastic limit predicted through the AT1 model are different.

Chapter 4. Phase-field approach to fracture of composite laminates • For short bars, we observe a single crack breaking the bar in two segments of equal length L/2.

Results with the isotropic surface energy

• For sufficiently long bars, we observe a non-uniform crack spacing with either one or two cracks, breaking the bar in one long and one or two short segments.

• In agreement with the prediction of Chapter 3, there is only one long segment and a short segment and they are always at the boundaries.

• In the configuration with two cracks, they are symmetrically placed. This configuration can be regarded as the symmetrical version of the case with a single crack breaking the bar in a short and long segment, as sketched in Figure 4.15.

• The internal length η plays an important role, affecting the critical length for the bifurcation from the solution with uniform and non-uniform crack spacings. This is not in agreement with the results of Chapter 3, which ignores the effect of the internal length.

• For long bars, the delamination is present at the first nucleation event. We will ignore this case with spurious delamination in our analysis. The evolution can be with no cracks (η large), with a smooth growth of a crack, or with a brutal crack nucleation (η small).

Appendices

A

Calculus of shape derivatives

In Chapter 3 we discuss solutions to a maximization problem of the total elastic compliance (3.15) and the stability given from the Hessian matrix 3.16. The purpose of this appendix is to introduce an alternative method that has been developed in this work to calculate first and second order shape derivatives of the elastic/potential energy and compliance with respect to the domain length. In this appendix, we introduce the analytical passages to calculating the first and second shape derivatives of the energy functional and its inverse, the symbolic calculus applied on FEniCSx and a numerical application to a problem with known analytical solution.

A.1 Change of variables

In order to calculate the shape derivatives of the potential energy function P, we introduce a change of variables, allowing us to describe the stretch of the domain along the abscissa's axis. Let's consider a reference body Ω, where we use a vector field ξ → θ(ξ) continuously differentiable describing the virtual extension, we introduce a change of variables s.t.:

Where h is the increment translating vector θ(ξ). We express the gradient of the transformation (A.1) and its inverse following Einstein's notation in the new set of coordinates and obtain:

To express the potential energy in the new coordinates, we need to express the strain tensor ε ij and the determinant in the new coordinates. As ε ij is a function of G jk = F -1 jk , we look for its expression first:

We derive once Equation (A.3), with respect to variable h. Here henceforth derivatives with respect to h will be denoted:

We derive once Equation (A.4), with respect to variable h:

By evaluating (A.4) and (A.5) at h = 0 we obtain:

Finally, we write the new expression of the volume element:

The strain tensor writes in the new set of coordinates:

Finally, the potential energy P h writes in the new domain:

where w( 1 2 u i,k G kj + 1 2 u j,k G ki ) denotes the elastic energy density.

A.2 A Gθ-like derivative: first order

The first directional derivative of Equation (A.10) can be found deriving the potential energy with respect to h and evaluating it at h = 0.

This derivative can be seen as an application of the Gθ method, introduced by Destuynder et al. [1984], with the difference of the θ function. In the Gθ method, θ is a function used to calculate the virtual advancing of a crack.

Thanks to the equilibrium condition (A.12) and expressions (A.6), Equation A.11

like derivative: second order

The second derivative can be found deriving twice with respect to h equation (A.10) and evaluating the expression at h = 0, doing so we obtain:

Thanks to the equilibrium condition,

dh 2 cancels out. However, first order derivatives of u, du i,k dh cannot be simplified. The complexity relies in solving this derivative numerically. We show a possible application in the following section by recasting the problem through symbolic calculus which can be easily implemented in FEniCS and FEniCSx.

A.4 Symbolic calculus of domain derivatives

We denote with P the potential energy. Any subscript denotes the directional derivative with respect to the sub-scripted variable and all directions are denoted with a "hat" symbol.

The potential energy is re-written as a function of u and θ, like so:

One delicate aspect is that the displacement field u depends on θ as well. Equation A.11 can then be written as:

And equation A.14 as:

(A.17) We therefore need the solution uθ . The problem can be reformulated starting from the equilibrium equation:

Find

This condition should be valid ∀θ s.t.

Calculus of shape derivatives

This is a linear system that allows us to compute uθ = u θ (θ)( θ) as solution of

Simplifications Some simplifications due to the equilibrium condition at order one and two can be made:

Finally, Equation A.16 becomes: .23) and Equation A.17 becomes :

In the composite bar that we are modelling, we choose a θ that corresponds to a longitudinal stretch of the bar, hence:

A.5 Plane Stress Cylinder -Validation of G-theta Method at order 2

The purpose of this section is to look at a simple elastic problem where the analytical solution of the displacement field is known, so that its first and second derivative with respect to the external radius can be calculated analytically. Doing so, we can verify the validity of the numerical solutions of the G-θ like approach at first and second order.

We will be considering a linear elastic cylinder with an internal radius R in , and external radius R out . The cylinder is subjected to an internal pressure P in and an external pressure P out . The material properties depend on the Young modulus E and the Poisson ration ν. Thanks to the structure's symmetry, the 2D cross-section corresponds to a quarter of a hollow cylinder under plane-stress assumptions. For axisymmetric and plane-stress conditions, the displacement field is radial u = u(r)er. The problem can be simplified thanks to the loading, geometric and material symmetries.

1. Loading symmetry: the loading is applied according to the radial axis 2. Geometric symmetry: symmetrical according to the radial axis 3. Material symmetry: the material used is linear, homogeneous and isotropic

A.5.1 Variational Formulation

We hereby consider a quasi-static problem that can be solved with the equilibrium equation and its respective boundary conditions, see system A.25.

Equilibrium:

The potential energy of the problem writes:

We specify the affine space of the trial functions: Ĉ0 = { v(x) regular so that v(0) = 0} The variational formulation yields to Equation A.27.

The non dimensional form of the potential energy is expressed in equation (A.29), thanks to the normalized parameters introduced herebelow: Appendix A. Calculus of shape derivatives

A.5.2 Analytical Solution

Thanks to the symmetries the displacement field is radial: u = u(r)er , the strain and stress tensors are further simplified.

The displacement solution: u(r) = Ar + B r where A and B are constants (A.30)

For a practical purpose we rename the following variable: P int = p a , P out = p b , R int = a, R out = b. The stain and strain and stress are defined in the the systems A.31 and A.32.: ε = u r,r (er er) + u r r (e θ e θ ) σ = σ rr (er er) + σ θθ (e θ e θ )

Since σ rr (r = a) = -p a and σ rr (r = b) = -p b , we find the constants A and B:

Finally, the displacement solution and strains simplify to: Where, F, I, ∇θ are second order tensors indicating respectively the gradient of the transformation, an identity matrix and the gradient of θ. N is a vector and n the normal vector to the external surface. For a circular crown, we define the stretching function θ(x, y) to be the following one: θ(x, y) = ( r-R int R out -R int cos(α) er, r-R int R out -R int sin(α) e θ ), r = x 2 + y 2 , α = arctan( y x )

If we apply the change of coordinates to Equation A.29, we obtain the following Equation: