Experimental and numerical investigation of interface damage in composite L-angle sections under four-point bending

Curved laminates in aero-structures, such as the L-angle sections where webs and flanges meet, are prone to delamination due to high interlaminar stresses in these regions. Some efforts to investigate delamination in these structures can be found in the literature but commonly structures are limited to unidirectional layups or modelling approaches are constrained to the cohesive element based methods. In this work, multi-directional L-angle laminates were manufactured using unidirectional prepregs and tested under four-point bending load conditions to examine the interface damage. Acoustic emission technique was used to assist the capture of damage initiation and propagation. Three interface modelling strategies for predicting delamination, namely cohesive element, cohesive surface and perfectly bonded interface were used in the numerical study. The interface damage behaviour was successfully predicted by the simulation methods and differences among the strategies were compared.


Introduction
Composite materials have increasingly been used to meet the growing need for lightweight structures, primarily in aerospace, wind energy, automotive and construction sectors, due to their high specific stiffness and strength along with high energy absorption. The large volume usage of fibre reinforced composites in recent aircraft such as Airbus A350, A320 and Boeing B787 has demonstrated the maturity and acceptance of these advanced materials. Due to their absence of through thickness reinforcement, laminated composite materials are prone to delamination, which reduces the material properties and is considered as one of the most catastrophic failure mechanisms for composite structures. 1 This is particularly obvious for some curved parts connecting different flat panels, 2-5 e.g. L-angle sections where flanges and webs meet in large aircraft structures, as the primary loading conditions for these regions are folding/unfolding forces and/or moments. Characterisation and modelling and of the interlaminar damage behaviour of such structural elements is essential in designing of aircraft composite structures. Tensile testing in the through-thickness direction to waisted samples is the basic method to determine the interlaminar properties. In order to avoid premature damage near the clamps, such tests are only suitable for very thick samples. 6 Some inplane uniaxial tests can also be used to determine the out-of-plane properties but they are only applicable to some samples with specific layups, which can induce delamination at the free edges under tension. 7,8 Without losing generality, the four-point bending tests on L-angle samples 9 present less complexity in the mechanical tests and coupon preparations, which have been standardised 10 for engineering reference.
Consequently, in this work the four-point bending test was selected to determine the interlaminar damage behaviour of the laminate L-angle sections.
Delamination in composite materials is a progressive process 11 and therefore different methods have been developed for modelling of damage initiation and propagation. The Virtual Crack Closure Technique (VCCT) based on fracture mechanics calculates energy release rate, 12 assuming that the energy needed to separate a pair of surfaces is the same as the energy needed to close the same surfaces. It is only applicable for modelling of delamination propagation due to the prerequisite which is the need for the presence of pre-existing cracks. Another technique to simulate the delamination in composites is the cohesive zone model (CZM), which relates interfacial tractions to separations through a softening law. 12,13 CZM is able to predict both damage initiation and propagation without the limitation faced by the VCCT technique and therefore is widely employed in modelling of composite delamination. The CZM method has been successfully employed to model structures of different complexities, ranging from simple beams for model validation [14][15][16] to applications on complex composite structures. [17][18][19] The CZM is currently the mostly widely used method to model delamination in composite structures due to its widespread implementation in commercial FE solvers. However, care should also be taken when using the CZM method because of its multiple model parameters 20 and mesh dependency effect. 15,21 The accuracy of the CZM depends on having a sufficient number of elements within the cohesive zone. 15 Using a coarse mesh might not be sufficient to capture the accurate stress distribution along the cohesive zone and might lead to unreliable results. Mesh design should be based on the consideration of cohesive zone length governed by several material properties as equations presented by Turon et al. 21 and Yang and Cox. 22 Typically, the cohesive zone length for an aerospace carbon/epoxy laminate is approximately 0.3 to 1 mm as suggested by Cox and Yang. 23 To minimise meshdependency, appropriately reducing the interface strength in CZM was found necessary when a large structure is meshed with an element length slightly greater than the above typical cohesive zone length, because a lower interface strength leads to a longer cohesive zone length which could incorporate more interface elements to capture the cohesive zone stress distribution. 15,21 Inaccurate results were also observed as a result of specifying an excessively low interface strength due to earlier damage initiation at the interface leading to overly compliant structural behaviour. 15 The cohesive zone method can be implemented based on either a cohesive element or a cohesive surface (a form of contact), and a significant number of applications of both types have been reported in previous research. The cohesive element method is able to model the behaviour of finite-thickness bond lines, such as debonding in adhesively bonded joints. 24 The thickness of a constituent interface in composite materials is usually negligibly small. When using the cohesive element, it is necessary to reduce its thickness to a meshing-allowed minimum to avoid spurious compliance, or to use a zero-thickness cohesive element. However, it is difficult to generate zero-thickness cohesive elements or even finite-thickness cohesive elements at constituent interfaces for complex structures, hence the surface-based cohesive modelling method has been more widely employed for modelling delamination in complex composite structures.
The prediction of delamination in L-angle composite sections has attracted broad research interest in recent years. 16,25,26 Nguyen et al. 16 experimentally and numerically investigated the failure behaviours of L-angle sections under four-point bending. 3D cohesive elements at the interface between two adjacent layers were used to predict the delamination of the specimen. Good agreement of predicted load-displacement curves with experiment was obtained after a parametric study on cohesive parameters. A similar study was conducted by Cao et al. 25 based on cross-ply laminates. Due to the symmetry in the adopted lay-ups (0 8 /90 8 only), zerothickness 3D cohesive element FE models were employed to represent a quarter of the L-angle arc. As cohesive elements were also embedded at the intra-laminar interfaces to account for matrix cracking, the mesh was extremely refined so that only a section of about 10 mm of the leg of the L-shape specimen was kept in the model to apply moment. The predicted failure bending moment (no load-displacement response due to model simplification) and failure locations were in good agreement with experiment. A plane strain FE model with 2D cohesive elements was used by Geleta et al. 26 to study the delamination behaviour of L-angle laminates but 3D effects were not considered. It is concluded from the above literature that although curved structures have been investigated in terms of delamination prediction, there is still a lack of a detailed comparison between the use of the 3D cohesive element method and the cohesive surface method to predict delamination in multi-directional composite L-angle sections. In addition, damage initiation is sometimes of interest at the design stage and the validation of using a computationally efficient method to predict delamination initiation has been ignored in the literature.
In this paper, experimental study is carried out to characterise the interface damage behaviour in a multidirectional L-angle composite section under four-point bending. Different interface modelling strategies for delamination modelling will be investigated and compared. This includes the use of the cohesive element method and the cohesive surface method to predict the progressive interface damage behaviour as well as using a perfectly bonded interface to predict delamination initiation. With the help of acoustic emission and digital image correlation techniques, both interface damage initiation and propagation can be accurately captured in the experiment and used for model validation. The predicted results from different modelling strategies will be compared and validated against experimental results.

