Verifying material nonlinear frame element formulations is pretty difficult. There are only a handful of analytical solutions that do not rely on bilinear uniaxial constitutive response. You’re much better off verifying geometric nonlinearity.
However, an example based on a relatively straightforward biaxial stress resultant plasticity model has intrigued me over the years. The example, described in more detail here, is of the two-span continuous beam shown below.
The transverse load is a sinusoidal function, , while the axial load increases linearly with time, .
The ensuing interaction of axial force and bending moment along the beam is described by a stress resultant plasticity model with an elliptical yield surface, described here. The model parameters include elastic stiffness (k1 and k2) and yield force (Fy1 and Fy2) for each stress resultant as well as an isotropic hardening parameter, .
When loaded along direction 1, the force-deformation response is bilinear with isotropic hardening modulus . If you load along the other axis, the isotropic hardening modulus is . Under combined loading, the isotropic hardening should fall somewhere in between these two values. But those details are not important.
This section model, including an extension for kinematic hardening, is implemented in OpenSees as
Elliptical. Along with its less complicated, older, and documented sibling
Bidirectional, this model can couple any two stress resultants in a single section, e.g., moment-shear as described here. The input commands allow you to specify any two modes of section response to be coupled, i.e., P-M is not hard wired. Side note: These models are also appropriate for base isolators in zero length section elements.
Returning to the beam example at hand, let’s use the elliptical stress resultant plasticity model with axial response k1=EA=3000 and Fy1=Ny=20, flexural response k2=EI=3000 and Fy2=My=40, and isotropic hardening parameter =0.2. The beam is discretized in to four force-based elements, each with 10 Gauss-Lobatto points. The excessive number of integration points gives better resolution on the distribution of section deformations.
The global load-displacement response for t=[0, 5.2] is shown below.
And the distributions of curvature and axial deformation along the beam at the end of the analysis (t=5.2), are shown below.
Neither the global nor local response matches the results described in the paper, but they capture the essence. So, nothing’s really been verified. And I’m OK with that. It just feels good to put these OpenSees results out there.