Approximation of Column Buckling

1. Introduction

This discussion will explore then dynamics of equilibrium states of a column or beam as it buckles under compressive loading or phenomena which alter the stability of a structure. For instance, during heat waves the buckling of rails has been a cause of train derailments.

Buckling due to glaze icing Source: University of Cape Town http://www.mth.uct.ac.za/~myers/	Buckling of rails Source: Brock University http://www.physics.brocku.ca/courses/1p23/Heat/rail.html

Consider the system d²y/dt² = - w ²y, where w is a constant.

General solution: d²y/dt² = - w ²y    Þ   y = C₁ sin w t + C₂ cos w t .

We know there is a family of oscillating solutions consisting of the sinusoidal functions, with each solution having some amplitude and phase, and each following elliptical orbits in phase space. In the solution of this second-order linear differential equation, the sine and cosine functions are eigenfunctions of the linear system. Keeping this system in mind, the attempt will be to generalize it to nonlinear equations.

2. Equilibrium States

Along the axis of a column (or a wire), force can be applied in either direction. If the force is directed outward, tension is created, with periodic orbits in time. For example, a string under tension can produce acoustic vibrations. This phenomenon has been studied and found to include nonlinear vibrations, such as those found in guitar strings (Tufillaro, 1991). In such a system under tension, there is only equilibrium state once kinetic energy is lost.

In contrast, if the force is compressive, it can be shown that natural vibrations of a beam will decrease in frequency as load is applied. The frequency drops to zero along a square root function, and provides a way to solve for the critical buckling load (James et al, p.104). As more weight is added and buckling occurs, new equilibrium states emerge as the load is increased.

Instead of using functions in the time domain, y(t), equilibrium solutions will live in the spatial domain, y(s). Solutions will be presented as a bifurcation at a critical load, using the displacement as a metric of some loading force, P, applied to a vertical column.

These are sketches of equilibrium states of columns. Each column supports a mass, m, and may buckle in different modes, k. A series of pitchfork bifurcations will be shown to exist as the mass is increased. When stability changes and a column abruptly moves from one mode to another, this is called mode jumping. For analysis purposes, a column is taken as length p, and graphs will be shown with columns in a horizontal position (along the x-axis), with the y-axis being the displacement.

Placing one endpoint at the origin and allowing the other to slide along the x-axis, it is clear that the upper limit of x is not a constant. Rather than sticking with Cartesian coordinates, the function can be specified as y(s), where s is the position along the arc length. This function y(s) is the displacement of the column for any s between 0 and p, and this independent variable provides initial conditions y(0) = y(p) = 0 .

3. The Euler-Bernoulli Law

Analogous to Hooke's law for springs, the displacement function can be related to the applied force. By imagining that the column contains small springs within its structure, small amounts of torque are applied to each segment ds by the force P with a lever arm being equal to the displacement y(s). The Euler-Bernoulli law relates this torque to the local curvature of the function at s, which is f'(s), where f is angle of the tangent line. The use of f' is natural, since it relates to the severity of the column buckling, and (unlike y'') it preserves the isometric character of the local forces acting to straighten the column.

The angle f relates to the slope as f = arctan(dy/dx).
The two torques will balance and there is equilibrium when
P y(s) + M f'(s) = 0, where M is the modulus of flexural rigidity (a measure of the stiffness of the column).

4. Changing to the second order

We can apply the following to get the problem down to a single normalized parameter with a second-order differential equation.

By the Euler-Bernoulli law, P y(s) + M f'(s) = 0. This is equivalent to -f'(s) = (P/M) y(s).
With l = P/M, the equation is parameterized as -f'(s) = l y(s). Taking the derivative of both sides, -f''(s) = l y'(s)

Looking at the right triangle, y'(s) = dy/ds = sin f(s). So by substitution, -f''(s) = l sin f(s).
This is a second-order differential equation involving only the variables f and s, along with some parameter l.

Recalling that y(0) = y(p) = 0, and that -f'(s) = l y(s), we have initial conditions f'(0) = f'(p) = 0.

Our goal is to solve: -f''(s) = l sin f(s) subject to f'(0) = f'(p) = 0.

5. What happens if we use sin f » f ?

We can attempt to linearize the second-order problem: -f''(s) = l sin f(s), subject to f'(0) = f'(p) = 0.

With the approximation sin f » f, we now have the problem f''(s) = -l f(s).
But this is just like the spring problem, with solutions f(s) = C₁ sin l^1/2 s + C₂ cos l^1/2 s .
So the angle of the column is varying sinusoidally along the arc length.

