Minimal Plate Buckling Example

OpenSees is not built to perform linear buckling analysis. But a few years ago, Luigi Caglio shared a workaround described in this post.

In the post, the example application is a frame model, but there’s no reason the approach cannot work for shell models. So, here’s a minimal working example.

Consider a rectangular steel plate with simple boundary conditions and compressive axial load. The elastic modulus is E=29,000 ksi and the Poisson ratio is 0.3. The plate length, width, and thickness are L=50 inch, h=8 inch, and t=1 inch, respectively.

A mesh of geometrically nonlinear ShellNLDKGT elements is defined below using the mesh command shown in this post.

import openseespy.opensees as ops
import numpy as np
import scipy.linalg as slin

# Plate dimensions
L = 50.0
h = 8.0
t = 1.0

# Material properties
E = 29000
v = 0.3

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

ops.node(1,0,-h/2,0)
ops.node(2,0,h/2,0)
ops.node(3,0,h/2,L)
ops.node(4,0,-h/2,L)

c = h/4 # Mesh size

ops.mesh('line',1,2,*[1,2],0,6,c)
ops.mesh('line',2,2,*[2,3],0,6,c)
ops.mesh('line',3,2,*[3,4],0,6,c)
ops.mesh('line',4,2,*[4,1],0,6,c)

ops.section('ElasticMembranePlateSection',1,E,v,t)

ops.mesh('quad',13,4,*[4,3,2,1],0,6,c,'ShellNLDKGT',1)

ops.fixZ(0,1,1,1,0,0,0) # bottom nodes, Z=0
ops.fixZ(L,1,1,0,0,0,0) # top nodes, Z=L

ops.fix(1,0,0,0,1,0,0)

The last fix command restrains one drilling DOF in the model in order to avoid artificially low eigenvalues.

Before applying any loads, we can get the material stiffness for the plate using the printA command.

ops.integrator('LoadControl',0.0)
ops.system('FullGeneral')
ops.analysis('Static','-noWarnings')
ops.analyze(1)

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

We can now apply a unit reference load, get the stiffness matrix again accounting for geometric nonlinearity, then obtain the geometric stiffness matrix by subtracting off the material stiffness.

endNodes = ops.getNodeTags('-mesh',3) # tag 3 = line mesh across top of plate
Pref = 1
dP = Pref/len(endNodes)

ops.timeSeries('Constant',1)
ops.pattern('Plain',1,1)
for nd in endNodes:
    ops.load(nd,0,0,-dP,0,0,0)

ops.analyze(1)

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

Kgeo = Kmat-Kgeo

Finally, we call the eig() function from the scipy.linalg package, then find the minimum positive eigenvalue (the critical buckling load) and the index of the corresponding eigenvector which will correspond to the buckled shape.

lam, x = slin.eig(Kmat,Kgeo)
lamPos = lam[lam>0]
xPos = x[:,lam>0]

lamMin = min(lamPos)
lamMinIndex = lamPos.argmin()

For this analysis, the computed buckling load (lowest eigenvalue) is 76.7 kip while the exact flexural buckling load is =76.3 kip where I=ht³/12.

After some low-level gymnastics mapping equation numbers from the critical eigenvector back to nodes of the OpenSees model, we can view the buckled shape.

Looks correct to me!