Buckling of composite beams with two enveloped delaminations: Lower and upper bounds

Computers and Structures 86 (2008) 2155–2165

Contents lists available at ScienceDirect

Computers and Structures journal homepage: www.elsevier.com/locate/compstruc

Buckling of composite beams with two enveloped delaminations: Lower and upper bounds M.S.R. Parlapalli a, Dongwei Shu b,*, Gin B. Chai b a b

Production Technology, Department of Mechanical Engineering (CTW), University of Twente, P.O. Box 217, 7500 AE, Enschede, The Netherlands School of Mechanical and Aerospace Engineering, Nanyang Technological University, 50 Nanyang Avenue, Singapore 639798, Singapore

a r t i c l e

i n f o

Article history: Received 14 June 2007 Accepted 5 June 2008 Available online 19 August 2008 Keywords: Buckling Delamination Composites Finite element analysis Mode shapes

a b s t r a c t Lower and upper bounds of the buckling load of a composite beam with two enveloped delaminations are obtained from newly developed analytical models. The characteristic equation, governing the delamination buckling is derived by using Euler–Bernoulli beam and classical lamination theory, performing proper linearization and by imposing the appropriate continuity and boundary conditions. Linear and non-linear finite element analyses carried out to verify the accuracy of the developed model. S-shaped buckling mode dominates more than a single hump buckling mode in the case of upper bounds, whereas, for lower bounds of the buckling load, mixed and/or local buckling modes are observed. Ó 2008 Elsevier Ltd. All rights reserved.

1. Introduction Fibre debonding [1], broken fibers, delaminations [2] and micro cracks of the matrix [3] are some of the failure mechanisms that contribute to the property degradation of the composite structures. Out of these failure mechanisms, delaminations play a critical role when the structure is under compressive loading. Delaminations arise at the interlaminar region by impact loading, fatigue and/or poor manufacturing process [3]. Under compressive loading, a delaminated composite laminate may buckle and possibly undergo propagation of delamination. This propagation can leads to lowering of the buckling load of the laminate and to an unexpected structural failure at loads below the design level. The delamination buckling analyses become further intricate as the delaminations vary in number (single or multiple), types (enveloping, enveloped, etc.), locations (spanwise or thicknesswise) and/or shapes (embedded or through-width). Chai et al. [4] studied the buckling analysis of delaminated composite laminates by thin film approximation where local buckling of a surface delamination was studied analytically. The local delamination growth, stability and arrest were governed by a fracture mechanics-based energy release rate. The energy-release rate associated with delamination growth was computed by numerical differentiation of the total potential energy with respect to delamination length. Thin film approximation neglects any bending deformation in the base plate. Bottega and Maewal [5] analyzed * Corresponding author. Tel.: +65 67904440; fax: +65 67911859. E-mail address: [email protected] (D. Shu). 0045-7949/$ - see front matter Ó 2008 Elsevier Ltd. All rights reserved. doi:10.1016/j.compstruc.2008.06.008

the axisymmetric buckling and growth in a compressively loaded two-layer symmetric isotropic clamped circular plate with a concentric penny-shaped delamination. The influence of imperfections in the form of transverse loads were also examined and found that the imperfections do not have any noticeable effect on the delamination process. They developed an analytical model based on asymptotic analysis of post-buckling behavior. They found that delamination growth based on a Griffith-type fracture criterion might have occurred following the delamination buckling, provided that sufficient bending energy was produced at the delamination crack front. Simitses et al. [6] developed a onedimensional model for predicting buckling loads of homogenous, elastic, laminated clamped and simply-supported beams. The delamination was assumed to exist prior to the applying the compressive load and it divides the column (or plate) into four regions, each having such dimensions that Euler–Bernoulli beam theory could apply. They observed that buckling load decreases as the delamination length increases, for given delamination positions. Yin et al. [7] presented an elastic buckling and post buckling analysis of an axially loaded beam-plate with through-width delamination, symmetrically located at an arbitrary depth. They found that for a relatively short and thick delamination, the buckling load of the delaminated plate is a lower bound of the ultimate axial load capacity. Kassapoglou [8] analyzed the buckling and postbuckling composite laminates having elliptical delamination by a perturbation method. von Karman non-linear strain-displacement equations were used in the classical laminated plate theory. The governing differential equations were reset into three partial differential equations in three displacements and solved as an

2156

M.S.R. Parlapalli et al. / Computers and Structures 86 (2008) 2155–2165

eigenvalue problem for buckling load. He proposed critical delamination length, which was defined as the length below which the panel structure under compression was not affected by the presence of delamination. He pointed out that the critical delamination length is important in design, since it represents the maximum delamination size that can be tolerated in a panel without lowering its compression strength. His method has yielded buckling loads with large errors. Jane and Yin [9] developed an analytical procedure for computing the buckling loads and the postbuckling solutions of cross-ply and angle-ply laminates containing the elliptical delaminations by thin film approximation. Rayleigh–Ritz method and von Karman’s non-linear theory of plates have been used in their analysis. The transverse shear effect appreciably reduces the buckling load and significantly increases the energy release rate. Although higher-order solutions are desirable for improved accuracy, the computing storage and time requirement to evaluate the integrals in the minimum potential energy expression are large. Shu and Mai [10] modified the classical beam theory with the relaxation to the plane-section assumption by introducing rigid connector (plane sections at the ends of the delaminations are assumed to be remain plane and perpendicular to midplane of the beam during buckling) and soft connector (plane sections at the ends of the delaminations are assumed to be remain plane). They found that the soft connector analysis gives lower bound values where as the rigid connector analysis gives upper bound values. They have not discussed in details the intermediate load values. Kutlu and Chang [11] investigated analytically and experimentally the compression response of laminated composites containing multiple through-width delaminations. A non-linear finite element code based on the updated Lagrangian formulae was developed. Experiments were carried out on Fiberite T300/976 graphite/epoxy composites to verify the model and numerical calculations. The interface friction between the upper and lower surfaces of delaminations is not considered in this analysis. Lee et al. [12] studied the buckling of axially loaded composite beam plates with multiple through-width delaminations by developing a displacementbased, one-dimensional finite element model. Governing equations were derived from the principle of the stationary value of potential energy. The generalized displacements were expressed over each element as a linear combination of the one-dimensional Lagrangian interpolation functions. Both symmetric and antisymmetric buckling modes were obtained using a half-model of the composite. Lim and Parsons [13] developed a mathematical model to predict the buckling load of a composite beam with multiple through-width double equal delaminations. Symmetrical delaminations were considered and Rayleigh–Ritz method has been employed. Lagrangian multipliers were used to enforce constraints arising from the boundary conditions and the kinematic continuity requirements between the different sections of the beam. Antisymmetric s-shaped buckling mode was observed for clamped single delaminated beam at the delamination lengths of 0.35L and 0.4L. They observed that for the single delamination case, Euler buckling dominates if the delamination is sufficiently short. For longer, shallow delaminations, the thin film buckling dominates. The contact between laminate has not been considered in this analysis. Adan et al. [14] developed an analytical model for buckling of multiple delaminated composite beams. The buckling equations were formulated by the principle of virtual displacements and considering von Karman strains, which were solved by a perturbation method. They observed that when the delamination boundary is located over the mid-span of the other delamination, the decrease in the buckling load is maximal and the effect of interaction decays, as the overlapping zone between delaminations reduces. The theoretical singularity at the delaminations, plasticity effects interlaminar penetrations (contact effects) and end effects, which may occur,

