ELSEVIER

Contents lists available at ScienceDirect

Materials and Design

journal homepage: www.elsevier.com/locate/matdes

A simplified layered beam approach for predicting ply drop delamination in thick composite laminates

Khong Wui Gan a'*, Giuliano Allegrib, Stephen R. Hallettc

a Faculty of Engineering and the Environment, University of Southampton Malaysia Campus (USMC), KotaIlmu Educity @ Iskandar, 79200 Nusajaya, Johor, Malaysia b Department of Aeronautics, Imperial College London, South Kensington Campus, London SW72AZ, UK c Advanced Composites Centre for Innovation and Science, University of Bristol, University Walk, Bristol BS8 1TR, UK

CrossMark

HIGHLIGHTS

GRAPHICAL ABSTRACT

1 A novel global-local approach is proposed to predict ply-drop delamination in thick tapered composites.

1 The ply-drop is represented locally by an assembly of Timoshenko beams with loading obtained from a global FE model.

1 The relative accuracy given by the approach can be used to efficiently rank delamination hazards of multiple ply-drops.

1 Its applicability to thick tapered composite specimens is demonstrated.

ARTICLE INFO

Article history:

Received 28 April 2016

Received in revised form 23 June 2016

Accepted 24 June 2016

Available online 7 July 2016

Keywords: Composite laminate Ply drop-off Delamination Numerical analysis

ABSTRACT

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 finite 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://

creativecommons.org/licenses/by/4.0/).

1. Introduction

1.1. Strength prediction for tapered laminates

The reasons why through-thickness tapering is commonly adopted in lightweight structures are twofold: 1) achieving a pre-

* Corresponding author. E-mail address: K.W.Gan@soton.ac.uk (K.W. Gan).

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. dropping-off) 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

http://dx.doi.org/10.1016/j.matdes.2016.06.105

0264-1275/© 2016 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).

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-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.

Therefore, several analytical and numerical techniques for the predicting the strength of tapered laminates have been proposed in the literature. These approaches can be broadly classified in [1]: 1) strength-based [3,6-14]; 2) fracture mechanics-based [15-24]. Strength approaches rely upon calculating the stress distribution at a ply drop-off location and then apply interactive stress criteria (e.g. quadratic) for de-lamination onset [3,10]. Analytical strength models are primarily based on shear lag approximations [8,10,14], while numerical techniques rely either on displacement-based [6] or hybrid finite element analysis [9]. Owing to the singular nature of tractions at a ply termination, point stresses or average stresses should be considered [6]. Fracture mechanics methods involve estimating the energy release rate (ERR) associated with delamination emanating from a ply drop-off and then employ the Griffith criterion for propagation. This is done either in the context of analytical beam/plate modelling [18,21-23] or via finite element analysis (FEA). The latter can be coupled either with the virtual crack closure technique [15-17,19,24] or cohesive zone modelling [25]. Cohesive zone FEA applied to ply-drop analysis has the advantage of being able to predict both delamination onset and growth, unifying strength and energy-based integrity assessments [25].

1.2. 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.

13. 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 shear-undeformable 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

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].

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.

2. Modelling methodology

2.1. 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 de-formable) 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 Timo-shenko 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 mesh-dependent. 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.

2.2. 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, Sa [27]. The overall ERR, or G, can be partitioned into pure mode I and mode II components, i.e. G = Gi + GH. Fracture occurs when G = Gc, where Gc is the critical mixed mode ERR. The ERR expressions associated with axial

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.

Thick section Thin section

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

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.

c. Beam elements are used in the local models

Beam elements cannot capture the stress gradients in the through-thickness 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.

2.3. 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 de-laminations 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].

3. Methodology validation

3.1. Virtual testing of single ply drop-off coupons

The objective of this section is to evaluate the accuracy of using one-dimensional 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_0. The first three letters a, b and c, refer to the number of unidirectional plies in the "core", "drop" and "belt" sub-laminates, respectively; 9 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:

Fig. 4. Global-local modelling framework for ply-drop delamination in tapered laminates.

Table 1

Mechanical properties for T300/914C unidirectional tape and 914C resin.

T300/914C ply properties (1 = fibre direction) [24]

E1 (MPa) E2 = E3 (MPa) Gi2 = Gi3 (MPa) G23 (MPa) v„ = V13 V23

135,200 9000 5000 3300 0.3 0.5

914C resin properties [24]

E (MPa) G (MPa) V

4000 1481 0.35

T300/914C inter-ply cohesive properties (I = mode I) [28]

GIc GIIc = GIIIc (N/mm) (N/mm) °Imax °IImax — °IIImax KI KII = KIII (MPa) (MPa) (N/mm3) (N/mm3)

0.17 0.494 75 80 4.0e6 1.5e6

i) A detailed two-dimensional plane stress FE model (Fig. 5a) with cohesive elements (COH2D4) placed at the interfaces between the sub-laminates and the resin pocket. The meshes are regular and comprise 2 elements per ply thickness. This FE model serves as a numerical benchmark, from which reference results such as the delamination initiation location and load are obtained for validation purposes. The delamination initiation load is then applied as the far-field boundary condition in the homogenised global FE model described in Model (iii);

ii) A detailed two-dimensional plane stress FE model similar to Model (i), but with a void instead of a resin pocket. Cohesive elements are embedded along interfaces between the sub-laminates as in Model (i). The purpose of this model is to investigate the effect of the hollow resin pocket assumption made in the analysis;

iii) A coarse homogenised global FE model, meshed as shown in Fig. 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 (<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;

Table 2

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

Tapered Cohesive FE models WITH resin Cohesive FE models WITHOUT

Laminate (Model (i)) resin (Model (ii))

Pirntintioit Delamination Pinitiation Delamination

(N/mm) initiation location (N/mm) initiation location

2 4 2 27° 257 Thin section 237 Thin section

4_4_4_27° 420 Thin section 434 Thin section

4_4_4_14° 761 Thin section 568 Thin section

4 4 4 7° 938 Thick section 859 Thin section

8_4_8_27° 809 Thin section 838 Thin section

16_8_8_27° 1044 Thin section 844 Thin section

16_4_8_27° 1206 Thin section 1170 Thin section

20_4_8_27° 1415 Thin section 1354 Thin section

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, Gc is calculated according to the BK criterion [28] in Abaqus:

Gc = Gic + (Glk - GiJbr-^) (1)

where Glc and Gllc respectively are the material mode I and mode II fracture toughness values and n is an empirical exponent. For T300/ 914C, Ref. [28] reports GIc = 0.170 N/mm, G„c = 0.494 N/mm and n = 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

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. 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.

beam models, a failure index, FI, in terms of the load applied, can be defined as a function of ERR:

where G = Gi + Gn corresponds to the mixed-mode ERR extracted from the local beam models using the far-field delamination initiation load Pinitiation 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 < 0, the prediction is non-conservative; if FI > 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.

3.2. Results

The delamination initiation loads, Pinitiation, 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 Pinitiation 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 delamina-tion 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 > 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

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. FI at the thick (left) and thin (right) sections calculated by different local FE models for various tapered laminates under uniaxial tension. FI > 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.)

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.

4. 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 Ti-moshenko beam models, run the analysis and import the results in the global FE model for post-processing purposes.

4.1. 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 ply-drop 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

Fig. 9. A global FE model and three possible locations to place the ply-drops.

1400 ■ _

| 1200

I™ 1 800 | 600

I™ II II

24_4_4_37* 20.4 3.27'

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

nearer to the free surface, which is in agreement with the experimental results presented by Weiss et al. [30].

4.2. 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 GIc = 0.235 N/mm, GIIc = 0.775 N/ mm and n = 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 > 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.

Fig. 11. Stacking sequence of the symmetric MPD specimen in the tapered section. The high-fidelity FE mesh was generated from the scanned image of an actual specimen [31 ].

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 > 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 delamina-tion 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 delamina-tion 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.

5. 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 Timo-shenko 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.

Acknowledgements

The presented work was funded by the EPSRC through a University ofBristol Impact Acceleration Award 'A Tool for the Multi-Scale Damage Tolerant Design of Tapered Composite Laminates' (Grant no. EP/ K503824/1). The support of Rolls-Royce plc through the Composites University Technology Centre (UTC) at the University ofBristol is also gratefully acknowledged. Supporting data may be requested from Prof. S.R. Hallett. Access to supporting data may be granted, subject to

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.

consent being requested and granted from the original project participants.

Appendix AA.1. Derivation of ERR mode partition due to shear forces in thick section

The overall ERR due to shear in the thick section can be written as [27]:

Gt v> = 3 f I _ Vthick |