Taking the derivative, f'(s) = C₁ l^1/2 cos( l^1/2 s ) - C₂ l^1/2 sin( l^1/2 s ).
Along with our initial conditions, this tell us C₁ = 0, and that either C₂ = 0 or sin( l^1/2 p ) = 0.
When C₂ is nonzero, l^1/2 must be an integer, so all solutions have the form f(s) = 0 or f(s) = f₀ cos k s , where k is any integer and l = k².

Since -f'(s) = l y(s), the solution is y(s) = (-1/l) (-f₀ k sin k s) = (f₀ / k) sin k s which achieves a maximum magnitude of Y = f₀ / k . Because f₀ is arbitrary, possible bifurcations can be graphed as a sequence of vertical lines.

The function cos k s is an eigenfunction, and the perfect squares are eigenvalues. Buckling occurs only when l is a perfect square. In this solution, the angle f₀ can be any arbitrary real number. For example, nothing prevents the column from collapsing completely at l = 1. But paradoxically, it supports much heavier loads, say at l = 1,000,000,000.

6. What went wrong?

The higher-order terms are critical here. We can get a better approximation by adding terms from the Taylor expansion, such as sin f » f - f ³ / 3! . This leaves us with f'' = -l ( f - f ³ / 3! ) = -(l - lf ² / 6 ) f, which involves a cubic term that becomes difficult to solve analytically. In general, this type of problem where we are given u''(x) as a some function of u(x), along with some a pair of initial conditions for u'(x), is called a nonlinear Sturm-Liouville problem.

Rather than attempting a rigorous analytic solution, we can take a closer look at the troublesome lf ² / 6 term, which was basically taken to be zero in the linearization, and let this term be some parameter e in the interval 0 £ e < l . We are left with the linear system f'' = -( l - e ) f under the same initial conditions. Using theorems from nonlinear bifurcation theory, it has been shown that (under certain constraints) linear approximations of this type can produce qualitatively similar buckling solutions (Brown, 1993).

7. Improvements to the bifurcation diagram

f(s) = f₀ cos k s , where k is any integer and l - e = k²
Under this constraint, any value of l = k² + e can cause the column to buckle.
However, the steepest angle f₀ can no longer take on any real value, but is limited by the parameters.
If the smaller parameter e is nonzero, then e = lf ² / 6 > k²f ² / 6.
This tells us also that |f| £ (6e) ^1/2 / k . So there are pairs of approximate solutions for f₀ = ± (6l - 6k²) ^1/2 / k for every positive perfect square k² £ l .

This solution also tells us, quite intuitively, that the column will buckle more severely as the load is increased. But this solution allows f₀ to exceed p/2, displaying an inaccurate graph as e becomes significant. This happens when:
    (6l - 6k²) ^1/2 > kp/2
Þ 6l - 6k² > k² p²/4
Þ 6l > k² (6+p²/4)
Þ l > k² (1+p²/24) » 1.4 k².
The accuracy of graphs will tend to be limited by this condition.
Looking at y(s) = (-1/l) (-f₀ k sin ks) = ± (1/l) (6l - 6k²) ^1/2 sin ks ,
the maximum magnitude is Y = (1/l) (6l - 6k²) ^1/2.

8. Animating different modes in Maple

For graphing purposes, a Cartesian function C(x) can be estimated using an arc length approximation,
sin ks = sin kx(s/x) » sin(kx sec f₀) , and graphs will be plotted from 0 to p.
(An even better solution would integrate ds = (ds/dx) dx along the arc length)
Animations are shown with modes k=1 and k=2, and use changing values of l = L :

f0 = ± (6l - 6k2) 1/2 / k	y(s) = (k/l) f0 sin ks	C(x) = (k/l) f0 sin(kx sec f0)
in red 0 < s < p	in green 0 < s < p	in yellow 0 < x < p cosf0

with(plots):
k := 1;
Phi0 := sqrt(6*L-6*k^2)/k;
animate(plot,[[Phi0*cos(k*s),Phi0*k/L*sin(k*s),Phi0*k/L*sin(k*s/cos(Phi0)) *`if`(s/3.14 < cos(Phi0),1,0)],s=0..Pi],L=1.05*k^2..0.95*k^2);

Mode 1 Animation	Mode 2 Animation

References

[1] S. Strogatz, Nonlinear Dynamics and Chaos, Perseus Books Publishing, 1994.

[2] R. Brown, A Topological Introduction to Nonlinear Analysis, Birkhäuser Boston, 1993.

[3] N. Tufillaro, T. Abbott, J. Reilly, An Experimental Approach to Nonlinear Dynamics and Chaos, Addison-Wesley Publishing Company, 1991.

[4] M. James, G. Smith, J. Wolford, P. Whaley, Vibration of Mechanical and Structural Systems (Second Edition), HarperCollins College Publishers, 1994.