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 `ops_vis`

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

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.

Mode 1 Mode 2 Mode 3

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.

Dear PD,

thank you for the very interesting post.

I thought it might be interesting to share that using the rigid diaphragm also the following distributions of masses give the same periods:

ops.mass(10,0,0,0,0,0,mz)

mz_edge = – m/4*( (L/2)**2+(L/2)**2 )

ops.mass(11,m/4,m/4,0,0,0,mz_edge)

ops.mass(12,m/4,m/4,0,0,0,mz_edge)

ops.mass(13,m/4,m/4,0,0,0,mz_edge)

ops.mass(14,m/4,m/4,0,0,0,mz_edge)

as well as:

mz_edge = mz/4 – m/4*( (L/2)**2+(L/2)**2 )

ops.mass(11,m/4,m/4,0,0,0,mz_edge)

ops.mass(12,m/4,m/4,0,0,0,mz_edge)

ops.mass(13,m/4,m/4,0,0,0,mz_edge)

ops.mass(14,m/4,m/4,0,0,0,mz_edge) #no rotational mass in node 10

Where the term “- m/4*( (L/2)**2+(L/2)**2 )” is mass times square of the distance, with a negative sign.

I am puzzled by what the meaning of negative rotational mass of a single node might be. But in global terms, it is reducing the rotational inertia of the system.

Best

LikeLiked by 1 person

That’s why FEM isn’t easy. Different people produce different result.

LikeLiked by 1 person

what about using joint offset of column d/2, as the column and beam are modeled as the centreline.

LikeLike

Do you mean offsets at each end of the beams equal to column d/2? I don’t think that will matter because the beams are “rigid”.

LikeLike

No, I mean column joint offset at top -d/2. I have done that in most of the cases which makes structure flexible. You have considered -d

LikeLike

One beam depth is what’s stated in the ETABS verification example.

LikeLike