24–15 Solutions to Exercises

Homework Exercises for Chapter 24

Solutions

EXERCISE 24.1
Because K = K0 + λK1 is diagonal, the null eigenvector (buckling mode) z for the lowest critical load and
the incremental load vector q are    
1 0
z= , q= (E24.9)
0 1
Obviously z and q are orthogonal.

EXERCISE 24.2
The linear equations for element (1) are the same as in §9.6, whereas for element (2) we have
    
c2 sc −c2 −sc u X2 0
 sc s2 −sc −s 2   u Y 2   0 
k (2)  2 = (E24.10)
−c −sc c2 sc   u X 3   0 
−sc −s 2 sc s2 uY 3 0

where s = sin ϕ and c = cos ϕ. Assembling and imposing boundary conditions we find that the LPB stress
distribution is the same as in the case ϕ = 0◦ of §9.6, that is, bar (1) takes all the load. This happens because
the structure is still statically determinate.
The LPB eigensystem is
 λq 
k (2) c2 − k (2) sc    
 L (1)  u X2 = 0 .
0
λq (E24.11)
k (2) sc k (1) + k (2) s 2 − uY 2 0
L (1)
0

The characteristic equation obtained by setting detK = 0 is

λq λq
(k (2) c2 − )(k (1) + k (2) s 2 − ) − k (2) s 2 c2 = 0. (E24.12)
L (1)
0 L (1)
0

Dividing through by (k (1) )2 and calling

q
α=λ , r = k (2) /k (1) , (E24.13)
k (1) L (1)
0

we get the quadratic equation

(c2 r − α)(1 + r s 2 − α) − r 2 s 2 c2 = 0, (E24.14)
or
c2 r − α(1 + r s 2 ) + α 2 = 0. (E24.15)
The roots are
α1,2 = 1
2
1 + r s2 ± (1 + r s 2 )2 − 4c2 r . (E24.16)

We are interested in the smallest root α1 obtained by taking the − sign. This root has the small-r expansion

α1 = c2 r + O(r 2 ) (E24.17)

24–15
Chapter 24: BIFURCATION: LINEARIZED PREBUCKLING I 24–16

Retaining only the first term we obtain as r → 0

k (1) L (1) k (2) L (1)

λcr = α1 0
= 0
c2 (E24.18)
q q

which agrees with the crtical value in (24.21) for ϕ = 0◦ . In a similar manner one obtaines that the limit
buckling mode is the same as (E24.2) thus verifying orthogonality.

EXERCISE 24.3
The axial load in bar (1) is p = −λcu 2Y 2 . Consequently the LPB eigensystem is
 
(2) λcu 2
k − (1)Y 2 0
 L0 
det 
λcu 2
=0 (E24.19)
0 k (1) − (1)Y 2
L0

This does not fit within the generalized eigenvalue problem (24.11) in which both matrices are constant. It
must be categorized as a nonlinear eigenvalue problem, which pertains to the general type (24.3).

EXERCISE 24.4

(a) The two translational rigid body motions along the x and y directions are

u x = cx , v = 0, θz = 0, u = 0, u y = cy , θz = 0 (E24.20)

where cx and c y are arbitrary constants. The corresponding nodal displacement vector is

urig = [ cx cy 0 cx cy 0 ]T (E24.21)

The check is carried out by showing that the associated nodal forces are identically zero for any cx and
cy :
K0 urig = 0, K1 urig = 0 (E24.22)
This can be verified by simple algebra. Alternatively one may postmultiply the two RBM basis vectors
[ 1 0 0 1 0 0 ]T and [ 0 1 0 0 1 0 ]T . By inspection the sum of columns 1 and 4 and
columns 2 and 5 of both K0 and K1 are zero.
Note that infinitesimal rigid-body rotation modes such as

utilt = [ 0 − 12 vr vr /L 0 v
1 r
2
vr /L ]T (E24.23)

which corresponds to the beam “tilting” about its midpoint, are not zero-force modes for the geometric
stiffness unless the axial force N vanishes. The product of K1 by (E24.23) produces a system of
equilibrium-restoring (if N > 0) or equilibrium-upsetting (if N < 0) nodal moments.

(b) The LPB eigenproblem is (E24.4) with the matrices given by (E24.2) and (E24.3). Upon applying the
boundary conditions u X 1 = u Y 1 = u X 2 = 0 this 6 × 6 system reduces to a 3 × 3 eigenproblem:
   
4E I /L 0 2E I /L λP 4L 2 0 −L 2
E A/L 0 + 0 0 z=0 (E24.24)
symm 4E I /L 30L symm 4L 2

24–16
24–17 Solutions to Exercises

in which  
θz1
z= u X2 (E24.25)
θz2

(c) The second equation in (E24.24) is

(E A/L)u X 2 = 0 (E24.26)
Since E A/L = 0 we must have u X 2 = 0. We can therefore treat this as a displacement boundary
condition and remove u X 2 from (24.23). This leaves us a 2 × 2 eigenproblem with z containing only the
two end rotations:
  λP  2   
4E I /L 2E I /L 4L −L 2 θ
+ z = 0, z = z1 (E24.27)
2E I /L 4E I /L 30L −L 2
4L 2 θz2

which can be further simplified to

 
EI 4 2
 λP L 
4 −1
+ z=0 (E24.28)
L 2 4 30 −1 4

This may be solved directly as a quadratic equation in λ, which yields the roots λ1 = −12 and λ2 =
−60. For hand solution it is more instructive, however, to split the solution into a symmetric and an
antisymmetric mode as suggested.

(d) For the assumed symmetric buckling mode,

   
θ 1
z1 = =θ (E24.29)
−θ −1

