Silvia Regalado Olivares

Roberto López López

Manuel Collados

Carlos Quintero Noda

Maria Jesús Martínez González

Sara Esteban Pozuelo

Bernard Gelly

Optical design of an image slicer-based Integral Field Unit for the THEMIS solar telescope

Keywords: IFU, IFS, Image Slicer, Spectroscopy, Spectro-polarimetry, Solar Observations, Telescopes, Diffraction, Optical design, Spectrograph, Multi-slit

THEMIS, a 90 cm solar telescope, offers spectropolarimetric observations in the VIS and NIR spectral ranges through a long-slit spectrograph. Now, an image slicer-based Integral Field Unit (IFU) is being developed as an evolution of the long-slit set-up. The IFU will allow observing a 2D Field-of-View (FoV) and several spectral regions simultaneously. Integral Field Spectroscopy mitigates very well the loss of information generated by the atmospheric seeing over a single long-slit and is more suitable for 2D evolutionary studies of solar phenomena. The IFU optical design inherits the concept of the IFU currently installed on the GRIS spectrograph at the GREGOR telescope and has been adapted to THEMIS. Also, we took advantage from the image slicers manufacturing process evolution since the GRIS IFU was created. First, the spatial resolution is increased by reducing the width of the slicer mirrors from 100 to 80 µm. Second, the FoV is increased by changing the slices number from 8 to 16. In terms of optomechanical integration, the slicer mirror unit is placed at the telescope focal plane. The available volume is limited and has constrained the IFU elements position and the mirrors focal length. The design strategy generates 16 different pupil planes that are superimposed and keeps the optical path constant to allow interchangeability between the new IFU and the long-slit mode. The most challenging condition was keeping the pupil plane position unaltered since the input beam is not telecentric. This contribution presents the IFU optical design and describes the challenges faced in design phase.

INTRODUCTION

THEMIS, 1 the 90-cm optical solar telescope, offers spectropolarimetric observations in the VIS and NIR spectral ranges through a long-slit spectrograph. Currently, the possibility of updating THEMIS with advanced instrumentation is being explored. In this paper, we present a prototype of a slicer-based Integral Field Unit (IFU) designed for THEMIS. This IFU system will offer the simultaneous acquisition of 2D images in several spectral regions. Moreover, an advantage of this design is its interchangeability between long-slit and IFU modes. This IFU optical design heritages the concept of the IFU that is currently installed on the GRIS spectrograph 2 at the GREGOR solar telescope and is also based on image slicers. [START_REF] Dominguez-Tagle | First light of the integral field unit of gris on the gregor solar telescope[END_REF] The IFU system designed for THEMIS has been adapted to the specifications of this telescope. Since the image slicer's manufacturing process has evolved since the GRIS IFU was created, presently it is possible to produce a system with better performance. First, the reduction of the width of each mirror slice from 100 µm to 80 µm has led to an increase of the spatial resolution. Second, the number of slices grows from 8 to 16, augmenting the FoV covered by the new Integral Field System (IFS).

The long-slit is placed after a transfer optics, at the focal plane of the telescope, F2. Then, the beam reaches the pre-disperser of the THEMIS spectrograph 4 but first, at 210 mm just after the slit, a field lens reimages the pupil providing telecentrism to the predisperser since the input is not telecentric. At IFU mode, the Integral Field Unit will substitute the long-slit.

One of the limitations of the optical design is the reduced space available for placing the IFU. The requirement to fit the IFU in this cavity constrains the mirror focal length and the element distribution of the system. The design of the IFU system maintains the optical path invariant in comparison to the long-slit mode and also avoids placing elements outside the cavity. Thus, the presented design allows a straightforward interchangeability between long-slit and IFU modes. To compensate for optical path differences, a final extra mirror sends the light to the spectrograph. The pupil position and size after the IFU have been maintained because the lens that provides telecentrism cannot be removed. Since the management and alignment of all pupil planes formed in the IFU system are crucial, working in non-telecentric conditions for the IFU pupil management and overlapping has posed a challenge. As a solution, the IFU provides a telecentric exit pupil and then, a final mirror reimages the pupil at the corresponding position. In the next section, the optical design of the IFU is described.

OPTICAL DESIGN DESCRIPTION

Starting requirements

The next table shows the starting requirements for the optical design. The IFU will be placed at the focal plane of the telescope, F2. At that level, the input pupil is at a distance of 2626.26 mm from F2. Then, there is a lens, called LCS, that provides a telecentric pupil as an input for the predisperser. This lens is kept in its place on IFU mode, so, for reproducing identical pupil conditions, the IFU must provide an exit pupil at the same distance from LCS with the same diameter as when the long-slit is used. The first element of the predisperser is the collimator which forms a pupil at 6298.87 mm with diameter size of 126 mm.

General

The image slicer is composed by 16 flat mirrors of 80 µm width and 1.73 mm length each, forming 16 mini-slits. At the image plane, all the mini-slits are reorganised as one long-slit with no magnification on the mini-slits. The elements that form the system can be seen on the side and top view in Figure 2a and Figure 2b. The image-slicer is the first element of the IFU and is placed at the focal plane F2. At the slicer, the FoV is divided into sixteen sub-fields. Each slicer mirror sends the light of each sub-field to a collimator mirror with 85 mm of focal length. The collimator mirrors are placed at their focal length from the slicer and at 11 mm from the vertical Y-axis, which ensures that they are not blocking the light from the telescope that reach the slicer mirrors. Figure 3 displays the front view of the collimator mirrors with their configuration number from the OpticStudio software. As can be seen, the collimator mirrors are arranged in two rows, where they are alternated. In order to simplify the pupil management at all stages, the sixteen image pupils formed by the collimator mirrors have been overlapped at the same plane generating a single IFU pupil. Figure 2 shows the location of the IFU pupil, which is at 88 mm after the collimator mirrors and at 7.6 mm on the Y-axis. In Figure 4, the footprint of the IFU shows the good overlapping of all the sixteen pupils. At this level, a mask can be placed to avoid scattered light. The pupil diameter is 1.52 mm. After the IFU pupil, the sixteen collimated beams reach the next element, the camera mirrors array. These mirrors have the same focal length as the collimator mirrors, leading to a 1:1 magnification in the system. They also have an alternate distribution, oppositely arranged in comparison to the collimator mirrors distribution since the beams cross each other at the pupil plane. With the alternate distribution, the individual alignment of the output mini-slits can be optimised (this issue is described later in this section). The camera mirrors are located at 85 mm after the IFU pupil, at their focal distance, giving a telecentric pupil as a result, which allows the alignment of the sixteen IFU pupils. To get the required pupil position, a Pupil Compensator Mirror (PCM, described later in this section) has been added.

Optical design solution

Each camera mirror redirects the beam down to an individual mirror of the folding mirrors array, which, in turn, sends the beam down to the Y-axis and also compensate for the differences in focal plane location for each mini-slit at the output plane. Figure 5 shows a rear view of the folding mirrors where the individual displacement can be seen. After the folding mirrors, the IFU first image plane is generated. It consists of sixteen mini-slits arranged in a single continuous line as can be seen in the general overview portrayed in Figure 6. Figure 7 shows a close view of the IFU first output mini-slits. Finally, the PCM mirror is the last optical element. It is a monolithic mirror with a 3552.76 mm radius of curvature whose functions are diverse. On one hand, it redirects the light to the predisperser without deviations from the beam coming from the telescope. On the other hand, it generates a virtual exit pupil simulating the same conditions in location and diameter as when the long-slit is placed. Additionally, the PCM mirror reimages the mini-slits providing the final IFU image plane, the entrance to the predisperser. The image is virtually placed at the exact position of the image-slicer compensating then for the optical path differences.

In Figure 8, a side general view of the last elements of the system is displayed. Figure 9 shows a top view of the PCM mirror and the final mini-slits. Due to the curvature of the PCM mirror, the final output mini-slits have a small decentering on the Y-axis from their theoretical horizontal axis, which is more significant at the edges where it reaches up to 44.8 µm.

Furthermore, one of the reasons for the disposal of the elements on the system is to avoid the mini-slit tilt. This issue was also considered for the design of the 35 µm image-slicer IFU for GREGOR, [START_REF] Regalado Olivares | Optical and mechanical updates in the GREGOR Infrarred Spectrograph for simultaneous spectral observations[END_REF] which is already in development status. Ideally, the mini-slits should be aligned with their theoretical horizontal axis. However, they are tilted because the inclined optical elements of the system affect the mini-slit alignment. The design has been optimized in order to obtain a reasonable tilt, leading to a maximum value of 0.8 • .

