A Solution, Just Not The Solution

Force-based elements satisfy equilibrium in strong form, even with member loads. However, this does not mean force-based elements always get the exact solution.

Consider a simple prismatic, linear-elastic beam with a point load at mid-span.

Using a single force-based element with a single point load applied to the element using the eleLoad command.

E = 1; I = 1; A = 1
ops.section('Elastic',1,E,A,I)
Np = 5
ops.beamIntegration('Lobatto',1,1,Np)
ops.geomTransf('Linear',1)
ops.element('forceBeamColumn',1,1,2,1,1)

P = 40
ops.timeSeries('Constant',1)
ops.pattern('Plain',1,1)
ops.eleLoad('-ele',1,'-type','beamPoint',-P,0.5)

With this model, we do obtain the exact equilibrium solution for the vertical reactions (P/2) and the internal bending moment.

However, because the numerical integration cannot handle the discontinuity caused by the point load, the end rotation is not exact. The computed end rotation is 17394.478 while the exact end rotation is =16000–a relative error of about 9%. The curvature, virtual moment, and product of these two functions–which is what’s integrated numerically in the force-based element state determination–are shown below.

That change in at mid-span causes problems for numerical integration based on smooth polynomials with continuous derivatives. However, as we use more integration points, the error reduces, e.g., if we used 9 Lobatto points instead of 5, the relative error for the end rotation reduces to about 2%. But the error will never go to zero.

Not really a big deal in my opinion. Besides, the model is statically determinate, so we don’t need no stinkin’ compatibility equations to get the equilibrium solution anyway.

But what happens when the model is indeterminate, e.g., by fixing the rotation at the left support?

Because of the error in compatibility, which we now need in order to find equilibrium, the internal moments are slightly off from the exact solution.

Similarly, the vertical reactions are slightly off from the exact solutions of 11P/16 and 5P/16, and the end rotation is 8697.239 instead of the exact =8000.

Among the infinite equilibrium solutions for this indeterminate model, we found one that was close, but not equal, to the unique solution that also satisfies compatibility. You could fix this “problem” by using composite integration, e.g., composite Simpson integration, on each side of the point load.

But instead of going to such great lengths with integration points and weights, you could define two elements, make the point load a nodal load, and save yourself some trouble. And what if the point load moves during the analysis? You’re not going to re-define elements or integration points on the fly. Just accept the error and move along. And with material nonlinear response, this integration error will be the least of your concerns.