S4 and B31 ignore shell and beam temperature gradients (the expanded C3D8I takes the element's mean temperature): no thermal bending, no message

Analysis issues
1

Summary

For shells, the second value of a *TEMPERATURE data line is the temperature gradient through the thickness. For beams, the second value is the gradient in the 2-direction and the third value the gradient in the 1-direction (data-line format of *TEMPERATURE, section 7.132 of the 2.23 manual). *INITIAL CONDITIONS, TYPE=TEMPERATURE (section 7.77) takes the same values. On S8, S8R, B32 and B32R the gradient bends the element as expected. On S4 and B31 it has no effect; only the mean temperature is applied:

  • A free 2 x 1 plate with a gradient of 100 should curl by 0.1 at its far corner. With S8 or S8R it does (-0.1000000); with S4 the corner stays at 0.
  • A cantilever with a gradient of 100 in the 2-direction should deflect 0.05 at the tip (beam theory). With B32 or B32R it does (-0.05000029); with B31 the tip deflection is 0.

The manual promises gradients for these elements:

  • section 6.2.14 (S8 and S8R), to which section 6.2.12 (S4 and S4R) refers “for more information on shell elements”: “In addition to a temperature for the reference surface of the shell, a temperature gradient in normal direction can be specified on the *TEMPERATURE card.”
  • section 6.2.33 (B32 and B32R), to which section 6.2.32 (B31 and B31R) refers: “In addition to a temperature for the reference surface of the beam, a temperature gradient in 1-direction and in 2-direction can be specified on the *TEMPERATURE.”
  • *TEMPERATURE (section 7.132) and *INITIAL CONDITIONS (section 7.77) say the same, without exceptions.

ccx does spread the gradient onto the expanded nodes of S4 and B31 (gen3dtemp.f). The C3D8I element they are expanded into then gives all its integration points the element’s mean temperature (resultsmech.f). For solid C3D8 and C3D8I elements this averaging may well be intended. For expanded shells and beams it removes the gradient that the manual promises.

S4R and B31R (documentation note). They give no thermal bending either (corner 8.081101E-16, tip 4.577773E-16). There it follows from the element: both are expanded into C3D8R, whose single integration point at mid-thickness sees only the mean temperature. Section 6.2.32 says this about bending in general: “Since the C3D8R element has only one integration point in the middle of the element, bending effect cannot be taken into account. Therefore, the B31R element should not be used for bending.” Section 6.2.12 has no such sentence for S4R. A line saying that S4R and B31R ignore temperature gradients (in 6.2.12, 6.2.32 and at the gradient sentences of *TEMPERATURE and *INITIAL CONDITIONS), or a warning when a gradient is given for them, would help.

Version

ccx 2.23 (macOS arm64). Source: master 9ba5d76 (16 Oct 2026).

Already reported or fixed?

Not reported: I found no matching GitHub issue (open or closed) or pull request in Dhondtguido/CalculiX, and no calculix.discourse.group topic (searched 21 Oct 2026). Not fixed: resultsmech.f is unchanged on master (9ba5d76).

Minimal input

s4_gradient.inp: a flat 2 x 1 plate of two S4 elements (thickness 0.1, E 1000, nu 0, alpha 1E-3), supported statically determinate. The step sets a mean temperature of 10 and a gradient of 100.

*NODE, NSET=NALL
1, 0, 0, 0.
2, 1, 0, 0.
3, 2, 0, 0.
4, 0, 1, 0.
5, 1, 1, 0.
6, 2, 1, 0.
*ELEMENT, TYPE=S4, ELSET=EALL
1, 1, 2, 5, 4
2, 2, 3, 6, 5
*MATERIAL, NAME=M
*ELASTIC
1000., 0.
*EXPANSION
1.E-3
*SHELL SECTION, ELSET=EALL, MATERIAL=M
0.1
*BOUNDARY
1, 1, 3
2, 2, 3
4, 3, 3
*INITIAL CONDITIONS, TYPE=TEMPERATURE
NALL, 0.
*STEP
*STATIC
*TEMPERATURE
NALL, 10., 100.
*NODE PRINT, NSET=NALL
U
*END STEP

s4r_gradient.inp: the same with TYPE=S4R.

s8_gradient.inp: the same plate with two S8 elements; also run with TYPE=S8R.

s8_gradient.inp 
*NODE, NSET=NALL
1, 0, 0, 0.
2, 0.5, 0, 0.
3, 1, 0, 0.
4, 1.5, 0, 0.
5, 2, 0, 0.
6, 0, 0.5, 0.
7, 1, 0.5, 0.
8, 2, 0.5, 0.
9, 0, 1, 0.
10, 0.5, 1, 0.
11, 1, 1, 0.
12, 1.5, 1, 0.
13, 2, 1, 0.
*ELEMENT, TYPE=S8, ELSET=EALL
1, 1, 3, 11, 9, 2, 7, 10, 6
2, 3, 5, 13, 11, 4, 8, 12, 7
*MATERIAL, NAME=M
*ELASTIC
1000., 0.
*EXPANSION
1.E-3
*SHELL SECTION, ELSET=EALL, MATERIAL=M
0.1
*BOUNDARY
1, 1, 3
3, 2, 3
9, 3, 3
*INITIAL CONDITIONS, TYPE=TEMPERATURE
NALL, 0.
*STEP
*STATIC
*TEMPERATURE
NALL, 10., 100.
*NODE PRINT, NSET=NALL
U
*END STEP

b31_gradient.inp: a B31 cantilever along z (length 1, two elements, RECT 0.1 x 0.1), clamped at z = 0, n1 = x (so n2 = t x n1 = y), mean temperature 10, gradient 100 in the 2-direction (second value), 0 in the 1-direction (third value).

*NODE, NSET=NALL
1, 0., 0., 0
2, 0., 0., 0.5
3, 0., 0., 1
*ELEMENT, TYPE=B31, ELSET=EALL
1, 1, 2
2, 2, 3
*MATERIAL, NAME=M
*ELASTIC
1000., 0.
*EXPANSION
1.E-3
*BEAM SECTION, ELSET=EALL, MATERIAL=M, SECTION=RECT
0.1, 0.1
1., 0., 0.
*BOUNDARY
1, 1, 6
*INITIAL CONDITIONS, TYPE=TEMPERATURE
NALL, 0.
*STEP
*STATIC
*TEMPERATURE
NALL, 10., 100., 0.
*NODE PRINT, NSET=NALL
U
*END STEP

b31_gradient.inp is also run with TYPE=B31R.

b32_gradient.inp: the same beam as two B32 elements; also run with TYPE=B32R.

b32_gradient.inp 
*NODE, NSET=NALL
1, 0., 0., 0
2, 0., 0., 0.25
3, 0., 0., 0.5
4, 0., 0., 0.75
5, 0., 0., 1
*ELEMENT, TYPE=B32, ELSET=EALL
1, 1, 2, 3
2, 3, 4, 5
*MATERIAL, NAME=M
*ELASTIC
1000., 0.
*EXPANSION
1.E-3
*BEAM SECTION, ELSET=EALL, MATERIAL=M, SECTION=RECT
0.1, 0.1
1., 0., 0.
*BOUNDARY
1, 1, 6
*INITIAL CONDITIONS, TYPE=TEMPERATURE
NALL, 0.
*STEP
*STATIC
*TEMPERATURE
NALL, 10., 100., 0.
*NODE PRINT, NSET=NALL
U
*END STEP

How to run

ccx s4_gradient
ccx s8_gradient
ccx b31_gradient
ccx b32_gradient

Expected

Plate: with nu = 0 a linear temperature through the thickness bends the free plate stress-free to the curvature alpha g = 0.1 in both directions. With the supports at (0, 0), (1, 0) and (0, 1): w = -(alpha g / 2)(x² - x + y² - y), so w = -0.1 at (2, 0) and (2, 1). The mean temperature stretches it by alpha T = 0.01: u = 0.02 at x = 2.

Beam: curvature alpha g = 0.1 about the 1-axis, tip deflection alpha g L²/2 = 0.05 (in -y, away from the hot side), axial alpha T L = 0.01.

Actual

All runs exit 0 without a warning.

deck element node in-plane U1 (plate) / axial U3 (beam) bending deflection: U3 (plate) / U2 (beam)
s8_gradient S8, S8R 5 at (2, 0) 2.000000E-02 -1.000000E-01
s4_gradient S4 3 at (2, 0) 2.000000E-02 3.756002E-15
b32_gradient B32, B32R 5, tip 1.000000E-02 -5.000029E-02
b31_gradient B31 3, tip 1.000000E-02 -3.774921E-16

The mean temperature is applied; the gradient is lost. S4R (8.081101E-16) and B31R (4.577773E-16) behave the same; see the documentation note in the Summary.

With *NODE FILE, OUTPUT=3D and NT, the expanded nodes of the S4 plate carry -5 and +5 for a step line NALL, 0., 100., so the gradient does reach the expanded C3D8I nodes.

Where it goes wrong

  • gen3dtemp.f lines 214-262 (and 69-184 for initial values) give the expanded nodes of shells (‘L’) and beams (‘B’) the temperature plus or minus the gradient times the distance from the reference surface, for linear and quadratic elements alike.
  • S4 and B31 are expanded into C3D8I (sections 6.2.12 and 6.2.32). resultsmech.f computes the internal forces, and with them the thermal load of a linear step (linstatic.c line 484, b[k]=fext[k]-f[k]). There every integration point of an 8-node brick takes the element’s mean temperature:
              if((lakonl(4:5).eq.'8 ').or.
     &             (lakonl(4:5).eq.'8I')) then
                do i1=1,8
                  t0l=t0l+t0(konl(i1))/8.d0
                  t1l=t1l+t1(konl(i1))/8.d0
                enddo

(resultsmech.f lines 848-853; lines 868-873 for coupled analyses.) The mean of the bottom and top nodes is the reference temperature, so no thermal bending is left. For solid C3D8 and C3D8I elements a constant temperature per element may be deliberate. Expanded elements can be told apart by the flag gen3dtemp.f uses, lakon(7:7) = ‘L’ (shell) or ‘B’ (beam).

Fix

In resultsmech.f, skip the averaging for expanded shell and beam elements (lakonl(7:7) ‘L’ or ‘B’) and interpolate the temperature at each integration point with the shape functions, as the else branch does, at least through the thickness. Solid C3D8 and C3D8I can keep the mean if that is intended.

For S4R and B31R, a sentence in the manual or a warning (see the documentation note above).

Workaround: S8, S8R, B32 or B32R where gradients matter.