Using nonlinear elements, particularly the forceBeamColumn element, with elastic sections is just as good as, if not better than, using the elasticBeamColumn element for many reasons. Not only do force-based elements with elastic sections make the transition to material nonlinearity easy, they also facilitate debugging your model.
Another reason I like force-based elements is you can record the section response at the integration points. But, how many integration points are required to get the correct linear-elastic solution? I will focus on Gauss-Lobatto integration, but there are many other options.
For linear-elastic, prismatic frame elements without member loads, you have to integrate the quadratic polynomials that arise from the produce of linear curvature with linear weighting functions. Gauss-Lobatto quadrature integrates simple polynomials up to order . Therefore, to integrate
, you need
, or
. (I know, the algebraic solution is
, but you can’t have half an integration point).
Consider the indeterminate beam shown below with L=20 ft, E=29,000 ksi and I=800 in4. The applied moment is M=400 kip-ft.

And here is a snippet of code to create the element.
E = 29000
A = 10 # Not important
I = 800
ops.section('Elastic',8,E,A,I)
ops.beamIntegration('Lobatto',12,8,3)
ops.geomTransf('Linear',5)
ops.element('forceBeamColumn',1,1,2,5,12)
The exact bending moment distribution is linear and the end rotation is . Using one force-based element with three Gauss-Lobatto points, we get the following bending moment and end rotation, 0.01241379 rad, which matches the exact solution.

With force-based elements, you can include member loads in the element equilibrium relationship. For the common case of a uniform distributed load, the curvature distribution will be quadratic while the weighting function is linear, so you need to integrate . Fortunately, three Gauss-Lobatto points still get the job done.
Now consider the same beam but with a uniform distributed load w=1 kip/ft.

The exact bending moment distribution is quadratic and the end rotation is . Using three Gauss-Lobatto points in a single force-based element, we obtain the following distribution of internal moment. The computed end rotation, 0.00103448, matches the exact solution.

For both examples shown in this post, you can use more Gauss-Lobatto points and get the same, exact, results.
Note that if your member loads are linearly varying, e.g. hydrostatic pressure on a vertical member, you will need to add another integration point to get the elastic solution.
Also, member loads that cause discontinuities along an element, e.g., point loads or partial distributed loads, will give an integration error no matter how many integration points you use in a linear-elastic, prismatic element. You will also always get an error for non-prismatic elements no matter what integration you throw at it.
Just remember that with material nonlinearity you will want to use more integration points to capture the spread of inelasticity, or a regularization technique to simulate localized response.

Nice blog.
I have one confusion. As Force Based BeamColumn assumes the linear distribution of Forces (Moments). How does it get the exact result in case of the uniformly distributed beam, where bending moment distribution is non-linear.
Regards,
Zeeshan
LikeLike
Hello Zeeshan,
Have a look at this paper: https://doi.org/10.1061/(ASCE)0733-9445(1997)123:7(958)
PD
LikeLiked by 1 person
Dear Prof. Scott,
Thank you for this insightful post. While reproducing your results, I noticed that the Plain constraint handler introduces a significant displacement error (approx. 20%). Curiously, this error appears sensitive to the number of analysis steps (N) in LoadControl.
Switching to Transformation constraints immediately yields the exact solution. It seems the Plain handler struggles to enforce boundary equilibrium for a single force-based element when eleLoad is present.
I’ve attached a code link below illustrating this “step-dependency” and the discrepancy between handlers. I would greatly appreciate your insights on why the constraint choice and step size impact a linear elastic element so significantly.
https://drive.google.com/file/d/1DurWjM7JJCS7-_NrKx1DMlpJ8Yr1sdIM/view?usp=drive_link
LikeLiked by 1 person
Hello Fei,
Thanks for providing the example! The force-based element handles member loads without adding nodal loads, so the global solution will see zero change in nodal displacement and converge. For your script with ‘test NormDispIncr’, the norm is zero on the first iteration, but the norm of the residual is >> 0. Should work if you change the test to ‘NormUnbalance’. It’s not a constraint handler issue!
Michael
LikeLike
Thank you, Prof. Scott! That makes perfect sense. Switching to NormUnbalance fixed it immediately.
LikeLike