Compressive buckling of laminates with an embedded delamination

Composites Science and Technology 59 (1999) 1247±1260

Compressive buckling of laminates with an embedded delamination Ning Hu a,*, Hisao Fukunaga b, Hideki Sekine b, Kouchakzadeh Mohammad Ali b a

Department of Engineering Mechanics, Tsinghua University, Beijing 100084, People's Republic of China b Department of Aeronautics and Space Engineering, Tohoku University, Sendai 980-77, Japan Received 8 January 1998; received in revised form 10 August 1998; accepted 1 October 1998

Abstract In this paper, a buckling analysis of laminates with an embedded delamination has been conducted by employing a ®nite-element method based on the Mindlin plate theory. An e€ective solution method has been put forward to deal with the contact problem in the buckling mode. In this method, an iterative updating process incorporating ®rst-order sensitivity analysis and quadratic programming technique has been proposed for computing ®ctitious forces in the contact area, which can be used to calculate sti€ness parameters of some arti®cial springs. Finally, the penetration between two delaminated layers in the buckling mode can be e€ectively prevented by inserting these arti®cial springs into the overlapped area in the buckling mode. Numerical examples indicated that this method is very ecient for solving the contact problem in buckling analysis. Its accuracy and convergence speed are very high. The numerical results demonstrate that the contact analysis is very important for buckling analysis in some cases. The e€ects of various sizes, shapes and positions of the delamination and the ®ber angle of the sublaminate on the buckling load have also been investigated. # 1999 Elsevier Science Ltd. All rights reserved. Keywords: C. Buckling; C. Delamination; C. Laminates; C. FEA; Contact

1. Introduction Delamination is one of the most common failure modes in composite laminates, and may be formed as a consequence of various impact events, poor fabrication processes and fatigue. As is well known, the compressive strengths of structures made from laminated composite materials may be severely reduced by the presence of this delamination damage. The region bounded by a delamination and the laminate free surface is liable to buckle locally (commonly referred to as delamination buckling) under compressive in-plane loads, thereby creating conditions conducive of delamination growth, and consequent global failure of the structure. This laminate failure mode is a critical design concern for many composite structures. From the ®rst research on this problem by Chai et al. [1] in 1981, extensive research work has been conducted in this ®eld. A common problem occurring in the buckling analysis of delaminated plates is layer overlap, which means that delaminated layers may penetrate each other. The contact phenomena may appear for a number of reasons. First, the generally non-symmetric nature of the * Corresponding author.

sublaminate (thin delaminated layer) layup may tend to make it de¯ect in an asymmetric mode which is, however, prevented by the balance of the parent substrate (the remaining part of the plate). Second, a sublaminate with a small out-of-plane thickness ratio may be liable to de¯ect in a local mode, and at the same time the parent substrate deforms in a global one. In most of the previous researches on buckling analysis of delaminated laminates, widely used models which are essentially one-dimensional have been investigated [1±9]. Except for that indicated by Chai et al. [1] where contact occurs under conditions of some bending rigidity in the parent substrate and multiple contacts between di€erent layers caused by multiple delaminations shown in the literature, it was found that usually no contact occurs in most one-dimensional models which may include some multiple delamination cases [7± 9]. It may be explained that the governing partial differential equations in one-dimensional cases are reduced to ordinary di€erential equations which may imply the eigenmode corresponding to the lowest eigenvalue to be of one sign. Hence, the contact analyses usually have not been considered in most one-dimensional cases. For more complex situations, e.g. two- or three-dimensional problems, buckling analyses have been conducted [10±15].

0266-3538/99/$ - see front matter # 1999 Elsevier Science Ltd. All rights reserved. PII: S0266 -3 538(98)00166 -3

1248

N. Hu et al. / Composites Science and Technology 59 (1999) 1247±1260

Most of these analyses considered the e€ects of contact conditions among the delaminated layers, since it was found that contact frequently occurs in two- or threedimensional models. For example, Peck and Springer [11] and Whitcomb [13] used the concept of a ®ctitious force to deal with the contact problem. However, their solution methods [11,13] can only solve contact problems arising from the equilibrium problem under practical compressive loads, in this case the critical buckling load is calculated from the load/strain history. Actually, the critical buckling load can be obtained e€ectively from the solution of an eigenvalue problem. However, many previous researches for obtaining the buckling load of delaminated laminates were focused on the linear eigenvalue analysis where no contact is considered [2,4±9]. This linear eigenvalue analysis may result in some physically impossible or inadmissible buckling modes. Especially when the sublaminate is very thin and ¯exible, or with non-symmetric layup, its local or asymmetric buckling mode may occur in the delamination area to cause the penetration between two delaminated layers. Unlike usual contact problems in the equilibrium problem under the known applied load [11,13], there are two unknowns desired in the present eigenvalue problem, i.e. buckling load and buckling mode. Therefore, the solution of this contact problem becomes more dicult. Recently, the work conducted by Suemasu et al. [15] was concerned with contact conditions in the buckling mode of a free-edge delamination problem. They employed a kind of imaginary spring, which can provide the resistance force and moment to prevent the penetration between two delaminated layers. This method may be classi®ed into the well-known penalty function approach for dealing with contact problems. However, no detailed algorithm to estimate the sti€ness parameters of the spring has been proposed [15]. In fact, there is no clear rule to choose these parameters in the ordinary penalty function method and it depends on particular problems considered. Moreover, an e€ective algorithm for dealing with contact problems in the buckling analysis of delaminated laminates was proposed by the ®rst author of this paper, based on the concept of a ®ctitious force and an arti®cial spring [16]. Also, the sensitivity analysis and quadratic programming technique were incorporated into this iterative algorithm. However, the prime contents in Ref. [16] were concentrated on the investigation of some important properties of this algorithm and veri®cation of this algorithm through some examples in the literature. Hence, there is insucient evidence: (a) E€ect of contact analysis on the buckling analysis has not been investigated within a wide range of case presentations. (b) The in¯uence of various sizes, shapes and positions of embedded delaminations and the ®ber angle of the sublaminate on the buckling load has not been studied in detail using this method.

Our work is designed to move the present state of knowledge forward by producing new evidence. This paper is organized as follows. Section 2Section 3Section 4 will present the overall theories used in numerical analyses. This will keep the paper structured. Section 5 addresses some important veri®cation examples. Section 6 illustrates and discusses the results of our research and Section 7 presents the conclusions. 2. Governing equations Firstly, for the composite laminated plate, eightnoded isoparametric Mindlin plate elements including the in¯uence of the shear deformation are used. The shear correction factor of the Mindlin plate theory used is 5/6. Due to non-symmetric layup in delaminated layers caused by delamination, the in-plane deformation should be taken into account. Five degrees of freedom, i.e. {ui, vi, wi, fxi, fyi} on each node are used. Furthermore, the local deformations in the vicinities of the delamination front are not considered, because they are thought to have little in¯uence on the global behavior of the plate. The plate can be divided into three parts: sublaminate, parent substrate and the undelaminated portion. As shown in Fig. 1, the continuity conditions of displacements on the delamination front are: wi ˆ w0 ;

…1a†

xi ˆ x0 ;

…1b†

yi ˆ y0 ;

…1c†

ui ˆ u0 ‡ zi x0 ;

…1d†

i ˆ 0 ‡ zi y0 ;

…1e†

Fig. 1. Continuity conditions of the displacements at the delamination front.

N. Hu et al. / Composites Science and Technology 59 (1999) 1247±1260

1249

where i denoting two delaminated layers is from 1 to 2, and 0 represents the mid-plane of the undelaminated portion. The governing equation of the buckling analysis from the ®nite-element method can be expressed by the following equation, which can be obtained from the wellknown Tre€tz criterion [4]: Ku ˆ l K u;

…2†

where K is the global elastic sti€ness matrix of the laminated plate, u the buckling mode, ls the critical buckling load, Ks the global geometric sti€ness matrix of the plate. The matrices K and Ks can be obtained from the ordinary ®nite-element method. Moreover, the simultaneous iteration algorithm proposed by Corr and Jennings [17] is adopted to analyze this eigenvalue problem.

It is expected that this ®ctitious force vector F in Eq. (3) can overcome the penetration between two delaminated layers. Usually, the contact solution of the above equation can be transferred into the following more general optimization model of a complementarity problem: … minimize G ˆ Fus dc …4a† c

3. Contact analysis As explained above, in some cases the sublaminate may penetrate into the substrate to lead to the physically inadmissible buckling modes. Hence, it is important to deal with this contact problem. However, as aforementioned the diculty arises from the two unknowns in the eigenvalue problem of buckling analysis, i.e. buckling load and buckling mode. To simplify the analysis, the contact problem is considered by ignoring the friction and sliding between the two delaminated layers. Therefore, the contact conditions among possible contact pairs only depend on the relative lateral de¯ections of the two delaminated layers in the buckling mode. For the contact solutions in general equilibrium problems with known external loads, most of the present existing methods can be classi®ed into three major categories which include the Lagrange multiplier method [18], Penalty function method [19] and the Mathematical programming method [20,21]. Naturally, these methods cannot be directly employed to solve the contact problems in the eigenvalue analysis due to two unknowns. In this section, we will describe our method for solving this special contact problem.

subject to F50

…4b†

us 50;

…4c†

where C is the contact boundary and us the gap vector of two delaminated layers in the lateral direction. There are two unknowns F and us in Eqs. (4a)±(4c). Hence, it is very dicult to solve it directly. Here, our objective can be changed to minimize the summation of gaps in the contact area and then Eqs. (4a)±(4c) can be modi®ed into the following least-square form: …5a† minimize G ˆk us k2 subject to F50:

As shown in Fig. 2, a ®ctitious force is assumed at the ith contact pair which possesses the node number iu on the upper delaminated layer, il on the lower delaminated layer. Therefore, a ®ctitious force vector F which includes the ®ctitious forces at all possible contact pairs can be formed in the contact area and Eq. (2) can be modi®ed as follows: …3†

…5b†

For enforcing the contact complementarity condition in Eq. (4a), the components in us vector should include those gaps at the following two kinds of contact pairs, which are called the contacting pairs: uis 40; uis 50;

…6a† and

Fi > 0:

…6b†

In case of no rigid-body displacement, the gap function us can be expressed as: us ˆ us …F; p†;

3.1. QP optimum model for contact problem

Ku ˆ l K u ‡ F:

Fig. 2. Fictitious force between two overlapped layers.

…7†

where p is the position vector of contacting pairs. Moreover, the gap function in iteration (k + 1) can be expanded into the form of a Taylor series with respect to contact force vector F:  uks ‡ uk‡1 s

@uks
Fk ˆ uks ‡ Sk
Fk @Fk

…8†

into Eq. (5a)Eq. (5b) leads to the Substitution of uk‡1 s following equations:

1250

N. Hu et al. / Composites Science and Technology 59 (1999) 1247±1260

minimize

1 T T G ˆ
Fk Hk
Fk ‡ Ck
Fk ‡ Dk 2

…9a†

subject to

Fk ‡
Fk 50;

…9b†

where T

…9c†

T

…9d†

Hk ˆ S k S k ; Ck ˆ Sk uks ; 1 T Dk ˆ uks uks : 2

…9e†

In the above optimization model, there is only one unknown vector F. The above quadratic programming problem can be solved by employing some ecient algorithms, e.g. Lemke algorithm. Here, the method proposed by Goldfarb and Idnani [22] is adopted directly. After obtaining the increment of the ®ctitious force vector from the above computation, we can update F vector in the next iteration as follows: Fk‡1 ˆ Fk ‡
Fk :

…10†

Due to the constraint condition in Eq. (9b), the components in Fk + 1 vector in Eq. (10) can be guaranteed to be positive, which denotes the resistance forces within the contacting pairs as shown in Fig. 2. 3.2. Arti®cial spring

k‡1

i ;

i ˆ 1; . . . nc ;

k‡1

where ukis … i ˆ 1; . . . nc †is the gap at contacting pair i in iteration k. Therefore, Eq. (12) can be changed into the following form: k‡1 Fk‡1 …13† Kpi ˆ ik ;  where is a parameter to accelerate the convergence speed. k‡1 Inspection of Eq. (13) reveals that Kpi … i ˆ 1; . . . nc † at di€erent contacting pairs are determined by the values of the corresponding components in F k+1 vector with the same ratio, since and ak are the same for all contacting pairs. Furthermore, it can be found that k‡1 is always positive, hence it possesses correct phyKpi sical meaning. By employing these sti€ness parameters, the element sti€ness matrix Ki of the ith arti®cial spring can be constructed [16]. Finally, the sti€ness matrix of the system can be updated by: nc X k‡1 Ki ; …14† Kk‡1 ˆ K0 ‡ iˆ1

After obtaining the ®ctitious force vector from the QP computation, the next step is to ®nd a way to solve Eq. (3). Obviously, Eq. (3) is di€erent from the usual eigenvalue problem due to the presence of F vector. Therefore, we cannot solve Eq. (3) directly by employing all the existing algorithms for analyzing eigenvalue problems. Note that F vector is used to overcome the negative gaps among all contacting pairs. Naturally, it can be thought that the components in F vector are related to some arti®cial springs through the gaps as shown in Fig. 3, which leads to the following equation: ˆ Kpi Fk‡1 i

of the gap at the same contacting pair after inserting the spring, nc the number of all contacting pairs in iteration k. However, we do not know anything about ai before k‡1 is obtained. For simplicity, a parameter (average Kpi gap in iteration k) is chosen for all contacting pairs: Pnc kis j ; …12† k ˆ iˆ1 nc

…11†

is the sti€ness parameter of the arti®cial where Kpi spring at contacting pair i in iteration k, ai the variation

0

where K is the original sti€ness matrix of the system. In the above iterative procedure, there are two approximations made in Eqs. (8) and (12). Actually, the convergence cannot be guaranteed theoretically. However, for the ®rst approximation, the basic idea is almost the same as that of the Newton±Raphson method for general non-linear problems. For the second one, the computational error in the sti€nesses of arti®cial springs may be actually caused in the ®rst several iterations. However, as the iteration proceeds, the gaps at di€erent contacting pairs tend to be the same small value, the average gap can be used to calculate the sti€nesses of the arti®cial springs at di€erent contacting pairs. Furthermore, this error can also be removed partially by employing the following step to modify the ®ctitious forces. 3.3. Modi®cation of ®ctitious force vector

Fig. 3. Arti®cial spring between two overlapped layers.

By employing the updated sti€ness matrix, the eigenvalue problem in Eq. (2) can be resolved and a new gap vector can be obtained. The ®ctitious force vector should be modi®ed at this step, since the e€ect of the ®ctitious force vector is considered by approximately using the arti®cial springs from Eq. (13) where ak signi®es the average gap at the previous iteration. The ®ctitious force vector is modi®ed as follows:

N. Hu et al. / Composites Science and Technology 59 (1999) 1247±1260

Fk‡1 ˆ ÿKPi i

k‡1

uk‡1 is ;

i ˆ 1; . . . nc

…15†

is the new gap at contacting pair i. where uk‡1 is From Eq. (15), it can be found that the ®ctitious force will become negative when the new gap is positive. Therefore, the judgment for contacting pairs described in Eq. (6) should be changed into the following alternative conditions: uis < 0 uis 50

…16a† and

Fi < 0

…16b†

The negative ®ctitious force is unrealistic, since it denotes the tension within a contacting pair. It can be removed in the optimum step from Eq. (9b). 3.4. Sensitivity analysis The sensitivity matrix S can be calculated by employing the forward di€erence method. Assume the k di€erence step of the ®ctitious force vector as
iF for contacting pair i, Eq. (3) can be written as: k

Kk u ˆ lk K u ‡
iF ;

…17† k

where the di€erence step vector
iF can be expressed as:  k …18†
iF ˆ 0


0 ik 0


0 ÿ ik 0


0 ; where the positions of ik and ÿik are determined by the equation numbers of the degrees of freedom of de¯ections of nodes iu and il on upper and lower layers shown in Fig. 2, respectively. De®ne: Xk ˆ lk K uk ;

…19†

where uk is the buckling mode in the kth optimum iteration step, and ik is suggested as follows: ÿ
…20† ik ˆ km j ‡kn j =100; where Xkm and Xkn are the components in the vector Xk. Their positions are identical to those of ik and ÿik , respectively. It can be found that …Kk u ÿ lk K u† is rank-de®cient. Hence it is very dicult to obtain u by employing some direct solution methods to calculate the gap vector under the di€erence step. An iterative procedure is suggested here:

1251

The convergence check on this iteration is performed using the following criterion:  ˆ j…
uj‡1 ÿ
uj †=
uj j4";

…22†

where
uj‡1 denotes the di€erence between two iterations: p …23†
uj‡1 ˆ k uj‡1 ÿ uj k2 : The convergence tolerance in Eq. (22) is assumed as 0.01 in our computation. It should be noted that the convergence criterion in Eq. (22) is di€erent from the usual convergence criterion, which may use
uj‡1 4" directly. However, in our problem it is unsuitable to do k so, since the di€erence step vector
iF of ®ctitious force vector is comparatively smaller by observing Eq. (20) and consequently
uj‡1 can reach a very small value quickly and easily before suciently accurate data is obtained. Therefore, in this case, it may be very hard to select the convergence tolerance. To avoid this problem, we utilize Eq. (22) to perform the convergence check, which is comparatively more strict. For the usual problems, the tolerance 0.01 may not be small enough. However, for the convergence criterion in Eq. (22), it was found from our numerical experiences [16] that this value is sucient for obtaining the accurate sensitivity in general. After convergence of the above sub-iteration, we can get the gap vector for contacting pair i under the ®ctitious force di€erence step. Finally, the sensitivity matrix S can be calculated from the gap vectors for all contacting pairs under the di€erence step, original gap vector obtained in Eq. (2) and the di€erence step. Detailed description of this process and the computational amount of this method can be referred to Ref. [16]. 3.5. Iteration strategy By summarizing the algorithm described previously, the sequence of calculations in iteration k can be stated in the following manner:

…21†

1. Decompose the sti€ness matrix Kk and perform the eigenvalue analysis. 2. Calculate the corresponding gap vector. 3. Check for the convergence in the optimum updating iterations. For the sake of convenience, the following convergence criterion is used:

4; …24†

where j is the iteration number for obtaining the sensitivity before convergence. Since the sti€ness matrix Kk has been already decomposed at the previous solution of the eigenvalue problem, the computation amount here is only the backsubstitution of the system equation which is usually small.

where is the smallest absolute gap value within all contacting pairs and  the convergence tolerance. 4. When k > 1, modify the ®ctitious force according to Eq. (15) and determine the contact states in possible contact pairs from Eqs. (16a) and (16b). 5. Use the decomposed sti€ness matrix in step 1 and the gap vector in step 2 to calculate the sensitivity matrix.

k

Kk uj‡1 ˆ lk K uj ‡
iF ;

1252

N. Hu et al. / Composites Science and Technology 59 (1999) 1247±1260

Fig. 4. Geometry for a through-the-width delaminated plate.

6. Evaluate the matrices in Eqs. (9c)±(9e), and solve the QP problem to get the increment vector of the ®ctitious force
Fk. 7. Update the ®ctitious force vector and sti€ness matrix of the system from Eqs. (10) and (14), then go back to step 1. In general, the procedure described above can be thought of as a hybrid-version of the Mathematical programming method and the Penalty function method. 4. Numerical test problems A FORTRAN program, based on the theories described above, has been developed and veri®ed by some one-dimensional [2,3,6,9] and two-dimensional examples [14,15] which also include the experimental results [14]. A good consensus was obtained. In this section, some selected examples are illustrated and discussed. Other detailed validations of our method can also be found in Ref. [16]. 4.1. Through-the-width delamination Firstly, a specially orthotropic composite laminate containing one centrally located mid-plane delamina-

tion shown in Fig. 4 is considered, which was studied by Simitses et al., [2], Chen [6] and Lee et al. [9] Material properties are listed as follows: E11 ˆ 1:81  102 GPa;

E22 ˆ 1:03

1

 10 GPa; G12 ˆ G13 ˆ G23 ˆ 7:17 GPa; ˆ 0:28: Plate speci®cations are taken as:

12 …25†

h=L ˆ 1=400;

…26a†

h ˆ H=2:

…26b†

The clamped boundary conditions at the two ends are applied. Comparison of the present results and those in Refs. [2,6,9] is shown in Table 1 for various delamination lengths. Normalized buckling loads are obtained by dividing the buckling loads of the delaminated plate to that of a plate without delamination. From Table 1, it can be found that all results of the present method are in excellent agreement with the other three solutions [2,6,9]. There is one interesting point about the results. Buckling mode is global symmetric for most cases of a/ L, however it changes into a global antisymmetric mode when a/L=0.4. This mode cannot be obtained if symmetry assumptions are applied along the axial direction as done in Ref. [6]. The symmetry model in Ref. [6] will lead to a little higher result (shown in Bold Italic in Table 1). The present result con®rms this phenomenon.

Table 1 Comparison of normalized buckling loads for various lengths of a through-the-width mid-plane delamination a/L

Ref. [2]

Ref. [6]

Ref. [9]

Present

0.3 0.4 0.5 0.6 0.7 0.8 0.9

0.9638 0.8481 0.6896 0.5411 0.4310 0.3514 0.2923

0.9638 0.8561 0.6896 0.5411 0.4310 0.3514 0.2933

0.9639 0.8482 0.6898 0.5413 0.4311 0.3515 0.2934

0.9606 0.8445 0.6860 0.5384 0.4288 0.3497 0.2919

N. Hu et al. / Composites Science and Technology 59 (1999) 1247±1260

Fig. 5. Comparison of buckling loads by Kardomateas and Schmueser [3] and the present method for a through-the-width delaminated plate with various delamination lengths and depths.

1253

Fig. 7. Comparison of buckling loads by Lee et al. [9] and present FEM solutions for three delaminations.

The boundary conditions are identical to those in the above example. The present results and those obtained by Kardomateas and Schmueser [3] from the theory including transverse shear deformation are displayed in Fig. 5. A good consensus is exhibited between the two kinds of results.

(h/H) on the buckling load of laminates. Three delaminations are considered. Two delaminations symmetrically located with respect to the mid-plane and the last one is at the mid-plane. Non-dimensional buckling load vs. non-dimensional delamination length is shown in Fig. 7. For clear observation, only the results in Ref. [9] and the present results are plotted. It can be found that the present results are in good agreement with those by Lee et al. [9] Also, they agree well with the results by Wang et al. [7] except for the case of three delaminations with h/ H=0.5, the reason can be found in Ref. [9]. In the above multiple delamination example, it should be noted that none of the critical buckling load results show physically inadmissible buckling mode. This phenomenon was also found and pointed out by Lee et al. [9].

4.2. Multiple delaminations

4.3. Embedded delamination

For comparison purposes, a multiple delamination problem in Fig. 6 solved by Wang et al. [7] and Lee et al. [9] is further examined. The material used in the study is SMC-R50, which possesses the following material properties:

A laminate with an elliptical embedded delamination [14] as shown in Fig. 8 is analyzed. This plate consists of eight plies with stacking sequence (+/ÿ)2s. Material properties and speci®cations of the plate and delamination are

Furthermore, to check the validation of the in¯uence from shear deformation in thicker sublaminate, a through-the-width delamination which can be located anywhere through the thickness of the plate [3] is considered. Material properties and plate speci®cations are E11 ˆ 7:0  101 GPa; E22 ˆ 4:5 GPa; ˆ 2:5 GPa; 12 ˆ 0:35; L=H ˆ 1=200:

E11 ˆ 1:090  101 GPa;

G12 …27a† …27b†

E22

ˆ 7:586 GPa; G12 ˆ 2:483 GPa; ˆ 0:22; 13 ˆ 0:31:

12 …28†

This example demonstrates the e€ects of non-dimensional length (a/L) and through-the-thickness location

Fig. 6. Three symmetrically located delaminations through the thickness.

Fig. 8. Model of a laminated plate with an elliptical delamination by Yeh and Tan [14].

1254

N. Hu et al. / Composites Science and Technology 59 (1999) 1247±1260

4.4. Free-edge delamination For veri®cation of the present contact algorithm, a plain woven carbon ®ber-reinforced composite plate with a free-edge delamination shown in Fig. 10 is considered [15]. Material properties and the speci®cations of the plate and delamination are: E11 ˆ E22 ˆ 6:17  101 GPa; ˆ 5:17 GPa; 12 ˆ 0:046;

G12

L ˆ 12:0 cm;

H ˆ 0:348 cm;

ˆ 1=8: Fig. 9. Comparison of buckling loads by Yeh and Tan [14] and present FEM solutions for an embedded delamination.

E11 ˆ 9:63  101 GPa; ˆ 3:2 GPa; 2a=L ˆ 0:61; ˆ 1:72;

E22 ˆ 6:7 GPa;

G12

12 ˆ 0:27; L ˆ 6:0 cm;

H ˆ 0:096 cm;

…29a† W ˆ 4:0 cm; h=H ˆ 3=8:

a=b …29b†

The same boundary conditions as those in the above examples are used. Yeh and Tan [14] carried out the experiment of this model and compared it with the FEM solutions. In their FEM approach, a degenerated shell element with large displacement and rotation was employed and the buckling load was obtained from the load/strain history. Also, the updated Lagrangian formulations were used to deal with the large deformation. Comparison of the present results and those obtained by Yeh and Tan [14] for various ®ber orientations is shown in Fig. 9. It can be found that the present results coincide with the numerical and experimental results obtained by Yeh and Tan [14] with permissible accuracy. This comparison validates the e€ectiveness of the present method when sublaminates are not with single ®ber orientation.

Fig. 10. Geometry for a composite plate with a free-edge delamination by Suemasu et al. [15].

W ˆ 5:5 cm;

…30a† h=H …30b†

Three edges including two edges with applied loads and the undelaminated edge are simply supported. In our computation, the convergence tolerance  in Eq. (24) is de®ned as 10ÿ4 which is same as that used in Ref. [15]. The parameter in Eq. (13) is arbitrarily chosen as 100. Comparison of the present results and those obtained by Suemasu et al. [15] is illustrated in Fig. 11. The buckling load of the delaminated plate is normalized using the buckling load of the intact plate. Two kinds of results by Suemasu et al., [15] which denote the physically admissible solutions with contact analysis and physically inadmissible ones without consideration of contact, are plotted. From Fig. 11, it can be seen that the present results agree very well with the physically admissible solutions of Suemasu et al. [15] Also, it can be found that the buckling load increases due to the introduction of contact analysis. This increasing tendency is very obvious when the delamination width is large. Regarding the convergence speed of contact analysis, for this free-edge delamination problem, usually the present technique needs only four or ®ve iterations to obtain the converged results, whereas it needs about ten iterations in Ref. [15]. However, direct comparison of the two methods is quite dicult owing to completely di€erent theoretical backgrounds. The detailed description

Fig. 11. Comparison of results by Suemasu et al. [15] and present FEM solutions for a free-edge delamination.

N. Hu et al. / Composites Science and Technology 59 (1999) 1247±1260

1255

Fig. 12. Geometry for a laminated plate with an elliptical delamination at the center.

of the computational amount of the present method was already given in Ref. [16]. Also, detailed investigation of the contact analysis in the buckling analysis will be addressed in the following section for an embedded delamination problem. 5. Model for an embedded delamination problem For investigation of the laminates with an embedded delamination, a standard plate is de®ned in this section. Our analytical models in Section 6 will be constructed based on this standard plate by introducing some small modi®cations. A simply supported ten-plied [90/0/90/0/0]s laminated square plate shown in Fig. 12 is analyzed. The material constants of single lamina are identical to those in Eq. (25). Dimensions of this square plate are given as: L=20.0 cm and H=0.3 cm. An elliptical delamination existing between the 1st and the 2nd laminae is assumed to be located at the center of the laminate. Two radii of the delamination are de®ned as: a=b=4.5 cm. As stated in Section 4.1, though there is symmetry in this problem, at least a half part of this plate should be modeled in our numerical analyses by noting that there may be an anti-symmetric buckling mode in analysis. The FEM mesh is illustrated in Fig. 13, which includes 518 eight-noded isoparametric plate elements and 1613 nodes. Convergence of this mesh has been checked by more re®ned meshes. It was found that this mesh is enough to tackle this problem. 6. Results and discussions In this section, we will investigate the e€ect of the contact analysis on the buckling load and discuss the

Fig. 13. FEM mesh for a square composite plate with an embedded delamination.

elimination process of overlap ®rst. Also, the e€ects of various sizes, shapes and positions of the delamination and the ®ber angle of the sublaminate on the buckling load are studied in detail. 6.1. Contact analysis The standard plate as described above is analyzed to show the elimination process of the penetration between two delaminated layers. Two parameters and  are identical to those in Section 4.4. In the ®rst step of computation, the sectional de¯ection distribution along P±Q line in the buckling mode is displayed in Fig. 14(a). It can be found that there exists obvious overlap between the upper and lower layers. The buckling mode of the upper thin layer is a complex local symmetric shape, but the thick lower layer has no obvious deformation. From Fig. 14(b)±(d), it can be observed that the

1256

N. Hu et al. / Composites Science and Technology 59 (1999) 1247±1260

Fig. 14. De¯ection distribution along P±Q line in the buckling mode within the iteration process.

overlap between two delaminated layers is eliminated as contact analysis proceeds. At the same time, the deformation of the lower layer becomes signi®cant. Acceptable convergence is obtained in the fourth iteration. Hence, the convergence speed of this method is very high. By comparing the shapes of the initial buckling mode and the ®nal buckling mode, it can be found that the buckling mode changes drastically. It changes from an overlapped local mode to a global one with the small local deformation of the upper thin delaminated layer. It should be mentioned that the buckling mode presented in the individual iteration in Fig. 14 does not represent the buckling history. It is, however, introduced by an iterative procedure designed to lead to the correct solution at the buckling point. In other words, the intermediate results presented in Fig. 14 represent only the iteration history, which demonstrates the e€ectiveness of the proposed method. Furthermore, it should be pointed out that the overlapped buckling mode (see Fig. 14(a)) will not be present without imposed perturbation. Moreover, if contact buckling mode appears as shown in Fig. 14(a) due to perturbation, it may be usually unstable as reported by Suemasu

et al. [15]. They found that the buckling loads obtained from experiments in this case are usually scattered. Actually, in the post-buckling process, the upper layer will de®nitely move upwards and then there will be no contact.

Fig. 15. Convergence phenomenon of the buckling load within the iteration process.

N. Hu et al. / Composites Science and Technology 59 (1999) 1247±1260

The iterative history of the buckling load is shown in Fig. 15. As can be expected, the buckling load increases as the iteration proceeds. The value of the buckling load is from 1023.92 (N/cm) to 1252.17 (N/cm) with 22% increase. Our numerical experience has shown [16] that the convergence speed can be raised if a relatively large value of is chosen. Moreover, the stability and accuracy of computation are not a€ected. It was also revealed [16] that the contact analysis can lead to considerable increase of the buckling load only when the buckling mode changes drastically within the contact iteration process. Usually, the drastic changes in the buckling mode are liable to occur when there is signi®cant overlap in the initial analysis. In contrast, the in¯uence of the contact analysis on the buckling load is very small when there is no obvious overlap in the initial analysis. 6.2. E€ect of delamination size on buckling load

1257

mode and local buckling mode. The representative signs of three buckling modes are put on the appropriate portions of the graph for a better understanding of the deformation behavior of the plate. From Fig. 16, it can be found that the plate buckles in a global mode for small delamination diameters. In this case, the delamination has very little in¯uence on the buckling load of the plate. Moreover, no serious overlap occurs and the e€ect of contact analysis is negligible. With the increase of the delamination size, the upper thin delaminated layer tends to buckle independently which results in mixed modes. In this case, the e€ect of overlap turns out to be more visible. Also, the buckling load starts to decline from this stage. Finally, the local buckling mode with very low buckling loads appears and the e€ect of contact analysis is very important. 6.3. E€ect of delamination shape on buckling load

In order to investigate the e€ect of the delamination size on the buckling load, a circular delamination existing between the ®rst and second layers is considered. All of the laminated plate speci®cations are the same as the standard plate stated previously except for the delamination diameter, which is used as a variable. Fig. 16 shows the e€ect of the delamination size on the buckling load. In Fig. 16, x-axis denotes the ratio of the delamination area to the plate area, and y-axis denotes the normalized buckling load which is obtained by dividing the present buckling load to that of the intact plate. This normalized buckling load will be used in the following sections. Furthermore, the buckling load with overlap in the buckling mode (buckling load at the beginning of iteration) is shown along with the buckling load of the physically admissible mode (buckling load at the end of iteration) to emphasize on the e€ect of the overlap. To facilitate the following discussion, the buckling mode can be classi®ed into three main types as many other authors [3,9], i.e. global buckling mode, mixed buckling

To study the e€ect of the delamination shape on the buckling load, an elliptical delamination existing between the ®rst and second layers of the laminate is considered and other speci®cations are described in Section 5. Two kinds of models are established by selecting two radii a and b of the ellipse as variables, respectively. In one model, the radius b which is perpendicular to the applied load direction is taken as a variable and the radius a is considered as 4 cm invariably. In another one, b is set to be 4 cm and a varies during computation process. Two kinds of results are shown in Fig. 17. In this ®gure, it can be found that the delamination enlarged in the direction perpendicular to the load direction reduces the buckling load of the laminate signi®cantly. In this case, the buckling mode changes very quickly from a global mode with high buckling load into a local one with very low buckling load. Another important phenomenon is that the e€ect of the contact analysis becomes smaller and smaller gradually with the increase of the delamination size. In this regard, the little e€ort

Fig. 16. E€ect of delamination size on the buckling load.

Fig. 17. E€ect of delamination shape on the buckling load.

1258

N. Hu et al. / Composites Science and Technology 59 (1999) 1247±1260

Fig. 18. E€ect of delamination position on the buckling load for a laminate consisting of 20 plies with 90 ®bers.

or small sti€nesses of arti®cial springs may account for this small increase, since the delamination with a very large area makes the upper thin delaminated layer become very soft, and the overlap can be overcome easily. Fig. 17 also demonstrates that the delamination enlarged on the load direction little in¯uences the buckling load. In this case, the buckling mode changes smoothly from a global mode into a mixed one. 6.4. E€ect of delamination position on buckling load The e€ect of the delamination position through the thickness of the laminate on the buckling load is also investigated. A laminate of 20 plies with only 90 ®bers is analyzed. A circular delamination of radius 4 cm is assumed, and other speci®cations are as those stated in Section 5. The analyzing results are displayed in Fig. 18. Inspection of Fig. 18 reveals that the plate buckles with an apparent local mode in the upper thin delaminated

Fig. 20. Buckling load of a laminate [90/0/90/04/90/0/90] for various delamination sizes and positions.

layer when the delamination is near the surface of the plate. The buckling load is very low in this case. As the delamination goes deeper inside the plate, the buckling load increases rapidly and the mixed and global modes appear in succession. When the thickness of the upper thin delaminated layer is a quarter of the thickness of the laminate, the buckling load reaches its highest level. However, after entering the global buckling mode, the buckling load decreases when the delamination approaches the mid-plane of the plate since the global bending sti€ness of the plate becomes lower in this case. Furthermore, a laminate consisting of 20 plies with stacking sequence (02/902/452/ÿ452)s has also been studied. There is a circular delamination of radius 4 cm in the plate. The results are displayed in Fig. 19. As the delamination moves to the mid-plane, the increasing tendency of the buckling load can be identi®ed from Fig. 19, which is almost the same as that in Fig. 18. Also, there is a considerable jump in the buckling load when the delamination between the fourth and ®fth layers goes to the next position. This can be explained from the fact that the buckling mode changes from the local mode to the mixed one. 6.5. E€ect of ®ber angle in sublaminate on buckling load

Fig. 19. E€ect of delamination position on the buckling load for a laminate consisting of 20 plies with stacking sequence (02/902/452/ ÿ452)s.

The in¯uence of the ®ber angle in the sublaminate (upper thin delaminated layer) on the buckling load is analyzed and discussed here. Two stacking sequences of the laminated plate: [90/0/90/04/90/0/90] and [0/0/90/04/ 90/0/0] are considered. The diameter of a circular delamination existing between di€erent layers is used as a variable in this analysis and other speci®cations are identical to those in Section 5. The in¯uence of the delamination position through the thickness of the laminate is also investigated. The results of the two models are illustrated in Figs. 20 and 21, respectively, where // denotes the position of the delamination.

N. Hu et al. / Composites Science and Technology 59 (1999) 1247±1260

Fig. 21. Buckling load of a laminate [0/0/90/04/90/0/0] for various delamination sizes and positions.

In Fig. 20, there is an interesting phenomenon. By comparing the results of two cases [90//0/90/04/90/0/90] and [90/0//90/04/90/0/90], it can be observed that the buckling load decreases in the second case. The local buckling mode is liable to occur in this case although the thickness of sublaminate increases. This phenomenon is completely di€erent from those in Figs. 18 and 19, which was also found by Marshall et al. [3]. It can be explained from the fact that, after inclusion of 0 ®bers the sublaminte with non-symmetric layup may buckle easily and the e€ect of highly loaded 0 plies in the sublaminate is to lower the buckling stress [10]. The buckling load increases as the delamination moves to the mid-plane of the plate except for the case mentioned above. By inspecting Fig. 21 for case [0/0/90/04/90/0/0], it can be found that the buckling load increases monotonically as the delamination approaches the mid-plane of the plate. As the delamination size increases, Figs. 20 and 21 show that the buckling mode starts with a global buckling mode, changes smoothly into a mixed mode and ends ®nally into a local one. Hence, the variation of the buckling mode can be clearly divided into two stages. When the delamination is close to the upper free

Fig. 22. E€ect of ®ber angle in sublaminate on the buckling load.

1259

surface (thin sublaminate), the local buckling mode is liable to occur and the buckling load decrease dramatically with the increase of the delamination size. This phenomenon is clearly reasonable. To show the e€ect of ®ber angle in the sublaminate on the buckling load more clearly, the results of cases [90//0/90/04/90/0/90] and [0//0/90/04/90/0/0] are drawn in Fig. 22. From Fig. 22, it can be found that the buckling mode of the model with 0 ®bers in the sublaminate changes from the global mode into the local mode very quickly even when the delamination size remains very small. The buckling load drops signi®cantly. This phenomenon can be explained from the fact that 0 ®bers carry an extensive in-plane load and buckle easily [10]. Also, unlike the case of the sublaminate with 90 ®bers, it was found in our analysis that the sublaminate with 0 ®bers generally does not penetrate the remaining part of the laminates. Hence, there is no need to perform the contact analysis in this case. 7. Conclusions The compressive buckling of the laminates with an embedded delamination was analyzed in this study. To avoid the overlap between two delaminated layers in the buckling mode, a solution method for contact analysis in eigenvalue problems was proposed. In this method, the ®rst-order sensitivity analysis and quadratic programming technique were employed to compute ®ctitious forces in the contact area ®rst. Then the sti€ness parameters of some arti®cial springs can be determined through these ®ctitious forces. Finally, the overlap between two delaminated layers in the buckling mode can be eliminated e€ectively by introducing the contribution of these arti®cial springs into the structural sti€ness matrix. For the laminates with an embedded delamination, the e€ect of the introduced contact analysis on the buckling load was analyzed in detail. Furthermore, the e€ects of various sizes, shapes and positions of the delamination and the ®ber angle of the sublaminate on the buckling load have also been investigated. The following conclusions can be drawn: 1. Inclusion of contact analysis to prevent overlaps increases the buckling load. Particularly, when the buckling mode changes signi®cantly after including contact analysis, the buckling load may increase to a great extent and contact analysis is vital in this case. 2. Buckling load decreases as delamination size increases. The variation of the buckling mode can be generally divided into two stages in this process. In the ®rst stage, the buckling mode varies from a global mode to a mixed one with only a small reduction in the buckling load. However, in the

1260

3.

4.

5.

6.

N. Hu et al. / Composites Science and Technology 59 (1999) 1247±1260

second stage, the buckling mode changes from a mixed mode to a local one with considerable reduction in the buckling load. The reduction in the buckling load caused by an elliptical delamination with the short radius in the load direction is more considerable than that caused by a delamination of the same size but with the long radius in the load direction. As the position of the delamination approaches the mid-plane of the laminates, the buckling load increases signi®cantly initially due to the change in the buckling mode from a local mode to a mixed mode or a global one. Thereafter, the buckling load decreases due to the reduction of the global bending sti€ness of the laminates and no obvious change is visible in the buckling mode. A sublaminate with 0 ®bers buckles more easily than that with 90 ®bers. Therefore, 0 plies should be situated far from the outer surface of the laminate, as also reported by Chai and Bacock [10]. In case of a sublaminate with 0 ®bers, one can perform the buckling analysis directly without considering the contact problem, since generally there is no obvious overlap or penetration in the delamination area. However, overlaps are liable to occur when there are 90 ®bers in the sublaminate.

References [1] Chai H, Babcock CD, Knauss WG. One delamination modelling of failure in laminated plates by delamination buckling. Int J Solids Struct 1981;17:1069±83. [2] Simites GJ, Sallam S, Yin WL. E€ect of delamination of axially loaded homogeneous laminated plates. AIAA J. 1985;26:1437±44. [3] Kardomateas GA, Schmueser DW. Buckling and postbuckling of delaminated composites under compressive loads including transverse shear e€ects. AIAA J. 1988;26:337±43. [4] Marshall RD, Sandor€ PE, Lauraitis KN. Buckling of a damaged sublaminate in an impacted laminate. J Comp Tech Res 1988;10:107±13. [5] Shaw D, Tsai MY, Analysis of delamination in compressively loaded laminates. Comp Sci Tech 1989;34:1±17.

[6] Chen HP. Shear deformation theory for compressive delamination buckling and growth. AIAA J 1991;29:813±9. [7] Wang SS, Zahlan NM, Suemasu H. Compressive stability of delmainated random short-®ber composites. II-experimental and analytical results. J Comp Mater 1985;19:317±33. [8] Kutlu Z, Chang FK. Modelling compression failure of laminated composites containing multiple through the width delaminations. J Comp Mater 1992;26:350±87. [9] Lee J, GuÈrdal Z, Grin Jr. OH. Layer-wise approach for the bifurcation problem in laminated composites with delaminations. AIAA J. 1993;31:331±8. [10] Chai H, Bacock C. Two-dimensional modelling of compressive failure in delaminated laminates. J. Comp. Mater. 1985;19:67±98. [11] Peck SO, Springer GS. The behavior of delaminated composite plates. J Comp Mater 1991;25:907±29. [12] Barbero EJ, Reddy JN. Modelling of delamination in composite laminates using a layer-wise plate theory. Int J Solids Struct 1991;28:1069±83. [13] Whitcomb JD. Analysis of a laminate with a postbuckled embedded delamination including contact e€ects. J Comp Mater 1992;26:1523±35. [14] Yeh MK, Tan CM. Buckling of elliptically delaminated composite plates. J Comp Mater 1994;28:36±52. [15] Suemasu H, Gozu K, Hayashi K. Compressive buckling of rectangular composite plates with a free-edge delamination. AIAA J 1995;33:312±9. [16] Hu N, Buckling analysis of delaminated laminates with consideration of contact in buckling mode. Int J Numer Mech Engng (in press). [17] Corr RB, Jenning A. A simultaneous iteration algorithm for symmetric eigenvalue problems. Int J Numer Mech Engng 1976;10:647±63. [18] Bathe KJ, Chaudhary A. A solution method for planar and axisymmetric contact problems. Int J Numer Mech Engng 1985;21:65±88. [19] Oden JT. Exterior penalty methods for contact problems in elasticity. In: Wunderlich W et al., editors. Nonlinear Finite Element Analysis in Structural Mechanics. New York: Springer, 1981. [20] Kong XA, Gakawya A, Cardou A, Cloutier L. A numerical solution of general frictional contact problems by the direct boundary element and mathematical programming approach. Comp Struct 1992;45:95±101. [21] Simunovic S, Saigal S. Quadratic programming contact formulation for elastic bodies using boundary element method. AIAA J 1995;33:325±31. [22] Goldfarb D, Idnani A. A numerically stable dual method for solving strictly convex quadratic programs. Mathematical Progr 1983;27:1±33.