). A Simplified Layered Beam Approach for Predicting Ply Drop Delamination in Thick Composite Laminates

The prediction of delamination onset is a challenging task in the design of thick tapered composite laminates, where multiple ply terminations ( “ drop-offs ” ) are present. This paper addresses the development of a global-local ﬁ nite element-based design approach for tapered laminates, whereby layered Timoshenko beam models are employed to predict delamination initiation from individual drop-offs. This modelling strategy provides a fast and conservative method for evaluating the strength of tapered composite laminates. Parametric test cases are presented in order to validate the methodology and understand its limitations. Finally, the application of the tool to a relatively thick tapered composite test specimen comprising multiple ply-drops is demonstrated. © 2016 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://


Strength prediction for tapered laminates
The reasons why through-thickness tapering is commonly adopted in lightweight structures are twofold: 1) achieving a pre-defined geometrical shape (e.g. a given aerofoil shape for a helicopter rotor blade); 2) minimizing structural weight while retaining sufficient stiffness and strength.In composite laminates, tapering always involves shedding thickness via "terminating" (i.e.droppingoff) plies at locations determined by the target geometry.Due to the ensuing geometrical and material discontinuity, ply drop-offs act as interlaminar stress risers, causing the onset and propagation of delamination [1,2].This represents the primary failure mode of tapered laminates, usually causing a severe strength knockdown compared to the nominal load carrying capability of the constituent composite material.
The design of tapered composite laminates is usually based on rules of thumb [3][4][5], which have been deduced both from experimental testing, strength/buckling optimisation and/or elementary fracture mechanics.These design guidelines involve: 1) having a minimum number of plies (thickness) dropped at any given station; 2) keeping the stagger distance between ply terminations to at least three times the dropped thickness; 3) terminating plies in order, starting from the stiffest (0°) and ending with the most compliant (90°); 4) keeping the laminate symmetric and balanced while plies are dropped, in order to avoid membrane/bending and bending/twisting coupling.However, the rules of thumb summarised above are often conflicting; for example, it is difficult to maintain symmetry and balance while dropping plies in order according to their relative stiffness.Moreover, depending on the specifics of the problem considered (e.g.overall geometry and material employed), some rules of thumb may be more important than others in dictating the final strength.Finally, the design guidelines discussed above do not allow for a strength assessment, even in an approximated fashion.

Open challenges
Despite the vast literature devoted to predicting the effect of tapering on composite strength, significant challenges still arise when trying to estimate the load carrying capability of thick laminates comprising multiple ply drop-offs, particularly in the context of preliminary design.This is due to the fact that "high-fidelity" models of tapered components require meshes at sub-millimetric scale (typically one ply thickness as characteristic element dimension), in order to resolve the stress field and the associated ERR at individual ply terminations.An illustrative example of a "high-fidelity" model for a severely tapered composite component [25] is provided in Fig. 1.The model (Fig. 1b) comprises one solid element per ply through the thickness, with zero-thickness cohesive elements on each interface to model delamination onset.The mesh was constructed from scanned images of an actual coupon (Fig. 1a), representing in detail the curvature of individual plies as well as the geometry of the resin pockets associated with individual drop-offs.These features were found to strongly influence the actual strength of tapered laminates [3,6,10,14,15,20,23,25], hence they need to be included in the models for the sake of accuracy.A model as that presented in Fig. 1 requires 1 man-day to be set up and several hours to run on a multi-core high performance computer.Clearly modelling at such level of detail is unfeasible during preliminary design, when multiple design alternatives must be evaluated.This is the primary reason why "global-local" FEA approaches have been proposed in the literature for the analysis of delamination onset/growth at/from ply drop-offs [3,13,26].A "global-local" FEA strategy may streamline the modelling procedure, but the computational cost of local ply-by-ply high-fidelity models is still prohibitive.Moreover, the actual features associated with individual ply drops (i.e.local ply curvature, resin pocket geometry) are heavily influenced by the manufacturing process and are not known a priori.Predicting the aforementioned features would require modelling of the manufacturing process before any virtual strength testing can be carried out, further aggravating the computational burden.Hence, any predictive method for delamination prediction from ply drop-offs must be robust enough to cope with manufacturing uncertainties and process variability.

