Do It Your Self-Weight

Original Post - 05 Nov 2023 - Michael H. Scott

Visit Structural Analysis Is Simple on Substack.

Most solid elements in OpenSees don’t do body forces very well, if at all. Among elements that have body forces, the implementation and execution are inconsistent. With surface loads, the story is even more convoluted.

However, OpenSees is pretty good at dynamics, so all solid elements handle mass density. Yes, the implementations are inconsistent because some elements have mass density as an element property while other elements pull the mass density from their materials. But, either way, you can count on solid elements in OpenSees providing an element mass matrix.

Then, logically, if elements can provide mass, you should be able to apply self-weight as mass times gravitational acceleration, g, as nodal loads, external to the elements.

With PRs #1265 and #1266, I made applying self-weight to solid elements an easy DIY project. The PRs implemented the printA function, previously restricted to FullGeneral, for the Diagonal system of equations.

Combined with the GimmeMCK integrator and some easy scripting, you can use printA from the Diagonal system to grab the mass at each node, then apply self-weight as a nodal load in a pattern with constant time series. Or, you can ramp up the self-weight in a linear time series and then use the loadConst command.

A simple test case is a 50x50x20 box of linear-elastic tetrahedral elements (plotted with opsvis below). The mesh size is 10 and the elements have mass density \(\rho\). The specific weight is then \(\rho g\) and the total weight of the box is 50,000\(\rho g\). The numbers aren’t all that important for this demonstration.

Box of tetrahedral elements

First, use GimmeMCK to get the mass, which will be a one-dimensional array of mass values. Be sure to use the lumped option on the Diagonal system, otherwise you won’t get all the mass if the elements return a consistent mass matrix. Then, loop over the nodes and apply mass times gravity in the direction of gravity load.

#
# Define your model
#

NDF = 3 # DOFs per node
gravDOF = 2 # 0 = X, 1 = Y, 2 = Z

ops.integrator('GimmeMCK',1,0,0)
ops.system('Diagonal','lumped')
ops.analysis('Transient')
ops.analyze(1,0.0)

M = ops.printA('-ret')

ops.timeSeries('Constant',1)
ops.pattern('Plain',1,1)
for nd in ops.getNodeTags():
    
    dofs = ops.nodeDOFs(nd)
    gravload = [0]*NDF

    # Skip constrained DOFs
    if dofs[gravDOF] < 0:
        continue
    else:
        gravload[gravDOF] = -g*M[dofs[gravDOF]]

    ops.load(nd,*gravload)
    
ops.fixZ(0,1,1,1)

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

Be sure to impose boundary conditions after applying the self-weight; otherwise, you will under count the total weight of your model because the loop ignores constrained DOFs which have equation number -1. The check on negative equation numbers is not entirely necessary if you haven’t yet imposed boundary conditions, but the check prevents you from getting unexpected results if you have imposed boundary conditions. Python uses index -1 as the last list entry.

Finally, perform static analysis for the gravity load. Be sure to not use the Diagonal system for the static analysis. You’ll get results–very incorrect results!

To verify this DIY self-weight approach is correct, sum the vertical reactions and compare to the total weight of the box.

ops.reactions()

Ftotal = 0
for nd in ops.getNodeTags():
    Ftotal += ops.nodeReaction(nd,gravDOF)    
print(Ftotal)

The printA implementation for the Diagonal system is available in OpenSeesPy version 3.5.1.8 and later as well as OpenSees.exe version 3.6.0 and on.

The self-weight approach shown in this post will also work for beam-column and shell elements. For some cases with the beam-column elements, this approach may be easier and more general than using the eleLoad command.

Also, beware that when you have massless rotational DOFs in shell and frame models, you’ll get a solver failed message when forming the mass matrix via the Diagonal/GimmeMCK combo. You can ignore that message.

OpenSees Cloud

OpenSees AMI

Do It Your Self-Weight