Final analysis

At the final image plane, the image quality has been analyzed. Figure 10 displays the spot diagram for the sixteen configurations for the 0.850 µm wavelength at nine FoV. showing that the system is diffraction-limited. Figure 11a shows the footprint of the sixteen pupils overlapped after the PCM mirror, which has the same dimensions as the input pupil. Figure 11b displays the obtained footprint for the overlapped pupil after the collimator mirror of the predisperser, which shows also the same size as when the long-slit mode is used. The good overlapping of the sixteen pupils is confirmed.

CONCLUSIONS

A preliminary design has been presented for an Integral Field Unit based on image slicers for the THEMIS spectrograph. This IFU fits into the available free space and fulfills the initial requirements. This system allows an easy interchangeability between the long-slit and IFU modes, while keeping the optical path length and the pupil size and position unaltered at the entrance of the predisperser.

Figure 1

 1 Figure 1 shows a 3D model of the IFU general layout

Figure 2 :

 2 Figure 2: IFU Side view (a) and top view (b).

Figure 3 :

 3 Figure 3: Rear view of the image-slicer and front view of the collimator mirrors array with their configuration number from OpticStudio.

Figure 4 :

 4 Figure 4: Footprint of the IFU pupil from OpticStudio

Figure 5 :

 5 Figure 5: Rear view of folding mirrors array, the IFU pupil and front view of the camera mirrors array.

Figure 6 :

 6 Figure 6: General view of the camera mirrors, folding mirrors, IFU first output mini-slits and PCM.

Figure 7 :

 7 Figure 7: Folding mirrors array and IFU first output mini-slits.

Figure 8 :

 8 Figure 8: Side view of the last elements of the IFU.

Figure 9 :

 9 Figure 9: Top view of the PCM mirror and the final output virtual mini-slits.

Figure 10 :

 10 Figure 10: Spot diagram for the 0.850 µm at the final output slit for the sixteen configuration at nine FoV.

Figure 11 :

 11 Figure 11: Footprint of the sixteen overlapped pupils after the PCM mirror (a) and after the collimator mirror of the predisperser (b).

Table 1 :

 1 Starting requirements for the THEMIS IFU optical design

	Magnification	1
	Wavelength range	0.550 µm -0.850 µm
	Illumination type	Not Telecentric
	Entrance pupil distance to F2	2626.26 mm
	Entrance pupil diameter	46 mm
	Image-slicer
	Number of slicers mirrors	16
	Slices width	80 µm
	Slices length	1.73 mm
	Total slicer unit size	1.73 mm x 1.28 mm
	Output slits
	Number of total output slits	1
	Slits arrangement	Aligned and continuous
	(Non-illuminated) Distance in-between mini-slits	0.273 mm

ACKNOWLEDGMENTS

This project/development has been funded by the European Union's Horizon 2020 Research and Innovation programme under Grant Agreement No. 824135 (SOLARNET). SEP acknowledges the funding received from the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation programme (ERC Advanced Grant agreement No. 742265).


H Girard
email: hugo.girard@insa-lyon.fr

A Doitrand

B Koohbor

R G Rinaldi

N Godin

D Long

J Bikard

Influence of nearby fiber on fiber-matrix debonding: Coupled Criterion prediction and debonding shape determination

Keywords: Fiber-matrix debonding, Finite Fracture Mechanics, Coupled Criterion, Linear Elastic Fracture Mechanics, Neighboring fiber

Fiber-matrix interface debonding in two-fiber specimens under remote tensile loading is studied both experimentally and numerically by means of a coupled stress and energy criterion. Depending on its relative position, the neighboring fiber induces a perturbation of both stress and energy fields at the reference fiber interface which results in asymmetrical debonding initiation and propagation. The determination of the debonding initiation and propagation shape is addressed based on either i) stress isocontours, ii) energy isocontours or iii) the Coupled Criterion (CC). It was found that the debonding initiation configuration can be determined based on stress (respectively energy) isocontours for small (respectively large) enough interface brittleness number. For intermediate brittleness number, the debonding initiation configuration cannot be obtained using neither the stress nor the energy isocontours, but requires a coupling of both aspects. Despite different initiation debonding configurations, the corresponding initiation remote stresses do not differ much, which results in similar debonding configurations after unstable crack propagation following initiation.

Introduction

Interface fracture is a challenging problem to address from both experimental and numerical points of view. First, the very definition of an interface may not be straightforward. In numerical models, it is generally defined either as a zero-thickness surface that links two different materials or possibly considering a finite thickness when its dimension is comparable or not negligible with respect to the other material dimensions. Depending on the materials involved and the manufacturing process, the transition from one material to another might not be just a surface separating two homogeneous and uniform materials, but a transition zone of a finite thickness may also be present [START_REF] Mäder | Interphase characterization in polymer composites -monitoring of interphasial behavior in dependence on the mode of loading[END_REF][START_REF] Karger-Kocsis | Recent advances in fiber/matrix interphase engineering for polymer composites[END_REF][START_REF] Fritah | Influence of Water on an Epoxy/Amine-Metal Interphase: A Combined DFT and Mixing Calorimetry Approach[END_REF]. Then, interface fracture is another challenge as it is generally described in models as the separation of the two lips of the surface between both materials, including or not a process zone. In experiments, either adhesive, cohesive, or a combined adhesive-cohesive cracking may be encountered depending on the interface bonds and homogeneity (possible presence of local defects for instance). As a consequence, the viewpoint of interface fracture modeling can be regarded as a homogenized vision that gathers all these mechanisms. The direct consequence is that interface fracture properties (e.g., strength, fracture toughness) may not be straightforward. Nevertheless, numerical interface fracture models such as, e.g., Cohesive Zone Models (CZM) [START_REF] Barenblatt | The formation of equilibrium cracks during brittle fracture. General ideas and hypotheses. Axially-symmetric cracks[END_REF][START_REF] Dugdale | Yielding of steel sheets containing slits[END_REF] or the Coupled Criterion (CC) [START_REF] Leguillon | Strength or toughness? A criterion for crack onset at a notch[END_REF] are useful tools to describe material behavior observed experimentally. For instance, the CC enables assessing debonding initiation at the fiber-matrix interface. This approach was introduced by Leguillon [START_REF] Leguillon | Strength or toughness? A criterion for crack onset at a notch[END_REF] to extend Linear Elastic Fracture Mechanics (LEFM) [START_REF] Griffith | The phenomena of rupture and flow in solids[END_REF] to crack initiation.

Mantič [START_REF] Mantič | Interface crack onset at a circular cylindrical inclusion under a remote transverse tension. Application of a coupled stress and energy criterion[END_REF] first applied the CC to assess debonding at the interface of a circular inclusion embedded within an elastic matrix subjected to a far field tensile loading. The Hutchinson and Suo [START_REF] Hutchinson | Mixed Mode Cracking in Layered Materials[END_REF] relation was used to describe the interface mode mixity leading to the initiation of debonding triggered either by both stress and energy conditions or by the energy condition exclusively. García et al. [START_REF] García | Debonding at the fibre-matrix interface under remote transverse tension. One debond or two symmetric debonds?[END_REF] extended the embedded inclusion problem to symmetric debonding at the two opposite poles and showed that non-symmetric initiation was more favorable than symmetric crack initiation due to a larger energy released per crack surface unit. Several extensions of this approach were then proposed for bi-axial loading, for instance, by Mantič and García [START_REF] Mantič | Crack onset and growth at the fibre-matrix interface under a remote biaxial transverse load. Application of a coupled stress and energy criterion[END_REF] and Muñoz-Reja et al. [START_REF] Muñoz-Reja | Crack onset and propagation at fibre-matrix elastic interfaces under biaxial loading using finite fracture mechanics[END_REF] as well as comparison with other numerical approaches such as Cohesive Zone Model (CZM) [START_REF] Gentieu | Size effect in particle debonding: Comparisons between finite fracture mechanics and cohesive zone model[END_REF][START_REF] García | Fiber-size effects on the onset of fiber-matrix debonding under transverse tension: A comparison between cohesive zone and finite fracture mechanics models[END_REF][START_REF] Távara | Revisiting the Problem of Debond Initiation at Fibre-Matrix Interface under Transversal Biaxial Loads -A Comparison of Several Non-Classical Fracture Mechanics Approaches[END_REF]. Besides, the CC can be an efficient method to perform an inverse identification of the interface properties, as demonstrated by Doitrand et al. [START_REF] Doitrand | Determination of interface fracture properties by micro-and macro-scale experiments in nacre-like alumina[END_REF]. A range of interface shear and opening critical energy release rates and strengths were also determined by Girard et al. [START_REF] Girard | Numerical simulation of fiber-matrix debonding: Inverse identification of interface properties[END_REF] for a sample with a single glass fiber embedded into an epoxy matrix.

Previous studies have mainly focused on a single fiber in an infinite matrix. However, a more realistic composite microstructure involves numerous fibers with varying relative positions (which can be for instance quantified by inter-center distances and angles with respect to the loading direction). Several authors have investigated the influence of neighboring fibers on the debonding process. Muñoz-Reja et al. [START_REF] Muñoz-Reja | Crack onset and propagation at fibre-matrix elastic interfaces under biaxial loading using finite fracture mechanics[END_REF][START_REF] Muñoz-Reja | Influence of a neighbour fibre on the onset and growth of a fibre-matrix debond under biaxial loading. A study by Finite Fracture Mechanics at linear elastic interfaces[END_REF] assessed debonding initiation in a two-fiber configuration under tensile loading, varying the inter-center distance between 2.25 and 2.5 times the fiber radius and compared two different inter-center angles. Asymmetric debonding locations with respect to the loading direction were obtained. It was shown that the smaller the inter-center distance, the smaller the remote tensile load required for the debonding to initiate. Sandino et al. [START_REF] Sandino | Numerical analysis of the influence of a nearby fibre on the interface crack growth in composites under transverse tensile load[END_REF] studied the influence of angle and distance between fibers center on debonding propagation. Intercenter angle influenced the debonding arrest angle and remote loading required for propagation. Moreover, the influence of the nearby fiber vanishes after a certain distance. Within the same framework, Velasco et al. [START_REF] Velasco | Interaction between fibres in the transverse damage in composites[END_REF] and Zhuang et al. [START_REF] Zhuang | Transverse crack formation in unidirectional composites by linking of fibre/matrix debond cracks[END_REF] proposed a similar analysis in the presence of a pre-existing debonding and found that the debonding of another fiber may either be delayed or favored by the pre-existing debonding depending on its location. It was also found that the debonding kink into the matrix was affected by the neighboring fiber and was either damped or reinforced depending on whether there was a pre-existing debonding. Kushch et al. [START_REF] Kushch | Numerical simulation of progressive debonding in fiber reinforced composite under transverse loading[END_REF] used CZM for the prediction of debonding initiation with a two-fiber configuration and focused on interface stresses and debonding locations that vary with the inter-center angle, either at the top or bottom poles. However, most of the aforementioned studies do not focus on a configuration that is asymmetric with respect to the loading direction, so that the debonding location and propagation remain symmetric. Besides, there is a lack of experimental confrontation to validate the numerical approach employed.

The asymmetric debonding shape leads to some difficulties in the implementation of the CC which requires the a priori definition of the crack shape. Such shape of the debonding can be prescribed symmetric in the case of a single inclusion [START_REF] Mantič | Interface crack onset at a circular cylindrical inclusion under a remote transverse tension. Application of a coupled stress and energy criterion[END_REF][START_REF] París | Kinking of Transversal Interface Cracks Between Fiber and Matrix[END_REF]. However, multi-fiber configurations induce a perturbation of the fields at the interface and the debonding no longer remains symmetric with respect to the loading direction [START_REF] Muñoz-Reja | Influence of a neighbour fibre on the onset and growth of a fibre-matrix debond under biaxial loading. A study by Finite Fracture Mechanics at linear elastic interfaces[END_REF][START_REF] Kushch | Numerical simulation of progressive debonding in fiber reinforced composite under transverse loading[END_REF]. To overcome this shape identification problem, some authors estimated the crack shape based on stress isocontours [START_REF] Doitrand | 3D application of the coupled criterion to crack initiation prediction in epoxy/aluminum specimens under four point bending[END_REF][START_REF] Doitrand | Comparison between 2D and 3D applications of the coupled criterion to crack initiation prediction in scarf adhesive joints[END_REF][START_REF] Carrere | Theoretical study based on 2D assumptions of the influence of small pores on crack initiation in adhesively bonded joints[END_REF]. This solution led to an efficient way of determining the crack shape. However, it does not take into account any energy aspects which, in some cases, can drive the initiation of debonding [START_REF] Mantič | Interface crack onset at a circular cylindrical inclusion under a remote transverse tension. Application of a coupled stress and energy criterion[END_REF][START_REF] Girard | Numerical simulation of fiber-matrix debonding: Inverse identification of interface properties[END_REF][START_REF] Doitrand | Damage onset modeling in woven composites based on a coupled stress and energy criterion[END_REF][START_REF] Martin | Competition between deflection and penetration at an interface in the vicinity of a main crack[END_REF]. In this regard, García et al. [START_REF] García | A numerical study of transverse cracking in cross-ply laminates by 3D finite fracture mechanics[END_REF] and Doitrand et al. [START_REF] Doitrand | Damage onset modeling in woven composites based on a coupled stress and energy criterion[END_REF] proposed an iterative approach by comparing the CC solution based on a parameterized crack geometry in the vicinity of the stress concentration which, therefore, considers the energy aspects. The two above-mentioned approaches for crack shape determination do not consider all possible crack geometries. As a consequence, it is yet not established which initiation crack shape among all possible crack shapes is the most favorable, i.e. results in the minimum imposed loading at initiation.

The present study aims at determining the CC solution in the case of an asymmetric debonding by simultaneously considering all possible crack shape configurations. The configuration under investigation and related experiments are presented in Section 2. The remote loading required for a debonding to initiate is therefore evaluated with nearby two-fiber specimens having different inter-center angles and distances, as detailed in Section 3. In the following, the CC solution is compared to the stress-and energy-based debonding shape solutions to quantify the accuracy of the stress isocontours assumption in terms of debonding angles and remote loading. The CC solutions are then confronted to experimental observations [START_REF] Livingston | Characterizing fiber-matrix debond and fiber interaction mechanisms by full-field measurements[END_REF] in Section 4. In details, the location and size of the debonding are compared, as well as the remote loading required for the debonding to initiate.

Experiments

Sample geometries and testing protocol

Two-fiber specimens under uniaxial tensile loading are considered to study the influence of a nearby fiber on fiber-matrix debonding initiation and propagation. Experimental results are taken from [START_REF] Livingston | Characterizing fiber-matrix debond and fiber interaction mechanisms by full-field measurements[END_REF][START_REF] Girard | Experimental characterization of inplane debonding at fiber-matrix interface using single glass macro fiber samples[END_REF]. Figure 1 describes specimen geometry and dimensions, sized accordingly to ASTM D638 with a sample thickness of 6 mm. Thirty samples are subjected to uniaxial tensile loading along the (Oy) direction, leading to a far-field stress (σ ∞ ). The inter-center angle (α) from the (Ox) direction and distance (d) are varied in the range 0, 30, 45, 60 and 90 degrees and 2.5r, 4r and 5r, respectively, with r the fiber radius equal to 1 mm. It is worth noting that both fiber absolute positions do change between the configurations as they are rotated with respect to the inter-center midpoint. Samples are manufactured by embedding two glass macro-fibers in a room temperature curable thermoset epoxy resin whose isotropic elastic properties are detailed in Table 1.

In the following, the fibers are referred to as respectively reference and neighboring fibers and all quantities will be referring to the reference fiber in the sequel (see Figure 1). The angular position at the interface of the reference fiber is denoted θ f . In order to study the field perturbation near the fiber and the debonding process, a speckle pattern was applied on the samples surfaces to enable performing 2D Digital Image Correlation (DIC). A camera equipped with a high magnification

Property

Epoxy

Glass Fiber Young Modulus [GPa] 2.36 ± 0.10 63 Poisson's ratio [-] 0.40 ± 0.04 0.30

Table 1: Linear elastic properties of the epoxy matrix and glass fiber. Fiber properties are provided by the manufacturer and matrix properties are based on 3 in-house measurements from [START_REF] Livingston | Characterizing fiber-matrix debond and fiber interaction mechanisms by full-field measurements[END_REF].

lens was used to capture debonding initiation and propagation in situ. The camera recording was synchronized with the load cell for a similar extraction time. Figure 1: Specimen geometry and dimensions including two fibers (radius r), where d is the inter-center distance and α is the angle between the (Ox) axis and the direction passing through both fiber centers. The reference fiber interface is divided into a top half circle (T) and a bottom half circle (B) where θ f denotes the angular position.

Experimental measurements

The numerical simulation of fiber-matrix debonding initiation and propagation presented in the following enables determining the following physical quantities:

• Debonding initiation location (top half circle, bottom half circle or both fiber half circles);

• Range of debonding initiation angles along the interface (θ d );

• Remote applied stress (σ ∞ ) associated with debonding initiation.

On the one hand, the extraction of the debonding locations and angles is based on DIC residuals [START_REF] Tomičević | Evaluation of fatigue crack network growth in cast iron for different biaxial loading paths via full-field measurements[END_REF]. The debonding leads to the appearance of a new black area at the fiber-matrix interface. Since the emerging black pixel area was not present in the reference image, high correlation residuals can be detected and thus related to the debonding location and size. On the other hand, the determination of the debonding initiation loading is based on previous work by Livingston et al. [START_REF] Livingston | Characterizing fiber-matrix debond and fiber interaction mechanisms by full-field measurements[END_REF] and Girard et al. [START_REF] Girard | Experimental characterization of inplane debonding at fiber-matrix interface using single glass macro fiber samples[END_REF] where virtual gauges are placed across the fiber-matrix interface for local displacement and strain measurements. The DIC solution is used to extract the gauge length variation as a function of the applied loading. Debonding initiation thus results in a sudden change in slope which allows the determination of the far-field load (and stress by simply dividing the load by the pristine transverse cross section area). All the experimental results obtained for each configuration are detailed in Section 4.

Debonding shape prediction

Coupled Criterion implementation

The initiation of debonding is predicted using the CC. A 2D plane strain Finite Element (FE) model is implemented under the assumption of linear elasticity and small deformations. As a consequence, the stress is proportional to the applied load whereas the elastic strain energy is proportional to the square applied load. In the sequel, the calculations were performed with a σ ∞ = 6.85 MPa remote stress. The FE model consists of two 2 mm diameter fibers embedded in the matrix. The material properties presented in Table 1 are assigned to each component. A 0.02 mm mesh size is chosen at the fiber-matrix interface, which ensures that the influence on the stress and elastic strain energy is smaller than 1% for a finer mesh. The geometry of the model is limited to 50 mm and 13 mm along the longitudinal (Oy) and transverse (Ox) directions with limited nodal displacements so that the configuration is close to the experimental setup.

Debonding initiation is assessed using the CC, which combines two conditions to provide the initiation loading and debonding angle range. On one hand, the stress condition must be satisfied at any location over the entire debonding path Γ before initiation. As a consequence, the stress requirement is not a point stress but a nonlocal condition since it must be fulfilled all over the debonding path before initiation. It can thus be employed even in the case of a nonuniform stress field. On the other hand, the Incremental Energy Release Rate (IERR) must be larger than the interface average critical Energy Release Rate (ERR). These two conditions are summarized in Equation [START_REF] Mäder | Interphase characterization in polymer composites -monitoring of interphasial behavior in dependence on the mode of loading[END_REF].

σ eq = σ 2 nn + 1 µ 2 τ 2 nt ≥ σ c ∀ x ∈ Γ, G inc (θ d , σ ∞ ) ≥ G c (θ d ).
(

) 1
The parameter µ denotes the τ c to σ c ratio, which are respectively the shear and tensile strengths. The IERR is denoted by G inc and G c is the average critical ERR. The normal and shear stresses of the interface are σ nn and τ nt , respectively. The calculation of all these quantities is detailed in the sequel. Contrary to the CC application to single-fiber configuration where the initiation debonding location is a priori known and symmetric with respect to (Oy) direction [START_REF] Mantič | Interface crack onset at a circular cylindrical inclusion under a remote transverse tension. Application of a coupled stress and energy criterion[END_REF][START_REF] García | Debonding at the fibre-matrix interface under remote transverse tension. One debond or two symmetric debonds?[END_REF][START_REF] Girard | Numerical simulation of fiber-matrix debonding: Inverse identification of interface properties[END_REF], the initiation debonding configuration may now depend on the nearby fiber (intercenter distance and angle, see Figure 1). As a consequence, solving the CC requires to determine not only the remote imposed stress causing debonding initiation but also the initiation debonding configuration, which is addressed in next sections.

Methodology

Previous studies have shown that for multi-fiber samples, the interface stress fields are influenced by their respective neighboring fibers [START_REF] Sandino | Numerical analysis of the influence of a nearby fibre on the interface crack growth in composites under transverse tensile load[END_REF][START_REF] Velasco | Interaction between fibres in the transverse damage in composites[END_REF][START_REF] Kushch | Numerical simulation of progressive debonding in fiber reinforced composite under transverse loading[END_REF]. They affect both the stress level at the interface and the elastic strain energy released during debonding initiation. Consequently, the debonding trajectory can no longer be considered symmetric with respect to the loading configuration and must now be determined.

Based on the CC, debonding can initiate provided the stress and energy criteria given in Equation 1 are fulfilled. As a consequence, it is likely that the potential debonding configurations do not lie far from the ones maximizing either the equivalent stress or the elastic strain energy variation.

Figure 2 shows the variation of the equivalent stress (σ eq ) and elastic strain energy variation (-∆W ) computed as a function of the angular position along the fiber-matrix interface (θ f ). The equivalent stress variation as a function of θ f is obtained based on a calculation without debonding. The evaluation of -∆W is performed by releasing one interface node at a given angular position θ f . A constant mesh size along the interface is used to ensure that the angular location maximizing -∆W also maximizes the IERR. For a single-fiber configuration under remote tensile loading, the equivalent stress and -∆W profiles are symmetric and both their maxima are attained for the same angle (90 deg. angular position, see Figure 2a). Asymmetry in the equivalent stress and elastic strain energy variation nevertheless occurs for a two-fiber configuration, as shown in Figure 2b in a case where α = 45 deg. and d = 2.5r. More precisely, asymmetric fields prevail and both quantities reach a maximum at a similar 100 deg. angular position, different from the single fiber configuration. In the following θ i is used to defines the angular position θ f maximizing the stress or the elastic strain energy variation, where debonding could most likely initiate from. Table 2 lists the angles maximizing the strain energy released (θ -∆W i ) or the equivalent stress criterion (θ σeq/σc i ) for different inter-center angles at a fixed inter-center distance of 2.5r, distance for which the fields are most affected. The ratio between the two angles are presented in Table 2. All ratios remain close to 1 whatever the configuration meaning that both angles are close to each other, even for highly perturbed configurations. Single fiber 2: θ i variation for various inter-center angles where d = 2.5r and µ = 1.5, evaluated at the reference fiber top half circle. The upperscripts refer to the energy or stress conditions used for the angular position determination. Ratio between the two angles are added to highlight the small discrepancy between the two approaches.