Paper overview
This paper presents a novel global-local FEA framework for the strength assessment of thick composite laminates.The underlying approach is based on considering coarse meshes at the global scale, where the presence of ply drop-offs is not explicitly modelled, homogenised mechanical properties are considered and relatively few elements are employed through the thickness.The global models are linked with local FEA analyses, where individual ply drop-offs are represented as assemblies of shear-deformable Timoshenko beams.The main objective is to provide a computationally cheap method for estimating the strength of tapered laminates.The methodology proposed here is an extension of that presented in Refs.[23,24], which relied on shearundeformable Euler-Bernoulli beam assemblies.Regarding robustness, the method is formulated considering worst-case scenarios for the geometrical arrangement of each individual ply drop-off, i.e. those that correspond to the largest strength knockdowns.This is done in order to yield a conservative strength prediction.The paper is organised as follows: the beam assembly representing individual ply drops is presented in Section 2, together with the methodology for the estimation of the ERR associated with the delaminations emanated from the ply termination.Section 3 addresses the validation of the simplified modelling methodology for single ply drop-offs, which is carried out via comparing the predicted strength with that obtained from high-fidelity cohesive zone models in Abaqus/Standard.Finally, Section 4 provides two illustrative examples of application of the simplified modelling approach proposed here; these include: 1) the optimal "placement" of a ply drop off through the thickness; 2) the analysis of a thick symmetrically tapered coupon.

Problem statement
Consider a cross-section of an asymmetric internal ply-drop in a generic composite laminate, shown in Fig. 2a.The resin pocket associated with the ply-drop is assumed to have a right triangular shape, with base length L. In the tapered section, the composite is split into two sub-laminates; the "upper" segment going over the resin pocket is denoted as "belt", while the lower sub-laminate is indicated as "core" [23,24].The "drop" sub-laminate is sandwiched between the "belt" and the "core" segments.Together, the "belt", "core" and "drop" sub-laminates form the "thick" section, while the "belt" and the "core" segments, once reconnected at tip of the resin pocket, constitute the "thin" section.The fundamental assumption adopted here is that all the aforementioned sub-laminates can be described as Timoshenko (i.e.first order shear deformable) beams.Each sub-laminate is modelled with homogenised material properties corresponding to its stacking sequence.
The resin pocket associated to each individual ply drop-off is here considered as a void.This assumption is supported by experimental evidence [6,23,24], which show that ply-drop resin pockets usually are highly porous and crack at very low loading; hence, their presence can be neglected for conservative strength assessments [21].As a consequence, the "drop" sub-laminate is free from loading.Both ends of the "belt" and "core" beams are rigidly connected, allowing the transfer of an axial force P, a shear force V, and a bending moment M from the thick section to the thin one.The internal forces and bending moments of this statically indeterminate beam system, shown in Fig. 3, can be easily calculated, either from elementary Timoshenko beam theory or FE analysis.Fracture mechanics is then applied, with the assumed void at resin pocket constituting the initial crack.From the internal axial/ shear forces and moments in the beams, the ERR associated with the thick and thin section delaminations are computed by means of the approach proposed by Williams for cracked laminates [27].In the thick section, the delaminations emanated from the tips of the resin pocket are here assumed to have equal length and hence represented as a "double" interlaminar crack surrounding the dropped plies.This assumption for the thick section delamination was originally proposed in Ref. [22] following experimental observations on asymmetric tapered specimens subjected to tension and also adopted in Refs.[23,24].
The Timoshenko beam model discussed above, which takes into account the geometrical details of a ply-drop configuration, is adopted as the local model (Fig. 2c) in the global-local analysis strategy presented here.The through-thickness shear deformation accounted for by Timoshenko beam theory becomes important when the sub-laminates are relatively thick.
The applied boundary conditions, i.e. the sectional forces/moment of the thick section, for the local beam model are extracted at the ply-drop location in the global FE model shown in Fig. 2b.The mesh density of the global FE is irrelevant, as the sectional forces/moments are not meshdependent.This is because the equilibrium of forces is always guaranteed in an implicit FE solution, implying that a coarse linear elastic global FE model can be here employed without significant restrictions.The P, V and M per unit width at any cross-section in the global FE model can be calculated by integrating the elemental stresses along that section.The corresponding positive sign conventions for forces/moments are illustrated in Fig. 2b.
Both the global FE model and the local beam model are computationally cheap, making them suitable for preliminary design iterations.The global-local modelling strategy adopted here is summarised in Fig. 4; this procedure may be fully automated to account for multiple or sequences of ply-drops in the same component.In this case, one needs to run only one global FE model and as many local beam models as the number of ply-drops in the composite component.

ERR calculation
In a linear elastic fracture mechanics framework, the ERR is computed as the derivative of the potential energy per unit width with respect to an infinitesimal increment of delamination length, δa [27].The overall ERR, or G, can be partitioned into pure mode I and mode II components, i.e.G = G I + G II .Fracture occurs when G = G c , where G c is the critical mixed mode ERR.The ERR expressions associated with axial load and bending moment [23,27], as well as their corresponding mode partitioning for both the thick and thin sections, are recalled in Appendix A. However, since Timoshenko beams are considered in this work to represent the ply drop-off sub-laminates, the ERR contributions due to through-thickness shearing need to be accounted for.The expression of the ERR due to shear forces for a single delamination [27] in the thin section is also recalled in Appendix A. On the other hand, the ERR expression due to shear forces for the "double" interlaminar crack at the thick section has not been presented in the literature before and it is here derived in Appendix A.

Remarks on the conservative nature of the local models
The idealisation of a complex tapered laminate comprising multiple ply-drops into an assembly of beams with homogenised linear elastic material properties makes the analysis cheap to run.Moreover, the solution obtained is expected to be conservative, i.e. the ERR for the delaminations is overestimated, for the following reasons: a.The geometry of the resin pocket is idealised The resin pockets are idealised to be right-triangular-shaped, as sketched in Fig. 2a.A right-triangular resin pocket represents a worst-case scenario, i.e. it is more prone to delamination than other triangular shapes (especially in the thin section) [23,24].This shape idealisation also helps minimising the number of parameters needed to describe the ply-drop geometry.b.The resin pocket is considered as a void Considering the resin pocket as a void increases the stress intensity factors in the local models.This promotes thin section delamination, while the onset and growth of interlaminar cracks in the thick section are both unaffected by the material/geometrical properties of the resin pocket [18,22,24].
c. Beam elements are used in the local models Beam elements cannot capture the stress gradients in the throughthickness direction.Hence, in the local models the singular stress fields at the delamination tip can only manifest themselves as nodal force jumps at the beam neutral axes.The strain energy integrated through the sub-laminate thickness therefore gives a large delamination driving ERR.In contrast, a continuum elastic model with a cohesive zone formulation allows only the strain energy in the vicinity of the delamination tip to be dissipated when an interlaminar crack initiates.

Virtual testing of single ply drop-off coupons
The objective of this section is to evaluate the accuracy of using onedimensional beam elements to model the sub-laminates at a ply-drop location and to understand the limitations of such an approach.For the sake of simplicity, only unidirectional laminates are considered.A parametric study is performed on 8 different tapered laminates, varying the overall thicknesses, tapering angle and dropped thickness.The T300/914C composite system, having mechanical properties given in Table 1, is considered in this study.The tapered laminates will be named hereafter with the following notation: a_b_c_θ.The first three letters a, b and c, refer to the number of unidirectional plies in the "core", "drop" and "belt" sub-laminates, respectively; θ denotes to the taper angle.The 8 tapered laminates considered are: 2_4_2_27°; 4_4_4_27°, 4_4_4_14°, 4_4_4_7°, 8_4_8_27°, 16_8_8_27°, 16_4_8_27°and 20_4_8_27°.For each laminate, a total of six FE models have been set up and run in Abaqus/Standard, namely:  2b . The resin pocket is filled with the surrounding composite material.The local material axes of the elements in the tapered section are defined such that the fibre direction follows the bottom edge of the elements.This approximation might not give accurate stress distributions in the tapered section, but the difference in sectional forces/bending moments integrated from elemental stresses immediately outside the tapered section has been found to be small (b 6%).The difference becomes smaller as the laminate gets thicker.The sectional forces/bending moment at the thick section immediately to the left of the resin pocket are extracted from the model and applied as the boundary conditions at the thick section in the local beam models described in Models (iv) and (v); iv) A local beam model comprising only Euler-Bernoulli beam elements (B23).The nodal forces/bending moments of the individual beams representing the sub-laminates, as well as the reaction forces/bending moments at the thin section, are computed in Abaqus; v) A local beam model using Timoshenko beam elements (B21).
The model is the same as Model (iv), but with a different beam formulation.The ERRs computed from Model (iv) and Model (v) are compared to assess the accuracy of the two local beam models; vi) A 2D continuum plane stress FE model (Fig. 5b) with a similar mesh as in Model (i), but without the resin (i.e. a void) and cohesive elements.The delamination initiation load obtained from Model (i) is again applied as its 'far-field' boundary condition.Sectional forces/moments of the sub-laminates and of the thick and thin sections (at the locations highlighted in purple and red respectively in Fig. 5b) are computed by integrating the elemental stresses to find the delamination ERRs.
The results are then compared to those obtained from all the other approaches discussed above.
The critical mixed-mode ERR for delamination, G c is calculated according to the BK criterion [28] in Abaqus: where G Ic and G IIc respectively are the material mode I and mode II fracture toughness values and η is an empirical exponent.For T300/ 914C, Ref. [28] reports G Ic = 0.170 N/mm, G IIc = 0.494 N/mm and η = 1.62.The other relevant properties of the cohesive elements used in the models are given in Table 1.For details on the cohesive element formulation available in Abaqus/Standard, the reader is referred to [29].
Only uniaxial tension loading is considered in this parametric study.In a linear elastic analysis, the magnitude of ERR is proportional to the square of the applied load.To gauge the accuracy of different local   beam models, a failure index, FI, in terms of the load applied, can be defined as a function of ERR: where G = G I + G II corresponds to the mixed-mode ERR extracted from the local beam models using the far-field delamination initiation load P initiation obtained from the "benchmark" Model (i).If the computed FI of the beam model is found equal to one, then delamination prediction agrees with the cohesive zone model.So if FI b 0, the prediction is nonconservative; if FI N 0, the delamination onset load estimation is conservative.The FI may also be used to rank different ply-drop configurations, as will be demonstrated in Section 4.2, since a higher FI implies a higher propensity to delaminate.

