One-dimensional finite-difference formulation

This page is a draft public technical description of the current FDM 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-difference method (FDM) solves the transient one-dimensional thermal problem in two workpiece regions, \(\sideR{R}\) and \(\sideS{S}\), using a cell-centred spatial support. The local coordinate in each region points from the welding interface toward the distant end of that workpiece.

The spatial method is independent of the temporal integrator and of the calculation engine. The active GNU Octave implementation uses the project-wide time-integration contract, and an independent C++ transcription exists for cross-engine work. This page describes the numerical formulation, not the Octave or C++ software architecture.

The current production FDM operator identifies itself as:

fdm-1d-cell-centred-positive-logkappa-flux-form-quadratic-interface
revision 2.1.0

when the default quadratic-second-order interface closure is active.

2. Continuous thermal problem

For each workpiece role \(\xi\in\{\sideR{R},\sideS{S}\}\), the current thermal model is written in one-dimensional flux form as

\rho_\xi(T_\xi)c_{p,\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).

Here \(\rho\), \(c_p\) and \(\kappa\) are resolved constitutive properties, while \(q'''_v\) represents resolved distributed source/sink contributions. Depending on the active process stage, these may include volumetric electrical heating and reduced lateral heat exchange. Physical distant-end and welding- interface conditions are defined independently at the model layer and are only then converted into FDM boundary fluxes.

The thermal problem therefore separates:

physical/model laws
        -> resolved material/process data
        -> FDM spatial closure
        -> semi-discrete ODE system
        -> selected time integrator

3. Cell-centred spatial support

For a region of physical 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 stored temperature \(T_{\xi,\fdx{i}}\) is therefore associated with the centre of a physical cell. In particular, the first stored temperature is half a cell from the welding interface and the last stored temperature is half a cell from the distant physical surface.

This distinction is essential: a stored endpoint value is not automatically identical to a physical surface temperature. The half-cell distance is handled by the relevant numerical boundary/interface closure.

4. Conservative shared-face form

The active production operator constructs a single shared sequence of conductive face quantities

\Phi_{\xi,\fdx{1/2}}, \Phi_{\xi,\fdx{3/2}}, \ldots, \Phi_{\xi,\fdx{N_\xi+1/2}}

and defines the cell-centred conductive contribution by the exact discrete face difference

\mathcal D_{h,\xi}(T)_{\fdx{i}} = \frac{ \Phi_{\xi,\fdx{i+1/2}} -\Phi_{\xi,\fdx{i-1/2}} }{\fdx{\Delta x_\xi}}.

The first and last entries of the same face sequence are the physical boundary fluxes supplied by the welding-interface and distant-end closures. They are not post-processing values added after an interior operator has already been formed.

This is also the invariant exercised by the permanent conservation test: D = diff(face_flux)/dx exactly in the implementation.

5. Interior face-gradient approximation

For an ordinary internal face between cells \(\fdx{i}\) and \(\fdx{i+1}\), the temperature gradient is evaluated using the fourth-order centred face formula

\left. \frac{\partial T_\xi}{\partial \fdx{x}} \right|_{\fdx{i+1/2}} \approx \frac{ T_{\xi,\fdx{i-1}} -27T_{\xi,\fdx{i}} +27T_{\xi,\fdx{i+1}} -T_{\xi,\fdx{i+2}} }{24\fdx{\Delta x_\xi}}.

Near the first and last internal faces, the implementation uses the approved one-sided counterparts so that no exterior ghost temperature is introduced merely to complete the stencil. For the first internal face the derivative is

\left. \frac{\partial T_\xi}{\partial \fdx{x}} \right|_{\fdx{3/2}} \approx \frac{1}{\fdx{\Delta x_\xi}} \left( -\frac{11}{12}T_{\xi,1} +\frac{17}{24}T_{\xi,2} +\frac{3}{8}T_{\xi,3} -\frac{5}{24}T_{\xi,4} +\frac{1}{24}T_{\xi,5} \right),

with the right-side expression obtained by the corresponding mirrored stencil.

6. Positive conductivity reconstruction

The current implementation does not interpolate \(\kappa\) directly at an internal face. Instead, it reconstructs \(\log\kappa\) and exponentiates the result. For an ordinary internal face,

\log \kappa_{\xi,\fdx{i+1/2}} \approx \frac{ -\log\kappa_{\xi,\fdx{i-1}} +9\log\kappa_{\xi,\fdx{i}} +9\log\kappa_{\xi,\fdx{i+1}} -\log\kappa_{\xi,\fdx{i+2}} }{16},

followed by

\kappa_{\xi,\fdx{i+1/2}} = \exp\!\left( \log \kappa_{\xi,\fdx{i+1/2}} \right).

This is a current implementation detail that supersedes any older explanatory formula that interpolates conductivity itself. For finite positive nodal conductivities, the reconstruction is designed to preserve positivity unless floating-point range is exceeded, in which case the runtime fails explicitly.

The conductive face quantity is then

\Phi_{\xi,\fdx{i+1/2}} = \kappa_{\xi,\fdx{i+1/2}} \left. \frac{\partial T_\xi}{\partial \fdx{x}} \right|_{\fdx{i+1/2}}.

7. Welding-interface closure

The physical interface law is defined at the model layer. The current production cell-centred closure converts that law into the default quadratic-second-order interface flux. Defining \(\sideinterface{q''_0}>0\) as heat transferred from \(\sideR{R}\) to \(\sideS{S}\),

\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) }.

The corresponding reconstructed interface-side temperatures are

\sideR{T_R^I} = \frac{9\sideR{T_{R,1}}-\sideR{T_{R,2}}}{8} - \frac{3\fdx{\Delta x_R}\sideinterface{q''_0}} {8\kappa_{R,1}},

\sideS{T_S^I} = \frac{9\sideS{T_{S,1}}-\sideS{T_{S,2}}}{8} + \frac{3\fdx{\Delta x_S}\sideinterface{q''_0}} {8\kappa_{S,1}}.

The same thermal transfer is applied to the two regions with opposite signs. There is no fictitious interface temperature used as a ghost node to complete the FDM stencil.

Interface electrical dissipation, when present, is a separate source term. The current runtime incorporates the allocated interface electrical power into the region boundary-flux bookkeeping rather than redefining the thermal contact law.

A former two-point first-order interface closure remains available only for legacy reproduction/regression. It is not the current production default.

8. Distant-end Robin closure

The active physical distant-end condition is Newton/Robin exchange,

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

where \(T_{s,\xi}\) is the physical surface temperature. Because the final stored temperature is located half a cell inside that surface, the active cell-centred closure eliminates the un-stored surface value and uses

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 adiabatic limit is obtained for \(h_{end,\xi}=0\). The endpoint degree of freedom remains part of the evolving thermal state; it is not reset after each time step to a prescribed ambient temperature.

9. Semi-discrete system

After spatial closure, each regional nodal balance can be written as

\rho_{\xi,\fdx{i}}c_{p,\xi,\fdx{i}} \frac{dT_{\xi,\fdx{i}}}{dt} = \mathcal D_{h,\xi}(T)_{\fdx{i}} +Q_{\xi,\fdx{i}}(T,t),

where \(Q\) denotes the resolved distributed source/sink contribution supplied by the active physical/process model.

Equivalently,

\frac{d\mathbf T_\xi}{dt} = \mathbf F_\xi\!\left( t,\sideR{\mathbf T_R},\sideS{\mathbf T_S};\mathcal P(t) \right).

The FDM solver supplies this semi-discrete right-hand side to the shared project-level temporal integrator. The spatial method itself does not define a fixed time-stepping algorithm.

10. Current temporal configuration

The temporal architecture is method-neutral. At present, the principal explicit configuration is Adams--Bashforth of order two (AB2), with Forward Euler as the configured starter/restart method.

For a right-hand side \(\mathbf F^n\), Forward Euler is

\mathbf T^{n+1} = \mathbf T^n +\Delta t\,\mathbf F^n,

while AB2 is

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

This is a current reproducible configuration, not a definition of FDM. Changing the time integrator does not redefine the physical stage-transfer law or the spatial finite-difference operator.

At a process event, physical state and scientific process history follow the model-level transfer contract. Private multistep integrator memory may continue, restart or be reconstructed only according to the temporal-integrator contract.

11. Stability reference

For the constant-coefficient fourth-order five-point diffusion operator, the most negative Fourier eigenvalue used by the current runtime reference is

\lambda_{min} = -\frac{16\alpha}{3\fdx{\Delta x}^2}, \qquad \alpha=\frac{\kappa}{\rho c_p}.

For AB2, whose current negative-real stability radius is one in this normalisation, this gives the reference Fourier-number limit

\mathrm{Fo} = \frac{\alpha\Delta t}{\fdx{\Delta x}^2} \leq \frac{3}{16}.

The runtime scales the stored spatial reference by the selected temporal integrator's negative-real stability radius. This bound is a practical constant-property reference. Temperature-dependent coefficients, physical boundary closures, interface coupling and event-driven changes still require refinement and verification; exceeding the configured reference is treated as an error rather than silently advancing.

12. Accuracy statement: local high order does not imply global fourth order

The current FDM implementation retains fourth-order centred and one-sided ingredients for face gradients and related internal reconstruction. However, the project's manufactured-solution test for the active variable- conductivity conservative cell-centred operator explicitly accepts an observed global order of approximately two, with an acceptance floor of 1.9.

Accordingly, the current public claim is:

The solver uses high-order/fourth-order spatial ingredients in its internal flux construction, while the approved observed global convergence order of the present variable-conductivity cell-centred operator is approximately second order.

The complete solver must not be presented as globally fourth-order unless a later accepted verification record changes this status.

13. Verification evidence represented by the current test suite

The current source tree contains dedicated verification logic for, among other items:

exact shared-face difference structure and boundary-flux conservation;
constant-conductivity reconstruction and positivity of reconstructed face

conductivity;

exact equilibrium of a linear temperature field with constant conductivity;
manufactured-solution spatial refinement for the active variable-

conductivity operator;

cell-centred Robin closure and steady Robin refinement;
interface-source convergence and interface-source energy conservation;
positive-conductivity production gates;
temporal discretisation and restart behaviour;
compatibility/equivalence of historical wrappers with the canonical API.

The existence of these tests shows what the project checks. A future public verified or PASS badge must still be tied to an accepted run/CI record for the actual revision published by MPX rather than inferred from the mere presence of test files.

14. Current assumptions and limits

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

the complete multiphysics platform is one-dimensional;
every future FDM formulation must use the same cell-centred support;
AB2/Forward Euler is the only permitted temporal scheme;
the current endpoint coefficient represents a fully resolved physical

jaw/electrode model;

fixed-grid FDM and FEM results must be identical;
numerical verification constitutes experimental validation.

The current model keeps the physical law, numerical method and software engine as separate architectural concerns so that any of these may evolve without silently redefining the others.

15. Scientific notation used on this page

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

\\fdx{...} — finite-difference/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 must remain sufficient to identify the same roles when colour is unavailable.