Buckling of cross-ply laminates subject to linearly varying compressive loads and in-plane boundary restraints

Abstract

In-plane prebuckling stresses arise under axial edge loading as a result of a plate boundary being restrained from deformation in the plane of the plate. These stresses exhibit a uniform distribution only under a special set of boundary conditions and in general such uniform in-plane stress distributions do not exist under most boundary conditions encountered in practice. In the present study, in-plane prebuckling stresses are investigated for various in-plane boundary constraints and the stress contour lines are given for rectangular cross-ply laminated plates under linearly varying edge loads. The second part of the article involves the computation of the optimal layer thickness to maximize the buckling load under various combinations of in-plane boundary restraints. The numerical solutions are obtained by finite elements based on first-order shear deformable plate theory which include the in-plane deformations as nodal degrees of freedom. Numerical results show that the exclusion of the in-plane restraints may lead to errors in stability calculations and consequently in the design of laminated plates.

Introduction

In this study, the influence of in-plane boundary restraints on in-plane stresses and the buckling of rectangular cross-ply laminates are investigated under linearly varying compressive edge loads. It is observed that presence of in-plane restraints on the deformation of a plate leads to nonuniform distribution of in-plane stresses which in turn affect the buckling load.

The effect of various in-plane boundary restraints on buckling of plates has been investigated by Harris, ¹ Sherbourne and Pandey, ² Bedair,^3,4 and Nemeth. ⁵ Such restraints are known to lead to thermal buckling as studied by Nemeth^6,7 and Jones. ⁸ These boundary conditions arise in many applications and lead to a nonuniform distribution of in-plane stresses which require a numerical approach for their study. Such stresses, for example, arise in the web sections of transversely loaded composite beams and may lead to premature failure if not taken into account in the design phase. Restraints against in-plane movement typically also arise in aerospace structures where adjacent panels and stiffeners cause these constraints. ⁵

Another complicating factor arises when the plate is subjected to linearly varying in-plane loads leading to buckling. In this case, prebuckling deformations have to be taken into account in the computation of the buckling load.^9–11 This problem has been studied by Papazoglou et al., ¹² Nemeth,^13–15 Zureick and Shih, ¹⁶ Leissa and Kang, ¹⁷ Kang and Leissa, ¹⁸ Wang et al., ¹⁹ and Zhong and Gu.^20,21 Papazoglou et al. ¹² investigated the buckling of symmetric laminates under linearly varying in-plane loads using the Rayleigh-Ritz method based on classical lamination theory. The infinitely long laminates under linearly varying axial loads were studied by Nemeth.^14,15 Leissa and Kang ¹⁷ obtained exact solutions for the free vibration and buckling of rectangular plates subjected to a linearly varying compressive load with two opposite edges simply supported (SS) and the other two edges clamped (C). They extended these results to other boundary conditions as reported elsewhere. ¹⁸ Wang et al. ¹⁹ employed a differential quadrature method to study the vibration and buckling of an SS-C-SS-C plate subject to linearly varying in-plane loads. The effect of various in-plane boundary restraints on the optimization of laminates under linearly varying buckling loads has been studied by Adali and Cagdas. ²² Experimental results on the buckling of anisotropic laminates under linearly varying loads were discussed by Romeo and Ferrero. ²³ Apart from the linearly varying in-plane loads, buckling of laminates under other nonuniform edge loads have been investigated by Kam and Chu, ²⁴ Bert and Devarakonda, ²⁵ Devarakonda and Bert, ²⁶ Jana and Bhaskar, ²⁷ Wang et al.,^28,29 and under in-plane moments by Kang and Leissa. ³⁰ Devarakonda and Bert ²⁶ explained the stress diffusion phenomena and showed that under nonuniform in-plane loads, the membrane force distribution is not constant inside the plate, and presented explicit results for buckling loads in comparison with numerical results obtained using finite element method (FEM).

In many of the previous studies stated above, the buckling loads were calculated assuming the stress distribution in the plate was constant and also the in-plane boundary conditions were not fully described. ² The boundaries were generally defined as either SS or C; however, specifying the in-plane boundary conditions often models a structure more precisely. For example, stiffeners also change the stress distribution considerably in the plates as shown by El-Ghazaly et al. ³¹ In this study, only the effect of in-plane/membrane boundary conditions was investigated and the in-plane stress distributions under external loads were calculated before performing a linearized stability analysis by FEM. It was observed that the exclusion of the in-plane restraints leads to errors in the stability calculations. Reasons for excluding in-plane boundary restraints in many studies include the facts that analytical solutions in the presence of in-plane boundary constraints do not exist or not accurate enough. Similarly, the use of finite elements which do not include the in-plane deformations as nodal degrees of freedom leads to inaccurate computations. To compare the present results with those available in the literature, the boundary conditions where the in-plane stress distribution can be taken as constant are also studied in the present work.

Buckling analysis of laminated composite plates using the first-order shear deformation plate theory (FSDT) was given by Khdeir and Librescu, ³² but this work was restricted to uniform loads. Zhong and Gu ²⁰ developed an analytical solution for the buckling of SS rectangular Reissner–Mindlin plates subjected to linearly varying in-plane loading. Buckling loads obtained by the classical laminated plate theory can have significant errors increasing with an increase in the degree of anisotropy and thickness of the plate.^32,33 Thus, finite element formulations based on classical theory are of limited use in the analysis of composite structures due to the transverse shear effects. ³⁴ In view of this, in the present work the plate-finite element is based on the FSDT. Daripa and Singha ³⁵ studied composite laminates under localized in-plane loads using shear deformable plate-finite elements.

The present work investigates the in-plane stress distribution of cross-ply laminates under linearly varying edge loads and in-plane boundary restraints using finite elements based on the FSDT. An outline of the finite element formulation and the description of the finite elements are discussed in the next section. Next, the finite element solution is applied to a number of verification problems to assess its accuracy. Subsequently the prebuckling stresses are studied under linearly varying edge loads. A related problem is the computation of the optimal layer thicknesses to maximize the buckling load. The results for the optimal design problem are given for linearly varying and uniform edge loads and the efficiencies of the designs are obtained.

Finite element formulation

Next, the formulation for the parabolic isoparametric plate element used in the calculations is discussed. One of the reasons for selecting the 8-noded elements is to avoid the parasitic shear or shear locking problem which arises when linear isoparametric elements are used. Using the curved elements and making use of reduced integration, the shear locking problem can be avoided. The second reason is the ease in modeling complex geometries. Note that, the finite element code used in this study is written in Matlab^® environment.

First-order shear deformable plate theory

The displacement field corresponding to FSDT for composite laminated plates is as follows:

$\begin{array}{l}u(x,y,z)=u_0(x,y)+z{\text{θ}}_x(x,y)\\ v(x,y,z)=v_0(x,y)+z{\text{θ}}_y(x,y)\\ w(x,y,z)=w_0(x,y)\end{array}$

1

where u ₀, v ₀, and w ₀ are the mid-surface displacement components and θ _x , and θ _y are the rotations of the transverse normal about the y and x axes. The x–y plane coincides with the mid-surface of the plate, with the plate thickness denoted by t. The strain–displacement relations are given by

$\mathbf{\chi }=\left\{\begin{array}{c}\begin{array}{l}{\text{χ}}_x\\ {\text{χ}}_y\\ {\text{χ}}_{xy}\end{array}\end{array}\right\}=\left\{\begin{array}{c}\begin{array}{l}{\text{θ}}_{x,x}\\ {\text{θ}}_{y,y}\\ ({\text{θ}}_{x,y}+{\text{θ}}_{y,x})\end{array}\end{array}\right\}{\epsilon }^0=\left\{\begin{array}{c}\begin{array}{l}{\text{ε}}_x^0\\ {\text{ε}}_y^0\\ {\text{γ}}_{xy}^0\end{array}\end{array}\right\}=\left\{\begin{array}{c}\begin{array}{l}{u}_{0,x}\\ v_{0,y}\\ (u_{0,y}+v_{0,x})\end{array}\end{array}\right\}\mathbf{\phi }=\left\{\begin{array}{l}{\text{φ}}_x\\ {\text{φ}}_y\end{array}\right\}=\left\{\begin{array}{l}{\text{θ}}_x+w_{0,x}\\ {\text{θ}}_y+w_{0,y}\end{array}\right\}$

2

where φ _x and φ _y are the rotations due to shear about the y and x axes; and χ _x , χ _y , and χ _xy and ${\text{ε}}_x^0$ , ${\text{ε}}_y^0$ , and ${\text{ε}}_{xy}^0$  are the curvatures and the mid-surface strains, respectively. A comma denotes partial differentiation with respect to the following variable. These equations can be expressed in matrix form as

$\mathbf{\text{χ}}={\mathbf{\text{B}}}_{\mathbf{\text{χ}}}\mathbf{\text{δ}},{\mathbf{\text{ε}}}^0={\mathbf{\text{B}}}_{\mathbf{\text{ε}}}\mathbf{\text{δ}},\mathbf{\text{φ}}={\mathbf{\text{B}}}_{\mathbf{\text{φ}}}\mathbf{\text{δ \hspace{1em} where \hspace{0.17em} δ}}={\left\{\begin{array}{ccccc}{u}_0& v_0& w_0& {\text{θ}}_x& {\text{θ}}_y\end{array}\right\}}^{\mathrm{T.}}$

Moments, normal forces, and shear forces are given in equations (3)–(5), as follows:

$\mathbf{\text{M}}=\left\{\begin{array}{c}{M}_x\\ M_y\\ M_{xy}\end{array}\right\}=\mathbf{\text{D χ}}+\mathbf{\text{B}}{\mathbf{\text{ε}}}^0=\mathbf{\text{D}}\left\{\begin{array}{c}{\text{χ}}_x\\ {\text{χ}}_y\\ {\text{χ}}_{xy}\end{array}\right\}+\mathbf{\text{B}}\left\{\begin{array}{c}\begin{array}{l}{\text{ε}}_x^0\\ {\text{ε}}_y^0\\ {\text{γ}}_{xy}^0\end{array}\end{array}\right\}$

3

$\mathbf{\text{N}}=\left\{\begin{array}{c}{N}_x\\ N_y\\ N_{xy}\end{array}\right\}=\mathbf{\text{B χ}}+\mathbf{\text{A}}{\mathbf{\text{ε}}}^0=\mathbf{\text{B}}\left\{\begin{array}{c}{\text{χ}}_x\\ {\text{χ}}_y\\ {\text{χ}}_{xy}\end{array}\right\}+\mathbf{\text{A}}\left\{\begin{array}{c}\begin{array}{l}{\text{ε}}_x^0\\ {\text{ε}}_y^0\\ {\text{γ}}_{xy}^0\end{array}\end{array}\right\}$