Results
The delamination initiation loads, P initiation , which correspond to a kink (for stable propagation) or a significant load drop (for unstable propagation) in the load-displacement curves are obtained from Model (i) and Model (ii) for each of the tapered laminates considered here.The initiation locations are reported in Table 2.It is found that delamination initiates from resin pocket in the thin section in all laminates, except for the case of 4_4_4_7°, i.e. for shallowest tapering angle.In the full cohesive FE models with resin, the interfaces between the resin pocket and the sublaminates failed first.This was very quickly followed by delamination in the thin section, as shown in Fig. 6.Once the interfaces have failed, the resin carries no load and its presence is no longer of importance.Therefore, the initiation loads predicted by the two cohesive FE models, with/without considering the presence of the neat resin in the pocket, do not differ significantly in most cases.This justifies the removal of neat matrix pocket from the local beam model.However, the differences become more significant when the resin pockets are relatively elongated, as in the case of 4_4_4_14°, 4_4_4_7°and 16_8_8_27°laminates.In these cases, a considerable extent of interfacial failure between the resin pocket and the "belt" sub-laminate must occur before a thin section delamination can initiate.
For the case of 4_4_4_7°laminate, a "double" delamination is seen to occur in the thick section in the cohesive FE model with resin (Fig. 7).Although the delaminated lengths are different, both delaminations initiate at the same time.On the other hand, in the cohesive FE model without resin, thin section delamination occurs instead of thick section "double" delamination and this occurs at a slightly lower load.The assumption to ignore the resin in the analysis should therefore be treated with caution when the resin pocket is long and shallow as it can affect the resulting failure mode and delamination size/locations.
The P initiation obtained from the benchmark cohesive Model (i) was applied to the corresponding global FE model (Model (iii)).The sectional forces and bending moment extracted at the ply drop-off location were then applied to the three local FE models (Euler-Bernoulli beam -Model (iv), Timoshenko beam -Model (v) and plane stress continuum -Model (vi)).
The values of ERR components at both thin and thick sections for the tapered laminates are detailed in Table A1, Appendix B. The delamination in the thick section occurs predominantly in mode II, while the thin section delamination propagates in mixed-mode regime.Fig. 8 presents the resultant FIs given by different local models.An FI N 1 indicates that the computed ERR is sufficient to initiate delamination, and a further increase from unity implies a more conservative prediction.The FIs at the thin section computed from the plane stress continuum models are consistently closer to unity compared to the beam models.The former correctly identified the delamination onset location for all cases compared to the cohesive FE models, except for 4_4_4_7°laminate, where delaminations propagate in both the thick and thin sections.The generally more conservative results in the thin section are due to the absence of resin in the models and the reasons discussed in Section 2.3.The difference between the results calculated using the plane stress continuum models and the beam models is more evident as the laminate gets thicker.This is clearly because the continuum models allow accounting for stress concentration at the delamination tips and the variation of stresses across thickness.Therefore, the sublaminate forces/bending moments integrated with respect to the neutral axes can be different from the nodal forces/moment directly obtained from the beam solutions, demonstrating a limitation associated with a full beam idealisation.Although the continuum models give better  results, the main drawback is that they require a time-consuming meshing effort and a significant running time.
The Timoshenko beam models perform nearly as well as the continuum models and predict correctly the delamination onset location for all cases, albeit in a more conservative fashion.Fig. 8 clearly shows that more conservative predictions are obtained for thicker laminates, especially for mode I delamination occurring at the thin section.By comparing laminates 4_4_4_27°and 4_4_4_14°which both delaminate in the thin section, it can be observed that the results are more conservative in laminate 4_4_4_14°, due to the larger area of resin pocket being neglected, although the thickness is the same in both cases.On the other hand, the Euler-Bernoulli beam models, which do not account for through-thickness shear deformation, only work well in cases where the sub-laminates are thin.In those cases, the ERR values computed using the Euler-Bernoulli beams do not deviate much from those using the Timoshenko beams.For other cases, the results from the Euler-Bernoulli beam solutions are erroneous, as they fail to predict delamination onset when the sub-laminates are thick.

Applications
This section demonstrates the application of the global-local beam approach as a potential design strategy for studying ply-drop problems.The analysis methodology has been implemented in Abaqus/Standard via Python scripts.Once the global FE model has been set up, the user is required to specify the ply drop-off locations, drop thickness and resin pocket length; the Python scripts automatically set up the local Timoshenko beam models, run the analysis and import the results in the global FE model for post-processing purposes.

Positioning a ply-drop
A hypothetical problem requires one to taper asymmetrically a unidirectional laminate from 32 plies in the thick section to 28 plies in the thin section.It is of interest to find out the best location to place the plydrop when the laminate is subjected to uniaxial tensile loading.The taper angle is held constant at 27°.Three positions of the ply-drops are proposed (Fig. 9), i.e. at 4 plies below the top surface (near to the free surface), 8 plies below the top surface and 12 plies below the top surface (near to the neutral axis).
A global homogenised FE model is built for the problem, as shown in Fig. 9. 20 elements are used in the through-thickness direction of the model.The composite system being considered is T300/914C.The local material axes of the elements in the tapered section are defined such that the fibre direction follows the bottom edge of the elements.An arbitrary 1400 N/mm (it can be any value) tensile load per unit width is applied at the end of the thin section.As the analysis is linear elastic, the delamination load predicted by the beam model can be obtained simply by dividing the applied load by the failure indices (Eq.( 2)).It is worth mentioning that only one global FE model is needed, from which sectional forces/bending moment are extracted and imposed as the boundary conditions in the three local beam models, which correspond to the three different positions of ply-drop.To assess the accuracy of the local Timoshenko beam approach, detailed FE models of those three arrangements of ply-drop with cohesive elements embedded along potential delamination paths are also run.Thin section delamination occurs in all three models.The resulting delamination initiation loads predicted by the cohesive and beam models for the three positions of ply-drop are illustrated in Fig. 10.
It can be observed that the present beam models are generally 2 to 3 times more conservative that the cohesive FE models.Unfortunately, there is no simple relationship or ratio between the two approaches regarding the level of conservatism, as it depends on the relative thickness of both "core" and "belt" sub-laminates.Nonetheless, the local models predict the trend of delamination onset load correctly, allowing selection of the "best" configuration (i.e.laminate 24_4_4_27°out of the three), and, most importantly, they run in very short time, i.e. few seconds.The results indicate that it is advantageous to place the ply-drop nearer to the free surface, which is in agreement with the experimental results presented by Weiss et al. [30].

Symmetric multiple ply-drop specimen
Symmetric multiple ply-drop (MPD) coupons have been manufactured and tested at the University of Bristol [31].They were made of IM7/8552 carbon fibre-reinforced epoxy and were tapered from 32 plies in the thick section to 16 plies in the thin section.The nominal ply thickness was 0.25 mm.The fracture properties for the BK criterion applied to IM7/8552 are G Ic = 0.235 N/mm, G IIc = 0.775 N/ mm and η = 2.6645 [32].Fig. 11 shows a high fidelity FE mesh generated from the scanned image of an actual specimen, illustrating the stacking sequence and the position of ply-drops in the tapered section of the coupon.This comprises 8 ply-drops on each side of the symmetry plane, with ID number 1 to 8 from the thin section.The taper angle of the resin rich area is very shallow, about 7.9°.Results from both experiment and high-fidelity FE analysis (with cohesive elements along the interfaces) show that the interlaminar crack initiated as a "double" delamination at the thick section of the first ply-drop and later propagated towards the thick section, before extending into the thin one.
Two global FE meshes are built to demonstrate the mesh independency of the global-local approach.Both of these are slice models 1 element wide, with applied generalized plane strain boundary conditions.The first model comprises 8 elements (hexahedral elements with reduced integration, C3D8R) in the through-thickness direction, whilst the other one has 16.The outer "mould" lines of the models are kept the same as in the high fidelity FE mesh.As there is same proportion of 0°and ±45°plies in both the thick and thin sections of the laminate, homogenised material properties from classical laminate theory are used for the whole specimen.The average experimental failure load of 65 kN is applied at the end of the thin section.The FI predicted by the model will therefore be defined based on this reference experimental load.As before, FI N 1 implies a conservative prediction.The elemental stresses on the sections where the ply-drops are located are integrated to obtain the sectional forces/bending moments.Due to the symmetry of the problem, only ply-drops in the top half of the models are considered.In the local beam models, the thicknesses of the individual beams are those of the sub-laminates.Both the global and local models employ the same homogenised properties.
The delamination FIs in the thick and thin sections of the 8 ply-drops in the top half of the symmetric MPD specimen are computed by the beam models and shown in Fig. 12.There is no difference in the results obtained from the coarse and fine global FE models, once again proving that the proposed global-local approach is mesh independent.The results are shown to be very conservative for thin section delamination (all have FI N 1, mode I dominated) due to the reasons already discussed in Section 2.3.The first ply-drop has the highest FI.The ply-drops nearer to the mid-plane of the specimen (those with odd number) are more inclined to delaminate than those nearer to the free surface.On the other hand, the model predicts that thick section delamination occurs only in the first ply-drop with an FI = 1.1 (mode II dominated).This delamination prediction agrees with the experimental observations and the results obtained from the high-fidelity FE analysis [31].
Despite the simplification and homogenisation procedures adopted in the analysis, the beam model captures mode II dominated delamination quite accurately.This is indeed the main fracture mode in the thick section of ply-drops under uniaxial tension.This is also observed in the case of 4_4_4_7°laminate in Section 3.1.Even without foreknowledge of the actual failure mode (thick or thin section), one can use the more conservative thin section delamination load prediction for the relative ranking of different stacking sequences.The model correctly predicts which ply-drop delaminates first and also predicts that inner ply drop-offs are more likely to delaminate than the outer ones.

Conclusions
The proposed global-local approach to predict delaminations emanating from the resin pocket associated with a ply-drop requires a global FE model to compute sectional forces/bending moments for the local model.The global force/moments are mesh independent, thus allowing using of a very coarse global FE mesh.On the other hand, the local Timoshenko beam model takes into account the geometry of the resin pocket and it is able to capture the fracture onset at both the thick and thin sections with minimal computational resources.Comparative studies with cohesive FE analysis have demonstrated that the local beam approach is able to correctly identify both the failure location and the relative ranking of laminates.The approach proposed also consistently provides conservative predictions of the ERR, especially for a mode I dominated fracture.This is because the transverse forces calculated by the local beam model tend to give a high mode I component at the thin section for thick laminates (see Table A1).
However, caution should also be exercised when dealing with very shallow taper angles, as the removal of the resin pocket from the analysis can affect both the failure mode and location.Although the absolute accuracy of the global-local approach (i.e. the prediction of the actual failure load) is limited due to the idealisation of the problem, the proposed method is more than adequate for comparative studies between different ply-drop configurations.Overall, the global-local approach proposed here is easy to implement as an automated analysis tool for quick preliminary design studies, when there is a need to compare many alternative possibilities in terms of the through-thickness location of ply terminations, drop thickness and the overall order of ply terminations.Future work will address the implementation of the methodology developed in this paper into an optimisation framework for the preliminary design of tapered composite laminates.The overall ERR due to shear in the thick section can be written as [27]: In order to split it into pure mode I and mode II contributions, the load partitioning scheme in Fig. A1 is adopted.Since the 'double' delamination configuration creates three sub-laminates, pure mode II fracture occurs when the shear strain γ in the three arms is the same, i.e.: The same homogenised interlaminar shear modulus G xz is assumed for all the sub-laminates.Hence, the identities (Eq.(A2)) lead to: where ξ is the ratio of the corresponding sub-laminate thickness t to the thick section laminate thickness.The subscripts A, B and C refer to the 'belt', 'drop' and 'core' sub-laminates respectively.Mode I opening can occur independently in the top and bottom delaminations of the thick section.This is accounted for in Fig. A1 by introducing the mode I component V IA and V IC .Considering the overall shear force resultant for the "belt" and "core" sub-laminates and the free end assumption for the "drop" sub-laminate, one gets: Eq. (A4) leads to three linear equations in 3 unknowns, i.e.V IA , V IC and V II .The two mode I components V IA and V IC are not independent from the mode II component, V II .Solving for the three unknowns, one obtains: Moreover, the shear force equilibrium condition dictates that: Substituting the first and third of Eq. (A4) into Eq.(A6) leads to: From the second of Eq. (A4) and exploiting the identity ξ A + ξ B + ξ C =1, Eq. (A7) can be simplified to: Finally, by substituting Eqs.(A4) and (A8) into Eq.(A1) and after some further simplifications, the expression of the ERR due to shear at the thick section becomes: It is obvious that the mode I and mode II components in Eq. (A9) are decoupled.The ERR components are therefore: Note that Eq. (A10) holds only when V IA N 0 and V IC N 0, i.e. when the shear forces are acting in the crack-opening direction in the "belt" and "core" sub-laminates, respectively.Moreover, G II , thick (V) vanishes when t drop → 0.

