Summary
Section 7.66 (*GAP CONDUCTANCE) of the 2.23 manual: “The conductance coefficient can be defined as a function of the contact pressure and the mean temperature of slave and master surface.” Section 7.130 (*SURFACE BEHAVIOR) offers EXPONENTIAL, LINEAR, TABULAR, TIED and HARD pressure-overclosure laws.
The thermal contact spring routines compute the contact pressure themselves, from the clearance, to look up the conductance. In three of the six pair/law combinations they get that pressure wrong, so the conductance comes from the wrong point of the table. The mechanical contact pressure of the same run is right. In the deck below every law gives a contact pressure of 1.0E4 at the prescribed overclosure, so the conductance should be 100 and the heat flow 3333.333 in every case:
| pair | law | CPRESS (*CONTACT PRINT CSTR) | heat flow, ccx 2.23 and master | conductance this implies |
|---|---|---|---|---|
| NODE TO SURFACE | LINEAR | 1.000E+04 | 3750.000 | 150 (last table point) |
| NODE TO SURFACE | TABULAR | 1.000E+04 | 2500.000 | 50 (first table point) |
| NODE TO SURFACE | EXPONENTIAL | 1.000E+04 | 3333.333 | 100 |
| SURFACE TO SURFACE | LINEAR | 1.000E+04 | 3333.333 | 100 |
| SURFACE TO SURFACE | TABULAR | 0 (see PR 205) | 2505.885 | 50.24 |
| SURFACE TO SURFACE | EXPONENTIAL | 1.000E+04 | 3333.333 | 100 |
Exit 0 and no message in every case. There are two causes, one in the node-to-surface routines and one in the face-to-face routines (below). With both corrected, all six cases give 3333.333.
Version
ccx 2.23 (macOS arm64) and current master (537aa81, 2 Nov 2026): the same numbers to all printed digits.
Already reported or fixed?
Not reported: I found no GitHub issue or pull request in Dhondtguido/CalculiX about the thermal contact pressure or the gap conductance, and no matching calculix.discourse.group topic title (searched 3 Nov 2026).
Related: springforc_f2f, springstiff_f2f: use the 802-entry plconloc layout by oursland · Pull Request #205 · Dhondtguido/CalculiX · GitHub (open) corrects the same plconloc(81) read in the mechanical face-to-face routines springforc_f2f.f and springstiff_f2f.f (a TABULAR face-to-face pair transmits no force, hence CPRESS 0 in the table). It does not touch the thermal twins springforc_f2f_th.f and springstiff_f2f_th.f. With the change of PR 205 applied to master, the face-to-face TABULAR pair’s CSTR becomes 1.000000E+04, but its heat flow is still 2505.885.
Minimal input
gapk.inp (node-to-surface, LINEAR):
** Two unit cubes stacked at z = 1, conductivity 100 (each cube is a thermal
** resistance of 0.01). Every node of the upper cube is moved 0.001 down, so
** the overclosure is 0.001 everywhere and the contact pressure is 1.E4.
** Gap conductance 50 at pressure 0 and 150 at 2.E4, so 100 at 1.E4.
** Heat flow from T = 100 (top face) to T = 0 (bottom face):
** q = 100/(0.02 + 1/k): k = 100 -> 3333.333, k = 150 -> 3750, k = 50 -> 2500
*NODE, NSET=NALL
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.
9, 0., 0., 1.
10, 1., 0., 1.
11, 1., 1., 1.
12, 0., 1., 1.
13, 0., 0., 2.
14, 1., 0., 2.
15, 1., 1., 2.
16, 0., 1., 2.
*ELEMENT, TYPE=C3D8, ELSET=LOW
1, 1, 2, 3, 4, 5, 6, 7, 8
*ELEMENT, TYPE=C3D8, ELSET=UP
2, 9, 10, 11, 12, 13, 14, 15, 16
*NSET, NSET=NLOW
1, 2, 3, 4, 5, 6, 7, 8
*NSET, NSET=NUP
9, 10, 11, 12, 13, 14, 15, 16
*NSET, NSET=BOT
1, 2, 3, 4
*NSET, NSET=TOP
13, 14, 15, 16
*MATERIAL, NAME=M
*ELASTIC
1.E6, 0.
*CONDUCTIVITY
100.
*SOLID SECTION, ELSET=LOW, MATERIAL=M
*SOLID SECTION, ELSET=UP, MATERIAL=M
*SURFACE, NAME=LTOP
LOW, S2
*SURFACE, NAME=UBOT
UP, S1
*SURFACE INTERACTION, NAME=SI
*SURFACE BEHAVIOR, PRESSURE-OVERCLOSURE=LINEAR
1.E7, 1.E-3
*GAP CONDUCTANCE
50., 0.
150., 2.E4
*CONTACT PAIR, INTERACTION=SI, TYPE=NODE TO SURFACE
UBOT, LTOP
*INITIAL CONDITIONS, TYPE=TEMPERATURE
NALL, 0.
*BOUNDARY
NLOW, 1, 3
NUP, 1, 2
*STEP, NLGEOM
*COUPLED TEMPERATURE-DISPLACEMENT, STEADY STATE
1., 1.
*BOUNDARY
NUP, 3, 3, -0.001
BOT, 11, 11, 0.
TOP, 11, 11, 100.
*NODE PRINT, NSET=TOP, TOTALS=ONLY
RFL
*CONTACT PRINT
CSTR
*END STEP
The other five rows of the table replace the *SURFACE BEHAVIOR card and its data line, and/or TYPE=NODE TO SURFACE by TYPE=SURFACE TO SURFACE. All three laws give a pressure of 1.E4 at an overclosure of 0.001:
*SURFACE BEHAVIOR, PRESSURE-OVERCLOSURE=TABULAR
0., 0.
1.E7, 1.
*SURFACE BEHAVIOR, PRESSURE-OVERCLOSURE=EXPONENTIAL
1.E-3, 100.
(EXPONENTIAL: p = p0·100^(overclosure/c0) = 100·100 = 1.E4.)
How to run
ccx gapk
Expected
total heat generation for set TOP 3.333333E+03 for every law and pair type, as the EXPONENTIAL rows and the face-to-face LINEAR row give.
Actual
ccx 2.23 and master 537aa81, last increment:
LINEAR, NODE TO SURFACE: 3.750000E+03
TABULAR, NODE TO SURFACE: 2.500000E+03
TABULAR, SURFACE TO SURFACE: 2.505885E+03
The other three cases print 3.333333E+03. Every run ends with “Job finished” and exit 0. Apart from the *INFO lines for NLGEOM and for *GAP CONDUCTANCE (below), the only message is the warning “no conductivity constants were assigned to material SI”, which all six runs print. *CONTACT PRINT CSTR shows a pressure of 1.000000E+04 (9.999999E+03 for LINEAR node-to-surface) at every slave node or integration point, except for the face-to-face TABULAR pair, where it shows 0 (the mechanical bug of PR 205).
Where it goes wrong
1. Node-to-surface: the law is chosen by the sign of elcon(1). src/springforc_n2f_th.f (master line 111) and src/springstiff_n2f_th.f (line 114):
if(elcon(1,1,imat).gt.0.d0) then
! exponential overclosure
...
pressure=dexp(-beta*val+dlog(alpha))
else
! linear overclosure
pi=4.d0*datan(1.d0)
eps=elcon(1,1,imat)*pi/elcon(2,1,imat)
pressure=-elcon(2,1,imat)*val*
& (0.5d0+datan(-val/eps)/pi)
endif
surfacebehaviors.f stores the law in elcon(3) (its header, lines 26-29: 1 exponential, 2 linear, 3 tabular, 4 tied). For LINEAR, elcon(1) is σ∞, which surfacebehaviors.f forces to be positive (default: the first elastic constant / 70000). So LINEAR always takes the exponential branch, with α = K and β = σ∞: the pressure comes out as about K = 1.E7, and the conductance is taken at the last point of its table. TABULAR leaves elcon(1) = elcon(2) = 0, so it takes the linear branch with eps = 0·π/0 = NaN. The pressure is NaN, ident returns 0, and the conductance is taken at the first point. The mechanical routine springforc_n2f.f (lines 269, 288, 331) and both face-to-face thermal routines switch on int(elcon(3,1,imat)) instead.
2. Face-to-face TABULAR: the curve length is read from plconloc(81). src/springforc_f2f_th.f (line 193) and src/springstiff_f2f_th.f (line 181) take niso=int(plconloc(81)) after materialdata_sp, which stores the length in plconloc(801). These two routines already declare plconloc(802), but they copy the curve into xiso(20), yiso(20) (the node-to-surface routines use 200, the length plcopy allows). This is the same mistake that PR 205 corrects in the mechanical routines.
I checked both causes with a patched build of master (PR 205’s change plus two more): switching the two node-to-surface routines on int(elcon(3,1,imat)), with a TABULAR branch as in springforc_n2f.f, and reading plconloc(801) into 200-entry arrays in the two face-to-face routines. All six cases then give 3.333333E+03. With PR 205’s change alone they give the numbers of the table.
Fix
springforc_n2f_th.f,springstiff_n2f_th.f: select the law withint(elcon(3,1,imat))(1 exponential, 2 linear, 3 tabular), asspringforc_n2f.fdoes. The TABULAR branch needspliconandnplicon, which these two routines do not receive yet; both callers (resultstherm.fline 319,e_c3d_th.fline 250) have them, and pass them to the face-to-face routines already.springforc_f2f_th.f,springstiff_f2f_th.f:niso=int(plconloc(801))andxiso(200),yiso(200), as in PR 205.
Also (smaller, same routines)
Without NLGEOM on the *STEP card, the node-to-surface thermal routine works with the undeformed positions (springforc_n2f_th.f lines 44-60: “for geometrically linear static calculations the undeformed position is taken”). It then evaluates the pressure at a clearance of 0, while the mechanical routine always uses the deformed positions (springforc_n2f.f lines 50-58: “otherwise no contact force”). In the deck above without NLGEOM, node-to-surface EXPONENTIAL gives 2.512438E+03 (conductance 50.5, i.e. the pressure p0 = 100 at zero clearance) although CSTR shows 1.000000E+04. Face-to-face EXPONENTIAL still gives 3.333333E+03. The run prints
*INFO reading *GAP CONDUCTANCE: nonlinear geometric
effects are turned on
but gapconductances.f (lines 64-65) sets only iperturb(1)=2; the line iperturb(2)=1 is commented out. Either the message or the geometry should change.