have not considered in this work. Kyoung et al. [15] studied the buckling and postbuckling analysis of multiple through-width delaminated composites by a non-linear finite element analysis. The updated Lagrangian description and modified arc-length method were adopted in finite element formulation. Eight-node degenerated shell elements located along the mid-plane of these divided regions were used. First order shear deformation theory was implemented. The interference of the layers was prevented by employing the contact node pair (node-to-node contact was assumed in postbuckling behavior) defined by use of virtual beam element. They observed that multiple through-width and the multiple embedded delaminations lowers the buckling load of the cross-ply laminates due to the reduction of bending stiffness. Huang and Kardomateas [16] developed an analytical method for predicting the critical buckling loads of a composite beam-plate with double enveloped delaminations. They observed that the deeply located delamination affects the critical load significantly, when it is close to the nearest-to-the-surface delamination or when it is longer. When the two delaminations are too close to each other, the middle subplate buckles initially. This model has not studied delaminated layers made of different materials. Shu [17] performed an exact buckling analysis for beams with double equal delaminations. He observed that the clamped delaminated beams withstand lesser buckling loads then the simply-supported delaminated beams. Lachaud et al. [18] carried out experimental and numerical studies of delaminated thermoplastic carbon-fiber composites and thermoset composites. The compression test samples were cutout from T300/914 carbon epoxy and AS4/PEEK carbon (poly ether ether ketone) plates. Circular Teflon films (aluminum for the thermoplastic) of 20 lm thickness (12 lm for aluminum) were introduced during the fabrication, between plies of different orientation in order to create a delamination. They observed that the load at which local buckling appears, increases with the difference in angle between the buckled plies and the loading axis as well as with thickness of the buckled plies. Hu [19] analyzed the buckling analysis of a delaminated laminate with FEM based on the Mindlin plate theory. A sensitivity-based updating algorithm incorporating the quadratic programming technique was presented to solve the fictitious forces in contacting area. The fictitious forces were transferred into artificial springs, which prevent the penetration between two delaminated layers. The buckling load increases when the effect of the contact is considered. When the buckling mode changes significantly during the contact iteration, the buckling load also varies to a great extent. Gu and Chattopadhyay [20] studied the mechanics and mechanisms of delamination buckling and postbuckling of composites. Compression tests were carried out on HYE-3574 OH Graphite/epoxy with built-in-delaminations in order to evaluate the critical load and the actual postbuckling load-carrying capacity. Three types of flat-plate test specimen with single centered delamination locating at midplane and near surfaces with varying ply sequences were used. They found that the delamination-buckling mode was closely related to the location and length of the delamination. For composites with thinner and longer delaminations, the ultimate load obtained through the experiment is as high as three times their critical load. Sekine et al. [21] investigated the buckling analysis of elliptically delaminated composite laminates under in-plane compressive loads by taking into account of partial closure of delamination. A penalty function method was proposed to deal with the contact problem for the partial closure of the delamination. Fictitious springs were inserted iteratively at all of the overlapped points. They observed that the buckling load of a laminate with delaminations oblong in the load direction is higher than the buckling load of a laminate with delaminations oblong in the direction perpendicular to the load direction. They found that the buckling load of a laminate containing small number of delamina-

2157

M.S.R. Parlapalli et al. / Computers and Structures 86 (2008) 2155–2165

tions is higher than the buckling load of a laminate containing a hole. Hwang and Mao [22] used non-linear buckling and postbuckling analyses to study the through-width delamination buckling of interply hybrid composite plates under compression. The total strain-energy release rate and interlaminar-stress criterion were used to predict the delamination growth load. ANSYS 5.3 FEM software was used to calculate the buckling load. For global buckling mode, the buckling loads are suggested to represent the failure loads of delaminated composites. For mixed and local buckling modes, the failure of the delaminated composite is attributed to the growth of delamination. An increase in the number of glassfiber layers within an interply hybrid composite results in a decrease in the buckling load and the delamination failure load. More recently, the current authors studied the delamination buckling of two-and three-layer delaminated beams in terms of newly defined non-dimensionalized axial and bending stiffnesses and effectiveslenderness ratio [23,24]. They observed that the effective-slenderness ratio is found to be a controlling parameter of the type of delamination buckling mode configuration (global, mixed or local). Various researchers studied some important topics on delamina-

tion such as non-linear phenomena [25,26] as well as dynamic effect [27,28] and wave propagation [29]. Analytical models to obtain lower and upper bounds of the buckling load of composite beams with two enveloped delaminations have not yet been reported. These two bounds will be useful to gauge the working range of the bridging [30] and give guidelines for practical applications [31]. The objective of this work is to present analytical models based on the Euler–Bernoulli beam theory to analyze the buckling of composite beams with two enveloping delaminations. In addition, a parametric study on the influence of the lengths and locations of the enveloped delaminations on the buckling load is carried out. The present work also includes the ‘constrained-mode’ and the ‘free mode’ models. The effects of the differential-stretching and the bending–extension coupling are considered. 2. Formulation Fig. 1 shows the geometry of a composite beam with length L, thickness H and having two enveloped delaminations. In Fig. 1a,

SubC beam 6

P Virgin-beam 7

Sub-beam 5 D

E Sub-beam 3

A x1

H5

B

Virgin-beam 1

x w(x)

V

Sub-beam 4

Subbeam 2

P

x4 IV

x2

d2

F x6

x5

b

H1=H7=H

II

H3

I

a d1 III

x3

H4

a

Y3

Y4

Sub-beam 4 Sub-beam 2

Virgin-

B

a

Sub-beam 5

P beam 1

C Subbeam 6

b

D

Y2 Y1=Y7

b

Virginbeam 7

