[go: up one dir, main page]

Next Article in Journal
Synthesis and Characterisation of Alginate-Based Capsules Containing Waste Cooking Oil for Asphalt Self-Healing
Previous Article in Journal
Predicting Children with ADHD Using Behavioral Activity: A Machine Learning Analysis
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Sensitivity Analysis for Transient Thermal Problems Using the Complex-Variable Finite Element Method

by
Juan-Sebastian Rincon-Tabares
1,
Juan C. Velasquez-Gonzalez
1,
Daniel Ramirez-Tamayo
1,2,
Arturo Montoya
1,3,
Harry Millwater
1 and
David Restrepo
1,*
1
Department of Mechanical Engineering, University of Texas at San Antonio, San Antonio, TX 78249, USA
2
Advanced Computing, Mathematics and Data Division, Pacific Northwest National Laboratory, Richland, WA 99354, USA
3
Department of Civil and Environmental Engineering, University of Texas at San Antonio, San Antonio, TX 78249, USA
*
Author to whom correspondence should be addressed.
Appl. Sci. 2022, 12(5), 2738; https://doi.org/10.3390/app12052738
Submission received: 9 February 2022 / Revised: 2 March 2022 / Accepted: 2 March 2022 / Published: 7 March 2022
(This article belongs to the Section Applied Thermal Engineering)
Figure 1
<p>(<b>a</b>) 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>b</b>) Transient Thermal ZFEM methodology for sensitivity analysis.</p> ">
Figure 2
<p>(<b>a</b>) Mesh perturbation function <math display="inline"><semantics> <mrow> <mi>P</mi> <mrow> <mo>(</mo> <mi>x</mi> <mo>)</mo> </mrow> </mrow> </semantics></math> for nodal coordinates required for geometrical sensitivities using Transient Thermal ZFEM, where <math display="inline"><semantics> <mi>p</mi> </semantics></math> represents a fraction of the mesh that was perturbed. (<b>b</b>) Duplicated node approach for the CR representation for a 2D finite element. (<b>c</b>) Detail of the local coordinate system and quadrature points for the 4 real 4 imaginary nodes 2D element.</p> ">
Figure 3
<p>Flowchart of Transient Thermal ZFEM and Abaqus User Element Subroutine. All complex operations are performed inside the UEL and the <span class="html-italic">A</span> matrix and the <span class="html-italic">RHS</span> vector are decomposed into Cauchy-Riemann form before exiting the UEL.</p> ">
Figure 4
<p>(<b>a</b>) Convergence analysis for the heat transfer coefficient per unit of thickness where the optimal value that minimizes the NRMSE is <math display="inline"><semantics> <mrow> <mn>735</mn> <mrow> <mo>(</mo> <mrow> <mi mathvariant="normal">W</mi> <mo>/</mo> <mrow> <mo>(</mo> <mrow> <msup> <mi mathvariant="normal">m</mi> <mn>2</mn> </msup> <mi mathvariant="normal">K</mi> </mrow> <mo>)</mo> </mrow> </mrow> <mo>)</mo> </mrow> <mo>/</mo> <mi mathvariant="normal">m</mi> </mrow> </semantics></math>. (<b>b</b>) The dimensionless temperature history of fin tip: solid black line, Equation (23), yellow star, the real part of Transient Thermal ZFEM solution, blue square, Abaqus built-in element solution, thin red dashed line, steady-state solution (Equation (24)).</p> ">
Figure 5
<p>(<b>a</b>) Step size convergence analysis comparing Transient Thermal ZFEM and FD for the UEL and the selected analyzed variables. (<b>b</b>) Convergence study for mesh perturbation fraction <math display="inline"><semantics> <mi>λ</mi> </semantics></math>. The insets show <math display="inline"><semantics> <mrow> <mi>P</mi> <mrow> <mo>(</mo> <mi>x</mi> <mo>)</mo> </mrow> </mrow> </semantics></math> for the rule of thumb of the literature and <math display="inline"><semantics> <mrow> <mi>λ</mi> <mo>=</mo> <mn>1</mn> </mrow> </semantics></math> in the top left and bottom right corners, respectively.</p> ">
Figure 6
<p>Dimensionless temperature sensitivity history at the fin tip to (<b>a</b>) thermal conductivity; (<b>b</b>) room temperature; (<b>c</b>) fin length. Solid black line, analytical sensitivity using analytical differentiation in Equation (23), yellow star, sensitivity using Transient Thermal ZFEM, blue square, sensitivity using 1% FD and Abaqus built-in element solutions, thin red dashed line, sensitivity using analytical differentiation for steady-state using Equation (24).</p> ">
Figure 7
<p>The heating rate in time for fin tip: solid black line, rate using analytical differentiation in Equation (23), yellow star, rate using Transient Thermal ZFEM, blue square, FD using Abaqus built-in element solution.</p> ">
Figure 8
<p>(<b>a</b>) Temperature-dependent specific heat of Ti64. Experimental data obtained from [<a href="#B90-applsci-12-02738" class="html-bibr">90</a>]. The inset shows the original and the perturbed property in the <math display="inline"><semantics> <mrow> <mi>T</mi> <mo>−</mo> <msub> <mi>C</mi> <mi>P</mi> </msub> <mo>−</mo> <mi>i</mi> </mrow> </semantics></math> space; (<b>b</b>) dimensionless temperature history of fin tip for temperature-dependent specific heat; (<b>c</b>) sensitivity of <math display="inline"><semantics> <mi>θ</mi> </semantics></math> to specific temperature-dependent specific heat. Yellow star, sensitivity using Transient Thermal ZFEM, blue square, sensitivity using <math display="inline"><semantics> <mrow> <mn>1</mn> <mo>%</mo> </mrow> </semantics></math> FD and Abaqus built-in element solution, black dashed line, steady-state solution.</p> ">
Figure 9
<p>(<b>a</b>) Time history of fin <math display="inline"><semantics> <mi>ε</mi> </semantics></math>, (<b>b</b>) time history of fin <math display="inline"><semantics> <mi>η</mi> </semantics></math>. Black line, Equation (13), yellow star, Transient Thermal ZFEM solution, blue square, Abaqus built-in element solution, black thin dashed line, black dashed line, the steady-state solution.</p> ">
Figure 10
<p>Thermal effectiveness sensitivity history in fin tip to (<b>a</b>) thermal conductivity; (<b>b</b>) room temperature; (<b>c</b>) fin length. Solid black line, analytical sensitivity using analytical differentiation in Equation (26), yellow star, sensitivity using Transient Thermal ZFEM, blue square, sensitivity using <math display="inline"><semantics> <mrow> <mn>1</mn> <mo>%</mo> </mrow> </semantics></math> FD and Abaqus built-in element solutions, thin red dashed line, sensitivity using analytical differentiation for steady-state using Equation (26). (<b>d</b>) Steady-state influence of fin parameters in <math display="inline"><semantics> <mi>ε</mi> </semantics></math> based on relative sensitivities. In the pie chart, the missing parameters (<math display="inline"><semantics> <mrow> <mi>b</mi> <mo>,</mo> <mi>ρ</mi> <mo>,</mo> <msub> <mi>C</mi> <mi>P</mi> </msub> <mo>,</mo> <msub> <mi>T</mi> <mo>∞</mo> </msub> <mo>,</mo> <msub> <mi>T</mi> <mn>0</mn> </msub> </mrow> </semantics></math>) have an influence of less than <math display="inline"><semantics> <mrow> <mn>1</mn> <mo>%</mo> </mrow> </semantics></math> compared to the sum of all relative sensitivities. Moreover, the <math display="inline"><semantics> <mrow> <mrow> <mo>(</mo> <mo>−</mo> <mo>)</mo> </mrow> </mrow> </semantics></math> in the pie sections implies an inverse relationship between the output variable and the input parameter, while a <math display="inline"><semantics> <mrow> <mrow> <mo>(</mo> <mo>+</mo> <mo>)</mo> </mrow> </mrow> </semantics></math> parameter has a direct relationship with the output.</p> ">
Figure 11
<p>Thermal efficiency sensitivity history in fin tip to (<b>a</b>) thermal conductivity; (<b>b</b>) room temperature; (<b>c</b>) fin length. Solid black line, analytical sensitivity using analytical differentiation in Equation (27), yellow star, sensitivity using Transient Thermal ZFEM, blue square, sensitivity using <math display="inline"><semantics> <mrow> <mn>1</mn> <mo>%</mo> </mrow> </semantics></math> FD and Abaqus built-in element solutions, thin red dashed line, sensitivity using analytical differentiation for steady-state using Equation (27). (<b>d</b>) Steady-state influence of fin parameters in <math display="inline"><semantics> <mi>η</mi> </semantics></math> based on relative sensitivities. In the pie chart, the missing parameters (<math display="inline"><semantics> <mrow> <mi>ρ</mi> <mo>,</mo> <msub> <mi>C</mi> <mi>P</mi> </msub> <mo>,</mo> <msub> <mi>T</mi> <mo>∞</mo> </msub> <mo>,</mo> <msub> <mi>T</mi> <mn>0</mn> </msub> </mrow> </semantics></math>) have an influence of less than <math display="inline"><semantics> <mrow> <mn>1</mn> <mo>%</mo> </mrow> </semantics></math> compared to the sum of all relative sensitivities. Moreover, the <math display="inline"><semantics> <mrow> <mrow> <mo>(</mo> <mo>−</mo> <mo>)</mo> </mrow> </mrow> </semantics></math> in the pie sections implies an inverse relationship between the output variable and the input parameter while a <math display="inline"><semantics> <mrow> <mrow> <mo>(</mo> <mo>+</mo> <mo>)</mo> </mrow> </mrow> </semantics></math> parameter has the opposite relationship with the output.</p> ">
Versions Notes

Abstract

:
Solving transient heat transfer equations is required to understand the evolution of temperature and heat flux. This physics is highly dependent on the materials and environmental conditions. If these factors change with time and temperature, the process becomes nonlinear and numerical methods are required to predict the thermal response. Numerical tools are even more relevant when the number of parameters influencing the model is large, and it is necessary to isolate the most influential variables. In this regard, sensitivity analysis can be conducted to increase the process understanding and identify those variables. Here, we combine the complex-variable differentiation theory with the finite element formulation for transient heat transfer, allowing one to compute efficient and accurate first-order sensitivities. Although this approach takes advantage of complex algebra to calculate sensitivities, the method is implemented with real-variable solvers, facilitating the application within commercial software. We present this new methodology in a numerical example using the commercial software Abaqus. The calculation of sensitivities for the temperature and heat flux with respect to temperature-dependent material properties, boundary conditions, geometric parameters, and time are demonstrated. To highlight, the new sensitivity method showed step-size independence, mesh perturbation independence, and reduced computational time contrasting traditional sensitivity analysis methods such as finite differentiation.

