## Abstract

Approximate expressions for buckling loads of long rectangular plates with flexural/twist anisotropy and simply supported edges are derived. The novelty of the current approach lies in the use of bounded non-dimensional quantities to reduce the complexity of formulation to such a state that closed form solutions are possible. Such expressions are validated with respect to detailed finite element analysis and numerical solutions found in the literature. A number of studies are conducted using a composite material with extremely large ply stiffness in the fibre direction and whose layers are arranged into extreme stacking sequences. Such an arrangement maximizes anisotropic effects and so challenges the validity of the solutions in a most testing manner. Despite such extreme studies excellent agreement is found.

## 1. Introduction

Composite materials continue to find increasing use in primary structural parts of large aircraft. Two current examples include wings of the large military transporter, Airbus A400M and both fuselage and tailplane of the superjumbo, Airbus A380. As weight is often a primary driver, internal structures such as spars, ribs and stringers are formed from very thin-walled plates. These structures are subject to a combination of flexural, compression and shear loads, and as a consequence, buckling behaviour is often a design driver. While finite element analysis (FEA) is invaluable for detailed design there remains a strong need for reliable simple design formulae for initial sizing and optimization studies. These exist for isotropic and orthotropic plates and most are contained in the classic texts of Timoshenko (1936) and Lekhnitskii (1968). However, most composite plates exhibit some form of anisotropy, which renders orthotropic buckling formulae inaccurate at best and non-conservative at worst. Unfortunately, balanced symmetric plates that contain angle-plies, generally, exhibit flexural/twist anisotropy to some extent and no closed form simple design formulae exist for buckling loads of laminated plates that are flexurally anisotropic. Furthermore, it is well known that flexural/twist anisotropy reduces buckling loads of laminated plates subjected to compression loading (Chamis 1969; Nemeth 1986). Nemeth (1992*a*, 1995) has produced a set of comprehensive and validated design charts to assist designers basing the work on a numerical solution to an analytical model. Therefore, the objective of the current work is to present approximate closed form solutions that prove useful for preliminary design purposes. In this paper, the effects of uniform compression buckling of long flexurally anisotropic plates with simply supported edges is considered. By presenting such results engineers may immediately assess the impact of flexural/twist anisotropy on buckling loads.

Firstly, the buckling analysis is derived and approximate closed form expressions for the buckling load in terms of anisotropic parameters are presented. At the appropriate stages of development, relevant accomplishments by other authors are discussed in context with the present work. Secondly, a discussion of the implications for design is given.

## 2. Derivation of model

### (a) Equilibrium equations

The governing expression for static equilibrium of a flat plate subjected to compression loading, as given, for example, by Timoshenko (1936) is(2.1)where *M*_{x},*M*_{y} and *M*_{xy} are bending and twisting moments, respectively, and *N*_{x} represents the compression load.

To introduce anisotropy into the problem, the constitutive behaviour of the material is considered.

### (b) Constitutive equations

The laminate is assumed to be balanced and symmetric such that(2.2)hold. The *D*_{ij} terms are flexural rigidities given in terms of the material invariants, *W*_{1}–*W*_{5}, and lamination parameters, *ξ*_{1}–*ξ*_{12} (e.g. Onoda 1985; Fukunaga 1991), as(2.3)The lamination parameters may be calculated from the following integrals:(2.4a)where *h*_{i} is the distance of a particular ply surface from the mid plane, *θ* is the fibre angle of a particular ply and *t* is the thickness of the laminate. Note that the flexural/twist coupling is given in terms of *D*_{16} and *D*_{26}.

The material invariants are linear functions of the ply stiffness properties, *Q*_{ij}, as given, for example, by Jones (1999),(2.4b)The curvatures are(2.5)The immediate goal is to reduce the governing equations into our unknown displacement *w*, by substituting equation (2.5) into equations (2.2) and (2.1). The governing equation becomes(2.6)where an expression of this kind was first published in the dissertation by Gehring (1860), written in Latin and apparently independently, in English, by Lord Kelvin (1867). One method of obtaining a buckling load is to represent the transverse displacement, *w*, as a complete series of linearly independent basis functions of variable amplitude. A Fourier series is a widely used example. However, one reason an accurate closed form solution has not previously been produced before is the reliance on using many (more than three) terms in the Fourier series for representing transverse plate deflections. Once, more than a few terms occur in such a series; the resulting expressions become unwieldy and intractable in a closed form. Having observed FEA buckling patterns, a one-term expression for transverse deflection that appears to capture the majority of deformation behaviour in a useful way for conceptual design purposes, is proposed. The effect of exactly satisfying boundary conditions is assessed by comparison with FEA and the charts obtained by the exact solutions for the infinitely long plate (Nemeth 1992*a*, 1995). It transpires that exactly satisfying boundary conditions are secondary in effect for thin and long plates.

The proposed transverse displacement for a plate of length *a* and width *b*, at the onset of buckling, is(2.7)and is effectively a one-term displacement series that has the distinct advantage over an infinite series in that an approximate closed form solution becomes viable. The terms *m* and *n* reflect the number of longitudinal and transverse displacements while the *k* term reflects the influence of flexural/twist anisotropy on buckling patterns, i.e. converts a nodal line that is perpendicular to the long edges into one that is slanted (figure 1). The maximum amplitude of the mode is *w*_{0}. It is now convenient to digress slightly and discuss the origins of using equation (2.7) as a waveform in buckling analysis. Timoshenko (1913) was the first person to propose such an expression and used it for approximating the buckling mode of a web of a I section that is subject to shear loading. He realized that such a mode could not exactly satisfy the boundary condition of simple support along the longitudinal edges, since the zero-bending moment condition is not exactly satisfied. However, he appeared to realize that this was not necessarily important and continued his analysis. He obtained a simple expression for the shear-buckling load of an infinitely long isotropic plate and published the result in 1913. Another 10 years passed before Skan & Southwell (1924) calculated an exact expression from successive approximations until convergence was obtained. Timoshenko's result was 6% too great, yet the simplicity of the result justifies the approximate solution in design.

Later, Sekerzh-zen'kovich (1931), following Timoshenko's approach, and using equation (2.7) as the modeform, determined an approximate shear buckling load for an orthotropic plate in terms of undetermined wavenumbers *m* and *n*. Balabuch (1937) extended Sekerzh-zen'kovich solution to anisotropic plates and considered combined loading. Although, the framework was set up for developing closed form solutions, buckling loads were once again left in terms of unknown wavenumbers *m* and *n*. Independently of Balabuch's work, Thielemann (1950) derived similar equations for buckling loads of anisotropic plates. Surprisingly, these works have not received great attention by other researchers. Perhaps it is because they do not approximate their solutions further to a closed form that they have not received greater use by practising engineers and the fact that the boundary conditions are not obeyed exactly may have hindered acceptance by the academic community. Finally, the work of Green & Hearmon (1945) is of merit. They attempted to solve a similar problem using exponential functions, as done earlier by Skan and Southwell on shear buckling of isotropic plates, and once again used the method of successive approximations to obtain numerical results. Convergence was found to be relatively slow.

### (c) Development of model

Initially, a similar approach to that taken by Balabuch and Thielemann is followed before simplifying procedures using knowledge of material parameters is applied. Substituting equation (2.7) into equation (2.6) and making the necessary differentiations and reductions gives(2.8)with and .

Note that to obtain equation (2.8) above the following non-conventional reduction method was used. After substituting equation (2.7) into equation (2.6) and gathering similar trigonometric variables leaves coefficients of just two terms:Conventional practice would now adopt either Rayleigh–Ritz or Galerkin procedures, both of which involve integration. As integration is an averaging process in itself, it appears simpler, in this instance, to make use of the average values of the trigonometric term over the buckling wave. Dividing through by the double sine term leaves a double cotangent term that has an average value of zero and as such, terms with this as a coefficient have been eliminated1 from equation (2.8). Furthermore, equation (2.8) simplifies to the standard expression, originally developed by Huber (1929), and given, for example, by Jones (1999) for special orthotropy. In general, the flexural/twist anisotropy terms, *D*_{i6} that appear in the skew parameter, *k* and coefficients *δ* and *γ* disappear under conditions of special orthotropy. It is possible to determine the minimal buckling load with respect to *λ*_{m}, *λ*_{n} and *k*. Obviously, *n*=1, since equation (2.8) is monotonic in *n* for all known materials. To determine values for *m* and *k* such that a closed form solution to equation (2.8) becomes viable is trickier. However, using a little insight simplifies the solution procedure. Firstly, equation (2.8) is re-cast in non-dimensional form as(2.9)using the non-dimensional parameters recommended by Nemeth (1992*b*),(2.10)and also with(2.11)The material parameters, *α*, *β*, *γ* and *δ* are not independent but they do capture the effects of orthotropy and flexural anisotropy on buckling resistance. This is particularly evident when it is realized that(2.12)for known composite materials (Nemeth 1992*a*,*b*; Weaver 2003). These values provide limiting bounds to the values of coefficients in equation (2.9) that is essential in rationalizing the curtailment of small terms in the solution process, and hence facilitate the development of approximate closed form solutions.

It is possible to obtain a measure of the angle of skew (angle of nodal lines) due to anisotropy (angle *ψ* in figure 1). This is found by setting the transverse displacement, *w* to zero in equation (2.7) and finding an expression for *ψ*. Firstly, note that *ψ* is given by(2.13)and that by setting *w*=0 gives(2.14)where *n*′ is an integer. Differentiating equation (2.14) with respect to *x* and rearranging gives(2.15)as an expression to determine the angle of skew, *ψ*. As such, *k*′ may be viewed as a non-dimensional measure of the angle of skew.

To obtain an approximate closed form solution it is noted that , and that the differential d*K*_{x} is given by(2.16)To obtain a buckling load, d*K*_{x}=0 or(2.17)Considering the second term, since then as a necessary condition. However, this leads to as a necessary condition in the first term. By differentiating equation (2.9) and ensuring that and leaves(2.18)with(2.19)These solutions are soluble but leave the expressions for *K*_{x} in an inelegant form. It is easier to drop all small terms in the coefficients *a*–*i*. Noting that the terms *γ* and *δ* have values relatively close to zero, it is proposed to drop terms of 0(*δ*^{2}, *γ*^{2}) or higher as a first approximation. Furthermore, as it is anticipated that *k*′ is a function of *δ* and *γ* (i.e. the nodal slope is expected to be a function of flexural anisotropy), the value of the *c* approximates to unity. Then(2.20)is a solution. This is the expected result for an orthotropic plate as originally found by Huber (1929). Therefore, we have assumed that flexural/anisotropy does not directly affect the value of the longitudinal wavelength. These assumptions are validated later in §3. Furthermore, the expression for *k*′ is found by neglecting terms of *k*′^{2} or higher in equation (2.19), and is(2.21)noting that *k*′ is a linear function of *δ* and *γ*. A simple expression for *K*_{x} is found by substituting equations (2.20) and (2.21) back into equation (2.9) to give:(2.22)The orthotropic solution, originally found by Huber (1929) is regained by setting *γ*, *δ* =0. As the buckling load is an even function of flexural/twist anisotropy (*γ* and *δ*) it is apparent that this form of anisotropy always reduces axial compression buckling loads and confirms the conclusion of other researchers (Chamis 1969; Nemeth 1986; Grenestedt 1991). Furthermore, it appears that *δ* is three times more effective at reducing buckling loads as is *γ*. This has important consequences for lay-up selection and optimization studies.

If it is deemed necessary, it is also possible to obtain a more accurate second approximation for the buckling coefficient. Firstly, combining equations (2.18*a*) and (2.9) gives a simpler form of *K*_{x} as(2.23)Also, note that the last four terms in *c* are small compared to the first term (unity) when *k*′ is given by equation (2.21), because they are second order in the anisotropy parameters and are also divided by (*β*+3). Thus, a binomial expansion can be used to determine an expression for *η*^{2} that is obtained by taking the square root of minus *c*. This gives(2.24)noting that if *η* is required, equation (2.18*a*) is preferable rather than taking the square root of the less accurate equation (2.24). Next, by combining equations (2.23) and (2.24), and simplifying gives a second approximation to the buckling coefficient as(2.25)A second approximation for the angle of nodal lines is obtained by substituting the value of *η*^{2} from equation (2.24) into the second expression in equation (2.18) and solving the resulting cubic expression to determine a second approximation for *k*′ (see appendix A). If necessary, a third approximation for the buckling coefficient is found by using the second approximation for *k*′ and the third approximation for *η*. As will be shown later, this solution has converged for the entire range of elastic properties and gives results that compare to an accuracy of greater than 0.3% with FEA.

The above model is expected to predict the buckling response of flexurally anisotropic plates for finite lengths in which . At such lengths the cusps on the buckling curve (*K*_{x} versus ) converge to constant value (e.g. Stavksy & Hoff 1969). For shorter rectangular plates boundary conditions at the loaded edge is expected to raise buckling loads. Shear deformation is expected to increasingly depress buckling loads for *b*/*t* or *l*/*mt* ratios less than 100.

## 3. Validation

Equations (2.22) and (2.25) for the buckling coefficient reduce to the familiar expression for the orthotropic laminate (*δ*=*γ*=0):(3.1)(Huber 1929). Since equations (2.22) and (2.25) are even functions of *δ* and *γ*, and these in turn, are odd functions of *θ* it is expected that replacing *θ* by −*θ* in the stacking sequence does not affect the magnitude of the buckling coefficient. This is to be expected from flexurally anisotropic laminates. However, the sign of the angle of nodal lines does change since, *k*′ is an odd function of ply angle, *θ*. In general, it is expected that the current analysis will result in an overestimation of buckling loads. This is partially due to the non-satisfaction of boundary conditions on bending moment along the plate's longitudinal edges and the incomplete satisfaction of equation of motion (equation (2.6)). As such, according to Rayleigh's Principle (e.g. Southwell 1936) the model has less degrees of freedom than necessary and is, in a sense, stiffer than reality. It is now necessary to quantify this induced error. As it is known that the induced error vanishes for orthotropic materials and is maximal for lay-ups with greatest flexural anisotropy, a qualitative measure of error is expected to be given by the magnitude of angle of skew or *k*′, to be more exact. In turn, *k*′ is a monotonic function of *δ* and *γ*, themselves direct measures of flexural anisotropy. It follows that accuracy is expected to decrease for increasing magnitudes of *δ* and *γ*. The present validation study is, therefore, made on a material with large amounts of flexural/anisotropy. For balanced, symmetrical laminates made from existing materials it may be shown that *δ* and *γ* never exceed a magnitude of 0.7 and then only for a [±*θ*]_{s} laminate made from P100/AS3501 (*E*_{11}=369 GPa, *E*_{22}=5.03 GPa, *G*_{12}=5.24 GPa and *ν*_{12}=0.31) and with *θ* close to 40° or 60°. This material also gives the largest value of *β* with a value of *β*=2.75 for *θ*=45°.

The models were compared with Nemeth's buckling charts (Nemeth 1992*a*) and with respect to FEA using ABAQUS Standard (Hibbitt, Karlsson and Sorenson Inc. 2003). Eight noded isoparametric shells (S8R) proved reliable elements and gave accuracies of better than 0.2% when compared with orthotropic solutions. Boundary conditions were such that the loaded edge remained straight and were free to rotate and the longitudinal edges were also kept straight and allowed to rotate and translate laterally (due to Poisson's ratio effects). To eliminate mesh effects, sufficient care was taken to ensure that at least 10 elements per half wavelength were used. This was done by modifying the length of the plate accordingly, for each lay-up. For instance, a plate with zero-degree fibres along the plates major axis had longer wavelengths than a plate with fibre angle at 90°. To ensure that the FEA gave solutions for a long plate, lengths that gave at least eight longitudinal half wavelengths were used. Once a viable solution was found two neighbouring solutions were also found, each for plates that were 20% longer and 20% shorter than the original. The solution was accepted if these neighbouring results were within 1% of the original. Finally, to minimize transverse shear deformation effects *b*/*t* ratios were chosen to exceed values of 100.

Firstly, the FEA analysis is validated with respect to Nemeth's buckling charts as the latter are already validated (Nemeth 1992*a*). Both methods give practically identical values for the buckling coefficient. That is, except for , when FEA gives approximately 2% lower values of *K*_{x}. This small difference probably arises from the curtailment of the infinite Fourier series at a limited number of terms, in the analysis that gives those buckling charts.

## 4. Results

In accordance with usual practice the buckling coefficient is plotted as a function of ply angle for balanced, angle-ply laminates in figures 2 and 3 for the P100/AS3501 material. These results represent those buckling calculations obtained from the third approximation or FEA as they are both indistinguishable from each other. However, *K*_{x} is not an appropriate measure for optimizing buckling loads because *K*_{x} is itself a function of lay-up, via *D*_{11} and *D*_{22}. Therefore, a more objective measure is required. To accomplish this one further non-dimensional parameter is necessary,(4.1)so that *N*_{x}/*N*_{x(iso)} is used as a more informative measure of buckling performance. It is a relative term as it normalizes the buckling load with respect to the quasi-isotropic lay-up of the same material. As such, its value reflects the relative performance of elastic tailoring with respect to the quasi-isotropic lay-up:(4.2)This represents a practical measure for design purposes, since the quasi-isotropic lay-up is frequently used and *N*_{x}/*N*_{x(iso)} gives a direct measure of differing lay-ups on buckling behaviour. It is given by(4.3)noting that(4.4)for isotropic materials, noting that, for this case, *β*=1, *δ*, *γ*=0 are substituted into equation (3.1).

The effect of ply orientation and stacking sequence on *N*_{x}/*N*_{x(iso)} is shown in figure 3 and results are compared directly with FEA in tables 1–3. Unidirectional laminates have the largest values of *δ* and *γ* and as such provide the most testing lay-ups to model buckling accurately. Their results are shown in table 1. Comparing buckling coefficients obtained by the approximate methods developed herein with FEA gives very good correlation. Indeed, the third approximation estimates have converged to those solutions obtained by FEA and as such may be considered as exact in the sense that it is accurate to 3 significant figures. For the second approximation (equation 2.25), only overestimates buckling loads by a maximum of 2.5%. Indeed, for the vast majority of laminates and material systems results for the second approximations are indistinguishable (less than 1% in error) from exact results (FEA or Nemeth's charts). Furthermore, for laminates with , results from the first approximation are within 2% of those of the second approximation. Balanced laminates are more commonly used in design and have smaller magnitudes of *δ* and *γ* as reflected by the data in table 2. Finally, orthotropic angles ply laminates are detailed in table 3 where all data is exact since *δ*=γ=0.

The magnitude of the skew angle for the various laminates are also listed in tables 1–3. The greater the difference from 90° represents the greater effect of flexural anisotropy. Unidirectional laminates with fibre angles close to 50° have the greatest skewed nodal lines that reflect the largest magnitudes of *δ* and *γ*. Figure 4*a* shows the contours of superimposed transverse deflection for the unidirectional 45° laminate and as such represents a plate with a high degree of flexural anisotropy. Note the large degree of skewing of nodal lines under compression loading. Figure 4*b* shows a small longitudinal section of the same plate highlighting the mesh and noting that only contours of magnitude up to 10% of the maximum transverse displacement are shown. This shows the straightness of the nodal line and hence, the validity of equation (2.7) in modelling the buckling pattern. Also note that the region of plate that is not satisfied by this waveform is very small. This is shown by the very small region close to the long edges, where the nodal lines sharply changes direction to meet the boundary perpendicularly so as to satisfy the zero-boundary moment.

## 5. Conclusions

Approximate closed form expressions have been developed to assess the effect of flexural anisotropy on buckling loads of long rectangular plates with simply supported sides subject to compression. The models were compared with detailed FEA for plates of extreme orthotropy and flexural anisotropy. The third level of approximate calculations converged to results that are indistinguishable from FEA. The buckling load expressions from the first and second approximations have broad ranging applicability to most common laminates in use. The expressions have been validated for design purposes and may assist in lay-up optimization procedures and the design of supporting structures such as the angular placement of ribs and stringers.

## Acknowledgments

The author wishes to thank Mike Nemeth for his encouragement and technical insight.

## Appendix

There is a closed form solution for the cubic polynomial (equation (2.18*b*)) in . It is found from taking the real root from(A1)with(A2)noting that there are six *z* roots in three sets of complex conjugate pairs.

An alternative approximate solution to *k*′ is found as follows.

A more accurate value of *k*′ is found by assuming that *k*′ differs from by a small amount, such that(A3)with *e*≪1. Then for all integers, *p*,(A4)is a reasonable approximation. Substituting equation (A 3) into the second expression of equation (2.18), neglecting nonlinear (i.e. small) terms in *e* and rearranging givesEquations (A 3) and (A 4) combine to give a second approximation for *k*′. From which a third approximation for *η* may be from equation (2.18*a*). Finally, a third (and converged) solution for *K*_{x} is found directly from equation (2.9).

## Footnotes

↵Note that the same result could have been obtained using Galerkin's method directly as done by Balabuch and Thielemann. However, the present way appears to be more physically intuitive to the designer.

- Received November 17, 2004.
- Accepted July 14, 2005.

- © 2005 The Royal Society