th'ck 10Gxz ytbelt fcore ("thick J

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' delam-ination configuration creates three sub-laminates, pure mode II fracture occurs when the shear strain y in the three arms is the same, i.e.:

VII 4>V dropII 4>va

Gxzfbelt Gxztd

'xz*- drop

Gxzt co

The same homogenised interlaminar shear modulus G^ is assumed for all the sub-laminates. Hence, the identities (Eq. (A2)) lead to:

— tdrop=tbelt — ^core — tcore/tbelt — §C/^A

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 VIA and VIC. 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:

Vtl = -Via + V„, 0 = Via-Vic + , Vbl = Vic + (A4)

Eq. (A4) leads to three linear equations in 3 unknowns, i.e. VIA, VIC and VII. The two mode I components VIA and VIC are not independent from the mode II component, VII. Solving for the three unknowns, one obtains:

Via = -(£,b + §cWtl + £,aVbl, Vic = (£,a + £,bWbl-£,cVtl, Vu

= §a(Vtl + Vbl) (A5)

Moreover, the shear force equilibrium condition dictates that:

V thick = Vtl + VBL (A6)

Substituting the first and third of Eq. (A4) into Eq. (A6) leads to:

V thick — _Via + Vic + (1 + ^W«

Fig. A1. Mode partitioning scheme for the thick section delamination due to shear forces.

From the second of Eq. (A4) and exploiting the identity ^ + + §C =1, Eq. (A7) can be simplified to:

V thick = vii /£,a

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:

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]:

g(m) MRMTR-MBR\ 2(J2 J2!

'■thin 2 E»A 1 + ) !vt3e,t + t3oj

G(V) _

thick _

iA | V2C | tdrop V11 tbelt tcore

G(M) _LJ_

11 2 Ex;

Mtr + Mba Y 12 12^R\-

1 + J Ui

12 Mt2h

It is obvious that the mode I and mode II components in Eq. (A9) are decoupled. The ERR components are therefore:

(Vk+VL)

tbelt tcore

where ipR = ( tcore/tbelt)3. On the other hand, the mode I and mode II ERR due to the bending moment at the thick section are [23]:

1 1 Î12MI 12m2L

/.thick

4 M t3elt t

u//,thick

10Gxzt2elt

G(M) _ 1 J_ u//,thick _ 4 E

,/12(1 + 2^), 12^"^ 12M2hi

tthick

Note that Eq. (A10) holds only when VIA > 0 and VIC > 0, i.e. when the shear forces are acting in the crack-opening direction in the "belt" and "core" sub-laminates, respectively. Moreover, GiVhick vanishes when

tdrop ^ 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]:

where ^ = (tdrop/tbeit)3, <Pc =(tcore/tdrop)3 and: M Mtl^a(1 + <ïc)-Mbl m MtlMc-Mbl(1 + M

MIA =-, M , , N , -,-, MIC = -

^a(1 + >Pc ) + 1 Mtl + Mbl

<M1 + ^c ) + 1

^a(1 + <ïc ) + 1

A.4. ERR contribution from shear (transverse) force (for thin section)

//.thin _ 2 E

(p2 p2 F2 \ I 1 TR | 1 BR — 1 thin |

V^belt tcore tthin J

At the thick section, the mode II ERR is:

1 1 / PTL ! PBL — Pthick | 4 Exx ^belt tcore tthick J

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]:

G/.thin " 5 Gxz

tbelt tcore

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.

V/ = £Vbr-(1-©VTR

and § = tbelt/tthin. VI < 0 implies that the shear forces are acting in the crack-opening direction and so Eq. (A20) gives a non-zero result.

Appendix B

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 > 1 implies delam-ination failure (highlighted in bold).

Tapered laminate FE model Thick section Thin section

G/ G// n G/ G// n

2 4 2 27° Euler-Bernoulli* 0.0516 0.1332 0.7161 0.3595 0.4542 1.6579

Timoshenko** 0.0234 0.1085 0.5701 0.2801 0.5704 1.5827

Plane stress continuum*** 0.0352 0.1025 0.6096 0.1735 0.4445 1.3104

4_4_4_27° Euler-Bernoulli 0.0154 0.1554 0.6178 0.4873 0.2381 1.8016

Timoshenko 0.0016 0.0877 0.4292 0.4415 0.4787 1.8050

Plane stress continuum 0.0388 0.0882 0.6031 0.1915 0.3589 1.2876

4_4_4_14° Euler-Bernoulli 0.0242 0.5032 1.0593 0.6294 0.6749 2.1526

Timoshenko 0.0123 0.3961 0.9244 0.6008 0.8707 2.1841

Plane stress continuum 0.0669 0.3314 0.9851 0.3543 0.6492 1.7439

(continued on next page)

/.thick

Table A1 (continued)

Tapered laminate FE model Thick section Thin section

Gi Gii FI Gi Gii FI

4_4_4_7° Euler-Bernoulli 0.0238 0.7219 1.2502 0.2929 0.8564 1.7600

Timoshenko 0.0201 0.6792 1.2088 0.2883 0.9211 1.7881

Plane stress continuum 0.0555 0.5055 1.1248 0.2081 0.4496 1.3810

8_4_8_27° Euler-Bernoulli 0.0537 0.1475 0.7417 0.0448 0.1298 0.6868

Timoshenko 0.0014 0.1340 0.5264 0.8403 0.6205 2.4121

Plane stress continuum 0.0483 0.0971 0.6554 0.2261 0.2940 1.3197

16_8_8_27° Euler-Bernoulli 0.0232 0.1268 0.6001 1.7151 0.1060 3.2395

Timoshenko 0.0169 0.1373 0.5933 0.9220 0.3775 2.4646

Plane stress continuum 0.0828 0.1210 0.8119 0.3937 0.2862 1.6494

16_4_8_27° Euler-Bernoulli 0.1081 0.0964 0.8770 0.0846 0.0385 0.7489

Timoshenko 0.0073 0.1777 0.6254 1.0622 0.3984 2.6392

Plane stress continuum 0.0527 0.1033 0.6811 0.2901 0.1583 1.3958

20_4_8_27° Euler-Bernoulli 0.1182 0.0942 0.9093 0.0892 0.0255 0.7601

Timoshenko 0.0063 0.1914 0.6437 1.0974 0.3461 2.6714

Plane stress continuum 0.0504 0.1108 0.6823 0.3046 0.1263 1.4172

* Global Model (iii) + local Model (iv). ** Global Model (iii) + local Model (v). *** Model (vi).

References

[1] K. He, S.V. Hoa, R. Ganesan, The study of tapered laminated composite structures: a review, Compos. Sci. Technol. 60 (2000) 2643-2657.

[2] S.V. Hoa, B.L. Du, T. Vu-Khanh, Interlaminar stresses in tapered laminates, Polym. Compos. 9 (5) (1988) 337-344.

[3] A. Mukherjee, B. Varughese, Design guidelines for ply drop-off in laminated composite structures, Compos. Part B 32 (2001) 153-164.

[4] D.B. Adams, L.T. Watson, Z. Gurdal, C.M. Anderson-Cook, Genetic algorithm optimization and blending of composite laminates by locally reducing laminate thickness, Adv. Eng. Softw. 35 (2004) 35-43.

[5] F.-X Irisarri, A Lasseigne, F.-H. Leroy, R. Le Riche, Optimal design of laminated composite structures with ply drops using stacking sequence tables, Compos. Struct. 107 (2014) 559-569.

[6] J.C. Fish, S.W. Lee, Delamination of tapered composite structures, Eng. Fract Mech. 34 (1) (1989)43-54.

[7] E.A Armanios, L. Parnas, Delamination Analysis of Tapered Laminated Composites under Tensile Loading, in: T.K O'Brien (Ed.), Composite Materials: Fatigue and Fracture, ASTM STP1110, 3rd vol., 1991, pp. 340-358.

[8] AJ. Vizzini, Strength of laminated composites with internal discontinuities parallel to the applied load, AIAAJ. 30 (6) (1992) 1515-1520.

[9] P.N. Harrison, E.R. Johnson, A mixed variational formulation for interlaminar stresses in thickness-tapered composite laminates, Int. J. Solids Struct. 35 (16) (1996) 2377-2399.

[10] A.J. Vizzini, Shear-lag analysis about an internally-dropped ply, J. Reinf. Plast. Compos. 16 (1) (1997) 73-85.

[11] F. Mortensen, O.T. Thomsen, A simple approach for the analysis of embedded ply drops in composite and sandwich laminates, Compos. Sci. Technol. 59 (1999) 1213-1226.

[12] B.R. Vidyashankar, A.V. Krishna Murty, Analysis of laminates with ply drops, Compos. Sci. Technol. 61 (2001) 749-758.

[13] S.C. Her, Stress analysis of ply drop-off in composite structures, Compos. Struct. 57 (2002) 235-244.

[14] K. He, R. Ganesan, S.V. Hoa, Modified shear-lag model for analysis of a composite laminate with Drop-off plies, Compos. Sci. Technol. 63 (10) (2003) 1453-1462.

[15] S.I. Salpekar, I.S. Raju, T.K. O'Brien, Strain-energy-release rate analysis of delamination in a tapered laminate subjected to tension load, J. Compos. Mater. 25 (1991) 118-141.

[16] G.B. Murri, S.A. Salpekar, T.K. O'Brien, Fatigue Delamination Onset Prediction in Unidirectioal Tapered Laminates, in: T.K. O'Brien (Ed.), Composite materials: Fatigue and Fracture, ASTM STP 1110, 3rd. vol., 1991, pp. 312-339.

[17] G.B. Murri, T.K O'Brien, C.Q. Rousseau, Fatigue life methodology for tapered composite flexbeam laminates, J. Am. Helicopter Soc. 43 (2) (1998) 146-155.

[18] W. Cui, M.R. Wisnom, M.I. Jones, An experimental and analytical study of delamination of unidirectional specimens with cut central plies, J. Reinf. Plast. Compos. 13 (1994) 722-739.

[19] W. Cui, M.R. Wisnom, M.I. Jones, Effect of through thickness tensile and compressive stresses on delamination propagation fracture energy, J. Compos. Technol. Res. 16 (4) (1994) 329-335.

[20] W. Cui, M.R. Wisnom, M.I. Jones, Effect of step spacing on delamination of tapered laminates, Compos. Sci. Technol. 52 (1994) 39-46.

[21] W. Cui, M.R. Wisnom, M.I. Jones, New model to predict static strength of tapered laminates, Composites 26 (2) (1995) 141-146.

[22] Z. Petrossian, M.R. Wisnom, Parametric study of delamination in composites with discontinuous plies using an analytical solution based on fracture mechanics, Compos. Part A 29A (1998) 403-414.

[23] G. Allegri, M.R. Wisnom, S.R. Hallett, A simplified approach to the damage tolerance design of asymmetric tapered laminates. Part I: methodology development, Compos. Part A 41 (10) (2010) 1388-1394.

[24] G. Allegri, M.R. Wisnom, S.R. Hallett, A simplified approach to the damage tolerance design of asymmetric tapered laminates. Part II: methodology validation, Compos. Part A 41 (10) (2010) 1395-1402.

[25] S.R. Hallett, J.K. Lander, M.I. Jones, L.F. Kawashita, M.R. Wisnom, Testing and Modelling of a Severely Tapered Composite Specimen, 5th International Conference on Composites Testing for Model Identification (CompTest), Lausanne, 2011.

[26] A. Mukherjee, B. Varughese, A ply drop-off element for inclusion of drop-off in the global analysis of layered composite structures, Comput. Struct. 54 (5) (1995) 865-870.

[27] J.G. Williams, On the calculation of energy release rates for cracked laminates, Int J. Fract. 36 (2) (1988) 101-119.

[28] M. König, R. Krüger, K. Kussmaul, M.V. Alberti, M. Gadke, Characterizing Static and Fatigue Interlaminar Fracture Behavior of a First Generation Graphite/Epoxy Composite, in: S. Hopper (Ed.), Composite materials: Testing and Design, AsTm STP 12421997, vol. 13, 1997, pp. 60-81.

[29] Abaqus 6.13 Analysis User's Guide. Section 32.5.

[30] A. Weiss, L. Michel, S. Mahdi, E. Cros,J.J. Barrau, Influence of Ply-Drop Position in Thickness Direction on Static and Fatigue Loading Behaviour of Carbon Fibre Epoxy Laminates, Proceedings of ICCM-17, Edinburgh, 2009.

[31] L.F. Kawashita, M. Jones, S. Giannis, S.R. Hallett, M.R. Wisnom, High fidelity modelling of tapered laminates with internal ply terminations, 18th International Conference on Composite Materials (ICCM), Jeju Island, 2011.

[32] G. Allegri, M.R. Wisnom, S.R. Hallett, A new semi-empirical law for variable stressratio and mixed-mode fatigue delamination growth, Compos. Part A 48 (2013) 192-200.