Materials and fabrication of L-angle specimens
One L-angle panel (then cut into three specimens) was prepared for experimental characterisation and validation of the numerical models. The material used for the manufacturing was EasyComposites XPreg XC130 which is a unidirectional prepreg consisting of T700SC carbon fibre and a proprietary epoxy resin. The prepreg has a nominal cured ply thickness of 0.3 mm and nominal fibre volume fraction of 67%. All laminates were manufactured by hand layup with room temperature debulking after the 1st, 4th, 8th and 12th ply for 15 minutes. An aluminium male tool with radius of 6.35 mm was used to manufacture L-angles and a flat aluminium tool was used for UD specimens. The prepreg layups were cured in the autoclave with the cure cycle recommended by the manufacturer: heating to 120 C with a 2 C/min ramp, 1 hour hold, ramp to 130 C with 0.3 C/min followed by a hold for 2 hours. The panel was cured under 6 bar pressure with full vacuum applied. The manufactured panel was cut into three specimens of 25 mm width using a diamond saw. Dimension of individual L-angle specimens are given in Table 1.

Mechanical testing methodology
The three specimens were tested in four-point bending using the setup shown in Figure 1. Due to the absence of testing standards for multi-directional laminates, experiments followed ASTM D6415 standard for unidirectional laminates where possible (distance between the rollers was reduced to ensure that failure happens at a displacement lower than 5 mm). The testing rig comprised loading rollers of 5 mm radius with span between the pairs of rollers being equal to 60 mm and 35 mm. Testing speed was 1 mm/min. Two acoustic emission (AE) sensors were attached to the bottom roller supports to detect damage events occurring in the specimens. The AE system from Physical Acoustics Corporation was used to acquire signals from the sensors using 0 dB gain and 45 dB threshold for signals. AEWin software was used to collect results.
Experimental force-displacement curves along with the absolute energy of damage events as detected by AE are given in Figure 2. In addition to the curved samples, two flat plaques were tested to determine Young's moduli of the material. The layups of the plaques were [0 ] 6 and [90 ] 12 , which were used for testing in fibre and transverse directions, respectively. Nominal thicknesses of the cured panels were 1.9 mm and 3.8 mm, respectively. Five specimens of each orientation were prepared and tested. Experiments followed BS EN ISO 527-1 standards where possible (thickness of samples was greater than recommended). The nominal width of specimens was 25 mm. Aluminium tabs were bonded to the ends of specimens such that testing grips could be attached, with the span between the tabs being 150 mm. The specimens were tested at a speed of 1 mm/min. Strains were measured using a clip-on extensometer. It was found that the Young's modulus in the longitudinal direction was 114.06 AE 1.72 GPa and the Young's modulus in the transverse direction was 7.34 AE 0.03 GPa.

