Resources Numerical accuracy Calculation policy
Solves the temperature field of a part with fixed temperatures, heat flux, power, convection and radiation, and reports temperatures, heat flux and the heat balance.
Temperatures, matrices and loads are Float64Array. The linear system is solved by preconditioned conjugate gradients to a relative residual of 1e-10, radiation is re-linearized until the temperature changes by less than 1e-6 K, and a transient run steps in time with the θ-method, so the mesh, the iteration tolerances and the time step set the accuracy. Study values arrive in SI units and are converted with powers of ten.
20 calculations in 1 source file. The code is read from the released app when this page loads. It is the calculation people run, not a copy kept for this page.
Tests tools/partcraft-test/test-thermal.mjs
The radiation constant in the solver's units, W/(mm²·K⁴).
σ = 5.670374419 × 10⁻¹⁴
Written as a literal with 10 significant digits, already scaled from m² to mm². Radiation results cannot carry more digits than this constant.
export const SIGMA_SB = 5.670374419e-14; // Stefan–Boltzmann, W/(mm²·K⁴)The factors that convert conductivity, flux, film coefficient, heat generation and heat capacity from SI to the solver's millimeter units.
k × 10⁻³; flux and h × 10⁻⁶; generation and ρc × 10⁻⁹
The factors are exact powers of ten in decimal but not in binary, so each conversion is one multiplication with one rounding. No rounded physical factor is involved.
export const SI = { k: 1e-3, flux: 1e-6, h: 1e-6, gen: 1e-9, rc: 1e-9 };The conduction, capacity and thermal-load matrices of a unit cube, scaled for every voxel element.
K = Σ_g ∇Nₐ·∇N_b ÷ 8; M = Σ_g NₐN_b ÷ 8
Integrated once with 2 × 2 × 2 Gauss points at ±1 ÷ √3 into Float64Array. Every voxel then multiplies these by its side length, so all elements share the same rounding.
const UNIT = (() => {
const K = new Float64Array(64), M = new Float64Array(64), TH = new Float64Array(24 * 8);
for (let g = 0; g < 8; g++) for (let a = 0; a < 8; a++) {
for (let b = 0; b < 8; b++) {
K[a * 8 + b] += (GG[g][a][0] * GG[g][b][0] + GG[g][a][1] * GG[g][b][1] + GG[g][a][2] * GG[g][b][2]) / 8;
M[a * 8 + b] += GN[g][a] * GN[g][b] / 8;
for (let i = 0; i < 3; i++) TH[(3 * a + i) * 8 + b] += GG[g][a][i] * GN[g][b] / 8; // ×E α h²/(1−2ν)
}
}
return { K, M, TH };
})();The share of a thermal load carried by the element's internal bubble modes, and their thermal stress.
f = E α h² (TH ÷ (1 − 2ν) + Aᵀ Fᵢ) ΔT, A = −Kᵢᵢ⁻¹ Kᵢc
The 9 × 9 bubble stiffness is solved with solveSmall from fea.js and the result is cached per Poisson ratio. The factor 1 ÷ (1 − 2ν) is not guarded, so ν = 0.5 divides by zero.
function bubbleThermal(nu) {
if (bubCache.has(nu)) return bubCache.get(nu);
const D = elasticD(nu), c2 = 1 / (1 - 2 * nu);
const bRow = (dN, nf) => { const B = new Float64Array(6 * 3 * nf), c = 3 * nf; dN.forEach(([dx, dy, dz], a) => { B[3 * a] = dx; B[c + 3 * a + 1] = dy; B[2 * c + 3 * a + 2] = dz; B[3 * c + 3 * a] = dy; B[3 * c + 3 * a + 1] = dx; B[4 * c + 3 * a + 1] = dz; B[4 * c + 3 * a + 2] = dy; B[5 * c + 3 * a] = dz; B[5 * c + 3 * a + 2] = dx; }); return B; };
const dNi = (x, y, z) => [[-4 * x, 0, 0], [0, -4 * y, 0], [0, 0, -4 * z]];
const Kii = new Float64Array(81), Kic = new Float64Array(9 * 24), Fi = new Float64Array(9 * 8);
GPTS.forEach(([x, y, z], g) => {
const Bi = bRow(dNi(x, y, z), 3), Bc = bRow(GG[g], 8);
const DBi = new Float64Array(54), DBc = new Float64Array(144);
for (let r = 0; r < 6; r++) for (let k = 0; k < 6; k++) { const d = D[r * 6 + k]; if (!d) continue; for (let j = 0; j < 9; j++) DBi[r * 9 + j] += d * Bi[k * 9 + j]; for (let j = 0; j < 24; j++) DBc[r * 24 + j] += d * Bc[k * 24 + j]; }
for (let i = 0; i < 9; i++) {
for (let r = 0; r < 6; r++) { const b = Bi[r * 9 + i]; if (!b) continue; for (let j = 0; j < 9; j++) Kii[i * 9 + j] += b * DBi[r * 9 + j] / 8; for (let j = 0; j < 24; j++) Kic[i * 24 + j] += b * DBc[r * 24 + j] / 8; }
const m = (Bi[i] + Bi[9 + i] + Bi[18 + i]) * c2;
for (let b = 0; b < 8; b++) Fi[i * 8 + b] += m * GN[g][b] / 8;
}
});
const A = solveSmall(Kii, 9, Kic, 24).map((v) => -v), Gi = solveSmall(Kii, 9, Fi, 8);
const TH = new Float64Array(24 * 8);
for (let r = 0; r < 24; r++) for (let b = 0; b < 8; b++) { let t = UNIT.TH[r * 8 + b] * c2; for (let i = 0; i < 9; i++) t += A[i * 24 + r] * Fi[i * 8 + b]; TH[r * 8 + b] = t; }
// stress per unit Eα at points: D·Bi(p)·Gi (6×8)
const at = (pts) => { const out = new Float64Array(pts.length * 48); pts.forEach(([x, y, z], p) => { const Bi = bRow(dNi(x, y, z), 3); for (let r = 0; r < 6; r++) for (let b = 0; b < 8; b++) { let t = 0; for (let k = 0; k < 6; k++) { const d = D[r * 6 + k]; if (d) for (let j = 0; j < 9; j++) t += d * Bi[k * 9 + j] * Gi[j * 8 + b]; } out[p * 48 + r * 8 + b] = t; } }); return out; };
const out = { TH, gauss: at(GPTS), points: at(XI.map((_, a) => [XI[a], ETA[a], ZETA[a]]).concat([[0, 0, 0]])) };
bubCache.set(nu, out);
return out;
}Conduction, capacity, heat load, gradient, thermal load, stresses and geometric stiffness of one voxel element.
K = k h K_unit; C = ρc h³ M_unit (lumped: ρc h³ ÷ 8); ∇T = Σ G₀ Tₐ ÷ h
The side h is read from the element's own coordinates as a difference of two node coordinates, then used as a scale; a part far from the origin loses digits in that difference. The thermal load and thermal stress divide by 1 − 2ν without a guard.
function hexFor(type) {
return {
type, npe: 8, nip: 8,
conductivity(X, k, out) { const s = k * side(X); for (let i = 0; i < 64; i++) out[i] = s * UNIT.K[i]; return out; },
capacity(X, rc, out, lumped) { const h = side(X), v = rc * h * h * h; if (lumped) out.fill(v / 8, 0, 8); else for (let i = 0; i < 64; i++) out[i] = v * UNIT.M[i]; return out; },
heatLoad(X, Q, out) { const h = side(X); out.fill(Q * h * h * h / 8, 0, 8); return out; },
gradient(X, Te, out) { const h = side(X); out[0] = out[1] = out[2] = 0; for (let a = 0; a < 8; a++) for (let i = 0; i < 3; i++) out[i] += G0[a][i] * Te[a] / h; return out; },
thermalLoad(X, E, nu, alpha, dT, out) {
const h = side(X), bub = type === 'h8i', s = E * alpha * h * h / (bub ? 1 : 1 - 2 * nu), M = bub ? bubbleThermal(nu).TH : UNIT.TH;
for (let r = 0; r < 24; r++) { let t = 0; for (let b = 0; b < 8; b++) t += M[r * 8 + b] * dT[b]; out[r] = s * t; }
return out;
},
stressPoints(X, E, nu, ue, out, th = null) {
const h = side(X), S = hexElement(nu, type).S, sc = new Float64Array(48), f = E / h;
for (let c = 0; c < 8; c++) for (let r = 0; r < 6; r++) { let t = 0; const b = c * 144 + r * 24; for (let j = 0; j < 24; j++) t += S[b + j] * ue[j]; sc[c * 6 + r] = f * t; }
for (let g = 0; g < 8; g++) {
let dTg = 0;
for (let r = 0; r < 6; r++) { let t = 0; for (let c = 0; c < 8; c++) t += GN[g][c] * sc[c * 6 + r]; out[g * 6 + r] = t; }
if (th) {
for (let c = 0; c < 8; c++) dTg += GN[g][c] * th.dT[c];
const st = E * th.alpha * dTg / (1 - 2 * nu); out[g * 6] -= st; out[g * 6 + 1] -= st; out[g * 6 + 2] -= st;
if (type === 'h8i') { const C = bubbleThermal(nu).gauss; for (let r = 0; r < 6; r++) { let t = 0; for (let b = 0; b < 8; b++) t += C[g * 48 + r * 8 + b] * th.dT[b]; out[g * 6 + r] += E * th.alpha * t; } }
}
}
return out;
},
geometricStiffness(X, sig, out) {
const h = side(X);
out.fill(0, 0, 64);
for (let g = 0; g < 8; g++) {
const s = sig.subarray ? sig.subarray(6 * g, 6 * g + 6) : sig.slice(6 * g, 6 * g + 6);
const m = [s[0], s[3], s[5], s[3], s[1], s[4], s[5], s[4], s[2]]; // 3×3, Voigt xy yz zx
const G = GG[g];
for (let a = 0; a < 8; a++) {
const t0 = m[0] * G[a][0] + m[3] * G[a][1] + m[6] * G[a][2], t1 = m[1] * G[a][0] + m[4] * G[a][1] + m[7] * G[a][2], t2 = m[2] * G[a][0] + m[5] * G[a][1] + m[8] * G[a][2];
for (let b = 0; b < 8; b++) out[a * 8 + b] += (t0 * G[b][0] + t1 * G[b][1] + t2 * G[b][2]) * h / 8;
}
}
return out;
},
face: {
npf: 4,
mass(Xf, out) { const A = quadArea(Xf); for (let i = 0; i < 16; i++) out[i] = A * QUAD_M[i] / 36; return out; },
load(Xf, out) { out.fill(quadArea(Xf) / 4, 0, 4); return out; },
},
};
}The area of a four-node boundary facet.
A = |d₁ × d₂| ÷ 2 (d₁, d₂ the diagonals)
Half the length of the cross product of the two diagonals, with Math.hypot. Exact in form for a planar quad, which voxel faces are.
function quadArea(Xf) { // planar quad: half the cross product of its diagonals
const d1 = [Xf[6] - Xf[0], Xf[7] - Xf[1], Xf[8] - Xf[2]], d2 = [Xf[9] - Xf[3], Xf[10] - Xf[4], Xf[11] - Xf[5]];
return Math.hypot(d1[1] * d2[2] - d1[2] * d2[1], d1[2] * d2[0] - d1[0] * d2[2], d1[0] * d2[1] - d1[1] * d2[0]) / 2;
}The mesh volume and the area of every exposed voxel face.
volume = elements × h³; facet area = h²
The volume is a count times h cubed, one product with no accumulation. Node coordinates are origin + index × h, so they do not drift across the grid.
const nn = (grid.nx + 1) * (grid.ny + 1) * (grid.nz + 1), xyz = new Float64Array(3 * nn);
const NX = grid.nx + 1, NXY = NX * (grid.ny + 1);
for (let n = 0; n < nn; n++) { const i = n % NX, j = Math.floor(n / NX) % (grid.ny + 1), k = Math.floor(n / NXY); xyz[3 * n] = grid.origin[0] + i * grid.h; xyz[3 * n + 1] = grid.origin[1] + j * grid.h; xyz[3 * n + 2] = grid.origin[2] + k * grid.h; }
const active = VX.activeNodes(grid, solidAt);
const surf = grid.surf;
const quads = surf ? surf.quads : new Int32Array(0);
return {
type, element: ELEMENTS[type], grid, elems: Int32Array.from(ids), nn, xyz, active, ne, npe: 8, conn, volume: ne * grid.h ** 3,
facets: { npf: 4, nodes: quads, area: new Float64Array(quads.length / 4).fill(grid.h * grid.h) },Factors the conduction matrix on its own pattern to speed up the temperature solve.
A pivot that is not above 0 restarts the factorization with the diagonal scaled by 1 + shift (1e-3, then × 4, 8 attempts). If all fail, the preconditioner falls back to dividing by the diagonal. The shift changes only the preconditioner, not the system solved.
/** Incomplete Cholesky IC(0) on the lower part of a CSR matrix → apply(r, z); a diagonal shift is added on breakdown. */
export function ic0(P, val) {
const { n, rowPtr, cols } = P;
const diag = new Int32Array(n);
for (let r = 0; r < n; r++) diag[r] = pos(P, r, r);
let shift = 0;
for (let attempt = 0; attempt < 8; attempt++) {
const L = new Float64Array(val.length); // lower entries (k ≤ diag) used
let ok = true;
for (let i = 0; i < n && ok; i++) {
for (let p = rowPtr[i]; p <= diag[i]; p++) {
const j = cols[p];
let s = val[p] * (j === i ? 1 + shift : 1);
// s −= Σ_{k<j} L_ik L_jk (sorted merge of rows i and j below j)
let a = rowPtr[i], b = rowPtr[j];
while (a < p && b < diag[j]) { const ca = cols[a], cb = cols[b]; if (ca === cb) { s -= L[a] * L[b]; a++; b++; } else if (ca < cb) a++; else b++; }
if (j === i) { if (!(s > 0)) { ok = false; break; } L[p] = Math.sqrt(s); } else L[p] = s / L[diag[j]];
}
}
if (ok) {
const y = new Float64Array(n);
return {
shift,
apply(r, z) {
for (let i = 0; i < n; i++) { let s = r[i]; for (let p = rowPtr[i]; p < diag[i]; p++) s -= L[p] * y[cols[p]]; y[i] = s / L[diag[i]]; }
z.set(y);
for (let i = n - 1; i >= 0; i--) { z[i] /= L[diag[i]]; const zi = z[i]; for (let p = rowPtr[i]; p < diag[i]; p++) z[cols[p]] -= L[p] * zi; }
return z;
},
};
}
shift = shift ? shift * 4 : 1e-3;
}
const inv = Float64Array.from(diag, (d) => 1 / val[d]); // Jacobi as the last resort
return { shift: Infinity, apply(r, z) { for (let i = 0; i < n; i++) z[i] = r[i] * inv[i]; return z; } };
}Solves the heat system for the temperature of every free node.
stop when ‖r‖ ÷ ‖b‖ < 10⁻¹⁰
The residual is updated recursively and tested against the norm of the right-hand side; default tolerance 1e-10, iteration cap max(1000, 100 × round(∛n)). The loop also leaves when p·Ap is not positive, which is then reported as not converged.
/** PCG on CSR → { x, iterations, residual, converged }. */
export function csrPcg(P, val, b, { x0 = null, tol = 1e-10, maxIter = 0, prec = null, cancelled = null } = {}) {
const n = P.n, M = prec || ic0(P, val);
const x = x0 ? Float64Array.from(x0) : new Float64Array(n), r = new Float64Array(n), z = new Float64Array(n), p = new Float64Array(n), q = new Float64Array(n);
let bn = 0; for (let i = 0; i < n; i++) bn += b[i] * b[i];
bn = Math.sqrt(bn);
if (bn === 0) return { x: new Float64Array(n), iterations: 0, residual: 0, converged: true };
csrMul(P, val, x, q);
for (let i = 0; i < n; i++) r[i] = b[i] - q[i];
M.apply(r, z); p.set(z);
let rz = 0; for (let i = 0; i < n; i++) rz += r[i] * z[i];
const cap = maxIter || Math.max(1000, 10 * Math.round(Math.cbrt(n)) * 10);
let it = 0, res = 1;
for (; it < cap; it++) {
let rn = 0; for (let i = 0; i < n; i++) rn += r[i] * r[i];
res = Math.sqrt(rn) / bn;
if (res < tol) break;
if (cancelled && (it & 31) === 0 && cancelled()) throw new Error('Cancelled');
csrMul(P, val, p, q);
let pq = 0; for (let i = 0; i < n; i++) pq += p[i] * q[i];
if (!(pq > 0)) break;
const a = rz / pq;
for (let i = 0; i < n; i++) { x[i] += a * p[i]; r[i] -= a * q[i]; }
M.apply(r, z);
let rz2 = 0; for (let i = 0; i < n; i++) rz2 += r[i] * z[i];
const beta = rz2 / rz; rz = rz2;
for (let i = 0; i < n; i++) p[i] = z[i] + beta * p[i];
}
return { x, iterations: it, residual: res, converged: res < tol };
}Turns a total power on picked faces into a flux per unit area.
q = P ÷ A (meshed area); flux input: q × scale
The area is the sum of the facets' load integrals. A power is divided by that meshed area, so the total applied equals the power; a given flux is multiplied by the ratio of true to meshed area. A zero area gives a flux of 0 instead of a division by zero.
const q = bc.type === 'flux' ? bc.q * s : (A > 0 ? bc.P / A : 0);The equivalent film coefficient of a radiating facet at its current temperature.
h_r = ε σ (T² + T∞²)(T + T∞) × scale, T in kelvin (°C + 273.15)
Temperatures are shifted by 273.15, which is the exact definition, and the product is formed in 64-bit from the facet's mean nodal temperature. The factored form avoids subtracting T⁴ − T∞⁴.
const updateHr = (T) => rad.forEach((i, j) => { const bc = bcs[i], Ta = bc.Tamb + KELVIN; bc.facets.forEach((f, m) => { const Tf = facetMeanT(T, f) + KELVIN; hr[j][m] = bc.eps * SIGMA_SB * (Tf * Tf + Ta * Ta) * (Tf + Ta) * (bc.scale || 1); }); });Repeats the temperature solve with updated radiation coefficients until the temperatures settle.
stop when max |T_new − T_old| ≤ 10⁻⁶ K, at most 60 passes
The stop test is an absolute temperature change of 1e-6 K over all active nodes. Each pass is a full conjugate-gradient solve started from the previous answer. The number of passes, the last change and whether it settled are returned.
for (; it < 60 && change > 1e-6; it++) {
updateHr(Tcur);
const val = Float64Array.from(Kv), rhs = Float64Array.from(fBase);
withRad(val, rhs);
const sol = csrPcg(P, val, rhs, { x0: x, tol, cancelled });
iterations += sol.iterations; residual = sol.residual; converged = converged && sol.converged;
const Tn = full(sol.x);
change = 0; for (let n = 0; n < nn; n++) if (active[n]) change = Math.max(change, Math.abs(Tn[n] - Tcur[n]));
Tcur = Tn; x = sol.x;
// i18n-ignore: progress labels never reach the screen (the workers pass on only the fraction)
if (onProgress) onProgress(Math.min(0.95, (it + 1) / 15), `Radiation: pass ${it + 1}, change ${change.toExponential(1)} K`);
}
radInfo = { iterations: it, converged: change <= 1e-6, change };Advances the temperature field in time and stores the frames shown on the time slider.
(C ÷ Δt + θK) T⁺ = (C ÷ Δt − (1 − θ)K) T + f, Δt = t ÷ steps, time of step s = s × Δt
Frame times are computed as step number × Δt, not by adding Δt repeatedly, so they do not drift. The system matrix is rebuilt each step only when radiation is present. Each step is a conjugate-gradient solve; the worst residual over the steps is returned.
const { T0, t, steps, theta = 1 } = transient, dt = t / steps, maxFrames = opts.frames || 51;
const every = Math.max(1, Math.ceil(steps / (maxFrames - 1)));
// frameSteps: the step numbers to store (0 = the start) instead of every `every`-th (thermal-history.js: one per stress solve)
const want = opts.frameSteps ? new Set(opts.frameSteps) : null, keep = (s) => (want ? want.has(s) : s % every === 0 || s === steps);
const last = want ? Math.min(steps, Math.max(0, ...want)) : steps; // no need to step past the last stored frame
let x = new Float64Array(nfree).fill(T0);
T = full(x); if (!want || want.has(0)) { frames.push(Float64Array.from(T)); times.push(0); }
let prevT = T, A = null, Ka = null, prec = null;
const C = new Float64Array(nfree); for (let n = 0; n < nn; n++) if (free[n] >= 0) C[free[n]] = Cl[n];
const diagPos = new Int32Array(nfree); for (let r = 0; r < nfree; r++) diagPos[r] = pos(P, r, r);
const tmp = new Float64Array(nfree);
radInfo = rad.length ? { iterations: steps, converged: true, linearised: 'about the previous step' } : null; // i18n-ignore: a data note, never shown
let fa = fBase;
for (let s = 1; s <= last; s++) {
if (cancelled && cancelled()) throw new Error('Cancelled');
if (!A || rad.length) {
Ka = Float64Array.from(Kv); fa = Float64Array.from(fBase);
if (rad.length) { updateHr(T); withRad(Ka, fa); }
A = new Float64Array(Ka.length);
for (let i = 0; i < Ka.length; i++) A[i] = theta * Ka[i];
for (let r = 0; r < nfree; r++) A[diagPos[r]] += C[r] / dt;
prec = ic0(P, A);
}
csrMul(P, Ka, x, tmp);
const rhs = new Float64Array(nfree);
for (let r = 0; r < nfree; r++) rhs[r] = C[r] / dt * x[r] - (1 - theta) * tmp[r] + fa[r];
const sol = csrPcg(P, A, rhs, { x0: x, tol, prec, cancelled });
iterations += sol.iterations; residual = Math.max(residual, sol.residual); converged = converged && sol.converged;
prevT = T; x = sol.x; T = full(x);
if (keep(s)) { frames.push(Float64Array.from(T)); times.push(s * dt); }The heat entering through each boundary condition, generated, stored, and what is left over.
imbalance = inflow − outflow − stored; fixed temperature: Q = (K T)ₙ − other heat into node n
All terms are 64-bit sums over nodes and facets. The imbalance subtracts totals that are nearly equal for a converged solution, so it shows the solve's residual rather than 15 digits of either total.
/** Heat into the part per boundary condition (W), generation, storage, and what is left over. */
function balanceOf(Tb, stored) {
const items = bcs.map(() => 0), KT = new Float64Array(nn), into = new Float64Array(nn); // into: all non-fixed heat entering each node
for (let e = 0; e < ne; e++) {
mesh.elemX(e, X); el.conductivity(X, k, Ke);
for (let a = 0; a < npe; a++) { let s = 0; for (let b = 0; b < npe; b++) s += Ke[a * npe + b] * Tb[mesh.conn[e * npe + b]]; KT[mesh.conn[e * npe + a]] += s; }
}
let Qgen = 0; for (let n = 0; n < nn; n++) { Qgen += genN[n]; into[n] += genN[n]; }
const facetIn = (f, hh, Tamb, i) => {
facetX(mesh, f, Xf); el.face.mass(Xf, Mf); el.face.load(Xf, Lf);
for (let a = 0; a < npf; a++) {
const na = mesh.facets.nodes[npf * f + a];
let q = hh * Tamb * Lf[a]; for (let b = 0; b < npf; b++) q -= hh * Mf[a * npf + b] * Tb[mesh.facets.nodes[npf * f + b]];
into[na] += q; items[i] += q;
}
};
bcs.forEach((bc, i) => {
if (itemLoad[i]) { for (let n = 0; n < nn; n++) if (itemLoad[i][n]) { into[n] += itemLoad[i][n]; items[i] += itemLoad[i][n]; } }
else if (bc.type === 'convection') for (const f of bc.facets) facetIn(f, bc.h * (bc.scale || 1), bc.Tamb, i);
});
rad.forEach((i, j) => bcs[i].facets.forEach((f, m) => facetIn(f, hr[j][m], bcs[i].Tamb, i)));
// fixed temperatures supply what the conduction there needs beyond the rest (storage at fixed nodes is not counted: they are held)
for (const n of fixedNodes) items[fixOwner[n]] += KT[n] - into[n];
let St = 0; if (stored) for (let n = 0; n < nn; n++) St += stored[n];
let inflow = Math.max(0, Qgen), outflow = Math.max(0, -Qgen);
for (const q of items) { if (q > 0) inflow += q; else outflow -= q; }
return { items, gen: Qgen, stored: St, inflow, outflow, imbalance: inflow - outflow - St };
}The heat flux of every element, its average at the nodes and the largest element value.
q = −k ∇T; |q| = √(q_x² + q_y² + q_z²)
Element values and nodal averages are Float64Array. The maximum is the largest element-center magnitude, taken from the 64-bit value. Nodes off the mesh get NaN.
/** Element-centre heat flux q = −k∇T (W/mm²) and its nodal average (3 per node; NaN off the mesh). */
export function heatFlux(mesh, T, k) {
const { ne, npe, nn, el = mesh.element } = mesh;
const X = new Float64Array(3 * npe), Te = new Float64Array(npe), g = new Float64Array(3), elem = new Float64Array(3 * ne), node = new Float64Array(3 * nn), cnt = new Uint16Array(nn);
let max = 0, maxAt = -1;
for (let e = 0; e < ne; e++) {
mesh.elemX(e, X);
for (let a = 0; a < npe; a++) Te[a] = T[mesh.conn[e * npe + a]];
el.gradient(X, Te, g);
for (let i = 0; i < 3; i++) elem[3 * e + i] = -k * g[i];
const m = Math.hypot(elem[3 * e], elem[3 * e + 1], elem[3 * e + 2]); if (m > max) { max = m; maxAt = e; }
for (let a = 0; a < npe; a++) { const n = mesh.conn[e * npe + a]; for (let i = 0; i < 3; i++) node[3 * n + i] += elem[3 * e + i]; cnt[n]++; }
}
for (let n = 0; n < nn; n++) for (let i = 0; i < 3; i++) node[3 * n + i] = cnt[n] ? node[3 * n + i] / cnt[n] : NaN;
return { elem, node, max, maxElem: maxAt };
}Adds every element's thermal load into the nodal force vector of a thermal-stress study.
f = Σₑ ∫ Bᵀ D ε_th dV
Element loads are added node by node into a Float64Array in element order.
/** The structural load of a temperature rise: f (3 per node) = Σ ∫ Bᵀ D ε_th dV. dT: nodal rise over the stress-free temperature (nn). */
export function thermalLoadVector(mesh, { E, nu, alpha }, dT, scale = null) {
const { ne, npe, nn, element: el } = mesh, f = new Float64Array(3 * nn);
const X = new Float64Array(3 * npe), te = new Float64Array(npe), fe = new Float64Array(3 * npe);
for (let e = 0; e < ne; e++) {
mesh.elemX(e, X);
for (let a = 0; a < npe; a++) te[a] = dT[mesh.conn[e * npe + a]];
el.thermalLoad(X, E * (scale ? scale[e] : 1), nu, alpha, te, fe);
for (let a = 0; a < npe; a++) { const n = mesh.conn[e * npe + a]; f[3 * n] += fe[3 * a]; f[3 * n + 1] += fe[3 * a + 1]; f[3 * n + 2] += fe[3 * a + 2]; }
}
return f;
}Converts the study's heat inputs to solver units and scales them for the difference between true and meshed face area.
scale = true face area ÷ meshed area (1 when either is 0)
The true area is the sum of the faces' recorded areas, the meshed area the sum over facets. Flux and film coefficient are multiplied by 1e-6. The scale guards both areas against zero.
/**
* A heat bc list from a job (face indices, SI values) → solver bcs on this mesh. exactArea: regen face areas (scale = true ÷ mesh area,
* as structural pressures do, so a curved face's staircase does not inflate the heat it exchanges).
*/
export function heatBCs(mesh, heat, exactArea = null) {
return heat.map((it) => {
const facets = mesh.facetsOn(it.faces);
if (!facets.length) throw new Error(t('study.thermal.bc_no_faces.' + it.type));
const A = facetsArea(mesh, facets), At = exactArea ? it.faces.reduce((s, fi) => s + (exactArea[fi] || 0), 0) : 0;
const scale = At > 0 && A > 0 ? At / A : 1;
const bc = { type: it.type, facets, scale, meshArea: A, area: At || A };
if (it.type === 'temperature') bc.T = it.T;
if (it.type === 'flux') bc.q = it.q * SI.flux;
if (it.type === 'power') bc.P = it.P;
if (it.type === 'convection') { bc.h = it.h * SI.h; bc.Tamb = it.Tamb; }
if (it.type === 'radiation') { bc.eps = it.eps; bc.Tamb = it.Tamb; }
return bc;
});
}The lowest and highest temperature with their locations, and the largest heat flux in W/m².
max flux in W/m² = max |q| ÷ 10⁻⁶
The minimum and maximum are read from the 64-bit temperature array. The flux is converted back to SI by dividing by 1e-6, one rounding.
/** A solved temperature field's summary: min / max with where, the flux maximum (W/m²). */
export function heatSummary(mesh, T, flux) {
let lo = Infinity, hi = -Infinity, loN = -1, hiN = -1;
for (let n = 0; n < mesh.nn; n++) { if (!mesh.active[n]) continue; if (T[n] < lo) { lo = T[n]; loN = n; } if (T[n] > hi) { hi = T[n]; hiN = n; } }
const at = (n) => (n >= 0 ? [mesh.xyz[3 * n], mesh.xyz[3 * n + 1], mesh.xyz[3 * n + 2]] : null);
return { minT: { value: lo, at: at(loN) }, maxT: { value: hi, at: at(hiN) }, maxFlux: flux ? flux.max / SI.flux : 0 };
}Conductivity and volumetric heat capacity of the material in the solver's units.
k × 10⁻³; ρc = density × cp × 10⁻⁹
Two or three 64-bit multiplications of the material's stored values. A missing conductivity or specific heat throws an error naming it.
/** The material's thermal values in solver units, or throws naming what is missing. */
export function thermalMaterial(m, { transient = false, expansion = false } = {}) {
const mat = readMaterial(m);
const miss = [!(mat.k > 0) && t('study.thermal.mat.k'), transient && !(mat.cp > 0) && t('study.thermal.mat.cp'), expansion && !(mat.alpha >= 0) && t('study.thermal.mat.alpha')].filter(Boolean);
if (miss.length) throw new Error(t('study.thermal.mat.missing', { what: miss.reduce((a, b) => t('study.thermal.mat.or', { a, b })) }));
return { mat, k: mat.k * SI.k, rc: mat.cp > 0 ? mat.density * mat.cp * SI.rc : 0 };
}The meshed volume and the mass shown in the study summary.
mass in kg = volume (mm³) × density (kg/m³) × 10⁻⁹
The volume is the voxel mesh's own volume (element count × h³, set in voxelMesh) and the mass is two multiplications. Both describe the voxel mesh, not the exact solid.
const summary = { elements: h.mesh.ne, nodes, dofs: h.res.nfree, h: grid.h, nx: grid.nx, ny: grid.ny, nz: grid.nz, oddRows: grid.oddRows, volume: h.mesh.volume, mass: h.mesh.volume * h.mat.density * 1e-9 };The arithmetic is 64-bit, but the method or the data sets the accuracy: a solver tolerance, a mesh, a step size, a datum shift, a table. The limit is stated where the result is shown.