Numerical analysis of the crack paths produced by fibre–matrix interface failure in cross-ply LFRP laminates

Composite Structures 284 (2022) 115222 Available online 8 January 2022 0263-8223/© 2022 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). Numerical analysis of the crack paths produced by fibre–matrix interface failure in cross-ply LFRP laminates Laura Moreno a , Luis T´ avara b , * , Elena Correa b a Departamento de aeronaves y vehículos espaciales, Escuela T´ ecnica Superior de Ingeniería Aeron´ autica y del Espacio, Universidad Polit´ ecnica de Madrid, Pza. del Cardenal Cisneros 3, 28040 Madrid, Spain b Grupo de Elasticidad y Resistencia de Materiales, Escuela T´ ecnica Superior de Ingeniería, Universidad de Sevilla, Camino de los Descubrimientos s/n, 41092 Sevilla, Spain ARTICLE INFO Keywords: Composite materials Inter-fibre failure Laminate Interface crack Debonding LEBIM ABSTRACT The so-called matrix failure, also named inter-fibre failure, is one of the main failure mechanisms in long fibre reinforced polymers (LFRPs). This mechanism is characterized by the debonding of some fibres when transverse loads are diving the failure process. In the present study, the prediction of cracks initiation and growth within fibre–matrix interfaces of the 90◦layer in [0,90 n ] s laminates is aimed. To this end, the behaviour of the fibre–matrix interface is characterized using a Linear Elastic Brittle Interface Model (LEBIM) that has been implemented in the commercial FEA code Abaqus ® and applying a solution algorithm built in Python called the Sequentially Linear Analysis (SLA). The results point out the capability of the debonding process to represent the generation of transverse cracks, showing how the onset and growth of interface cracks within a laminate are influenced by different factors such as the local distribution of fibres, the value of the fibre–matrix interfacial fracture properties and the thickness of the 90◦layer. 1. Introduction The so-called inter-fibre/matrix failure, that appears in laminates subjected to transverse loads, has been the object of previous studies by the authors [1–4], and other researchers [5–8]. According to some authors [9,10], the initiation of this failure could be linked to matrix cavitation or the appearance of a debond at the fibre–matrix interface. In these previous studies [1,4], the authors made the assumption that the failure starts with the appearance of small debonds at the interface and identified the following stages of the damage mechanism at micromechanical level in an isolated fibre, Fig. 1: Stage 1: Crack nucleation at the interface in a position aligned with the applied load. This result is based on the assumption that the interfacial radial stress between the fibre and matrix controls the first debond appearance. Stage 2: Unstable and symmetrical growth of the crack up to a semidebonding angle of 60-70◦. Stage 3: Kinking: The crack at the interface suddenly changes its direction of growth and advances towards the matrix in a direction perpendicular to the applied load The use of multi-fibre models has been performed by several authors to study the behaviour of composite laminates [11–15]. The present work presents a multi-fibre model and concentrates on the first and second stages of the inter-fibre failure. In particular, the objective of this paper is to present a tool that is able to predict the initiation and growth of fibre–matrix debonds of the 90◦layer of [0,90 n ] s laminates by modelling the whole laminate. This study puts special emphasis on the analysis of the crack pattern caused by several fibre–matrix debonds that may lead to a crack path (transverse crack in the 90◦layer) or to the appeareance of isolated debonds. Although the transverse crack in the 90◦layer is caused by the coalescense of the debonds by a kinking of the interface cracks into the matrix, in the present study only the debonds are allowed (no matrix failure is included). This assumption is coherent with the microstructure that can be frequently observed in current laminates where the V f achieved leads to touching fibres, and thus no kinking is present in all fibres, see Fig. 2. In addition, previous numerical studies that analyse the generation of transverse cracks considering both debonding and matrix cracks show that, for high V f values, more than 90% of the damage [16] and residual damage [17] are associated with fibre–matrix debonding. To this end, the fibre–matrix interfaces are modelled using the Linear * Corresponding author. E-mail address: [email protected] (L. T´ avara). Contents lists available at ScienceDirect Composite Structures journal homepage: www.elsevier.com/locate/compstruct https://doi.org/10.1016/j.compstruct.2022.115222 Received 29 October 2021; Received in revised form 13 December 2021; Accepted 4 January 2022 Composite Structures 284 (2022) 115222 2 Elastic Brittle Interface Model (LEBIM). LEBIM can be considered as a limit case of the well-known cohesive zone model (CZM) formulation [18]. LEBIM has proven to be a robust solution when compared to traditional CZMs and additionally requires less parameters than CZMs. Moreover, LEBIM is capable of including several failure criteria in a straightforward manner (i.e. Hutchinson and Suo [19], Benzeggagh and Kennane [20], among others, see [21] for a review). Furthermore, several studies have proven that the shape of the CZM (including LEBIM as limit case) does not affect the results when the overall stiffness of the problem under study is driven by the adherents (i.e. a change in the interface/adhesive stiffness does not affect the global load–displacement curve) [18,22,23]. In addition, a linear solution algorithm named Sequentially Linear Analysis (SLA) has also been applied. This algorithm is able to solve problems in the presence of large instabilities. One of its main advantages is that SLA is able to get information along the unstable zone(s) of the problem under study (i.e. local damage, crack grow at several crack tips, etc.). This kind of information is relevant when the crack path is analysed [18,24]. Thus the analysis of the crack patterns in the present research includes: (i) The study of the stress state in zones close to the crack pattern; the aim is to check that a group of aligned debonds behaves in a similar way than a transverse crack (ii) a parametric study of the interface properties and their effect on the crack patterns (iii) a study of the effect of the thickness of the 90◦layer on the crack pattern, in addition to the effect that could cause the proximity of the fibres and their disposition within the matrix. The structure of the paper is detailed next. In Section 2, the numerical tools used are presented. In Section 3, the model under study is described and, in Section 4, the results obtained are shown. Finally, in Section 5 the conclusions are summarized. 2. Numerical tools The main features of the numerical tools used and implemented in Abaqus ® are presented next, for further details see [18,25,26]. 2.1. LEBIM +Abaqus In the present study, the fibre–matrix interfaces are modelled using the Linear Elastic Brittle Interface Model (LEBIM), representing a continuous spring distribution. LEBIM was used satisfactorily to model the behaviour of interfaces within composites [25,26,27]. This kind of behaviour is included in Abaqus ® [28] by means of a user subroutine UMAT written in Fortran. The UMAT is designed to work at integration point level of the elements included along the interfaces. Then, two kinds of states could be found at each integration point: Undamaged interface point { σ (x) = kn ε n(x) τ (x) = kt ε t(x)(1) Damaged interface point { σ (x) = kn〈 ε n(x)〉− τ (x) = 0(2) where σ (x) and τ (x) are the normal and shear components of the stress vector; k n and k t denote the normal and tangential stiffness of the springs distribution, and ε n and ε t the normal and shear components of the strain vector at an interface point. Also, 〈〉±represents the positive or negative part of a function. Additionally, an energetic criterion is adopted so that the failure at a point x of the interface occurs when G(x) ≥G c ( ψ (x)), G(x) being the Energy Release Rate (ERR) [19,21] and G c ( ψ (x)) the fracture toughness Fig. 1. Stages of initiation and growth of the crack in the interface under uniaxial tension at the micromechanical level. Fig. 2. Transverse crack in a carbon-epoxy specimen tested under transverse tension. L. Moreno et al. Composite Structures 284 (2022) 115222 3 that depends on the fracture mode mixity ψ (x) and follows the proposal of Hutchinson and Suo [19]: Gc( ψ ) = Gc( ψ ) GIc =1+tan2[(1−λ) ψ ](3) Where G Ic is the Mode I ERR critical value for the interface and λ is the sensitivity mode parameter that characterizes interfaces with moderately strong fracture mode dependence [19]. The typical range is 0.2 ≤λ ≤0.3. With reference to G(x) calculation, it follows (see [21] and references therein): G(x) = GI(x) + GII (x) = kn〈δn(x)〉2 + 2+ktδ2 t(x) 2(4) 2.2. SLA +Abaqus ® The interface crack growth at a particular fibre–matrix pair has proven to be a very unstable process [21]. When this event appears simultaneously in more than one fibre this whole process can hardly be detected with usual solution algorithms, as those built-in on commercial FEA codes as Abaqus ® [28]. Several authors have used the Sequentially Linear Analysis (SLA) to obtain a solution in problems with the presence of large instabilities [18,21,24]. In the present research SLA is implemented through a Python script linked to Abaqus ® [18,25]. In the present SLA implementation, at the beginning of the problem solution, every interface element (and their corresponding integration points) is pristine (undamaged). Then, a sequence of solution steps is performed. Each step includes: (i) a linear elastic solution of the system, (ii) a post-process of the results which allows to determine those integration points that are most likely to fail (based on the used failure criterion). (iii) The creation or modification of a text file where the chosen (damaged) integration points and their corresponding elements are stored. 3. Details of the models The analysed models are developed in Abaqus ® [28] and represent [0, 90 n ] s laminates where the three main layers are clearly distinguished, see Fig. 3. The 90◦layer is constituted by an epoxy matrix where a random distribution of carbon fibres is embedded. The upper and lower 0◦layers are modelled as a homogeneous orthotropic material. The dimensions of the models (B, H and e) are similar to those selected in Arteiro et al [13], although slight modifications are included. In particular, B =200 μ m, H =75 μ m and three different values of e have been considered (20, 40 and 60 μ m). On the other hand, for every model in [13], obtained varying e, two models are analysed in the present manuscript, i.e. B in the present study is half of the B value used in [13]. Moreover, symmetry boundary conditions are used at one of the sides, see Fig. 3. This option was chosen for two reasons: (i) to study the influence of a different fibre distribution and (ii) to reduce the computational cost of the models. Analysed models are denoted by L (left part) and R (right part), shown in Figs. 3, 4 and 5, additionally, plane strain hypothesis (2D) is assumed. Every model includes linear - quadrilateral elements with full integration, i.e. including 4 integration points (C3P4 in Abaqus ®). Fibres are assumed to be circular and equal and the interface is modelled as an annulus of small thickness, t, see white zone in Fig. 6(b). The mesh is designed to be the same at each fibre and interface zone. Conforming Fig. 3. Boundary conditions applied to the models with e =20 μ m (a) L: left part of the laminate in [13] , L. (b) R: right part of the laminate in [13]. Fig. 4. Distribution of fibres for e =40 μ m (a) L: left part of the laminate in [13], L. (b) R: right part of the laminate in [13]. L. Moreno et al. Composite Structures 284 (2022) 115222 4 meshes are used in the fibre-interface and interface-matrix zones. A detail of the mesh used is shown in Fig. 6 for the left 20 μ m case. Regarding the epoxy matrix (red zone in Figs. 3-6), presented in the 90◦layer, it is modelled using a linear elastic isotropic constitutive law. The Young’s modulus, E m , considered is 2.79GPa. and the Poissons ratio, υ m , is 0.33. The number of elements at the matrix zone and the fibre volume fraction for each model are shown in Table 1 On the other hand, carbon fibres (green zones in Figs. 3-6) are also modelled following a linear elastic isotropic constitutive model, using a radius equal to 2.60 µm and E f =15.032GPa and υ f =0.25273. The number of elements at each fibre is 1356, details of the number of fibres included at each model are compiled in Table 1. The behaviour of the interface between the fibre and the matrix is modelled using LEBIM (described in Section 2.1) implemented in Abaqus ® through a UMAT subroutine. 92 equally distributed elements are used along the interface (element size is ~ 3.91◦) with a constant thickness, h, equal to 0.01 μ m. Baseline LEBIM parameters were taken from literature [21,25]: σ c =90 MPa, G Ic =2 J/m 2 , k t /k n =0.25 and λ = 0.25. Finally, the upper and lower 0◦layers (blue zones in Fig. 5) follow a homogeneous linear elastic orthotropic law. 1400 elements are used at each layer. The elastic properties and number of elements used are shown in Table 2. 4. Results The usual initial failure mechanism in cross-ply laminates (0◦/90◦) is the appearance of a transverse crack which is produced by the nucleation and coalescence of aligned debonds. Although the analysed models allow the presence of multiple debonds at any position, their coalescence is not allowed, as it was already explained and justified in the Introduction. Moreover, the models are able to capture, at the local level, that debonding occurring at a single fibre–matrix interface follows the aforementioned stages I and II (initiation and unstable growth) [18]. As an example of the model capabilities, in Fig. 7, the applied strain versus the number of damaged integration points is depicted for the model left 20 μ m. Applied strain is calculated as the imposed (uniformly Fig. 5. Distribution of fibres for e =60 μ m (a) L: left part of the laminate in [13], L. (b) R: right part of the laminate in [13]. Fig. 6. Detail of the mesh for the left 20 μ m case. Table 1 Number of elements in the matrix zone, fibre volume fraction (V f ) and number of fibres included at each model. e (µm) N◦elements in the matrix zone V f (%) N◦of fibres Left 20 19241 56.4 52 40 36722 53.5 101 60 57843 56.7 159 Right 20 18686 51.3 48 40 40392 53.5 101 60 72528 57.1 160 Table 2 Properties of the orthotropic layer. E 1 (GPa) E 2 =E 3 (GPa) υ 12 = υ 23 υ 13 G 12 (GPa) G 13 =G 23 (GPa) 171.42 9.08 0.32 0.4 5.29 3.92 L. Moreno et al. Composite Structures 284 (2022) 115222 5 Fig. 7. Strain, ε (%) vs. number of damaged integration points for the 20 μ m left case. Fig. 8. Damage variable (debond) produced at several fibres for the 20 μ m left case. Deformed shape for ε =2.213% and scale factor 3. Steps: (a) 50, (b) 100, (c) 150, (d) 200, (e) 250 and (f) 300. L. Moreno et al. Composite Structures 284 (2022) 115222 6 distributed) displacement on the lateral side divided by the initial length B/2. It is interesting to recall that when SLA is applied, each step provokes the failure of 2 integration points within an interface element (i.e. half element is “broken” ~1.96◦), thus the number of damaged integration points can be directly related to the total debond length (sum of debonds at several fibres). In the plot each valley is associated with a debond within a single fibre, then it can be easily detected from the plot the number of debonded fibres. As it can be noticed from Fig. 7, SLA is able to capture the whole instable process. It is interesting to notice that from a physical point of view, the first critical strain (provoking the first debond) also causes debonds in other fibres. In the specific case shown in Fig. 7, the first critical strain, 2.213% (depicted with a discontinuous blue line), causes the appearance of debonds in 5 fibres (which complete determines the Fig. 9. Damage variable (debond) produced at several fibres for the 20 μ m left case. Deformed shape for ε =2.361% and scale factor 3. Steps: (a) 350, (b) 400, (c) 450, (d) 500 and (e) 550. Fig. 10. Strain, ε (%) vs. number of damaged integration points for the 20 μ m right case. L. Moreno et al. Composite Structures 284 (2022) 115222 7 path of a complete transverse crack). After that, a new maximum strain, 2.361% (depicted with a discontinuous blue line), provokes the debonding in other set of 4 fibres (forming another transverse crack path). In Fig. 8, the damage variable is shown for steps 50, 100, 150, 200, 250 and 300, the deformed shape is plotted for a 2.213% strain and a scale factor of 3. While in Fig. 9, the damage variable is shown for steps 350, 400, 450, 500 and 550, the deformed shape is plotted for a 2.361% strain and a scale factor of 3. It is interesting to recall that the numerical tool predicts the appearance of aligned debonds (leading to a crack path), in this particular case and as mentioned, the second path is formed once the first one is completed. Once the numerical tool capabilities are shown, 3 analyses are carried out in the following: (i) The study of the stress state in zones close to a crack pattern (aligned debonds) with the aim to check that a group of aligned debonds behaves in a similar way than a transverse crack. (ii) A parametric study of the interface properties and their effect on the crack patterns. (iii) A study of the effect of the thickness of the 90◦layer on the crack pattern, in addition to the effect that could cause the proximity of the fibres and their disposition within the matrix. 4.1. Stress state in the neighbourhood of aligned debonds The aim of this section is to check if the presence of a set of aligned debonds leads to a similar stress state as the one produced by a transverse crack. i.e. stress relaxation close to a transverse crack. Although the analysed model corresponds to the 20 μ m right case, similar behaviours in terms of stresses are also observed for the rest of the cases. In Fig. 10. the applied strain versus the number of damaged integration points is depicted for the 20 μ m right model. It is interesting to recall that in this case, isolated debonds appear in contrast to what was presented in the previous section for the 20 μ m left model. First, the first critical strain, 2%, provokes just 1 debond, then, a strain equal to 2.21% causes the appearance of 2 debonds. Next, a set of 6 aligned debonds appear for a 2.27% strain value. Finally, a strain equal to 2.347% provokes 3 additional debonds. σ I has been selected here as the representative of the stress state aiming to evaluate the possible stress relaxation due to the relative position of the debonds versus the external load and under the assumption that the maximum principal stress criterion could be in this case quite representative of the process. Therefore, in Fig. 11, σ I distribution in the Fig. 11. σ I distribution in MPa for the 20 μ m right case. Deformed shape with scale factor 3. (a) Step 1, ε =2%, (b) Step 57, ε =2%, (c) Step 171, ε =2.21%, (d) Step 500, ε =2.27% and (e) Step 711, ε =2.347%. Fig. 12. Damage variable (debonds) produced at several fibres for the 20 μ m right case. Deformed shape with scale factor 3, step 711, ε =2.347%. L. Moreno et al. Composite Structures 284 (2022) 115222 8 90◦ply zone is shown for several steps. It is interesting to notice that, in Fig. 11(a) the stress distribution is quite uniform. A similar behaviour is observed in Fig. 11(b) and (c), except in the zone close to where isolated debonds appear. On the other hand, in Fig. 11 (d) and (e), where aligned debonds appear, the stress distribution in the zones close to the aligned debonds is clearly lower than in the zones far from them. This fact may indicate that aligned debonds have a similar behaviour as an “actual” transverse crack. Fig. 12 shows the damage variable for step 711 with an applied strain equal to 2.347%, where the numbers inside the fibres indicate the order of the debonding mechanism occurrence. The colour of the number indicates a set of debonds produced by the same applied strain (see Fig. 10). In order to see how σ I value changes while the debonding mechanism occur, in Fig. 13 the evolution of σ I at every load step for 3 elements is plotted. Specifically, the 3 chosen elements are the ones indicated in Fig. 12 and denoted as A, B and C. As expected, the stress state on Element A, clearly decreases after the appearance of the sixth debond, while the stress state on element C only decreases after the twelfth debond. On the other hand, the stress state on Element B (quite far from the damaged zones) remains at a quasi-constant level as it is not affected by the appearance of any debond. 4.2. Effect of the interface properties on the crack pattern The interface crack (debond) behaviour could be directly affected by the chosen interface properties [21]. Specifically in this section fracture mode I parameters needed in LEBIM are analysed, i.e. the interface strength, σ c , and the fracture toughness, G Ic . The presented results correspond with the 20 μ m right case (R20), although similar conclusions can be obtained with the rest of the cases. 4.2.1. Interface strength Five different values of the interface strength, σ c , are used for the present parametric analysis: 30, 60, 90, 120 and 150 MPa. All the chosen values are within the range of values found in literature. The rest of the LEBIM parameters remains constant in this study (G Ic =2 J/m 2 , k t /k n = 0.25 and λ =0.25). Table 3 presents, for each σ c case, the values of the strain applied along the fracture process until two crack paths are formed. The number of debonds generated at each strain increment are also compiled. Notice that the number of steps required for the formation of the two crack paths (also detailed in Table 3) do not necessarily coincide in all σ c cases and neither the length of the debonds. The first conclusion derived from the results presented in Table 3 is that as σ c increases the value of the strain required for the appearance of the first debond decreases. This result is coherent with the fact G Ic is forced to be constant and thus the effect in the LEBIM law of increasing σ c (maintaining the same value of the area under the curve, i.e. G Ic ) is to increase the stiffness of the interface, k n . As a consequence, δ nc and thus the strains, lowers. In fact, in the other extreme, it is remarkable that the first strain value for the 30 MPa case, i.e. 2.68%, represent the highest final strain value of the rest of σ c cases, which might indicate the singular character of this first case associated with a very small value of the interface stiffness. Fig. 14 presents the debonding situation in the last step specified in Table 3 for each σ c case. The debonding sequence is also detailed through the numbers in the fibres and the colours employed, as it was already done in Fig. 12. From the observation of Fig. 14, it can be deduced that for medium–high values of σ c (Figs. (c), (d) and (e)) the results are similar: the first debond appears isolated and the first crack path is located in the left side of the model (the central part of the complete model) and the second crack path in the right side. Therefore, the selection of σ c in the 90–150 MPa range would not alter the conclusions obtained. For low values of σ c the situation changes: the order of appearance of the transverse paths is the opposite to the one exposed for the medium–high cases and even changes its location for the 30 MPa extreme Fig. 13. σ I evolution at elements A, B and C indicated in Fig. 12. Table 3 ε (%) and number of debonds for σ c =30, 60, 90, 120 and 150 MPa cases associated with the R20 model. σ c (MPa) [step] 30 [450] 60 [500] 90 [700] 120 [550] 150 [700] ε (%) [number of debonds] 2.68 [1] 2.17 [1] 2 [1] 1.88 [1] 1.87 [1] 2.73 [1] 2.19 [1] 2.21 [2] 2.14 [4] 2.11 [3] 2.82 [4] 2.25 [3] 2.27 [6] 2.28 [4] 2.19 [1] 2.92 [2] 2.26 [5] 2.35 [3] 2.23 [1] 2.94 [1] 2.26 [1] 2.96 [3] 2.37 [4] L. Moreno et al. Composite Structures 284 (2022) 115222 9 case. As mentioned before this case shows a singular behaviour that is reinforced by the results presented in Fig. 14(a). Specifically, the transverse paths are closer than in the rest of the cases and a large interface deformation is also predicted by the model. This deformation is represented by the circumferential blue marks that appear around many fibres. None of these results seems to reproduce a very actual situation, from an experimental point of view, which might discard 30 MPa as a reliable value for σ c . 4.2.2. Fracture toughness Three different values of the fracture toughness, G Ic , are used for the present parametric analysis: 2, 5, and 10 J/m 2 . All the chosen values are within the range of values found in literature. The other LEBIM parameters remain constant ( σ c =90 MPa, k t /k n =0.25 and λ =0.25). Table 4 presents, for each G Ic case, the values of the strain applied along the fracture process until two crack paths are formed. The number of debonds generated at each strain increment as well as the number of steps required for the formation of the two crack paths are also compiled. It is clear from the results that as G Ic increases the strain required to initiate the fracture process also does. This result, of direct physical sense, is coherent with the fact that σ c is forced to be constant and thus the effect in the LEBIM law of increasing G Ic is to reduce the stiffness of the interface, k n . As a consequence, δ nc and thus the strain, increases. This fact is clearly reflected in the 10 J/m 2 extreme case where the first strain value, i.e. 5.37%, is higher than all final strain values of the rest of G Ic cases, which might be an indicative of the singular character of this situation associated with a very small value of the interface stiffness. The 5 J/m 2 case also shows faraway strain values with reference to the 2 J/ m 2 case, therefore pointing out that those G Ic could lead to very different situations. Fig. 15 presents the debonding situation in the last step specified in Table 4 for each G Ic case. The debonding sequence is also detailed through the numbers in the fibres and the colours employed, as it was already done in Figs. 12 and 14. From the observation of Fig. 15, it can be deduced that starting from the fracture toughness reference value 2 J/m 2 multiplying it by a 2.5 factor changes the order or appearance of the transverse paths and avoid the presence of isolated debonds. The 10 J/m 2 extreme value maintains the crack paths location of the 5 J/m 2 case. The singular character of the 10 J/m 2 case (and even 5 J/m 2 case) is again pointed out by Fig. 15(c) (and (b)) where large values of the strain at the fibre–matrix interface of many fibres is predicted by the model (represented by circumferential blue marks) which could lead to the dismissing of this value as an accurate representation of the experimental reality. 4.3. Effect of the 90◦ply thickness on the crack pattern The six different configurations (two options for each thickness, 20, 40 and 60 μ m) described in section 3 are used herein to see the influence Fig. 14. Debonding sequence for the 20 μ m right case. Deformed shape with scale factor 3. (a) σ c =30 MPa, ε =2.96% (b) σ c =60 MPa, ε =2.26%, (c) σ c =90 MPa, ε =2.347%, (d) σ c =120 MPa, ε =2.28% and (e) σ c =150 MPa, ε =2.74%. Table 4 ε (%) and number of debonds for G Ic =2, 5 and 10 J/m 2 cases associated with the R20 model. G Ic (J/m 2 ) [step] 2 [700] 5 [497] 10 [411] ε (%) [number of debonds] 2.00 [1] 3.46 [1] 5.37 [1] 2.21 [2] 3.49 [1] 5.40 [1] 2.27 [6] 3.57 [4] 5.53 [5] 2.35 [4] 3.59 [4] 5.55 [1] 5.72 [2] L. Moreno et al.