Finite element modelling strategies
In this section the FE modelling strategies for the L-angle sample are presented, along with the details of the FE models. The models are arranged in three categories with decreasing fidelity levels. The most complicated model features the 3D interlaminar damage model, which is typically implemented using the cohesive zone model, while the simplest model uses single layered shell elements without any interfaces, and thus only provides global stiffness information for the structure. For the sake of simplicity, in all models, the rollers which support or apply loads to the L-angle are assumed to be rigid bodies. Frictionless contact is defined between the L-angle and the rigid rollers. All simulations are based on the finite deformation theory. Figure 1 (right) shows the modelling configuration. For convenience and following the precedent of the models in the literature, the orientation is inverted from the experiment. The specimen dimensions and geometry and other parameters for the experiment are reported in Table 2. Note that the thickness and width of the L-angle in the FE models are the averages of those of the three specimens.
The material properties used for the plies are from two sources: the in-plane Young's moduli are set according to the experimental data, while the shear modulus and Poisson's ratio are chosen from the literature, 27 viz. G 12 ¼ 4:23 GPa, 12 ¼ 0:309 and 23 ¼ 0:4. For all the modelling strategies introduced below, the material is modelled using a linear elastic behaviour.

Full 3D model with cohesive interfaces
The full 3D model incorporates all the laminae, covering their individual ply properties and orientations. In this case, all elements modelling the plies are   continuum elements which are based on 3D continuum mechanics. The ply interfaces are modelled with cohesive zone models. Two typical 3D implementations provided by Abaqus are the cohesive elements (COH3D8 for hexahedra) and the cohesive contact surfaces. Both implementations use the same cohesive zone modelling principle based on the Crisfield theory. 28,29 A general introduction on the cohesive zone model used in this study is given in the Appendix. The bilinear traction-separation law, which is geometrically the most simple form to implement, has been used for delamination analyses. Results for global load displacement response are relatively insensitive to the exact shape of the traction-separation curve, provided that the correct interfacial strength and fracture toughness are applied. 22 In the absence of experimental data for the interface properties, values (Table 3) from the literature for carbon/epoxy interface are adopted. It is noted that the interface stiffness is a numerical value introduced by the CZM and the selection of this value requires extra care so that it doesn't cause spurious compliance or convergence issues. In addition, the predicted damage initiation loads are not sensitive to the exact value of the numerical interfacial strength values if they are in the right range. As this study deals with a large composite structure, a parametric study on the cohesive properties is computationally difficult to perform. The readers are referred to literature 15,19 for detailed discussions on the effects of cohesive parameters on structural response.
Using the cohesive element, a finite thickness of the cohesive layer can be properly modelled, however, such a thickness could reduce the global stiffness of the structure when a problem-dependent threshold is reached. In this study, the thickness of laminate is about 5 mm, and the cohesive element thickness is chosen to be 0.01 mm. A local view of this configuration with two plies is illustrated in Figure 3. It is also possible to create zero thickness cohesive elements when extra constraint on conformality between the meshes of neighbouring plies is imposed.
In practice, the nodes on opposing interlaminar interfaces may not necessarily be conformal due to the difference in size between cohesive and lamina elements (as demonstrated in Figure 3). This is because the mesh size of the cohesive element in the plane parallel to the shell midsurface is required to be sufficiently small to ensure the validity of simulation, and such size is usually smaller than the element size of the neighbouring lamina. In such a case, when inserting the additional of cohesive elements, there are in total four surfaces forming the two interfaces between the cohesive layer and the two connected plies. If the nodes on these surfaces do not match, additional constraints must be introduced to ensure the connection of the layers. The tied constraint is used here for its simplicity. For non-conformal interfaces, the tied constraint in Abaqus will retain only the lower number of DOFs for the nodes on either interface, and thus effectively reduces the total DOF number of the whole model. In particular, 1 088 924 nodes are defined by the mesh, while another 187 856 nodes are created by Abaqus for the roller contact surfaces.
The cohesive interface can also be modelled in a contact-like manner, which is known as "cohesive surface" or "cohesive contact" approach. For this approach, no specific element has to be created for the cohesive layer but a zero thickness contact pair. The utilities for establishing the contact pairs available in Abaqus/CAE can be used to ease the modelling process. During the solution process, the cohesive behaviour is actually implemented in a contact formulation, which will introduce additional internal DOFs. The number of those internal DOFs is proportional to the size of the cohesive surfaces as well as their mesh densities. Consequently, the total number of DOFs of the whole model will increase accordingly. In this study, there are in total 15 cohesive interfaces between the 16 plies, and the DOF number is roughly doubled during this modelling approach. In particular, there are 127 204 nodes defined by the geometrical discretisation,  while Abaqus generated additional 126 546 internal nodes to account for the cohesive surfaces and the roller contact surfaces. Due to the high computational cost for the high fidelity 3D models, the reduced integration hexahedral element C3D8R is used with enhanced hourglass controls to suppress the zero-state energy modes.