1. Introduction

The physics of transient heat transfer allows one to evaluate the dynamics of energy movements for bodies subjected to changes in temperature, enabling the analysis of the heat flux, temperature history, and temperature rate of different processes. Transient heat transfer problems have significant applications in food processing [1,2], alternative energies [3,4], oil and gas [5], metallurgical processes [6,7], metallography [8,9], infrastructure materials [10,11], and electronics and semiconductor cooling [12,13,14], among others. These applications require analysis of multiple parameters, e.g., time-dependent boundary conditions [15], temperature-dependent thermo-mechanical properties [16,17], and evolving geometries [18] to guide the design of components.
Experimental, analytical, and numerical approaches have been proposed to model transient heat transfer processes. In the experimental approach, level values of the inputs are set to generate random and independent experimental runs. This approach allows the empirical models to capture the dynamic characteristics of a thermal problem regardless of its complexity [19]. It has been used to optimize operational parameters of radiation heat exchange materials [20], induction heating [21], heat treatments [22], sensor design and calibration [23,24], cooling devices [25,26], and electronics refrigeration [27]. However, when the processes of interest involve many parameters, and the factors are prone to uncontrollable changes (e.g., uncontrolled room conditions), the experimental approach becomes time-consuming and expensive, limiting its application.
Approximations for specific processes can be obtained using simplified analytical solutions for the heat diffusion equation for plates [28], walls [29], cylinders [30,31], or spheres [32] under temperature, heat flux, or periodic heating [29]. For instance, predictions in the thermal profile of a plate due to a moving heat source can be approximated by this method [33]. However, explicit analytical solutions are uncommon when systems are large, complex, and time-dependent. The approximations based on simplified models lead to significant errors. As a result, numerical models such as the finite element method (FEM), finite-difference methods, and boundary element methods are popular for solving transient heat transfer problems due to their flexibility and versatility. These numerical methods have been used to approximate the solution of the heat diffusion equation in large and complex problems such as the analysis of fractal piping design [34], phase change materials [35], enhanced surfaces [36], heat exchangers [37], heat sinks [38] and additive manufacturing [39,40], among others. In addition, these and other numerical techniques have been used to study the rate of change of temperature to time in crystallization [2,41], microstructural evolution [42], and superconducting devices [43].
Regardless of the modeling technique, the more elaborate the model, the more challenging it is to track the effects of the input variables on the outputs of the models. Sensitivity analysis is a technique to determine the interactions among the input parameters in a model. For instance, an output could be insensitive to a specific input but highly sensitive to a combination of that input with other parameters in the model. Sensitivity analysis is also used to quantify the influence of these parameters on the model response [44]. In sensitivity analysis, each of the inputs is represented by sensitivity coefficients given by the partial derivative of the output of interest with respect to each input [19,23,24,45,46,47,48,49,50]. The traditional methods to obtain sensitivities are hand-coded analytical differentiation [44], symbolic differentiation [51], automatic differentiation [52], and finite differentiation (FD) [44]. Hand-coded analytical and symbolic differentiation are similar methods as both are exact and rely on analytical models. However, in the hand-code differentiation, the sensitivities are obtained by hand procedures and typed directly in the code. These make the method error-prone, time-consuming to develop and code, and challenging for implicit models. In contrast, symbolic differentiation generates explicit symbolic expressions of explicit analytical models and uses computational symbolic packages to obtain the sensitivities automatically [53]. Automatic differentiation or algorithmic differentiation takes advantage of the chain rule and function composition to sequentially decompose a function in the root variables and its operations in the function. Then the recurrent evaluation of the decomposition in the chain rule of the function allows the algorithm to recover the sensitivity using forward, backward, or adjoint schemes. The formulation of the method makes the algorithm slow and memory demanding [52,53,54]. Another traditional method for obtaining the sensitivities of a model is the numerical FD. This technique has been widely applied in sensitivity analysis in many fields because it is easy to implement as it only requires multiple model evaluations with a small step size for the analyzed variable. Nevertheless, when the model becomes nonlinear, obtaining high precision sensitivities using FD is challenging because the method is subject to computational subtraction cancellation errors and truncation errors [55,56]. In addition, the step size for each input to the model is problem dependent and usually requires a convergence analysis to select it for each of the inputs. Moreover, multiple evaluations of the model are required to calculate each parameter sensitivity of interest, making the technique disk memory intensive and time demanding [54].
An alternative approach for the numerical calculation of sensitivities is the Complex Taylor Series Expansion Method (CTSE) [57]. Unlike finite differencing, CTSE does not require the difference of two analyses. Instead, the derivatives are computed using a small perturbation along the imaginary axis of the input parameters, evaluating the function, extracting the imaginary part of the outputs, and dividing by the perturbation magnitude (more details shown in Section 2.2). In this context, since no subtraction operations are needed, no subtractive error is generated, and the computed sensitivity is highly accurate [56]. CTSE can be incorporated in partial differential equation solvers based on techniques such as the finite-differences method [58], the boundary element method [59], and the FEM [60,61].
The combination between CTSE and analytical modeling has been successfully applied to transient heat transfer analysis. For instance, Jayapragasam et al. [53] showed a sensitivity study for a 1D transient heat conduction problem. The authors compared analytical differentiation, FD, and CTSE methodologies for calculating sensitivities and concluded that CTSE provides certain advantages over other methods, such as near analytical accuracy, step size independency, reduction of computing time, and elimination of the subtractive cancelation problem. Similarly, CTSE has been applied to calculate sensitivities to solve inverse heat transfer problems using the finite-differences method. In this application, the sensitivities calculated using CTSE were used to guide the optimization algorithms to recover the values of temperature-dependent thermal properties and thermal boundary conditions [62,63,64]. Moreover, some authors have reported improved computational times and convergence rates in the minimization algorithms when CTSE is used over other numerical differentiation algorithms [65,66,67].
This work presents and verifies a methodology that integrates CTSE and the transient heat transfer FEM. We call this method the transient heat transfer complex-variable finite element method or Transient Thermal ZFEM. This methodology allows the calculation of faster and more accurate sensitivities of heat flux and temperature history with respect to constant and temperature-dependent thermo-physical properties, boundary conditions, initial conditions, time, and geometry than FD. Moreover, although the essence of the method takes advantage of complex variable operations to calculate sensitivities, it is demonstrated how this method can also be applied in real-variable solvers. A sensitivity analysis for a fin problem under 1D transient thermal conditions was performed to demonstrate our methodology. This device, as represented in Figure 1a, was selected due to the availability of an analytical solution with constant parameters and the relevance of the problem in thermal systems. Fins are present in almost all heat transfer applications such as engines, computers and electronics, and cooling systems. It is worth highlighting that ZFEM has been verified in a number of application areas such as elastic analysis [60], fatigue [68,69], fracture mechanics [61,70,71,72,73], plasticity [74], residual stresses [75], creep [76], variance estimates [77], reliability-based design optimization [78], steady-state thermo-elasticity [79], structural dynamics [80], and periodic material design [81]. However, this is the first demonstration of ZFEM in transient heat transfer problems.
The outline for the rest of the paper is as follows: Section 2 presents the background related to the transient heat transfer FEM and CTSE. Section 3 presents our new methodology that integrates FEM and CTSE to form the Transient Thermal ZFEM methodology. Section 4 demonstrates the application of the methodology on the analysis of the fin problem considering thermal transient conditions for constant parameters (e.g., constant material properties and boundary conditions, see Figure 1). The results are compared against analytical equations. Subsequently, we analyzed the case of temperature-dependent specific heat, and the results were compared against FD. Finally, the conclusions and future work are presented in Section 5.

2. Background