4

$\mathbf{\text{Q}}=\left\{\begin{array}{l}{Q}_x\\ Q_y\end{array}\right\}={\mathbf{\text{C}}}_{\mathbf{\text{s}}}\left\{\begin{array}{l}{\text{φ}}_x\\ {\text{φ}}_y\end{array}\right\}$

5

where for a laminated plate the matrices A, B, D, and C_s are given by Reddy. ³⁶ Note that the shear correction factors are considered as equal to 5/6, in this study.

Linear static analysis using FEM

The interpolation function of the displacement field is defined as follows:

${\left\{\begin{array}{ccccc}{u}_0& v_0& w_0& {\text{θ}}_x& {\text{θ}}_y\end{array}\right\}}^T={\displaystyle \sum _{i=1}^8{S}_i{\mathbf{\text{δ}}}_i},$

6

where S_i is the shape function at node i, and the displacement vector at node i is ${\mathbf{\text{δ}}}_i={\left\{\begin{array}{ccccc}{u}_{0i}& v_{0i}& w_{0i}& {\text{θ}}_{xi}& {\text{θ}}_{yi}\end{array}\right\}}^T$ . The shape functions are given by Cook. ³⁷ The strain displacement relations for an element are given as follows:

${\mathbf{\text{ε}}}^0=[{\mathbf{\text{B}}}_{\text{ε}}^1,\dots ,{\mathbf{\text{B}}}_{\mathbf{\text{ε}}}^8]{\mathbf{\text{δ}}}_e$

7

$\mathbf{\text{χ}}=[{\mathbf{\text{B}}}_{\mathbf{\text{χ}}}^1,\dots ,{\mathbf{\text{B}}}_{\mathbf{\text{χ}}}^8]{\mathbf{\text{δ}}}_e$

8

$\mathbf{\text{φ}}=[{\mathbf{\text{B}}}_{\mathbf{\text{φ}}}^1,\dots ,{\mathbf{\text{B}}}_{\mathbf{\text{φ}}}^8]{\mathbf{\text{δ}}}_e$

9

where δ _e is the element displacement vector. Using these equations, the strain energy of element e, denoted by U_e , can be written as

$U_e=\frac{1}{2}{\displaystyle \underset{A}{\int }({\mathbf{\text{M}}}^T\mathbf{\text{χ}}+{\mathbf{\text{Q}}}^T\text{φ}+{\mathbf{\text{N}}}^T{\text{ε}}^0)}dA$

10

The total potential energy can be minimized by differentiating U_e with respect to δ _e which results in the element stiffness matrix given by

${\mathbf{\text{K}}}_e={\displaystyle \underset{-1}{\overset{1}{\int }}{\displaystyle \underset{-1}{\overset{1}{\int }}\left({\mathbf{\text{B}}}_{\mathbf{\text{χ}}}^T\mathbf{\text{D}}{\mathbf{\text{B}}}_{\mathbf{\text{χ}}}+{\mathbf{\text{B}}}_{\mathbf{\text{ε}}}^T\mathbf{\text{B}}{\mathbf{\text{B}}}_{\mathbf{\text{χ}}}+{\mathbf{\text{B}}}_{\mathbf{\text{φ}}}^T{\mathbf{\text{C}}}_s{\mathbf{\text{B}}}_{\mathbf{\text{φ}}}+{\mathbf{\text{B}}}_{\mathbf{\text{ε}}}^T\mathbf{\text{A}}{\mathbf{\text{B}}}_{\mathbf{\text{ε}}}+{\mathbf{\text{B}}}_{\mathbf{\text{χ}}}^T\mathbf{\text{B}}{\mathbf{\text{B}}}_{\mathbf{\text{ε}}}\right)\left|J\right|d\text{ξ}d\text{η}}}$

11

where |J| denotes the determinant of the Jacobian matrix and (ξ,η) denotes the coordinates of any point in the curvilinear coordinate system.

Linearized stability analysis using FEM

The steps of the solution can be briefly explained as follows ³⁷ :

• Calculate the membrane forces using a standard plane-stress analysis.

• Calculate the initial stress or geometric stiffness matrix using the membrane stresses calculated in step (i).

• Calculate the buckling load parameter λ, which is defined as the ratio of actual buckling load to the applied forces at which the plate buckles, by solving the eigenvalue problem which is solved using the Matlab^® function as follows:

$\text{det( K}-\text{λ}{\mathbf{\text{K}}}_{\mathbf{\text{G}}}\text{)}=0$

12

where K is the global stiffness matrix and K_G is the global geometric stiffness matrix. Note that the solution is only valid for symmetric laminates where there is no bending–extension coupling. Next, the derivation of the element geometric stiffness matrix, ${\mathbf{\text{K}}}_{\mathbf{\text{G}}}^e$ , is discussed.

Element geometric stiffness matrix

The element geometric stiffness matrix is derived considering that the plate element is subjected to a prescribed in-plane stress system leading to buckling. The nonlinear terms of the mid-surface strains associated with the von Karman plate theory are given by Reddy ³⁶

${\mathbf{\text{ε}}}_{NL}^0={\left\{\begin{array}{ccc}\frac{1}{2}{w}_{0,x}^2& \frac{1}{2}{w}_{0,y}^2& w_{0,x}{w}_{0,y}\end{array}\right\}}^T$

13

The potential energy of the applied stresses denoted by V_e arising from the action of the in-plane stresses on the corresponding second-order strains is given in

$V_e={\displaystyle \underset{A}{\int }\left(\left\{\begin{array}{ccc}{N}_{xx}& N_{yy}& N_{xy}\end{array}\right\}{\mathbf{\text{ε}}}_{NL}^0\right)}dA$

14

where N_xx , N_xy , and N_yy are the in-plane forces at the Gauss points which are calculated by a primary static analysis. For this purpose, a linear static analysis should be performed prior to linearized stability analysis. Thus, the calculation of the buckling load parameter involves two stages, the first stage being the linear static analysis and the second stage being the linearized stability analysis. Rewriting the terms inside the integral yields:

$V_e=\frac{1}{2}{\displaystyle \underset{A}{\int }{\mathbf{\text{δ}}}_e^T({\mathbf{\text{G}}}^T\mathbf{\text{S G}}){\mathbf{\text{δ}}}_e}dA\text{\hspace{1em} where, \hspace{0.17em} G}{\mathbf{\text{δ}}}^e=[{\mathbf{\text{G}}}_1,\dots ,{\mathbf{\text{G}}}_8]{\{{\mathbf{\text{δ}}}_1\dots {\mathbf{\text{δ}}}_8\}}^T={\left\{w_{0,x}{w}_{0,y}\right\}}^T$

15

${\mathbf{\text{G}}}_i=\left[\begin{array}{ccccc}0& 0& N_{i,x}& 0& 0\\ 0& 0& N_{i,y}& 0& 0\end{array}\right]\text{\hspace{1em} and \hspace{1em} S}=\left[\begin{array}{cc}{N}_{xx}& N_{xy}\\ N_{xy}& N_{yy}\end{array}\right]$

16

The element geometric stiffness matrix can be obtained by differentiating V_e with respect to δ _e and can be expressed as

${\mathbf{\text{K}}}_{\mathbf{\text{G}}}^e={\displaystyle \underset{-1}{\overset{1}{\int }}{\displaystyle \underset{-1}{\overset{1}{\int }}{\mathbf{\text{G}}}^T\mathbf{\text{S G}}\left|\mathbf{\text{J}}\right|d\text{ξ}d\text{η.}}}$

17

Stress evaluation-smoothing technique

The Gauss points are the best points to calculate stresses; however, especially for comparison purposes, the stresses at the element nodes are also required. To solve this problem, a technique described by Barlow ³⁸ and in Hinton and Owen, ³⁹ which is a bilinear extrapolation of the 2 × 2 Gauss point stress values, is used.

Verification problems

Two verification problems are solved in order to validate the finite element code discussed in the section on finite element formulation.

Buckling of a rectangular plate subject to linearly varying loads

Stability analysis results for a SS homogeneous isotropic plate of dimensions a × b under the action of a linearly varying axial load are presented in comparison with the results presented by Timoshenko. ⁴⁰ The geometry and loading details are shown in Figure 1 and the edge load, N_xx , is defined in Eq. (18). The critical value of stress, σ _cr , is defined as ${\text{σ}}_{cr}=k({\text{π}}^{2}D)\text{/}(b^{2}H)$ , where $D=E{H}^3\text{/}12(1-{\nu }^2)$  represents the bending rigidity of the plate and H the plate thickness.

Figure 1.

Plate geometry and linearly varying edge load.

As shown in Table 1, the numerical results obtained using a 16 × 8 finite element mesh are in good agreement with the reference results. However, the finite element results obtained are slightly lower than the results presented by Timoshenko ⁴⁰ due to the influence of shear deformation taken into account here.

Table 1.

Numerical values of k in comparison with the analytical results with H/b = 0.01.

α	a/b=	0.4	0.5	0.6	2/3	0.75	0.8	0.9	1.0	1.5
2	a	29.100	25.600	24.100	23.900	24.100	24.400	25.600	25.600	24.100
	b	29.024	25.487	24.096	23.683	24.099	24.459	25.571	25.512	24.112
4/3	a	18.700	–	12.900	–	11.500	11.200	–	11.000	11.500
	b	18.824	15.012	12.956	12.123	11.467	11.224	10.988	11.004	11.472
1	a	15.100	–	9.700	–	8.400	8.100	–	7.800	8.400
	b	15.110	11.603	9.729	8.698	8.359	8.123	7.862	7.806	8.361
4/5	a	13.300	–	8.300	–	7.100	6.900	–	6.600	7.100
	b	13.274	10.043	8.340	7.654	7.105	6.891	6.650	6.590	7.107
2/3	a	10.800	–	7.100	–	6.100	6.000	–	5.800	6.100
	b	12.205	9.171	7.590	6.950	6.441	6.243	6.018	5.959	6.443
0	a	8.410	6.250	5.140	–	–	4.200	4.040	4.000	–
	b	8.389	6.238	5.130	4.688	4.335	4.198	4.041	3.997	4.337

^a Timoshenko. ⁴⁰

^b Current study using a 16 × 8 mesh.

Buckling of an antisymmetric cross-ply laminate under uniaxial compression