3D Model with perfectly bonded interfaces
It is well-known that the 3D models with cohesive zones are extremely expensive because of the progressive interface damage model, which limits their general usage in large composite structures. From the point of view of engineering design, a 3D model with perfectly bonded interfaces can be used for the evaluation of interface damage initiation if this is of interest at the design phase.
Two approaches are used in this study to create the perfectly bonded interface models: the "tied" interface approach uses the TIE constraint provided in Abaqus, and the shared node approach simply reduces the interlaminar interface to element interfaces. In this category, both models use the C3D8 element which provides stresses at 8 integration points per element.
In using "tied" (or surface to surface constraint) conditions to approximate lamina to lamina interactions, the requirement for coherent meshes between neighbouring lamina can be avoided. Using this approximation, interfaces are modelled as infinitely stiff and infinitesimally thin. Laminae may be modelled individually (as separate "parts", see Figure 4(a)) and meshed as such. Surface tied conditions may then be defined between neighbouring laminae. Although the application of a tied interface prevents damage initiation and evolution laws being implemented it is possible to evaluate the chosen delamination criterion (here given by equation (4)) a posteriori. To achieve this, the stress state at the interface must be approximated. In the present work, this is done by taking nodal averaged stresses at the surface where the tied condition is applied (see Figure 4(b)). It is important to note that, as averaging is applied on a part by part basis, two evaluations of the stress state (and consequently the evaluations of f, made by equation (4)) at the interface may be made (relating to the upper and lower laminae). As equation (4) does not include stress terms which may be discontinuous at the interface, it is to be expected that results relating to the upper and lower lamina interfaces are similar ( Figure 5 bears this out and suggests the use of a suitable mesh). Linearity is assumed in the interface stress states, with observations being scaled in order to estimate the time instant at which the delamination criterion expressed by equation (4) is satisfied. Evaluations have been made at multiple time instances in analyses, and linearity in the interface stress states (for tied model conditions) has been verified.
The tied interface approach is especially suitable for the interface between laminae with non-matching meshes. If it is possible to create conformal meshes for the bonded lamina, the shared node approach is more preferable due to its simplicity and relatively lower computational cost. To generate the conformal  Figure shows (a), the general tied interface modelling procedure, wherein individual lamina are modelled as separate parts (without coherent meshes) and (b), the stress projection procedure, wherein information is taken from either the lamina above or below the interface. Element integration points are marked by crosses. Projections are made using information from above the interface (the "upper ply") and below the interface (the "lower ply"). Note that a cylindrical coordinate system is adopted out of convenience, where h ¼ 0 is located on the plane of rotational symmetry for the L-angle component. mesh in Abaqus/CAE, different plies will be included in the same part, but with appropriate partitions exactly at the interfaces between the plies. Thus elements between the interface will share the same nodes in the mesh topology. In addition, the postprocessing of the interlaminar stress field and then the evaluation of failure criterion can be easily performed with Abaqus/ CAE or by using a Python script.
For the purpose of further reducing unnecessary computational costs, the 3D models in this category are simplified to half models by exploiting their rotational symmetry around the y-axis. As shown in Figure 6, the displacement of an arbitrary node a lies on the symmetric plane follows the constraint with its symmetric node b where the rotationally symmetric operator R is defined by It is trivial to verify that R is symmetric, unitary and orthogonal, namely For further details about the rotationally symmetric boundary conditions, the readers are kindly referred to Li and Reid. 30

