# Frame Buckling Analysis in OpenSeesPy

### Authors: Luigi Caglio and Michael H. Scott

June 2021

This notebook describes the elastic frame buckling computations developed by Luigi Caglio.
The original idea was posted as a reply to a post in the OpenSees Facebook group and the steps were then summarized in this blog post.

## Computational Procedure

As described in the blog post, the following steps are required in order to compute buckling loads.

### 1. Define the model

Define a frame model with geometric nonlinearity. This can be either a corotational mesh of elastic elements (either elasticBeamColumn or disp/forceBeamColumn with elastic sections) or elements with geometric nonlinearity "inside the basic system" (forceBeamColumnCBDI, dispBeamColumnNL, mixedBeamColumn, etc. with elastic sections).

### 2. Obtain the material stiffness matrix

Before applying any loads to the model, use printA to get the stiffness matrix.

### 3. Obtain the geometric stiffness matrix

Now apply reference axial loads to the model, use printA again to get the geometric stiffness matrix.

### 4. Eigenvalue analysis

Solve the generalized eigenvalue problem using numpy.slin.eig(). The returned eigenvalues give the critical buckling loads and the eigenvectors give the buckled shapes.

### Example Frame Structure¶

The computational steps will be demonstrated for the frame shown below.

In [56]:
import openseespy.opensees as ops
import openseespy.preprocessing.DiscretizeMember as dm

import matplotlib.pyplot as plt
import openseespy.postprocessing.ops_vis as opsv

import numpy as np
import scipy.linalg as slin

#plt.rc('text',usetex=True)
#plt.rc('font',family='serif')

%matplotlib notebook


### Units and Material Properties

In [57]:
inch = 1
kip = 1

ft = 12*inch
ksi = kip/inch**2

E = 29000*ksi


### Main OpenSees code¶

In [58]:
# Creating a figure
plt.figure()

# Loop through two load cases
for lc in [1,2]:

ops.wipe()
ops.model('basic','-ndm',2,'-ndf',3)

L = 12*ft
H = 30*ft
Lc = L
Lg = H

ops.node(1,0,0); ops.fix(1,1,1,1)
ops.node(2,0,L)
ops.node(3,H,L)
ops.node(4,H,0); ops.fix(4,1,1,1)

# W14x90
A = 26.5
Ic = 999
ops.section('Elastic',1,E,A,Ic)
ops.beamIntegration('Legendre',1,1,2)

# W18x76
A = 22.3
Ig = 1330
ops.section('Elastic',2,E,A,Ig)
ops.beamIntegration('Legendre',2,2,2)

# Corotational transformation
ops.geomTransf('Corotational',1)

# Number of elements/member
Nele = 8

# Columns
dm.DiscretizeMember(1,2,Nele,'dispBeamColumn',1,1,10,10)
dm.DiscretizeMember(3,4,Nele,'dispBeamColumn',1,1,30,30)

# Beam
dm.DiscretizeMember(2,3,Nele,'dispBeamColumn',2,1,20,20)

ops.numberer('RCM')
ops.system('FullGeneral')
ops.constraints('Plain')
ops.test('NormUnbalance',1e-8,10,0)
ops.algorithm('Newton')
ops.analysis('Static')

ops.analyze(1)

# Get number of equations
N = ops.systemSize()

# Get material stiffness
Kmat = ops.printA('-ret')
Kmat = np.array(Kmat)
Kmat.shape = (N,N)

ops.timeSeries('Constant',1)
ops.pattern('Plain',1,1)
if lc == 1:
if lc == 2:

ops.analyze(1)

# Get geometric stiffness
Kgeo = ops.printA('-ret')
Kgeo = np.array(Kgeo)
Kgeo.shape = (N,N)
Kgeo = Kmat-Kgeo

# Eigenvalue analysis, and sorting
lam, x = slin.eig(Kmat,Kgeo)
lamSort = np.sort(lam[lam>0])

# for plotting mode shapes later in the notebook
if lc == 1:
x_lc1 = x
lam_lc1 = lam
if lc == 2:
x_lc2 = x
lam_lc2 = lam

#
# Now get "exact" solution for this load case
# We are using direct assembly with exact stability functions for the columns
#
Kmat = np.zeros((3,3))

if lc == 1:
l1 = 1
l3 = 1
if lc == 2:
l1 = 2
l3 = 0

Pplot = []
Kplot = []

