Buckling of beam structures

01

Energy formulation

Although the displacements are assumed to be small, the second-order terms in the strain–displacement relations are retained because they are essential for the stability analysis. These terms generate the geometric stiffness matrix. [1]

Buckling is a bifurcation phenomenon. The structure follows the small-displacement equilibrium path up to a critical load, at which this equilibrium state loses stability and a new equilibrium path associated with lateral deformation becomes possible. Along the post-buckling path, large displacements may subsequently develop.

The strain energy of a beam element is:

The following developments retain all the terms required to obtain both the elastic stiffness matrix and the geometric stiffness matrix.

02

Euler–Bernoulli beam model

For the Euler–Bernoulli beam model, the constitutive matrix and the generalized strain vector are:

When the second-order contribution produced by the transverse displacement is retained, the axial strain contains the additional term (see Chapter 6):

After substitution in the strain energy and neglecting terms of order higher than two, the expression can be separated into an elastic part and a part that depends on the axial force:

The generalized strains and the slope of the beam axis are written as:

With $\xi=\bar x/L$, the displacement interpolation matrix is:

The element strain energy becomes:

For the two-node Euler–Bernoulli element, integration gives:

03

Timoshenko beam model — isoparametric linear element

The Timoshenko beam model includes shear deformation. The axial strain, shear strain and curvature are:

The linear isoparametric shape functions are:

The slope of the beam axis is:

Since $ds/d\bar x=2/L$, the element matrices are:

With one-point Gauss integration, the elastic stiffness matrix used in the program is:

04

Timoshenko beam model — third-degree shape function

A third-degree interpolation of the transverse displacement improves the behavior of the Timoshenko element. The axial displacement remains linear, while the transverse displacement and section rotation are interpolated consistently.

The slope of the beam axis, $\varphi=d\bar v/d\bar x$ (see Section 5.2), is obtained by differentiating this interpolation.

The generalized strain matrix keeps the form:

Let's check the above expressions. If $GA_0\rightarrow\infty$, the Timoshenko beam model becomes the Euler–Bernoulli one. The middle row of $\boldsymbol{B}$ becomes zero and disappears, as in Chapter 6:

The slope matrix also becomes the third row of the interpolation matrix given in Chapter 6:

The elastic and geometric matrices are again obtained from:

The elastic stiffness matrix is:

For the geometric stiffness matrix, it results:

The local matrices are transformed to the global coordinate frame:

05

Buckling eigenvalue problem

The element matrices are assembled into the global elastic and geometric stiffness matrices. After applying the boundary conditions, the buckling problem is:

or, equivalently:

With the sign convention used here, the axial force $N$ is negative in compression. The geometric stiffness matrix of a compressed member is therefore assembled using $-N$.

The eigenvalue $\lambda$ is a load multiplier applied to the reference loading used in the preliminary linear analysis. Since a unit compressive load is used in the present examples, the numerical value of the smallest positive eigenvalue is equal to the critical buckling load in newtons. For a structure subjected to several loads, the same eigenvalue multiplies the complete prescribed reference loading pattern.

06

MATLAB programs

The main program has a very simple structure:

% main
clear
gen
stiff
S=K\F;
stress
geomK

The gen subprogram generates the input data. The stiff subprogram generates the global stiffness matrix and applies the boundary conditions. The linear system is then solved. The stress subprogram computes the axial forces in the elements, and geomK generates the geometric stiffness matrix and solves the eigenvalue problem.

The matrices used here follow the developments of Section 5.1, Section 5.2, Section 5.3, and Chapter 6.

07

Example 1

A 1 m long beam with a square $20\times20\ \mathrm{mm}$ cross-section is considered. Young's modulus is $E=200000\ \mathrm{MPa}$, and Poisson's ratio is $\nu=0.25$. A reference compressive load of $1\ \mathrm{N}$ is applied. The left end is clamped, while the right end is simply supported in the transverse direction.

The Euler critical load is:

The following table compares the first buckling load obtained with the three beam elements.

Timoshenko-based finite elements lead to a slightly smaller critical load because they also account for shear deformation. For slender beams, this effect becomes negligible and the Timoshenko solution approaches the Euler–Bernoulli result. The convergence of the two-node isoparametric element is slower than that of the third-degree element.

08

Example 2

The same beam is subjected to an axial compressive force and a small transverse force $H$.

This is not a buckling eigenvalue problem. It is a large-displacement analysis solved using the Updated Lagrangian formulation.

The purpose of this example is to verify the linear buckling prediction through a geometrically nonlinear equilibrium path. As the lateral load $H$ tends to zero, the equilibrium curves approach the critical load obtained from the eigenvalue buckling analysis.

09

Example 3

The structure consists of two beams with square $20\times20\ \mathrm{mm}$ cross-sections. Young's modulus is $E=200000\ \mathrm{MPa}$, and Poisson's ratio is $\nu=0.25$. The geometry, supports and loading are shown below.

10

Reference

[1] Timoshenko, S.P., Gere, J.M., Theory of Elastic Stability, McGraw-Hill, 1985.