## Abstract

An analytical model based on variational principles for a thin-walled stiffened plate subjected to axial compression is presented. A system of nonlinear differential and integral equations is derived and solved using numerical continuation. The results show that the system is susceptible to highly unstable local–global mode interaction after an initial instability is triggered. Moreover, snap-backs in the response showing sequential destabilization and restabilization, known as cellular buckling or snaking, arise. The analytical model is compared with static finite element (FE) models for joint conditions between the stiffener and the main plate that have significant rotational restraint. However, it is known from previous studies that the behaviour, where the same joint is insignificantly restrained rotationally, is captured better by an analytical approach than by standard FE methods; the latter being unable to capture cellular buckling behaviour even though the phenomenon is clearly observed in laboratory experiments.

## 1. Introduction

The buckling of thin plates that are stiffened in the longitudinal direction represents a structural instability problem of enormous practical significance [1–4]. Since they are primarily used in structures where mass efficiency is a key design consideration, stiffened plates tend to comprise slender plate elements that are vulnerable to a variety of different elastic instability phenomena. In this work, the classical problem of a panel comprising a uniform flat plate stiffened by several evenly spaced blade-type longitudinal stiffeners under axial compression, made from a linear elastic material, is studied using an analytical approach. Under this type of loading, wide panels with many stiffeners can be divided into several struts, each with a single stiffener. These individual struts are primarily susceptible to a global (or overall) mode of instability namely Euler buckling, where flexure about the axis parallel to the main plate occurs once the theoretical global buckling load is reached. However, when the individual plate elements of the strut cross section, namely the main plate and the stiffener, are relatively thin or slender, elastic local buckling of these may also occur; if this happens in combination with the global instability, then the resulting behaviour is usually far more unstable than when the modes are triggered individually [5].

In this work, the development of a variational model is presented that accounts for the interaction between the global Euler buckling mode and the local buckling mode of the stiffener, such that the perfect and imperfect elastic post-buckling response can be evaluated. A system of nonlinear ordinary differential equations subjected to integral constraints are derived and solved using the numerical continuation package Auto-07p [6]. It is indeed found that the system is highly unstable when interactive buckling is triggered; snap-backs in the response show a sequence of destabilization and restabilization derived from the mode interaction and the stretching of the plate surface during local buckling, respectively. This results in a progressive spreading of the initial localized buckling mode. This latter type of response has become known in the literature as *cellular buckling* [7] or *snaking* [8] and it is shown to appear naturally in the current numerical results. The effect is particularly strong where the rotational restraint provided at the joint between the main plate and the stiffener is negligible. As far as the authors are aware, this is the first time this phenomenon has been found analytically in stiffened plates undergoing global and local buckling simultaneously. Similar behaviour has been discovered in various other mechanical systems such as in the post-buckling of cylindrical shells [9], the sequential folding of geological layers [10], the buckling of thin-walled I-section beams [11] and columns [12].

In previous work, cellular buckling has been captured in physical tests for some closely related structures that suffer from local–global mode interaction [11,13]; it is revealed by the buckling elements exhibiting a progressively varying deformation wavelength. However, by increasing the rotational stiffness of the aforementioned joint, the snap-backs are moderated and the equilibrium path shows an initially smoother response with the snap-backs appearing later in the interactive buckling process. For these cases, the results from a purely numerical model, formulated within the commercial finite element (FE) package Abaqus [14], are used to validate the analytical model and exhibit very good correlation, leading to the conclusion that the physics of the system is captured accurately by the analytical model.

## 2. Analytical model

### (a) Modal description

Consider a thin-walled simply-supported plated panel that has uniformly spaced stiffeners above and below the main plate, as shown in figure 1 with the associated coordinate system. The panel length is *L*, and the spacing between the stiffeners is *b*. It is made from a linear elastic, homogeneous and isotropic material with Young's modulus *E*, Poisson's ratio *ν* and shear modulus *G*=*E*/[2(1+*ν*)]. A portion of the panel that is representative of its entirety can be isolated as a strut, shown in figure 1*b*,*c*, because the transverse bending curvature of the panel during buckling would be relatively small, particularly if *L*≪*n*_{s}*b*, where *n*_{s} is the number of stiffeners in the panel. The main plate (width *b* and thickness *t*_{p}) has upper and lower stiffeners (heights *h*_{1} and *h*_{2}, respectively, with equal thickness *t*_{s}) connected to it. The strut is loaded by an axial force *P* that is applied to the centroid of the cross section. Presently, however, although the model is formulated for the general case, the numerical results focus on the practically significant case where the stiffeners are only connected to one side of the panel, i.e. where *h*_{2}=*t*_{p}/2.