Shell model without interfaces
This category uses the simplest FE modelling strategy, which is widely adopted for large-scale structures in engineering, especially in the sizing phase of industrial products such as automobiles, aircraft and ships. The FE models are based on shell theories, including the Kirchhoff-Love theory and Reissner-Mindlin theory, depending on the geometry related structural behaviour. In this work, the thickness of L-angles is not negligible as required for the proper use of the Kirchhoff-Love theory, thus the Reissner-Mindlin shell elements are used.
Two main classes of shell elements are provided by Abaqus: the conventional shell elements which represents a shell structure being an assembly of 2D geometric primitives immersed in the 3D space, and the continuum shell elements which are degenerated 3D elements. The element type S4 is used for the conventional shell model, which is a full model with all 4 rollers with single layer of mesh. The continuum shell model has, in fact, exactly the same mesh as the 3D element model. Therefore, it can be derived from previously built models after necessary modifications on element formulations and material section assignments. The advantage of using continuum shells is they are able to be stacked to lift the straight normal restriction from Kirchhoff-Love or Reissner-Mindlin shell theories, thus effectively avoiding the shear locking even in case that the in-plane mesh is relatively coarse. In this work, we derived the continuum shell model from the shared nodes 3D model by changing the element type from C3D8R to SC8R, thus it is also a half model with rotationally symmetric boundary conditions applied on the symmetry plane as introduced previously.
The shell models cannot provide the out-of-plane shear stress components directly, thus in this study they are only used to assess the global stiffness of the model for comparison with other models.

Results and discussion
Results of the four-point bending experiments are consistent between samples with all three of the tested specimens having almost matching initial stiffness. AE data from each of the specimens were used to detect damage initiation. The average displacement and force at which the first acoustic events happened were 3.35 AE 0.10 mm and 696.7 AE 33.9 N, respectively.
For all tested specimens, the damage initiation was almost coincident with the first small drop in the loaddisplacement curve which is seen in Figure 2. The first two consecutive AE events correspond to the damage initiation and the failure, respectively. After the first event, the occurrence of superficial cracks are observed by the recorded video. Although the two AE events, as seen in Figure 2, clearly happen at different displacements, the occurrence of superficial cracks are in fact almost simultaneous. Some of the experiments indicated that the first superficial crack to appear is the crack between (approximately) plies 3 and 4 which is then followed by the crack between (approximately) plies 12 and 13, both shown in Figure 7. These two cracks were followed by a crack between (approximately) plies 8 and 9 and a crack between (approximately) plies 10 and 11. However, it seems that the exact sequence of the crack appearance was different for each sample. Besides, it is not certain that the sequence of occurrence of the superficial cracks represents the truth, because the internal cracks were impossible to observe, while they could occur first without being seen and propagate to the surface to appear as the second or third crack. Potentially, use of a high-speed camera could provide more information about the crack sequence. A number of superficial transverse intra-ply cracks are also visible from Figure 7 after the detected damage initiation displacement (3.35 mm) by the AE. At the end of the loading, some apparent transverse cracks on the inner layer can be observed from the damaged specimens Figure 1 (left).
Unfortunately no scanning electron microscope (SEM) was performed to confirm the occurrence of transverse cracks in different layers.
A summary for the size of all run FE models are listed in Table 4. The models are run on different computer platforms according to their size. Therefore, the CPU times are reported for the readers' reference, but not for quantitative comparisons. The damage initiation and ultimate strength predicted by the FE models are compared with the experiment results as summarised in Table 5.