Y5

L

P

E

Sub-beam 3

A

F L

c

1 2

z0 z2 z1

midplane

x, u

zk-1

Hi k

zk zni-1

ni

zni

z, w Fig. 1. A composite beam with two enveloped delaminations: (a) under axial compressive loads; (b) equivalent beam model (shown with neutral axes) and (c) the ith beam laminate.

2158

M.S.R. Parlapalli et al. / Computers and Structures 86 (2008) 2155–2165

the longer delamination envelops the shorter delamination, which will be referred as ‘enveloping delaminations’. The two pre-existing through the width enveloped delaminations are denoted as delamination-I (of length a) and delamination-II (of length b) and are located at offset distances of d1 and d2 from the center of the beam, respectively. Considering one-dimensional through the width delaminations conveniently defines a unit width of the composite beam. Both the ends of the beam are assumed to be clamped and an external load, P, is applied along the neutral axis of the beam. Due to the presence of the two enveloped delaminations, the composite beam is analyzed as seven interconnected beams, i.e., virgin beams 1 and 7 (non-delaminated portion), and subbeams 2, 3, 4, 5 and 6 (delaminated portion) as shown in Fig. 1a and b. The virgin and sub-beams have thicknesses of Hi, lengths of Li (i = 1–7) and neutral axial positions of Y1 = Y7 = H/2, Y2 = Y6 = H2/2, Y3 = H4 + H5 + (H3/2), Y4 = H4/2 and Y5 = H4 + (H5/2), respectively. These sub- and virgin-beams can be treated as classical Euler–Bernoulli beams, provided that L  H [32]. In Fig. 1a, notations A and F denote the clamped ends and B, C, D and E denote the ends of the delaminations. The coordinate axes for the sub-beams are shown in Fig. 1a and x is measured from the left end of the composite beam. wi(x) denotes the small elastic deflection of the neutral axis of the ith (i = 1–7) beam segment. To simplify the analysis it is assumed that the compressive load is uniform and uniaxial, the enveloped delaminations exist before the application of the load and the interference between the laminated layers is neglected. Lower and upper bounds of the buckling loads will be obtained by means of ‘free’ and ‘constrained’ modes, respectively [10]. For a ‘free-mode’, it is assumed that the delaminated layers deform ‘freely’ without touching each other. For a ‘constrained-mode’, the delaminated layers are assumed to be in touch along their whole length all the time, but are allowed to slide over each other. Further, it is assumed that the interlaminar bond is good such that buckling precedes delamination propagation/ extension. The following sections present the formulations of the ‘freemode’ and the ‘constrained-mode’ models for the buckling of a composite beam with two enveloping delaminations. 2.1. Free-mode model: lower bounds of the critical buckling load The equations governing buckling of composite beam with two enveloped delaminations using the Euler–Bernoulli beam theory are 0000

Di wi ðxÞ þ Pi w00i ðxÞ ¼ 0 ði ¼ 1—7Þ

ð1Þ

where Di is the reduced bending stiffness of the ith beam segment, 0000 2 4 w00i ðxÞ ¼ d wi ðxÞ=dx2 , wi ðxÞ ¼ d wi ðxÞ=dx4 , Pi is the axial load (P = P1 = P7). The mechanical properties of the composite beams are determined using the classical laminate theory [33]. Thus

Di ¼ Di11 

ðBi11 Þ2

ð2Þ

Ai11

where

Ai11 ¼ b Bi11 ¼ Di11

b 2

ni  X k¼1 ni X

Q k11

 k

ðzk  zk1 Þ

  ðQ k11 Þk z2k  z2k1

k¼1

  4 2 Q k11 ¼ Q k11 cos4 / þ Q k22 sin / þ 2 Q k12 þ 2Q k66 cos2 / sin / E11 ; 1  m12 m21

Q 22 ¼

E22 ; 1  m12 m21

wi ðxÞ ¼ Ai cosðki xÞ þ Bi sinðki xÞ þ C i x þ Di

ði ¼ 1—7Þ

ð3Þ

pffiffiffiffiffiffiffiffiffiffiffi where ki ¼ Pi =Di and Ai, Bi, Ci, and Di are 28 integration constants. The critical buckling load of a composite beam with two enveloped delaminations will be obtained by solving for the integration constants (Ai, Bi, Ci, and Di) as an eigenvalue problem by imposing proper boundary conditions and appropriate kinematic continuity conditions. The lowest eigenvalue gives the critical buckling load. For clamped supported ends at x = x1 and x = x6, the appropriate boundary conditions are wi = 0 and w0i ¼ 0, where i = 1 and 7, and w0i ¼ dwi =dx. The continuity conditions for deflection and slope at x = x2 are given by

w1 ðxÞ ¼ w2 ðxÞ;

w1 ðxÞ ¼ w3 ðxÞ

ð4Þ

w01 ðxÞ ¼ w02 ðxÞ;

w01 ðxÞ ¼ w03 ðxÞ

ð5Þ

The continuity for shear and bending moments at x = x2 are given by

V 1 ðxÞ ¼ V 20 ðxÞ þ V 30 ðxÞ

ð6Þ

M1 ðxÞ ¼ M 20 ðxÞ þ M 30 ðxÞ  DP2 ðY 1  Y 2 Þ þ DP3 ðY 1  Y 3 Þ

ð7Þ

where Mi ðxÞ ¼ Di w00i ðxÞ, V i ðxÞ ¼ Di w000 i ðxÞ, Vi0 and Mi0 are the shear force and bending moment at the delamination ends, DP2 and DP3 are the incremental axial forces in the sub-beams 2 and 3, arising from the axial extension/compression [24], respectively. DP2 and DP3 are solved from the compatibility between the stretching/ shortening of the delaminated layers and axial equilibrium [24]. In Eq. (7), the third and fourth terms on the right hand side represents the contribution of the differential stretching between the sub-beams 2 and 3 to the bending moment [10]. In the same process, continuity conditions can be obtained at x = x3 between the sub-beams 2, 4 and 5; at x = x4 between the sub-beams 4, 5 and 6; at x = x5 between the sub-beams 3, 6 and the virgin beam 7. The boundary conditions and the continuity conditions provide 28 homogeneous equations to solve for the 28 integration constants (Ai, Bi, Ci, and Di, i = 1–7). The coefficients of the 28 integration constants from the above equations (Eqs. (4)–(7)), equations from the boundary conditions and the kinematic continuity conditions at the delaminations fonts x = x3, x = x4 and x = x5 are written in a matrix form. For a non-trivial solution to exist the determinant of the coefficient matrix must vanish and the lowest eigenvalue gives the non-dimensionalized buckling load, Pð¼ P=Pcr ; Pcr ¼ 4p2 EI=L2 for clamped endsÞ. 2.2. Constrained-mode model: upper bounds of the critical buckling load

ni    bX ¼ Q k11 z3k  z3k1 k 3 k¼1

Q 11 ¼

where Ai11 is the extensional stiffness, Bi11 is the coupling stiffness, Di11 is the bending stiffness, Q k11 is the coefficient stiffness of the kth lamina, b is the width, ni is the number of plies, E11 and E22 are the longitudinal and transverse Young’s moduli, respectively. G12 is the in-plane shear modulus, m12 and m21 are the longitudinal and transverse Poisson’s ratios, respectively. / is the angle of the kth lamina and zk and zk1 are the locations of the kth lamina with respect to the midplane of the ith beam (Fig. 1c). The second term on the right hand side of Eq. (2) represents the reduction of stiffness due to the bending–extension coupling [33]. The general solutions of Eq. (1) are

Q 66 ¼ G12 ;

m21 ¼

m12 E22 E11

For ‘constrained-mode model’ it is assumed that the delaminated layers are ‘constrained’ to have the same transverse deformations [10]. Therefore, the sub-beams 2 and 3 have the same transverse deformation (x2 6 x 6 x3). Similarly, sub-beams 4, 5 and 3(x3 6 x 6 x4), and sub-beams 6 and 3(x4 6 x 6 x5) have the same transverse deformation, due to which the composite beam is modelled as five sub-beams, which are denoted as I, II, III, IV and V in Fig. 1b. Sub-beam II is made of sub-beams 2 and 3. Simi-

M.S.R. Parlapalli et al. / Computers and Structures 86 (2008) 2155–2165

larly, sub-beam III is made of sub-beams 4, 5 and 3 and sub-beam IV is made of sub-beams 6 and 3. In the case of the constrainedmode model the equation governing buckling of virgin beam I (x1 6 x 6 x2) is 0000

D1 w1 ðxÞ þ P1 w001 ðxÞ ¼ 0

ð8Þ

For sub-beam II (x2 6 x 6 x3), we have 0000

ðD2 þ D3 Þw2 ðxÞ þ ðP2 þ P3 Þw002 ðxÞ ¼ 0

ð9Þ

For sub-beam III (x3 6 x 6 x4), we have 0000

ðD3 þ D4 þ D5 Þw4 ðxÞ þ ðP3 þ P4 þ P5 Þw004 ðxÞ ¼ 0

ð10Þ

For sub-beam IV (x4 6 x 6 x5), we have 0000

ðD3 þ D6 Þw6 ðxÞ þ ðP3 þ P6 Þw006 ðxÞ ¼ 0

ð11Þ

and for virgin-beam V (x5 6 x 6 x6), we have 0000

D7 w7 ðxÞ þ P7 w007 ðxÞ ¼ 0

ð12Þ

The generalized solutions for the ‘constrained-mode model’ are identical in form to the ‘free-mode model’. From the general solutions of the virgin- and sub-beams, the integration constants are reduced to 20 (Ai, Bi, Ci, and Di, i = 1, 2, 4, 6, and 7) and it requires 20 equations to solve for these integration constants. These twenty equations will be obtained from the 4 boundary and 16 continuity conditions (at x = x2, x = x3, x = x4 and x = x5), which are discussed below. The continuity conditions for deflection, slope, shear and bending moments at x = x2 are

w1 ðxÞ ¼ w2 ðxÞ

ð13Þ

w01 ðxÞ ¼ w02 ðxÞ 000 D1 w000 1 ðxÞ ¼ ðD2 þ D3 Þw2 ðxÞ 00 00 D1 w1 ðxÞ ¼ ðD2 þ D3 Þw2 ðxÞ

ð14Þ ð15Þ ð16Þ

Similarly, the continuity conditions at x = x3, x = x4 and x = x5 can be derived. The coefficients of the integration constants from the above equations (Eqs. (13)–(16)), 4 equations from the boundary conditions and 12 continuity conditions at the delaminations fonts x = x3, x = x4 and x = x5 are written in a matrix form and is solved as an eigenvalue problem.

2159

G12 = 6.2 GPa, m12 = 0.23 [22]. Linear and non-linear finite element analyses are carried out by using ANSYS 9.0 software. The composite beam is modeled with PLANE42 elements by considering extra shape functions options that could strongly influence the delaminated mode shapes. The PLANE42 element can be used either as a plane element (plane stress or plane strain) or as an axisymmetric element. The PLANE42 element is defined by four nodes having two degrees of freedom at each node: translations in the nodal X and Y directions. The element has plasticity, creep, swelling, stress stiffening, large deflection, and large strain capabilities. For nonlinear FEA, contact elements, CONTA171 and TARGE169 are used to form a contact pair in the enveloped delaminations region to prevent the physically impermissible penetration. CONTA171 is used to represent contact and sliding between 2-D ‘‘target” surfaces (TARGE169) and a deformable surface, defined by this element. Contact occurs when the element surface penetrates one of the target segment elements (TARGE169) on a specified target surface. Coulomb and shear stress friction is allowed. The correct node ordering of the contact element is critical for proper detection of contact. The nodes must be ordered such that the target must lie to the right side of the contact element when moving from the first contact element node to the second contact element node. The contact elements themselves overlay the solid elements describing the boundary of a deformable body and are potentially in contact with the target surface, defined by TARGE169. This target surface is discretized by a set of target segment elements (TARGE169) and is paired with its associated contact surface via a shared real constant set. This real constant set includes all real constants for both the target and contact elements. Each target surface can be associated with only one contact surface, and vice-versa. Similar to the non-linear FEA of Yap and Chai [35] for the buckling of composite laminates with a single delamination, three kinds of small deformations as shown in Fig. 3a–c are input to three identical finite element models before the axial compression to obtain the desired three types of buckling mode shapes, i.e., local, s-shaped and global. The s-shaped buckling mode for enveloped delaminations

a

3. Finite element analysis The finite element model of a composite beam with two enveloped delaminations is illustrated in Fig. 2. The clamped composite beam is 280 mm in length and 3 mm in thickness. Each ply thickness is 0.15 mm hence the composite beam is made of 20 unidirectional plies (6 layers of carbon/epoxy, 6 layers of glass/epoxy and 8 layers of carbon/epoxy). The material properties of glass/epoxy plies are: E11 = 37.9 GPa, E22 = 14.3 GPa, G12 = 5.6 GPa, m12 = 0.29, and that of carbon/epoxy plies are: E11 = 121 GPa, E22 = 9.4 GPa,

b

c

Fig. 2. Finite element model of a composite beam with two enveloped delaminations.

Fig. 3. Three types of delamination buckling modes: (a) local; (b) s-mode and (c) global.

2160

M.S.R. Parlapalli et al. / Computers and Structures 86 (2008) 2155–2165

was reported in [24]. Then, the three buckling loads are compared and the lowest load is the critical buckling load and the corresponding buckling mode will be the actual one. A parametric study (in terms of the lengths and locations of the delaminations) has been carried out to validate the developed analytical models. 4. Results and discussion For the verification of the accuracy of the present model, comparisons with published results on enveloped delaminations [16] and separately carried out linear and non-linear finite element analysis using ANSYS 9.0 are made.

Table 1 Comparison between results (P/Pcr) of the present analytical model and Huang and Kardomateas [16] for homogenous beams having two enveloped delaminations (H3 = 0.75H, H4 = 0.125H) with clamped ends b/L

0.2 0.4 0.6 0.7 1.0

a/L = 0.1

a/L = 0.3

a/L = 0.5

Present

Ref. [16]

Present

Ref. [16]

Present

Ref. [16]

0.9627 0.3754 0.1709 0.0969 0.0625

0.9628 0.3754 0.1709 0.0970 0.0625

0.1730 0.1508 0.1191 0.0852 0.0602

0.1731 0.1508 0.1191 0.0867 0.0602

0.0624 0.0624 0.0571 0.0497 0.0431

0.0624 0.0624 0.0571 0.0497 0.0431

4.1. Verification 1: homogeneous delaminated beam with two enveloped delaminations The accuracy of the present model is verified by comparing the computed non-dimensionalized buckling loads (P/Pcr) with Huang and Kardomateas [16] for the case of a homogeneous beam with two enveloped delaminations. The two enveloped delaminations are located at H3 = 0.75H, H4 = 0.125H, and d1 = d2 = 0. Table 1 shows an excellent agreement between the two results. 4.2. Verification 2: comparisons with linear and non-linear finite element analyses results After successful verification of the present analytical model for a homogeneous beam cases, non-dimensionalized buckling loads are computed for a three-layer beam and compared with linear and non-linear FEA results. The delamination-II is located at H3 = 0.4H whereas the delamination-I is at H4 = 0.3H. A close agreement between analytical and linear FEA results is observed as shown in Table 2. Whereas the buckling loads obtained from the non-linear FEA are higher than those from the analytical model and linear FEA. This is because using the linear FEA and the present analytical model, one can only get the critical buckling loads in which the sub-laminate buckles as shown in Fig. 4a. It is not cable of obtaining the actual buckling load when the entire laminate including the sub-laminate buckles. On the other hand, the non-

Table 2 Comparison between the non-dimensionalized buckling loads (P/Pcr) of present analytical model and FEA (linear and non-linear) for a composite beam with two enveloped delaminations, H3 = 0.4H and H4 = 0.3H b/L

a/L = 0.1 Present (free-mode)

0.1 0.3 0.5 0.6 0.7 0.8 0.9

0.9999 0.7568 0.4695 0.3377 0.2539 0.1980 0.1589

a/L = 0.3 FEA

Present (free-mode)

Linear

Non-linear

0.9994 0.7508 0.4620 0.3330 0.2508 0.1959 0.1574

0.9994 0.7803 0.5609 0.4074 0.3190 0.2608 0.2209

0.7480 0.5324 0.4562 0.3370 0.2538 0.1980 0.1589

a/L = 0.5 FEA Linear

Non-linear

0.7228 0.5273 0.4557 0.3324 0.2508 0.1959 0.1574

0.8784 0.5545 0.4231 0.3763 0.3075 0.2502 0.2100

Present (free-mode)

FEA Linear

Non-linear

0.2844 0.2833 0.2796 0.2534 0.2313 0.1972 0.1588

0.2792 0.2791 0.2774 0.2519 0.2302 0.1953 0.1573

0.6200 0.4977 0.3464 0.2726 0.2521 0.1982 0.1741

Fig. 4. Three types of buckling modes from the non-linear FEA: (a) local (a = 0.5L, b = 0.7L); (b) s-mode (a = 0.5L, b = 0.3L) and (c) global (a = 0.1L, b = 0.5L).

2161

M.S.R. Parlapalli et al. / Computers and Structures 86 (2008) 2155–2165

00

local or s-mode buckling shapes as shown in Fig. 4b–c. It is noted that the longer delamination of the two enveloped delaminations dominates the buckling behavior.

H4

Interface 1

900 Interface 2 00 Interface 3 900 Interface 4

4.3. A composite beam with two enveloped delaminations

900

The most useful feature of the analytical models described here is their usefulness to study rapidly the linearized buckling of the composite beams with two enveloped delaminations. As it is rather difficult to consider all the possible combinations of the lengths and locations (along the spanwise and thicknesswise) of the delaminations (delamination-I and -II) and to study their effect on the lower and upper bounds of the buckling load, selected results have been presented here. The composite beam is made of AS/ 3501 graphite/epoxy layers having a stacking sequence of [00/ 900]2s and material properties for the lamina as: E11 = 181 GPa, E22 = 10.3 GPa, G12 = 7.17 GPa and m12 = 0.28 [12]. The 8-ply beam has a thickness of 1.016 mm. The thicknesswise locations of the delaminations are shown in Fig. 5.

midplane

900 Interface 5 00 Interface 6 Interface 7

00

H3

Fig. 5. Thicknesswise locations of the two enveloped delaminations for a [00/900]2s AS/3501 graphite/epoxy composite laminate.

linear FEA is able to simulate buckling of the entire laminate and obtain the desired overall buckling load, resulting in either global,

a

1.0

Constrained model 1. a=0.2L 2. a=0.4L 3. a=0.6L Free model 4. a=0.2L 5. a=0.4L 6. a=0.6L

Free model Constrained model

0.8

1 2 4

P/Pcr

0.6

3

0.4 H4=0.125H a b

0.2 0.0 0.2

0.3

5 6

H3=0.125H

0.4

0.5

0.6

0.7

0.8

0.9

a = 0.2L, b = 0.2125L (Free-mode model)

(ii) a = 0.2L,b = 0.6L(Free-mode model)

Normalized modeshape

b (i)

Normalized modeshape

Delamination-II length (b/L)

1.0 3 0.5

4 2 1

0.0 0.0

6

5

7 0.5

1.0

1.0 3 0.5

12 0.0 0.0

4

5 0.5

6 7 1.0

Normalized distance (x/L)

Normalized distance (x/L)

Normalized modeshape

(iii) a = 0.2L, b = 0.6L (Constrained-mode model)

III

1.0 II

IV

0.5 I 0.0 0.0

V 0.5

1.0

Normalized distance (x/L)

Fig. 6. (a). Influence of the lengths of the two enveloped delaminations on P, H4 = 0.125H = H3, d1 = d2 = 0.0; (b) buckling mode shapes.

2162

M.S.R. Parlapalli et al. / Computers and Structures 86 (2008) 2155–2165

Fig. 7. Influence of the lengths of the enveloped delaminations on P, H4 = 0.375H = H3, d1 = d2 = 0.0.

buckling mode is observed as shown in Fig. 6b (iii). As the two enveloped delaminations occur at different interfaces, the two sides of the delamination may be totally separated (free-mode) or in full contact (constrained-mode) or in partial contact (in between free- and constrained-modes) with each other during the buckling, these two bounds give the guidelines for practical applications [31]. Fig. 7 shows the effect of the lengths of the two enveloped delaminations on the non-dimensionalized buckling loads ðPÞ, when H4 = 0.25H = H3. Though a similar trend to the case of H3 = H4 = 0.125H (Fig. 6a) is observed here, the gap between the lower and upper bounds of the buckling load decreases as the two enveloped delaminations are located deeply inside of the composite beam. In the case of the lower bounds of the buckling load, not all the curves have merged into a single line indicating a strong influence of the shorter delaminations on the buckling loads. This observation is different from the case discussed earlier (Fig. 6a). Hence, as the thicknesswise locations of the two enveloped delaminations vary, the buckling load is also strongly depends on the lengths of the shorter delamination.

4.3.1. Effect of the lengths of the two enveloped delaminations Fig. 6a shows the effect of the lengths of the two enveloped delaminations on the lower and upper bounds of the non-dimensionalized buckling loads ðPÞ of a composite clamped beam which are obtained from the ‘free-mode’ and ‘constrained-mode’ models, respectively. The two enveloped delaminations are located symmetrically at the center of the beam, i.e., d1 = d2 = 0 and H4 = 0.125H = H3. Pð¼ P=Pcr Þ is the ratio of the buckling load of a composite beam with two enveloped delaminations and the buckling load of a sound composite beam (no delaminations) of the same size. It is observed from Fig. 6a that as the length of the two enveloped delaminations increases the upper and lower bounds of the buckling load decreases, which is due to the decreasing of the stiffness of the composite beam. In the case of the lower bounds of the buckling load all the curves are merged into a single line indicating that the longer delamination of the two enveloping delaminations influences the buckling load more than that of the shorter delamination. In the case of the upper bounds of the buckling load, a monotonic variation is observed and not all the curves have merged into a single line. This difference between the buckling loads of the ‘constrained-mode’ and ‘free-mode’ model is largely depends on the opening of the mode shapes of the subbeams. Further, in the case of the ‘free-mode’ model the sub-beams are allowed to deform freely whereas in the case of the ‘constrained-mode’ model the sub-beams are constrained to deform together due to which they withstand higher buckling loads. The gap between the upper and lower bounds of the buckling loads is higher for shorter enveloped delaminations. As the lower and upper bounds of the buckling load are the minimum and the maximum loads that a composite beam can with stand, these bounds can be useful to know the effect working range of the bridging [30] that will be a possible future work. Fig. 6b shows the buckling mode shapes for different delaminations lengths. In Fig. 6b, X-axis contains the normalized spanwise locations of the two enveloped delaminations, Y-axis contains the normalized deflections of the neutral axes of the sub-beams, and the vertical lines on the curves indicate the lengths of the individual sub-beams. The deflections of all the sub-beams are non-dimensionalized by dividing with the maximum deflection of a sub-beam. At shorter lengths of the delamination-I (a = 0.2L and b = 0.2125L), mixed buckling mode is observed (Fig. 6b (i)), at the longer lengths of the delamination-I (b = 0.6L) local buckling of the sub-beam 3 is noted (Fig. 6b (ii)). In the case of the upper bounds of the buckling load, single hump

4.3.2. Effect of the thicknesswise locations of the two enveloped delaminations Fig. 8a shows the effect of the thicknesswise locations of the delamination -I (H4 = 0.125H, 0.25H, 0.375H and 0.5H) on the non-dimensionalized buckling loads ðPÞ for various lengths of the delamination-II. The delamination-II is located at H3 = 0.125H and d1 = d2 = 0. The length of the delamination-I is 0.2L. It is observed that the lower and upper bounds of the buckling loads reduces as the length of the delamination-II increases for given locations of the delamination-I. In the case of the upper bounds of the buckling load as the delamination-I is located deep inside, i.e., H4 = 0.375H and H4 = 0.5H, and b < 0.45L, s-shaped buckling mode as shown in Fig. 8b (i) dominates the buckling behavior than a single hump buckling mode as shown in Fig. 8b (ii). As the length of the delamination-II increases further i.e., b > 0.45L, single hump buckling mode as shown in Fig. 8b (iii) dominates the buckling behavior. Whereas, in the case of a composite beam with two separated delaminations s-shaped buckling modes have not observed [34]. The thicknesswise locations of the delamination-I strongly influence the lower and upper bounds of the buckling load. As the thicknesswise locations of the enveloped delaminations vary, the effects of bending-stretching coupling also plays a dominate role since the sub-beams consists of unsymmetrical laminae. When the delamination-II is located at H3 = 0.125H and when the delamination-I is located at H4 = 0.5H the lower bounds of the buckling load are minimum. Whereas the upper bounds of the buckling load are minimum when the delamination-I is located at H4 = 0.375H when b < 0.23L and at H4 = 0.5H when b > 0.23L, respectively. Similarly, the maximum values of the lower bounds of the buckling load are observed when the delamination-I is located at 0.125H. Whereas the upper bounds of the buckling load are maximum when the delamination-I is located at H4 = 0.125H when b < 0.46L and at H4 = 0.25H when b > 0.46L, respectively. It is interesting to note that in the case of upper bounds of the buckling load, the load distribution is proportional to individual subbeam’s flexural rigidities as the sub-beams 2 and 3 deflect together. Whereas for the case of the lower bounds of the buckling load, the load distribution is proportional to individual sub-beam’s axial stiffnesses as sub-beams 2 and 3 deflect independently and a possible local, global and/or mixed delamination buckling occurs [6,17,23]. Fig. 9 shows the effect of thicknesswise locations of the delamination-I (H4 = 0.125H, 0.25H, 0.375H and 0.5H) on the lower and upper bounds of the non-dimensionalized buckling load ðPÞ for various lengths of the delamination-II when the delamination-II

1.0 0.8 4

0.6

1

Constrained model Free model 1. a=0.2L Constrained model 2. a=0.4L 3. a=0.6L Free model 4. a=0.2L 5. a=0.4L 6. a=0.6L

P/Pcr

2

0.4

3 5

0.2 0.0 0.2

H4=0.375H a b

6

H3=0.375H

0.3

0.4

0.5

0.6

0.7

0.8

0.9

Delamination-II length (b/L)

2163

M.S.R. Parlapalli et al. / Computers and Structures 86 (2008) 2155–2165

a

1.0 0.8

P/Pcr

0.6

1

2

3

Constrained model 1. H4=0.125H 2. H4=0.25H 3. H4=0.375H single hump 4. H4=0.500H single hump 5. H4=0.375H s-shape 6. H4=0.500H s-shape

4

7

5 6

98 10

Free model Constrained model

0.4 0.2

H4

0.0 0.2

Free model 7. H4=0.125H 8. H4=0.250H 9. H4=0.500H 10. H4=0.375H

a b

0.3

H3=0.125H

0.4

0.5

0.6

0.7

0.8

0.9

Delamination-II length (b/L)

b

II

1.0 0.5 0.0 0.0 -0.5

I

III 0.5

(ii) Single hump (a = 0.2L and b = 0.3L)

Normalized modeshape

Normalized modeshape

(i) S-shape (a = 0.2L and b = 0.3L)

1.0 7

IV -1.0 Normalized distance (x/L)

III

1.0

IV

II 0.5 I

V 0.0 0.0 0.5 1.0 Normalized distance (x/L)

Normalized modeshape

(iii) Single hump (a = 0.2L and b = 0.6L)

III

1.0

0.5

IV

II I

V 0.0 0.0 0.5 1.0 Normalized distance (x/L)

Fig. 8. (a) Influence of the thicknesswise locations of the delamination-I on P, H3 = 0.125H, a = 0.2L d1 = d2 = 0.0; (b) buckling mode shapes of ‘constrained-mode’ model for H4 = 0.5H.

0.75

1 2

0.60

P/Pcr

is located at H3 = 0.25H. The length of the delamination-I is kept as 0.2L and d1 = d2 = 0. It is observed that the lower and upper bounds of the buckling loads reduces as the length of the delamination-I increases for given locations of the delamination-II. When the delamination-II is located at H3 = 0.25H, the lower and upper bounds of the buckling load are of minimum in magnitude when the delamination-I (a/L = 0.2) is located at H4 = 0.375H and H4 = 0.5H. The maximum values of the lower buckling load are observed when the delamination-I is located at 0.125H. Whereas the maximum values of the upper bounds of the buckling load are observed when the delamination-I is located at 0.25H when b < 0.725L and at 0.125H when a > 0.725L, respectively. From these observations, we noted that the lower and upper bounds of the buckling load are strongly dependent on the thicknesswise locations as well the lengths of the two enveloped delaminations. The difference between the lower and upper bounds of the buckling load is dependent on the buckling mode shapes of the subbeams as discussed previously.

65 7

0.45

Free model Constrained model Constrained model 1. H4=0.250H 2. H4=0.125H 3. H4=0.375H 3 4. H4=0.500H 4

Free model 5. H4=0.125H 6. H4=0.250H 7. H4=0.500H 8. H4=0.375H

8

0.30 0.15 0.00 0.2

H4

a b

0.3

H3=0.25H

0.4

0.5

0.6

0.7

0.8

0.9

Delamination-II length (b/L) Fig. 9. Influence of the thicknesswise locations of the delamination-I on P, H3 = 0.25H, a = 0.2L d1 = d2 = 0.0.

2164

M.S.R. Parlapalli et al. / Computers and Structures 86 (2008) 2155–2165

1.0

Free model Constrained model

0.8

2

3

0.6

P/Pcr

Constrained model 1. a=0.2L, b=0.4L 2. a=0.3L, b=0.5L 3. a=0.4L, b=0.6L

1

a

0.4

d1

Free model 4. a=0.2L, b=0.4L 5. a=0.3L, b=0.5L 6. a=0.4L, b=0.6L

H4=0.125H

b

H3=0.125H

d2

0.2

4

5 6

ment is observed between the results obtained from the analytical model and the finite element analysis, and between the analytical model and the literature. Detailed parametric studies are conducted which show that the lower and upper bounds of the buckling load are strongly depend on the thicknesswise and spanwise locations and lengths of the two enveloped delaminations. As the lengths of the two enveloped delaminations increases, the lower and upper bounds of the buckling load decreases. In the case of the upper bounds of the buckling loads, for certain thicknesswise locations of the enveloped delaminations, s-shaped buckling mode dominates the buckling behavior than a single hump buckling mode. Whereas in the case of the lower bounds of the buckling load mixed and/or local buckling mode dominates the buckling phenomena.

0.0 -0.2

-0.1

0.0

0.1

0.2

Spanwise locations of delaminations(d1/L=-d2/L) Fig. 10. Influence of the spanwise locations of the enveloped delaminations on P, H3 = 0.125H = H4.

4.3.3. Effect of the spanwise locations of the two enveloped delaminations Fig. 10 shows the influence of the spanwise locations of the two enveloped delaminations (d1/L and d2/L) on the lower and upper bounds of the non-dimensionalized buckling load ðPÞ for various lengths of the enveloped delaminations. The two enveloped delaminations are located at H3 = 0.125H = H4. Based on the lengths of the two enveloped delaminations, the spanwise location of the delaminations is restricted to 0.275L as shown in Fig. 10. The spanwise location of the delamination is zero i.e., d1/L = 0 means that the delamination-I is located symmetrically at the center of the beam. Similarly, d2/L = 0 means that the delamination-II is located symmetrically at the center of the beam. From Fig. 10, it is observed that the spanwise locations of the two enveloped delaminations has a strong influence on the upper bounds of the buckling load for all the cases of the lengths of the two enveloped delaminations. However, a nominal effect on the lower bounds of the buckling load is observed. The upper bounds of the buckling load are maximum for a/L = 0.2 and b = 0.4L when the delamination-I is located at the center, for a/L = 0.3 and b = 0.5L when d1 = d2 is 0.0125L and for a/L = 0.4 and b = 0.6L when d1 = d2 is 0.025L. The upper and lower bounds of the buckling load are minimum when the two enveloped delaminations are located close to the ends of the beam that suggests that the effect of the clamped boundary condition becomes a significant factor. Similar observation was noted for vibration of composite beam with enveloped delaminations by Shu and Della [36]. From these observations, it is noted that besides the lengths of the delaminations, the locations (along the thicknesswise and spanwise) of the delaminations strongly influence the upper and lower bounds of the buckling load and a detailed investigation of these parameters are essential to understand the buckling phenomena. The present models can be further improved by using Timoshenko first shear deformation theory or zig-zag plate theories [37] instead of using the Euler–Bernoulli beam theory. Readers can find more details about the zig-zag plate theories for multilayered plates and shells in the review paper by Carrera [37]. 5. Conclusions By considering the effects of differential stretching and the bending–extension coupling, the lower and upper bounds of the buckling load of a composite beam with two enveloped delaminations are obtained by developing analytical models. A good agree-

References [1] Reifsnider KL, Talug A. Analysis of fatigue damage in composite laminates. Int J Fatigue 1980;3(1):3–11. [2] Adams DF, Zimmerman RS, Odom EM. Frequency and load ratio effects on critical strain energy release rate Gc thresholds of graphite/epoxy composites. Philadelphia: ASTM; 1987. STP937. [3] O’Brien , Zimmerman RS, Odom EM. Mixed-mode strain energy release rate effects on edge delamination of composites. ASTM ; 1984. STP836, pp. 125–42. [4] Chai H, Babcock CD, Knauss WB. One-dimensional modeling of failure in laminated plates by delamination buckling. Int J Solids Struct 1981;17(11): 1069–83. [5] Bottega WJ, Maewal A. Delamination buckling and growth in laminates. J Appl Mech 1983;50:184–9. [6] Simitses GJ, Sallam S, Yin WL. Effect of delamination of axially loaded homogeneous laminated plates. AIAA J 1985:1437–44. [7] Yin WL, Sallam SN, Simitses GJ. Ultimate axial load capacity of a delaminated beam-plate. AIAA J 1986;24(1):123–8. [8] Kassapoglou C. Buckling, post-buckling and failure of elliptical delaminations in laminates under compression. Compos Struct 1988;9:139–59. [9] Jane WL, Yin WL. Refined buckling and post-buckling analysis of twodimensional delamination – I. Analysis and validation. Int J Solids Struct 1992;29(5):591–610. [10] Shu D, Mai YW. Effect of stitching on interlaminar delamination extension in composite laminates. Compos Sci Technol 1993;49(2):165–71. [11] Kutlu Z, Chang FK. Modeling compression failure of laminated composites containing multiple delaminations. J Compos Mater 1992;26(3):350–85. [12] Lee J, Gurdal Z, Griffin OH. Layer-wise approach for the bifurcation problem in laminated composites with delaminations. AIAA J 1993;31(2):331–8. [13] Lim YB, Parsons ID. The linearized buckling analysis of a composite beam with multiple delaminations. Int J Solids Struct 1993;30(22):3085–99. [14] Adan M, Sheinman I, Altus E. Buckling of multiply delaminated beams. J Compos Mater 1994;28(1):77–90. [15] Kyoung WM, Kim CG, Hong CS. Modeling of composite laminates with multiple delaminations under compressive loading. J Compos Mater 1998;32(10):951–68. [16] Huang H, Kardomateas GA. Buckling of orthotropic beam-plates with multiple central delaminations. Int J Solids Struct 1998;35(13):1355–62. [17] Shu D. Buckling of multiple delaminated beams. Int J Solids Struct 1998;35(13):1451–65. [18] Lachaud F, Lorrain B, Michel L, Barriol B. Experimental and numerical study of delamination caused by local buckling of thermoplastic and thermoset composites. Compos Sci Technol 1998;58:727–33. [19] Hu N. Buckling analysis of delaminated laminates with consideration of contact in buckling mode. Int J Numer Methods Eng 1999;44:1457–79. [20] Gu H, Chattopadhyay A. An experimental investigation of delamination buckling and postbuckling of composite laminates. Compos Sci Technol 1999;59:903–10. [21] Sekine H, Hu N, Kouchakzadeh MA. Buckling analysis of elliptically delaminated composite laminates with consideration of partial closure of delamination. J Compos Mater 2000;34(7):551–76. [22] Hwang SF, Mao CP. Failure of delaminated interply hybrid composite plates under compression. Compos Sci Technol 2001;61:1513–27. [23] MSRao P, Shu D. Buckling of two-layer beams with an asymmetric delamination. Eng Struct 2004;26(5):651–8. [24] Rao PMS, Wenge T, Shu D. Buckling analysis of tri-layer beams with enveloped delaminations. Compos Part B: Eng 2005;36(1):33–9. [25] Miklashevich A. Delamination of composites along the interface as buckling failure of the stressed layer. Mech Compos Mater 2004;40(4):279–86. [26] Aniello R, Francesco S, Pierluigi P. Influence of contact phenomena on embedded delaminations growth in composites. AIAA J 2003;41(5):933–40. [27] Yin WL, Jane KC. Vibration of a delaminated beam-plate relative to buckled states. J Sound Vibr 1992;156(1):125–40. [28] Luo H, Hanagud S. Dynamics of delaminated beam. Int J Solids Struct 2000;37:1501–19.

M.S.R. Parlapalli et al. / Computers and Structures 86 (2008) 2155–2165 [29] Ostachowicz W, Krawczuk M, Cartmell M, Gilchrist M. Wave propagation in delaminated beam. Comput Struct 2004;82:475–83. [30] Parlapalli M, Shu D. Buckling analysis of two layer delaminated beams with bridging. Eur J Mech – A/Solids 2006;25(5):834–53. [31] Ting-Shun WJ, Chien-Chang L. Engineering analysis of buckling of delaminated beam plates. Compos Struct 1996;34:397–407. [32] Shu D, Fan H. Free vibration of a bimaterial split beam. Compos Part B: Eng 1996;27B:79–84. [33] Jones RM. Mechanics of composite materials. 2nd ed. Philadelphia: Taylor and Francis; 1999.

2165

[34] Parlapalli M, Shu D. Effects of bridging on buckling analysis of tri-layer beams with separated delaminations. AIAA J 2006;44(9):2108–17. [35] Yap CW, Chai GB. Analytical and numerical studies on the buckling of delaminated composite beams. Compos Struct 2007;80:307–19. [36] Shu D, Della CN. Vibrations of multiple delaminated beams. Compos Struct 2004;64:467–77. [37] Carrera E. Historical review of zig-zag theories for multilayered plates and shells. Appl Mech Rev 2003;56(3):287–303.