Moreover, the geometries are chosen such that global (Euler) buckling about the *x*-axis is the first instability mode encountered by the strut, before any local buckling in the main plate or stiffener occurs. The formulation for global buckling is based on small deflection assumptions, since it is well known that it has an approximately flat (or weakly stable) post-buckling response [15]. The study is primarily concerned with global buckling forcing the stiffener to buckle locally shortly after or even simultaneously. It has been shown in several works [11,12,16,17] that shear strains need to be included to model the local–global mode interaction analytically. For thin-walled components, rather than, for instance, soft core materials used in sandwich structures [17], Timoshenko beam theory gives sufficiently accurate results [11]. This can be applied to global flexural buckling within the framework of two degrees of freedom, known as ‘sway’ and ‘tilt’ in the literature [16]. These are associated with functions describing the lateral displacement *W* and the angle of inclination *θ*, respectively, defining the appropriate kinematics, as shown in figure 1*c*. From linear theory, it can be shown that *W*(*z*) and *θ*(*z*) can be represented by the following expressions:
*q*_{s} and *q*_{t} are the generalized coordinates of the sway and tilt components of global buckling, respectively. The corresponding shear strain *γ*_{yz} during bending is given by the following expression
*γ*_{yz} would be zero and *q*_{s} would be equal to *q*_{t}. The initial global buckling displacement is assumed to put the stiffener into extra compression and it therefore becomes vulnerable to local buckling. The direction of global buckling is crucially important; if the stiffener goes into extra compression and buckles locally, then it can lead to excessive deflection. This can induce plasticity within the stiffener, leading to the *tripping* phenomenon that can lead to catastrophic failure of panels [2,18]; currently, however, linear elasticity is assumed throughout.

The stiffener is taken to be a flat plate (blade-type) and the local kinematics of it require careful consideration. A linear distribution in *y* for the local in-plane displacement *u*(*y*,*z*) is assumed owing to the application of Timoshenko beam theory. The tips of the stiffeners have free edges, but a rotational spring with stiffness *c*_{p} is included to model the resistance to rotation about the *z*-axis that the joint provides between the stiffeners and the main plate. Hence, if *c*_{p}=0, then the stiffeners are effectively pinned with the main plate. In the generic case (*c*_{p}>0), however, the assumed out-of-plane displacement along the width of the stiffener *w*(*y*,*z*) can be approximated by a function that is a summation of trigonometric and polynomial terms. Figure 2 shows the form of the functions *w* and *u* with the algebraic expressions being given thus
*B*_{0}, *B*_{1}, *B*_{2}, *B*_{3} and *B*_{4} being constants evaluated by solving for the appropriate boundary conditions for the stiffener. The choice of shape for *f* reflects the possibility of the slope at the junction between the main plate and the stiffener being small in the less compressed zone of the stiffener. For *w*(*y*,*z*) is found to be

In this study, both perfect and imperfect panels are considered. A displacement *W*_{0} corresponding to *W* and a rotation *θ*_{0} corresponding to *θ* are introduced; the respective expressions are
*q*_{s0} and *q*_{t0} being the amplitudes defining the initial global imperfection. A local imperfection *w*_{0}(*y*,*z*)=*f*(*y*)*w*_{0}(*z*), see figure 2*a*, is also introduced, the form of which is derived from a first-order approximation of a multiple-scale perturbation analysis of a strut on a softening foundation [19]:
*z*=[0,*L*] and *w*_{0} is symmetric about *z*=*η*. This form for *w*_{0} has been shown to be representative for local–global mode interaction problems in the literature [20]. It enables the study of periodic and localized imperfections, the latter being shown to be relatively more severe for related problems [21,22].

### (b) Total potential energy

The governing equations of equilibrium are derived from variational principles by minimizing the total potential energy *V* of the system. The total potential energy comprises the contributions of global and local strain energies of bending, *U*_{bo} and *U*_{bl} respectively, the strain energy stored in the ‘membrane’ of the stiffener *U*_{m} arising from axial and shear stresses, and the work carried out by the external load *W*_{0} and *w*_{0}, are stress-relieved. This implies that the elemental moment *M* and hence the bending energy both drop to zero in the unloaded state (figure 2*d*). The global bending energy accounts for the main plate and is hence given by
*z* and *x*-axis. The local bending energy, *U*_{bl}, accounting for the stiffener only, is given thus
*y*, *D* is the plate flexural rigidity given by *Et*^{3}_{s}/[12(1−*ν*^{2})] and
*F* is an example function. The compressive side of the panel only contributes to the local bending energy once global buckling occurs. The membrane energy *U*_{m} is derived from considering the direct strains (*ε*) and shear strains (*γ*) in the stiffener. The global buckling distribution for the longitudinal strain *ε*_{z} can be obtained from the tilt component of displacement from the global mode, thus
*ε*_{zp}=−Δ. The combined expressions for the direct strains for the top and bottom stiffeners *ε*_{zt} and *ε*_{zb}, respectively, including local and global buckling are
*U*_{m} is given by
*U*_{d} is the contribution from direct strains, which is given by the terms multiplied by the Young modulus *E*, whereas *U*_{s}, the contribution arising from the shear strains, is given by the terms multiplied by the shear modulus *G*. The transverse component of the strain *ε*_{y} is neglected because it has been shown that it has no effect on the post-buckling stiffness of a long plate with three simply-supported edges and one free edge [23]. The total direct strain energy *U*_{d} is therefore
*U*_{s} requires the shear strain *γ*_{yz}, which is also modelled separately for the compression and the tension side of the stiffeners. The general expression for the shear strain *γ*_{yzi} in the stiffeners is thus
*U*_{s} is
*F* is evaluated at *V* is the work done by the axial load *P*, which is given by
*V* is given by the summation of all the strain energy terms minus the work done, thus

### (c) Variational formulation

The governing equilibrium equations are obtained by performing the calculus of variations on the total potential energy *V* following a well-established procedure that has been detailed in Hunt & Wadee [16]. The integrand of the total potential energy *V* can be expressed as the Lagrangian (*V* is given by
*V* must be stationary, hence the first variation *δV* must vanish for any small change in *w* and *u*. Since *w* and *u*; these comprise a fourth-order and a second-order nonlinear differential equation for *w* and *u*, respectively. To facilitate the solution of the equations within the package Auto-07p, the variables are rescaled with respect to the non-dimensional spatial coordinate *w*/*L* and 2*u*/*L*, respectively. Note that the scalings exploit symmetry about the midspan and the equations are hence solved for half the strut length; this assumption has been shown to be acceptable for cases where global buckling is critical [20]. The non-dimensional differential equations are thus
*ψ*=*L*/*h*_{1} and *q*_{s} and *q*_{t} leading to the three integral equations in non-dimensional form
*q*_{s} and *q*_{t} before any interactive buckling occurs, i.e. when *w*=*u*=0. This link is assumed to hold also between *q*_{s0} and *q*_{t0}, which acts as a simplification by reducing the number of global imperfection amplitude parameters to one; this relationship is given by
*y*.

Linear eigenvalue analysis for the perfect strut (*q*_{s0}=*q*_{t0}=*A*_{0}=0) is conducted to determine the critical load for global buckling *P*^{C}_{o}. By considering that the Hessian matrix *V*_{ij}, thus
*q*_{s}=*q*_{t}=*w*=*u*=0, the critical load for global buckling is obtained
*P*^{C}_{o} reduces to the expected classical Euler column buckling load.

## 3. Numerical study

The full system of equilibrium equations is difficult to solve analytically. The continuation and bifurcation software Auto-07p is thus used; it has been shown in previous work [11] to be an ideal tool to solve the equations for this kind of mechanical system. The solver is adept at locating bifurcation points and tracing branching paths as model parameters are varied. To demonstrate this, an example set of geometric and material properties is chosen thus: *L*=5000 mm, *b*=120 mm, *t*_{p}=2.4 mm, *t*_{s}=1.2 mm, *h*_{1}=38 mm, *h*_{2}=*t*_{p}/2, *E*=210 kN mm^{−2}, *ν*=0.3. The global critical load *P*^{C}_{o} can be calculated using equation (2.39), whereas the local buckling critical stress *σ*^{C}_{l} can be evaluated using the well-known formula *σ*^{C}_{l}=*k*_{p}*Dπ*^{2}/(*b*^{2}*t*), where the coefficient *k*_{p} depends on plate boundary conditions; limiting values being *k*_{p}=0.426 and *k*_{p}=1.247 for a long stiffener connected to the main plate with one edge free and the edge defining the junction between the stiffener and the main plate being taken to be pinned or fixed, respectively. For the panel selected, table 1 shows that global buckling is critical, and that the stiffener could be the next in line to buckle locally, because its critical stress is much less than that of the main plate.

### (a) Analytical model results

To solve the governing equations, the principal parameters used in the numerical continuation procedure in Auto-07p were interchangeable although generally *q*_{s} was varied for computing the equilibrium paths for the distinct buckling modes. However, the load *P* was used as the principal continuation parameter for computing the interactive buckling paths. The perfect post-buckling path was computed first from *P*^{C}_{o} and then a sequence of bifurcation points were detected on the weakly stable post-buckling path; the value of *q*_{s} closest to the critical bifurcation point (*C*) was identified as the secondary bifurcation point (*S*), see figure 3, and denoted as *q*^{S}_{s}. This triggered interactive buckling, and a new equilibrium path was computed. Typically, this reduced the load *P* with *w* and *u* becoming non-trivial as the interaction advanced. On the interactive buckling path, a progressive sequence of destabilization and restabilization, exhibited as a series of snap-backs in the equilibrium paths was observed—the signature of cellular buckling [7,10,11,12]. In figure 3, the first such snap-back is labelled as *T* with

For the example panel being considered, figure 4 shows equilibrium plots of the normalized axial load *p*=*P*/*P*^{C}_{o} versus the normalized buckling amplitudes of (*a*) the global mode of the strut and (*b*) the local mode of the stiffener. The graph in (*c*) shows the relative amplitudes of the global and local buckling modes. Finally, (*d*) shows the relationship between the sway and tilt amplitudes of global buckling, which are almost equal; indicating that the shear strain from the global mode is small but, importantly, not zero. Figure 5 shows a selection of three-dimensional representations of the deflected strut that include the components of global buckling (*q*_{s} and *q*_{t}) and local buckling (*w* and *u*) for *S* and cells *C*_{4}, *C*_{6}, and *C*_{8}, respectively. Figure 6 illustrates the corresponding progression of the numerical solutions for the local buckling functions *w* and *u*. The results for the example clearly show cellular buckling with the spiky features in the graphs in figure 4*a*–*c*. Moreover, the buckling patterns, as seen in figure 5, clearly show an initially localized buckle gradually spreading outwards from the panel midspan with more peaks and troughs.

In the literature [7,10], the Maxwell load (*P*^{M}) is calculated for snaking problems, which represents a realistic lower bound strength for the system. Presently, however, it is more complex to determine such a quantity, because the system axial load *P* does not oscillate about a fixed load *P*^{M} as the deformation increases. This is primarily due to the fact that, unlike systems that exhibit cellular buckling, such as cylindrical shells and confined layered structures [7], there are two effective loading sources as the mode interaction takes hold: the axial load *P*, which generally decreases and the sinusoidally varying (in *z*) tilt generalized coordinate *q*_{t}, which represents the axial component of the global buckling mode that generally increases. The notion of determining the ‘body force’, discussed in Hunt & Wadee [16], could be used as way to calculate the Maxwell load, but this has been left for future work.

Hitherto, the results have been presented for the case where the joint between the main plate and the stiffener is pinned (*c*_{p}=0). Figure 7 shows how the introduction of the joint stiffness affects the response. First of all, the value of *σ*^{C}_{l,s} is increased, because rotation at the joint is more restrained. This, in turn, increases the value of *q*^{S}_{s} and *q*^{T}_{s} at different rates. It can be seen that as *c*_{p} is increased the initial interactive mode that emerges from the secondary bifurcation *S* is more localized with a smaller wavelength (figure 7) and resembles the more advanced cells shown in figure 6 directly. Moreover, it is demonstrated that the snap-backs occur later in the interactive buckling process.

## 4. Validation

The commercial FE software Abaqus was selected to validate the results from Auto-07p where appropriate. Four-noded shell elements with reduced integration were used to model the structure. Rotational springs were also used along the length to simulate the boundary conditions at the joint of the stiffener with the main plate. An eigenvalue analysis was used to calculate the critical buckling loads and eigenmodes. The nonlinear post-buckling analysis was performed with the static Riks [24] method with the aforementioned eigenmodes being used to introduce the necessary geometrical imperfection to facilitate this. For the numerical example, the strut described in §3 has a rotational spring stiffness *c*_{p}=1000 *N* mm mm^{−1}. As discussed at the end of the previous section, this basically increases the gap between triggering the interactive buckling mode and the onset of cellular buckling behaviour because the initial local buckling mode is naturally more spread through the length; this smoothing and hence simplification of the response potentially allows for a better comparison with the analytical model.

Linear buckling analysis shows that global buckling is the first eigenmode, and the corresponding critical load and critical stress are *P*^{C}_{o}=1.64 kN and *W*_{0}) where *q*_{s0}=0.001 is combined with a local imperfection (*w*_{0}) where *A*_{0}=0.12 mm, *α*=8.0, *β*=35 and *η*=*L*/2, which matches the initial imperfection from Abaqus for the analytical model such that a meaningful comparison can be made. The graphs in (*a*,*b*) show the normalized axial load *p* versus the buckling amplitudes of the global mode and the local deflection of the stiffener, respectively. The graph in (*c*) shows the local versus the global buckling modal amplitudes; all three graphs show excellent correlation in all aspects of the mechanical response.