Full 3D model with cohesive interfaces
The full 3D models with CZM can simulate the actual delamination process in the L-angle sample and its overall response. However, this is extremely expensive in terms of computational cost. This is because, when the damage has initiated, the significant stiffness change will cause convergence to become difficult, and one has to either 1) use a much more refined mesh or 2) decrease the minimal allowed increment length. Both ways were tested for the full 3D models. The cohesive element model used the former while the cohesive surface model used the latter. The loaddisplacement curve in Figure 8 shows that the failure occurred at almost the same displacement of around 4 mm.
At the damage initiation point, the contours for the initiation function f (equation (4)) are shown in Figure 9. The damage initiates from the lateral side between plies 3 and 4 in the cohesive element model, while it initiates from the centre of symmetry between plies 7 and 8 in the cohesive surface model. This difference, however, vanishes when the damage evolves through the two models. In the finely meshed cohesive element model, the lateral side of the damaged area stops propagation and new damage initiates from the centre of symmetry between plies 7 and 8 as it does in the cohesive surface model. Eventually, the first full failure occurs at the same location for both models. The contours for the damage variable at the first full failure point (where it is seen as a drastic drop of the load) are shown in Figure 10.
The cohesive surface model is based on a special contact algorithm, which will introduce extra DOFs for the contact implementation. This required increase Avg. of 3 samples of DOF is demanding when the laminate has a large number of plies. In this study, the 16-ply laminate has 15 interlaminar surfaces, which requires an equal number of cohesive surfaces. The additional DOFs are introduced by Abaqus by inserting internal nodes. It can be seen in Table 4 that the DOF increase reaches 99.48% for the cohesive contact model, and it is the largest increase among all the 3D models. Such a significant increase prevents further refinement on the mesh regarding a likely limitation of computational costs. In fact, the cohesive contact model ranks the most expensive one among all models in this study, since it requires the use of a contact algorithm during the FE solution, including contact pair searching and contact status determination. Nevertheless, the rotationally symmetric boundary condition does not work well with the contact constraints in Abaqus/Standard, and thus prevents the use of a half model. To achieve a compromised solution, the load increment had been decreased to capture the peak load. However, as shown in Figure 8, a spiky load was observed due to numerical instabilities. This could be corrected during postprocessing, or by further investing on the computational cost.
Comparing with the experimental results, it appears that the predicted damage initiation location at interface 3/4 (see Figure 9) matches the first AE event (see first column in Figure 7), while the predicted failure location at interface 7/8 (see Figure 10), which starts internally and will later propagate to the surface, corresponds to the crack in the middle (see second column in Figure 7).
In this category, both models are of high-fidelity and thus are computationally expensive. For the cohesive element model, the separation of interfaces is represented by the deformation of the cohesive elements,  while for the cohesive surface model, the separation of interfaces must be computed within the contact algorithm which requires a significantly large number of additional DOFs. Therefore, although in this study the cohesive element model possesses the largest number of DOF (3 548 568), it takes less than half of the CPU time than the model using cohesive contact (with 571 443 DOFs).

3D Model with perfectly bonded interfaces
For the shared node model, the contour for the initiation function f (equation (4)) at possible damage initiation locations is shown in Figure 11, both the computed half model and its mapped rotationally symmetric counterpart are plotted. It should be noted that the damage is only calculated and meaningful at the interface but it is unavoidably shown in the laminate elements when creating the contour on a 3D model due to a limitation in interpolation control. The predicted initiation location at interface 7/8 agrees with the cohesive surface model, but the quadratic criterion f only reached 0.744 at location interface 3/4. Moreover, the location at interface 10/11 has f ¼ 0.962, which could correspond to the lower crack occurred in the experiment (see Figure 7). Similarly, the contour for the initiation function f for the tied interface model is plotted in Figure 12. The predicted damage initiation contour agrees approximately with the that from the shared node model, however there are obvious discrepancies in the damage locations of high f values. This is likely due to the inaccurately recovered nodal stresses of the tied interface model, which were obtained from the Abaqus built-in functionality. As a result, the predicted damage initiation load and displacement from the tied interface model were less accurate than the results from other modelling strategies. Therefore it is recommended that further investigation on the algorithms for nodal stress recovery is required if using the tied interface approach.
In terms of the computational cost, the percentage of increased internal nodes for the tied interface model is 32.39%, which is significantly larger than that of the shared nodes model (11.60%). However, due to the fact that the tie constraint removes the DOF of the slave nodes, keeping only that of the master nodes, the total DOF 378 132 is even less than the shared nodes model's 420 393. In practice, the two models are run on computers with different configuration in parallel, thus the CPU time for the tied constraint is significantly less than than shared nodes model.
In this work, the material behaviour was considered as linear elastic as delamination is the dominant failure mode under this load case, although some transverse intra-ply cracks were observed. Unfortunately no SEM was performed immediately after the 1st AE event to confirm if transverse crack occurred during the 1st AE event. In the simulation, the outstanding transverse stress level for the laminae at the 1st AE event (3.4 mm) is around 40-50 MPa (with the highest appears in the inner layer), which is below the transverse strength for such materials. 27 However, future works are suggested to perform SEMs immediately after 1st AE event and employ a material damage model in the simulation if computational resources allowed to determine the exact failure event sequence and achieve high-fidelity simulation results.
In general, this work employed both cohesive element and surface methods to successfully simulate the  Figure 11. Damage initiation function f contour from the shared node model at displacement 3.96 mm for the L-angle section. Only the region of interest is shown, along with its rotationally symmetric half. Note that the damage initiation function f is only evaluated at the interface but here is also shown on the laminate elements due to the unavoidable interpolation when creating the contour on a 3D model. progressive interface damage in multi-directional curved laminates and presented a detailed benchmark between the two approaches in terms of prediction accuracy and computational cost, which were not fully accounted for in the literature. In addition, the use of perfectly bonded interface model to evaluate damage initiation in curved laminate has been performed and validated against experimental results, showing this is a sensible method for quicker queries at the design stage as the cohesive models are computationally expensive.

Conclusions
In this work, the interface damage behaviour in multidirectional L-angle laminates under four-point bending was experimentally and numerically investigated. Experimental tests were carried out under the ASTM D6415 standard with acoustic emission technique to capture damage initiation and propagation. An exercise has been performed among various interface modelling approaches for predicting delamination in multi-directional L-angle laminates, including the cohesive element method, the cohesive surface method and the perfectly bonded interface method. The high-fidelity models, featuring the inclusion of both cohesive element and surface methods, successfully simulated the progressive damage of the interlaminar interfaces, and good agreement in the predicted forcedisplacement curve up to ultimate strength was observed. Due to its significantly high computational costs to obtain such results, it is demonstrated that for laminates with multiple interfaces the cohesive surface approach is less preferable than the cohesive element approach, as the former method requires a huge amount of additional DOFs internally for the adopted contact algorithm. In case that only the damage initiation is of interest, the simplified 3D model with perfectly bonded interfaces (shared node method) can be used for quicker queries. Although it did not accurately match the experimental results, one could identify several plausible damage initiation locations from the simulation results, and perform further subsequent investigations to achieve sensible judgements in industrial applications.

Declaration of Conflicting Interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding
The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This project is funded by the Clean Sky Figure 12. Damage initiation function contour from the tied interface model at displacement 3.01 mm for the L-angle section. Only the region of interest is shown, along with its rotationally symmetric half. Note that the damage initiation function f is only evaluated at the interface but here is also shown on the laminate elements due to the unavoidable interpolation when creating the contour on a 3D model.