U1 beam element: nodal forces have wrong signs (axial at node 1, shear at node 2), so RF does not balance the loads

Analysis issues
3

And another related issue:

Summary

U1 is documented as a 3D Timoshenko beam (manual §6.2.46). The code cites Luo (2008) for its consistent shape functions. With such shape functions, one element should give the exact Timoshenko deflection of a cantilever under an end load.

Instead, the shear parameter φ = 12EI/(kGAL²) enters U1’s shape functions and stiffness with the wrong sign:

one-element tip deflection
ccx 2.23 PL³/(3EI) · (φ² + φ + 4) / (4(1 + φ))
Timoshenko PL³/(3EI) · (1 + φ/4)

So for φ < 3 the element is stiffer than Euler–Bernoulli, and at φ = 1 its stiffness matrix is singular.

Already reported or fixed?

Not reported: there is no matching GitHub issue (open or closed) or open pull request in Dhondtguido/CalculiX, and no calculix.discourse.group topic (searched 17 Oct 2026). Not fixed: the routine involved is unchanged on master, and conda-forge’s latest 2.23 build (pl5321hd17ca8b_8, which carries about 40 extra patches) gives the same output. Related: forum topic 600 (2021–2023) found another sign error taken from the same paper’s Appendix A, fixed in 2.21.

Minimal input (u1_shear.inp)

** Deep cantilever, one U1 element: L = 5, square section 1 x 1
** (A = 1, I = 1/12), E = 1000, nu = 0.25 (G = 400), k = 5/6,
** tip load P = 1 along n1.
** Timoshenko: v = P L^3/(3 E I) + P L/(k G A) = 0.5 + 0.015 = 0.515
*NODE, NSET=NALL
1, 0., 0., 0.
2, 5., 0., 0.
*USER ELEMENT, TYPE=U1, INTEGRATION POINTS=2, MAXDOF=6, NODES=2
*ELEMENT, TYPE=U1, ELSET=EBEAM
1, 1, 2
*MATERIAL, NAME=MAT
*ELASTIC
1000., 0.25
*BEAM SECTION, ELSET=EBEAM, MATERIAL=MAT, SECTION=GENERAL
1., 0.083333333, 0., 0.083333333, 0.83333333
0., 1., 0.
*BOUNDARY
1, 1, 6
*STEP
*STATIC
*CLOAD
2, 2, 1.
*NODE PRINT, NSET=NALL
U
*EL PRINT, ELSET=EBEAM
S
*END STEP

How to run

ccx u1_shear

Expected

Timoshenko beam theory (φ = 0.12):

  • tip deflection v = PL³/(3EI) + PL/(kGA) = 0.5 + 0.015 = 0.515, compared with Euler–Bernoulli 0.5;
  • tip rotation θz = PL²/(2EI) = 0.15;
  • section forces: shear force sxy = k·G·A·(v′ − θz) = +1 along the whole element; bending moment magnitude 5 at the root and 0 at the tip.

Actual

 displacements (v(i),i=1..ndof) for set NALL and time  0.1000000E+01

         1  0.000000E+00  0.000000E+00  0.000000E+00  0.000000E+00  0.000000E+00  0.000000E+00
         2  0.000000E+00  4.614286E-01  0.000000E+00  0.000000E+00  0.000000E+00  1.500000E-01

 stresses (elem, integ.pnt.,sxx,syy,szz,sxy,sxz,syz) for set EBEAM and time  0.1000000E+01

         1   1  0.000000E+00 -4.464286E+00  0.000000E+00 -7.857143E-01  0.000000E+00  0.000000E+00
         1   2  0.000000E+00 -5.357143E-01  0.000000E+00 -7.857143E-01  0.000000E+00  0.000000E+00
  • v = 0.4614, which is 10.4 % below Timoshenko and below the Euler–Bernoulli value 0.5.
  • The moments are 4.46 at the root and 0.54 at the tip, instead of 5 and 0.
  • The shear force is −0.786 instead of +1.
  • The rotation is right.

Refining the mesh converges only slowly to 0.515:

elements 1 2 4 8 16 128
tip v 0.4614286 0.4744595 0.4944521 0.5080876 0.5131085 0.5149695

The test suite shows the same effect. test/userbeam.dat.ref has a tip deflection of 1.022509E-01. Its own data give 1.024003E-01 for Euler–Bernoulli and 1.024502E-01 for Timoshenko (φ = 0.00195).

With 12EI = kGAL² (φ = 1), for example A = 3, I1 = 0.25, I2 = 2.25, k = 0.625, L = 2 with the material above, the run stops with exit code 255 right after “Factoring the system of equations using the symmetric spooles solver”. spooles.out says “matrix found to be singular”, and no error is printed.

Where it goes wrong

src/resultsmech_u1.f, lines 427–430, set up the shape-function derivatives (Luo’s Eqs. 19–20):

      aly=12.d0*e*xi11/(xk*um*a*dl*dl)
      alz=12.d0*e*xi22/(xk*um*a*dl*dl)
      bey=1.d0/(1.d0-aly)
      bez=1.d0/(1.d0-alz)

Together with the polynomials that use aly and alz, these are the standard exact Timoshenko shape functions with φ replaced by −φ. The exact functions have 1/(1 + φ) where these have 1/(1 − φ), and φ has the opposite sign in every polynomial.

The matrix in src/e_c3d_u1.f (lines 479–529) is exactly the strain energy of these shape functions; I checked this symbolically. For example, line 480:

            y2=(12.d0*e*xi11-xk*um*a*dl*dl)**2

This vanishes at φ = 1.

Flipping the sign of φ gives the standard Timoshenko matrix (12EI/((1 + φ)L³), (4 + φ)EI/((1 + φ)L), …). The changes are:

  1. In resultsmech_u1.f: aly=-12.d0*e*xi11/(xk*um*a*dl*dl), and alz likewise.
  2. In e_c3d_u1.f:
    • y2=(12.d0*e*xi11+xk*um*a*dl*dl)**2, and z2 likewise;
    • 15.d0 instead of 3.d0 in y3 and z3;
    • 6.d0 instead of 30.d0 in s(5,11) and s(6,12).

I re-implemented both routines in Python with these changes. The result is 0.515 for any number of elements, in both bending planes. Combined with the assembly fix from the separate U1 nodal-force issue, the nodal forces built from the section forces then equal K_e·u_e. The changes have not been built and tested in ccx itself.

Workaround

  • A large k makes the displacements those of an exact Euler–Bernoulli beam (k·A ≫ 12EI/(GL²) for every element). The S shear forces then have the right size but still the wrong sign.
  • For slender beams the error is small (−0.2 % in userbeam).
  • Otherwise, refine the mesh.