Find any indeterminate beam, frame, or truss problem from a structural analysis textbook, and you can make OpenSees solve it. But sometimes, replicating the basics is not so easy. Take, for instance, an axial-moment (P-M) interaction diagram of reinforced concrete (RC) sections.

The typical approach advocated with OpenSees is to use repeated moment-curvature analyses over a range of axial loads. For each axial load applied on a zero length section, use displacement control with a reference moment–or proportional reference moments for biaxial bending–until reaching a target fiber strain, section curvature, or some other limit state.

However, this approach doesn’t match the textbook procedure for analyzing RC sections:

- Assume the extreme concrete compression fiber is at 0.003 strain
- Obtain a linear strain profile by assuming a neutral axis depth or a reinforcing steel strain
- Compute stresses in the Whitney stress block (WSB) and the steel bars
- Compute stress resultants from the WSB and steel bars
- Determine the equivalent force-couple system acting at the section centroid

Repeat these steps over a range of neutral axis depths (or reinforcing steel strains) and you will generate a P-M interaction diagram. These steps are easily coded in MATLAB, Python, and even Excel.

The **first point of this post** is that we can make OpenSees do exactly this textbook WSB procedure. As usual, define your concrete and steel materials, put them in a fiber section, and then put that section in a zero length section element. Here we use *Concrete01* that reaches its peak compressive strength, , at 0.002 strain and spalls at 0.006 strain. The steel is EPP.

```
ops.node(1,0.0,0.0)
ops.node(2,0.0,0.0)
ops.fix(1,1,1,1)
ops.fix(2,0,1,0)
fc = 4.0 # ksi
fy = 60.0 # ksi
E = 29000.0 # ksi
ops.uniaxialMaterial('Concrete01',concreteTag,fc,0.002,0,0.006)
ops.uniaxialMaterial('Steel01',steelTag,fy,E,0.0)
#
# Define your section (secTag)
#
ops.element('zeroLengthSection',1,1,2,secTag)
```

Then, for any strain in the “bottom” reinforcing steel, compute the strain profile. For reference, a rectangular RC section is shown below along with a strain profile.

Using the strain profile, we can use `sp`

(single point) constraints to impose the axial deformation and curvature at the section centroid.

```
# e.g., at the balance point
epss = fy/E # Tension steel yield strain
kappa = (epss+epsc)/d
c = epsc/kappa # Depth to NA
ops.timeSeries('Constant',2)
ops.pattern('Plain',2,2)
ops.sp(2,1,-epsc + kappa*h/2) # Strain at the centroid
ops.sp(2,3,kappa) # Curvature
```

The axial load and bending moment are then the “reactions” on the node with the imposed `sp`

constraints. Be sure to use `'Transformation'`

constraint handler and also use the `'UmfPack'`

solver as it’s the only solver that won’t rain on your parade for having all DOFs constrained.

```
ops.constraints('Transformation')
ops.system('UmfPack')
ops.analysis('Static')
ops.analyze(1)
ops.reactions()
P = -ops.reaction(2,1)
M = ops.reaction(2,3)
```

Put all of this in a loop over the steel strain, `epss`

, from -0.003 (pure compression) to something large (say five times the steel yield strain) and you’ll get the P-M interaction curve shown below along with a straight-up Python implementation of the textbook WSB procedure.

I’m not too worried about the OpenSees interaction curve taking a circuitous route to pure compression. The top of the interaction curve will get chopped off at 85% anyway due to “accidental moments”, the rumored sequel to a 2007 Judd Apatow movie. If you’re curious, the roundabout path is due to the descending branch of *Concrete01* from 0.002 to 0.003 compressive strain. Use *Concrete23* and you’ll see a different detour to pure compression.

However, there’s still a difference between the WSB approach and OpenSees for axial loads below 85% of the pure compression load.

So, the **second point of this post** is that we can make OpenSees match the textbook WSB approach. To this end, we need a concrete stress-strain model that has compressive stress equal to for compressive strains exceeding . This strain is based on similar triangles with the depth, , of the WSB, where is the WSB parameter which depends on .

Nothing fits the bill, at least not directly. I can think of a couple kludgy workarounds, the least complicated of which is the *ElasticMultiLinear* material where we can define points that characterize the stress block.

```
fc = 4.0 # ksi
epsc = 0.003
beta1 = 0.85 # for 4 ksi concrete, smaller for higher strength concrete
strain = [-epsc,-epsc*(1-beta1),-epsc*(1-beta1),0]
stress = [-0.85*fc,-0.85*fc,0,0]
ops.uniaxialMaterial('ElasticMultiLinear',concreteTag,'-strain',*strain,'-stress',*stress)
```

Sure, with these inputs, the tangent of this material model is either zero of infinity, so I don’t recommend you use this model for general FEA. But this model is OK for our zero length section analysis because we’re not solving any nodal equilibrium equations–just imposing displacements and calculating forces.

As shown below, faking the WSB with *ElasticMultiLinear* instead of using *Concrete01* gives a P-M interaction curve that agrees much better with the straight-up Python implementation of the WSB.

Due to the concrete fiber discretization, the OpenSees solution wiggles around the balance point (I used 20 fibers). These wiggles go away as you use more concrete fibers. Try it out yourself!

You read this far, so let’s see if you paid attention. How many unique factors of 0.85 does this post contain? Explain your response and the different 0.85 factors in the comments section below.

I would say three unique 0.85 factors used.

1. 85% of pure axial load to account for accidental moments.

2. 85% of compressive strength to estimate the equivalent stress block parameters

3. 85% of neutral axis depth, again to estimate the equivalent stress block parameters

LikeLiked by 2 people