Pmax = 55000*kip
Nsteps = 1000
dP = Pmax/Nsteps

P = 0
while P <= Pmax:

# Element 1
psi = Lc*(l1*P/(E*Ic))**0.5
if psi > 0:
B1 = (psi*np.sin(psi)-psi**2*np.cos(psi))/(2*(1-np.cos(psi))-psi*np.sin(psi))
C1 = (psi**2-psi*np.sin(psi))/(2*(1-np.cos(psi))-psi*np.sin(psi))
else:
# Regular 4EI/L and 2EI/L from linear structural analysis
B1 = 4
C1 = 2

# Element 3
psi = Lc*(l3*P/(E*Ic))**0.5
if psi > 0:
B3 = (psi*np.sin(psi)-psi**2*np.cos(psi))/(2*(1-np.cos(psi))-psi*np.sin(psi))
C3 = (psi**2-psi*np.sin(psi))/(2*(1-np.cos(psi))-psi*np.sin(psi))
else:
# Regular 4EI/L and 2EI/L from linear structural analysis
B3 = 4
C3 = 2

# Direct assembly of stiffness matrix for this level of axial loads
Kmat[0,0] = 2*(B1+C1)*E*Ic/Lc**3 + 2*(B3+C3)*E*Ic/Lc**3 - l1*P/Lc - l3*P/Lc
Kmat[1,0] = (B1+C1)*E*Ic/Lc**2
Kmat[2,0] = (B3+C3)*E*Ic/Lc**2

Kmat[0,1] = Kmat[1,0]
Kmat[1,1] = B1*E*Ic/Lc + 4*E*Ig/Lg
Kmat[2,1] = 2*E*Ig/Lg

Kmat[0,2] = Kmat[2,0]
Kmat[1,2] = Kmat[2,1]
Kmat[2,2] = B3*E*Ic/Lc + 4*E*Ig/Lg

Pplot.append(P)
Kplot.append(np.linalg.det(Kmat)) # When determinant is zero, we know the value of P is critical

P += dP

plt.subplot(2,2,lc)
plt.plot(Pplot,Kplot,'b-',markersize=1,label='Direct')
Nmodes = 3
plt.plot(lamSort[:Nmodes],[0]*Nmodes,'ko',label='Eigen')
plt.grid()
plt.xlim([0,Pmax])
plt.xlabel('$P$ (kip)')
plt.ylim(-1.2*Kplot[0],1.2*Kplot[0])
if lc == 1:
plt.ylabel('$\det({\\bf K})$')
plt.legend()

C:\ProgramData\Anaconda3\lib\site-packages\numpy\core\_asarray.py:85: ComplexWarning: Casting complex values to real discards the imaginary part
return array(a, dtype, copy=False, order=order)
C:\ProgramData\Anaconda3\lib\site-packages\numpy\core\_asarray.py:85: ComplexWarning: Casting complex values to real discards the imaginary part
return array(a, dtype, copy=False, order=order)

In [59]:
# draw model:
import eSEESminiPy
eSEESminiPy.drawModel()


# Sort mode shapes for plotting¶

In [60]:
# function to sort lambda and modeshape vectors accordingly

def sortModes(lam,x):

#sort lambda and x
lambda_sorted = np.sort(lam)
lambda_inds = lam.argsort()
x_Sort = np.zeros_like(x)

for count, index in enumerate(lambda_inds):
x_Sort[:,count]=x[:,index]

#only modes with positive lambda
index = 0
while index<len(lambda_sorted) and lambda_sorted[index]<0:
index+=1
lamSort = lambda_sorted[index:]
x_Sort  = x_Sort[:,index:]

return lamSort,x_Sort

In [61]:
lam_lc1, x_lc1 = sortModes(lam_lc1,x_lc1)
lam_lc2, x_lc2 = sortModes(lam_lc2,x_lc2)


# Function to plot mode shapes¶

In [62]:
# function to plot the mode shapes
# from https://github.com/zhuminjie/OpenSeesPyDoc/blob/master/src/ops_vis/ops_vis.py

# default settings

# fmt: format string setting color, marker and linestyle
# check documentation on matplotlib's plot

# continuous interpolated shape line
fmt_interp = 'b-'  # blue solid line, no markers

# element end nodes
fmt_nodes = 'rs'  # red square markers, no line

# deformed model
fmt_defo = 'b-'  # green dashed line, no markers

# undeformed model
fmt_undefo = 'g:'  # green dotted line, no markers

# section forces
fmt_secforce = 'b-'  # blue solid line

# figure left right bottom top offsets
fig_lbrt = (.04, .04, .96, .96)

# azimuth and elevation in degrees
az_el = (-60., 30.)

# figure width and height in centimeters
fig_wi_he = (16., 10.)

#sfac=False

nep=17
unDefoFlag=1
fmt_undefo='g:'
interpFlag=1
endDispFlag=0
fmt_interp='b-'
fmt_nodes='rs'
Eo=0,
az_el=(-60.0, 30.0)
fig_wi_he=(16.0, 10.0)
fig_lbrt=(0.04, 0.04, 0.96, 0.96)

def plot_buckling_mode(mode,lambda_,vec,scalefac=10):  # modified from   _plot_defo_mode_2d    function
sfac=scalefac

ele_tags = ops.getEleTags()
#plt.figure()
for ele_tag in ele_tags:
nd1, nd2 = ops.eleNodes(ele_tag)

DOFs_nd1 = ops.nodeDOFs(nd1)
DOFs_nd2 = ops.nodeDOFs(nd2)
# element x, y coordinates
ex = np.array([ops.nodeCoord(nd1)[0],
ops.nodeCoord(nd2)[0]])
ey = np.array([ops.nodeCoord(nd1)[1],
ops.nodeCoord(nd2)[1]])

#if modeNo:
if DOFs_nd1[0] != -1:
ed1 = -vec[DOFs_nd1[0],mode]
else:
ed1 = 0
if DOFs_nd1[1] != -1:
ed2 = -vec[DOFs_nd1[1],mode]
else:
ed2 = 0
if DOFs_nd1[2] != -1:
ed3 = -vec[DOFs_nd1[2],mode]
else:
ed3 = 0
if DOFs_nd2[0] != -1:
ed4 = -vec[DOFs_nd2[0],mode]
else:
ed4 = 0
if DOFs_nd2[1] != -1:
ed5 = -vec[DOFs_nd2[1],mode]
else:
ed5 = 0
if DOFs_nd2[2] != -1:
ed6 = -vec[DOFs_nd2[2],mode]
else:
ed6 = 0

ed = np.array([ed1,
ed2,
ed3,
ed4,
ed5,
ed6])

if unDefoFlag:
plt.plot(ex, ey, fmt_undefo)

# interpolated displacement field
if interpFlag:
xcdi, ycdi = opsv.beam_defo_interp_2d(ex, ey, ed, sfac, nep)
plt.plot(xcdi, ycdi, fmt_interp)

# translations of ends
if endDispFlag:
xdi, ydi = opsv.beam_disp_ends(ex, ey, ed, sfac)
plt.plot(xdi, ydi, fmt_nodes)

plt.axis('equal')
plt.title("Mode "+str(mode+1)+'; lambda = '+str(int(np.real(lambda_[mode]))), x=0.5, y=0.8)
#plt.show()


# Plot mode shapes for load case 1¶

Notice the symmetry of the buckled shapes resulting from the symmetry of the loading

In [63]:
scale_factor = 50
lam = lam_lc1; x = x_lc1
plt.figure()
plt.subplot(2, 2, 1)
plot_buckling_mode(0,lam,x,scale_factor)   #mode 1
plt.subplot(2, 2, 2)
plot_buckling_mode(1,lam,x,scale_factor)   #mode 2
plt.subplot(2, 2, 3)
plot_buckling_mode(2,lam,x,scale_factor)   #mode 3
plt.subplot(2, 2, 4)
plot_buckling_mode(3,lam,x,scale_factor)   #mode 4
plt.show()


# Plot mode shapes for load case 2¶

Notice the asymmetry of the buckled shapes resulting from the asymmetric loading

In [64]:
scale_factor = 50
lam = lam_lc2; x = x_lc2
plt.figure()
plt.subplot(2, 2, 1)
plot_buckling_mode(0,lam,x,scale_factor)   #mode 1
plt.subplot(2, 2, 2)
plot_buckling_mode(1,lam,x,scale_factor)   #mode 2
plt.subplot(2, 2, 3)
plot_buckling_mode(2,lam,x,scale_factor)   #mode 3
plt.subplot(2, 2, 4)
plot_buckling_mode(3,lam,x,scale_factor*4)   #mode 4
plt.show()