Figure 8*d* shows the local out-of-plane deflection profiles at the respective locations (i)–(iii) shown in figure 8*a*–*c*, the comparisons being for the same values of *q*_{s}. These, again, show excellent correlation between the results from Auto-07p and Abaqus, particularly near the peak load. A visual comparison between the three-dimensional representation of the strut from the analytical and the FE models is also presented in figure 9. The correlation with Abaqus, for a relatively high value of *c*_{p}, is seen to be excellent, but when *c*_{p} is decreased, the analytical model exhibits the characteristic snap-backs earlier (figures 4 and 7). However, the Abaqus solution for static post-buckling, which relies on assuming an initially imperfect geometry, does not (and indeed possibly cannot) obtain the progressive cellular buckling solution. This is not entirely surprising because the wavelength of the Abaqus solution seems to be fixed throughout the loading history. This demonstrates that for situations where buckling wavelengths are likely to change, the conventional method for analysing instability problems with static FE methods, where initial imperfections affine to a combination of eigenmodes are introduced and analysed using a nonlinear solver (e.g. the Riks method), may not be able to capture some experimentally observed physical phenomena. As an alternative, the authors attempted to model the problem as a dynamical system within Abaqus such that the snap-backs may be replicated more readily. However, computational issues in terms of slow and unreliable convergence stalled this study and it has been left for the future. Nevertheless, a recent study of a closely related structure [12] clearly showed that the analytical approach matches the physical response significantly more accurately than static FE methods, in terms of the load versus displacement behaviour, and the observed change in the post-buckling wavelength. The observation of precisely the cellular buckling response has been reported for the flanges of thin-walled I-section struts experimentally [25], which is very similar to the current case where *h*_{1}=*h*_{2} (figure 1). Hence, the authors are of the opinion that because cellular buckling is an observed physical phenomenon for those related structural components, the current results from the analytical model are also valid.

## 5. Concluding remarks

This work identified a form of cellular buckling in axially compressed stiffened plates. This arised from a potentially dangerous interaction of local and global modes of buckling, with a characteristic sequence of snap-back instabilities in the mechanical response. It is also found that the effect of the instabilities can be reduced somewhat by increasing the rotational stiffness of the joint between the main plate and the stiffeners. However, the mode interaction persists and the local buckling profile still changes wavelength, as has been found in recent related studies [11–13]. The increased rotational stiffness of the joint allows the model to be validated from a purely numerical formulation in FE software. However, earlier studies suggest that for smaller values of the aforementioned stiffness, the analytical model would be superior in modelling the actual physical response. This is in terms of the equilibrium response, the load-carrying capacity and the deformed shape, which is known from earlier research to exhibit the change in the local buckling wavelength. The latter point is crucial because the static FE solution seems to be unable to capture the changing wavelength of the post-buckling mode that has been observed in physical experiments. It is likely that in an actual scenario, the subcritical response in conjunction with the snap-backs are likely to lead to buckling cells being triggered dynamically. This has been observed during the *tripping* phenomenon [2] and the postulate is that the elastic interactive buckling found currently, if combined with plasticity, would promote this kind of sudden and violent collapse of a stiffened panel.

Currently, the authors are extending the analytical model by including the possibility of buckling the main plate in combination with the stiffener. In this case, the load-carrying capacity is likely to be reduced and it would be representative of a scenario where the cross section has uniform thickness. Moreover, an imperfection sensitivity study is being conducted to quantify the parametric space for designers to avoid such dangerous structural behaviour.

- Received February 3, 2014.
- Accepted May 1, 2014.

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.