OpenSees Cloud

OpenSees AMI

Right Under Your Nose

Original Post - 29 May 2021 - Michael H. Scott

Visit Structural Analysis Is Simple on Substack.

I long ago accepted that buckling analysis would never be implemented in OpenSees.

Although there is a getGeometricTangentStiff() method in the Element interface, only PFEM fluid elements use it. Implementing this method for frame elements, assembling the geometric stiffness, and solving the generalized eigenvalue problem would require several updates to the innards of OpenSees.

Then Luigi Caglio, a Ph.D. student in civil engineering at the Technical University of Denmark (DTU), came up with a clever solution, posted in the OpenSees Facebook group.

I will do my best to explain Luigi’s idea, but the gist is to use the printA function twice to capture the material and geometric stiffness matrices of a model, then do eigenvalue analysis in your script (Python makes it easy, but you could use Tcl as well).

Define a frame model with geometric nonlinearity. Typically, a corotational mesh of geometrically linear elements, but you could also use geometrically nonlinear elements. The elements should be linear-elastic, e.g., elasticBeamColumn or forceBeamColumn with elastic sections.

First, perform a static analysis with no loads so that the stiffness matrix is assembled. Then obtain the stiffness matrix using printA. Remember to use the FullGeneral system of equations. This gives the material stiffness of the model, i.e., \({\bf K}_A={\bf K}_{mat}\).

Next, apply reference loads and perform a second static analysis. Then obtain the stiffness matrix again. This gives the material stiffness of the model minus the geometric stiffness, i.e., \({\bf K}_B={\bf K}_{mat}-{\bf K}_{geo}\). As you apply larger reference loads, the geometric stiffness will increase.

Solving for the material and geometric stiffness matrices, we obtain \({\bf K}_{mat} = {\bf K}_A\) and \({\bf K}_{geo}={\bf K}_A-{\bf K}_B\).

Now we can perform generalized eigenvalue analysis with {\bf K}{mat} and \({\bf K}_{geo}, i.e., solve {\bf K}_{mat}{\bf x}=P{\bf K}_{geo}{\bf x}\) where _P is the eigenvalue (buckling load) and x is the eigenvector (buckled shape).

Because the analyses are static, you don’t need to use the GimmeMCK integrator to get the stiffness matrix–LoadControl with \(\Delta t=0\) will suffice.

import numpy as np
import scipy.linalg as slin

#
# Define your model
#

ops.system('FullGeneral')
ops.integrator('LoadControl',0)
ops.analysis('Static')

ops.analyze(1)

N = ops.systemSize()
Kmat = ops.printA('-ret')
Kmat = np.array(Kmat)
Kmat.shape = (N,N)

ops.timeSeries('Constant',1)
# Reference loads for buckling
ops.pattern('Plain',1,1)
ops.load(2,0,-1.0,0)
ops.load(3,0,-1.0,0)

ops.analyze(1)

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

# Eigenvalue solution
lam, x = slin.eig(Kmat,Kgeo)

# Sort the positive eigenvalues
lamSort = np.sort(lam[lam>0])

Below is an example that I use in my nonlinear structural analysis course and that was shown to me when I took nonlinear structural analysis.

In a nod to reality, the columns in this frame are not slender enough for elastic buckling to precede material yielding–not by a long shot. But we’ll keep going with elastic buckling.

Assuming axial rigidity of all members, we can form the 3×3 stiffness matrix (combined material and geometric) for this frame using direct assembly with exact stability functions based on the column axial loads in each load case. The gory details of that assembly are omitted from this post.

Shown below is the determinant of that 3×3 matrix as a function of the load, P, along with the eigenvalues obtained from Luigi’s method.

The results match very well! Differences can be attributed to the inclusion of axial stiffness and the number of corotational elements per member in the OpenSees model, as well as, for the higher modes, the inherent linearization of eigenvalue problems. For Load Case 2, the determinant switches from positive to negative infinity, so you can ignore the sharp vertical line at about P=27,000 kip.

Although not shown here, you can check that the eigenvectors returned by slin.eig match the buckled shape of the frame.

One of the nice things about this approach is you don’t have to use initial sway to trick the columns in to buckling. In addition, although I only showed one story frames in this post, there’s no reason this approach will not work for multi-story frames.

So, the solution for elastic buckling analysis has been in OpenSees all along, hiding in plain sight.

Grazie mille, Luigi Caglio, for finding it.

Thank you to Luigi Caglio and Silvia Mazzoni for comments and feedback on an initial draft of this post.