2.1. 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 [ W / mK ] , is presented in Equation (1), where ρ [ kg / m 3 ] corresponds to the density, C P [ J / kg   K ] the specific heat, Q V ˙ [ W / m 3 ] is the volumetric internal heat generation rate, I corresponds to the identity matrix, and is the nabla differential operator.
ρ C P T ( x , y , z , t ) t = k I T T ( x , y , z , t ) + Q V ˙
T ( x , y , z , 0 ) = T 0
T ( x S , y S , z S , t ) = T S = T ¯
T ( x S , y S , z S , t ) = T S = q ¯
T ( x S , y S , z S , t ) = T S = h U ( T T S ) h U = h c o n v + h r a d h r a d = ϵ σ ( T 2 + T S 2 ) ( T + T S )
This equation allows one to find a solution for the temperature T ( x , y , z , t ) [ K ] at time t [ s ] and any position ( x , y , z ) [ m ] in a geometric domain. To solve Equation (1), a time-dependent temperature condition, e.g., the initial temperature ( T 0 [ 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 ¯ [ K ] , in Equation (3)) and predefined heat fluxes ( q ¯ [ W / m 2 ] , in Equation (4)). For example, Figure 1a represents a fin with a predefined wall temperature ( T ¯ = T W ), and insulation at the tip ( 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, ϵ is the material’s surface emissivity, σ = 5.67 × 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 [ W / m 2 K ] is used to represent the contribution from both heat transfer mechanisms, as shown in Equation (5). Here, T [ K ] is the room temperature, h c o n v [ W / m 2 K ] is the convection heat transfer coefficient, and h r a d is the radiation heat transfer. When T 0 and T are similar, h r a d [ W / m 2 K ] 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 [29,32].
Using the finite element method, the representation of the nodal temperatures, T e , in the element is specified by the shape functions, Φ e , given by:
T ( x , y , z , t ) = Φ e T e = j = 1 n T j e Φ j e ( x , y , z , t )
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 ). In this equation, M e [ J / K ] is the element capacitance matrix, K c e   [ W / K ] is the symmetric element conductivity matrix, K h c o n v e [ W / K ] and K h r a d e [ W / K ] are the asymmetric convective and radiative flux matrixes, respectively. The right-hand side terms include the heat flux vectors f q ¯ e [ W ] , f c o n v e   [ W ] , f r a d e   [ W ] , and the volumetric heat generation vector f Q V e   [ W ] for the element e .
M e T ˙ e + ( K c e + K c o n v e + K r a d e ) T e = f q ¯ e + f c o n v e + f r a d e + f Q V e
To obtain the temperature history from Equation (7), the generalized Crank-Nicholson method for the time integration recurrence scheme is applied. This integration scheme is obtained by a trapezoidal approximation of the elemental temperature rate T ˙ e and solving for T t s + 1 e . Using the substitutions K e = K c e + K c o n v e + K r a d e and f e = f q ¯ e + f c o n v e + f r a d e + f Q V e , the time integration scheme is given by [82]:
( 1 Δ t M e + β   K e ) T t s + 1 e = ( 1 Δ t M e ( 1 β ) K e ) T t s e + ( 1 β ) f t s e + β   f t s + 1 e
where Δ t is the time increment and the temperature of the time t s + 1 is obtained based on the values of M e , K e , T e , and f t s e , and the next value of f t s + 1 e . For linear problems, f q ¯ e , f c o n v e , f r a d e , f Q V e are constant, then the flux vectors f t s e = f t s + 1 e , therefore ( 1 β ) f t s e + β   f t s + 1 e = f t s e . Setting β = 1 in Equation (8) minimizes numerical integration problems, prevents spurious oscillations, and ensures stability for the solution [82,83]. With this condition, the system of equations becomes the fully implicit system A e T t s + 1 e = R H S e , where the coefficient matrix for the time increment ( A e [ W / K ] ) and the right-hand side of the equation ( R H S e [ W ] ) are defined in Equations (9) and (10), respectively.
A e = ( 1 Δ t M e + K e )
R H S e = 1 Δ t M e T t s e + f t s e
It is noted that when 2D conditions are assumed, the FEM formulation for 2D elements considers that the temperature gradient in the out-of-plane direction is null. For instance, if the problem is modeled in the x y plane, then T / z = 0 . This implies that the heat flux in the z -direction is zero, and the element is insulated along with the front and back faces. To overcome this issue, Morville et al. [84] and La Batut et al. [85] proposed the use of the volumetric sink term ( Q V in Equation (1)) to enhance the heat transfer in a 2D model. Q V is defined in Equation (11) with U V [ ( W / ( m 2 K ) ) / m ] corresponding to the heat transfer coefficient per unit of thickness. In this equation, the temperature-dependent heat sink dissipates more heat as the temperature difference increases. This is similar to having surface convection/radiation heat transfer in the out-of-plane direction. In addition, U V should be tuned to minimize the error of the temperature response of the FEM and a reference solution.
Q V = U V ( T T )

2.2. Numerical Differentiation by CTSE

CTSE is used to calculate the numerical sensitivities of any function g ( ϕ ) with respect to any input parameter ϕ . To obtain the sensitivity g ( 1 ) ( ϕ ) = g / ϕ , where the superscript indicates the order of differentiation, ϕ is perturbed by a small amount h (i.e., h 10 10 [61]) in the imaginary axis. Then, the evaluation of the function with the complex input, g ( ϕ + i h ) , is obtained by considering a Taylor’s series expansion as shown in Equation (12). Truncating the infinite series to the first-order term of h and solving for the first-order sensitivity g ( 1 ) ( ϕ ) , Equation (13) presents the formulation for the numerical sensitivity of the function g for the parameter of interest ϕ [56]. If h is small enough, the truncation of the higher-order terms is accurate because the powers of h of order two and higher will be approximately zero. Moreover, the original response of the solution variable is returned unchanged in the real component.
g ( ϕ + i h ) g ( ϕ ) + g ( 1 ) ( ϕ ) i h 1 ! + g ( 2 ) ( ϕ ) ( i h ) 2 2 ! +
g ( 1 ) ( ϕ ) = g ( ϕ ) ϕ I m [ g ( ϕ + i h ) ] h
The sensitivities obtained from Equation (13) have the units of the output parameter ( g ) divided by the perturbed parameter ( ϕ ). To simplify the analysis of sensitivities, it is recommended to transform them into relative coefficients [45]. The calculation of these relative coefficients can be achieved by the application of Equation (14), where the relative sensitivities of g with respect to the parameter of interest ϕ or S ϕ g , have the same units of g .
S ϕ g ( x , y , z , t ) = | ϕ | g ( x , y , z , t ) ϕ
On the other hand, when analyzing sensitivities regarding T , T 0 , and T ¯ in heat transfer problems, it is common to normalize these sensitivities using the maximum temperature difference in the system ( Δ T m a x [ K ] ) instead of the values of T , T 0 , or T ¯ . For example, in the extended surface of Figure 1a, | ϕ | = | Δ T m a x | = | T W T | [32]. This is commonly used to normalize temperature differences in thermal systems because the sensitivity to those temperatures can be explained as the change of the output g due to a small change Δ in T , T 0 , or T ¯ .

3. Sensitivity Analysis in Transient Heat Transfer Problems Using Transient Thermal ZFEM

To obtain sensitivity information for transient heat transfer problems, we propose a methodology that combines the solution of the transient heat transfer equation using the finite element method with the Complex Taylor Series Expansion. We call this methodology Transient Thermal ZFEM, and the main concepts are presented in Figure 1b. To begin, all of the components in the FEM formulation, such as nodal coordinates, material properties, and boundary conditions, become complex-valued variables [79]. Then, to obtain sensitivities of T or any other output of the Transient Thermal ZFEM model with respect to the input parameters ϕ , different perturbation schemes for ϕ should be followed depending on the kind of input, as follows:
  • Constant inputs: to obtain sensitivities with respect to the input variables in the model that are constant, the procedure presented in Equation (13) can be applied with only a small imaginary perturbation h applied to the variable of interest. For instance, to obtain sensitivities with respect to T , we evaluated the model is evaluated with a value of T * = T ( 1 + h i ) .
  • Geometry sensitivities: to obtain sensitivities with respect to geometric shapes, perturbations to the geometry are required. In Transient Thermal ZFEM, the domain shape is represented by the nodal coordinates. As previously mentioned, the nodal coordinates are converted to complex-variable. For reference, Figure 2a shows a mesh where the nodal coordinates are in the plane xy, and the out of plane represents the imaginary component I m . Moreover, Figure 2b shows the representation of a planar element with real and imaginary coordinates. To obtain sensitivities with respect to the geometry, the imaginary component of the nodal coordinates should be perturbed according to the following guidelines: First, the direction of the perturbation should be selected. For instance, considering sensitivities with respect to the length b [ m ] in the fin problem (see Figure 2a), as b is parallel to the x - axis and its boundaries normal to that axis, the perturbation should be applied by translating the imaginary nodal coordinates along the direction of positive I m axis. Second, a fraction of the nodes ( λ ) should be selected in the mesh to be perturbed. As presented in Figure 2a, λ can take values from zero to one. The selected value for this fraction depends on the balance between precision and computational time. This parameter is problem-dependent and could be determined with convergence analysis. However, previous works have shown that considering a value of λ large enough to perturb at least two elements normal to the boundary is enough for most problems [60,71,86,87]. Next, a perturbation function P ( x ) is selected. In Figure 2a a linear P ( x ) is considered. This function has a minimum value of zero and a maximum value of one. The maximum and minimum values of P ( x ) are located along the geometric boundaries of dimension λ b . Finally, the fraction of the imaginary coordinates of the nodes are shifted in the selected imaginary direction, a magnitude h × P ( x ) , and the sensitivity is obtained following Equation (13). An example of this perturbation scheme is presented in Figure 2a, where the wireframe mesh presents the real components of the nodal coordinates and surface pointing out of the plane, representing the perturbation on the imaginary axis.
  • Temperature-dependent properties: to obtain sensitivities with respect to temperature-dependent properties, the properties should be defined for any given temperature; therefore, fitting or interpolating from experimental values is necessary. After fitting, the perturbation is applied by shifting h units in the imaginary axis of the resulting temperature-dependent curve, and the perturbed function is used to perform the calculations as usual. This procedure is analogous to the perturbation of a constant parameter. Finally, Equation (13). is used to obtain the numerical sensitivities. An example of the fin problem is presented in Section 4.3.
  • Time sensitivities: to obtain sensitivities with respect to time, the time increment Δ t is perturbed by a value h in the imaginary component. Then the standard solution method is applied. After the model is solved, the time sensitivities are calculated using Equation (15). The time function in the denominator of Equation (15) compared to Equation (13) is a consequence of the numerical time integration of the imaginary perturbation associated Δ t where # i n c is the total number of converged times increments in the solution. A convergence study should be performed on the time increment to determine the correct balance between computational time and accuracy.
    g t I m [ g ( Δ t + h i ) ] h × t t m a x × # i n c

Solution of Complex-Valued System of Equations: Cauchy-Riemann Matrix Notation

As most numerical packages do not support complex variable linear algebra operations, it is possible to take advantage of the multicomplex algebra isomorphism to represent any complex number as a matrix with all real numbers [79]. This representation of complex numbers is known as the Cauchy-Riemann (CR) matrix notation and allows for the use of matrix functions and arithmetic operations as homologous of the operations between complex numbers [88].
Consider the complex-variable A e matrix and R H S e vector in Equations (9) and (10), respectively, these quantities can be represented as the CR matrices following Equation (16) and CR vectors following Equation (17). This representation implies the addition of extra degrees of freedom in the Transient Thermal ZFEM model. Consequently, a duplicated mesh approach is used to consider the extra degrees of freedom provided by the sensitivity information in the imaginary components. In Figure 2b, a 2D Transient Thermal ZFEM element with the duplicated mesh approach is presented. In this case, the original mesh of the FEM (nodes 1 to 4) will carry the real information, while the secondary or imaginary mesh (nodes 11 to 14) contains the sensitivity of the output with respect to the perturbed input. A four-node—four-quadrature point complex 2D isoparametric element is shown in Figure 2c.
A e = A R e e + i A I m e = [ A R e e A I m e A I m e A R e e ]
R H S e = R H S R e e + i R H S I m e = [ R H S R e e R H S I m e ]

4. Numerical Example: Transient Fin Problem

We implemented the method in 2D User Element Subroutine (UEL) within the commercial finite element analysis software Abaqus to demonstrate the performance of the Transient Thermal ZFEM method. Abaqus cannot solve systems with complex variables directly because its built-in equation solver is made for real-valued components. As a result, the complex-variable right-hand side of the element ( R H S e ), and the complex-variable coefficient matrix for the element ( A e ) are transformed to the Cauchy-Riemann form of all real numbers before returning the results to Abaqus. Then, the global system of equations is assembled internally by Abaqus and solved for each time increment. This process is presented in the flowchart of Figure 3, and the full version of the UEL is presented in the Supplementary Material.
The method was verified by analyzing the heat transfer behavior for the thin fin shown in Figure 1a. The geometry of the fin is characterized by the length b , height δ , and thickness L , and it is made from a material with conductivity k , density ρ , and specific heat C P . The results computed using ZFEM were compared with the analytical solution using symbolic differentiation and Abaqus simulation with 2D DC2D4 elements and finite differentiation. The errors were measured using the Normalized Root Mean Square Error (NRMSE) metric as described in Equation (18). In this equation, o corresponds to the order of the sensitivity of g , g ( o ) ¯ is the reference results, and # i n c is the total number of sampling points which is coincident with the total number of increments in the transient thermal problem.
NRMSE = 1 max ( g ( o ) ¯ ) min ( g ( o ) ¯ ) ( g ( o ) ¯ g ( o ) ) 2 # i n c
The analytical solution for the temperature as a function of time for the fin problem was obtained by considering the dimensionless parameters for the position X in Equation (19), the time τ in Equation (20), N 2 in Equation (21), λ n in Equation (22), the temperature θ in Equation (23), and by adding the steady-state θ s s ( Χ ) and the transient solutions θ τ ( Χ , τ ) as shown in Equations (24) and (25), respectively [28]. For this solution, it was also assumed that the axis x and X were aligned with the length of the fin, and the fin had an adiabatic tip.
Χ = x b
τ = t k b 2 ρ C P
N 2 = 2 h U b 2 k δ
λ n = 2 n 1 2 π
θ = T T T W T = θ s s ( Χ ) + θ τ ( Χ , τ )
θ s s ( Χ ) = cosh ( N Χ ) cosh ( N )
θ τ ( Χ , τ ) = 2 n = 1 ( 1 ) n + 1 ( θ 0 λ n λ n 2 N 2 + λ n 2 ) c o s ( λ n Χ ) e ( N 2 + λ n 2 ) τ
In this problem, the fin was in thermal equilibrium with the environment at t = 0   s , meaning that the room temperature ( T ) was equal to the initial temperature ( T 0 ). Then, a prescribed temperature T w was applied at the base of the fin, i.e., x = b . As soon as the fin started to heat, energy dissipation occurred due to the film conditions in the top, bottom, back, and front surfaces determined by T and the global heat transfer coefficient h U . In the case of ZFEM and Abaqus built-in models, combined convection and radiation boundary conditions were applied on the top and bottom edges, and T w was applied to the left edge as depicted in Figure 1a. In addition, enhanced dissipation was considered to obtain the effects of out-of-plane dissipation. The convergence analysis for U V is presented in Figure 4a, where the optimal was found at 735 ( W / ( m 2 K ) ) / m and is highlighted with the blue marker. The temperature and sensitivity information were extracted from the node located in the centerline of the fin’s tip. For reference, the analytical solution at that point is presented as a solid black line in Figure 4b. The mesh was discretized along with the height, δ [ m ] , using 10 equally spaced elements, and 1000 elements along the length, b , and the time increment for the solver was set to Δ t = 1   s . We used convergence analyses to set the values of discretization in space and time. We varied the number of elements in the domain by one order of magnitude up to 10 5 elements. A similar approach was followed for the time discretization up to a maximum of 10 4 time increments. Full details about the convergence analyses are presented in the Supplementary Material.
Using the parameters in Table 1 (obtained from [28]), the analytical dimensionless temperature history, ( θ ) , at the fin’s tip ( X = 0 ) is presented in Figure 4b. In this figure, the analytical θ history is presented as a solid black line, θ S S as a dashed red line, and the results using the Transient Thermal ZFEM UEL and the built-in Abaqus functions are overlaid using a yellow star and a blue square, respectively. This convention will be preserved for the rest of the results presented in this work. The sampling rate for the figures was every 14 time steps for visualization purposes. From the results presented in Figure 4b, it is noted that as the temperature approaches the steady-state (approximately at t = 350   s ), and the simulations using the Transient Thermal ZFEM UEL and the built-in Abaqus functions match the temperature history of the analytical solution.
Although the Transient Thermal ZFEM method for obtaining sensitivities can be applied to any input variable in the models, out of briefness, the next sections will focus on obtaining sensitivities for the transient behavior of a material property ( k ), a boundary condition ( T ), a geometric parameter ( b ), a temperature dependent-property ( C P ( T ) ), and for the calculation of heating rates ( θ / t ), which will cover all the perturbation cases described in Section 3. The results for the remaining parameters affecting the fin’s thermal transient behavior are relegated to the Supplementary Material.

4.1. Sensitivity Analysis of Temperature History

As stated in Section 3, the first step to obtain sensitivities using Transient Thermal ZFEM is to apply a perturbation, h, to the variable of interest. While many other references have shown the step size independence for CTSE when applied in other fields, a step size analysis was conducted for transient thermal analysis to verify previous findings. In addition, the optimal step size for FD was also computed for comparison purposes.
Considering the analytical equations as a reference for the calculation of the error, the NRMSE for the responses as a function of the step size for Transient Thermal ZFEM and FD calculated with built-in Abaqus functions is presented in Figure 5a, using continuous and dashed lines, respectively. All the FD results were obtained by using a forward differentiation scheme and first-order accuracy, which is analogous to applying a positive imaginary perturbation in the CTSE method integrated into Transient Thermal ZFEM. In agreement with previous literature about numerical sensitivities [89], for values above h × ϕ with h = 10 4 , the truncation error completely cancels the precision of the sensitivity response in the FD method. However, the Transient Thermal ZFEM is completely insensitive to the step size for the analyzed input parameters. The change in the trend for the dashed red and dashed green lines in FD is related to a shift in the dominance of the truncation error to round-off error. This analysis highlights one of the main advantages of the Transient Thermal ZFEM over the regular FD in that CTSE is dramatically less prone to step size inaccuracies. Based on the perturbation step analysis study, the rest of the results presented herein were calculated using h = 10 10 ϕ for the Transient Thermal ZFEM and h = 1 % ϕ for FD.
For geometrical perturbations, it is required to perturb a fraction of nodes in the mesh. Besides using the rule of thumb for this parameter, a convergence for study for λ was performed. Figure 5b, shows the NRMSE for the sensitivity with respect to b, S b θ , as a function of λ . In this figure, the top insert shows the rule of thumb suggested λ , and the bottom inset shows a completely perturbed mesh. For the transient thermal fin problem, the precision of S b θ is insensitive to λ , as shown by the almost flat line of Figure 5b. Therefore, the rule of thumb and the case when λ = 1 have almost the same precision. For illustration purposes in the fin problem, we picked λ = 1 for the sensitivity analysis of the fin length. This means that all the imaginary coordinates of the mesh have non-zero values as shown in the bottom right inset in Figure 5b.
The sensitivities history of θ concerning the thermal conductivity, S k θ , the room temperature, S T θ , and the fin’s height, S b θ , are shown in Figure 6a–c, respectively. A negative relative sensitivity predicts a decreasing trend for the temperature with an increment of the input parameter and vice versa. Overall, a good agreement between analytical differentiation, FD, and Transient Thermal ZFEM is observed. The solutions from the simulations using Transient Thermal ZFEM had a better match than built-in Abaqus functions with FD match when compared with the analytical values, as could be seen by error measurements in Table 2. These errors are on the same order as the NRMSE for θ ; therefore, we attribute the good agreement among methodologies to the enhanced dissipation implemented to counter the lack of heat transfer in the out-of-plane direction in the 2D FEM and Transient Thermal ZFEM elements.
From the time-dependent sensitivity analysis in Figure 6, it can be stated that T is only influential in the transient state, since S T θ 0 at steady state while k and b are not influential at the beginning, but they are influential in steady-state given that S k θ 0 and S b θ 0 at 425 . Based on the information provided by Figure 6a, a less conductive fin will reduce θ given that S k θ ( t ) 0 . On the other hand, Figure 6b,c show that longer fins exposed to higher room temperatures will reduce the fin’s tip temperature. This is explained by the inverse relation of θ , S b θ ( t ) , and S T θ ( t ) . In the current application, reducing the fin’s tip is desired because it means that more energy can be dissipated in the system before reaching the tip.

Computational Efficiency of Transient Thermal ZFEM

For the numerical solutions to the transient fin problem, the computational system used consisted of an Intel Xeon W-2265 at 3.50 GHz, 128 Gb of RAM, Abaqus 2021 Hot Fix 5, and Ubuntu 20.04 Long Term Support. In this specific case, the total CPU time of three pairs of sequential and independent simulation runs made with ZFEM and FD were conducted and averaged. The results showed that on average, ZFEM required 45 % less total CPU time to solve the Transient Thermal ZFEM model and obtain sensitivities compared to using built-in Abaqus functions and FD. This contrasts with that reported in the ZFEM literature for more complicated physics, where 2.5 to 4 times longer computational times were obtained for equivalent real variable runs [71,74,87]. This behavior shows another advantage of using CTSE for the transient thermal linear problem as it requires only one evaluation of the model to obtain sensitivities. In contrast, at least two evaluations of the models are required for FD, assuming the optimal step size is known a priori.

4.2. Heating Rate

Figure 7 shows the results for the heating rate, θ / t , at the fin’s tip. This rate was obtained by directly evaluating the Transient Thermal ZFEM model and applying the methodology shown in Section 3. In contrast, one is required to post-process the built-in Abaqus solutions by computing FD using the time domain to obtain the heating rate. The NRMSE for Transient Thermal ZFEM and Abaqus built-in functions with FD is 0.012 and 0.010, respectively. The small deviations with respect to the analytical solutions can be attributed to the time increment used in the simulations. A detail about the influence of the time incrementation in the heating rate can be seen in the Supplementary Material.

4.3. Sensitivity Analysis of Temperature History for a Fin with Temperature-Dependent Specific Heat

When any coefficient of Equation (1) is temperature-dependent, the transient heat transfer problem becomes nonlinear and analytical equations are no longer available [82]. To test the capacity of Transient Thermal ZFEM to handle nonlinearities, a model containing temperature-dependent specific heat was created and solved. For this analysis, the specific heat was selected as the variable of interest because it is temperature-dependent in the model and it directly modifies the thermal capacitance of the material. Moreover, according to Equations (1) and (7), this property is expected only to affect the process in the transient state. The values for C P ( T ) used in the model are presented in Figure 8a and correspond to Ti64 alloy [90]. Due to the lack of analytical equations, the response obtained from built-in Abaqus functions combined with FD was used as a reference to calculate the NRMSE. To calculate the reference response for S C p θ with FD, one model was run with the nominal parameters from the curve shown in Figure 8a, and for the second model, all the values in the curve defining C P ( T ) were shifted 5.5   J / Kg   K which corresponds to 1% of the minimum value of the specific heat. The convergence analysis for these models is presented in the Supplementary Material. For the calculation of the sensitivities using Transient Thermal ZFEM, the procedure described in Section 3 was followed as: (a) first interpolating the experimental values for C P ( T ) in the temperature range of the solution, and then (b) applying an imaginary perturbation, h = 10 10 to all values as shown in the inset of Figure 8a. Finally, to normalize the sensitivity results, we used | ϕ | = 580   J / Kg   K for Equation (10) which corresponds to the specific heat reported in Table 1.
The θ response for Transient Thermal ZFEM and built-in Abaqus functions is reported in Figure 8b, and the sensitivity response is presented in Figure 8c. A good agreement was found between both models with an NRMSE of 0.003 for θ , and 0.007 for S C p θ . These results show that the Transient Thermal ZFEM can accurately determine sensitivities considering material nonlinearities and temperature-dependent parameters.

4.4. Sensitivity Analysis of the Heat Flux History: Defining the Fin’s Performance

Using the heat flux output, two measurements of performance were calculated for the fin: heat transfer effectiveness, ε , calculated using Equation (26), and the thermal efficiency, η , defined in Equation (27). The heat transfer effectiveness corresponds to the ratio of the heat flux ( q ) from a wall with a fin to the heat flux of the same wall without a fin, and the efficiency is the relationship of effective heat transfer from the fin to the room regarding the heat transfer of an ideal fin kept at the initial temperature θ 0 [24]. The sensitivity analysis of these two performance indices represents useful information for improving the design of fins since it can identify which variables positively affect their performance. In this regard, the main idea is to reduce the values of parameters with negative sensitivities and increase the importance of parameters with positive sensitivities.
ε ( t ) = q ( θ ( Χ = 1 , t ) ) q ( Χ = 1 )
η ( t ) = Q ˙ ( θ ( Χ , t ) ) Q ˙ ( θ ( Χ ) = θ 0 ) = ε ( t ) δ 2 b
To calculate the performance parameters using Transient Thermal ZFEM and built-in Abaqus functions, the heat flux was obtained at the wall-fin interface ( x b ). The evolution of ε and η are presented in Figure 9a,b, respectively, showing good agreement between the Transient Thermal ZFEM and the Abaqus built-in element response and with NRMSE values of 0.052 and 0.052, respectively, when compared with the analytical equations. η and ε quickly reach the steady-state value with an asymptotic behavior tending to positive infinity in t = 0   s . This behavior is explained by the discontinuity of the temperature profile at the wall-fin interface at that time. In addition, this is similar to the behavior of the transient heat transfer coefficient studied in similar problems [30,31].
A comprehensive sensitivity analysis of the fin performance parameters involving all of the parameters affecting the model was performed. The time history evolutions for all sensitivities are shown in the Supplementary Material, and only the results for k , T , and b , are included in Figure 10 for ε and Figure 11 for η . The maximum error for the sensitivities corresponded to 0.048 for S ε and 0.048 for S η when compared to the analytical solution. As observed in the time history for the sensitivities, the values quickly stabilized to the steady-state. Therefore, pie charts were developed showing the values of the contribution of the individual sensitivities to the total sum of the sensitivities for all variables in the model after steady-state conditions have been reached. Figure 10d and Figure 11d show the sensitivities for ε and η , respectively. In the pie charts, the ( + ) symbol inside each of the colored sections implies a direct relationship between the output and the input variable while the ( ) implies an inverse relationship.
From Figure 10d, it is possible to identify that δ , h U , and k are the most influential parameters to modify ε . Each of these variables accounts for nearly a third of the total sensitivity of the effectiveness. Based on Equation (14), if δ , h U , or k are duplicated, an equivalent thermal effectiveness increment (or reduction) in the fin at steady-state will occur. In this regard, the ( ) inside the pie sections for δ and h U means that they should be reduced to increase the effectiveness. The other parameters, e.g., b, ρ , C P , T , and T 0 influence the effectiveness less than the 1% compared to the sum of all relative sensitivities. This means that an increment or reduction of these five variables will not produce a significant change in the effectiveness of the fin.
Similarly, from Figure 11d, it is possible to identify that δ , b , h U , and k are the most influential parameters to modify the thermal efficiency. Therefore, increments in δ and k will increase this index. For the efficiency, the weights are not equally split as in the case of the effectiveness. For instance, the fin length is almost twice more influential than h U , and k . Moreover, the fin thickness is approximately 2.4 times less influential in η than the length of the fin. This means that larger geometries will reduce the thermal efficiency of the fin while enhancing the heat transfer with more conductive materials larger heat transfer coefficients are beneficial for η .
From these results, some guidelines for the fin’s design can be stated: To cool a heated wall, it is preferred to have a fin with a low aspect ratio made of a material with high conductivity and low heat transfer coefficient. This boundary condition requirement is counterintuitive because larger heat transfer coefficients typically improve the heat exchange from the external surfaces to the room. However, this also means that the heat will be transferred only at the beginning of the fin. Therefore, according to the first-order sensitivity information, a shorter and thicker fin made of more conductive materials is preferred to increase the dissipation performance.
It is worth highlighting that although we show the Transient Thermal ZFEM methodology in a simple thermal problem, the methodology is general. For instance, we have included an additional verification of the capabilities of the Transient Thermal ZFEM in the Supplementary Material for a 2D problem corresponding to a rotating disk under nonconstant initial temperature.

5. Conclusions

In this work, the complex Taylor series expansion method was integrated with the transient heat transfer FEM to create the Transient Thermal ZFEM methodology. This methodology enables the calculation of high precision sensitivities regarding geometrical parameters, boundary conditions, time, and material properties. Although the method takes advantage of complex algebra to obtain the sensitivities, we do not require complex variable solvers. Instead, the system of equations is solved using the Cauchy-Reimann matrix representation of a complex number. This transformation allows the solution of the FEM equations in any real-variable solver facilitating the implementation of the new method within commercial finite element software.
We verified the Transient Thermal ZFEM method against the analytical solution for a fin with an insulated tip. Our methodology accurately computed temperature and temperature sensitivity history regarding specific heat, thermal conductivity, ambient temperature, and fin geometry. Good agreement with the reference solutions was obtained for the numerical example. Similarly, the thermal performance of the fin based on the fin’s heat flux presented a good agreement with the analytical solution. The sensitivities using Transient Thermal ZFEM were compared with FD using built-in Abaqus elements. Overall, we found a reduction of 45% in the computational time compared with FD. Moreover, we demonstrated that the Transient Thermal ZFEM is step-size independent and mesh perturbation independent for the transient thermal linear problem. On the other hand, the time derivative obtained by Transient Thermal ZFEM is highly dependent on the number of increments of time but only limited by the computing hardware.
Implementing the Transient Thermal ZFEM in finite element analysis programs will enable the sensitivity analysis of heat transfer models with a fraction of the computational cost and higher precision than sensitivity analysis using the traditional finite differentiation method. Moreover, Transient Thermal ZFEM can improve the convergence rates and computational times on gradient-based optimizations, especially in inverse heat transfer problems. These improvements are possible because the Transient Thermal ZFEM methodology allows one to obtain a more precise and computationally efficient calculation of partial derivatives than the traditional methods.
Given the advantages of Transient Thermal ZFEM over the traditional sensitivity analysis methodologies, our methodology is especially attractive for complicated systems in alternative energies, petroleum engineering, chemical, agriculture, materials, and metallurgical industries. Accurate sensitivity information in transient thermal problems allows one to understand the system’s thermal behavior better. For instance, in the numerical example, using Transient Thermal ZFEM, the variables can be sorted by their influence in the transient state and steady states. Then, it was possible to effectively filter the variables that affect the processes after reaching the steady-state, even if the parameters of the systems are temperature dependent.
As future work, this method will be extended to calculate higher-order sensitivities and to explore other multiphysics problems such as thermo-mechanical, microstructural evolution, crystal plasticity, heat treatments, and thermo-electrical. Moreover, we will explore the use of sensitivities to perform a derivative-based approximation of functions that enable predictions in complicated problems or to evaluate models with output uncertainty related to the input parameters for uncertainty quantification. In addition, these sensitivities could guide semiempirical techniques that rely on simulations to recover basic parameters of the process as the Nusselt numbers, Biot numbers, or Prandtl numbers.

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/app12052738/s1, Supplementary Material 1: Transient Thermal ZFEM user element for ABAQUS; Supplementary Material 2: Convergence analysis and detailed sensitivity analysis for θ ; Supplementary Material 3: Convergence study for temperature dependent specific heat model; Supplementary Material 4: Detailed sensitivity analysis of fin performance; Supplementary Material 5: numerical sensitivities of a rotating disk under nonconstant initial temperature.

Author Contributions

Conceptualization, J.-S.R.-T. and D.R.; Methodology, J.-S.R.-T., D.R., D.R.-T. and A.M.; software, J.-S.R.-T., D.R.-T., A.M. and H.M.; validation, J.-S.R.-T. and J.C.V.-G.; writing—original draft and editing, J.-S.R.-T. and D.R.; visualization, J.-S.R.-T. and J.C.V.-G.; writing—review, A.M., H.M. and D.R.; supervision, A.M., H.M. and D.R.; project administration, D.R.; Funding acquisition, H.M. and D.R. All authors have read and agreed to the published version of the manuscript.

Funding

This work was funded by the Department of Defense of the United States of America [grant number W911NF2010315].

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Acknowledgments

This work was supported by the Department of Defense of the United States of America [grant number W911NF2010315].

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Joardder, M.U.H.; Hasan Masud, M.; Joardder, M.U.H.; Masud, M.H. Possible solution of food preservation techniques. In Food Preservation in Developing Countries: Challenges and Solutions; Springer International Publishing: Berlin, Germany, 2019; pp. 199–218. [Google Scholar]
  2. Mo, J.; Groot, R.D.; McCartney, G.; Guo, E.; Bent, J.; van Dalen, G.; Schuetz, P.; Rockett, P.; Lee, P.D. Ice crystal coarsening in ice cream during cooling: A comparison of theory and experiment. Crystals 2019, 9, 321. [Google Scholar] [CrossRef] [Green Version]
  3. Zou, H.; Pei, P.; Wang, C.; Hao, D. A numerical Study on heat transfer performances of horizontal ground heat exchangers in ground-Source heat pumps. PLoS ONE 2021, 16, e0250583. [Google Scholar] [CrossRef] [PubMed]
  4. Kant, K.; Biwole, P.H.; Shukla, A.; Sharma, A.; Gorjian, S. Heat transfer and energy storage performances of phase change materials encapsulated in honeycomb cells. J. Energy Storage 2021, 38, 102507. [Google Scholar] [CrossRef]
  5. Huan, H.; Liu, L.; Huan, B.; Chen, X.; Zhan, J.; Liu, Q. A theoretical investigation of modelling the temperature measurement in oil pipelines with edge devices. Meas. J. Int. Meas. Confed. 2021, 168, 108440. [Google Scholar] [CrossRef]
  6. Kotrbacek, P.; Bellerova, H.; Luks, T.; Raudensky, M. Heat transfer correlations for secondary cooling in continuous casting. Steel Res. Int. 2021, 92, 2000465. [Google Scholar] [CrossRef]
  7. Stieven, G.d.M.; Soares, D.d.R.; Oliveira, E.P.; Lins, E.F. Interfacial heat transfer coefficient in unidirectional permanent mold casting: Modeling and inverse estimation. Int. J. Heat Mass Transf. 2021, 166, 120765. [Google Scholar] [CrossRef]
  8. Inyushkin, A.V.; Taldenkov, A.N.; Ralchenko, V.G.; Bolshakov, A.P.; Khomich, A.V. Isotope effect in thermal conductivity of polycrystalline CVD-Diamond: Experiment and theory. Crystals 2021, 11, 322. [Google Scholar] [CrossRef]
  9. Gu, C.; Lian, J.; Bao, Y.; Xiao, W.; Münstermann, S. Numerical study of the effect of inclusions on the residual stress distribution in high-Strength martensitic steels during cooling. Appl. Sci. 2019, 9, 455. [Google Scholar] [CrossRef] [Green Version]
  10. Thi, V.D.; Khelifa, M.; Oudjene, M.; El Ganaoui, M.; Rogaume, Y. Finite element analysis of heat transfer through timber elements exposed to fire. Eng. Struct. 2017, 143, 11–21. [Google Scholar] [CrossRef]
  11. Lee, H.; Seong, J.; Chung, W. Correlation analysis of heat curing and compressive strength of carbon nanotube–Cement mortar composites at Sub-Zero temperatures. Crystals 2021, 11, 1182. [Google Scholar] [CrossRef]
  12. Righetti, G.; Zilio, C.; Doretti, L.; Longo, G.A.; Mancin, S. On the design of phase change materials based thermal management systems for electronics cooling. Appl. Therm. Eng. 2021, 196, 117276. [Google Scholar] [CrossRef]
  13. Alhusseny, A.; Al-Fatlawi, A.; Al-Aabidy, Q.; Nasser, A.; Al-Zurfi, N. Dissipating the heat generated in high-Performance electronics using graphitic foam heat-Sinks cooled with a dielectric liquid. Int. Commun. Heat Mass Transf. 2021, 127, 105478. [Google Scholar] [CrossRef]
  14. Nakazawa, Y.; Imajo, S.; Matsumura, Y.; Yamashita, S.; Akutsu, H. Thermodynamic picture of dimer-Mott organic superconductors revealed by heat capacity measurements with external and chemical pressure control. Crystals 2018, 8, 143. [Google Scholar] [CrossRef] [Green Version]
  15. Cui, W.; Gawecka, K.A.; Taborda, D.M.G.; Potts, D.M.; Zdravković, L. Time-Step constraints in transient coupled finite element analysis. Int. J. Numer. Methods Eng. 2016, 106, 953–971. [Google Scholar] [CrossRef] [Green Version]
  16. Wcisło, B.; Pamin, J. Local and non-Local thermomechanical modeling of elastic-Plastic materials undergoing large strains. Int. J. Numer. Methods Eng. 2017, 109, 102–124. [Google Scholar] [CrossRef]
  17. Jahan, S.A.; El-Mounayri, H. A thermomechanical analysis of conformal cooling channels in 3D printed plastic injection molds. Appl. Sci. 2018, 8, 2567. [Google Scholar] [CrossRef] [Green Version]
  18. Jantos, D.R.; Hackl, K.; Junker, P. An accurate and fast regularization approach to thermodynamic topology optimization. Int. J. Numer. Methods Eng. 2019, 117, 991–1017. [Google Scholar] [CrossRef]
  19. Moges, T.; Yan, W.; Lin, S.; Ameta, G.; Fox, J.; Witherell, P. Quantifying uncertainty in laser powder bed fusion additive manufacturing models and simulations. In Proceedings of the Solid Freeform Fabrication 2018: Proceedings of the 29th Annual International Solid Freeform Fabrication Symposium—An Additive Manufacturing Conference, Austin, TX, USA, 13–15 August 2018; pp. 1913–1928. [Google Scholar]
  20. Pierre, T.; Jiménez-Saelices, C.; Seantier, B.; Grohens, Y. Transient pulsed technique to characterize the radiative and conductive properties of bio aerogels. Int. J. Therm. Sci. 2017, 116, 63–72. [Google Scholar] [CrossRef]
  21. Frivaldsky, M.; Pavelek, M.; Donic, T. Modeling and experimental verification of induction heating of thin molybdenum sheets. Appl. Sci. 2021, 11, 647. [Google Scholar] [CrossRef]
  22. Zhang, F.; Zhang, J.; Ni, H.; Zhu, Y.; Wang, X.; Wan, X.; Chen, K. Optimization of AlSi10MgMn alloy heat treatment process based on orthogonal test and grey relational analysis. Crystals 2021, 11, 385. [Google Scholar] [CrossRef]
  23. Alam, T.; Kumar, R. Stagnation point transient heat flux measurement analysis from coaxial thermocouples. Meas. J. Int. Meas. Confed. 2018, 128, 352–361. [Google Scholar] [CrossRef]
  24. Manjhi, S.K.; Kumar, R. Stagnation point transient heat flux measurement analysis from coaxial thermocouples. Exp. Heat Transf. 2018, 31, 405–424. [Google Scholar] [CrossRef]
  25. Xu, J.; Zhang, K.; Duan, J.; Lei, J.; Wu, J.; Xu, J.; Zhang, K.; Duan, J.; Lei, J.; Wu, J.; et al. Systematic comparison on convective heat transfer characteristics of several pin fins for turbine cooling. Crystals 2021, 11, 977. [Google Scholar] [CrossRef]
  26. Ye, B.; Rubel, M.R.H.; Li, H. Design and optimization of cooling plate for battery module of an electric vehicle. Appl. Sci. 2019, 9, 754. [Google Scholar] [CrossRef] [Green Version]
  27. Zhang, J.; Zhang, T.; Prakash, S.; Jaluria, Y. Experimental and numerical study of transient electronic chip cooling by liquid flow in microchannel heat sinks. Numer. Heat Transf. Part A Appl. 2014, 65, 627–643. [Google Scholar] [CrossRef]
  28. Annaratone, D. Transient heat transfer. SpringerBriefs Appl. Sci. Technol. 2011, 1–46. [Google Scholar] [CrossRef]
  29. Incropera, F.; Dewitt, D.; Bergman, T.; Lavine, A. Fundamentals of Heat and Mass Transfer, 6th ed.; John Wiley & Sons, Inc.: Hoboken, NJ, USA, 2007; ISBN 978-0-471-45728-2. [Google Scholar]
  30. Shevchuk, I.V. Unsteady conjugate laminar heat transfer of a rotating non-Uniformly heated disk: Application to the transient experimental technique. Int. J. Heat Mass Transf. 2006, 49, 3530–3537. [Google Scholar] [CrossRef]
  31. Indinger, T.; Shevchuk, I.V. Transient laminar conjugate heat transfer of a rotating disk: Theory and numerical simulations. Int. J. Heat Mass Transf. 2004, 47, 3577–3581. [Google Scholar] [CrossRef]
  32. Lienhard, J.H.; Lienhard, J.H. A Heat Transfer Textbook, 5th ed.; Phlogiston Press: Cambridge, MA, USA, 2020; ISBN 0486837351. [Google Scholar]
  33. Li, F.; Ning, J.; Liang, S.Y. Analytical modeling of the temperature using uniform moving heat source in planar induction heating process. Appl. Sci. 2019, 9, 1445. [Google Scholar] [CrossRef] [Green Version]
  34. Llano-Sánchez, L.E.; Domínguez-Cajeli, D.M.; Ruiz-Cárdenas, L.C. Thermal transfer analysis of tubes with extended surface with fractal design. Rev. Fac. Ing. 2018, 27, 31–37. [Google Scholar] [CrossRef] [Green Version]
  35. Levin, P.P.; Shitzer, A.; Hetsroni, G. Numerical optimization af a PCM-Based heat sink with internal fins. Int. J. Heat Mass Transf. 2013, 61, 638–645. [Google Scholar] [CrossRef]
  36. Sinha, K.K.; Ahirwar, P. Computational analysis of heat sink or extended surface. Int. Res. J. Eng. Technol. 2017, 4, 70–79. [Google Scholar]
  37. Rincón Tabares, J.S.; Perdomo-Hurtado, L.; Aragón, J.L. Study of gasketed-Plate heat exchanger performance based on energy efficiency indexes. Appl. Therm. Eng. 2019, 159, 113902. [Google Scholar] [CrossRef]
  38. Rahimi-Gorji, M.; Pourmehran, O.; Hatami, M.; Ganji, D.D. Statistical optimization of microchannel heat sink (mchs) geometry cooled by different nanofluids using RSM analysis. Eur. Phys. J. Plus 2015, 130, 1–21. [Google Scholar] [CrossRef]
  39. Viguerie, A.; Auricchio, F. Numerical solution of additive manufacturing problems using a two-Level method. Int. J. Numer. Methods Eng. 2021. [Google Scholar] [CrossRef]
  40. Galati, M.; Di Mauro, O.; Iuliano, L. Finite element simulation of multilayer electron beam melting for the improvement of build quality. Crystals 2020, 10, 532. [Google Scholar] [CrossRef]
  41. Fan, H.; Wang, R.; Xu, Z.; Duan, H.; Chen, D. Migration and enrichment behaviors of Ca and Mg elements during cooling and crystallization of boron-Bearing titanium slag melt. Crystals 2021, 11, 888. [Google Scholar] [CrossRef]
  42. Feng, Z.; Wen, Z.; Lu, G.; Zhao, Y. Influence of cooling scenarios on the evolution of microstructures in nickel-Based single crystal superalloys. Crystals 2022, 12, 74. [Google Scholar] [CrossRef]
  43. Kawaguchi, G.; Yamamoto, H.M. Control of organic superconducting field-Effect transistor by cooling rate. Crystals 2019, 9, 605. [Google Scholar] [CrossRef] [Green Version]
  44. Razavi, S.; Gupta, H.V. What do we mean by sensitivity analysis? The need for comprehensive characterization of “global” Sensitivity in earth and environmental systems models. Water Resour. Res. 2015, 51, 3070–3092. [Google Scholar] [CrossRef]
  45. Özisik, M.N.; Orlande, H.R.B. Inverse Heat Transfer; Routledge: London, UK, 2018. [Google Scholar]
  46. Du, W.F.; Zhuo, S.; Yang, L.; Zhao, R.C. Numerical simulation and parameter sensitivity analysis of coupled heat transfer by PCCS containment wall. Appl. Therm. Eng. 2017, 113, 867–877. [Google Scholar] [CrossRef]
  47. Wang, C.; Matthies, H.G. Non-Probabilistic interval process model and method for uncertainty analysis of transient heat transfer problem. Int. J. Therm. Sci. 2019, 144, 147–157. [Google Scholar] [CrossRef]
  48. Mirkoohi, E.; Ning, J.; Bocchini, P.; Fergani, O.; Chiang, K.-N.; Liang, S. Thermal modeling of temperature distribution in metal additive manufacturing considering effects of build layers, latent heat, and temperature-Sensitivity of material properties. J. Manuf. Mater. Process. 2018, 2, 63. [Google Scholar] [CrossRef] [Green Version]
  49. Lane, B.; Heigel, J.; Ricker, R.; Zhirnov, I.; Khromschenko, V.; Weaver, J.; Phan, T.; Stoudt, M.; Mekhontsev, S.; Levine, L. Measurements of melt pool geometry and cooling rates of individual laser traces on IN625 bare plates. Integr. Mater. Manuf. Innov. 2020, 9, 16–30. [Google Scholar] [CrossRef]
  50. Milani Shirvan, K.; Ellahi, R.; Mirzakhanlari, S.; Mamourian, M. Enhancement of heat transfer and heat exchanger effectiveness in a double pipe heat exchanger filled with porous media: Numerical simulation and sensitivity analysis of turbulent fluid flow. Appl. Therm. Eng. 2016, 109, 761–774. [Google Scholar] [CrossRef]
  51. Martins, J.R.R.A.; Hwang, J.T. Review and Unification of Methods for Computing Derivatives of Multidisciplinary Computational Models. Am. Inst. Aeronaut. Astronaut. 2013, 51, 2582–2599. [Google Scholar] [CrossRef] [Green Version]
  52. Baydin, A.G.; Pearlmutter, B.A.; Radul, A.A.; Siskind, J.M. Automatic differentiation in machine learning: A survey. J. Mach. Learn. Res. 2015, 18, 1–43. [Google Scholar]
  53. Jayapragasam, P.; Le Bideau, P.; Loulou, T. Computing sensitivity coefficients by using complex differentiation: Application to heat conduction problem. Numer. Heat Transf. Part B Fundam. 2018, 74, 729–745. [Google Scholar] [CrossRef]
  54. Margossian, C.C. A review of automatic differentiation and its efficient implementation. Wiley Interdiscip. Rev. Data Min. Knowl. Discov. 2019, 9, e1305. [Google Scholar] [CrossRef] [Green Version]
  55. Garza, J.E.; Millwater, H.R. Sensitivity Analysis in Structural Dynamics using the ZFEM Complex Variable Finite Element Method. In Proceedings of the 54th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, Boston, MA, USA, 8–11 April 2013; American Institute of Aeronautics and Astronautics: Reston, VA, USA, 2013; pp. 1–16. [Google Scholar] [CrossRef] [Green Version]
  56. Squire, W.; Trapp, G. Using complex variables to estimate derivatives of real functions. SIAM Rev. 1998, 40, 110–112. [Google Scholar] [CrossRef] [Green Version]
  57. Lyness, J.N.; Moler, C.B. Numerical differentiation of analytic functions. SIAM J. Numer. Anal. 1967, 4, 202–210. [Google Scholar] [CrossRef]
  58. Cui, M.; Gao, X.; Chen, H. A new inverse approach for the equivalent gray radiative property of a non-Gray medium using a modified zonal method and the complex-Variable-differentiation method. J. Quant. Spectrosc. Radiat. Transf. 2011, 112, 1336–1342. [Google Scholar] [CrossRef]
  59. Gao, X.W.; He, M.C. A new inverse analysis approach for multi-region heat conduction bem using complex-Variable-Differentiation method. Eng. Anal. Bound. Elem. 2005, 29, 788–795. [Google Scholar] [CrossRef]
  60. Voorhees, A.; Millwater, H.; Bagley, R. Complex variable methods for shape sensitivity of finite element models. Finite Elem. Anal. Des. 2011, 47, 1146–1156. [Google Scholar] [CrossRef]
  61. Millwater, H.; Wagner, D.; Baines, A.; Lovelady, K. Improved WCTSE method for the generation of 2D weight functions through implementation into a commercial finite element code. Eng. Fract. Mech. 2013, 109, 302–309. [Google Scholar] [CrossRef]
  62. Jaluria, Y. Solution of inverse problems in thermal systems. J. Therm. Sci. Eng. Appl. 2020, 12. [Google Scholar] [CrossRef]
  63. Daouas, N. An Alternative sensitivity method for a two-Dimensional inverse heat conduction-radiation problem based on transient hot-Wire measurements. Numer. Heat Transf. Part B Fundam. 2018, 73, 106–128. [Google Scholar] [CrossRef]
  64. Zhou, Z.F.; Xu, T.Y.; Chen, B. Algorithms for the estimation of transient surface heat flux during ultra-Fast surface cooling. Int. J. Heat Mass Transf. 2016, 100, 1–10. [Google Scholar] [CrossRef]
  65. Cui, M.; Gao, X.; Zhang, J. A new approach for the estimation of temperature-Dependent thermal properties by solving transient inverse heat conduction problems. Int. J. Therm. Sci. 2012, 58, 113–119. [Google Scholar] [CrossRef]
  66. Cui, M.; Yang, K.; Xu, X.L.; Wang, S.D.; Gao, X.W. A modified levenberg-Marquardt algorithm for simultaneous estimation of multi-Parameters of boundary heat flux by solving transient nonlinear inverse heat conduction problems. Int. J. Heat Mass Transf. 2016, 97, 908–916. [Google Scholar] [CrossRef]
  67. Cui, M.; Li, N.; Liu, Y.; Gao, X. Robust Inverse approach for two-Dimensional transient nonlinear heat conduction problems. J. Thermophys. Heat Transf. 2015, 29, 253–262. [Google Scholar] [CrossRef]
  68. Voorhees, A.; Bagley, R.; Millwater, H.; Golden, P. Application of Complex Variable Methods for Fatigue Sensitivity Analysis. In Proceedings of the 50th AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, Palm Springs, CA, USA, 4–7 May 2009; American Institute of Aeronautics and Astronautics: Reston, VA, USA, 2009; p. 2711. [Google Scholar]
  69. Voorhees, A.; Millwater, H.; Bagley, R.; Golden, P. Fatigue sensitivity analysis using complex variable methods. Int. J. Fatigue 2012, 40, 61–73. [Google Scholar] [CrossRef]
  70. Wagner, D.; Millwater, H. 2D weight function development using a complex taylor series expansion method. Eng. Fract. Mech. 2012, 86, 23–37. [Google Scholar] [CrossRef]
  71. Millwater, H.; Wagner, D.; Baines, A.; Montoya, A. A virtual crack extension method to compute energy release rates using a complex variable finite element method. Eng. Fract. Mech. 2016, 162, 95–111. [Google Scholar] [CrossRef] [Green Version]
  72. Wagner, D.; Garcia, M.J.; Montoya, A.; Millwater, H. A finite element-Based adaptive energy response function method for 2D curvilinear progressive fracture. Int. J. Fatigue 2019, 127, 229–245. [Google Scholar] [CrossRef]
  73. Ramirez-Tamayo, D.; Soulami, A.; Gupta, V.; Restrepo, D.; Montoya, A.; Millwater, H. A complex-Variable cohesive finite element subroutine to enable efficient determination of interfacial cohesive material parameters. Eng. Fract. Mech. 2021, 247, 107638. [Google Scholar] [CrossRef]
  74. Montoya, A.; Fielder, R.; Gomez-Farias, A.; Millwater, H. Finite-Element sensitivity for plasticity using complex variable methods. J. Eng. Mech. 2015, 141, 04014118. [Google Scholar] [CrossRef]
  75. Fielder, R.; Montoya, A.; Millwater, H.; Golden, P. Residual stress sensitivity analysis using a complex variable finite element method. Int. J. Mech. Sci. 2017, 133, 112–120. [Google Scholar] [CrossRef]
  76. Gomez-Farias, A.; Montoya, A.; Millwater, H. Complex finite element sensitivity method for creep analysis. Int. J. Press. Vessel. Pip. 2015, 132–133, 27–42. [Google Scholar] [CrossRef]
  77. Fielder, R.; Millwater, H.; Montoya, A.; Golden, P. Efficient estimate of residual stress variance using complex variable finite element methods. Int. J. Press. Vessel. Pip. 2019, 173, 101–113. [Google Scholar] [CrossRef]
  78. Chun, J. Reliability-Based design optimization of structures using complex-Step approximation with sensitivity analysis. Appl. Sci. 2021, 11, 4708. [Google Scholar] [CrossRef]
  79. Montoya, A.; Millwater, H. Sensitivity analysis in thermoelastic problems using the complex finite element method. J. Therm. Stress. 2017, 40, 302–321. [Google Scholar] [CrossRef]
  80. Garza, J.; Millwater, H. Multicomplex newmark-Beta time integration method for sensitivity analysis in structural dynamics. AIAA J. 2015, 53, 1188–1198. [Google Scholar] [CrossRef]
  81. Navarro, J.D.; Millwater, H.R.; Montoya, A.H.; Restrepo, D. Arbitrary-Order sensitivity analysis in phononic metamaterials using the multicomplex taylor series expansion method coupled with bloch’s theorem. J. Appl. Mech. 2022, 89, 1–43. [Google Scholar] [CrossRef]
  82. Taler, J.; Ocłoń, P. Finite element method in steady-State and transient heat conduction. In Encyclopedia of Thermal Stresses; Springer: Dordrecht, The Netherlands, 2014; pp. 1604–1633. [Google Scholar]
  83. Logan, D.L. A First Course in the Finite Element Method; Cengage Learning: Boston, MA, USA, 2011; ISBN 9781133169055. [Google Scholar]
  84. Morville, S.; Carin, M.; Peyre, P.; Gharbi, M.; Carron, D.; Le Masson, P.; Fabbro, R. 2D longitudinal modeling of heat transfer and fluid flow during multilayered direct laser metal deposition process. J. Laser Appl. 2012, 24, 032008. [Google Scholar] [CrossRef]
  85. de La Batut, B.; Fergani, O.; Brotan, V.; Bambach, M.; El Mansouri, M. Analytical and numerical temperature prediction in direct metal deposition of Ti6Al4V. J. Manuf. Mater. Process. 2017, 1, 3. [Google Scholar] [CrossRef] [Green Version]
  86. Ramirez Tamayo, D.; Montoya, A.; Millwater, H. Complex-Variable finite-Element method for mixed mode fracture and interface cracks. AIAA J. 2018, 56, 4632–4637. [Google Scholar] [CrossRef]
  87. Montoya, A.; Ramirez-Tamayo, D.; Millwater, H.; Kirby, M. A complex-Variable virtual crack extension finite element method for elastic-Plastic fracture mechanics. Eng. Fract. Mech. 2018, 202, 242–258. [Google Scholar] [CrossRef]
  88. Price, G.B. An Introduction to Multicomplex Spaces and Functions, 1st ed.; CRC Press: Boca Raton, FL, USA, 1991. [Google Scholar]
  89. Fike, J.; Alonso, J. The Development of Hyper-Dual Numbers for Exact Second-Derivative Calculations. In Proceedings of the 49th AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition, Orlando, FL, USA, 4–7 January 2011; American Institute of Aeronautics and Astronautics (AIAA): Reston, VA, USA, 2011; pp. 1–17. [Google Scholar]
  90. Parry, L.; Ashcroft, I.A.; Wildman, R.D. Understanding the effect of laser scan strategy on residual stress in selective laser melting through thermo-Mechanical simulation. Addit. Manuf. 2016, 12, 1–15. [Google Scholar] [CrossRef] [Green Version]
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.
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.
Applsci 12 02738 g001
Figure 2. (a) Mesh perturbation function P ( x ) for nodal coordinates required for geometrical sensitivities using Transient Thermal ZFEM, where p represents a fraction of the mesh that was perturbed. (b) Duplicated node approach for the CR representation for a 2D finite element. (c) Detail of the local coordinate system and quadrature points for the 4 real 4 imaginary nodes 2D element.
Figure 2. (a) Mesh perturbation function P ( x ) for nodal coordinates required for geometrical sensitivities using Transient Thermal ZFEM, where p represents a fraction of the mesh that was perturbed. (b) Duplicated node approach for the CR representation for a 2D finite element. (c) Detail of the local coordinate system and quadrature points for the 4 real 4 imaginary nodes 2D element.
Applsci 12 02738 g002
Figure 3. Flowchart of Transient Thermal ZFEM and Abaqus User Element Subroutine. All complex operations are performed inside the UEL and the A matrix and the RHS vector are decomposed into Cauchy-Riemann form before exiting the UEL.
Figure 3. Flowchart of Transient Thermal ZFEM and Abaqus User Element Subroutine. All complex operations are performed inside the UEL and the A matrix and the RHS vector are decomposed into Cauchy-Riemann form before exiting the UEL.
Applsci 12 02738 g003
Figure 4. (a) Convergence analysis for the heat transfer coefficient per unit of thickness where the optimal value that minimizes the NRMSE is 735 ( W / ( m 2 K ) ) / m . (b) The dimensionless temperature history of fin tip: solid black line, Equation (23), yellow star, the real part of Transient Thermal ZFEM solution, blue square, Abaqus built-in element solution, thin red dashed line, steady-state solution (Equation (24)).
Figure 4. (a) Convergence analysis for the heat transfer coefficient per unit of thickness where the optimal value that minimizes the NRMSE is 735 ( W / ( m 2 K ) ) / m . (b) The dimensionless temperature history of fin tip: solid black line, Equation (23), yellow star, the real part of Transient Thermal ZFEM solution, blue square, Abaqus built-in element solution, thin red dashed line, steady-state solution (Equation (24)).
Applsci 12 02738 g004
Figure 5. (a) Step size convergence analysis comparing Transient Thermal ZFEM and FD for the UEL and the selected analyzed variables. (b) Convergence study for mesh perturbation fraction λ . The insets show P ( x ) for the rule of thumb of the literature and λ = 1 in the top left and bottom right corners, respectively.
Figure 5. (a) Step size convergence analysis comparing Transient Thermal ZFEM and FD for the UEL and the selected analyzed variables. (b) Convergence study for mesh perturbation fraction λ . The insets show P ( x ) for the rule of thumb of the literature and λ = 1 in the top left and bottom right corners, respectively.
Applsci 12 02738 g005
Figure 6. Dimensionless temperature sensitivity history at the fin tip to (a) thermal conductivity; (b) room temperature; (c) fin length. Solid black line, analytical sensitivity using analytical differentiation in Equation (23), yellow star, sensitivity using Transient Thermal ZFEM, blue square, sensitivity using 1% FD and Abaqus built-in element solutions, thin red dashed line, sensitivity using analytical differentiation for steady-state using Equation (24).
Figure 6. Dimensionless temperature sensitivity history at the fin tip to (a) thermal conductivity; (b) room temperature; (c) fin length. Solid black line, analytical sensitivity using analytical differentiation in Equation (23), yellow star, sensitivity using Transient Thermal ZFEM, blue square, sensitivity using 1% FD and Abaqus built-in element solutions, thin red dashed line, sensitivity using analytical differentiation for steady-state using Equation (24).
Applsci 12 02738 g006
Figure 7. The heating rate in time for fin tip: solid black line, rate using analytical differentiation in Equation (23), yellow star, rate using Transient Thermal ZFEM, blue square, FD using Abaqus built-in element solution.
Figure 7. The heating rate in time for fin tip: solid black line, rate using analytical differentiation in Equation (23), yellow star, rate using Transient Thermal ZFEM, blue square, FD using Abaqus built-in element solution.
Applsci 12 02738 g007
Figure 8. (a) Temperature-dependent specific heat of Ti64. Experimental data obtained from [90]. The inset shows the original and the perturbed property in the T C P i space; (b) dimensionless temperature history of fin tip for temperature-dependent specific heat; (c) sensitivity of θ to specific temperature-dependent specific heat. Yellow star, sensitivity using Transient Thermal ZFEM, blue square, sensitivity using 1 % FD and Abaqus built-in element solution, black dashed line, steady-state solution.
Figure 8. (a) Temperature-dependent specific heat of Ti64. Experimental data obtained from [90]. The inset shows the original and the perturbed property in the T C P i space; (b) dimensionless temperature history of fin tip for temperature-dependent specific heat; (c) sensitivity of θ to specific temperature-dependent specific heat. Yellow star, sensitivity using Transient Thermal ZFEM, blue square, sensitivity using 1 % FD and Abaqus built-in element solution, black dashed line, steady-state solution.
Applsci 12 02738 g008
Figure 9. (a) Time history of fin ε , (b) time history of fin η . Black line, Equation (13), yellow star, Transient Thermal ZFEM solution, blue square, Abaqus built-in element solution, black thin dashed line, black dashed line, the steady-state solution.
Figure 9. (a) Time history of fin ε , (b) time history of fin η . Black line, Equation (13), yellow star, Transient Thermal ZFEM solution, blue square, Abaqus built-in element solution, black thin dashed line, black dashed line, the steady-state solution.
Applsci 12 02738 g009
Figure 10. Thermal effectiveness sensitivity history in fin tip to (a) thermal conductivity; (b) room temperature; (c) fin length. Solid black line, analytical sensitivity using analytical differentiation in Equation (26), yellow star, sensitivity using Transient Thermal ZFEM, blue square, sensitivity using 1 % FD and Abaqus built-in element solutions, thin red dashed line, sensitivity using analytical differentiation for steady-state using Equation (26). (d) Steady-state influence of fin parameters in ε based on relative sensitivities. In the pie chart, the missing parameters ( b , ρ , C P , T , T 0 ) have an influence of less than 1 % compared to the sum of all relative sensitivities. Moreover, the ( ) in the pie sections implies an inverse relationship between the output variable and the input parameter, while a ( + ) parameter has a direct relationship with the output.
Figure 10. Thermal effectiveness sensitivity history in fin tip to (a) thermal conductivity; (b) room temperature; (c) fin length. Solid black line, analytical sensitivity using analytical differentiation in Equation (26), yellow star, sensitivity using Transient Thermal ZFEM, blue square, sensitivity using 1 % FD and Abaqus built-in element solutions, thin red dashed line, sensitivity using analytical differentiation for steady-state using Equation (26). (d) Steady-state influence of fin parameters in ε based on relative sensitivities. In the pie chart, the missing parameters ( b , ρ , C P , T , T 0 ) have an influence of less than 1 % compared to the sum of all relative sensitivities. Moreover, the ( ) in the pie sections implies an inverse relationship between the output variable and the input parameter, while a ( + ) parameter has a direct relationship with the output.
Applsci 12 02738 g010
Figure 11. Thermal efficiency sensitivity history in fin tip to (a) thermal conductivity; (b) room temperature; (c) fin length. Solid black line, analytical sensitivity using analytical differentiation in Equation (27), yellow star, sensitivity using Transient Thermal ZFEM, blue square, sensitivity using 1 % FD and Abaqus built-in element solutions, thin red dashed line, sensitivity using analytical differentiation for steady-state using Equation (27). (d) Steady-state influence of fin parameters in η based on relative sensitivities. In the pie chart, the missing parameters ( ρ , C P , T , T 0 ) have an influence of less than 1 % compared to the sum of all relative sensitivities. Moreover, the ( ) in the pie sections implies an inverse relationship between the output variable and the input parameter while a ( + ) parameter has the opposite relationship with the output.
Figure 11. Thermal efficiency sensitivity history in fin tip to (a) thermal conductivity; (b) room temperature; (c) fin length. Solid black line, analytical sensitivity using analytical differentiation in Equation (27), yellow star, sensitivity using Transient Thermal ZFEM, blue square, sensitivity using 1 % FD and Abaqus built-in element solutions, thin red dashed line, sensitivity using analytical differentiation for steady-state using Equation (27). (d) Steady-state influence of fin parameters in η based on relative sensitivities. In the pie chart, the missing parameters ( ρ , C P , T , T 0 ) have an influence of less than 1 % compared to the sum of all relative sensitivities. Moreover, the ( ) in the pie sections implies an inverse relationship between the output variable and the input parameter while a ( + ) parameter has the opposite relationship with the output.
Applsci 12 02738 g011
Table 1. Parameters involved in the transient thermal fin problem.
Table 1. Parameters involved in the transient thermal fin problem.
Parameter TypeParameter of Interest ( ϕ ) ValueUnits
Material Properties k 7.1 W / ( m K )
C p 580 J / ( kg K )
ρ 4430 kg / m 3
Boundary Conditions h U 114 W / ( m 2 K )
T 283 K
T w 389 K
Initial Conditions T 0 283 K
Geometry b 51 mm
L 1 m
δ 4.75 mm
Table 2. NRMSE of θ and the selected sensitivities compared to the analytical solution.
Table 2. NRMSE of θ and the selected sensitivities compared to the analytical solution.
NRMSE
θ S k θ T S b θ
ZFEM0.0030.0110.0020.003
FD0.0030.0090.0090.019
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Rincon-Tabares, J.-S.; Velasquez-Gonzalez, J.C.; Ramirez-Tamayo, D.; Montoya, A.; Millwater, H.; Restrepo, D. Sensitivity Analysis for Transient Thermal Problems Using the Complex-Variable Finite Element Method. Appl. Sci. 2022, 12, 2738. https://doi.org/10.3390/app12052738

AMA Style

Rincon-Tabares J-S, Velasquez-Gonzalez JC, Ramirez-Tamayo D, Montoya A, Millwater H, Restrepo D. Sensitivity Analysis for Transient Thermal Problems Using the Complex-Variable Finite Element Method. Applied Sciences. 2022; 12(5):2738. https://doi.org/10.3390/app12052738

Chicago/Turabian Style

Rincon-Tabares, Juan-Sebastian, Juan C. Velasquez-Gonzalez, Daniel Ramirez-Tamayo, Arturo Montoya, Harry Millwater, and David Restrepo. 2022. "Sensitivity Analysis for Transient Thermal Problems Using the Complex-Variable Finite Element Method" Applied Sciences 12, no. 5: 2738. https://doi.org/10.3390/app12052738

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop