OpenSees Cloud

OpenSees AMI

A Marathon, Not a Sprint

Original Post - 21 Aug 2020 - Michael H. Scott

Visit Structural Analysis Is Simple on Substack.

When tasked with developing an OpenSees model for simulating the nonlinear dynamic response of, let’s say, a multi-story reinforced concrete frame, you may be tempted to go straight to force-based frame elements with fiber sections comprised of Concrete23 and Steel08 material models.

This sprint to the finish line will undoubtedly lead to an analysis that “fails”.

To increase the likelihood of achieving a successful analysis, I recommend you follow these steps when building a model in OpenSees. I assume you are not brand new to OpenSees and that you know the basic commands to perform a nonlinear dynamic analysis for an earthquake ground motion.

1 – Use variables to define units

Although one of my colleagues in Eastchester does not strongly back this approach, I like to define variables for the base units of force, length, and time, then define physical constants and derive other units and unit conversions from these three variables. There’s no need to do unit conversions off line and enter magic numbers–you’ll make a mistake.

# Base units
kip = 1
inch = 1
sec = 1

# Other units
ft = 12*inch
ksi = kip/inch**2
psi = ksi/1000

# Physical constants
g = 32.2*ft/sec**2

# Unit conversions
m = 3.281*ft
kN = 0.2248*kip
kPa = kN/m**2

Above is not an exhaustive list, but you get the idea. If you define the input parameters using the units variables, all analysis computations and output will be in the basic units. To make a unit-dependent calculation or to obtain output in a different unit, just divide by the variable for that unit.

Es = 29000*ksi
fc = 4000*psi

Ec = 57000*(fc/psi)**0.5 * psi
print(f'Ec = {Ec/kPa} kPa') # Python f-string

2 – Define your model using nonlinear elements with elastic sections

After you define nodes and boundary conditions, define the elements with elastic sections (E, A, I, etc.) in a nonlinear formulation (dispBeamColumn, forceBeamColumn, etc.) with linear geometric transformations, which are easily changed later if large displacements are important. If the model is 3D, make sure you get the vecxz input correct.

Plot the model (easy using opsvis plotting commands) to make sure the element connectivity is correct. Go ahead and define any multi-point constraints, e.g., rigid diaphragms and equalDOFs, in this step too.

3 – Apply gravity loads

Put gravity loads in a constant time series and perform one step of a static analysis using load control with \(\Delta t\)=0. Remember that member loads for frame elements are defined with respect to local element axes, not global axes.

ops.timeSeries('Constant',1)
ops.pattern('Plain',1,1)
ops.load(...) # nodal gravity load
ops.eleLoad(...) # element gravity load

ops.integrator('LoadControl',0.0)
ops.analysis('Static')
ops.analyze(1)
ops.reactions()

If the analysis fails because of the equation solver, you have a rigid body mode (missing boundary condition or unconnected nodal DOF) that needs to be resolved before you move on.

Print out the reactions and do some equilibrium checks. What goes up must equal what comes down. Plot the displaced shape to make sure there’s nothing bizarre.

For models with large gravity loads, you can define the gravity loads in a linear time series and perform multiple steps of a load-controlled static analysis, then freeze the loads and revert the domain to time 0.

ops.timeSeries('Linear',1)
ops.pattern('Plain',1,1)
ops.load(...) # nodal gravity load
ops.eleLoad(...) # element gravity load

Nsteps = 10
ops.integrator('LoadControl',1/Nsteps)
ops.analysis('Static')
ops.analyze(Nsteps)
ops.reactions()

ops.loadConst('-time',0)

Either way, this is the initial condition for your dynamic analysis.

4 – Define mass

Although some structural analysis software defines mass automatically based on gravity loads, OpenSees does not do this. Use that gravitational constant from step 1 to go between weight and mass.

5 – Do an eigenvalue analysis

An eigenvalue analysis is a good way to make sure the distribution of mass and stiffness is correct. Does the fundamental period make sense for the height of your building? Plot the mode shapes too. You can also use the eigenvalues to define either Rayleigh damping or modal damping.

Nmodes = 4 # or whatever
w = ops.eigen(Nmodes)
T1 = 2*np.pi/w[0]**0.5
print(f'Fundamental period, T = {T1} sec')

