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., .

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., . As you apply larger reference loads, the geometric stiffness will increase.

Solving for the material and geometric stiffness matrices, we obtain and .

Now we can perform generalized eigenvalue analysis with and , i.e., solve 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 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.

Update: Here is a discussion from the OpenSees Cafe on the frame buckling procedure.

Actually you gave this hint in one of your posts. But you dismissed off the idea saying that we will have a tough time mapping eigen vectors back.

https://portwooddigital.com/2020/11/13/ordinary-eigenvalues/

Really “Right under your nose” !

LikeLike

Hello ash,

I don’t think that’s what I said in the “Ordinary Eigenvalues” post. What I meant was perhaps you could fake a geometric stiffness by defining nodal mass in some clever way, then use the eigenvalue solver in OpenSees to get buckling loads. But what Luigi came up with is the ideal solution.

PD

LikeLike

Thank you. Now I understood !!!

LikeLike

An attractive blog post

LikeLiked by 1 person

Some time ago (~2009) I implemented this algorithm in our OpenSees fork. You can take a look here:

https://github.com/xcfem/xc/blob/master/verif/tests/solution/eigenvalues/linear_buckling_column01.py

If someone is interested in doing the backport to OpenSees I’ll be glad to assist.

LikeLike

Hello Mr. Scott. I was trying to perform the same procedure for buckling eigenvalue analysis for models with shell elements and when I have seen the article it has helped me a lot to settle my ideas. The triangular element ShellNLDKGT, for example, has a contribution to the stiffness matrix of geometric nonlinearities related to large deformations. I am not completely sure if this same procedure can be used in this case. What are your thoughts?. Thank you very much in advance. Kind Regards.

LikeLike

Hola Alfonso,

It should work if the element accounts for large displacements. Give it a try and report back if it works or not.

PD

LikeLike

Hello Mr. Scott. Thanks to this article I have noticed what I was doing wrong before. I was obtaining Kgeo as KB-KA instead of KA-KB, so the sign of lambda factors were wrong. In order to compare the eigenvalue buckling analysis I have made the same model of section 2.4.4 in the article “Development and Application of a High-Performance Triangular Shell Element and an Explicit Algorithm in Opensees for Strongly Nonlinear Analysis” which is the article of the implementation of ShellNLDKGT in OpenSees. With the correction of the geometric stiffness matrix, I have obtained a first critical load similar to that obtained in the article (~200 kN). I have also performed a eigenvalue buckling analysis of the same model in ansys and I have obtained slightly different results for the first critical load (~250 kN) but it must be taken into account that the ansys model has been made with quadrilateral elements and they are not the same as those used in OpenSees. In the absence of further testing and modeling different configurations, it seems that this procedure could also work for buckling eigenvalues analysis in models with shell elements. Thank you very much for everything and best regards.

LikeLike

Hola Alfonso,

Great to hear that it works! Please send me the input file so that I can share it more broadly.

PD

LikeLike