Exploring the influence of micro-structure on the mechanical properties and crack bridging mechanisms of tufts.

A constitutive model for tufts bridging a mode I delamination is presented. The tuft is modelled as a rod, laterally supported by an elastic medium and clamped at both ends. A fracture mechanics approach is introduced to describe the progressive debonding of the tuft from the embedding laminate. The debond- ing model requires the identiﬁcation of stiffness, strength and toughness properties, which depend both on the laminate/tuft architecture and the constituent materials. Such identiﬁcation is carried out via experimental data obtained from tensile tests on single tufts inserted in a pre-delaminated non-crimp fabric composite. The experimental results are complemented by micro-scale ﬁnite element analysis. The mode I bridging law obtained from the constitutive model is implemented into a meso-scale cohesive zone formulation. This formulation is applied to predict the response to delamination of tufted Double Cantilever Beam (DCB) coupons. The cohesive zone approach is validated by means of experimental data from DCB tests. It is shown that the proposed micro- to meso-scale modelling approach yields results in good agreement with the experiments.


Introduction
Through-the-thickness reinforcement (TTR) is applied to 2dimensional composites in order to control and suppress delamination. Most common TTR methods include Z-pinning [1], stitching [2] and tufting [3]. Tufting is the most recent among them and is performed by inserting carbon, glass or aramid threads through the thickness of a dry preform by means of a single needle. Neighbouring tufts are interconnected to each other by a seam on one side of the preform and form thread loops on the other. Once resin infused, tufts become integral parts of the preform architecture, making it locally 3-dimensional. Despite the proved potential of tufts to counteract the propagation of delamination in composite parts [4], a complete study of their crack bridging behaviour is not available in the open literature.
The aim of this paper is to identify and describe the influence of micro-structure on the mode I crack bridging response of tufts and use the observations made at the micro-scale as the basis for the development of a multi-scale modelling framework for tufted composites.
An analytical micro-mechanical model is proposed to simulate the mechanical response of tufts embedded in mode I delaminating composites. The governing equations and assumptions of the model are supported by experimental results obtained for singletuft coupons, complemented with observations of the tuft architecture, morphology and failure mode. The suitability of this model for the prediction of the mechanical behaviour of bridged interfaces has been assessed via its implementation into the finite element model of a DCB coupon. A cohesive zone approach [5] has been adopted for this purpose and experimental data have been used to validate the overall multi-scale modelling strategy presented.

Single-tuft tests
A set of pre-delaminated single-tuft coupons has been tested under mode I conditions in order to derive the bridging law of the tuft, i.e. the relation between the relative displacement of the surfaces of a bridged crack and the force exerted by the tuft to counteract it [6,7]. The specimens were made of four layers of biaxial carbon Non-Crimp Fabric (NCF) with an areal weight of 1010 g/m 2 , stacked in a symmetric [(0/90) s ] 2 layup. The stack was separated at the mid-plane by a thin release film. Each 0°/90°layer contained equally arranged 24k HTS carbon fibre tows from Tenax, held together by non-structural stitching. Each coupon was tufted with commercially available 2k HTA40 carbon fibre sewing thread, having a dry cross-section area of 0.077 mm 2 . Tufts were inserted orthogonally to the release film. After insertion, each tuft featured a free loop end 3-5 mm long. The tufted preform was injected with aerospace grade epoxy resin (MVR444, Advanced Composites Group) using a Vacuum Assisted Resin Transfer Moulding (VARTM) process. Injection was carried out at 70°C and 1 bar pressure, followed by cure at 160°C for 90 min at 4 bar. The cured panel was subjected to post-cure in the oven at 180°C for 120 min. The final thickness was 4 (AE0.01) mm, with resulting global fibre volume fraction of 56.5%. The single-tuft coupons, with dimensions 20 mm Â 20 mm Â 4 mm, were tested in out-of-plane tension, under displacement control, at a cross-head speed of 0.25 mm/ min. A Digital Image Correlation (DIC) system was used to monitor the relative opening of the testing fixtures. Testing conditions are illustrated in Fig. 1.

