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 these OpenSeesPy 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 . 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 very small 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.

You can have numerical problems in one unit, but do not in other. Maybe unit systems are also important

LikeLiked by 1 person

This is probably my favorite post!

LikeLiked by 1 person

Thanks! And this post is one of my favorites 🙂

https://portwooddigital.com/2021/05/29/right-under-your-nose/

LikeLiked by 1 person