-∆W 0 1 2 σ ∞ σ ∞ θ f σ eq /σ c [-] σ eq /σ c (
θ i Angular position (θ f ) [deg.] -∆W [N.mm] α = 45 deg. d = 2.5r -∆W 0.5 1 1.5 2 σ ∞ σ ∞ θ f σ eq /σ c [-] σ eq /σ c (b)
Figure 3 illustrates the field perturbation by pointing out the equivalent stress maximum location |θ i | for several fiber inter-center distances and angles. This angle is calculated at both top and bottom half circles of the reference fiber. The larger the difference from 90 degrees, the higher the perturbation of the fields with respect to the single fiber configuration. Noticeably, the stress fields at the fiber bottom half circle are only slightly influenced since the angles remain within the 90 deg. vicinity (+ or -3%) whatever the inter-center distance. On the contrary, at the top half circle of the fiber, the smaller the inter-center distance, the larger the perturbation. Additionally, θ i differs the most from 90 degrees for an inter-center angle close to 60 degrees. This configuration may thus result in strongly asymmetric crack initiation. In contrast, θ i reverts to 90 deg. when α = 90 deg. as the direction through the two fiber centers is parallel to the load direction. Starting from an initiation debonding angular position θ i , the debonding shape can then be described by increments of angles on each side of this initial location, denoted by ∆θ + and ∆θ -, as detailed in Figure 4. The sum of these two debonding increments corresponds to the total debonding angle θ d . In the sequel, we evaluate different approaches to determine the possible initiation debonding configurations, based on stress, energy or both stress and energy.

|θ i | [deg.] d = 2.5r B d = 3r B d = 3.5r B d = 4r B d = 2.5r T d = 3r T d = 3.5r T d = 4r T Single fiber

Debonding configuration based on stress criterion

A convenient crack shape determination approach consists in defining the possible crack configurations based on the stress isocontours [START_REF] Doitrand | 3D application of the coupled criterion to crack initiation prediction in epoxy/aluminum specimens under four point bending[END_REF][START_REF] Doitrand | Comparison between 2D and 3D applications of the coupled criterion to crack initiation prediction in scarf adhesive joints[END_REF][START_REF] Carrere | Theoretical study based on 2D assumptions of the influence of small pores on crack initiation in adhesively bonded joints[END_REF]. The stress isocontour-based crack definition consists in considering areas which strictly fulfill the stress criterion for a given imposed loading and specified interface strengths. As a consequence, it reverts to determining the debonding configuration described by θ i , ∆θ + and ∆θ -for which the stress criterion is fulfilled before initiation

σ ∞ σ ∞ θ i ∆θ - ∆θ + y x
Figure 4: Schematic of the debonding shape including the initiation debonding angular position θ i and debonding increments ∆θ + and ∆θ -on each side of the initiation location at the reference fiber top interface. The sum of both debonding increments correspond to the total debonding angle θ d .

(see Figure 4). The two debonding tips respectively correspond to angles θ i + ∆θ + and θ i -∆θ -. The advantage of this approach is that, for given shear and tensile strengths, the possible crack surfaces can be determined based on a single calculation without crack and that the stress criterion is strictly fulfilled for each possible crack. Figure 5 shows an example of the variation of the equivalent stress along the fiber-matrix interface at the reference fiber top half circle for a given sample geometry and varying strength ratios. Two configurations can be encountered depending on the tensile and shear strength magnitudes (that is a variation in the µ ratio):

• Tension-induced crack initiation leading to one equivalent stress maximum per half circle;

• Shear-induced crack initiation leading to two equivalent stress maxima per half circle.

More precisely, for µ larger than 1.5, the equivalent stress decreases monotonically from the maximum peak, inducing a single initiation site. On the contrary, for µ smaller than 1.5, the equivalent stress exhibits two local yet uneven maxima thus corresponding to two possible (yet non equiprobable) initiation sites. It is worth mentioning that experimental observations often show a single initiation site at one or both half circles of the fiber, which could indicate tension induced debonding initiation for the studied fiber-matrix interface. Nevertheless, we remind that without a high speed camera that would enable capturing either tension-or shear-induced crack initiation, it is not possible to draw a clear conclusion on this statement. Additionally, the many stochastic effects resulting from imperfections or heterogeneity could also influence debonding location. As a consequence, configurations with smaller µ, and thus shear-induced crack initiation, are also studied for comparison purposes.

Tension-induced crack initiation

In case of tension-induced crack initiation (µ > 1.5), the equivalent stress reaches a single maximum, which leads to a straightforward definition of the possible debonding configurations based on equivalent stress isocontours. Figure 6a shows equivalent stress to normal strength over the debonding angle corresponding to θ i + ∆θ + and θ i -∆θ -according to the stress condition introduced by Leguillon [6]. This curve is obtained as the minimum stress attained over the whole crack path prior to debonding. The equivalent stress to tensile strength ratio decreases with increasing both angle increments from θ i . The possible initiation debonding configurations based on the stress isocontours are depicted by the red line in Figure 6. The admissible initiation debonding angles correspond to configurations where σ eq is larger than σ c . Figure 6b shows the debonding shape based on the stress criterion isocontours for several values of µ larger than or equal to 1.5. ∆θ + TC indicates the transition angle between the areas undergoing tension or compression. Note that the range of µ presented in Figure 6b results in θ i differences smaller than 3 degrees. The debonding shape is thus slightly influenced by the interface shear to tensile strength ratio. Ultimately, the increase of µ results in a debonding shape close to that obtained with a normal stress criterion (corresponding to µ = ∞).

σ eq /σ c [-] α = 45 deg. d = 2.5r µ = 0.5 µ = 1 µ = 1.5 µ = 2
The average difference between the shapes for a normal stress criterion and the equivalent stress criterion for various µ values is evaluated. By comparing for fixed ∆θ + , the difference in ∆θ -between debonding shapes obtained using the equivalent stress (∆θ - σeq ) or the normal stress criterion (∆θ - σnn ) is evaluated using Equation 2.

diff =   1 ∆θ + TC ∆θ + TC 0 ∆θ - σeq -∆θ - σnn ∆θ - σnn   × 100 (2)
For all the configurations covered, the average difference turns out to be no larger than 6%. Only debonding angles located in the tensile zone are considered because a crack is considered not to initiate under compressive stress. All in all, whatever the shear to tensile strength ratio, the normal stress-based debonding configuration provides a good estimate of the stress-based debonding shape obtained based on the equivalent stress criterion. Importantly, the advantage of using the debonding shape based on the normal stress rather than on the equivalent stress is that it allows defining possible initiation configurations that do not depend on fracture properties. As a consequence, it is numerically more efficient for inverse identification approach for instance. The same study is performed for several angles and inter-center distances to ensure that this observation is valid for any fiber inter-center distances and angles and the maximum differences never exceeded 10%. It is worth mentioning that the difference vanishes for a symmetric configuration, i.e. intercenter angle equal to 90 degrees.

Shear-induced crack initiation

For µ smaller than 1.5, the equivalent stress shows two local maxima (Figure 5). This configuration could then lead to various initiation scenarios:

• Single debonding initiation at one stress maximum;

• Two debonding initiations at both maxima.

García et al. [START_REF] García | Debonding at the fibre-matrix interface under remote transverse tension. One debond or two symmetric debonds?[END_REF] showed that two simultaneous debonding initiations at both half circles of a single fiber were not favorable compared to a single initiation site. Similar approach is employed here to compare the IERR of a single or double initiation located on the same half circle of the fiber. The same isovalues are determined based on each stress maximum and the nodes are released on the corresponding angles. Table 3 The debonding configuration maximizing the IERR can be considered as the favorable initiation site. Thus, a single initiation site is energetically more favorable. Interestingly, it corresponds to the location where the equivalent stress is maximum (see Figure 5). Finally, in case of shearinduced crack initiation, possible debonding configurations can be determined following stress isocontours by considering only one debonding location that corresponds to the global equivalent stress maximum.

Debonding configuration based on energy criterion

The second approach consists in determining the debonding shape based on the configurations maximizing the IERR. To this end, nodes at the fiber-matrix interface are successively released on both sides of θ i in increments of ∆θ to compute the IERR for a given debonding angle (θ d ) using Equation 3, where W denotes the elastic strain energy of the model and r the fiber radius:

G inc (θ d , σ ∞ ) = W (0, σ ∞ ) -W (∆θ + + ∆θ -, σ ∞ ) r(∆θ -+ ∆θ + ) (3)
The local critical ERR is evaluated at both crack tips by means of the Hutchinson and Suo relation [START_REF] Hutchinson | Mixed Mode Cracking in Layered Materials[END_REF]. This relation links the interface opening (G IC ) and shear (G IIC ) critical energy release rates to the local critical ERR (G c ) using the mode mixity defined as ψ(θ) = arctan(|τ nt | /σ nn ) (Equation 4). The mode mixity is assessed at a distance of 0.015 mm from the two crack tips for each configuration according to [START_REF] Mantič | Interface crack onset at a circular cylindrical inclusion under a remote transverse tension. Application of a coupled stress and energy criterion[END_REF] and the maximum of the two G c is retained.

G c (ψ(θ d )) = G IC 1 + tan 2 [(1 -λ) ψ(θ d )] (4)
The parameter λ is linked to the ratio between G IC and G IIC . The average interface critical ERR, G c , corresponds to the average of the local interface critical ERR over the complete debonding shape, corresponding to θ d (Equation 5).

G c = 1 θ d θ i +∆θ + θ i -∆θ - G c (ψ(θ d ))dθ d (5)
Figure 7a shows the IERR to critical ERR ratio as a function of debonding increments on each side of the debonding initiation location. For a small debonding angle, the energy released is too small to satisfy the criterion. By increasing the debonding angle, the criterion is fulfilled and the surface exhibits a maximum, corresponding to an optimal angle maximizing the criterion. Therefore, this maximum corresponds to the favorable condition for the debonding to initiate, meaning that the loading to apply in order to fulfill the energy criterion is minimal. The debonding shape can be evaluated by maximizing the surface gradient until the maximum is reached, denoted by the blue marker in Figure 7.

The G IIC to G IC ratio, related to the parameter λ, influences the shape of the debonding as well as the configuration maximizing G inc /G c . Figure 7b shows the shapes obtained from three different values of λ for the studied configuration. The debonding initiation angle remains the same for a fixed value of λ. Setting λ = 1 results in G IC = G IIC and a smaller value of λ corresponds to G IC < G IIC . Decreasing λ leads to an increase in the critical ERR and G IIC for a fixed G IC and mode mixity. Therefore, no shape change is observed for a fixed value of λ as the ratio between G IC and G IIC remains the same.

Coupled stress and energy criterion

Coupling both stress and energy conditions allows the determination of the minimum loading to apply in order to initiate a debonding as well as the corresponding debonding angle. When both conditions are involved, the initiation debonding is sought as the angle for which both criteria are fulfilled for a minimum applied loading. As a consequence, it is likely that the initiation debonding configuration lies between the debonding shapes determined based on either the stress isocontours or maximizing the IERR. Figure 8 shows the applied remote stress required to fulfill both conditions (σ ∞ req ) for a range of angle increments. The loading required to fulfill the equivalent stress condition increases monotonically with increasing debonding angle whereas the one required to fulfill the energy condition reaches a minimum. Combining both criteria thus results in determining the debonding configurations and corresponding applied stress for which the CC is fulfilled. Among these debonding configurations, the initiation debonding configuration is determined as the one minimizing the applied stress.

The debonding configuration for which both criteria are fulfilled for a minimum loading (i.e. the initiation debonding configuration) is highlighted on Figure 8a and 8b with a green marker. The same approach is employed to determine the minimum required remote stress that follows either the stress-based or energy-based debonding configurations, also highlighted in Figure 8a and 8b. These minima differ slightly from the CC solution, both in terms of the required remote stress and debonding angle.

For given λ and µ, the debonding configuration minimizing the loading required to fulfill the CC actually depends on the interface brittleness number γ which was introduced by Mantič [START_REF] Mantič | Interface crack onset at a circular cylindrical inclusion under a remote transverse tension. Application of a coupled stress and energy criterion[END_REF] as follow (Equation 6):

γ = 1 σ c G IC E * r (6)
Each marker in Figure 8b corresponds to the CC solution with a different γ, where E * denotes the harmonic mean of the effective elastic moduli of the two constituents. For small enough γ, the debonding configuration obtained using the CC is close to the debonding shape predicted based on stress isocontours, which is also close the one determined by maximizing G inc /G c . How- ever, for large enough γ the debonding configuration retrieves the debonding shape determined by maximizing G inc /G c . Whatever the interface brittleness number, the debonding initiation configuration actually lies between the debonding shapes determined based on either stress isocontours or maximizing G inc /G c .

G inc /G c CC σ eq /σ c (b) θ i γ = 0.66 (c) θ i γ = 1.16 (d) θ i γ = 1.45 (e)
Figure 9 shows the stress-based, energy-based and CC-based debonding configurations obtained for inter-center angles of 90, 60 and 0 degrees with d = 2.5r. For the 90 degrees inter-center angle, a symmetric configuration is observed with respect to the loading direction. The stress-based, energybased and CC-based debonding shapes are similar since there is no asymmetric perturbation of the fields. Similarly to the previous example, for inter-center fiber angles different from 90 degrees (Fig. 9b and Fig. 9c), the debonding configuration predicted using the CC always lies between the stress-based and energy-based debonding configurations. Moreover, it is also consistent with the stress-based debonding shape for small γ and with the energy-based debonding shape for large γ. For intermediate γ, there is a transition zone for which the debonding configuration does not correspond to neither shapes. Figure 9b shows the predicted debonding configurations for a 60 deg. inter-center angle for which the transition from stress-based to energy-based occurs abruptly, i.e. all debonding configurations predicted by the CC are close either to the stress-based or the energy-based debonding configuration whatever γ. On the contrary, Figure 9c shows a smoother transition, where the debonding configuration predicted by the CC does not correspond to stressbased or energy-based debonding configuration for a significant range of γ (from 0.96 to 1.36). 2.

It is noteworthy that even in the case of a strongly influenced initiation location (when θ i is much larger or smaller than 90 deg.), the debonding initiation angle tends to recover a symmetric configuration with respect to the loading direction. Figure 9b highlights this last observation with a α = 60 deg. then the increment of angle ∆θ -is small compared to the increment of angle ∆θ + tending towards an initiation angle centered at about 90 degrees. The same observation was obtained by Mantič et al. who numerically [START_REF] Mantič | A linear elastic-brittle interface model: application for the onset and propagation of a fibre-matrix interface crack under biaxial transverse loads[END_REF] found that, wherever debonding initiation might occur, propagation will tend towards the loading direction as the amount of available elastic strain energy is larger. Stress isocontours allow efficient numerical prediction of the debonding shape because they can be derived from a purely linear elastic calculation without node release. The difference induced by the stress-based debonding shape from the favorable configuration can be evaluated. The stressbased solution is determined by locating the minimum on surface σ ∞ req /σ ∞ (Fig. 8a) that follows the stress-based debonding shape and directly comparing it to surface global minimum, corresponding to the favorable CC configuration. Figure 10a shows the variation of the relative difference in terms of remote loading between the two solutions as a function of γ for inter-center angles of 60 and 0 degrees. The relative difference shows increasing trends with γ. The maximum relative difference remains smaller than 5.1 % and the maximum is reached for the largest values of γ (γ = 1.54) when the configuration tends towards energy-based debonding shape. It is noticeable that the difference induced by the stress remains small and arguably not worth the computation time cost it requires. Therefore, basing the crack shape on the stress isocontours could be a good compromise to determine the debonding initiation remote loading. Figure 10b shows the location of the debonding based on the CC and the stress at initiation for α = 60 deg. and γ = 1.54. The debonding angles at initiation predicted by the two shapes are identical, but their locations differ slightly with respect to θ i . The location of the debonding based on the CC solution is equivalent to the energy solution in this case, tending to maximize the energy released by being symmetric with respect to the loading direction.

0.6 0.8 1 1.2 1.4 1.6 0 2 4 γ [-] Relative difference [%] d = 2.5r α = 60 deg. α = 0 deg. (a) σ ∞ σ ∞ σ eq θ d CC θ d θ i α = 60 deg. d = 2.5r γ = 1.54 (b)

Unstable crack propagation

After debonding initiation, it is possible to assess its propagation based on Griffith's criterion

G(θ d ) = -dW/dθ d ≥ G c (ψ(θ d ))
. The ERR G can actually be obtained from G inc based on the following relation (Equation 7), where S denotes the debonding surface i.e. rθ d or r (∆θ + + ∆θ -):

G(θ d ) = G inc (θ d ) + S × dG inc dS (7)
When

G inc /G c maximum is reached, d(G inc /G c )/dS = 0 so that G inc /G c = G/G c (see Appendix A)
. Therefore, two configurations can be encountered after the initiation of debonding. On the one hand, if d(G inc /G c )/dS ≥ 0, then G/G c is larger than 1 which implies further unstable debonding phase without increasing the applied load until G/G c = 1. On the other hand, d(G inc /G c )/dS = 0 for the initiation debonding configuration so that G inc /G c = G/G c = 1, resulting in no further propagation of debonding after initiation, as d/dS(G/G c ) < 0 (see Appendix B).

The debonding has two different crack ends, which may separately or simultaneously propagate. One approach to identify the unstable debonding shape is to consider the direction that maximizes the ERR to critical ERR ratio. This enables the determination of the arrest angle corresponding to an angle for which G/G c becomes smaller than 1.

Figure 11a shows the unstable propagation phase of the crack in the configuration driven by both stress and energy for γ = 1.1. The debonding propagates in an unstable manner from the debonding angle solution, at initiation, obtained using the CC to the angle where G/G c becomes smaller than 1, after unstable propagation phase. Figure 11b shows the unstable propagation shape of the debonding for several values of γ. Higher values of γ result in an energy-driven configuration where no further propagation of the debonding is possible after initiation (γ = 1.5). On the opposite, a large unstable debonding phase is encountered when γ decreases, leading to a larger arrest debonding angle (γ = 1). Although the initiation debonding angle differs from the location maximizing G/G c , the unstable debonding retrieves an identical shape regardless of the interface properties. Configuration with 0 deg. inter-center angle highlights a transition zone where neither stressbased nor energy-based debonding shapes describe the CC solution (see Figure 9c). By computing the favorable solution that minimizes the loading necessary for initiation based on either stress, energy or CC shape, several starting angles can be used to assess the unstable debonding phase. Figure 12 shows the three corresponding debonding shapes. As the required remote loadings for the three solutions have relative differences smaller than 1%, the ERR to critical ERR ratio surface presented is based on the CC solution for illustration purposes. However, the ERR calculation differs because the three solutions have different initiation remote loadings, so that the shape could be slightly different. All cases provide a similar debonding arrest angle with the same location regarding θ i even with a slightly different initiation remote loading. Moreover, all shapes retrieve the same trajectory that maximizes the ERR to critical ERR ratio. Regardless of the type of shape used (stress, energy, CC), a configuration driven by both stress and energy will propagate to approximately the same debonding arrest angle, as the remote loading required for both solutions are close. Some limitations exist for energy-driven configuration, where a stress-based debonding shape solution could lead to an unstable propagation phase while the energy-based solution leads to a stable propagation. However, the final solution can be obtained using the local ERR released for each configuration, that is computationally expensive to determine.

Nearby fiber influence on debonding initiation

Perturbation induced by the neighboring fiber influences the debonding location along the matrix-fiber interface. However, fields intensity are also influenced by the neighboring fiber, both in terms of stress and energy, depending on the geometrical parameters d and α. Stress fields at the fiber-matrix interface are therefore either intensified or damped and the same observation can be made concerning the elastic strain energy released by the debonding initiation. For similar interface properties, debonding initiation location (top or bottom) and debonding angle range vary as a function of both inter-center angle and distance.

Figure 13a shows the variation of the maximum normal to far field stress ratio evaluated at both top and bottom half circles of the reference fiber as a function of the neighboring fiber position (varying d and α). The bottom half circle is less affected than the top half circle. The latter is more affected by the neighboring fiber, highlighting a larger variation in the stress field. Also, with the distance between the top half circle and the neighboring fiber reduces, larger stress perturbations occur. The stress intensification reaches a maximum for 90 deg. configuration. For a configuration close to 45 deg., the top half circle stress decreases, resulting in a higher stress at the bottom half circle. Figure 13b shows the variation of the maximum IERR to critical ERR ratio which is influenced similarly to the stress fields. The latter is obtained by releasing debonding based on the tensile stress isocontours. Higher IERR to critical ERR ratio is observed for angles close to 90 degrees, where the stress intensity is maximum. However, for inter-center angle range from 20 to 60 degrees, IERR to critical ERR ratio appears to be higher at the bottom half circle. Based on the latest observations, it is more likely that debonding would occur at the top half circle for inter-center distances close to 90 degrees at a lower imposed loading. The smaller the distance between the centers, the smaller the load required for debonding to initiate. Conversely, for an inter-center angle close to 45 degrees, debonding could most likely occur at the bottom half circle. The perturbation of the stress and energy fields at the fiber-matrix interface therefore has a significant influence on:

• The debonding initiation location;

• The required loading for debonding to initiate.

In the following, the CC solution is compared to experimental observations for different geometrical configurations.

Comparison with experimental observations

The interface properties of a single fiber sample made of the same constituents were previously determined by inverse identification [START_REF] Girard | Numerical simulation of fiber-matrix debonding: Inverse identification of interface properties[END_REF]. The normal and shear strengths were selected within a range of acceptable strengths. Table 4 summarizes the identified interface fracture properties. An application of the CC for all configurations is carried out in the following. The debonding initiation stresses at the top and bottom fiber half circles obtained using the CC and debonding shapes based on stress isocontours are compared with experimental measurements. Furthermore, considering the relevance of the CC, the comparisons will allow to challenge or validate the accuracy of the interface properties determined previously [START_REF] Girard | Numerical simulation of fiber-matrix debonding: Inverse identification of interface properties[END_REF].

The debonding initiation sites (reference fiber top or bottom half circle) are compared with the experimental observations in are found to be favorable at the top half circle are denoted as "T" and at the bottom half circle as "B". Concerning the simulations, when the difference between the top and bottom half circles is less than 1% relative error, the location of the debonding is designated as "T or B" since the position of the debonding may take place at either half circle. The debonding locations show no particular trend as a function of the angle and distance inter-center. The experimental results show poor agreement with the CC solution, especially for an inter-center angle of 90 degrees. The difference observed may be associated with the idealized loading employed numerically whereas the experimental loading can differ from a simplified far-fields loading. Additionally, 2D DIC allows one single free edge to be observed at a time, while the initiation of debonding may occur on the other side due to geometric imperfections for instance. Figure 14 shows the load required to initiate debonding as a function of the inter-center angles. These values are evaluated experimentally and using the CC for three different inter-center distances. Uncertainties of the experimental results are obtained by tracking a sudden change of slope in the extension of the virtual gauges located at each fiber pole [START_REF] Girard | Experimental characterization of inplane debonding at fiber-matrix interface using single glass macro fiber samples[END_REF]. Two samples per configuration are studied, resulting in 8 different values. The results of each pole are averaged and the two extremes are then used for the estimation of the uncertainties.

For all angles and inter-center distances, the CC solutions are indicated by round and square markers for the top and bottom position, respectively. The CC solution that requires a load is the favorable configuration. The latter depends on the inter-center angle. For angles close to 0 degrees, both the top or bottom locations are favorable (difference smaller than 1%). For angles between 30 and 60 degrees, the favorable initiation site is located at the bottom half circle. Finally, for angles close to 90 degrees, the top half circle seems to be the most favorable initiation site. The higher the inter-center distance, the smaller the difference between the top and bottom solutions, due to lower field perturbation, see Figures 13 and14.

The overall comparison leads to a good agreement between the CC solution and the experimental observation. This observation attests to the accuracy of the interface properties determined previously. However, a larger difference is observed for an inter-center angle of 90 degrees for all three inter-center distances, possibly due to the aforementioned difficulties in the identification of the debonding locations. The experimental determination of the debonding angle at initiation after the potential unstable propagation is relatively difficult. The main difficulty lies in the identification of the locations of the two crack tips, possibly masked by the DIC speckles or considered as a process zone. Moreover, experimental capture of the unstable propagation phase remains impossible without the use of a rapid camera. In the sequel, experimental debonding angles are therefore compared to arrest angles obtained numerically. The CC solutions provide an almost identical debonding angle after initiation, close to 77 degrees, whatever the inter-center angle and distance. Figure 15a shows the variation of the debonding angles as a function of the inter-center angles. The experimental values correspond to the debonding angles averaged and their standard deviation observed at initiation for all inter-center distances with fixed inter-center angle. Only small variations are seen and the experimental angles are close to 62 degrees. The latter is close to the averaged CC solution, with 24% relative difference, as the CC solutions slightly overestimate the experimental observations. Figure 15b shows the location and size of the debonding angles predicted by the CC using the stress isocontours and the corresponding experimental results, for the smallest distance (d = 2.5r), extracted from a half circle of one of the two samples. The CC solution accurately describes the location of the debonding by providing an upper limit to the experimental angle as observed on the overall trend (Figure 15a). The CC solution and the experimental angle locations and sizes are in good agreement, endorsing the stress isocontour as an effective approximation to describe the debonding shape. Moreover, even in the presence of the neighboring fiber, the experimental debonding is approximately centered at 90 degrees. The same result is obtained with the CC solution which validates the previous numerical observation (see Section 3.5) concerning the debonding initiation angle obtained with a strongly influenced θ i .

Conclusion

Debonding initiation and propagation were assessed using both the CC and Griffith criterion on a two-fiber sample subjected to remote uniaxial tensile loading. The debonding shape was predicted on the basis of stress-and energy-based isocontours compared to the CC solution. Finally, the results of the CC solution were compared with the corresponding experimental results.

Using normal stress-based or equivalent stress-based debonding configurations lead to an almost identical initiation configuration. Normal stress-based debonding configurations can be efficiently extracted without introducing any fracture properties and therefore appear as a relevant compromise to describe the stress-based debonding configurations.

The debonding angle size and location predicted by the CC lies between possible debonding configurations based on energy or stress isocontours. Depending on the interface properties, e.g. the interface brittleness number γ, three different configurations can be encountered:

• Small γ: the stress-based debonding prediction provides a reasonable approximation for the possible initiation debonding shapes;

• Intermediate γ: the possible initiation debonding shapes lie in-between the stress-based and energy-based debonding predictions;

• Large γ: the energy-based debonding prediction provides a good approximation, as the solution tends to be energy driven.

Despite the difference observed between the three possibilities, similar initiation remote loadings are obtained using either the stress-based, the energy-based or the CC-based debonding configurations. The stress-based debonding configurations therefore appear to be a good compromise as it is numerically much more efficient since it can be extracted based on a single calculation without crack. However, it is not possible to accurately describe the unstable debonding phase based on stress isocontours because the stress-based and energy-based debonding configurations may differ widely for some configurations. Therefore, the final location of the angle described by the two shapes could be slightly different. However, basing the debonding configuration on stress isocontours and assessing the unstable debonding phase on energy isocontours leads to the same debonding arrest angle, but is computationally more expensive. derestimation of the critical shear ERR that strongly influences the debonding arrest angle. The CC results obtained for a 90 degree inter-center configuration differ from the experimental results, which could be explained by the idealized loading applied numerically or by an experimental intercenter angle slightly different from 90 degrees. It is also important to mention that the properties used were determined using a 2D simulation without considering the singularity acting at the free edge. Therefore, the 3D simulation could allow a better modeling of the problem, especially of the 90 geometrical configuration where the CC fails to accurately reproduce the experimental results.

Figure 2 :

 2 Figure 2: Equivalent stress and elastic strain energy variation as a function of the angular position along the fibermatrix interface for (a) single-fiber configuration (b) two-fiber configuration for a constant µ = 1.5.

  angle (α) [deg.]

Figure 3 :

 3 Figure 3: Angle θ i for which the maximum equivalent stress is attained at the top and bottom fiber half circles (T and B respectively) as a function of the fiber inter-center distance and angle for µ = ∞.

µ

  Angular position (θf ) [deg.]

Figure 5 :

 5 Figure 5: Equivalent stress to tensile strength ratio variation as a function of the angular position along the reference fiber top interface for various µ values.

Figure 6 :

 6 Figure 6: (a) Equivalent stress to tensile strength variation as a function of angle increments (b) Stress isocontours for several values of µ where tension zone (σ nn > 0) and compression zone (σ nn < 0) are highlighted. The dashed lightgray line corresponds to symmetric configuration.

  summarizes the IERR obtained for a randomly selected isovalue for µ = 0.5 where α = 45 deg. and d = 2.5r. Both maxima θ f 40 deg. maximum θ f 132 deg. maximum G inc [N/mm] 3.97 • 10 -3 4.13 • 10 -3 3.7 • 10 -3

Figure 7 :

 7 Figure 7: (a) IERR to critical ERR variation. (b) Debonding shape for several values of λ and the dashed lightgray line corresponds to symmetric opening around θ i .

σα

  eq /σ c G inc /G c α = 45 deg., d = 2.5r = 45 deg. d = 2.5r

Figure 8 :

 8 Figure 8: (a) Required loading to satisfy stress or energy conditions where solid lines indicate the configurations corresponding to either stress-based or energy-based possible debonding shapes and circles indicates the corresponding initiation configurations minimizing the imposed loading. (b) Isocontours of the loading required to fulfill the CC emphasizing the minimum indicated in Fig. 8a and minima obtained for various other properties range, the smaller the gray intensity, the higher γ is. Visualization of θ i 100 deg. and associated debonding angle at initiation, before any potential unstable growth, taken from Fig. 8b for several values of γ, (c) γ = 0.66 with θ d = 33 deg. (yellow circle), (d) γ = 1.16 with θ d = 85 deg. (green marker) and (e) γ = 1.45 with θ d = 102 deg. (orange circle).

  shape λ = 0.2 Stress-based shape µ = 1.5

Figure 9 :

 9 Figure 9: Comparison between debonding shapes for configuration where d = 2.5r and (a) α = 90 deg. (b) α = 60 deg. and (c) α = 0 degrees. The markers denotes the CC solutions for various γ values and dashed lightgray lines correspond to symmetric openings on both sides from θ i whose values are given in Table2.

Figure 10 :

 10 Figure 10: (a) Relative difference between CC and stress-based debonding shape solution as a function of γ. (b) Schematic comparison of the debonding location and angle at initiation predicted based on both stress (θ d 100 deg.) and CC (θ d 97 deg.) shapes with θ i 60 degrees.

Figure 11 :

 11 Figure 11: (a) ERR to critical ERR variation. (b) Isocontours of the critical ERR (G c ) and unstable debonding shapes obtained for several values of γ.

  deg. d = 2.5r γ = 1.16 CC Stress Energy θ d [deg.] CC Stress Energy Initiation 67.6 75.

Figure 12 :

 12 Figure 12: Unstable debonding shape obtained from the stress-based, energy-based and CC-based initiation debonding angle.Table shows the associated debonding initiation angle and arrest angle obtained after unstable propagation phase.

Figure 13 :

 13 Figure 13: (a) Maximum normal to far-field stress ratio variation and (b) maximum IERR to average critical ERR ratio variation as a function of the inter-center angle and distance, for an arbitrary remote loading.

Table 5 :

 5 Debonding location (T for top and B for bottom) for several inter-center distances (columns) and angles (lines). Within a single cell, the upper right term corresponds to the experimental result and the lower left value to the computed CC one.

Figure 14 :

 14 Figure 14: Loading required for debonding initiation as a function of the inter-center angle for (a) d = 2.5r (b) d = 4r (c) d = 5r. The scatter bars associated with the experimental data points represent variations between several independent tests.

  [-] 1.02 1.00 1.03 0.96 1.02

	α [deg.]	0	30	45	60	90
	θ i σeq/σc	[deg.]	87	93	98	63	90
	θ -∆W i	[deg.]	89	95	102	60	90
	θ -∆W i	σeq/σc i /θ
	Table

Table 3 :

 3 IERR for a given initiation configuration considering single initiation site separately or both sites.

  Table shows the associated debonding initiation angle and arrest angle obtained after unstable propagation phase.

Table 5 .

 5 For various geometrical configurations, debondings thatG IC [N/mm] G IIC [N/mm]

	λ	σ c [MPa] τ c [MPa]

Table 4 :

 4 Fiber-matrix interface properties from[START_REF] Girard | Numerical simulation of fiber-matrix debonding: Inverse identification of interface properties[END_REF].

	d	α	0 deg.	30 deg.	45 deg.	60 deg.	90 deg.

Blind application of the CC and comparison with experimental results showed good overall agreement both in terms of required remote load and debonding angle after propagation using previously identified interfacial shear and tensile strengths and critical opening ERR. The debonding angle after initiation is slightly overestimated by the CC, which could be explained by an un-

Recalling that G inc , G c , G and G c are depending on the debonding surface S. When, a debonding initiates for a debonding surface maximizing G inc /G c :

It derives from the right part:

The IERR can be expressed as:

Integrating A.4 into A.3:

The local critical ERR can be expressed as:

Replacing into A.8:

Remembering that:

So it can be stated that:

Finally, the following equality is obtained by replacing A.13 into A.11:

As demonstrated by Mantič [START_REF] Mantič | Interface crack onset at a circular cylindrical inclusion under a remote transverse tension. Application of a coupled stress and energy criterion[END_REF], dG dS < dGc dS for S = S min , so the following equality is obtained: