Modified shear-lag model for analysis of a composite laminate with drop-off plies

Modified shear-lag model for analysis of a composite laminate with drop-off plies

Composites Science and Technology 63 (2003) 1453–1462 www.elsevier.com/locate/compscitech Modified shear-lag model for analysis of a composite laminat...

430KB Sizes 11 Downloads 94 Views

Composites Science and Technology 63 (2003) 1453–1462 www.elsevier.com/locate/compscitech

Modified shear-lag model for analysis of a composite laminate with drop-off plies K. He*, R. Ganesan, S.V. Hoa Department of Mechanical Engineering, Concordia Center for Composites, Concordia University, Montreal, Quebec, Canada Received 31 October 2002; received in revised form 10 March 2003; accepted 11 March 2003

Abstract In this paper, a modified shear-lag model was developed and implemented to perform interlaminar stress analysis of the unidirectional laminate with drop-off plies in tension. The essential assumptions for the established model are that both plies and resin layers are treated to act as carriers of tensile stress and also to act as stress-transfer media. The model, therefore, overcomes the limitations exhibited in the conventional shear-lag and is capable of addressing delaminations in the laminated composite due to material and geometric discontinuities at the region close to the drop-off plies. The validity and efficiency of the model have been illustrated by comparison of results obtained with the current model and ones from the existing literature. # 2003 Elsevier Science Ltd. All rights reserved. Keywords: Modified shear-lag model; Drop-off plies; Interlaminar stresses; C. Delamination

1. Introduction The laminated composite formed by dropping off some of the plies within the structure, also named thicknesstapered laminate, has received wide applications in aircraft industry because of its elastic tailoring capacities and great potential for creating more significant weight savings than regular laminated composite components. A comprehensive review on the works related to the tapered laminated composite was presented in Refs. [1,2]. To better understand failure mechanisms induced by drop-off plies in the tapered construction, stress state in the vicinity of ply drop-offs should be examined first. Many approaches were considered in this regard for prediction of interlaminar stress in tapered composite structures, which include finite element methods [3–8], extentional and shear spring model [9], mixed formulation [10], and shear-lag analysis [11]. Throughout these efforts one objective has been to understand the load transfer mechanism caused due to material and geometric discontinuities at the ply drop locations. * Corresponding author at present address: 1700 Kirts Boulevard, No. 201, Troy, MI 48084, USA. Tel.: +1-248-838-5236; fax: +1-248838-5300. E-mail address: [email protected] (K. He).

Armanios and Parnas [9] developed a model composed of extensional and shear springs and shear springs and predicted the interlaminar stresses based on the tendency of the plies to align themselves to the applied load. Although the comparison between their model and a finite element solution is good for interlaminar shear, it fails to capture the tensile nature of the interlaminar normal stress at the ply drops. Harrison and Johnson [10] developed a mixed formulation that correlates well with their finite element model of the section. However, their model lacks a physical description of the load transfer mechanism. On the other hand, extensive work has been done to stress analysis of the tapered laminate using finite element methods. The finite elements employed by previous authors to model the tapered composite construction are in general formulated based on displacement assumption, generalized plane deformation theory, or hybrid theory. Full three-dimensional, quasi-three-dimensional or two-dimensional elements have been considered and used to predict the stresses in the drop-off region of the laminate in terms of analysis requirements. Elements are refined at ply dropped-off region and the intersection between thin section and taper as to capture the realistic influence of the ply discontinuities on interlaminar stresses. The smallest element size is

0266-3538/03/$ - see front matter # 2003 Elsevier Science Ltd. All rights reserved. doi:10.1016/S0266-3538(03)00166-0

1454

K. He et al. / Composites Science and Technology 63 (2003) 1453–1462

equal to one quarter of the ply thickness. The density of element (number of elements per unit area) in the ply drop location is about 43 times of that of other normal region. So, computational effort is prohibitive, particularly for three-dimensional analysis. From the stand of improving computational efficiency, the modified shearmodel, which is extended from the ordinary shear-lag theory, was developed and implemented to solve the tapered laminated composite problem. In the idealized shear-lag theory, the elastic fibers carry only simple tension and the elastic matrix in-between the fiber deforms only in simple shear. Vizzini [11] applied the ordinary shear model to load transfer mechanism analysis for tapered composite structures in 1997. The fundamental basis of his shearlag model to be developed is that of three elastic layers separated by two shear layers. It was shown that this ordinary shear-lag model could capture a majority of the load transfer mechanism for the composite structures only when Young’s modulus of resin layer is significantly lower than the Young’s modulus of unidirectional tape. Another limitation in the conventional shear-lag model is its inability to address the matrix cracking problems that are often seen in composite laminates. This constrains its wide applications in laminated composites where resin cracks appear to be one of most premature failure event that could lead to eventual delamination failure of the structure. The new model, however, overcoming the limitations of the model with ordinary shear-lag and capable of being used to any type of composite laminates, was developed based on the assumptions that both plies and resin layers are treated to act as carriers of tensile stress and also to act as stress-transfer media. It is therefore applicable to more complicated analysis of laminates, such as delamination behavior and any modulus ratios between resin and ply layers. Meanwhile, it was formulated based on a conceptual model with comparatively simple formations and could provide with less computing efforts than that required in finite element methods more accurate stress responses of tapered laminates than the ordinary shear-lag model. In general, simplified mechanical models for the interlaminar analysis of tapered laminates, as indicated in previously published works, provide more physical insight than FEM approaches, and reasonable results comparable to that predicted with FEM methods are achievable based upon physically appropriate assumptions.

2.1. Laminate without drop-off layers A schematic drawing of model with variable layer thickness is shown in Fig. 1. For equilibrium of the force in each layer in x direction, as indicated in Fig. 2, there exists: dP1 þ 2=1 ¼ 0 dx

ðBottom layer;

dPi þ iþ1=i  i=i1 ¼ 0 dx ðIntermediate layers; dPN  N=N1 ¼ 0 dx

ð1aÞ

ð1bÞ i ¼ 2; 3; . . . N  1Þ

ðTop layer;

i ¼ NÞ

ð1cÞ

where, Pi,  i/j are, respectively, the axial force per unit thickness of the ith layer and interfacial shear stress. (2i1)th and (2i)th layers are corresponding to plies and resin layers, respectively.

Fig. 1. Schematic drawing of model with variable layer thicknesses.

2. General formulation Based on the aforementioned assumptions, two types of laminates were considered in formulation of modified shear-lag model for analysis of laminated composites with drop-off plies.

i ¼ 1Þ

Fig. 2. Diagram of free body.

K. He et al. / Composites Science and Technology 63 (2003) 1453–1462

Under the conditions of linear elasticity, the following equations hold, considering the assumptions above: Pi ¼ Ei ti i=j ¼

dui dx

ði ¼ 1; 2;    NÞ

   2Gi Gj  ui  uj ¼ Ki=j ui  uj ti Gj þ tj Gi

ð2Þ ð3Þ

where, Ei—Elastic modulus, E2i1 for plies, and E2i for resin layers; Gi—Shear modulus, G2i1 for plies, and G2i for resin layers; ui—Displacement of the ith layer; ti— Thickness of layer as ti=ai+bi x/L (i=1, 2, 3, . . . N); and Ki/j—Interfacial modulus, derived using the definition of shear deformation as follows (refer to Fig. 3 for the derivation).     Gi ui  ui=j Gj ui=j  uj   ¼ ð4Þ i=j ¼ ðti =2Þ tj =2 Therefore, the final differential system can be obtained by the substitutions of the latter equations into Eq. (1).   dt1 du1 d 2 u1 E1 þ t1 ð5aÞ þ K2=1 ðu2  u1 Þ ¼ 0 dx dx dx2 

2



dti dui d ui Ei þti 2 þKiþ1=i ðuiþ1 ui ÞKi=i1 ðui ui1 Þ¼0 dx dx dx i ¼ 2; 3; 4;    N  1

ð5bÞ

  dtN duN d 2 uN þ tN EN  KN=N1 ðuN  uN1 Þ ¼ 0 ð5cÞ dx dx dx2 For the differential equation system in Eq. (5), it is impossible to obtain a closed form solution, so we can solve it in a numerical way. Suppose that the displacements in each layer are of the following form: 1  x n X ui ¼ Cni i ¼ 1; 2; 3;    N: ð6Þ L n¼0

for each layer. By setting the coefficients of like terms on both sides of the equations equal, a solution can be obtained with two initial values for each layer determined by boundary conditions. Derivation of the coefficients was given in Ref. [2]. Eq. (5) can be further simplified to model laminated composites with five layers/ five sublaminates, as implemented in this work. 2.2. Laminate with drop-off layers The simplest tapered laminate consists of a two-zone area close to the drop off location, as shown in Fig. 4. The first zone consists of four constant thickness plies and resin layer and one variable thickness resin layer. The top ply is inclined. In zone two, three resin layers are located between the top and bottom plies. Similarly to the derivation of equations for the regular plate without drop-off layer, the governing equilibrium equations for general laminate with total number of N layers and with drop-off layer at the top the laminate are given in the following. For the above equilibrium of the force in each layer in zone 1 of the typical drop off, as indicated in Fig. 5, there exists: Bottom layer (i=1): dP1 þ 2=1 ¼ 0 dx

ð8Þ

Intermediate layers (i=2, 3, . . . N2): dPi þ iþ1=i  i=i1 ¼ 0 dx

ð9Þ

Intermediate thickness variable resin layer (i=N1): dPN1 0 0 þ N=N1 þ tanN=N1  N1=N2 ¼ 0 dx 0 0 N=N1  tanN=N1 ¼ 0

ð10aÞ

or   dPN1 0 þ N=N1 1 þ tan2   N1=N2 ¼ 0 dx

Substituting Eq. (6) and formulation for ti(x) into Eq. (5) results in a series of recursive equations involving Cn

Fig. 3. Derivation for interfacial modulus.

1455

Fig. 4. Schematic drawing of typical ply drop region.

ð10bÞ

1456

K. He et al. / Composites Science and Technology 63 (2003) 1453–1462

Top inclined ply (i=N):

The governing equilibrium equations in each layer in zone 2 as shown in Fig. 6 are of the forms: Bottom layer (N=1):

dPN 0 cos2   N=N1 ¼0 dx 0 0 N=N1  N=N1 tan ¼ 0

dP1 þ 2=1 ¼ 0 dx

ð11Þ

ð12Þ

Intermediate layers (i=2, 3, . . . N3): dPi þ iþ1=i  i=i1 ¼ 0 dx

ð13Þ

Intermediate thickness variable resin layer (i=N2): dPN2 0 0 þ N1=N2 þ tanN1=N2  N2=N3 ¼ 0 dx ð14aÞ 0 0 N1=N2  tanN1=N2 ¼0 or   dPN2 0 þ N1=N2 1 þ tan2   N2=N3 ¼ 0 dx

ð14bÞ

Top inclined resin layer (i=N1):  dPN1  0 0 þ N=N1  N1=N2 dx   0 0 þ tan N=N1  N1=N2 ¼0 Fig. 5. Equilibrium conditions of ply drop at zone 1.

ð15aÞ

    0 0 0 0 N=N1  N1=N2  N1=N2  tan N=N1 ¼0 or   dPN1 0 0 cos2  þ N=N1  N1=N2 ¼0 dx

ð15bÞ

Top inclined ply (i=N): dPN 0 cos2   N=N1 ¼0 dx 0 0 N=N1  N=N1 tan ¼ 0:

In the governing equations above, Pi is a function of (E ui) which is of the same form as in Eq. (2), while  i/j is a function of (Gi, Gj, ti, tj, ui, uj) as in Eq. (3). The recursive formulation for the coefficients of polynomial was provided in Ref. [2]. The required coefficients are contained in a linear system. Coefficients in one layer are coupled with that in adjacent layer as well as the sequential dependence of former term coefficient on the latter in the same layer. If the boundary conditions such as C0i and C1i for displacements and stress/strain on one side of each layer are known, the total coefficients are i, ti,

Fig. 6. Equilibrium conditions of ply drop at zone 2.

ð16Þ

K. He et al. / Composites Science and Technology 63 (2003) 1453–1462

defined using the above recursive formulation. The convergence of the series depends on the number of terms chosen. In this model using the modified shear-lag analysis, we can also present its difference from the model using original shear-lag method. The comparison with FEM results is performed so as to determine the accuracy of the model.

The force, in each layer, can be determined by Fi ðxÞ ¼  ðxÞtðxÞ ¼ Ei "ðxÞtðxÞ ¼ Ei u0i ðxÞtðxÞ 1 X Cni n  x n1 ¼ Ei tðxÞ L L n¼1

Implementation and validation of the modeling consist of the consideration of boundary conditions/constraints, conversion of the global system to local system (tapered coordinate system) for obtaining correlation with FEM solutions, determination of interlaminar normal stress with force equilibrium considerations, and comparison of the current modeling results with others available in the literature. 3.1. Interlaminar stress distributions The actual problem of interest to be solved involves partitioning the region about the drop-off into two zones, as shown in Fig. 4. In each of these regions both continuous (belt and core) plies and resin layers are assumed to carry the extensional and shear stress. In region 2, the dropped sub-laminate is replaced with a resin layer with variable thickness, and the thickness values of the top and bottom resin layers are assumed to be equal. Each zone has ten free constants. The twozone model has 10 boundary conditions, related to five deflections at each end of the model. The additional constraints on the system are the continuity requirements at the boundary between zones 1 and 2, where displacement and force must be continuous. To solve the two-zone problem, a set of 20 equations is set up and the initial coefficients in each zone can be determined. Based on these initial coefficients, all of the coefficients can be determined. Thus, the displacements, strains, and stresses are known. All of the constraints can be reduced by the definitions of the assumed displacements. For example in layer 5 1  x n X u5 ðxÞ ¼ Cn5 ð17Þ L n¼0

C1i L 1 X nCni Fi ðLÞ ¼ Ei ti ðLÞ L n¼1

Fi ð0Þ ¼ Ei ti ð0Þ

ð20Þ

For a two-zone problem as shown in Fig. 4, total of 20 constraints must be applied. Each zone generates ten unknowns (C0i and C1i , i=1,2,0 5). The 20 constraints can be broken into constraints on displacements at the boundaries, continuity of displacements between zones, and continuity of forces between zones. For this problem, the 20 constraints are: Boundary Conditions (10). 1. Left boundary condition on each ply and resin layer (total 5). 2. Right boundary condition on each ply and resin layer (total 5). Displacement continuity conditions (5). 3. Displacement continuity between zone 1 and 2 in each ply and resin layer. Force continuity conditions (5). 4. Force continuity between zone 1 and 2 in each ply and resin layer. To demonstrate the capability of the shear-lag model developed, an application of the model to predict interlaminar stresses about the ply drop region in the uniaxially loaded tapered laminate was made. All plies in the model as shown in Fig. 7 are 0 unidirectional carbon/epoxy with material properties given in Table 1,

thus the displacements at the boundary for any layer i are given by ui ð0Þ ¼ C0i ð18Þ

n¼0

where i=1, 2, . . . 5.

ð19Þ

Thus the forces per unit width at the boundaries in each layer are given by

3. Model implementation and validation

1 X ui ðLÞ ¼ Cni

1457

Fig. 7. Schematic drawing of laminate with drop-off layers.

1458

K. He et al. / Composites Science and Technology 63 (2003) 1453–1462

Table 1 Material properties Glass/epoxy unidirectional tape

Resin

Extensional moduli, msi ET=1.8 EN=1.8 EL=6.4

Young’s modulus, msi E=0.57 msi

GLT=0.65

LT =0.29

Shear moduli, msi GTN =0.65 GNL=0.60 Poisson’s ratios LN=0.29

TN=0.50

Shear moduli, msi G=0.21 Poisson’s ratio =0.37 Tensile strength  0=9 ksi

which, for convenience of comparison and evaluation, are quoted from Refs. [11,12]. Also, thin (0.1tply thick) resin layers are included in the model to allow for the direct calculation of interlaminar stresses. In this problem, an eight-ply laminate is tapered to four plies in a symmetric fashion. The laminate is 50.8 mm (2 inches) long in the thick section, 50.8 mm (2 inches) long in the thin section, 25.4 mm (1 inch) wide and has a taper ratio of 10:1 (ply drop-off step space is 10 times ply thickness). Uniaxial tension is applied by fixing thicker end of the laminate and applying an axial displacement of 0.0254 mm (0.001 inch) at the other end. Convergence of the shear-lag model is controlled by the term at which the infinite series in Eq. (6) are terminated. Fig. 8 indicates the peak values of interlaminar shear stress in the top resin layer using 5, 10, 15 and 20 terms in the series. For this analysis the series will be truncated at 20 terms. Fig. 9 illustrates interlaminar shear peak values about the foremost ply drop region (Point C), where a critical load transferring occurs and, it further induces delamination failure at the interfaces (ABCDE) between dropped plies and continuous plies. Improvement of MSL model over OSL model is observed with direct comparison of estimations of stresses components against FEA solution. Predictions by both shear-lag

Fig. 8. Convergence for the shear-lag model.

Fig. 9. Comparison of different models for interlaminar shear peak values in the top resin layer (left) and bottom resin layer (right). FEM: finite element method; MSL: modified shear-lag model; and OSL: ordinary shear-lag model.

models agree well with finite element solution, with results obtained from MSL being closer to the stress state rendered by FEA. Over the range of dropped ply region as shown in Figs. 10 and 11, it is indicated that about 26–43% of improvement in the magnitudes of interlaminar stress components is found with the present model. 3.2. Simulation of progression of the initial delamination With the confidence in using the modified shear-lag theory shown in the analysis of stress distributions at the ply-drop location, the progressive delamination that initializes at this region is further simulated with the model. This model therefore overcomes the limitations of ordinary shear-lag model, in which the fracture of the resin layer due to a tensile load could not be addressed because of its assumption. Analysis configurations are shown in Fig. 12; where the thick lines show the progressive delamination about the ply drop region. In Fig. 12(a), named case (P), a transverse matrix crack first originates at the interface between the foremost ply drop (Point C) and the resin rich pocket that is commonly observed from previous works and experimental program detailed in Ref. [6]. This is expected since the foremost ply drop region is critical due to material discontinuities and overall thinnest thickness of the laminate. Furthermore, the interlaminar normal stress profiles at ply drop step of the laminate which are shown in Fig. 11 experience sudden change from compression at the left side of the ply drop end to tension at the other side, indicating that the interlaminar normal stress has less contribution to delamination initiation at the drop off location, while interlaminar shear stress component plays significant role in causing delamination failure at this critical region. After the initiation of the transverse cracking occurring at Point C under a lower load, the delamination further propagates into the thick section along ‘a’ and ‘b’, respectively. Delamination ‘a’ goes along the interface

K. He et al. / Composites Science and Technology 63 (2003) 1453–1462

1459

Fig. 10. Interlaminar shear stress distributions.

Fig. 11. Interlaminar normal stress distributions.

between the dropped ply and the resin layer below it, and ‘b’ goes along the tapering direction. In Fig. 12 (b), the transverse matrix crack extends vertically from the original fashion as shown in Fig. 12(a) up to the belt and core plies and then deflects with the both ends into the thickness section of the laminate, with ‘a’ being a length going along the interface between the resin layer and the core ply, and ‘b’ a length going along the interface between the resin layer and belt ply. In addition, in

order to perform simulation of delamination, a modification to the Eq. (3) is made by introduction of a new parameter ai/j as shown:    2Gi Gj  i=j ¼ i=j ui  uj ¼ i=j Ki=j ui  uj : ð21Þ ti Gj þ tj Gi The new parameter i/j is thus defined, which is equal to unity when the interface is bonded (or intact) and is zero when it is debonded.

1460

K. He et al. / Composites Science and Technology 63 (2003) 1453–1462

middle resin layer are set to zero. For case (R), in addition to the above conditions, strains at x=0 at the ends of top resin and bottom resin layers are also set to be zero. ii. Similarly, the forces at x=0 are also zero. iii. At certain debonded length x=a (or b) (say ai or bi), the displacements should be continuous; namely ui(a) and ui(b) of debonded layer/region should be equal to ui(a) and ui(b) of undebonded layer/region. iv. At certain debonded length x=a (or b) (say a or b), the interfacial (interlaminar) shear stresses are equal to interfacial (interlaminar) shear strength  int.

Fig. 12. Schematic representation of the delamination initiation configurations. Case (P): fractured ply; Case (R): fractured resin matrix.

The boundary conditions for delamination analysis by  c (=0.5  0) criterion are shown as: i. The strains at x=0 for both ends of drop ply and resin layers are zero. For case (P) as shown in Fig. 10, the strain at the right end of the dropped ply in zone 1 and the stain at the left end of the

For the simulation of the progress of delamination, first  i/j values that dominate delamination failure mechanisms of tapered structures studied in Refs. [11,12] were calculated for each bonded interfacial cell having a length two times ply thickness H. If one cell (say k) gave the maximum value among all interface cells, this cell was judged as the first interface to be debonded at the corresponding load and interfacial parameter was thus set to zero. Then the shear stress  i/j values were again calculated for all bonded interfaces, and the interface with the maximum value was searched, which was set to be the second broken interface. With repletion of this process, the broken interface was determined one after another. The sign (+, ) of  i/j is dependent on whether ui+1/j ui/j is positive or negative. This sign is reversed when the delamination direction is taken reversed. Hereafter, the absolute value is used for  i/j for convenience. Fig. 13 illustrates the progress of delamination along ‘a’ and ‘b’ shown in Fig. 11 for case (P). The closed circles in the Fig. 11 correspond to the delaminated cells at the start ‘a1’ and restart ‘b1’, ‘a3’ and ‘a10’ after increment of the applied stress, and the open circles to the delamination cells whose debonding occurs at the same stress level as that for the preceding debonding.

Fig. 13. Calculated normalized interlaminar shear stress against the order of delaminated cell number for the case (P).

K. He et al. / Composites Science and Technology 63 (2003) 1453–1462

1461

Fig. 14. Calculated relationship between the normalized interlaminar shear stress and the delaminated length.

Fig. 14 shows the calculation results for the relationship between the normalized interlaminar shear stress and debonded lengths along ‘a’ and ‘b’ for both cases. For each case, the debonded region grows as follows: in case (P), the initiation load for delamination along ‘a’ is higher than the load required for the delamination along ‘b’. The load increases further for both lengths of less than five cells (each cell length is equal to two times ply thickness). Then, with almost constant loading the cracks increase mutually for both lengths. At last, the increments of the load are required for further increase of the delamination lengths. A different trend is observed for case (R), where the initial loads for debonding along both directions is lower than those for case (P). For case (R), as the debonding progresses, the loading required decreases for both cracks. After the first five cell debondings, the loading values increase until unstable fracture occurs.

4. Conclusions Modified shear-lag model that shows significant advantages over conventional shear-lag theory, both in computational accuracy and in application capability, has been developed and implemented to conduct interlaminar stress analysis of the tapered laminate. Improvement provided by the current model over generally used shear-lag model is demonstrated in comparison with solutions obtained in the literature. This new model has been extended to perform the delamination initiation simulation using interlaminar shear criterion. This is based on the fact that the interlaminar shear

stress component dominates premature failure mechanisms of the initial delamination caused at the ply drop. Two types of cracks are considered for the purpose of identifying the delamination mechanisms. It is shown in the analysis that the load required to initiate the debonding in case (P) is larger than the one in case (R), indicating that the initial transverse resin matrix crack would decrease the overall strength of the laminate. The different debonding trends for the two cases are also observed, which again demonstrate whether or not the inclusion of resin layers at the critical locations of the laminate would significantly affect the simulation results. References [1] He K, Hoa SV, Ganesan R. The study of tapered laminated composite structures: a review. Journal of Composites Science and Technology 2000;60:2643–57. [2] He K. Interlaminar stresses and fracture behavior of thicknesstapered composite laminates. PhD thesis, Concordia University; 2002. [3] Kemp BL, Johnson ER. Response and failure analysis of graphite–epoxy laminate containing terminating internal plies. In: Proceedings of the AIAA/ASME/ASCE/AHS 26th Structures, Structural Dynamics, and Materials Conference, Pt. 1, AIAA, New York, April 1985. p. 13–24 [AIAA paper 85-0608]. [4] Curry JM, Johnson ER, Starnes Jr JH. Effect of dropped plies on the strength of graphite-epoxy laminates. AIAA Journal 1992; 30(2):449–56. [5] Hoa SV, Du BL, Vu-Khanh T. Interlaminar stresses in tapered laminates. Polymer Composites 1988;9(5):337–44. [6] Fish JC, Lee SW. Delamination of tapered composite structures. Engineering Fracture Mechanics 1989;34(1):43–54. [7] Llanos AS, Vizzini AJ. The effect of film adhesive on the delamination strength of tapered composites. Journal of Composite Materials 1992;26(13):1968–83.

1462

K. He et al. / Composites Science and Technology 63 (2003) 1453–1462

[8] Salpekar SI, Raju IS, O’Brien TK. Strain-energy-release rate analysis of delamination in a tapered laminate subjected to tension load. Journal of Composite Materials 1991;25:118–41. [9] Armanios EA, Parnas L. Delamination analysis of tapered laminated composites under tensile loading. In: O’Brien TK, editor. Composite materials: fatigue and fracture, vol. 3. ASTM STP 1110. Philadelphia: American Society for Testing and Materials; 1991. p. 340–58.

[10] Harrison PN, Johnson ER. A mixed variational formulation for interlaminar stresses in thickness-tapered composite laminates. Int J Solids Structures 1996;35(16):2377–99. [11] Vizzini AJ. Shear-lag analysis about an internally-dropped ply. Journal of Reinforced Plastics and Composites 1997;16(1). [12] He K, Hoa SV, Ganesan R. Stress analysis of tapered composite laminates using partial hybrid finite elements. Journal of Reinforced Plastics and Composites [submitted for publication].