An Efficient Method for Local Buckling Analysis of Stiffened Panels

WANG Xinwei，YUAN Zhangxian

1.State Key Laboratory of Mechanics and Control of Mechanical Structures，Nanjing University of Aeronautics and Astronautics，Nanjing 210016，P.R.China；

2.School of Aerospace Engineering，Georgia institute of Technology，Atlanta，GA 30332-0150，USA

Abstract: The local buckling of stiffened panels is one of possible failure modes and concerned by engineers in the preliminary design of lightweight structures. in practice，a simplified model，i.e.，a rectangular plate with elastically restrained along its unloaded edges，is established and the Ritz method is usually employed for solutions. To use the Ritz method，however，the loaded edges of the plate are usually assumed to be simply supported. An empirical correction factor has to be used to account for clamped loaded edges. Here，a simple and efficient method，called the quadrature element method（QEM），is presented for obtaining accurate buckling behavior of rectangular plates with any combinations of boundary conditions，including the elastically restrained conditions. Different from the conventional high order finite element method（FEM），non-uniformly distributed nodes are used，and thus the method can achieve an exponential rate of convergence. Formulations are worked out in detail. A computer program is developed. improvement of solution accuracy can be easily achieved by changing the number of element nodes in the computer program. Several numerical examples are given. Results are compared with either existing solutions or finite element data for verifications. it is shown that high solution accuracy is achieved. in addition，the proposed method and developed computer program can allow quick analysis of local buckling of stiffened panels and thus is suitable for optimization routines in the preliminary design stage.

Key words: local buckling;quadrature element method;elastically restrained;stiffened panel

0 introduction

Stiffened plate is one of the structural elements widely used in aeronautical and marine structures.Due to its thin-walled in nature，stiffened plate is susceptible to buckle. Thus the buckling behavior of stiffened plate is one of the major concerns by the structural designers［1-2］.

Depending on the geometry，the stiffened plate exhibits several different buckling modes. One of the common buckling modes in typical lightweight structures is the local buckling，i.e.，the skin buckling between stiffeners. in the preliminary design stage，a simplified model，i.e.，a rectangular plate with elastically restrained on its unloaded edges，is used and the Ritz method is then employed to obtain the local buckling load［3-6］. Either Saint Venant torsion bar［3］or infinite independent rotational springs［4］are employed to model the elastic constraint caused by the stiffeners. if only one unloaded edge of the rectangular plate is rotationally restrained，the model can be also used to obtain the local buckling load of the flange of H-beams［5］. With a properly selected plate width，the accuracy of the local buckling load can be further improved［6］.

Although the local buckling problem can be solved by conventional finite element method（FEM）without any difficulty，however，it is inconvenient to use the FEM during the preliminary de-sign stage，since the dimensions of the stiffened plate and stiffeners are not finalized and to be determined by the optimization design. This is perhaps the main reason why the Ritz method is widely used in practice［1-6］. To employ the Ritz method，however，the loaded edges of the rectangular plate are usually assumed to be simply supported. in experiment，the loaded edges are usually clamped. Thus a case-by-case dependent correction factor has to be used for comparisons between the Ritz solutions and experimental results. Therefore，alternative efficient methods should be resorted to solve the title problem.

The weak form quadrature element method（QEM）is one of the alternative efficient methods.The QEM was originally proposed by Striz et al.［7-8］. After further developments by Zhong et al.［9-10］，Xing and Liu［11］，and present authors'research group［12-13］，now the QEM has been projected by its proponents as a potential alternative to the conventional FEM［14-15］. QEM combines the generality of the conventional FEM with the accuracy of spectral techniques and thus can achieve an exponential rate of convergence［14］. The solution accuracy of the QEM can be easily adjusted by increasing the number of element nodes. Besides，performing a parametric study by the QEM is an easy task. Thus the method is especially suitable for local buckling analysis of stiffened plate during the preliminarily design stage.

According to literature survey［14-15］，the QEM has not been used to solve the buckling problem of thin rectangular plates with elastically restrained edges. Therefore，the main objective of present investigation is to present a new and efficient approach for the local buckling analysis of stiffened plates. To model the rotationally elastic constraint，novel torsion bar element is proposed. Besides，a rectangular thin plate element is developed. Formulations and solution procedures are worked out in detail. Numerical examples are given，and results are compared with either existing solutions or finite element data for verifications. Finally，some conclusions are drawn based on the results reported herein.

1 Formulations of the Torsion Bar and Thin Rectangular Plate Elements

1.1 Energy expressions

A stiffened plate under uni-axial compression is schematically shown in Fig.1（a）. Fig.1（b）shows the simplified model for local buckling analysis of the stiffened plate. The plate side lengths are denoted by a and b. A uniform plate thickness t is considered. Elasticity modulus and Poisson's ratio are denoted by E and μ. Cartesian coordinate system（x，y，z）is set at the middle plane of the plate，thus-t/2 ≤z ≤t/2. The origin of the coordinate system is located at the plate center. The elastic constraints along unloaded edges，caused by the stiffeners，are modeled by two torsion bars.

The strain energy of the thin rectangular plate is

where D is a 3×3 symmetric matrix，and κ the cur-vature vector defined as

Fig.1 Sketches of stiffened plate and simplified model

where subscripts“x”and“y”denote the partial derivatives with respect to the Cartesian coordinates x and y，i.e.，wxx=∂2w/∂x2，wyy=∂2w/∂y2，wxy=∂2w/∂x∂y，and w(x，y) is the deflection，respectively.

For the isotropic homogeneous plate，elements in matrix D are defined as

The work done by the axial compressive distributed force Nxis given by

Currently，two ways are available to model the elastic constraint. One way is using the Saint Venant torsion bar to model the rotational elastic constraint［3］. The strain energy of the torsion bar is given by

where G is the shear modulus，J the half of the torsion constant of the stiffener about the x axis［3］，and θ=∂w/∂y at y= -b/2 or y=b/2.

The other way is using infinite rotational springs to model the rotational elastic constraint［4］.The strain energy of the torsion bar is given by

where kris the rotational spring constant， θ=∂w/∂y at y= -b/2 or y=b/2. Note that the unit kris the same as the one of the force.

Since QEM is to be used，only essential boundary conditions are required and given.

Simply supported edge（S）

Clamped edge（C）

Rotationally elastic restrained edge（E）

1.2 Quadrature torsion bar element

Different from the conventional high order finite elements，non-uniformly distributed nodes are used to achieve an exponential rate of convergence.Several types of node are available［14］and Gauss-Lobatto-Legendre（GLL）nodes are used in this paper.GLL nodes are roots of the following equation

where PN-1(ξ) is the（N-1）th-order Legendre polynomial. Note that ξi∈[-1，1] and thus xi=aξi/2(i=1，2，…，N ).

The rotational angle θ of the N-node torsion bar element is assumed as

where lj(x) is the Lagrange interpolation function defined as

The stiffness matrix of the quadrature torsion bar can be explicitly given by using GLL quadrature and differential quadrature（DQ）law. if Eq.（5）is used，the elements in the stiffness matrix denoted by kb1ijare

where Hkis the weight of GLL quadrature，Aikand Akjare the weighting coefficients of the first-order derivative with respect to x. A short program is available in Ref.［14］to compute the abscissas and weights in GLL quadrature.

in Eq.（13），Aij(i，j=1，2，…，N ) can be explicitly computed by using the differential quadrature（DQ）law as［14-15］

if Eq. （6） is used，the stiffness matrix of quadrature torsion bar can be also explicitly given by using GLL quadrature and the DQ law. Note that the stiffness matrix is a diagonal matrix and its diagonal terms denoted by kb2iiare given by

1.3 Quadrature rectangular plate element

Let N be the number of node in either x or y direction. An N×N-node quadrature rectangular plate element is formulated. For assemblage considerations，GLL nodes are also used.

According to the criteria for the selection of displacement functions［16］，three different displacement functions are assumed for formulating the rectangular plate element，which is novel and different from the common way to formulate a C1compatible finite element.The three displacement functions are

where li(x) and lj(y) are the Lagrange interpolation functions and their definition is given by Eq.（11）；and hi(x) and hj(y) the Hermite interpolation functions and their definition can be found in Ref.［17］.

To satisfy the criteria for forming a C1compatible element［16］，Eq.（16）is used to compute wx，wyand wxy，Eq.（17）is used to compute wxx，and Eq.（18）is used to compute wyy. in this way，the quadrature thin plate element contains only (N2+4N ) degrees of freedom（DOFs），one DOF at all inner nodes，two DOFs at all boundary nodes，and three DOFs at four corner nodes.

With GLL quadrature and the DQ law，the plate stiffness matrix k can be explicitly given by

where Hi，Hjare the weights of GLL quadrature.

in Eq.（19），the strain matrix at an integration point(xi，yj)（i，j=1，2，…，N）can be explicitly given by

where δjkand δilare the Kronecker symbols，contains wij（i，j=1，2，… ，N），(wx)l1，(wx)lN，(wy)1kand (wy)Nk，(l，k=1，2，…，N )；Axiland Ayjkthe weighting coefficients of the first-order derivative with respect to x or y obtained by using Eq.（16），and their explicit formula is given by Eq.（14）；the weighting coefficients of the second-order derivative with respect to x obtained by using Eq.（17）；andthe weighting coefficients of the second-order derivative with respect to y obtained by using Eq.（18）. Explicit expressions for computingandare also available and can be found in Refs.［15，17］. it is worth noting that three different displacement equations are used to compute Axil，Ayjk，and.

The geometric matrix g is also obtained by using the GLL quadrature.its non-zero elements are

1.4 Solution procedures

For the problem to be solved，one N×N-node rectangular plate element and two N-node torsion bar elements are used. The assemblage procedures are exactly the same as the conventional FEM.

After assemblage and applying essential boundary conditions，following matrix equation is resulted.

where δecontains non - zero derivative DOFs at boundary nodes and is to be eliminated，wrcontains non-zero displacement DOFs and is to be retained.

After eliminating δe，Eq.（22）can be rewritten as

Solving Eq.（23）by a generalized eigen-solver yields the eigen-values. The lowest eigen-value is the local buckling load of the stiffened plate.

2 Numerical Results and Discussion

A FORTRAN computer program is developed. To verify the program and show the excellent performance of the QEM，the buckling of rectangular plates（a/b=1 and 5）with all edges simply supported and subjected to uni-axial compression is analyzed first. Results with N varying from 5 to 13 are listed in Table 1. Exact solution is also included for comparisons.

Table 1 Comparison of buckling coefficient λ

in Table 1，λ=Nxb2/(π2D0). it is seen that the rate of convergence of the proposed thin plate element is high.One 7 × 7-node element can yield accurate buckling load up to three decimal places for the square plate. For a rectangular plate with aspect ratio of 5，one 13 × 13-node element can also yield accurate buckling load up to three decimal places.Due to a higher wave number in the buckling mode，it is expected that larger number of node is required for the rectangular plate than the square plate since only one rectangular plate element is used in the analysis.

For demonstrations and comparisons，two different stiffened panel geometries are considered.One has four equally spaced stiffeners shown in Fig.1（a），and the corresponding simplified model，as the one shown in Fig.1（b），has a dimension of a =700 mm and b=280 mm. The other case has five equally spaced stiffeners. The corresponding simplified model，as the one shown in Fig.1（b），has a dimension of a=700 mm and b=140 mm.The skin thickness t is 1 mm in both cases.

The cross section of the stiffener is a rectangular plate with height（hw）20 mm and thickness（tw）3 mm. The material of both skin and stiffener is aluminum alloy. Elasticity modulus is 72 GPa and Poisson's ratio is 0.33.

The half of the torsion constant of the stiffener about the x axis is calculated as［18］

When the stiffened plate is in the state of local buckling，the stiffeners do not buckle and thus the unloaded edges can be treated as simply supported（S），clamped（C），or rotationally elastic restrained（E）ones depending on the cross-sectional shape of the stiffeners. For comparisons with existing data，the loaded edges are either simply supported or clamped，and the unloaded edges are rotationally elastic restrained. in other words，an SESE plate or a CECE plate is analyzed.

Fig.1（b）is modeled by one 13 × 13 -node quadrature plate element and two 13-node quadrature torsion bar elements. The local buckling stresses of the stiffened plates with four and five stiffeners，denoted by σcr（σcr=Ncrx/t），are listed in Table 2. Both the entire stiffened panel and the simplified model are analyzed by using finite element methods，and the corresponding results are denoted by σcrFulland σcrStrip. Existing data obtained by Ritz method are also included in Table 2 for comparisons. The percentage relative difference is calculated by

it is seen that present QEM yields the most accurate local buckling stress. Ritz method with different simplified model yields different results. Besides，the data obtained by the QEM are in excellent agreement with the ones obtained by the FEM if the same simplified model is analyzed. A model with 2 mm×2 mm mesh size is used in the finite ele-ment analysis to ensure the solution accuracy.

Table 2 Comparison of local buckling stress σcrof stiffened plates

Table 3 Local buckling stress σcrof stiffened plates with various boundary conditions

The local buckling mode is shown in Figs.2，3.The 700 mm×280 mm plate shown in Fig.2 has three half waves in x direction，and the 700 mm×140 mm plate shown in Fig.3 has seven half waves in x direction. They agree with the ones reported in Refs.［3，4］.

Fig.2 Buckling mode shape (700 mm×280 mm)

Fig.3 Buckling mode shape (700 mm×140 mm)

Similar to the conventional FEM，one of the advantages of the QEM can treat any boundary condition easily. The local buckling of stiffened plates with various boundary conditions is analyzed and the buckling stresses are listed in Table 3. For validations，the finite element data are also included. The edge numbers are shown in Fig.1（b）and the stress ratio α is defined as

where subscripts represent the boundary conditions of the loading edges，i.e.，edges 1 and 3.

From data listed in Table 3，it is seen that the smaller the spacing of the stiffeners，the closer the local buckling stress to the one of the plate with unloaded edges clamped. it is also observed that α depends not only on the boundary conditions of unloaded edges but also on the aspect ratio of the rectangular plate.

Fig. 4 shows the side view of the buckling mode shapes for the 700 mm×280 mm rectangular plate. The boundary conditions of loading edges are simply supported and clamped. Although both buckling mode shapes have three half waves in x direction，the amplitude of the waves are different. For the SESE plate，the amplitudes of two half waves close to the boundary are larger than the middle one.For the CECS plate，the amplitude of the middle half wave is higher than the other two and just opposite to the SESE plate.

Fig. 4 Side view of buckling mode shapes ( 700 mm×280 mm)

Fig. 5 shows the side view of the buckling mode shapes for the 700 mm×140 mm rectangular plate. Although the amplitudes have similar distributions as the ones shown in Fig.4，however，the numbers of half waves in x direction are quite differ-ent. The mode of the SESE plate contains seven half waves and is symmetric about y axis，but the mode of CECE plate contains only six half waves and is anti-symmetric about y axis.

Fig. 5 Side view of buckling mode shapes ( 700 mm×140 mm)

Since the boundary conditions in a real structural can never be ideal，Figs.6，7 show the variations of the local buckling stress with the torsion constant J for the stiffened plates with four and five stiffeners，respectively. it is worth noting that the scale of horizontal axis in Figs.6，7 is in logarithm.

Fig. 6 Variation of local buckling stress of the stiffened plates with four stiffeners

Fig. 7 Variation of local buckling stress of the stiffened plates with five stiffeners

From Figs.6，7，it is clearly seen that when J is less than 0.1 mm4，the elastic restrained boundary can be approximated as the simply supported boundary. On the other hand，the elastic restrained boundary can be approximated as the clamped boundary if J is greater than 100 000 mm4. Therefore，by adjusting the torsion constant J from 0 to infinity，it is able to simulate the real boundary，from the simply supported edges to the fully clamped edges. The method can be also used for treating the boundary conditions of loading edges if necessary. in summary，the presented method is suitable for analyzing local buckling of stiffened plates containing stiffeners with different cross-sectional shapes and with any combinations of boundary conditions.

it is also observed that the stress ratio α depends strongly on the aspect ratios of the rectangular plate and boundary conditions of the unloaded edges. Therefore，the stress ratio α should be determined case by case to compare the Ritz data obtained based on the simply supported loading edges with experimental results obtained using clamped loading edges. Since actual boundary conditions can be easily applied by the QEM，thus a correction factor α is not needed to compare the QEM data with experimental results. Besides，parametric study can be easily performed by the QEM.

it is worth noting that initial imperfections evidently exist in real structures and affect the buckling behavior. Several ways［19］to consider the initial imperfections exist in the finite element analysis and a static analysis，instead of an eigen-value analysis，is then performed. The same procedures can be adopted by the QEM since it is essentially a high order FEM. However，this is beyond the scope of the paper.

To use the simple eigen-value method for the buckling analysis，initial imperfections are not considered. The first or the second eigen-vector of the perfect stiffened plate multiplying a small factor may be used to introduce the initial imperfection for the buckling and post-buckling analysis of the stiffened panel［19］.

3 Conclusions

A new approach，called QEM，is presented for local buckling analysis of stiffened panels. A simpli-fied model，i.e.，a rectangular plate with elastically restrained along its unloaded edges，is established first and then analyzed by the QEM. Quadrature torsion bar and thin rectangular plate elements are developed. Several numerical examples are given. Numerical data are compared with either existing results or data obtained by FEM for verifications.Fine mesh is used in the finite element analysis to ensure the solution accuracy. A parametric study is also performed by the quadrature element method.

Numerical results show that the convergence rate of the QEM is high. One of the advantages of presented approach is that the improvement of the solution accuracy can be easily done by simply increasing the number of element nodes in the computer program. Besides，parametric studies can be also easily carried out by the written program.

Based on the results reported herein，one may conclude that the presented method can yield accurate local buckling load of the stiffened plates with minimum computation effort. The solution accuracy is higher than the existing ones. Therefore，the quadrature element method is suitable for optimization routines in preliminary design besides for structural analysis.