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)|
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”?
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?
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.
12 thoughts on “Verifying Ain’t Easy”
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:
mz_edge = – m/4*( (L/2)**2+(L/2)**2 )
as well as:
mz_edge = mz/4 – m/4*( (L/2)**2+(L/2)**2 )
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.
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.
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”.
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
One beam depth is what’s stated in the ETABS verification example.
thank you so much for cheking result between opensees and etabs, it can be very helpful for verification. i’m also has to do verification one story building with just 4 elastic columns and no beam. the model has irregulaty in mass distribution because of it calculating of rotational mass is more complicated than commuonmass distribution. i can verify my model in etabs and achive my modal frequences in resource with 0.33% error in first mode. in opensees community someone offer me your cheking. i built my model as you said but it was unssucsesfull in verification.
i do some research about how etabs evaluate rotational mass and found that etabs calculate rotational mass for each mass node by multipleing mass by node’s distance from center of mass (like dear prafulla said). so i tried build model with this pattern also i use fiber section and displacement element then i achive first mode frequences with 1% error and i think it’s acceptable. also i use 5 point for integration and rigid diaphragm. i don’t think other things you considerd may efiicint like offset, type of integration or the way you evaluate rotational mass.
I look forward to discussing this further if you intrested and also i can share my files with you.
with best regard
Thank you, Arsham, for the information about how ETABS calculates the rotational mass.