Buckling Analysis for Wind Turbine Tower Design: Thrust Load versus Compression Load Based on Energy Method

Abstract

Tubular steel towers are the most common design solution for supporting medium-to-high-rise wind turbines. Notwithstanding, historical failure incidence records reveal buckling modes as a common type of failure of shell structures. It is thus necessary to revisit the towers’ performance against bending-compression interactions that could unchain buckling modes. The present investigation scrutinises buckling performances of a cylindrical steel shell under combined load, by means of the energy method. Within the proposed framework, the differential equations to obtain dimensionless expressions showed the energy-displacement relations taking place along the shell surface. Furthermore, shell models integrated with initial imperfection have been embedded into finite element algorithms based on the Riks method. The results show buckling evolution paths largely affected by bending moments lead to section distortions (oval-shaped) that in turn change the strain energy dissipation routine and section curvature. The shell geometrical parameters also show a strong influence on buckling effects seemingly linked to a noticeable reduction of the shell bearing capacity during the combined loading scenarios.

1. Introduction

Tubular steel shell structures are the most practical configuration for wind turbine tower design, that, nevertheless, are experiencing increasing challenges largely due to the scale development of the wind industry [1,2]. Historical data reveal that tower collapse is a significant factor threatening wind energy harvesting, especially in extreme wind events, whereas most of the published research on wind turbine failures were focused on blade failures and power generation system malfunctions. According to the practical operation guidelines of wind turbines in extreme wind conditions, the wind turbines move to the parking mode for avoiding uncontrollable fast blades rotation that could destabilise the supporting tower structure through additional cyclic dynamic loads [1]. Fundamental cylindrical steel shell analyses [3,4] helps understanding collapse mechanisms therefore enables optimization designs of wind tower shell configurations, door openings, and further steel detailing. Noticeably, various buckling modes seem to develop during such load scenarios, ending up in severe tower collapse cases [5,6]. Tubular shells of wind turbine towers when standing still or operating are subject to combined compression, torsion, and bending from the thrust load associated to cross winds and among which, the compression force and bending moment are dominant.

Modern shell structures show remarkable efficiency particularly with regards to material performance, so that failure mechanisms mainly relate to uncontrolled geometrical distortions such as buckling [7,8,9,10]. Forensic studies of wind turbine tower collapse cases have revealed wave-shaped buckling patterns containing bulges and dents along the circumferential direction of the tubular shells, as well as creases and wrinkles along the longitudinal direction [11,12,13,14]. The number of waves both circumferential and longitudinal falls within specific ranges that have previously been observed in pure compression tests and theoretical models [15,16,17,18,19,20]. More broadly, buckling types seem controlled by specific parameters such as tower diameter, shell thickness, tubular shell section length, stiffening rings (type and location), door-opening area, type and size of connection components, and welding [1,2,5,6,11,12,13,14]. Those parameters could normalize for running parametrical analyses and further clarify the nature of buckling-related collapse cases of wind turbines recorded in the past. For that reason, the present investigation focuses on the study of the buckling phenomenon [7,8,9,10,15,16,17] of an unstiffened cylindrical steel shell subject to compression and bending with a particular focus on the progression of buckling that leads to the collapse of structural towers.

2. Buckling Theories and Equations

General

Buckling theory started to develop over 3 centuries ago, particularly across scientists drawn from engineering and mathematics. Notwithstanding, knowledge gaps remain, which is demonstrated by discrepancies between theoretical or numerical solutions and direct observations undertaken on physical structures. To an extent, those discrepancies relate to buckling and yielding (geometric and material softening, respectively) that follow ultimate load conditions—even though modern wind turbine design specifies material strengths that are well above the safe region for regular supporting structures while buckling tends to be limited by serviceability requirements for shell structures [21].

Buckling theories did not see a breakthrough in recent years, as research has mainly focused on narrowing the gaps for engineering application [22,23,24]. For example, Koiter [25] and Donnell [26] delivered remarkable works in classical buckling analysis, introduced eigenvalue calculations addressing the ‘two-surface’ model, and highlighted instability problems linked to geometrical imperfections. The classical buckling theory therefore aims to obtain the ultimate load before buckling instabilities occur, hence ensuring the material reaches the design strength. However, the progression of buckling is not totally unrelated to working design conditions, specifically the transient buckling stage. It is thus hypothesized that if such relationships between linear and nonlinear limit solutions integrate, the respective nonlinear response would be easier to predict, based on the linear result. Notwithstanding, this simplified approach depends largely on the specific structure under the nominated loading condition.

The buckling behaviour of cylindrical shells distinguishes between bifurcation and the limit buckling point, as the local compression increases. Figure 1 shows that the bifurcation tends to occur in the ideal elasticity state (curves 1 and 2) involving both linear and nonlinear prebuckling paths. In this figure, the vertical axis represents axial compression while ${\lambda }_{cr}$ is the ideal critical load factor of the structure. The horizontal axis represents the displacement, namely the end shortening of the shell structure. Bifurcation buckling paths could be symmetric or asymmetric and stable or unstable, which reflects in the secondary bifurcation buckling paths past the first bifurcation point, e.g., curves 2i, 2ii, and 2iii (and there are probably more subsequent bifurcation paths). Another key aspect of buckling is the limit point buckling or ultimate point along the smooth load–displacement curve, with or without snap-through behaviour. Shells with imperfections due to factory construction or on-site installation are prone to experiencing snap-through behaviour. That is a dynamic jump takes places after the limit point, which triggers a new equilibrium (perhaps more stable) path. The bifurcation path may emerge along the limit point buckling path such that the relation between the load and deformation curve fulfils a bifurcation point solution before reaching the limit point at load factor $m{\lambda }_{cr}$, as depicted along the curve 3–3ii. This is a problem of buckling equilibrium paths, based on buckling control equations, when considering large initial geometrical imperfections.

Before buckling fully manifests, the structure remains stable but as the critical buckling occurs, the load–displacement equilibrium path descends dramatically to a certain distance below the critical buckling load. This poses challenges to analytical or numerical solutions, since energy-method-based buckling control equations contain more than 3-ordered items, which, in light of the observations above, derive multiple equilibrium paths [27]. In classical approaches, initial imperfections feeds into to these equations, which can resolve with existing algorithms to provide a single equilibrium path. The present investigation follows such approximations to obtain the minimum point along the equilibrium path aiming at identifying the most realistic critical buckling load.