A.2. ERR contribution from axial force
Axial forces applied with respect to the sub-laminate neutral plane produce only a pure mode II ERR.At the thin section, the mode II ERR is given by [23,27]: At the thick section, the mode II ERR is: Note that the additional factor of 1/2 in Eq. (A13) comes from the 'double' delamination (i.e.double fracture surfaces) assumption made for the thick section delamination; t is the thickness, with the subscript denoting the corresponding laminate/sub-laminate.

A.3. ERR contribution from bending moment
The total ERR due to the bending moment can be partitioned into mode I and mode II components as follows.For delamination at the thin section, these are [23,27]: where ψ R =(t core /t belt ) 3 .On the other hand, the mode I and mode II ERR due to the bending moment at the thick section are [23]: where ψ A =(t drop /t belt ) 3 , ψ C =(t core /t drop ) 3 and: A.4. ERR contribution from shear (transverse) force (for thin section) It is assumed that the shear stress has a parabolic distribution in the laminate/sub-laminates.At the thin section, the shear forces acting in the 'crack-opening' direction give rise to mode I only and the associated ERR contribution is [27]: where and ξ = t belt /t thin .V I b 0 implies that the shear forces are acting in the crack-opening direction and so Eq.(A20) gives a non-zero result.

Fig. 1 .
Fig. 1.(a) A severely tapered carbon/epoxy composite laminate, (b) the corresponding high-fidelity FE half-model generated from the scanned image of an actual specimen [25].

Fig. 2 .
Fig. 2. (a) Tapered composite section with idealised ply-drop configuration; (b) global FE model with coarse mesh and homogenised material properties; (c) local equivalent beam assembly.

Fig. 3 .
Fig. 3.A two-dimensional ply-drop section reduced to a simpler assembly of Timoshenko beams (not to scale).

Fig. 5 .
Fig. 5. (a) Model (i): An FE continuum (plane stress) model with cohesive elements placed along the interfaces between the sub-laminates and the resin pocket.(b) Model (vi): An FE continuum (plane stress) model with a hollow resin pocket and without cohesive elements.

Fig. 6 .
Fig. 6.Damage variable, SDEG, of cohesive elements showing: (Left) Interfaces between the resin pocket and the sub-laminates fail first, followed by delamination along the thin section in a full cohesive FE model.(Right) Thin section delamination happens in the model without the resin pocket.

Fig. 7 .
Fig. 7. Delamination propagation (indicated by arrows) of the cohesive FE models of 4_4_4_7°tapered laminate, (left) with and (right) without the resin pocket modelled.

Fig. 8 .Fig. 9 .
Fig. 8. FI at the thick (left) and thin (right) sections calculated by different local FE models for various tapered laminates under uniaxial tension.FI N 1 (above the red lines) implies delamination failure.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 10 .Fig. 11 .
Fig. 10.The delamination onset loads predicted by the beam models and compared against results from the cohesive FE models.

Fig. 12 .
Fig. 12. Failure indices for delamination in the (left) thick and (right) thin sections of the 8 ply-drops in the top half of the MPD specimen as computed by the beam models.

Table 2
Delamination initiation loads and locations predicted by the cohesive FE models with and without modelling the resin under uniaxial tension.

Table A1
ERR components and failure indices (FI) at the thick and thin sections calculated by different local FE models for various tapered laminates under uniaxial tension.FI N 1 implies delamination failure (highlighted in bold).Global Model (iii) + local Model (iv).⁎⁎ Global Model (iii) + local Model (v).⁎⁎⁎ Model (vi).