Critical Load of Elastic Columns by the Energy Method

When the static method (§11-3) meets a column with variable cross-section or a distributed axial load, the deflection ODE becomes a variable-coefficient equation — usually without an analytical solution — or the stability equation's order grows too high to expand. In such cases, the energy method is highly effective.

As a column bifurcates from straight to bent equilibrium, the axial force stays constant. The bending strain energy is

Substituting $EIy'' = -M$ gives the form in terms of $y''$:

Load potential U_P

The external load $F_{\mathrm P}$ does positive work on the displacement $\Delta$; the potential is therefore negative:

With multiple concentrated loads $F_{\mathrm Pi}$ each moving $\Delta_{i}$:

Total potential · Ritz foundation

Combining $E_{\mathrm P} = U + U_{\mathrm P}$:

The stationary condition $\partial E_{\mathrm P}/\partial a_i = 0$ gives homogeneous algebraic equations in the generalized coordinates $a_i$.

Assume the buckling displacement can be written as a linear combination of basis functions with generalized coordinates:

Now $y(x)$ is fully determined by the $n$ coordinates $a_i$, reducing the infinite-DOF problem to the n-DOF stability problem of §11-2.

Guidelines for choosing the trial function

• Must satisfy the displacement boundary conditions (geometric constraints); satisfying force boundary conditions also helps but is not required;
• Should be close to the real buckled shape — the closer, the better;
• Adding more terms should improve the estimate; if $n+1$ terms give essentially the same answer as $n$, the estimate is near exact;
• Power series or trigonometric series are common (see table below).

Table 11-1 · Common series satisfying the displacement BCs

Using the prismatic cantilever of Animation 11.4.4-a, try four different trial functions to compute the critical load and analyze the results. Displacement boundary conditions: $y = 0,\ y' = 0$ at $x = 0$.

(1) First term of a cosine series · $y = a_{1}(1 - \cos \pi x/2l)$

Derivatives: $y' = (\pi a_{1}/2l)\sin(\pi x/2l)$, $y'' = (\pi^{2} a_{1}/4 l^{2})\cos(\pi x/2l)$.

From $\mathrm d E_{\mathrm P}/\mathrm d a_{1} = 0$ with $a_{1} \ne 0$:

(2) Cubic polynomial · $y = a_{1}(x^{2} - x^{3}/3l)$

Generic form of the deflection under a tip transverse force. Derivatives:

Stationary condition gives

(3) Two coordinates · $y = a_{1} x^{2} + a_{2} x^{3}$

$y' = 2 a_{1} x + 3 a_{2} x^{2}$, $y'' = 2 a_{1} + 6 a_{2} x$. Thus

Let $\alpha = F_{\mathrm P}\, l^{2}/EI$; the stationary conditions give

Setting the determinant to zero: $3\alpha^{2} - 104\alpha + 240 = 0$, smallest root $\alpha = 2.486$, so

(4) Quadratic · $y = a_{1} x^{2}$ (a warning)

Here $y'' = 2 a_{1}$ is constant — the curvature is constant, far from the real cantilever's curvature (which varies as $\cos$).

Applying the stationary condition gives

Accuracy comparison

Recompute Ex. 11-3 (portal frame) of §11-3 by the energy method — columns of equal $EI$, beam of $EA = \infty$, column tops loaded by $F_{\mathrm P}$, sway buckling mode. The static method gave $F_{\mathrm{Pcr}} = 2.706\, EI/l^{2}$; compare.

① Trial function (CD column)

During sway buckling (see Animation 11.4.5-a), the strain energy comes from bending of the compressed $CD$ and the unloaded $EF$; the load potential is due to the displacements at $B$ and $D$.

Take for $CD$ the deflection shape generated by a transverse force at D (same form as Ex. 11-5 variant (2)):

② Strain energy · EF via its lateral stiffness

The strain energy of $EF$ can be found by computing the external work on $F$'s sway $\Delta$. From the chosen trial function $\Delta = \tfrac{2}{3}\, l^{2}\, a_{1}$. $EF$'s lateral stiffness is $k = 3 EI/l^{3}$, so

Total strain energy of both columns:

③ Load potential · two contributions

At buckling $AB$ undergoes only rigid-body rotation of angle $\Delta/l$; thus the potential of $F_{\mathrm P}$ at $B$ is

The potential of $F_{\mathrm P}$ acting on $CD$ follows Eq. (11-17) (same form as Ex. 11-5 variant (2)):

④ Total potential · Stationary condition · Critical load

Total potential:

From $\mathrm d E_{\mathrm P}/\mathrm d a_{1} = 0$:

With $a_{1} \ne 0$ we solve

$F_{\mathrm{Pcr}}$ from the true buckled shape is the exact value; $F_{\mathrm{Pcr}}$ from any assumed shape is generally higher than the exact value. The reason: the assumed shape is only a subset of all possible buckled configurations — effectively imposing additional constraints, which raises the system's resistance to buckling.

Improving accuracy

1. Choose a trial function closer to the real shape — the cantilever's cosine first term is exact;
2. Add more terms — if $n+1$ terms give essentially the same answer as $n$, the estimate is near exact;
3. The trial function must satisfy the displacement BCs; satisfying force BCs is a bonus;
4. Ritz gives an upper-bound estimate of $F_{\mathrm{Pcr}}$ — unsafe in design; apply a safety factor.

Animation 11.4.7-a: a pinned-pinned prismatic column $AB$ under a uniformly distributed self-weight $q$ and an axial concentrated load $q\, l$ at point $C$ (at $2l/3$ from the bottom end). Use the energy method to find the critical self-weight intensity $q_{\mathrm{cr}}$.

① Trial function (satisfying pinned-pinned BCs)

Satisfying $y = 0$ at $x = 0$ and $x = l$ (see Animation 11.4.7-a), take

Derivatives:

② Strain energy

③ Potential of the distributed self-weight

From the figure, the axial shortening of an elemental length $\mathrm d x$ due to slope $\theta$ is $\tfrac{1}{2}(y')^{2}\,\mathrm d x$. Under the distributed self-weight $q$, the weight carried above the segment at $x$ is $q(l - x)$, so

Integrating over $[0, l]$:

④ Potential of the concentrated load $q l$ at $C$

The concentrated force $q l$ only travels through the axial shortening of the segment $[0, 2l/3]$ above it:

⑤ Total potential · Stationary condition · Critical load

From $\mathrm d E_{\mathrm P}/\mathrm d a = 0$ with $a \ne 0$: