Forum topic:
Summary
In a *DYNAMIC, EXPLICIT step, ccx sizes the increment of a quadratic beam with a circular section as if the section were √2 times thicker than it is. The increment it selects is therefore unstable:
- The displacements exceed 1e66 by increment 200 and are NaN from increment 800 on.
- The job still prints “Job finished” and exits with code 0.
A RECT section of the same width gets a √2 smaller increment and runs fine.
Version
ccx 2.23 (conda-forge build pl5321h484372f_2, macOS arm64). On master, calcstabletimeincvol.f has changed since 2.23 only for C3D8R, and gen3dfrom1d.f has not changed. The problem should therefore still be there.
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.
Minimal input
circ_explicit.inp:
** One B32R cantilever, length 40, CIRC section of diameter 20, steel (N, mm, s).
** Tip load 1000 in y applied at t = 0 and held, explicit dynamics, 2e-3 s.
*NODE, NSET=NALL
1, 0.0, 0.0, 0.0
2, 20.0, 0.0, 0.0
3, 40.0, 0.0, 0.0
*ELEMENT, TYPE=B32R, ELSET=EALL
1, 1, 2, 3
*MATERIAL, NAME=STEEL
*ELASTIC
210000.0, 0.3
*DENSITY
7.85e-9
*BEAM SECTION, ELSET=EALL, MATERIAL=STEEL, SECTION=CIRC
20.0, 20.0
0.0, 0.0, -1.0
*BOUNDARY
1, 1, 6
*STEP, INC=100000
*DYNAMIC, EXPLICIT
, 2.e-3
*CLOAD
3, 2, 1000.0
*NODE PRINT, NSET=NALL, FREQUENCY=200
U
*END STEP
How to run
ccx circ_explicit
Then cap the increment below ccx’s own by changing the *DYNAMIC data line to , 2.e-3, , 6.e-7, and run again.
Expected
The increment ccx selects should be stable, so the tip displacement stays bounded. For a load applied suddenly and then held, it oscillates between 0 and about twice the static deflection. Euler-Bernoulli gives PL³/(3EI) = 1000·40³/(3·210000·7854) = 0.0129.
With the increment capped at 6e-7 the run does exactly that. Over 3334 increments the tip u_y runs from 3.8e-6 to 0.0264, with a mean of 0.0131, and there is no NaN.
The manual (6.9.7) bases the estimate on ω_max ≈ 2c/h_min, “where … h_min is the minimum height of the element”. The expanded section here is 20 wide, so with h = 20 the formula gives 5.90e-7. That is the increment ccx itself selects for SECTION=RECT, 20.0, 20.0, and that run is stable.
Actual
stdout:
Explicit time integration: Volumetric COURANT initial stable time increment:8.341279e-07
SELECTED time increment (not considering penalty contact):8.341279e-07
...
Job finished
The exit code is 0. .dat:
INCREMENT 200
displacements (vx,vy,vz) for set NALL and time 0.1668256E-03
1 0.000000E+00 0.000000E+00 0.000000E+00
2 1.513252E+66 4.967292E+65 6.496154E+66
3 -7.070765E+66 -3.819730E+67 -4.990000E+68
...
INCREMENT 800
displacements (vx,vy,vz) for set NALL and time 0.6673024E-03
1 NaN NaN NaN
2 NaN NaN NaN
3 NaN NaN NaN
Every later output, up to increment 2398 at t = 2e-3, is NaN. B32 behaves the same way.
I bisected the stable limit by capping the increment over a step of 1e-2 s:
- This one-element deck is stable at 7.5e-7 for 13,333 increments and turns NaN at 7.6e-7.
- A cantilever of four such elements (4 x 40) is stable at 6.3e-7 and turns NaN at 6.4e-7.
In both cases ccx selects 8.34e-7, which is 1.11 and 1.32 times the limit.
Where it goes wrong
calcstabletimeincvol.f takes h as the smallest volume/face-area ratio. Both are evaluated with a one-point rule, at the element centre and at the face centres:
volume=weight*xsj
...
area=weight*dsqrt(xsj2(1)*xsj2(1)+xsj2(2)*xsj2(2)+
& xsj2(3)*xsj2(3))
hmin=min(hmin,(volume*volfac/area))
...
smscale(nelem)=critom*hmin*elemfac/(2.d0*wavspd)
For SECTION=CIRC, gen3dfrom1d.f places the C3D20(R) corner nodes at ±d/(2√2) (sc=.5d0/dsqrt(2.d0)), at 45 degrees on the circle, and the midside nodes at ±d/2. The two one-point quantities then come out as:
- The Jacobian at the centre is that of the box around the circle, so volume = l·d² (16000 here). The true volume of the expanded element is 12418.
- Each side face is measured by the chord between its corner nodes, so area = l·d/√2.
Together they give h = min(l, √2·d) = 28.3 instead of about d = 20. Then dt = 1.9172·28.28·0.3/(2·6.001e6)·0.8/1.3 = 8.34e-7, the value printed above. With the volume and the face areas integrated properly, h = 12418/625.0 = 19.9 and dt = 5.86e-7, which is stable in both models above. Taking h from the section size for CIRC beams would also work.
Separately, the explicit run carries on past NaN and reports success, so a diverged run looks like a successful one.
Workaround
Either of these avoids the NaN:
- Cap the increment with the fourth field of the
*DYNAMICdata line, at half of ccx’s own value for CIRC beams (4.17e-7 here, the estimate for a square of side d/√2). - Model the solid circle as
SECTION=PIPE,10.0, 10.0(wall thickness equal to the radius). PIPE is expanded as the square around the section, so ccx selects 5.90e-7. That run is stable, with a mean tip deflection of 0.0136 against 0.0131 for CIRC.