Tuft morphology
Micrographic analysis was carried out to assess the post-cure tuft morphology, as shown in Fig. 2a. The insertion of tufts in a preform causes a local disruption of the in-plane fibre architecture, resulting in the formation of resin-rich regions around the through-thickness reinforcement. This is consistent with what has been reported for other TTR types in the open literature [7,8]. The cross section of the tuft is modelled by both the fabric architecture and the preform layup. Micrographs have revealed resin-rich regions characterised by maximum and minimum diameters of 5.6 mm (Coefficient Of Variation (COV) = 14.3%) and 0.55 mm (COV = 14.5%), respectively. The average impregnated cross-sectional area of the tuft, measured at 25% and 75% thickness of the samples, was 0.27 mm 2 (COV = 11%). Sectioning of the specimens has shown further that tufts are characterised by curved profiles and a random arrangement of their constituting thread segments, as in Fig. 2b. Such complex inherent features render a topological definition of tufts very difficult, and help explain the large experimental scatter in the derived bridging laws, as in Fig. 3a. Fig. 1. Mode I test on single-tuft specimen. The test is carried out in displacement control with an Instron 5500R and a 5 kN load cell. The arrows at the fours corners of the T-tabs identify the monitored displacements. Two cameras, one on the front and one at the back of the specimen, have been used at this purpose. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Mechanical response of tufts
With reference to Fig. 3a, four stages can be identified in the crack bridging response of tufts. Initially the bridging load increases, approximately in a linear fashion, up to a critical load in the range 50-100 N. Beyond this load, the slope of the loaddisplacement curves decreases significantly, characterising the second stage of the bridging response. Micrographs of partially loaded specimens reveal debonding cracks developing along the tufts at this stage, as shown in Fig. 3b. After reaching a peak value, P f , the load suddenly drops. At this stage, tufts break due to tensile fibre failure, which in the majority of cases occurs in proximity of the delamination plane. However, if fibre rupture occurs at a point within the composite, frictional pull-out takes place and, after a first large drop, the load decreases gradually to zero. Due to the tufts curvature, fibres are not loaded uniformly and some fail before others, as suggested by the presence of multiple load drops during the debonding stage. This process reflects in the scatter that can be appreciated in Fig. 3a, affecting the load to failure of tufts.
The mechanical energy involved in the process of mode I delamination bridging operated by tufts, and denoted by W I , is therefore determined by the work done to deform, disbond and finally pullout the tufts. For the tested specimens, the average value of this energy is 35.6 mJ.

