One-dimensional finite-element formulation

This page is a public technical preview of the current FEM spatial formulation. It documents the present project state; it is not a republication of the 2018 dissertation and it does not make the dissertation the normative source of the model.

1. Scope and current status

The current finite-element method (FEM) solves the same transient one-dimensional thermal problem used by the FDM formulation in the two workpiece regions \(\sideR{R}\) and \(\sideS{S}\). It changes the spatial semi-discretisation while retaining the same model-level material laws, volumetric sources, welding-interface law, distant-end boundary condition and shared temporal-integration architecture.

The active FEM uses a cell-centred support with two-node linear Galerkin (P1) elements connecting consecutive thermal degrees of freedom. The active GNU Octave implementation evaluates the spatial residual and passes the resulting semi-discrete right-hand side to the project-wide time integrator. An independent C++ implementation also exists for resolved-case cross-engine work.

With the default quadratic-second-order interface closure, the current FEM operator identifies itself as:

fem-1d-linear-galerkin-lumped-mass-quadratic-interface
revision 1.1.0

This page describes the numerical formulation, not the Octave or C++ software architecture.

2. Continuous thermal problem

For either workpiece role \(\xi\in\{\sideR{R},\sideS{S}\}\), define \(\beta_\xi=\rho_\xi c_{p,\xi}\). The continuous one-dimensional thermal balance is

\beta_\xi(T_\xi) \frac{\partial T_\xi}{\partial t} - \frac{\partial}{\partial\fdx{x}} \left( \kappa_\xi(T_\xi) \frac{\partial T_\xi}{\partial\fdx{x}} \right) = q'''_{v,\xi}(T_\xi,t).

The formulation is written per unit cross-sectional area. Multiplying every term by a uniform represented area gives the equivalent total-power form.

The physical boundary and interface laws are defined outside FEM. The finite- element method is responsible only for constructing a spatially discrete representation consistent with those laws.

3. Cell-centred support and Galerkin skeleton

For a physical region of length \(L_\xi\) divided into \(\fdx{N_\xi}\) equal cells,

\fdx{\Delta x_\xi}=\frac{L_\xi}{\fdx{N_\xi}}, \qquad \fdx{x}_{\xi,\fdx{i}} = \left(\fdx{i}-\frac{1}{2}\right)\fdx{\Delta x_\xi}, \qquad \fdx{i}=1,\ldots,\fdx{N_\xi}.

The \(\fdx{N_\xi}\) thermal degrees of freedom lie at cell centres. The \(\fdx{N_\xi}-1\) linear elements therefore span only the centre-to-centre Galerkin skeleton

\widetilde\Omega_\xi = \left[ \frac{\fdx{\Delta x_\xi}}{2}, L_\xi-\frac{\fdx{\Delta x_\xi}}{2} \right].

The physical domain still includes the two omitted half-cells of width \(\fdx{\Delta x_\xi}/2\). They are not discarded and the P1 interpolation is not silently extended to the physical faces. Their capacity and volumetric source measures are completed explicitly, while their conductive transfer is handled by the welding-interface and distant-end flux closures.

4. Regional weak form

Let \(v_\xi\) be an admissible test function on the centre-to-centre skeleton. Multiplying the strong form by \(v_\xi\), integrating over \(\widetilde\Omega_\xi\), and integrating the conduction term by parts gives

\int_{\widetilde\Omega_\xi} \beta_\xi v_\xi \frac{\partial T_\xi}{\partial t} \,\mathrm d\fdx{x} + \int_{\widetilde\Omega_\xi} \kappa_\xi \frac{\partial v_\xi}{\partial\fdx{x}} \frac{\partial T_\xi}{\partial\fdx{x}} \,\mathrm d\fdx{x} = \int_{\widetilde\Omega_\xi} v_\xi q'''_{v,\xi} \,\mathrm d\fdx{x} + \left[ v_\xi\kappa_\xi \frac{\partial T_\xi}{\partial\fdx{x}} \right]_{\fdx{x}_{\xi,1}}^{\fdx{x}_{\xi,N_\xi}}.

The endpoint term is not closed by pretending that the first and last FEM unknowns are physical surface temperatures. The welding side receives the cell-centred contact/interface flux; the distant side receives the cell-centred Robin flux.

5. Two-node linear Galerkin element

