*FREQUENCY step leaves the last eigenmode's strain energy in the internal energy (ELSE), and a following *DYNAMIC or NLGEOM *STATIC step adds to it

Analysis issues
1

Summary

After a *FREQUENCY step, the element internal energy (*EL PRINT ELSE) holds the strain energy of the last computed eigenmode. The modes are mass-normalised, so this is half the eigenvalue: 5.0E+03 for a cube whose only mode has eigenvalue 1.0E+04. A following implicit *DYNAMIC step, or a *STATIC step with NLGEOM, accumulates its energy on top of that value. Its ELSE output is then 5.000000E+03 instead of a few 1E-07 (dynamic) or 5.0E-04 (static). With two modes it is 1.5E+04, half the second eigenvalue 3.0E+04. This happens whether the frequency step itself prints ELSE or not.

A *FREQUENCY step changes nothing in the model state: it has no load and no displacement. The same *DYNAMIC step without the frequency step before it starts from zero, and a linear *STATIC step after the frequency step gives the right 5.0E-04. arpack.c restores the state variables xstate at the end of each mode’s pass, so the mode shapes are not meant to stay in the model’s state, but it does not restore the energies.

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: the mode loop of arpack.c is unchanged on master (9ba5d76).

Related:

  • A related energy problem with a different cause: the internal energy at the start of an implicit *DYNAMIC step after a strained static step is doubled (reported separately).

Minimal input

freq_dyn.inp is one C3D8 cube. Only the z-displacements of the four top nodes are free (E = 1000, nu = 0, density 0.3). Step 1 is a frequency step with one mode, and it requests ELSE, because ccx needs energy requests in the first step. Step 2 is a dynamic step from rest under 4 x 0.25 on the top.

*NODE
1, 0., 0., 0.
2, 1., 0., 0.
3, 1., 1., 0.
4, 0., 1., 0.
5, 0., 0., 1.
6, 1., 0., 1.
7, 1., 1., 1.
8, 0., 1., 1.
*ELEMENT, TYPE=C3D8, ELSET=EALL
1, 1, 2, 3, 4, 5, 6, 7, 8
*NSET, NSET=BOT
1, 2, 3, 4
*NSET, NSET=TOP
5, 6, 7, 8
*MATERIAL, NAME=M
*ELASTIC
1000., 0.
*DENSITY
0.3
*SOLID SECTION, ELSET=EALL, MATERIAL=M
*BOUNDARY
BOT, 1, 3
TOP, 1, 2
*STEP
*FREQUENCY
1
*EL PRINT, ELSET=EALL
ELSE
*END STEP
*STEP
*DYNAMIC, DIRECT
2.51327412287E-4, 2.51327412287E-3
*CLOAD
TOP, 3, 0.25
*END STEP

Variants:

  • dyn_only.inp: step 1 removed, and the *EL PRINT, ELSET=EALL / ELSE lines moved into the dynamic step.
  • freq2_dyn.inp: 2 modes instead of 1 on the *FREQUENCY data line.
  • freq_nlgeom.inp: step 2 replaced by *STEP, NLGEOM / *STATIC / 0.25, 1. / *CLOAD / TOP, 3, 0.25 / *END STEP.
  • freq_static.inp: the same as freq_nlgeom.inp without NLGEOM and without the 0.25, 1. line.
  • static_freqS_dyn.inp: a *STATIC step first (*CLOAD / TOP, 3, 0.25, *EL PRINT, ELSET=EALL / ELSE), then the frequency step with *EL PRINT, ELSET=EALL / S instead of ELSE, then the dynamic step with *EL PRINT, ELSET=EALL / ELSE.

How to run

ccx freq_dyn (and the variants the same way)

Expected

In step 2 of freq_dyn.inp, ELSE is the strain energy of the dynamic response alone, as in dyn_only.inp: 4.985689E-11 after the first increment and 4.933924E-07 at the end of the step. freq_nlgeom.inp should give about k u^2 / 2 = 5.0E-04, as freq_static.inp does.

Actual

freq_dyn.dat: the frequency step prints eigenvalue 0.1000000E+05. Every increment of the dynamic step then prints

 internal energy (element, energy) for set EALL and time  0.2513274E-03

         1  5.000000E+03

up to the end of the step. That is 1.0E+04 / 2, the strain energy of the mass-normalised mode (phi^T K phi / 2 = omega^2 / 2), plus the dynamic step’s own 1E-07. ccx’s energy balance on stdout for the dynamic step starts from the same value:

 initial energy (at start of step) = 5.000000e+03
 ...
 internal energy = 5.000000e+03
 kinetic energy = 3.157227e-07

Results per deck:

deck eigenvalues ELSE in the last step
dyn_only.inp 4.933924E-07
freq_dyn.inp 1.0E+04 5.000000E+03
freq2_dyn.inp 1.0E+04, 3.0E+04 1.500000E+04
freq_static.inp 1.0E+04 5.000000E-04
freq_nlgeom.inp 1.0E+04 5.000000E+03
static_freqS_dyn.inp 1.0E+04 5.000000E+03

Where it goes wrong

arpack.c loops over the modes (line 1046). For each mode it calls results() with the mode shape (line 1130, or 1154 when iperturb!=0) to compute the mode’s stresses for output. It passes the model’s energy array ener, together with enerini, emeini and stiini, which are freshly allocated as zero for each mode (lines 1121-1126). With iout=1, resultsmech.f sets ener = enerini + (emec-emeini)*(stre+stiini)/2, i.e. the mode’s strain energy, so after the loop ener holds the last mode’s value. The state variables get restored at the end of each mode’s pass, inside the loop (lines 1218-1222; the loop closes at line 1223):

    if(*iperturb!=0){
      for(k=0;k<*nstate_*mi[0]*(ne0+maxprevcontel);++k){
	xstate[k]=xstateini[k];
      }	  
    }

ener does not. nonlingeo.c starts the next step with enerini = ener (line 855), so every later energy is offset by the mode’s energy. linstatic.c starts from a zeroed enerini (line 576), which is why the linear static step is right. arpackbu.c (*BUCKLE) has the same kind of per-mode results() calls with ener (lines 790, 813), not checked here.

Fix

In arpack.c, keep a copy of ener before the mode loop and copy it back at the end of each mode’s pass, next to the xstate restore (or once after the loop). Alternatively, pass a scratch array to the per-mode results() calls. Not built or tested here.

Workaround: none inside one deck. Run the frequency analysis as a separate job.