A mode I constitutive model for tufts
Micro-mechanical models of through-thickness reinforcements include those of Jain et al. [9,10] and Cox et al. [11,12], which apply both to stitches and to short rods, and those of Allegri et al. [13] and Bianchi and Zhang [5], developed specifically for z-pins. Structural similarities between stitches and tufts make the first two models particularly relevant for the purposes of this paper. The aim of these models is the prediction of the bridging function characteristic of a specific z-reinforcement. Once determined, this bridging function, or bridging law, can be represented at the meso-scale by employing discrete non-linear springs, homogenised cohesive zone models or combinations of the two [14][15][16][17][18].
Jain et al. [9,10] proposed a simple analytical model for a single stitch bridging a crack under mode I loading conditions. The TTR unit was assumed perfectly cylindrical and in frictional contact with the surrounding laminate. Cox et al. [11,12,19], extended this model to inclined tows, considering mixed-mode and mode II loading conditions. In these models the pre-debonding stage for the stitch is neglected. With focus on mode I loading, a frictional slip zone ðl s Þ is considered to start developing immediately along the stitch under tensile load. For interconnected stitches, when l s equals the embedded length of the stitch L, the load and displacement keep increasing until the thread reaches its tensile strength and breaks. In reference [11], Cox and Sridhar also accounted for the anchoring action of the horizontal thread segments of the stitch. When l s = L, a further load increase is considered to be absorbed by these thread segments, which are pulled into the composite contributing to the axial displacement of the vertical thread.
Here, we present an analytical micro-mechanical model developed to capture the experimentally determined bridging mechanisms for tufts, as discussed in Section 2. The tuft is represented as a rod of cross-sectional area A, embedded in an elastic medium. It is inserted vertically and bridges a delamination at the mid-plane of a laminate. Neglecting, at this stage, the complexities of the tuft geometry and the consequences these may have on its mechanical response, the same damage mechanisms are considered to develop in both halves of the tuft. Under symmetric loading, this assumption allows to model only half of the studied domain (see Fig. 4). Unlike in the existing models [9][10][11][12], here the z-reinforcement is assumed to have an initial elastic stage, during which it is fully bonded to the composite. Where bonded, the tuft is capable of transferring load to the material it is embedded in. This load transfer causes a localised perturbation of the stresses and strains of the composite surrounding the tuft. Adopting the fundamental assumption of Cox's shear-lag theory [20], we define the shear  force at the tuft interface as a linear function of the relative displacement between the tuft and the non-perturbed composite: p s ¼ k z ½w T ðzÞ À w C ðzÞ, where k z is a parameter defined by the elastic properties of the material around the tuft, w T ðzÞ is the axial displacement of the tuft and w C ðzÞ the through-thickness displacement of the undisturbed composite. As soon as a characteristic load threshold is reached, the tuft starts to debond. This process of progressive debonding is treated as a crack propagation problem, whereby a cylindrical fracture grows along the tuftcomposite interface [21][22][23]. Along its debonded length, the tuft stretches and slides, opposed by an assumed uniformly distributed friction force, p 0 [9,11,12,19]. This fracture mechanics-based approach to describe the debonding process is one of the main differences between the presented model and those available in the open literature, where the energy dissipated for crack propagation is generally neglected. Fully inserted tufts, as those considered in this paper, are constrained to the top and bottom surfaces of the composite they are embedded in, in the same way as stitches. Neglecting any compliance of the composite surfaces, this constraint, exerted by the horizontal thread segments on one side and the thread loops on the other, is represented in the model as an encastre support for the tuft. For simplicity, any possible interaction between neighbouring tufts has been neglected in the formulation of the model.

Governing equations
The response of the tuft has been divided into two stages called bonded regime and progressive debonding regime, respectively.
In the bonded regime, the axial equilibrium of the tuft requires where E is the equivalent Young's modulus of the tuft modelled as a cylindrical rod, A its cross-sectional area, w T ðzÞ is the tuft axial displacement and w C ðzÞ the through-thickness displacement of the composite volume unaffected by the presence of the tuft. Since no forces are applied to the composite at the delamination plane, that is a free surface, the through-thickness deformation of the unperturbed composite is zero ðdw C ðzÞ=dz ¼ 0Þ and Eq. (1) can be rewritten in terms of the relative displacement between the tuft and the composite, wðzÞ, as Having established that wð0Þ ¼ 0, since the tuft is constrained to the external surfaces of the composite by means of the surface seams and loop, and EA dwðzÞ dz where P is the force that needs to be applied to the tuft at the delamination plane in order for half of the system to be in equilibrium (see Fig. 4), the solution to Eq. (2) is where a ¼ ½k z =ðEAÞ 1=2 . Eq. (4) defines the axial displacement of the tuft as a function of the axial coordinate z and the external load P (see Fig. 4). At the delamination plane ðz ¼ LÞ, this displacement becomes which corresponds to half of the relative displacement between the surfaces of a bridged crack. If, on the other hand, we assume the tuft to be debonded over a generic length l d , the axial equilibrium of the tuft is defined by the following system of equations: where the subscripts b and d refer to the bonded and debonded portions of the tuft, respectively. Eq. (6a) must satisfy the geometric and natural boundary conditions defined as with b P ¼ P À p 0 l d , following from simple equilibrium considerations. The tuft axial displacement in the domain z ¼ ½0; L À l d , is The solution to Eq. (6b) can now be obtained by imposing: whereŵ follows from Eq. (8) for z = L À l d . Satisfaction of the boundary conditions implies . Schematic of the problem described by the analytical model. Despite the irregular shape of the tuft and the fact that often failure occurs at a certain, but small, distance from the delamination plane, for simplicity only half of the system is studied, under the assumption of symmetric load and geometry.
which defines the z-displacement of the tuft for L À l d 6 z 6 L, in the regime of progressive debonding. At the delamination plane, this displacement is It can be noted that for a zero debonding length (l d = 0), Eq. (11) correctly coincides with Eq. (5). On the other hand, in the case of a fully debonded tuft, the displacement at the delamination plane becomes which is independent from the elastic foundation stiffness, k z .

Debonding criterion
When the process of debonding starts, a crack is considered to propagate along the lateral surface of the cylindrical tuft. The energy required for an infinitesimal increment of the debonded length is G C;int l p dl d , where G C;int is referred to as 'fracture toughness of interface' and l p is the cross-section perimeter of the zreinforcement. Since friction is assumed to act along the debonded length of the tuft, for the principle of stationarity of the total energy, the energy available for fracture growth is where G int is the equivalent mode II strain energy release rate of the interface, W ext is the total work done by conservative forces on the system, U fr is the energy dissipated by friction and U e is the elastic strain energy. It follows that the energy available for a unit increment of the crack surface is The numerator of Eq. (14) is the energy difference between a state in which the debonded length is l d and a second state in which the debonded length is l d þ dl d , as illustrated in Fig. 5. Following this criterion, and with few manipulations (see Appendix A), Eq. (14) can be re-written as By substituting the expressions of dw Eq. (15), one obtains For l d ¼ 0 and G int ¼ G C;int , Eq. (16) provides the critical load marking the onset of debonding The critical load depends on both the fracture toughness of interface and the friction force acting along the debonded portion of tuft. To initiate debonding, the applied load needs to be high enough to provide sufficient energy for the crack to propagate and to equilibrate the forces, including p 0 , acting along the TTR.
During debonding, G int ¼ G C;int and P ¼ P C þ DP. The substitution of G C;int in place of G int in Eq. (16) provides a relation between load and debonded length during the regime of progressive debonding. For a given l d , the corresponding load value is determined:

Identification of model parameters
The proposed model requires the evaluation of a series of geometrical and physical input parameters. The cross-sectional area and perimeter of the tuft, A and l p , and the composite thickness, 2L, can be determined from micrographs. The equivalent Young's modulus, E, the friction force, p 0 , the foundation stiffness, k z , and the equivalent mode II fracture toughness of interface, G C;int , are physical parameters requiring ad À hoc experiments to be determined.

Tufting thread characterisation
The most common architecture of commercially available continuous carbon fibre threads consists of a given number of fibre bundles, called yarns, interlaced together so that each yarn follows a helical path. The structure of the tufting thread resembles that of a rope, without a central core. When the thread is embedded into a preform, the preform architecture controls the final shape of its cross-section and the length of its perimeter, l p . For a specific preform, l p and A need to be measured from micrographs of tufted samples. As for the parameter defining the equivalent elastic modulus of the tuft in the model, E, this is determined by the elastic properties and volume fractions of the fibres and matrix constituting the tuft, as well as by the in situ geometry of the tuft. The latter is influenced both by the in-plane fibre architecture and stacking sequence of the preform and by the process of consolidation. This makes the identification of the final shape of the tuft, and consequently of the equivalent elastic modulus E, not trivial, which is the reason why E has been selected as one of the calibration parameters of the model.

Elastic foundation stiffness
As explained in Section 3, stress transfer between the tuft and the composite over the bonded length of the tuft has been taken into account by modelling the material around the tuft as a bed of tangential springs of stiffness k z [20]. A way of approximately quantifying this stiffness parameter is by means of Finite Element (FE) calculations. The model illustrated in Fig. 6 has been used to characterise k z as a function of the tuft equivalent elastic modulus E and the tuft embedded length, L. In this model, the tuft is described as a cylindrical rod of cross sectional area A (defined in Section 3.1), fully bonded to the surrounding medium. The FE analyses have been carried out using Abaqus 6.14. A static-general step has been implemented using C3D8 elements. A linearly varying displacement field has been applied at the nodes of the tuft, as shown in Fig. 6, and the reaction force on the top surface of the composite has been requested in output. Elastic transversely isotropic properties have been assigned to both the tuft and the composite. Fig. 7a shows that k z has a strong dependency on the thickness of the composite, which in the analysed models corresponds to the embedded length of the tuft, L. This result would suggest that the adoption of a constant k z value throughout the process of progressive debonding of the tuft is not valid, and it can only be accepted as a consequence of the low sensitivity of the model to k z in that stage of the bridging response, as demonstrated by the plot in Fig. 7b. k z varies not only with the embedded length of the tuft, but also with the boundary conditions of the tufted unit-cell. This means that its value would need to be re-evaluated for every new tufted structure under investigation. However, Fig. 7b shows that any variation of k z has little influence on the constitutive law of the tufted interface, within the limits of experimental repeatability. This implies that the curve calibrated on the results of the single-tuft tests has general validity, i.e. it can be used to predict the delamination behaviour of any tufted composite structure manufactured with the same materials and fibre architecture of the single-tuft specimens, independent of any variation of the boundary conditions. The effect of a different composite thickness can be predicted with the use of the analytical model only. Fig. 7c further shows that k z is also a function of the equivalent longitudinal modulus of the tuft, E; hence, in order to be able and determine the value of k z for a particular system, it is necessary to know the values of L and E for that specific system. However, since E is a calibration parameter, its value cannot be known a priori. This obstacle can be overcome using the empirical formula below, which establishes a relation between E and k z , based on the results of Fig. 7d: For the system described in Section 2, the coefficients a and b have been determined to have values 0.1571(AE0.0196) mm and 6976 (AE1528) N/mm, respectively. Eq. (20) has been used during the process of model calibration to determine which value of k z corresponded to the E value providing the best fit.

Friction force
The friction force, p 0 , acting along the debonded length of the tuft, is influenced by two different mechanisms. First, thermal residual stresses due to post-cure cool down affect the interfacial friction between the tuft and the composite. Secondly, the irregularity of the tuft longitudinal profile promotes mechanical interlocking between the tuft and the surrounding resin material. Since these irregularities have not been accounted for in the proposed model formulation, their effect reflects on the value of p 0 , as well as on the values of the other material parameters of the model. The quantification of p 0 , which in this study is assumed to be uniformly distributed [9,12,11], has been carried out by means of pull-out tests on single-tuft specimens in which the tuft loop has been machined off. The load-displacement curves obtained for the specimens described in Section 2 are reported in Fig. 3c. Given the average maximum force P PO in the tuft at pullout initiation, the tuft cross-section perimeter l p , and the effective pull-out length, L PO , the equivalent uniform friction stress opposing the sliding of the tuft can be calculated as [6] s 0 ¼ For the considered composite system, s 0 = 36 (COV = 16.7%) MPa and the equivalent friction force distributed along the tuft, defined as p 0 ¼ s 0 l p , is 80.6 N/mm.

Fracture toughness of interface
The process of debonding has been modelled as the propagation of a crack at the interface between the tuft and the composite, in pure mode II. This is a simplification derived from the assumption of the tuft as a straight rod. In reality, tufts often display curved profiles with irregular perimeters, which translates into a variable stress-state at the tuft-composite interface and a debonding crack propagating in mixed-mode. The mode mixity varies during crack growth, as therefore does G C;int . Depending on the in situ shape of the tuft, debonding may affect only a portion of its lateral surface, i.e the portion subjected to the highest tensile and shear stresses.
The assumption of a debonding crack propagating uniformly and symmetrically in the two halves of the delaminated tufted composite is expected to reflect inevitably on the input value of G C;int , which is therefore difficult to be determined a priori. Based on these considerations, the fracture toughness of interface has been selected, together with E, to be calibrated on experimental data sets obtained from single-tuft tests. gressive debonding. The critical load, P C , and corresponding critical displacement, w C , are determined by means of Eqs. (17) and (5). The bridging law in the bonded regime is obtained by discretising the displacement range up to w ¼ w C , and evaluating the corresponding values of P through Eq. (5). The load-displacement curve is linear up to P ¼ P C . The length of debonding, l d , is given as input to the second part of the script. Given the values of l d ; P can be determined from Eq. (18) and the corresponding displacement obtained by means of Eq. (11). In this regime, the bridging law is non-linear due to the energy dissipated by the system for increasing the debonded length of the tuft and for overcoming friction. Table 1 provides a summary of the model input parameters, with their experimentally determined values. The only two parameters requiring calibration are the equivalent axial elastic modulus of the tuft, E, and the fracture toughness of the tuft-composite interface, G C;int , as discussed in Section 4. These unknown parameters have been identified by means of the Genetic Algorithm (GA) routine in MATLAB [13]. The cost function to be minimised in the GA optimisation scheme is

Calibration of the model
where 2 Pw is the relative mean square error associated with the averaged load-displacement data points, and 2 W b and 2 W d b are the relative mean square errors of the work done by the external forces in the bonded and progressive debonding regimes (neglecting the work of partial frictional pull-out, if present), respectively. The optimisation has been carried out considering 50 combinations of populations of 50 individuals each. An initial range of values has been assigned to each calibration parameter: [0.1 240,000] MPa and [0 2] kJ/m 2 for E and G C;int , respectively. The upper bound for the tuft equivalent Young's modulus has been chosen as the tensile modulus of the impregnated carbon fibre yarns forming the tuft (Tenax Ò -J HTA40 H15 1K 67tex [24]), whereas G C;int is allowed to vary in a range of carbon-epoxy laminates common values of fracture toughness [25]. The values of the calibrated parameters are reported in Table 1 and the resulting calibrated bridging law is displayed in Fig. 8. As mentioned in Section 2.3, the failure load of the tuft is a function of the tuft in situ topology and cannot therefore be determined from independent tensile tests on impregnated tufting threads. Hence, in this paper, an experimentally-based failure criterion deduced from mode I testing of single-tuft coupons has been selected. In particular, the calibrated bridging law of Fig. 8 has been interrupted at the average measured deformation at failure of the tested carbon-fibre tufts (see Table 1). The debonded length corresponding to this deformation is l s = 1.77 mm, 88.5% of the embed-  ded length of the tuft. As for the energy involved in the process, the current model predicts that the energy loss due to friction and to debonding contribute to approximately 40% and 13% of the total energy of the system, respectively. Thus, although the energy dissipated by friction is dominant, the energy spent for crack propagation at the tuft-composite interface is not negligible. This supports the modelling strategy detailed in Section 3.
The determined E and G C;int parameters can now be used to predict the response to delamination of tufted meso-scale structures with the same characteristics, in terms of tufting thread material and preform material and architecture, of the single-tuft specimens the model has been calibrated on.

Mode I double cantilever beam test on tufted specimens
The new micro-mechanical analytical model has been validated by comparing meso-scale FE simulations to delamination tests of DCB specimens. Fig. 9 shows the schematics of sample dimensions and tuft row positioning used in the DCB tests. The longitudinal axis of the DCB specimen is parallel to the 0°ply orientation at the mid-plane of the specimen. The tufts were inserted normal to the laminate plane, in a square pattern with 5.6 mm Â 5.6 mm spacing, resulting in an areal density of tufts of 0.86%. The tufting seams were aligned orthogonal to the 0°ply orientation, with the first tuft row positioned 15 mm from the end of the release film, covering a total delamination length of 100 mm. A 10 lm thick PTFE film, at the mid-plane of the preform, extends 61 (AE0.5) mm from the edge producing an initial crack length of 50 mm in the consolidated panels. Each panel was injected using the VARTM process, resulting in a consistent laminate thickness of 4.0 (AE0.01) mm and a fibre volume fraction of 56%.
Delamination tests on the DCB specimens were performed on a universal test machine Instron 5500R, with 5 kN load cell, at a constant crosshead displacement rate of 1 mm/min. In absence of standardised procedures for testing through-the-thickness reinforced bi-directional laminates, mode I tests were performed following the guidelines developed for 'unreinforced' unidirectional composites [26]. Prior to testing, the initial pre-crack was propagated for 5 mm from the film insert. The new crack tip ensured the measured initiation values to be independent of the insert film thickness. Finally, for monitoring purposes, one side of each beam was coated with a thin layer of brittle white spray-on paint and marked with millimetre increments to enable visual monitoring of the crack tip position. On the opposite side, a black and white speckle pattern was applied for monitoring of crack propagation and for measurement of opening displacements with a Digital Image Correlation (DIC) system (Limess VIC 2D and VIC 3D).
To avoid premature flexural failure of the 2 mm thick beams of the tufted DCB specimens (see Fig. 10), 3 mm thick strips of cured unidirectional composite were bonded to both sides of the samples, as suggested by other authors [27,28]. The unidirectional laminate had 10 layers of carbon fibre uni-weave fabric (OCV Technical Fabrics™, 12k Grafil 34-700) with an areal weight of 310 g/m 2 , injected with MVR444 epoxy resin.

Modelling tufted interfaces at the meso-scale
Cohesive zone modelling is the approach selected to model the delamination behaviour of the tufted DCB specimens of Section 6.1. The application of this approach implies the use of cohesive elements to simulate the delamination behaviour of both the unreinforced interface and the discrete through-thickness reinforcements [5]. The former is governed by a bilinear traction-separation law in which the onset of damage is related to the interfacial strength, i.e., the maximum traction on the traction-separation jump relation, and the propagation of delamination follows a fracture mechanics criterion. When the area under the traction-separation jump relation is equal to the intrinsic fracture toughness of the resin, the tractions are reduced to zero and a new delamination front is formed [29,30]. As for the cohesive elements used to simulate the behaviour of the tufted regions, their constitutive law has been obtained using the micro-mechanical model presented in this paper (see Fig. 11). This constitutive law describes the initial localised elastic deformation of the tufts and material surrounding them (bonded regime), followed by progressive debonding of the tufts up to failure. It should be noted that the localised deformation taking place in the material around the bonded tufts can only be accounted for via the constitutive law of the interface elements, as the meshes of meso-scale models are generally too coarse to capture it. To avoid any further unrealistic through-thickness deformation of the material above and below the bridged interface elements, columns of rigid elements have been inserted in the arms of the DCB model at the tuft locations. The implementation of the mode I constitutive law of tufts into this finite element framework has been achieved via a user-defined constitutive formulation in which an arbitrary normalised traction-separation law can be prescribed.

Model description
The tufted DCB specimens of Section 6.1 have been simulated in the explicit finite element code, LS-DYNA v971-R7.1. 8-node selectively reduced-integration solid elements with hourglass control were employed to simulate both the composite NCF and UD beams with at least 6 elements in the through-thickness direction. This was done for both sub-laminates in order to capture the bending stiffness and laminate rotations. Symmetry conditions were exploited by modelling only one half of the specimen and taking into account the periodic arrangement of the tufts. In the tufted beam, tufts were modelled discretely with rigid elements and a tuft-to-tuft distance of 5.6 mm (see Fig. 12). The displacement rate was approximately 1 mm/s (defined initially by a smooth ramp rate followed by a constant velocity boundary condition imposed at the load application points), which made dynamic effects negligible. Mass scaling was used to reduce the solution time. A summary of the elastic homogenised properties assigned to the composite beams and a list of the cohesive parameters used to describe the mechanical responses of the un-tufted and tufted interfacial regions are provided in Tables 2 and 3, respectively. Due to the presence of the stiffening UD tabs, the average initiation fracture toughness measured for the tested DCB specimens has been used as the mode I matrix fracture toughness and its value is G C;resin I = 0.257 kJ/m 2 (COV = 5.9%). The fracture toughness of the tufted interface, G C;tuft I , has been calculated as the overall bridging work of the tuft, W I , divided by its cross-sectional area A. Following a mesh convergence study, in which the results were checked during the elastic and crack propagation phases, a minimum cohesive element length of 0.25 mm in the direction of crack propagation was adopted. This allowed the model to have three cohesive elements in the fully developed process zone ahead of the crack tip, in agreement with current meshing guidelines for cohesive zone models [31,30]. . It is clear that the delamination initially propagates through the untufted region, up to the first tuft row. Here, the localised bridging of the tufts causes a significant load increase, i.e. almost 130%. After failure of the first row of tufts, a large-scale bridging zone fully develops. For the laminate and DCB configuration considered here, the bridging region comprises two rows of intact tufts. With further opening displacement, the crack propagates unstably and is arrested by the presence of the localised reinforcements, causing the rows of tufts to progressively stretch and fail [15]. In this region of the load-displacement curve, the FE model underestimates the bridging response of the tufts, showing a more pronounced saw-tooth behaviour compared to the experimental curves. This could be corrected by considering the effect of fibre bridging, caused by the in-plane fibres of the fabric, on crack propagation (see Fig. 9).

Conclusions
A multi-scale modelling framework for predicting the mode I delamination behaviour of tufted composites has been presented and validated. It applies to the case of structures tufted through their entire thickness, with loops long enough to anchor the tufts and prevent them from pulling-out under load. The proposed multi-scale approach is based on a micro-mechanical model describing the mode I response of bridged interfaces, coupled with a meso-scale cohesive zone formulation for tufted structures. In the paper, this strategy has been applied successfully for the prediction of the delamination behaviour of 4 mm thick double cantilever beam specimens containing an array of tufts at 0.86% areal density.
The modelling approach presented has the advantage of requiring the calibration of only two parameters at the micro-scale (i.e. at single-tuft level), namely the equivalent Young's modulus of the cylindrical tuft and the equivalent mode II fracture toughness of the tuft/composite interface, and it captures the main bridging mechanisms observed for tufts under tensile loading, i.e. initial linear elastic deformation of the bonded tuft, followed by a progressive debonding event and frictional sliding of the debonded tuft, which are responsible for the non-linearity of its response. Despite the level of idealisation, the micro-scale model presented, when implemented into interface elements at the meso-scale, was capable of predicting the stiffness and first load peak of the tested DCB specimens within the limits of experimental repeatability.

Acknowledgments
The authors would like to acknowledge the Engineering and Physical Sciences Research Council and Rolls-Royce plc for their    support of this research through the EPSRC Centre for Doctoral Training in Advanced Composites for Innovation and Science (EP/ G036772/1) and the University Technology Centre (UTC) at the University of Bristol, UK. Original EPSRC studentship to Treiber (EP/F037937) is also gratefully acknowledged. All underlying data to support the conclusions are provided within this paper.
Appendix A. On the formulation of the debonding process Eq. (16) has been obtained by imposing the limit for Dl d ! 0 to the energy variation characterising the system described in Fig. 4, when the length of the debonding crack increases of Dl d . With reference to Fig. 5, in state 1 the energy of the system is in which N b ðL À l d Þ ¼ N d ðL À l d Þ for continuity.