For an element \(e=[\fdx{x}_a,\fdx{x}_b]\) of length \(\fdx{h_e}=\fdx{x}_b-\fdx{x}_a\), introduce the reference coordinate \(\zeta\in[-1,1]\). The P1 shape functions are

N_1(\zeta)=\frac{1-\zeta}{2}, \qquad N_2(\zeta)=\frac{1+\zeta}{2},

with isoparametric map and Jacobian

\fdx{x}(\zeta) =N_1\fdx{x}_a+N_2\fdx{x}_b, \qquad J_e=\frac{\fdx{h_e}}{2}.

Defining

\mathbf N= \begin{bmatrix}N_1&N_2\end{bmatrix}, \qquad \mathbf B = \frac{\partial\mathbf N}{\partial\fdx{x}} = \begin{bmatrix} -1/\fdx{h_e} & 1/\fdx{h_e} \end{bmatrix},

the trial and test fields are approximated by

T_{\xi,h}=\mathbf N\mathbf T_{\xi,e}, \qquad v_{\xi,h}=\mathbf N\,\delta\mathbf T_{\xi,e}.

No extra common thermal degree of freedom is introduced at the welding interface.

6. Element capacity, conduction and source terms

The element thermal-capacity matrix is

\mathbf M_{\xi,e}(\mathbf T_{\xi,e}) = \int_{\Omega_e} \rho_\xi(T_{\xi,h})c_{p,\xi}(T_{\xi,h}) \mathbf N^{\mathsf T}\mathbf N \,\mathrm d\fdx{x},

the conductive matrix is

\mathbf K_{\kappa,\xi,e}(\mathbf T_{\xi,e}) = \int_{\Omega_e} \kappa_\xi(T_{\xi,h}) \mathbf B^{\mathsf T}\mathbf B \,\mathrm d\fdx{x},

and the consistent volumetric-load vector is

\mathbf f_{v,\xi,e}(\mathbf T_{\xi,e},t) = \int_{\Omega_e} \mathbf N^{\mathsf T} q'''_{v,\xi}(T_{\xi,h},t) \,\mathrm d\fdx{x}.

For frozen coefficients that are constant over one element, \(\beta_{\xi,e}=\rho_{\xi,e}c_{p,\xi,e}\), the familiar reference forms are

\mathbf M_{\xi,e} = \frac{\beta_{\xi,e}\fdx{h_e}}{6} \begin{bmatrix} 2&1\\ 1&2 \end{bmatrix},

\mathbf K_{\kappa,\xi,e} = \frac{\kappa_{\xi,e}}{\fdx{h_e}} \begin{bmatrix} 1&-1\\ -1& 1 \end{bmatrix},

and, for constant \(q'''_{v,\xi,e}\),

\mathbf f_{v,\xi,e} = \frac{q'''_{v,\xi,e}\fdx{h_e}}{2} \begin{bmatrix}1\\1\end{bmatrix}.

These are reference limits. The active nonlinear thermal implementation reevaluates temperature-dependent material coefficients at the current thermal state rather than replacing them with constants.

7. Two-point Gauss quadrature and current assembly

The active Octave assembly uses a two-point Gauss rule on every P1 element. For \(\zeta=\pm1/\sqrt{3}\), the two shape-function weight pairs are

(a,b) = \left( \frac{1+1/\sqrt{3}}{2}, \frac{1-1/\sqrt{3}}{2} \right)

and \((b,a)\) at the opposite Gauss point. The code evaluates \(\rho c_p\), \(\kappa\), and the interpolated volumetric source at these states.

The current implementation vectorises this assembly across all elements. This is an algebraic re-expression of the same P1/two-point-Gauss formulation: it does not change the trial space, quadrature, lumping rule, boundary completion or constitutive evaluation points.

8. Row-sum mass lumping

The Galerkin capacity matrix is retained in the mathematical derivation, but the active explicit solver uses row-sum mass lumping,

\mathbf M_{L,\xi,e} = \operatorname{Diag}\!\left( \mathbf M_{\xi,e}\mathbf 1 \right).

For constant coefficients on a linear element,

\mathbf M_{L,\xi,e} = \frac{\beta_{\xi,e}\fdx{h_e}}{2} \begin{bmatrix} 1&0\\ 0&1 \end{bmatrix}.

Mass lumping makes the assembled capacity diagonal, so the explicit right-hand side can divide componentwise by the lumped capacity rather than solve a global capacity system at every evaluation. Diagonal capacity does not mean constant capacity: the current \(\rho c_p\) values still depend on temperature.

9. Half-cell completion of capacity and volumetric source

The centre-to-centre P1 skeleton covers only \((\fdx{N_\xi}-1)\fdx{\Delta x_\xi}\), while the physical region has length \(\fdx{N_\xi}\fdx{\Delta x_\xi}\). The two missing measures are therefore added directly to the endpoint rows.

For each boundary half-cell \(b\in\{1,N_\xi\}\), with \(w_{\xi,b}=\fdx{\Delta x_\xi}/2\), the active completion is

\Delta M_{L,\xi,b} = \beta_\xi(T_{\xi,b})w_{\xi,b},

\Delta f_{v,\xi,b} = q'''_{v,\xi}(T_{\xi,b},t)w_{\xi,b}.

For constant capacity and source this restores the exact represented domain measure,

\mathbf 1^{\mathsf T}\mathbf M_{L,\xi}\mathbf 1 = \beta_\xi L_\xi, \qquad \mathbf 1^{\mathsf T}\mathbf f_{v,\xi} = q'''_{v,\xi}L_\xi.

The dedicated domain-integral test checks these identities. Conductive stiffness is not completed by extending the P1 skeleton through the half- cells; conductive transfer there is represented by the support-aware physical flux closures.

10. Global assembly and residual form

Let \(\mathcal A\) denote the usual finite-element assembly operator. For each region,

\mathbf M_{L,\xi} =\mathcal A_e\left(\mathbf M_{L,\xi,e}\right), \qquad \mathbf K_{\kappa,\xi} =\mathcal A_e\left(\mathbf K_{\kappa,\xi,e}\right), \qquad \mathbf f_{v,\xi} =\mathcal A_e\left(\mathbf f_{v,\xi,e}\right),

with the endpoint half-cell capacity/source contributions included in \(\mathbf M_L\) and \(\mathbf f_v\).

The active implementation is most directly expressed as the regional residual balance

\mathbf M_L(\mathbf T)\dot{\mathbf T} = \mathbf f_v(\mathbf T,t) -\mathbf K_\kappa(\mathbf T)\mathbf T +\mathbf b_{end}(\mathbf T) +\mathbf b_{\Gamma}(\mathbf T) +\mathbf f_{0,J}(\mathbf T,t).

Here \(\mathbf b_{end}\) represents distant-end heat exchange, \(\mathbf b_\Gamma\) the equal-and-opposite thermal contact exchange, and \(\mathbf f_{0,J}\) any zero-thickness electrical-interface power assigned to the adjacent degrees of freedom.

11. Distant-end Robin closure

The physical distant-end law remains

q''_{out,\xi} = h_{end,\xi} \left(T_{s,\xi}-T_{\infty,end,\xi}\right).

Because the last FEM unknown is a cell-centred temperature located half a cell inside the physical boundary, the active runtime uses the same support-aware closure as the other cell-centred spatial methods,

q''_{out,\xi} = \frac{ T_{\xi,\fdx{N_\xi}}-T_{\infty,end,\xi} }{ \dfrac{\fdx{\Delta x_\xi}}{2\kappa_{\xi,\fdx{N_\xi}}} +\dfrac{1}{h_{end,\xi}} }.

The outward flux is subtracted from the corresponding endpoint residual. The classical nodal insertion K_NN += h and f_N += h*T_inf is a useful reference when the terminal FEM unknown lies on the physical surface, but it is not the active cell-centred boundary closure.

12. Welding-interface coupling

The physical R/S contact law is shared across spatial methods. With \(\sideinterface{q''_0}>0\) defined as heat transferred from \(\sideR{R}\) to \(\sideS{S}\), the active runtime evaluates the common cell-centred interface state and applies the same heat flux with opposite residual signs,

r_{R,1}\leftarrow r_{R,1}-\sideinterface{q''_0}, \qquad r_{S,1}\leftarrow r_{S,1}+\sideinterface{q''_0}.

For the default quadratic cell-centred closure,

\sideinterface{q''_0} = \frac{ 9\sideR{T_{R,1}}-\sideR{T_{R,2}} -9\sideS{T_{S,1}}+\sideS{T_{S,2}} }{ 8\sideinterface{R''_{th,0}} +3\left( \frac{\fdx{\Delta x_R}}{\kappa_{R,1}} + \frac{\fdx{\Delta x_S}}{\kappa_{S,1}} \right) }.

No additional FEM degree of freedom is introduced at the welding interface. The historical scalar interface temperature is reconstructed from the updated one-sided states after the time step.

Electrical power dissipated at a zero-thickness interface remains a separate source. The allocated power is added directly to the two interface-adjacent residual rows per unit represented area; it is not hidden inside an arbitrary element thickness.

13. Explicit semi-discrete right-hand side and time integration

After assembly and physical flux coupling, the active FEM supplies

\mathbf F_{\mathrm{FEM}}(\mathbf T,t) = \left[\mathbf M_L(\mathbf T)\right]^{-1} \mathbf r_{\mathrm{FEM}}(\mathbf T,t)

to the shared temporal layer. Because \(\mathbf M_L\) is diagonal, the runtime implements this inversion componentwise.

The temporal architecture is independent of FEM. The current principal explicit configuration is AB2 with a Forward-Euler starter/restart path,

\mathbf T^{n+1} = \mathbf T^n +\Delta t \left( \frac{3}{2}\mathbf F^n -\frac{1}{2}\mathbf F^{n-1} \right),

with Forward Euler used when compatible multistep history is unavailable. This is a current reproducible configuration, not a definition of FEM.

At process events the physical state follows the model-level transfer contract, while private integrator memory is independently continued, restarted or reconstructed according to the time-integrator contract.

14. Constant-property stability reference

For an infinite uniform mesh with constant \(\alpha=\kappa/(\rho c_p)\), P1 linear elements, row-sum lumped mass and no sources, the interior semi-discrete equation is

\beta\fdx{h}\,\dot T_i + \frac{\kappa}{\fdx{h}} \left(-T_{i-1}+2T_i-T_{i+1}\right) =0.

The corresponding Fourier eigenvalue is

\lambda_L(\theta) = -\frac{4\alpha}{\fdx{h}^{2}} \sin^2\!\left(\frac{\theta}{2}\right).

For the current AB2 negative-real-axis stability interval, this gives the reference

\mathrm{Fo}_{\mathrm{FEM},L} = \frac{\alpha\Delta t}{\fdx{h}^{2}} \leq \frac{1}{4}.

The runtime uses this value as a constant-property interior reference scaled by the selected temporal integrator's negative-real stability radius. It is not a nonlinear stability proof for the complete finite-domain problem with state-dependent coefficients, Robin exchange, interface coupling, electrical sources and process events.

15. Verification evidence and publication caution

The current source tree contains dedicated FEM checks for, among other items:

physical-domain capacity and source integrals under the cell-centred half-cell

completion;

equivalence of the vectorised assembly to the reference element-by-element

formulation;

the constant-property lumped-mass stability reference

\(\mathrm{Fo}_{crit}\approx1/4\);

Robin sign conventions and manufactured/reference boundary tests;
wrapper/canonical-path equivalence;
FDM/FEM cross-method diagnostics.

The FDM/FEM comparison is deliberately diagnostic: neither method is treated as an exact reference, and identical values on a fixed coarse mesh are not a correctness requirement. The meaningful scientific question is whether independent discretisations approach the same continuous problem under refinement while preserving the same physical laws and process inputs.

A repository validation note also records that, in the environment where that note was written, native GNU Octave execution still had to be performed in the author's Octave environment. Therefore this page does not convert the mere existence of test files into a public validated or PASS badge. Any such badge must be tied to an accepted run, CI record or explicitly reviewed result for the published revision.

16. Current assumptions and limits

This page presently documents only the active one-dimensional thermal FEM path. It does not claim that:

the complete Multiphysics platform is one-dimensional;
every future FEM formulation must use P1 elements or cell-centred support;
mass lumping is a universal FEM requirement;
AB2/Forward Euler is the only permitted temporal configuration;
a nodal Robin matrix insertion is the active cell-centred boundary model;
numerical cross-method agreement constitutes experimental validation;
the C++ development engine has production status.

The current architecture intentionally separates physical model, spatial method, temporal integration and calculation engine so that each may evolve without silently redefining the others.

17. Scientific notation used on this page

The technical website preserves the semantic equation colours inherited from the refactored dissertation:

\\fdx{...} — discrete spatial notation;
\\sideR{...} — workpiece role R;
\\sideS{...} — workpiece role S;
\\sideinterface{...} — welding-interface quantities.

These colours encode scientific meaning and are independent of the decorative website palette. Text labels and mathematical notation remain authoritative when colour is unavailable.