# Rayleigh damping
zi = 0.05 # Mode 1
zj = 0.03 # Mode 3
wi = w[0]**0.5
wj = w[2]**0.5
A = np.zeros((2,2))
A[0,0] = 1/wi; A[0,1] = wi
A[1,0] = 1/wj; A[1,1] = wj
b = np.zeros(2)
b[0] = zi
b[1] = zj
x = np.linalg.solve(0.5*A,b)

ops.rayleigh(x[0],0,0,x[1])

6 – Define your recorders

Give some thought to what you want to record, both globally (nodal displacements, base shear, etc.) and locally, e.g., moment-curvature. Also think about how you name your recorder output files.

7 – Run dynamic analysis with earthquake ground motion

Plot the response history. Does it make sense? Does the fundamental period of vibration shine through in the response?

8 – Change from elastic sections to fiber sections with elastic materials

After the initial definition of nodes (step 2), this is the second biggest step for your model. Change the elastic sections to fiber-discretized cross-sections with elastic uniaxial materials. Remember, you don’t need a lot of fibers. Repeat steps 3, 5, and 7, making sure you get the same response, but with minor changes due to fiber discretization.

If you want to record fiber response, define those recorders now.

9 – Switch the elastic materials to simple inelastic materials

To ease into material nonlinearity, switch the concrete to ElasticNoTension (ENT) and make the steel Steel01. Repeat steps 3, 5, and 7 and be able to explain why the response changes. You shouldn’t run into any convergence problems with this step. If you do, something’s up with your model.

Like mile 20 in a marathon, you are now at the mental half way point of your OpenSees analysis. Replenish those electrolytes, then get ready for the next 6.2 miles.

10 – Switch the inelastic materials from simple to moderate

Before making the leap to Concrete23 and Steel08, switch the concrete from ENT to Concrete01. Repeat steps 3, 5, and 7 and notice how the response changes. It’s very likely you’ll run into convergence problems after the switch to Concrete01. There’s nothing inherently wrong with Concrete01, it just introduces a lot of material nonlinearity. If you haven’t already done so, you may need to ramp up the gravity loads instead of using a constant time series as described in step 3.

11 – Diagnose convergence issues

If your analysis fails to converge, you can try different time steps, algorithms, constraint handlers, and equation solvers (if the system becomes non positive definite). It’s best to address any convergence issues with these “moderate” material models before moving on to the complex material models in the next step.

If you do have convergence problems, pay attention to the norms coming out of the convergence test. Is the norm converging to something or is the norm bouncing around? Whatever you do, don’t use the linear algorithm–that will not solve the convergence problem. Switch to a FullGeneral system of equations and obtain the dynamic tangent stiffness matrix and run some diagnostics in Python or MATLAB.

Nsteps = 2000
dt = 0.01*sec

for i in range(Nsteps):
   ok = ops.analyze(1,dt)
   if ok < 0:
      # Analysis failed, switch system and get KT
      ops.system('FullGeneral')
      ops.analyze(1,dt)
      KT = ops.printA('-ret')
      KT = np.array(KT)
      N = ops.systemSize()
      KT.shape = (N,N)
      # Dump KT to a file and then run some diagnostics
      break

You can also use the GimmeMCK integrator to get the individual mass, stiffness, and damping matrices.

12 – Switch the inelastic materials from moderate to complex

Now make the switch to Concrete23 and Steel08. Repeat steps 3, 5, and 7 and start working your way through the results. Are the results terribly different from step 10? If not, ask yourself why you wanted to use Concrete23 and Steel08 in the first place. If the displacements are large, change the column geometric transformations to PDelta or Corotational.

13 – Fine tuning

If you’re using fiber sections with Concrete23 and Steel08, there’s not much you can do to speed up your analysis besides using fewer fibers in each section. But, with a functioning model and analysis, you can play with solvers and look at other improvements in your model.

If you set up your script so that you can easily swap out nonlinear elements for their elastic counterparts, i.e., to move freely between steps 8, 9, 10, and 12, it will be much easier to diagnose convergence problems. The concepts and steps described here will be roughly the same for nonlinear static analysis (pushover) and for models of steel structures and soil models.