## 2 Background  
### Finite Element Method Formulation for Transient Heat Transfer Problems  
The governing equation for transient heat transfer for a solid body with unit length, isotropic thermal conductivity \(k[\mathrm{W/mK}]\), is presented in Equation (1), where \(\rho\left[\mathrm{kg/m^{3}}\right]\) corresponds to the density, \(C_{P}[\mathrm{J/kg}\mathrm{K}]\) the specific heat, \(\dot{Q}_{V}[\mathrm{W/m^{3}}]\) is the volumetric internal heat generation rate, \(I\) corresponds to the identity matrix, and \(\
abla\) is the nabla differential operator. \[\rho C_{P}\frac{\partial T(x,y,z,t)}{\partial t}=-kI\
abla^{T}\
abla T(x,y,z,t )+\dot{Q}_{V} \tag{1}\]  
\[T(x,y,z,0)=T_{0} \tag{2}\]  
\[T(x_{5},y_{5},z_{5},t)=T_{S}=\overset{-}{T} \tag{3}\]  
\[\
abla T(x_{5},y_{5},z_{5},t)=\
abla T_{S}=\overset{-}{q} \tag{4}\]  
\[\
abla T(x_{5},y_{5},z_{5},t)=\
abla T_{S}=h_{\mathrm{U}}(T_{\infty}-T_{S})\]  
\[h_{U}=h_{conv}+h_{rad} \tag{5}\]  
\[h_{rad}=\epsilon\sigma\Big{(}T_{\infty}^{2}+T_{S}^{2}\Big{)}(T_{\infty}+T_{S})\]  
This equation allows one to find a solution for the temperature \(T(x,y,z,t)[\mathrm{K}]\) at time \(t[\mathrm{s}]\) and any position \((x,y,z)[\mathrm{m}]\) in a geometric domain. To solve Equation (1), a  
Figure 1: (**a**) Schematic representation of the fin problem, including geometrical parameters and boundary conditions. The colormap represents the steady-state solution for a fin with an insulated tip subjected to room conditions with a predefined wall temperature. (**b**) Transient Thermal ZFEM methodology for sensitivity analysis. time-dependent temperature condition, e.g., the initial temperature (\(T_{0}[\text{K}]\)) in Equation (2), and surface boundary conditions are required. Surface conditions are denoted by the sub-index \(S\) in Equations (3)-(5). The surface boundary conditions are predefined temperatures (\(T[\text{K}]\), in Equation (3)) and predefined heat fluxes (\(\overset{\rightarrow}{q}\left[W/\text{m}^{2}\right]\), in Equation (4)). For example, Figure 1a represents a fin with a predefined wall temperature (\(T=T_{W}\)), and insulation at the tip (\(\overset{\rightarrow}{q}=0\) W/m\({}^{2}\)). Heat fluxes can be position-dependent, time-dependent, or surface temperature-dependent, as in convection and radiation conditions (Equation (5)). Here, \(\epsilon\) is the material's surface emissivity, \(\sigma=5.67\times 10^{-8}\) W/m\({}^{2}\)K\({}^{4}\) is the Stefan-Boltzmann constant. If there are multiple boundary conditions on the same surface, superposition is used for the analysis. For instance, if radiation and convection boundary conditions are applied over the same surface, the global heat transfer coefficient \(h_{U}\left[W/\text{m}^{2}\text{K}\right]\) is used to represent the contribution from both heat transfer mechanisms, as shown in Equation (5). Here, \(T_{\infty}[\text{K}]\) is the room temperature, \(h_{conv}\left[W/\text{m}^{2}\text{K}\right]\) is the convection heat transfer coefficient, and \(h_{rad}\) is the radiation heat transfer. When \(T_{0}\) and \(T_{\infty}\) are similar, \(h_{rad}\left[W/\text{m}^{2}\text{K}\right]\) can be treated as a constant by replacing the surface temperature (\(T_{S}\)) in Equation (5) by the average temperature among \(T_{0}\) and \(T_{\infty}\)[29; 32]. Using the finite element method, the representation of the nodal temperatures, \(T^{\epsilon}\), in the element is specified by the shape functions, \(\Phi^{\epsilon}\), given by:  
\[T(x,y,z,t)=\Phi^{\epsilon}T^{\epsilon}=\sum_{j=1}^{n}T_{j}^{\epsilon}\cdot\Phi _{j}^{\epsilon}(x,y,z,t) \tag{6}\]  
where \(n\) represents the number of nodes of the element [82]. Combining Equations (1)-(6), we obtain Equation (7), which represents the transient heat transfer FEM formulation for any time step (\(t\)s).