Inserting this into (E24.28) reduces it to the scalar eigenproblem

 
EI λP L
2+ 5 θ =0 (E24.30)
L 30

The expression in parenthesis must vanish for the root λ = λ1 and a simple calculation yields

EI EI
λ1 = −12 , pcr 1 = λ1 P = −12 . (E24.31)
P L2 L2

(e) For the antisymmetric mode,    

θ 1
z2 = =θ (E24.32)
θ 1
the same calculation scheme yields

EI EI
λ2 = −60 , pcr 2 = λ2 P = −60 . (E24.33)
P L2 L2

(f) The two lowest exact critical loads for the continuum problem are (E24.5). Since π = 3.14159265 we
have
EI EI
exact
pcr 1 = −9.869604 2 , exact
pcr 2 = −39.478417 2 . (E24.34)
L L

24–17
Chapter 24: BIFURCATION: LINEARIZED PREBUCKLING I 24–18

Comparing with (E24.31) and (E24.33) we find that the errors for the symmetric and antisymmetric
buckling cases are approximately 22% and 52%, respectively. Both eigenvalues are overestimated,
which is as it should be according to the Rayleigh-Ritz theory2 and the Courant-Fisher theorem.

(g) The sequence of calculations for the reduced-integration geometric stiffnesses (E24.6)–(E24.7) is identi-
cal. The results are summarized in the following table, which also includes the values previously obtained
for the exactly integrated geometric stiffness.

Case Coefficient K1 K∗1 K∗∗

1
Symmetric λ1 /(E I /P L 2 ) −12 −8 −12
Antisymmetric λ2 /(E I /P L 2 ) −60 ∞ ∞

A value of ∞ reports that the denominator in the expression for λ, which comes from the geometric
stiffness term, vanishes. This occurs with the one and two-point integrated geometric stiffnesses for the
antisymmetric mode3 . The exactly integrated K1 gives a finite, though inaccurate, value for λ2 . The
value for λ1 given by K∗1 is slightly better than the other two but still well outside engineering accuracy.
Note that this value of -8 underestimates the true value −π 2 because the Rayleigh-Ritz bound theorem
does not hold if one or both stiffness matrices are approximately integrated.
(h) A 2-element, 3-node model may be processed in two ways. Two elements of length 12 L are assembled.
The 9 degrees of freedom are reduced to 6 by imposing the 3 support constraints, and two axial freedoms
at nodes 2 and 3 are decoupled, reducing the determination of critical loads to a 4 by 4 eigenproblem.
This may be split into two 2 by 2 using symmetry and antisymmetry constraints.
Work can be reduced significantly (especially for hand computations) by imposing the symmetry and
antisymmetry conditions from the start. One can then work throughout with a single half-length element,
say 1-2, without need to do an assembly.
For the symmetric buckling we impose u x1 = u y1 = 0 from support conditions, u x2 = 0 from uncoupling,
and θz2 = 0 from symmetry, leaving two active freedoms: θz1 and u y2 . The eigensystem is
  λP  2   
8E I /L −12E I /L 2 8L −6L θz1
+ z = 0, z= (E24.35)
−12E I /L 24E I /L 3 30L −6L 72 u y2

Solving this with Mathematica gives the roots λ1 = −9.9438E I /(P L 2 ) √

and λ1 = −128.72E I /(P L 2 ),
where the coefficients are the 5-place numerical values of (−416 ± 64 2)/6. Hence the symmetric
critical load is
EI
pcr 1 = λ1 P = −9.9438 2 (E24.36)
L
The error is only 0.75%4

2 The exactly-integrated Hermitian element for the BE beam model fits within the framework of the Rayleigh-Ritz approx-
imation theory.
3 This mode is a cubic with odd variation about the beam midpoint, as depicted in Figure E24.3. One needs 3 Gauss points
to capture this polynomial variation.
4 The book I took this exercise from, as well as the original 1965 paper by H. C. Martin that presented this analysis for
the first time, report a coefficient of 10 instead of 9.9438. The small difference is not clear since they do not show the
source eigensystem. But I have checked that the elemental K0 and K1 are the same as they used.

24–18
24–19 Solutions to Exercises

The antisymmetric case can be also reduced to a 2 by 2 eigensystem by setting u x1 = u y1 = u x2 = u y2 =

0. The smallest critical load given by this system is

EI
pcr 2 = λ2 = −48 . (E24.37)
L2

which is in error by about 22% percent with respect to −4π 2 E I /L 2 .

EXERCISE 24.5

EXERCISE 24.6
W comes from the expansion of λP times the axial shortening of the column, expanded in Taylor series about
u 1 = u 2 = 0, and retaining only the quadratic terms.
The critical load factors are λ S = P/k, which corresponds to the symmetric buckling mode zTS = [1, 1]
and λ A = 3P/k, which corresponds to the antisymmetric buckling mode zTA = [1, −1]. Here is a solution
computed using Mathematica 2.2.

ClearAll[U,V,PE,u1,u2,k,lambda];
U=(k/2)*(u1^2+u2^2);
V=(lambda*P/2)*(u1^2+(u1-u2)^2+u2^2);
PE=U-V;
r={D[PE,u1],D[PE,u2]}; r=Simplify[r];
K={D[r,u1],D[r,u2]}; K=Simplify[K];
Print["K=",K];
KM=Coefficient[K,lambda,0];
KG=Coefficient[K,lambda,1];
A=-Inverse[KM].KG;
{eval,evec}=Eigensystem[A];
Print["eigenvalues =",eval];
Print["eigenvectors=",evec];

K={{k - 2 lambda P, lambda P}, {lambda P, k - 2 lambda P}}

P 3 P
eigenvalues ={-, ---}
k k
eigenvectors={{1, 1}, {-1, 1}}

24–19