Stability analysis results for a SS rectangular antisymmetric cross-ply laminated plate subject to a uniform axial compression are presented in comparison with the results presented by Reddy and Miravete. ⁴¹ The lamina thicknesses are equal and the material properties are taken as

E ₁ /E ₂ = 25, ν ₁₂ = 0.25, G ₁₂ /E ₂ = G ₁₃ /E ₂ = 0.6, G ₂₃/E ₂ = 0.25.

The results obtained are tabulated in Table 2, in comparison with the reference results, where the abbreviations 2Q8R, 2Q9R, and CFS stand for 8 node finite elements with reduced integration, 9 node finite elements with reduced integration, and closed form solution, respectively. The plate is modeled using a 16 × 8 regular mesh. The nondimensional buckling load, $\overline{N}$ is defined as $\overline{N}=(N_{xx}{b}^2)\text{/}(E_2{H}^3)$ . Note that the prebuckling stress distribution is assumed to be uniform and the prebuckling deformations are neglected in order to compare the results obtained here with the reference results. Excellent agreement with the reference results can be observed from Table 2.

Table 2.

Nondimensional buckling loads for antisymmetric cross-ply laminated plates under uniaxial compression, a/H = 10

Layers	Source	E 1/E 2=	3	10	20	30	40
2	CFSa		4.6948	6.1181	7.8196	9.3746	10.8170
	2Q9Ra		4.7718	6.2465	8.0423	9.7347	11.3553
	2Q8Rb		4.7437	6.1904	7.9667	9.6444	11.2501
4	CFSa		5.7138	9.0164	13.7430	17.7830	21.2800
	2Q9Ra		5.2543	9.2552	14.3320	18.8150	22.8060
	2Q8Rb		5.2401	9.2338	14.3025	18.7802	22.7677
6	CFSa		5.2670	9.6051	15.0010	19.6394	23.6690
	2Q9Ra		5.3430	9.7893	15.3940	20.2800	24.5770
	2Q8Rb		5.3317	9.7773	15.3797	20.2636	24.5598
10	CFSa		5.3159	9.9134	15.6690	20.6347	24.9640
	2Q9Ra		5.3884	10.0600	15.9270	21.0080	25.4500
	2Q8Rb		5.3786	10.0533	15.9206	21.0010	25.4429

^a Reddy and Miravete. ⁴¹

^b Current study using a 16 × 8 mesh.

Problem definition

A rectangular plate of lateral dimensions a × b and thickness H is subjected to a linearly varying compressive edge loading as shown in Figure 1. N_xx , N_yy , and N_xy are normal forces per unit length in the x and y directions and shearing force per unit length in the xy plane, respectively. The in-plane edge forces N_xx , N_yy , and N_xy along the edges x = 0 and x = a are given in Eq. (18) with the other edges being stress free, as follows:

$N_{xx}=-N_0\left(1-\alpha \frac{y}{b}\right),N_{yy}=0,N_{xy}=0$

18

where N ₀ is the intensity of N_xx at the edge y = 0 and α is a constant. By changing α, we can obtain different load cases, that is, α = 2, α < 2, and α > 2 corresponds to pure bending, a combination of bending and compression, and a combination of bending and tension, respectively. For all values of α, N_xx  = N ₀ at y = 0. In particular, N_xx  = 0 at y = b for α = 1, N_xx  = 0 at y = b/2 and N_xx  = −N ₀ at y = b for α = 2 which are the α values used in the numerical study of the problem with the case α = 0 corresponding to the uniform compression. The nondimensional buckling load parameter, λ, is defined as

$\text{λ}=\frac{N_0{b}^2}{{\text{π}}^2{D}_{11}}$

19

where D ₁₁ denotes the flexural rigidity of the plate defined by $D_{11}=\frac{E_1{H}^3}{12(1-{\nu }_{12}^2)}$ .

The entries of S matrix given in Eq. (16) are the in-plane forces at the Gauss points which are calculated by a primary static analysis.

Boundary conditions

The four different sets of boundary conditions are tabulated in Table 3 and schematically shown in Figure 2, where the sides AC and BD refer to x = 0 and x = a, respectively, and the sides AB and CD refer to y = 0 and y = b, respectively, as shown in Figure 1. Rotational restraints are the same for SS-I and SS-II boundaries and for C-I and C-II boundaries, the only differences being the in-plane boundary conditions. The boundary sets SS-I and SS-II are grouped as the SS plates and the boundary sets C-I and C-II are grouped as the C plates.

Figure 2.

Schematic representation of the boundary conditions simply supported (SS)-I, SS-II, clamped (C)-I, and C-II edges.

Table 3.

Boundary conditions; 1: free, 0: restrained

Edge	Boundary set	u	v	w	θ x	θ y
AC	SS-I	0	0	0	1	0
	SS-II	0	0	0	1	0
	C-I	0	0	0	0	0
	C-II	0	0	0	0	0
BD	SS-I	1	1	0	1	0
	SS-II	1	1	0	1	0
	C-I	1	1	0	0	0
	C-II	1	1	0	0	0
AB	SS-I	1	1	0	0	1
	SS-II	0	0	0	0	1
	C-I	1	1	0	0	0
	C-II	0	0	0	0	0
CD	SS-I	1	1	0	0	1
	SS-II	0	0	0	0	1
	C-I	1	1	0	0	0
	C-II	0	0	0	0	0

SS: simply supported edge, C: clamped edge.

Design problem

The design objective is to maximize the nondimensional buckling load of the symmetrically laminated rectangular plate of total thickness H, total layer number n by optimizing the ply thicknesses, H_i , where i = 1, 2, …,n. The design problem is defined below as:

Determine the ply thicknesses H_i , i = 1, 2, …,n of a symmetric cross-ply laminate of total thickness H such that the buckling load parameter,λ, will be maximized, viz.

$\underset{H_i}{\max }\text{\hspace{0.17em} λ \hspace{1em} such \hspace{0.17em} that \hspace{1em}}{\displaystyle \sum _{i=1}^n{H}_i}=H\text{\hspace{1em} where \hspace{0.17em}}i=1,2,\mathrm{...},n$

20

Numerical results

The numerical results are obtained using a 16 × 8 mesh, which consists of 8 noded parabolic isoparametric plate-finite elements. The numerical results obtained in the verification problems and a short convergence study show that this mesh is capable of accurately modeling the plates considered. The side-to-thickness ratio b/H was taken as equal to 25 and 100 in the numerical examples. A relatively high b/H ratio will minimize the contribution to shear deformation and a lower b/H ratio will increase the contribution to shear.

The material properties are taken as follows:

E ₁/E ₂ = 25, G ₁₂/E ₂ = G ₁₃/E ₂ = 0.5, G ₂₃/E ₂ = 0.3, v ₁₂ = v ₂₃ = v ₁₃ = 0.25

Analysis of in-plane stresses

We first analyze the effect of edge loading and the boundary conditions on the in-plane stress distribution. The applied edge loads are converted into nodal loads by numerical integration and a primary linear static analysis is performed in order to obtain the normal force distribution inside the plate. The normal forces at the Gauss points are extrapolated to the element corner nodes to plot the normal force distributions using the previously described stress evaluation-smoothing technique.

We consider a cross-ply laminate with four layers of equal thickness under axial compression N_xx given in Eq. (18), having a stacking sequence (0°/90°) _s and an aspect ratio of a/b = 1.

Figure 3 shows the in-plane force distribution for SS-I and C-I boundary conditions when there are no constraints at the three edges AB, BD, and CD in x and y directions, that is, when u(x,y) and v(x,y) are free. Figure 4 shows the in-plane stress distribution of the same laminate for the boundary conditions SS-II and C-II for α = 0, 1, 2. In this case, u(x,0) = u(x,b) = 0 and v(x,0) = v(x,b) = 0, that is, the in-plane deflections of the horizontal sides of the plate (Figure 1) are restrained. Moreover, u(0,y) = v(0,y) = 0, but u(a,y) and v(a,y) are free, that is the in-plane deflections are restrained only on the left side of the plate. The resulting distribution of N_xx is shown in Figure 4 which shows Poisson’s effect due to restrained horizontal edges as well as the nonsymmetric stress distribution due to the nonsymmetry of the boundary conditions. Figure 4a corresponding to a uniformly distributed axial load (α = 0) shows that even under this loading, the imposed boundary conditions lead to a nonuniform stress distribution. Figure 4 also shows that the distribution of in-plane forces depends not only on the boundary conditions but also on the in-plane loading.

Figure 3.

Distribution of normal force N_xx for the cross-ply laminate (0°/90°), for N ₀ = 1, simply supported (SS)-I and clamped (C)-I edges. (a) α = 0, (b) α = 1, and (c) α = 2.

Figure 4.

Distribution of normal force N_xx for the cross-ply laminate (0°/90°) _s , for N ₀ = 1, simply supported (SS)-II and clamped (C)-II edges. (a) α = 0, (b) α = 1, and (c) α = 2.

Design for maximum buckling load

Next, we investigate the effect of various problem parameters on the optimal layer thickness and the maximum buckling load. A four-layered symmetric cross-ply laminate with a stacking sequence (0°/90°/90°/0°) and aspect ratio a/b = 1 is considered. The nondimensional thicknesses of the 0° and 90° plies are given by h ₀ = H ₀/H and h ₉₀ = H ₉₀/H, respectively, where H ₀ and H ₉₀ are the thicknesses of 0° and 90° plies. h ₉₀ can be written in terms of h ₀ as h ₉₀ = 0.5 − h ₀ so that the unknown is only h ₀. Figure 5 shows the curves of λ plotted against h ₀ for boundary conditions SS-I and SS-II (Figure 5a) and C-I and C-II (Figure 5b) with α = 0, 1, 2 and b/H = 100. In this figure, λ values are normalized with respect to its value for the equal layer thickness laminate (h ₀ = 0.25). It is observed that pure bending, that is, α = 2, yields the highest and the lowest values of λ for all boundary conditions, emphasizing the importance of an optimal design. This case is also more sensitive to the changes in layer thicknesses compared to α = 0 and α = 1. The least sensitiveness to h ₀ are due to the cases for α = 0, 1 for the boundary conditions SS-I and C-I. Thus, an optimal design based on these boundary conditions will show little improvement over the laminate with equal thickness layers (h ₀ = h ₉₀ = 0.25). It is also observed that the optimal values of h ₀, in general, depend on the presence or absence of in-plane boundary restraints. Furthermore, the buckling load of cross-ply laminates can be increased considerably by optimal thickness design if the boundaries are restrained against in-plane movement, that is, for SS-II and C-II. For SS-II and α = 2, the increase in the buckling load is more than 50% as compared to a laminate with equal layer thicknesses (Figure 5a), and the corresponding value for C-II and α = 2 is 30%. However, the sensitivity of buckling load to the optimal thickness needs to be taken into account in a design situation.

Figure 5.

λ versus h ₀ for a cross-ply laminate (0°/90°) _s , a/b = 0, and b/H = 100.

The corresponding results for b/H = 25 are shown in Figure 6, which shows that the effect of shear deformation on the buckling load is more pronounced in the case of C plates.

Figure 6.

λ versus h ₀ for a cross-ply laminate (0°/90°) _s , a/b = 1, and b/H = 25.

Optimal layer thicknesses

Numerical results for the design problem stated in Eq. (20) are considered for a (0°/90°/90°/0°) laminate and can be stated as

$\underset{h_0}{\text{max}}\text{\hspace{0.17em} λ \hspace{1em} such \hspace{0.17em} that \hspace{0.17em}}2h_0+2h_{90}=1$

21

which is a one-dimensional optimization problem which is solved by the Golden section method. The optimization results for the thickness/height ratio b/H = 100 are shown in Tables 4–6 for α = 0, 1, 2, respectively, for different values of the aspect ratios a/b. In these tables, λ₀ corresponds to the case h ₀ = 0.25, that is equal layer thickness laminate. λ_max and h _opt denote the maximum value of the nondimensional buckling load λ and the optimal value of h ₀. Efficiency (η) is defined as $\text{η}=\frac{{\text{λ}}_{\text{max}}-{\text{λ}}_0}{{\text{λ}}_0}\times 100$ , that is the percentage increase in the buckling load compared to an equal layer thickness laminate.

Table 4.

Optimization results for α = 0 (uniform edge load), b/H = 100

a/b	Boundary set	λ0	h opt	λmax	η (%)
0.5	SS-I	3.37	0.50	3.79	12.6
	SS-II	3.39	0.48	3.79	11.8
	C-I	12.42	0.50	14.08	13.3
	C-II	12.44	0.50	14.05	12.9
1	SS-I	1.07	0.08	1.07	0.1
	SS-II	1.13	0.05	1.31	15.9
	C-I	3.92	0.50	3.94	0.5
	C-II	4.11	0.08	4.39	6.7
2	SS-I	0.90	0.07	1.07	18.5
	SS-II	1.10	0.09	1.54	40.1
	C-I	2.41	0.12	2.83	38.1
	C-II	3.18	0.09	4.39	38.3

SS: simply supported edge, C: clamped edge.

Table 5.

Optimization results for α = 1 (triangular edge load), b/H = 100

a/b	Boundary set	λ0	h opt	λmax	η (%)
0.5	SS-I	5.32	0.39	5.49	3.2
	SS-II	5.49	0.35	5.63	2.6
	C-I	18.76	0.50	19.34	3.1
	C-II	19.32	0.35	19.83	2.6
1	SS-I	2.03	0.04	2.11	4.0
	SS-II	2.19	0.05	2.59	18.3
	C-I	7.25	0.16	7.40	2.0
	C-II	7.84	0.08	8.62	10.0
2	SS-I	1.78	0.06	2.11	18.3
	SS-II	2.18	0.09	3.04	39.5
	C-I	4.58	0.13	5.44	18.8
	C-II	6.23	0.09	8.64	38.5

SS: simply supported edge, C: clamped edge.

Table 6.

Optimization results for α = 2 (pure bending edge load), b/H = 100

a/b	Boundary set	λ0	h opt	λmax	η (%)
0.5	SS-I	8.15	0.20	8.23	1.0
	SS-II	8.99	0.16	9.30	3.5
	C-I	27.37	0.20	27.48	0.4
	C-II	30.22	0.19	30.74	1.7
1	SS-I	5.20	0.12	7.83	50.7
	SS-II	7.19	0.10	11.03	53.4
	C-I	15.65	0.18	17.72	13.2
	C-II	21.49	0.09	27.78	29.2
2	SS-I	5.20	0.11	6.98	34.1
	SS-II	8.55	0.10	11.26	31.8
	C-I	11.57	0.12	14.24	23.1
	C-II	22.03	0.10	29.56	34.1

SS: simply supported edge, C: clamped edge.

It is observed that the aspect ratio a/b = 2 and α = 2, in general, gives higher design efficiency; however, the highest efficiencies (more than 50%) are obtained for the aspect ratio a/b = 1 and α = 2 (Table 6). The optimal values of the layer thickness h _opt depend on the specific boundary condition and in particular on the absence or presence of boundary restraints as can be seen by comparing the results for SS-I with SS-II and C-I with C-II. Thus, in applications where such restraints exist, neglecting them will lead to suboptimal designs. Comparing the buckling loads of laminates with SS-I and SS-II indicates that restrained boundaries (SS-II) lead to higher buckling loads. Moreover, the difference increases with increasing aspect ratio. Similar results hold for C-I and C-II. It is observed that in the case of pure bending (α = 2; Table 6), h _opt remains fairly uniform for a given aspect ratio.

Next the effect of the aspect ratio on the efficiency, optimal layer thickness, and the maximum buckling load is studied. Corresponding results for b/H = 25 are tabulated in Tables 7–9. Results for efficiency, h _opt, and λ_max versus the aspect ratio are given for α = 0 in Figure 7, for α = 1 in Figure 8, and for α = 2 in Figure 9 for SS-I, SS-II, C-I, and C-II (restrained boundaries) with b/H = 100. It is observed that efficiencies and the optimal h ₀ are not monotonic functions of the aspect ratio and h _opt fluctuates in a bigger range for α = 0 (Figure 7b) and a smaller range for α = 2 (Figure 9b). Optimization results for b/H = 25 are shown in Figures 10 –12. It is observed that the shear deformation has the largest effect on the optimal layer thicknesses of laminates under pure bending edge loads, that is α = 2 as a comparison of Figures 9 and 12 indicates.

Figure 7.

Efficiency, h _opt, and buckling load versus a/b, α = 0 (uniform edge load), and b/H = 100.

Figure 8.

Efficiency, h _opt, and buckling load versus a/b, α = 1 (triangular edge load), and b/H = 100.

Figure 9.

Efficiency, h _opt, and buckling load versus a/b, α = 2 (pure bending edge load), and b/H = 100.

Figure 10.

Efficiency, h _opt, and buckling load versus a/b, α = 0 (uniform edge load), and b/H = 25.

Figure 11.

Efficiency, h _opt, and buckling load versus a/b, α = 1 (triangular edge load), and b/H = 25.

Figure 12.

Efficiency, h _opt, and buckling load versus a/b, α = 2 (pure bending edge load), and b/H = 25.

Table 7.

Optimization results for α = 0 (uniform edge load), b/H = 25

a/b	Boundary set	λ0	h opt	λmax	η (%)
0.5	SS-I	0.64	0.49	0.73	14.6
	SS-II	0.64	0.49	0.73	13.9
	C-I	1.43	0.49	1.68	17.4
	C-II	1.43	0.49	1.67	17.0
1	SS-I	0.25	0.08	0.25	2.1
	SS-II	0.26	0.05	0.31	16.3
	C-I	0.76	0.19	0.77	0.7
	C-II	0.80	0.11	0.85	6.1
2	SS-I	0.22	0.08	0.25	14.6
	SS-II	0.27	0.10	0.36	34.1
	C-I	0.52	0.13	0.58	12.2
	C-II	0.69	0.10	0.87	27.2

SS: simply supported edge, C: clamped edge.

Table 8.

Optimization results for α = 1 (triangular edge load), b/H = 25

a/b	Boundary set	λ0	h opt	λmax	η (%)
0.5	SS-I	1.02	0.49	1.07	4.6
	SS-II	1.05	0.49	1.08	2.8
	C-I	2.19	0.50	2.36	7.8
	C-II	2.26	0.50	2.39	5.4
1	SS-I	0.47	0.06	0.50	5.5
	SS-II	0.51	0.05	0.60	18.5
	C-I	1.39	0.18	1.42	1.9
	C-II	1.51	0.11	1.63	7.8
2	SS-I	0.44	0.06	0.50	14.2
	SS-II	0.53	0.10	0.70	33.2
	C-I	0.98	0.14	1.10	12.4
	C-II	1.34	0.11	1.68	26.0

SS: simply supported edge, C: clamped edge.

Table 9.

Optimization results for α = 2 (pure bending edge load), b/H = 25

a/b	Boundary set	λ0	h opt	λmax	η (%)
0.5	SS-I	1.58	0.18	1.62	2.2
	SS-II	1.74	0.13	1.85	6.3
	C-I	2.86	0.31	2.98	4.1
	C-II	3.00	0.36	3.19	6.2
1	SS-I	1.21	0.15	1.61	33.5
	SS-II	1.63	0.12	2.10	28.4
	C-I	2.56	0.20	2.61	2.2
	C-II	3.20	0.27	3.21	0.3
2	SS-I	1.21	0.11	1.51	25.3
	SS-II	1.82	0.13	2.13	16.8
	C-I	2.21	0.16	2.41	8.9
	C-II	3.39	0.24	3.39	0.1

SS: simply supported edge, C: clamped edge.

Conclusions

The effect of in-plane boundary restraints was investigated using the in-plane prebuckling stresses of cross-ply laminates under linearly varying edge loads. Results were given for several combinations of boundary restraints in the form of counter lines. Three edge loads were studied, namely uniform compression for comparison purposes, triangular compressive load, and pure bending edge load. The constant stress distributions under these loads were also shown, which occur under specific boundary conditions that do not apply when the in-plane movement of the plate is constrained. After analyzing the buckling loads with respect to layer thicknesses (Figures 5 and 6), results for the optimal layer thicknesses are given in table and graphical form for various boundary conditions and edge loads. The results presented in the tables were obtained for rectangular plates with aspect ratios a/b = 0.5, 1, 2 under six different combinations of boundary conditions and include the efficiencies of the optimal designs. As there are no analytical solutions for the problems formulated in the present study due to the complicating factors, a finite element method was developed and implemented to calculate the buckling loads and prebuckling stresses. The solution is based on the FSDT and the influence of thickness, that is the contribution of shear deformation was examined by taking the b/H ratio as equal to 25 and 100.

Numerical results show that, as expected, the buckling loads are higher for C plates compared with SS plates, and the presence of in-plane restraints causes an increase in buckling load and, in general, the efficiency increases as a/b and nonuniformity of the edge loads, that is, increase in α.