I talked with a graduate student from Mechanical Engineering the other day. The student is learning OpenSees and successfully analyzed a truss. No, not that truss from Example 1.1.
After showing me the truss results, the student said something along the lines of “Deflections are fine and OpenSees does a good job, but I really want to know how to get stresses from an analysis. Like, can I analyze a cantilevered L-shaped rod with loads at the free end and get the stresses?” I knew exactly what the student was talking about and agreed that getting stresses out of OpenSees is not always straightforward.
I recall a similar problem from when I took strength of materials. The problem was Example 8-4 from Hibbeler’s 2nd edition Mechanics of Materials.
At point B, on the perimeter of the section, Hibbeler shows the normal stress is 21.4 ksi (due to axial force and bending moment) and the shear stress is 17.5 ksi (due to shear force and torsion).
This analysis should be doable in OpenSees.
First, the rod is statically determinate, so the bend radius doesn’t really matter. As a result, we only need four nodes: at the support, at section A-A, at the bend, and at the free end.
import openseespy.opensees as ops
inch = 1.0
kip = 1.0
lb = kip/1000
ksi = kip/inch**2
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,0,8*inch,0)
ops.node(3,0,18*inch,0)
ops.node(4,14*inch,18*inch,0)
The next step is the section definition. The NDFiber section is the way to go because we want to track axial and shear stresses. We’re not worried about deflections, so the material is not important, as long as the response remains elastic.
E = 29000*ksi
v = 0.3
ops.nDMaterial('ElasticIsotropic',1,E,v)
r = 0.75*inch
# Wedges and rings for the section discretization
Nwedge = 24
Nring = 10
dtheta = 0.5*360/Nwedge
ops.section('NDFiber',1)
ops.patch('circ',1,Nwedge,Nring,0,0,0,r,dtheta,360+dtheta)
The start and finish angles of dtheta and 360+dtheta ensure fiber centroids fall on the section y and z axes, as shown below. A fiber will be close to B, but no matter how many rings we use, no fiber will ever be exactly at B (more on that later).
For a round section, the circular patch discretization in OpenSees is not optimal because many fibers are crowded near the center of the section where not much happens in terms of flexural and torsional stresses.
Sure, you can define multiple patches as rings and DIY a mesh gradient from low density fibers in the center to higher density closer to the perimeter. But that’s a lot of work.
Instead, triangulation, e.g., with the sectionproperties library, gives a more uniform distribution of fibers over the cross-section.
from sectionproperties.pre.library import circular_section
from sectionproperties.post.fibre import to_fibre_section
A = 3.14159*r**2 # Area
geom = circular_section(d=2*r,n=40) # Diameter and number of perimeter points
geom.create_mesh(mesh_sizes=A/200) # Maximum element size
fiberfile = 'solidrod.txt'
to_fibre_section(geom, analysis_type="3D", save_to=fiberfile)
To read the exported fiber data into an OpenSees fiber section, see this post.
But again, no matter what we do with triangulation, we’ll never get a fiber centroid exactly at point B.
However you define the section, the next step is to define the elements. We need one of the element integration points to line up with the location of section A-A. I’ll use Gauss-Lobatto integration in the element between the support and section A-A so that an integration point falls right at the section of interest. There’s other ways to finagle the element definition to get an integration point at A-A.
Two-point Gauss-Legendre integration is fine for the other two elements.
ops.beamIntegration('Lobatto',1,1,3) # element 1
ops.beamIntegration('Legendre',2,1,2) # elements 2 and 3
ops.geomTransf('Linear',1,1,0,0) # elements 1 and 2
ops.geomTransf('Linear',2,0,-1,0) # element 3
ops.element('forceBeamColumn',1,1,2,1,1)
ops.element('forceBeamColumn',2,2,3,1,2)
ops.element('forceBeamColumn',3,3,4,2,2)
The vector in the x-z plane for elements 1 and 2 coincides with the global X-axis while the x-z vector points along the negative Y-axis for element 3. We could also use [1,-1,0] for the vector in the x-z plane of all three elements and let the cross products sort out the actual local y and z axes. Try it yourself.
With the model defined, we can now define the load at the free end and perform a static analysis.
ops.timeSeries('Constant',1)
ops.pattern('Plain',1,1)
ops.load(4,0,500*lb,-800*lb,0,0,0)
ops.analysis('Static','-noWarnings')
ops.analyze(1)
To record stress near, but not exactly at, point B of section A-A, we can use the eleResponse command with coordinates (y,z)=(0,r) for integration point (section) 3 of element 1.
sigma = ops.eleResponse(1,'section',3,'fiber',0,r,'stress')
When using the NDFiber section in three-dimensional models, three stresses are returned for each fiber: normal stress and two shear stresses.
For the section defined by a circular patch of fibers, the stresses at the fiber closest to B are:
- Normal stress,
= 20.5 ksi
- Shear stress,
= 15.7 ksi
- Shear stress,
= 0
As expected, the normal stress is lower than the exact value of 21.7 ksi because the fiber is not quite at the perimeter of the section. Also, the shear stress does not match the exact value of 17.5 ksi. Although most of the shear stress comes from torsion, the contribution from shear force is underestimated because the NDFiber section uses an average shear stress distribution through the section depth instead of the exact parabolic distribution. Using more fibers won’t get us any closer to the exact solution, but that’s an LBU for another day.
And for the section defined by triangulated fibers:
- Normal stress,
- Shear stress,
- Shear stress,
The shear stress
In both cases of section discretization, we can get a lot closer to the exact normal stress of 21.7 ksi by defining a fiber with zero area directly at point B. But that’s another LBU for yet another day. And still, we’ll always be slightly off from the exact solution due to the fiber discretization’s inherent error in approximating the second moments of section area.
Portwood Digital 