As mentioned above, algorithms to model buckling stability dominated research over the past decades [7,8,9,10,21,25,26,28,29,30,31]. The analysis and interpretation of buckling stability therefore underpins the analysis for real wind turbine tower. We note in that process that, based on the potential energy (classical) theory, energy is insensitive to the evolution path, and thus the stability equation derives from the total potential energy equilibrium, as depicted in Figure 2. In this chart, the origin point represents the state of the unloaded (hence undeformed) shell section. The following point is the initial stage field, as extreme to the vector $\overrightarrow{r_0}$. That state occurs past the application of gravity load and initial residual stress. Then, the basic stage field appears to describe the state under design loading, including upper structures installed during regular working periods. When external loads exceed the basic stage field, a trajectory following the vector $\overrightarrow{b_{+}}$, identifies the final stage. Noting that the final state accommodates arbitrary displacement or input energy components $\overrightarrow{r_0}$, $\overrightarrow{b}$, and $\overrightarrow{b_{+}}$ and perhaps other interval components. Each field therefore represents a set of possible stable points at a given stage, and stage evolution happens only when additional energy components exceed partial stability energy ranges.

For instance, one of the possible additional energy paths to the final stage field from the origin could be described as in Equation (1) while each vector possesses multiple set elements:

$$\overrightarrow{r}=\overrightarrow{r_0}+ \overrightarrow{b}+ \overrightarrow{b_{+}}$$

(1)

In the static analysis the kinetic energy component cancels, hence, the overall energy relates to the potential energy. In an arbitrary energy stage of a shell structure, the total potential energy ∏ can divide into internal energy and external energy:

$$\mathsf{\Pi }={\Pi }_{internal }+{\Pi }_{\begin{array}{c}external\\ \end{array}}={\Pi }_{i }+{\Pi }_{\begin{array}{c}e\\ \end{array}}$$

(2)

where, ${\Pi }_{internal }$represents the total potential energy stored by the shell elements, and the ${\Pi }_{external }$represents energy input from external load or environmental factors. Various possible energy paths have been applied in the above equations, and when δ∏ = 0, the neutral equilibrium state of the transition realises. When this principle translates to the eigenvalue problem involving the stiffness, geometries, and displacement matrices, in both its linear and nonlinear part, a typical simplified quadratic eigenvalue problem of the form unveils [32,33],

$$K_{T}V=K_e+{\lambda }^{*}\left(K_{gl}+K_{vl}\right)+{\lambda }^{*}{}^2\left(K_{gn}+K_{vn}\right)V=0$$

(3)

$$\left(K_e+{\lambda }^{*}{K}_{gl}\right)V=0$$

(4)

In Equations (3) and (4), $K_e$ represents the matrix of initial elastic stiffness, $K_{gl}$ represents the matrix of initial linear geometries, $K_{vl}$ represents the respective matrix of linear displacements, $K_{gn}$ represents the initial geometrical stiffness matrix in quadratic form, $K_{vn}$ represents the initial displacement matrix in the quadratic form, and $V$ is the nodal displacements of the selected element in the form of a vector. Equation (4) approximates the classical eigenvalue problem by ignoring the shell initial displacement. Incremental methods could solve the numerical problem through incremental load steps, whose equilibrium is checked iteratively. The solution of Equation (3) also relates to the critical load factor ${\lambda }_{cr}$ depicted in Figure 1; ${\lambda }^{*}$ is the eigenvalue solution of Equations (3) and (4), and is equal to ${\lambda }_{cr}$ when it gets to its minimum value; this factor could in turn establish a relationship with a membrane and bending components (${\lambda }_m \mathrm{and} {\lambda }_b$, respectively), as follows:

$${\lambda }_{cr}={\lambda }_b+m\times {\lambda }_m$$

(5)

The membrane part ${\lambda }_m$ is strongly linked to the total shell load bearing capacity, as shown in Figure 1. It is worth highlighting that, in that figure, the upper and lower load bounds delimit realistic buckling performances as an extension of linear-elastic bending strength and as part of membrane effects developed by cylindrical shells.

3. Load Analysis and Prebuckling Analysis of Wind Turbine Tower Shells

3.1. Axial Compressive Load

As previously stated, this research focused on the secondary bottom shell section of the wind turbine tower. The reason for this is that historical records and onsite inspection identify such a secondary bottom shell involved in the highest percentage of recorded buckling failures, apparently due to large induced bending and compression, which not always compare to predictions [6,11,14,30,34,35]. Buckling related failure distinguishes from other types of tower structural failure, such as bolt connection failure, flange failure [36], or blade-hitting damage, by having more dependency on initial imperfections or material non-homogeneity that causes uneven stress flow across the secondary bottom shell. Figure 3 shows a portion of the secondary bottom section, embedded within a cylindrical coordinate system, and all acting loads on its surface and boundaries.

In Figure 3, the circumferential load $N_{\theta }$ and moment $M_{\theta }$ do not vary with coordinate θ, because the circumferential displacement equals zero in this adopted linear model. Furthermore, by letting $u$, $\upsilon ,$ and $f$ be the displacements in the x, y, and z direction, respectively, one derives the following basic equilibrium equations [21,32,37]:

$$\{\begin{array}{l}\frac{d{N}_x}{dx}dx\cdot rd\theta +P_{x}rd\theta dx=0 \left(x-axis\right)\\ \frac{d{\theta }_x}{dx}dx\cdot rd\theta +N_{\theta }tan\theta d\theta +P_{r}dx\cdot rd\theta =0 \left(z-axis\right)\\ in which tand\theta \approx d\theta \\ \frac{d{M}_x}{dx}dx\cdot rd\theta -\left(Q_{x}rd\theta dx\right)=0 \left(y-axis moment\right)\end{array}$$

(6)

It follows that, by cancelling out similar items $dx$ and $rd\theta$ Equation (6) becomes:

$$\{\begin{array}{l}\frac{d{N}_x}{dx}+P_x=0\\ \frac{d\theta }{dx}+\frac{1}{r}{N}_{\theta }-P_r=0\\ \frac{d{M}_x}{dx}-Q_x=0\end{array}$$

(7)

$$N_x=\frac{Et}{1-{\nu }^2}\left(\frac{du}{dx}-\nu \frac{f}{r}\right)$$

(8)

where, $\nu$ represents the Poisson ratio (assumed as 0.3 for structural steel). Moreover, taking into consideration the stress–strain relationship as per Hook’s law while letting ε represent the strain, the following equations develops as:

$$N_{\theta }=\frac{Et}{1-{\nu }^2}\left({\epsilon }_{\theta }+\nu {\epsilon }_x\right)=-\frac{Et}{1-{\nu }^2}\left(\frac{f}{r}-\frac{du}{dx}\right)$$

(9)

Here we note that the relation between $M_x$ and displacement f can map from a plane slab once curved to a cylindrical shell:

$$M_x=-D\frac{d^{2}f}{d{x}^2}$$

(10)

$$M_{\theta }=\nu M_x$$

(11)

$$D=\frac{E{t}^3}{12\left(1-{\nu }^2\right)}$$

(12)

In which, D is the stiffness of the slab.

Whereas to eliminate the term $Q_x$ from Equation (7), we write,

$$\frac{d^2}{d{x}^2}\left(D\frac{d^{2}f}{d{x}^2}\right)+\frac{Eh}{r}f-\nu N_x-P_r=0$$

(13)

Equation (13) can be simplified by letting $P_r=0$, which seems sensible for a tower shell without radial pressure imposed to it:

$$D\frac{d^{4}f}{d{x}^4}+\frac{Eh}{r}f-\nu \frac{N_x}{r}=0$$

(14)

Equations (13) or (14) are basic differential equations for cylindrical shells, here representing a portion of the wind turbine tower, whose solution is [38]:

$$f=C_1{e}^{h_{1}x}+C_2{e}^{h_{2}x}+C_3{e}^{h_{3}x}+C_4{e}^{h_{4}x}$$

where the terms $C_i$ are constants of the integration while $h_i$ are roots to the following polynomial:

$$h^4+4{\beta }^4=0\hspace{1em}\beta =\frac{Et}{4D{r}^2}=\frac{3\left(1-v^2\right)}{r^{2}t}$$

(15)

If the dimension of $\beta$ is reciprocal to the length L, then we find:

$$\begin{gathered} h=\pm \beta \left(1\pm i\right)f=e^{-\beta x}\left(C_1{e}^{i\beta x}+C_2{e}^{-i\beta x}\right)+e^{\beta x}\left(C_3{e}^{i\beta x}+C_4{e}^{-i\beta x}\right) \\ f=e^{-\beta x}\left(C_{1}cos\beta x+C_{2}sin\beta x\right)+e^{\beta x}\left(C_{3}cos\beta x+C_{4}sin\beta x\right)+f\left(x\right) \end{gathered}$$

(16)

Equation (16) defines f as the displacement field in direction z and stands as the general solution to Equation (13).

At this end, strain energy has split in bending strain energy and stretching or membrane energy:

$$U=U_b+U_m$$

(17)

Based on the basic assumption of plane slab without the strain components ${\epsilon }_z$, ${\gamma }_{xz}$, and ${\gamma }_{yz}$, the accumulated energy of slab $U_s$ results in,

$$U_s=\frac{1}{2}{{\displaystyle \int }}_v\left({\sigma }_x{\epsilon }_x+{\sigma }_y{\epsilon }_y+{\tau }_{xy}{\gamma }_{xy}\right)dxdydz$$

(18)

where,

$$\begin{gathered} {\sigma }_x=-\frac{Ez}{1-{\nu }^2}\left(\frac{{\partial }^{2}f}{\partial x^2}+\upsilon \frac{{\partial }^{2}f}{\partial y^2}\right), {\epsilon }_x=-z\frac{{\partial }^{2}f}{\partial x^2} \\ {\sigma }_y=-\frac{Ez}{1-{\nu }^2}\left(\frac{{\partial }^{2}f}{\partial y^2}+\upsilon \frac{{\partial }^{2}f}{\partial x^2}\right), {\epsilon }_y=-z\frac{{\partial }^{2}f}{\partial y^2} \\ {\tau }_{xy}=-\frac{Ez}{1+{\nu }^2}\left(\frac{{\partial }^{2}f}{\partial x\partial y}\right), {\gamma }_{xy}=-2z\left(\frac{{\partial }^{2}f}{\partial x\partial y}\right) \end{gathered}$$

(19)

If the slab thickness is maintained constantly without thickness variation, the integration of the energy equation shows no variation with z, considering Equation (10), and introducing the curvature items to the equation, to describe cylindrical shells, the following equation yields:

$$\begin{gathered} U=\frac{D}{2}{{\displaystyle \iint }}_A \left\{{\left({\nabla }^{2}f\right)}^2-2\left(1-\nu \right)\left[\frac{{\partial }^{2}f}{\partial x^2}\frac{{\partial }^{2}f}{\partial y^2}-{\left(\frac{{\partial }^{2}f}{\partial x\partial y}\right)}^2\right]\right\}dxdy \\ =\frac{1}{2}{{\displaystyle \iint }}_A \left\{M_x\frac{{\partial }^{2}f}{\partial x^2}+2M_{xy}\frac{{\partial }^{2}f}{\partial x\partial y}+M_x\frac{{\partial }^{2}f}{\partial y^2}\right\}dxdy \\ =\frac{1}{2}{{\displaystyle \iint }}_A D\left\{{\left(K_x+K_y\right)}^2-2\left(1-\nu \right)\left[K_x{K}_y-K_{xy}{}^2\right]\right\}dxdy \end{gathered}$$

(20)

In Equation (20), A represents the shell surface area. Furthermore, by letting curvature components $K_x$,$K_y,$ and $K_{xy}$ be curvature variations represented by as $x_x$, $x_y,$ and $x_{xy}$, we could express the bending strain energy as:

$$U_b=\frac{D}{2}{{\displaystyle \iint }}_A\left\{{\left(x_x+x_y\right)}^2-2\left(1-\nu \right)\left[x_x{x}_y-x_{xy}^2\right]\right\}dxdy$$

(21)

The in-plane strain energy due to the membrane effect, the mid-surface stretching then becomes:

$$\begin{gathered} U_m=\frac{1}{2}{{\displaystyle \iint }}_A\left\{N_x{\epsilon }_{x0}+N_y{\epsilon }_{y0}+N_{xy}{\gamma }_{xy0}\right\}dxdy \\ =\frac{Et}{2\left(1-{\nu }^2\right)}{{\displaystyle \iint }}_A\left\{{\left({\epsilon }_{x\theta }+{\epsilon }_{y\theta }\right)}^2-2\left(1-\nu \right)\left[{\epsilon }_{x\theta }{\epsilon }_{y\theta }-\frac{1}{4}{\gamma }_{xy\theta }^2\right]\right\}dxdy \end{gathered}$$

(22)

In terms of the circumferential element displacement, as depicted in Figure 4, we used $ds$ to represent the original shell surface, $d{s}^{\prime }$ is the displaced shell mid-surface, and $df$ represents the displacement in z direction. In such terms, the distortion rate of the shell element is revealed as:

$$\frac{\partial \theta }{\partial s}=\frac{\partial \theta }{r\partial \theta }=\frac{1}{r}$$

(23)

$$\frac{\partial {\theta }^{\prime }}{\partial s^{\prime }}=\frac{d\theta +\left(\frac{{\partial }^{2}f}{\partial s^2}\right)ds}{\left(r-f\right)d\theta }$$

(24)

$$x_{\theta }=\frac{\partial {\theta }^{\prime }}{\partial s^{\prime }}-\frac{1}{r}\approx \frac{1}{r^2}\left(f+\frac{{\partial }^{2}f}{\partial {\theta }^2}\right)$$

(25)

$$x_{x\theta }rd\theta =-\left(\frac{{\partial }^{2}f}{\partial x\partial \theta }+\frac{\partial v}{\partial x}\right)d\theta$$

(26)

In Equation (26), $x_{x\theta }$ represents the distortion rate of the cross section of the shell’s surface, and $rd\theta$ is the arc length of the selected shell element.

Under the assumption that the shell mid-surface is inextensible, the shell surface displacement shows variations in the x-axis, along the radial direction and normal to radial direction. Ignoring the mid-surface strain change, the first simultaneous function of Equation (7) could thus be updated as:

$$\{\begin{array}{l}\frac{\partial u}{\partial x}=0={\epsilon }_x\\ \frac{\partial v}{r\partial \theta }-\frac{f}{r}=0\\ \frac{1}{r}\frac{\partial u}{\partial \theta }+\frac{\partial v}{\partial x}=0={\gamma }_{x\theta }\end{array}$$

(27)

All of the shell surface displacement functions could be expressed accurately with trigonometric functions embedded in relatively simple algorithms. We assumed the nodal displacement function of the tower shell approximate by a sum of two items as:

$$u=u_1+u_2, v=v_1+v_2, f=f_1+f_2$$

(28)

In that way, each displacement component is described as:

$$\{\begin{array}{l}{u}_1=0\\ u_2=-r{\displaystyle \sum _{n=1}^{\infty }}\frac{1}{n}\left(b_n\sin n\theta +d_n\cos \theta \right)\\ v_1=r{\displaystyle \sum _{n=1}^{\infty }}\left(a_n\cos n\theta -c_n\sin \theta \right)\\ v_2=r{\displaystyle \sum _{n=1}^{\infty }}\left(b_n\cos n\theta -d_n\sin \theta \right)\\ f_1=-r{\displaystyle \sum _{n=1}^{\infty }}n\left(a_n\sin n\theta +c_n\cos \theta \right)\\ f_2=-r{\displaystyle \sum _{n=1}^{\infty }}n\left(b_n\sin n\theta +d_n\cos \theta \right)\end{array}$$

(29)

The shell strain energy stays at a minimum under the hypothesis of inextensible shell performance. The bending strain energy $U_{b1}$, under the condition $\frac{{\partial }^{2}f}{\partial x^2}=x_x=0$, therefore follows:

$$U_{b1}=\frac{D}{2}{{\displaystyle \iint }}_A\left\{{x_{\theta }}^2+2\left(1-\nu \right)x_{x\theta }^2\right\}rd\theta dx$$

(30)

Equation (30) is distinguished from Equation (21). The former provides the integration in polar coordinates and is split from stretching strain energies. Moreover, by substituting Equation (27) into Equation (30), combined with displacement functions (28) enables quantifying the basic total strain energy of the unstretched shell section (the first element obliterated) as follows [21]:

$$U=\pi DL{\displaystyle \sum }_{n=2}^{\infty }\frac{{\left(n^2-1\right)}^2}{r^3}\left\{r^2{n}^2\left[a_n^2+{\left({\mathrm{c}}_n\right)}^2\right]+\frac{1}{3}{L}^2{n}^2\left[b_n^2+{\left({\mathrm{d}}_n\right)}^2\right]+2\left(1-\nu \right)r^2\left[b_n^2+{\left({\mathrm{d}}_n\right)}^2\right]\right\}$$

(31)

where, $a_n$, $b_n$, $c_n$, and ${\mathrm{d}}_n$ are constants that depend on loading conditions. Once defined, these lead to trigonometric functions that express initial assumed displacements. For example, when n = 1,

$$\{\begin{array}{l}{v}_1=r\left(a_{1}cos\theta -c_{1}sin\theta \right)\\ f_1=r\left(a_{1}sin\theta -a_{1}cos\theta \right)\end{array}$$

(32)

The longitudinal and circumferential buckling under pure axial compressive load, expressed by trigonometric functions, can be assumed as in Figure 5.

Figure 5 represents a typical surface displacement conformed by periodic functions with a half-wavenumber equal to 7; based on which then the circumferential deflection becomes,

$$f=-C_{0}sin\frac{n\pi x}{L}$$

(33)

where $C_0$ is a constant and n equals the half-wavenumber along the shell length L. The radical strain of the unit mid-surface circumference can thus be expressed as:

$$\begin{gathered} {\epsilon }_{xl}=-\frac{N}{Et} \\ {\epsilon }_{x}E=\frac{E}{1-{\nu }^2}\left({\epsilon }_{x_2}+\nu {\epsilon }_{{\theta }_2}\right) \\ {\epsilon }_{x_2}+\nu {\epsilon }_{{\theta }_2}=\left(1-{\nu }^2\right){\epsilon }_{x_l} \\ {\epsilon }_{{\theta }_2}=-\nu {\epsilon }_{x_l}+\frac{C_0}{r}sin\frac{n\pi x}{L} \\ {\gamma }_{xy}=x_y=x_{xy}=0 \\ {x}_x=C_0\frac{n^2{\pi }^2}{L^2}sin\frac{\pi nx}{L} \end{gathered}$$

(34)

In the above equations, ${\epsilon }_{x_2}$ and ${\epsilon }_{{\theta }_2}$ represent mid-surface axial strains after buckling. The relationship between virtual work and applied axial load N is therefore [21]:

$${\delta }_w=2\pi N\left[\nu {{\displaystyle \int }}_0^L{C}_{0}sin\frac{n\pi x}{L}dx+\frac{r}{2}{{\displaystyle \int }}_0^L{\left(C_0\frac{n^2{\pi }^2}{L^2}sin\frac{n\pi x}{L}\right)}^{2}dx\right]$$

(35)

To sum up the virtual works ${\delta }_{U_m}$ (membrane strain energy), and ${\delta }_{U_b}$ (bending strain energy), the equation above progresses to:

$$\begin{gathered} {\delta }_U={\delta }_{U_m}+{\delta }_{U_b} \\ =-2E\nu \pi t{\epsilon }_{x_l}{{\displaystyle \int }}_0^L{C}_{0}sin\frac{n\pi x}{L}dx+\frac{{\pi }^2{C}_0^{2}EtL}{2r}+C_0\frac{2{\pi }^{4}aLD\pi n^4}{2L^4} \end{gathered}$$

(36)

when ${\delta }_w={\delta }_U$.

Finally, by considering Equations (35) and (36) while letting $C_0\ne 0$, one identifies the following approximation to the critical load $N_{cr}$:

$$\left(\frac{\pi EtL}{2\gamma }+\frac{{\pi }^2\gamma n^{4}D}{2L^3}-\frac{{\pi }^{3}r{n}^{2}N}{2L}\right)C_0^2=0$$

(37)

$$N_{cr}=D\left(\frac{n^2{\pi }^2}{L^2}+\frac{E{L}^{2}t}{D{r}^2{n}^2{\pi }^2}\right)$$

(38)

3.2. Bending Moment and the Ovalisation

The study of transitional shell geometry into oval-shaped due to bending is essential for the analysis of a wind turbine tower under combined loading conditions. Oval-shaping occurs as cylindrical shell segments flatten to planes of bending as in Figure 6a and Figure 7b [39,40].

Along the longitudinal direction of tubular shells, pure bending wrinkles shell regions subject to the highest compression stress, as shown in Figure 6. The pure bending moment thus triggers prebuckling on the compressed wall while stretching the tension side, as in Figure 6a. The interaction between bending and axial stress is complex and has deserved ample study in the past [37,41,42,43,44,45,46]. The present investigation thus intended to shed light on the nature of the collapse mechanism derived from lateral loading affecting wind turbine towers, through the mathematical formulation exposed in the following paragraphs.

Brazier and Southwell [47] defined the non-linear moment-curvature of a section as:

$$M=EI\left(M\right)\frac{d^{2}u}{d{x}^2}$$

(39)

The terms $\frac{d^{2}u}{d{x}^2}$ representing the curvature of the shell section, from which the limit of the moment derived as [47]:

$$M_{lim}=2\sqrt{2}\frac{ER{t}^2\pi }{9}$$

(40)

In line with the previous section, the bifurcation moment for perfectly cylindrical shells turns out to be [48]:

$$M_{Bifurcation}=Er\pi \frac{t^2}{\sqrt{3(1-{\nu }^2)}}$$

(41)

As the cross section of the shell progresses into an oval-shape the buckling stress ${\sigma }_c$ with respect to the local minor axis $R_1$ (after Hutchinson [49]) turns out to be:

$${\sigma }_c=\frac{1}{\sqrt{3(1-v^2)}}\left(\frac{Et}{R_1}\right)$$

(42)

Figure 6c identifies an equivalent radius shell $R_1$. This, together with Figure 4 and Equation (33), enables the formulation of the local buckling equation linked to the minimum curvature strip. Furthermore, the moment-balance stress configuration depicted in Figure 7a enables devising the following extra control functions:

$$\begin{array}{l}N=F_{C_1}+F_{C_2}-F_t=F_y\left(Area{C}_A-Area{T}_A\right)\\ \hspace{1em}=F_y\left(2\left(\pi -\alpha \right)rt-2\alpha rt\right)\\ \hspace{1em}=F_y\left[2rt\left(\pi -2\alpha \right)\right]\end{array}$$

(43)

$$\begin{array}{l}M=F_t\left(r_t+r_{C_1}\right)-F_{C_1}{r}_{C_1}+F_{C_2}{r}_{C_2}\\ \hspace{1em}=F_y\left\{2\alpha rt\left(r_t+r_{C_1}\right)-\left[2\left(\frac{\pi }{2}-\alpha \right)r_t{r}_{C_1}\right]+\left(\pi rt\cdot r_{C_2}\right)\right\}\end{array}$$

(44)

As the ovalisation shown in Figure 7b manifests, the control equations become:

$$\begin{array}{l}N=F_{C_1}+F_{C_2}-F_t=F_y\left(Area{C}_A-Area{T}_A\right)\\ \hspace{1em}=F_y\left[2\pi R_2\left(\pi -2\alpha \right)t+2\alpha R_{1}t-2\alpha R_{1}t\right]\\ \hspace{1em}=F_y\left[2\pi R_{2}t\left(\pi -2\alpha \right)\right]\end{array}$$

(45)

$$\begin{array}{l}M=F_t\left(r_t+r_{C_1}\right)-F_{C_1}{r}_{C_1}+F_{C_2}{r}_{C_2}\\ \hspace{1em}=F_y\left\{2\alpha R_{1}t\left(r_t+r_{C_1}\right)-\left[R_2\left(\pi -\alpha \right)t\left(r_{C_1}\right)\right]+\left(\pi rt\cdot r_{C_2}\right)\right\}\end{array}$$

(46)

In summary, based on energy methods that incorporate linear elastic effects evolving into section ovalisation, induced by bending-axial interactions, we formulated a prebuckling system of linear equations to further study transient and post-buckling stages, via the numerical simulation discussed in the following paragraphs.

4. Transient and Post-Buckling Analysis

In the current approximation, the post- and prebuckling analysis are separate by a transient stage. Although static and quasi-static methods do not fully capture the dynamic nature of buckling, they are valid for simulating plate stresses under particular conditions. Linear analysis results are well founded with eigenvalue analysis that defines singular points along the shell deformation curve, as well as to identify bifurcation or snap through points.

The geometrical nonlinear static solutions of tubular shells typically contain a progression of bifurcation or snap-through points that derive from finite element static step-by-step analyses. The arc length method being one of the most efficiency methods for this purpose, particularly when merged with the modified Riks method [50,51,52] and embedded into a finite element algorithm like the one powered by ABAQUS (ABAQUS Inc., Palo Alto, CA, USA) [53].

Is worth noting that the arc length path in ABAQUS is insensitive to the stability of the structure. The short board of the modified Riks method in ABAQUS is unable to solve the bifurcation buckling behaviour of a perfect shell, because the increment step should be continuous as opposed to discreet throughout the post-buckling analysis. However, the initial imperfections can still feed into this process to find meaningful buckling modes, as external load increments approach the ideal critical load.

5. Finite Element Analysis

5.1. Modelling Details and Assumptions

A series of finite element models were developed in this study considering four different shell lengths and three different tube thicknesses, in all cases addressing a tower of 6 m diameter, fixed at the base. Tube segments were of 9 m, 12 m, 15 m, and 20 m length with 24 mm, 27 mm, and 36 mm thickness, respectively. This enables parameterising wind turbine towers with length–diameter (L/t) and diameter–thickness (D/t) ratios. Typically, D/t varied from 166.67 to 250 for R/t ratios in the range of 83.34–125, as shown in

Table 1. Model geometry size in the FEM analysis.

No.	Diameter (D)	Length (L)	Thickness (t)	L/D	D/t
1	6000	20,000	24	3.33	250
2	6000	15,000	24	2.5	250
3	6000	12,000	24	2	250
4	6000	9000	24	1.5	250
5	6000	20,000	27	3.33	222.22
6	6000	15,000	27	2.5	222.22
7	6000	12,000	27	2	222.22
8	6000	9000	27	1.5	222.22
9	6000	20,000	36	3.33	166.67
10	6000	15,000	36	2.5	166.67
11	6000	12,000	36	2	166.67
12	6000	9000	36	1.5	166.67

. L/D ratios, on the other hand, varied between 1.5 and 3.33, otherwise expressed as L/R falling in the range of 3–6.66.

Our study considered three loading conditions: axial compression, pure bending, and axial–bending interactions. The above parameterisation derived in 48 tower shell prototypes, amongst which 12 relate to linear buckling models and 36 specific loading configuration models, for the purpose of the respective eigenvalue analysis. The boundary conditions used for axial compression analyses consists of a fixed bottom edge and free to displace and rotate the top edge. For pure bending we restricted any horizontal displacement of the bottom edge and rotations with respect to the x and y axis, again leaving the top edge free to displace and rotate. The combined loading scenario replicates the boundary conditions set for pure bending. Keeping with the modelling, the shapes of the tube ends did not deform after assigning tie constraints along the circumferences and according to a reference point around the geometrical centre of the shell, which simulated the connection flange and ring stiffening inside the wind turbine tower; this in turn prevented the ovalisation taking place at both ends. The Finite Element Method software ABAQUS was therefore calibrated for this study by choosing S4R element types [53] whose four nodes reduce the integration of the shell element linearly by means of the Riks method [50] including linear buckling eigenvalue processing and iterations for numerical stability. The maximum number of increments of static Riks analysis was set as 1000, here considered adequate to approach the collapsed state of the tower shell, and fixes the cease for increments, which results in being particularly efficient when running models containing various structural configurations. The minimum arc length increment was 10⁻⁸ to accurately define the performance in the neighbourhood of bifurcation or snap-through points. The number of the thickness integration point was 7 and the element size of all models was approximately 200 mm × 200 mm with the curvature control at a maximum of 0.1, which provided additional post-buckling accuracy. The fixed diameter guarantees insensitivity to circumferential wavelength according to element numbers while the longitudinal element performance shows less inference on the simulation results. The elastoplastic parameterisation of material and stress–strain information map experimental data linked to S355 steel. This data constructed an 11-point curve characterised by a Poisson ratio of 0.3.

The geometric imperfections introduced in our modelling were reflected in modal shape imperfections calculated with a linear buckling analysis for axial compression scenarios. As shown in Figure 8, the first compression modal was approximately symmetric with respect to the z axis while the second and third modal were asymmetric. The initial loading condition of the wind turbine tower shell was assumed purely compressive hence without eccentricity generated from lateral motion of the upper segment of the structure, including nacelle and blades, during static load. The first three compressive modal shapes then reproduced complex initial imperfections at rates of 10% for the 1st modal, 5% for the 2nd modal, and 2.5% for the 3rd modal. This enabled us to determine the lower buckling bound effectively, which depicts a typical compression modal combination, as represented sin Figure 8.

The amount of the axial compressive load applied on the upper reference point during the Riks algorithm fed from the previous linear buckle analysis for each model, and the amount of bending moment was assigned at relatively large scales to exceed the plastic range of material stress as predicted by Equations (44) and (46). The Riks increments implemented in ABAQUS would thus approach the failure state until the fixed 1000 steps are completed otherwise the job is aborted to prevent a lack of convergence.

The knock down factor $\alpha$ for the imperfect shell buckling is given by [54] as:

$${\sigma }_1=\alpha {\sigma }_p$$

(47)

where ${\sigma }_p$ represents the critical compressive stress component of the perfect shell in the one direction. The experimental data concluded by NASA [55,56] is shown for a better reference in Figure 9.

In line with the above, the recommended uniaxial comprehension expression is given as [56]:

$$\alpha =1-0.901\left(1-e^{-\sqrt{\left(R/t\right)}/16}\right)$$

(48)

The buckling simulation models in this study founded in verifications of axial compression according to Equation (48) while the modelling results marked in Figure 10, in which the ${\sigma }_1$ was obtained by imperfect shell models and ${\sigma }_p$ derived from perfect shell models:

The three specific R/t ratios of wind turbine towers considered in our numerical simulations were 83.34, 111.11, and 125, see Table 1. This proved the basic settings of the buckling analysis fit experimental evidence of imperfect shell performance for similar R/t ratios. The net deviations obtained seemed acceptable with our simulation results lying above the experimental regression curve shown in Figure 9, as expected.

In terms of bending load scenarios, the simulation results validated with data provided by Dimopoulos [14], who investigated the bending performance of the wind turbine tower both experimentally and numerically. Figure 11 compares the corresponding load-displacement curves of clamped cylindrical shells using data from the referred sources. The R/t ratio in [14] was set to 50 while the L/D ratio to 3, therefore our model No. 9 whose R/t and L/D ratios were respectively 83.34 and 3.33 enabled validation, given the proximity of valued parameters with respect to experimental evidence. The concentrated load applied at the centre of the top edge in our model No.9 was equivalent to 34.91 kN, by considering the ratio of sectional resistance against acting moments. The progression of curves fitted logically while the normalised deviation would correspond to differences of different material quality and geometrical configurations.

5.2. The Typical Buckling Induced Failure under Combined Loading Conditions

For combined loading we selected for validation our model of 20 m long with 24 mm wall thickness (R/t = 125, L/R = 6.67). Figure 12 illustrates the structural performance that we observed. The sequence shown depicts a characteristic progression of failure. That simulation resulted from progressing increment steps of combined loading, which revealed the longitudinal half-wavenumber, emerging between steps 12 and 18, which aligns with the prebuckling displacement that is implicit in Equation (33), characterised by features proper of periodic trigonometry functions. At increment step 16, the flattening effect resulted in mid-length stress concentrations and the original longitudinal wavelength kept unchanged until incremental step 25. Then, the inward single buckle dent developed instantly within the time framed by steps 20–25 while buckling started expanding more actively from increment 38. Furthermore, the circumferential half-wave number developed between steps 25 and 135, gradually evolving from one single buckle point to eight inward dents (two other dents were on the side viewpoint at step 135). The number of buckle point development stopped in the post-buckling stage when the shell section underwent large displacements, which triggered a rapid decrease of the load factor while the displacement and rotation at the top simultaneously increased. Past this point, the wind turbine tower could not provide regular service requirements. The compressive side circumferential half-wave number in the range 6–8 is commonly observed in real historical collapse cases reported by onsite photos [6,34,35,57]. The lateral side view shown for step 200 shows the shell surface flattening on the tension side and expansion of buckling across the cylindrical shell. These modelling thus enabled us to accurately devise serviceability and strength limit states [1,58,59].

5.3. Buckling Progression under Different Loading Conditions

Figure 13 shows the deformed tower tube (20 m at 24 mm) when subject to uniform compression; Figure 14 shows the equivalent surface deformation paths under pure bending while Figure 15 reproduces these for combined effects, for the critical compression region. Figure 13 shows a longitudinal wave near the mid-length of the tubular shell, which derives from initial imperfections (eigenvalue results determined the highest gradient in the mid-length areas). It can be seen in that figure how the section curvature increased under pure bending with the loading incremental. At step 39 we counted 4 half waves while the number decreased to 3 at step 121, 1 at step 168, and back to 3 at step 187. The system became unstable past step 232. Classical buckling theory describes these changes as dissipation and rearrangement of strain energy on the shell surface, both through stretching ${\lambda }_m$ and bending ${\lambda }_b$. Under uniform compression, buckling also appeared at the fixed edge of the bottom tube segment, partially stimulated by the imposed motion constraints. Figure 14 shows asymmetric bending taking place at high incremental steps also as a result of initial imperfections, which revealed the failure mode was highly sensitive to such imperfections, as confirmed by experimental research [17,44,56,60,61]. This was captured during our interpretation of buckling equilibrium paths discussed above.

The results obtained for pure bending load are in line with previous modelling outcomes [40,62,63], which predict that initial waving develops at the middle span, followed by the flattening of both compressive side and tension sides, to induce the cross-sectional ovalisation illustrated in Figure 6 and Figure 7. The deformation patterns observed at step 39 and step 71 justified the use of a simple trigonometry function in Equation (33) to simulate longitudinal periodic displacement. Figure 14 also suggests that the curvature at mid-span unchained the buckling collapse framed between steps 71 and 153, as expected, the tension side remained flat during that process.

The performance of the shell under combined loading naturally reproduces compression and bending effects. Longitudinal wave shapes resemble those seen for axial compression, i.e., containing a long extension of prebuckling surface displacement due to initial imperfections and the joint loading. The bending moment partially flattened the compression side to then increase half-waves along the length of the segment [21]. Buckling emerged at step 39, with a large dent and two adjacent bulges. Buckling, along the longitudinal direction thus developed rapidly as the shell membrane effect vanished. In the corresponding numerical process the factor of ${\lambda }_m$ was eliminated to define the lower bound of bearing capacity, which fuelled the shell’s progression to the post-buckling branch. In summary, the combined loading condition accelerated the evolution of buckling via mid-length stress concentration initiated by pure bending and magnified by compressive effects, as shown in Figure 12.

5.4. Buckling Development for Different Shell Lengths

The previous results were enhanced with a series of shell models. These included a segment of 9 m length with 24 mm thickness (L/R = 3), which enabled comparing shell length effects (L/R ratio)—note the above results related to a 24 m-long one with L/R = 6.67, hence the R/t ratio remained equal to 125. Figure 16 shows the progression of stress and deformation through 10 characteristic configurations. It is noticeable that the buckling affected area differed for the long and short tube. The failure of the short tube was located near the bottom as opposed to the mid-length as observed in the long segment. The distribution of the first visible waves varied less in the sort segment, as predicted by Equation (33) where the assumed displacement function has a denominator of length L. The stress concentration observed at step 6 developed near the base, which as seen before relates to constraints imposed at the boundary (these represent the flange to flange connection of adjacent segments or stiffening rings that induce sudden stiffness changes). Three half waves appeared at step 16, which caused stress redistribution, local flattening, and curving hence inducing geometrical changes of adjacent shell elements.

Figure 17 shows stress and deformation associated to axial compression and pure bending. This figure illustrates how the axial compression led to symmetric buckling, both at the top and bottom of the shell with buckle bulges developing 800 mm from the edge, which we understood as a characteristic failure mode of cylindrical metal shells [17,45]. In practice, we could expect such features to appear on short water tanks and silos, even though these could have different R/t ratios with respect to the shell section studied here [21,64,65,66]. Figure 17b describes pure bending effects. There we could see how stress concentration took place at step 31, pinpointed at both edges; this followed initial periodic waves and flattening combined through the tube length. It thus seemed that the shorter shells are more susceptible to boundary conditions as opposed to the assumed initial displacement or geometrical imperfections.

Figure 18 and Figure 19 provide an energy perspective. These enabled tracking the dissipation of strain energy and the change of curvature change of two different shell lengths. In these figures, the vertical axis measures the strain energy magnitude in relation to the normalised longitudinal dimension and through a sequence of incremental steps.

In the early stage of the load incremental process, both the 9 m and 24 m shells were characterised by evenly distributed strain energy—with tiny periodic patterns along the length as the result of initial wavenumbers in the prebuckling stage. At that point, the strain-stress relationship was linear and according to Hook’s Law. As the load incremented, two different strain energy redistribution patterns revealed on each shell: the 9 m shell experienced a strain energy decline in the upper and mid-length and the accumulation of energy at the bottom end. In contrast, the 24 m shell underwent an increase of strain energy at its middle length followed by a large decrease at both edges.

Figure 20 addresses the evolution of strain energy distribution associated to pure bending for the 9 m long shell. In this case, bending, i.e., ovalisation, has an important contribution to buckling, seemingly defining strain and stress distribution and expansion, as shown in Figure 21.

Furthermore, Figure 21 shows the curvature of the 20 m shell under pure bending, which precedes local buckling. In this case, the curvature fluctuated throughout the length of the tube while becoming notorious at the middle span, which resembles the performance of a solid beam [21,27,67]. Figure 22 enhances the curvature appreciated at the lower scale in Figure 21. Notoriously, curvatures at steps 4, 8, and 12 were suitably described by periodic functions, with such configuration changing rapidly afterwards as the section underwent buckling—see also Figure 14. The longitudinal section curvature of the combined loading scenario was similar to the one observed for pure bending ones. This can be appreciated by scrutinizing Figure 21, Figure 22 and Figure 23.

On the other hand, the circumferential strain energy and evolution is shown in Figure 24. The transition from smooth to highly asymmetric curves seen along the length of the shell replicated along the circumference, however the distribution of strain varied radically between the short and long tubes. In general, the strain distribution seems step-like for the short tube segment (Figure 24a) while for the long segment, these appeared to concentrate in specific points (peaks).

5.5. Buckling Development for Different Shell Thickness

Further to the scrutiny of the evolution of stress and deformations for varying shell length, we present here the results collected for varying shell thickness. The models to compare include the 9 m-long segment with a thickness of 36 mm (R/t = 83.33) and 24 mm thickness (R/t = 125). These correspond to a L/R ratio of 3. Figure 25 shows the deformation and surface stress distributions via 10 characteristic configurations in relation to combined loading. The features shown in Figure 25 contrast with those shown in Figure 16, for the 24 mm shell. The 36 mm thick shell shows less dependency on boundary conditions, arguably due to the lower rate of variation of stiffness from the inner sections with respect to the edges. Stress augmentation became visible at steps 20–24, as a half-ring bulging emerged, leading to further circumferential wavenumber. The circumferential wavenumber passed from 1 to 5 then to 9, which contrasted with the more gradual evolution seen in the 24 mm thick shell.

According to these results, the lower R/t ratio provided better energy expansion for shell structures in buckling, which manifested at early load stages. The higher the amount of wavenumber during the prebuckling state, the better the energy dissipation, and the more the similarity to pure bending buckling (buckling was triggered at the middle length of the member). Based on Equations (31), (36), and (46), the shell thickness was proved as a controlling factor for higher energy storage capacity, which in turn allows multiple buckling equilibrium paths, as a basis the developed mathematical framework [68]. This mechanism seen in this investigation is consistent with those seen elsewhere [10,21,29,30,44,68,69].

6. Conclusions

In the present paper, the types of the general buckling behaviour of steel tubular wind turbine towers in different parameters of L/R, R/t, and under different loading conditions were investigated and the following conclusions were drawn:

• The ovalisation effect underpinned buckling derived from bending and axial interactions. This implies flattening the compression segment, longitudinally, and tapering the local curvature, circumferentially.
• The buckling point tended to appear at the mid-length of the shell section although the post-buckling scenarios were less redundant than axial compression ones, which reaffirms previous statements, as in [68,69]. It is worthy to note that multiple models rarely predicted bifurcation branching.
• Under the combined load, a lower L/R ratio showed strong dependency on boundary conditions, whereas higher L/R show patterns that resemble those observed under pure bending buckling. This could be due to a larger shell surface area enabling redistribution of energy that fuels buckling mechanisms.
• Similarly, lower R/t boosted energy dissipation and expansion, which in turn enhanced surface deformation rates hence buckling equilibrium paths.
• Strain energy appeared strongly dependent on local curvature, a property that feeds into linear functions that describe local bending and membrane strain. Furthermore, the arc length method fit any buckling analysis of shell structures, as it determines the realistic solution for instability problems, particularly when those derive from initial geometrical imperfections.
• Combined axial compression and bending reduced the critical buckling load of a cylindrical steel shell. This appeared in Figure 26 and Figure 27 shown below. On the other hand, a load factor corresponding to a 24 mm thickness shell peaked for the 15 m-long model, other than the 9 m or 20 m, which revealed an optimum solution for that load factor interval; noting that, the maximum load factor was 0.312, far less than the experimental recommended knock-down factor of 0.68, for the ideal theoretical buckling load ${\lambda }_{cr}$ [55,56].