Verifying Ain’t Easy

I’ve posted a few modeling challenges on frame analysis (strongback, Ziemian, and stability) and soil-structure interaction.

However, I recently accepted a challenge from George Chamosfakidis to see if OpenSees can give the same periods and mode shapes reported in the ETABS verification example shown below.

Published verification examples typically just show the “correct” result and not all the ways you can get an “incorrect” result. This post will show the “correct”–and several “incorrect”–approaches for this example using OpenSees.

The model is shown in the figure above, so I won’t redraw it; however, here are the modeling assumptions.

• Neglect shear and axial deformation in the columns
• Beams are infinitely rigid
• Columns have rigid joint offsets equal to the beam depth (36 inch)
• Columns are completely fixed at their base

With these assumptions, the modulus of elasticity, and the column sizes provided, the distribution of stiffness is pretty straightforward. In the code below, I multiply by 1000, where appropriate, to make elements “rigid”. For all members, I use nonlinear elements with elastic sections and two-point Gauss integration.

inch = 1.0
kip = 1.0
sec = 1.0

g = 386.4*inch/sec**2

ft = 12*inch
lb = kip/1000
ksi = kip/inch**2
psf = lb/ft**2

E = 3000*ksi
v = 0.3
G = 0.5*E/(1+v)
q = 150*psf

L = 30*ft
H = 15*ft

import openseespy.opensees as ops

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

ops.node(1,0,0,0); ops.fix(1,1,1,1,1,1,1)
ops.node(2,L,0,0); ops.fix(2,1,1,1,1,1,1)
ops.node(3,0,L,0); ops.fix(3,1,1,1,1,1,1)
ops.node(4,L,L,0); ops.fix(4,1,1,1,1,1,1)

ops.node(11,0,0,H)
ops.node(12,L,0,H)
ops.node(13,0,L,H)
ops.node(14,L,L,H)

# Beam depth
d = 36*inch

# Column C1, C2
b = 24*inch
A = b*b
I = b*b**3/12
J = 2*I

ops.section('Elastic',1,E,1000*A,I,I,G,J)
ops.beamIntegration('Legendre',1,1,2)

ops.geomTransf('Linear',1,1,0,0,'-jntOffset',0,0,0,0,0,-d)

ops.element('forceBeamColumn',1,1,11,1,1)
ops.element('forceBeamColumn',2,2,12,1,1)

# Column C3, C4
b = 18*inch
A = b*b
I = b*b**3/12
J = 2*I

ops.section('Elastic',2,E,1000*A,I,I,G,J)
ops.beamIntegration('Legendre',2,2,2)

ops.element('forceBeamColumn',3,3,13,1,2)
ops.element('forceBeamColumn',4,4,14,1,2)

# Beams
b = 18*inch # assumed
A = b*d
Iz = b*d**3/12
Iy = d*b**3/12
J = Iz+Iy

ops.section('Elastic',3,1000*E,A,Iz,Iy,1000*G,J)
ops.beamIntegration('Legendre',3,3,2)

# B1,B2
ops.geomTransf('Linear',3,0,1,0)
ops.element('forceBeamColumn',5,11,12,3,3)
ops.element('forceBeamColumn',6,13,14,3,3)

# B3,B4
ops.geomTransf('Linear',4,1,0,0)
ops.element('forceBeamColumn',7,11,13,4,3)
ops.element('forceBeamColumn',8,12,14,4,3)

The ETABS manual states “the example is a three degree-of-freedom system”, which means we should define a rigid diaphragm for the floor. Let’s define the primary node at the geometric center of the floor.

ops.node(10,L/2,L/2,H)
ops.fix(10,0,0,1,1,1,0)
ops.rigidDiaphragm(3,10,11,12,13,14)

The resulting OpenSees model, rendered with opsvis, is shown below for reference.

The distribution of mass is not so straightforward and there’s a few ways to go off the rails with the specified floor weight.

I’ll first show the “correct” distribution of mass–correct in the sense that the natural periods match the theoretical solution reported in the ETABS manual. Then I’ll show “incorrect” distributions of mass, along with a few other incorrect modeling assumptions George and I made on our way to the correct solution. It’s like reporting a crime scene investigation.

The correct approach is to assign the total floor mass to the horizontal DOFs of the primary rigid diaphragm node, along with rotational mass about the vertical DOF of this node.

W = q*L*L     # Total floor weight
m = W/g       # Total floor mass
mz = m*L**2/6 # Rotational mass about vertical axis
ops.mass(10,m,m,0,0,0,mz)

With this distribution of mass, the natural periods are summarized below.

Mode	Theoretical (sec)	ETABS (sec)	OpenSees (sec)
1	0.1389	0.1389	0.1385
2	0.1254	0.1254	0.1254
3	0.070	0.0702	0.06966

I’m not sure why the OpenSees model is a wee bit stiffer in torsion (modes 1 and 3) compared to the ETABS results. I am sure that I will not lose any sleep over it though.

Update January 2022 – this post shows how to better match ETABS.

Below are the mode shapes from OpenSees. Be sure to use the 'fullGenLapack' option for the eigen command so that you can get all three modes.

There is torsion in the first mode due to the different column sizes. The second mode is symmetric flexure and the third mode is torsion. You can verify that the eigenvectors from OpenSees are in direct proportion to those reported in the ETABS example.

Now let’s take a look at some incorrect modeling approaches. I encourage you to plot the mode shapes for each case to see the effect of modeling assumptions.

First, what happens if you keep the rotational mass, but assign mass to the top of each column based on tributary area?

#ops.mass(10,m,m,0,0,0,mz)
ops.mass(10,0,0,0,0,0,mz)
ops.mass(11,m/4,m/4,0,0,0,0)
ops.mass(12,m/4,m/4,0,0,0,0)
ops.mass(13,m/4,m/4,0,0,0,0)
ops.mass(14,m/4,m/4,0,0,0,0)

In this case the model becomes more flexible and you get periods 0.1701, 0.1254, and 0.1135 sec. Moving the masses away from the center of mass increases the rotational inertia of the floor.

And what happens if you remove the rotational mass, but keep the mass at the top of each column?

#ops.mass(10,m,m,0,0,0,mz)
#ops.mass(10,0,0,0,0,0,mz)
ops.mass(11,m/4,m/4,0,0,0,0)
ops.mass(12,m/4,m/4,0,0,0,0)
ops.mass(13,m/4,m/4,0,0,0,0)
ops.mass(14,m/4,m/4,0,0,0,0)

The model is still too flexible. The periods are 0.1564, 0.1254, and 0.1069 sec.

On a side node, let’s return to the correct mass distribution then see what happens with a DIY approach to the rigid diaphragm. Instead of using the rigidDiaphragm command, the DIY approach uses equalDOF to constrain the column nodes to have the same horizontal translations and rotation about the vertical axis as the primary node.

ops.fix(10,0,0,1,1,1,0); 
#ops.rigidDiaphragm(3,10,11,12,13,14)
ops.equalDOF(10,11,1,2,6)
ops.equalDOF(10,12,1,2,6)
ops.equalDOF(10,13,1,2,6)
ops.equalDOF(10,14,1,2,6)

Like my DIY projects around the house, something major goes awry. In this case, torsion essentially disappears because the equalDOF constraints prevent the floor from rotating. The periods are 0.1254, 0.1254, and 0.006524 sec.

Let’s return to the correct mass distribution and rigid diaphragm and look at a few modeling assumptions that affect the stiffness rather than the mass.

First, what if the beam stiffnesses are based on the nominal material properties, i.e., flexible, not “rigid”?

#ops.section('Elastic',3,1000*E,A,Iz,Iy,1000*G,J)
ops.section('Elastic',3,E,A,Iz,Iy,G,J)

In this case, the periods become 0.1743, 0.1694, and 0.09140 sec. A beam width was not stated in the example, so I assumed 18 inch.

Back to rigid beams, what happens if the columns are not axially “rigid”?

#ops.section('Elastic',1,E,1000*A,I,I,G,J)
#ops.section('Elastic',2,E,1000*A,I,I,G,J)
ops.section('Elastic',1,E,A,I,I,G,J)
ops.section('Elastic',2,E,A,I,I,G,J)

Columns are already stiff for axial response, so the natural periods change slightly to 0.1390, 0.1260, and 0.06973 sec.

What happens if we get overzealous and make the torsional response of the columns “rigid” as well?

#ops.section('Elastic',1,E,1000*A,I,I,G,J)
#ops.section('Elastic',2,E,1000*A,I,I,G,J)
ops.section('Elastic',1,E,1000*A,I,I,1000*G,J)
ops.section('Elastic',2,E,1000*A,I,I,1000*G,J)

We reduced torsional flexibility and, as expected, the periods become 0.1258, 0.1254, and 0.01560 sec.

Finally, what if we had everything else correct, but omitted the rigid joint offset at the top of each column?

#ops.geomTransf('Linear',1,1,0,0,'-jntOffset',0,0,0,0,0,-d)
ops.geomTransf('Linear',1,1,0,0)

The columns become more flexible and the natural periods increase to 0.1933, 0.1752, and 0.09684 sec.

This post covered three approaches to mass distribution and five approaches to element stiffness for a very simple frame. That’s 15 possible paths to follow, only one of which is “correct”. The number of paths increases rapidly if we include the DIY rigid diaphragm approach, other modeling assumptions, and the myriad other potential errors, e.g., due to units, element orientation, and boundary conditions.