Transient Thermal Spreading From a Circular Heat Source in Polygonal Flux Tubes

The ability to characterize thermal spreading resistance is an important factor for many applications, including aeronautics [1], tribology, micro-electronics, and marine icing [2] to name a few. Early works included experimental [3] and analytical studies [4]. The understanding of thermal spreading resistance has evolved through the examination of various aspects, such as transient spreading resistance, the influence of geometry, and the role of numerical modeling. All of the historical work can be found in the text by Muzychka and Yovanovich [5].

In micro-electronics and other fields, understanding and optimizing the thermal spreading resistance is crucial for effective thermal management of micro-electronic devices [6]. This resistance can be influenced by factors such as material properties, thickness, geometry of the materials involved, and temperature difference between the two surfaces. Despite its significance, many aspects of thermal spreading resistance have not yet been solved, leading to a focus on understanding this phenomenon under various physical and geometrical conditions in recent studies.

This study numerically analyzed the transient thermal spreading in polygonal flux tubes and the impact of different geometric conditions on the heat transfer process (see Fig. 1). The study of transient thermal spreading from a circular heat source in a polygonal flux channel involves the analysis of the conduction of heat from a circular heat source, such as a device, droplet, or component, to the surrounding material through an array of flux channels with polygonal cross-sectional shapes and symmetry boundary conditions on the lateral edges. The polygonal flux tube results from the arrangement of the heat sources on a larger domain as shown in Fig. 2.

Several researchers have studied the influence of geometry and edge cooling on the thermal spreading resistance. Fried [1] examined the effect of surface roughness, normal load, and temperature on the thermal contact conductance between two surfaces in contact, and found that the thermal contact conductance could be improved by reducing the surface roughness and increasing the normal load. Mikic et al. [7] analyzed the heat transfer between a droplet and solid surface during dropwise condensation, noting that the heat transfer was influenced by the thermal spreading resistance at the interface and the physical properties of the droplet and surface.

Yovonovich et al. [8] introduced an elaborate solution for a thermal system consisting of a circular heat source in perfect contact with a compound circular disk. The compound disk, comprised of two different isotropic materials each with unique thermal conductivities and thicknesses, is perfectly aligned. This disk's free surfaces are adiabatic, and the heat flux in the heat source area is either uniform or of specific patterns. Moreover, the bottom surface of the compound disk is consistently cooled through either uniform convective or contact conductance.

In their substantial contribution to the field, Nelson and Sayers [9] provided both numerical and analytical solutions for a scenario involving an isoflux circular source. Their study offered mathematical expressions for determining area-average and peak temperatures. Additionally, the study also introduced a concise, closed-form equation that achieves an accuracy level within 10% of the comprehensive solution. In a pair of associated studies [10] and [11], the authors introduced analytical solutions for the isoflux circular source on a finite thermal spreader. Song provided formulas for calculating both average area and peak temperatures, he also introduced a simple yet accurate expression, demonstrating precision within a 10% margin of the complete solution. Complementing this, Lee's research presents an analytical model specifically tailored to predict these resistances under varying cooling conditions in electronic components. The model offers a comprehensive view by providing dimensionless resistances across a spectrum of relative source sizes and inverse external resistances. They highlighted that their fully computed results aligned impressively well with the numerical values provided by Nelson and Sayers [9].

Additionally, numerous studies have investigated the influence of the Biot number, a key dimensionless parameter, on thermal spreading resistance. This number represents the ratio of internal to external thermal resistance within a body, crucial for assessing heat transfer dynamics. Muzychka et al. [12] underlined the crucial role of the Biot number in thermal spreading resistance modeling. They provide key insights into how this parameter influences the spreading parameter's behavior and the heat transfer characteristics within compound annular sectors. In a recent study, Jain [13] discusses the thermal spreading resistance from an isothermal source into a finite-thickness body, which is relevant for the design of heat sinks and thermal spreaders for micro-electronics and other heat-generating devices. The study highlights the importance of selecting a suitably high Biot number, like Bi_max = 1000, to accurately represent isothermal conditions in mixed boundary scenarios.

Bagnall et al. [14] developed an analytical solution for temperature rise in complex multilayer structures with discrete heat sources, analyzing the temperature distribution in a multilayer structure with a circular heat source and considering the effect of thermal resistance at the interfaces. Yovanovich et al. [15] presented a general expression for the spreading resistance of an isoflux rectangular heat source on a two-layer rectangular flux channel with convective or conductive cooling at one boundary. They also discussed the effect of heat flux distribution over strip sources on two-dimensional spreading resistances and provided tabulated values for the three flux distributions. Muzychka et al. [16] provided a detailed analysis of the thermal spreading resistance in multilayered contacts and its applications in thermal contact resistance.

In a separate study, Muzychka et al. [17] analyzed the thermal spreading resistance in isotropic and compound rectangular flux channels using a solution for an isotropic or compound circular flux tube. The research found that the circular disk solution, which requires only a single series summation, was more efficient than the solution for rectangular flux channels. Furthermore, the authors presented a new analytical solution that considered various limiting cases, including adiabatic edges, to examine the thermal spreading resistance in rectangular flux channels with edge cooling. This approach was proven to be effective for a broad range of dimensionless parameters, as shown in Ref. [18]. Bhatt and Rhee [19] also used numerical modeling to analyze thermal spreading resistance in square and rectangular regions and found that the results matched well with analytical solutions.

Muzychka et al. [20] developed general expressions for determining the spreading resistance of an eccentric rectangular heat source on different types of flux channels and demonstrated that the solution for a central heat source could be used to calculate the spreading resistance of corner and edge heat sources using images. In a separate study, Muzychka et al. [21] reviewed thermal spreading resistance theory in compound and orthotropic systems, presented solutions for various compound and isotropic systems and developed new solutions for orthotropic systems. Negus and Yovanovich [22] analyzed the transient temperature rise at a surface due to arbitrary contacts on half-spaces, whereas Negus et al. [23] examined the quantification and normalization of constriction resistance in heat flux tubes. Turyk and Yovanovich [24] analyzed heat transfer through various types of flux channels, including circular and noncircular shapes, and the impact of different parameters on constriction resistance.

Muzychka et al. [25] developed new solutions for thermal spreading resistance in compound orthotropic systems that include interfacial resistance. These solutions can be applied to both circular disks and rectangular flux channels with a central uniform heat flux source, and represent novel extensions of previous solutions for other systems. In addition, Yovanovich [26] presented dimensionless spreading resistances for two types of isoflux-planar contact areas: regular polygons and a hyperellipse. The spreading resistance was based on the transient centroid temperature rise, and closed-form expressions for the steady-state spreading resistance were provided for regular polygons and hyperellipses as functions of their shape parameter and aspect ratios. Yovanovich et al. [27] developed a mathematical model for transient conduction from isothermal convex bodies of arbitrary shape.

The model combines short- and long-time asymptotic solutions, which correspond to the half-space solution and the conduction shape factor, respectively. They also presented the values for the blending parameter as a function of both the aspect ratio and body shape, which provided the ability to develop correlations using short-time steady-state asymptotes. In thermal spreading resistance studies of compound annular sectors, the blending parameter indicates the relative contact size and is essential for calculating spreading resistance. It also helps evaluate the impact of coatings or deposits on larger tubes in these sectors.

In practical applications, the shape of the contact area is often unpredictable owing to the random nature of the contacting surfaces. To address this challenge, a model that can be applied to any arbitrarily shaped contact area is highly desirable. Previous studies have provided insights into the thermal spreading resistance in flux channels, including the effects of edge cooling and the accuracy of using circular disk solutions. However, there is still a need for further understanding of transient thermal spreading in polygonal geometries from an isothermal or isoflux circular heat source.

Previous studies have shown that a universal time function, taking into account the specific geometry and material properties of the flux tube, can accurately describe the transient thermal spreading resistance. This function can be derived by analyzing the asymptotic characteristics of the solutions obtained from these studies.

where $Ts¯$ is the mean source temperature and T_i is the constant sink temperature.

where R_s is the thermal spreading resistance, $Ts¯$(t) is the temperature of the heat source at time t, $T¯cp$(t) is the contact plane temperature at time t, and Q is the heat flow from the source. This equation takes into account the semi-infinite length of the flux channel and the specific geometry and material properties of the tube in order to describe the transient thermal spreading process.

Spreading resistance problems generally address either half-space or flux channel geometries when the depth of the substrate is large compared to the width of the source. The length scales involved determine which is applicable, Muzychka and Yovanovich [5]. When the length scale of the source is small compared to the length scale of their spacing, the problem can be framed as a single contact transferring energy to a half-space [28]. However, when the length scale of the sources is on the same order of magnitude as the spacing between them, the effects of the neighboring sources must be considered. This is done by dividing the substrate into an array of adjacent flux tubes or channels, one for each source, with symmetry boundary conditions between them (see Fig. 1). The shape of the flux tube or channel depends on the geometry of the source array. Figure 2 shows some possible configurations. It incorporates various factors that affect the heat transfer process, such as the material properties of the flux tubes, the geometry of the tubes, the temperature difference between the heat source and the surrounding material, and the constriction ratio of circular heat source to flux tube or channel area. Figure 1 illustrates the simulation domain of a hexagonal flux tube of infinite length with a circular heat source at the center of one end of the flux tube. For all polygonal geometries the constriction ratio (ϵ) is equal to $ϵ=AsAt$⁠. The initial temperature in the system is uniform at $Ti$⁠. At time t > 0, a constant temperature boundary condition is applied at the surface of the heat source region, which is located at the center of the hexagonal flux tube. The surface outside of the heat source region is assumed to be adiabatic.

for an isoflux heat source.

where $n$ is an inward directed normal.

The temperature distribution in the flux tubes was modeled using a finite volume method to solve the heat diffusion equation and boundary conditions. The time-derivative was calculated using the Euler scheme. The simulation was performed using the OpenFOAM framework [29], and in particular, the discretization scheme was based on standard Gaussian finite volume integration, where cell face values are calculated by summing interpolated values from cell centers. For this model, the standard finite volume discretization of Gaussian integration is chosen, which requires linear interpolation or central differencing from cell centers to face centers. In addition, the residual control was set to enforce an additional condition for the solution's accuracy. In this case, the simulation would continue until the residual for “T”—the difference between successive iterative solutions-fell below $10−5$⁠. This control aids in ensuring that the numerical solution converges to a sufficiently accurate result. The length of the flux channel was assumed to be five times the radius, which is a reasonable assumption for this model.

To confirm the robustness of the results and ensure that they were not dependent on the specific mesh used, a mesh independence study was conducted. Sensitivity analyses were performed using meshes with varying resolutions, and the results were compared to ensure consistency and independence of the mesh for all geometries. This confirms that the numerical solutions are robust and can be applied to other similar problems with confidence. The use of different meshes also helped to identify the optimal mesh resolution for this problem, providing an efficient numerical solution.

The constriction ratio of $ϵ=0.4$ was considered in the analysis, and the dimensionless thermal resistance was used as the comparison parameter. Fourier number of Fo = 0.0005 was considered and the triangular flux channel was studied for mesh independence as it was assumed to have a higher error compared to the circular tube flux channel.

Based on this mesh independence analysis, a mesh resolution of 1,100,000 as shown in Fig. 3 was chosen for the remainder of this study. To determine an appropriate time-step size, a time independence study was conducted using the same geometry. During the study, simulations were run for $ϵ=0.2$ with the time-step size being varied in each simulation. The results of the study, as shown in Fig. 4, suggest that a time-step size of 0.0001 s is suitable for this analysis.

Since we are investigating the effect of flux tube geometry, particularly for larger constriction ratios than the half-space limits, we will also compare the results with established solutions for the circular flux tube. For both an isothermal and isolfux heat source, simple models are reported in Muzychka and Yovanovich [30].

In the next section, we will show that the above models are able to predict the numerical results quite well for different flux channel shapes for both the isothermal and isoflux heat sources.

To validate numerical results, a comparison was made between our steady-state study and that of Negus et al. [23]. The latter study provided exact values of the constriction resistance parameters (ψ) for various contact area configurations, including circle-circle (ψ_cc), circle-square (ψ_cs), and square-square (ψ_ss).

In this study, we validated our results for two relative contact sizes (ϵ = 0.4, 0.6) and for the circle on circle (ψ_cc) and circle on square (ψ_cs) configurations. The results are presented in Tables 1 and 2, respectively.

To further validate our findings, we compared our results to those of a study by Yovanovich on the transient spreading resistance of arbitrary isoflux contact areas [26]. This study provides a comprehensive analysis of the transient behavior of heat spreading from a point source through various flux tubes, including the development of a universal short-time asymptote based on the solution for the isoflux circular contact area, which has an analytical solution for all dimensionless times. We used this study as a reference to confirm the validity of our own transient results by comparing our numerical results for a triangular flux tube with those presented in Yovanovich's study [26], as illustrated in Fig. 5. The comparison between our numerical results for a triangular flux tube and those presented in Yovanovich's study [26] showed a high degree of correlation, indicating the reliability and accuracy of our transient results. Overall, the validation process has confirmed the accuracy and reliability of our results for both transient and steady parts.

Figure 6 provides a top view of the isothermal lines for various polygonal flux tubes, offering a visual insight into the nuanced differences in thermal distribution among these geometric configurations. This figure is based on nondimensional temperature. Comparing all the representations, the circular cylinder distinctly exhibits the most uniform thermal distribution, as evidenced by its evenly spaced isothermal lines. As we transition from shapes with fewer sides, like triangles, to those with more, such as hexagons, the thermal distribution within these polygonal cylinders approaches the uniformity observed in the circular case. It's noteworthy that the triangle, with its sharp corners, diverges most from the ideal concentric isothermal pattern inherent to the circle.

Building upon the insights gained from Figs. 6 and 7 delves deeper into a more specific aspect: the comparative analysis of penetration depth. This analysis focuses on hexagonal and cylindrical flux channels at a constriction ratio of ϵ = 0.1 across three distinct Fourier numbers: $Fo=0.01,0.1,$ and 10. As the Fourier number increases, there's a marked enhancement in the penetration depth for both geometries, signifying a more profound heat penetration. The trend aligns with the larger Fourier numbers indicating longer durations and thus, more substantial heat diffusion. More importantly, Fig. 7 illustrates the subtle congruence in penetration depth between hexagonal and cylindrical flux channels at this constriction ratio across all Fourier numbers. It reveals that in terms of thermal performance, the two geometries have closely related behavior, emphasizing the flexibility of design choices, i.e., the layout of the heat sources.

As depicted in Fig. 8, five different constriction ratios (ϵ = 0.01, 0.2, 0.4, 0.6, 0.8) were incorporated into the analysis. The solid lines within the figure are the simple analytical model for a circular flux tube, as highlighted in the research by Lam and Muzychka [31]. The triangular flux channel consistently demonstrates the most significant deviation when compared to the simple model and other geometries. This deviation underscores the unique thermal challenges posed by triangular configurations in transient heat transfer scenarios. This multifaceted approach not only reaffirms the observations from the earlier figures but also paves the way for a more comprehensive understanding of the thermal behavior across varying geometric structures.

Turyk and Yovanovich [24] delivered solutions pertaining to the corresponding isoflux problem in a circular flux tube. They formulated the transient spreading resistance and used the dimensionless scheme presented in Ref. [27]. This not only elucidated the complexity of the problem but also provided a clear perspective on how thermal spreading could be modeled under dynamic conditions. Furthermore, in Fig. 9, we examine five unique constriction ratios, namely, ϵ = 0.01, 0.2, 0.4, 0.6, 0.8. It is important to underline that the derived analytical solutions here are based on Turyk and Yovanovich [24] solutions for the isoflux condition. The simple model presented earlier agrees well with the current numerical results for other flux tube shapes.

The results of the current study demonstrate that the thermal spreading resistance values are dependent on both the geometry of the flux tube and the contact area (or constriction ratio). As the geometry of the polygonal flux tube transitions from a triangle to a circle, the thermal spreading resistance approaches the resistance value of a circle-on-circle configuration. As shown in Table 3, the circle, with its uniform structure, consistently registers the lowest RMSE (root-mean-square-errors) values. The RMSE is calculated by comparing the exact simulation values (representing the actual observed data) with the predicted values obtained from the analytical models for circle-on-circle configuration. For each constriction ratio, this calculation involves assessing ten distinct data points, each corresponding to one of ten Fourier numbers. As polygonal shapes evolve from a triangle to a hexagon, they tend to emulate the circle's efficiency, evidenced by the hexagon's RMSE nearing that of the circle. Particularly at a constriction ratio of ϵ = 0.2, the triangle's performance starkly contrasts with the circle, producing a modest 3.48% RMSE disparity. However, the hexagon provides reasonable closeness to the circle, with an RMSE difference consistently under 0.38%. The triangular flux tube is further spotlighted by its consistently high RMSE values across all constriction ratios that were considered. In summary, the circular flux tube emerges as the optimal shape for heat distribution, while the data also emphasizes the potential of polygons like the hexagon, especially as they geometrically approach a circular form. Table 4 provides similar results for the isoflux source case.

Thermal spreading resistance plays a crucial role in various thermal engineering problems that involve variations in temperature and cross-sectional area. In this paper, a numerical model has been developed to simulate the transient thermal spreading in polygonal flux tubes from a circular heat source for both an isothermal and isoflux heat source. This study revealed that the effect of flux tube shape on the thermal spreading resistance is small when the results are appropriately nondimensionalzed using the square root of the source area as a characteristic length and by defining the constriction ratio as $ϵ=As/At$⁠. As the number of sides increases, the polygon approaches a circular configuration, and the thermal resistance approaches the circle-on-circle configuration. Furthermore, the root-mean-square error is calculated and it was shown that circle on triangle configuration has the largest error as compared with the simple model for the circle-on-circle case. Overall, this study provides valuable insight into the thermal spreading resistance in polygonal flux channels which can be applied to various thermal engineering problems in contact heat transfer problems. Good agreement between the numerical results and the analytical solution for the circular source on a circular flux tube is also shown to be better than 5 percent on average.

• Natural Sciences and Engineering Research Council of Canada and Mitacs (Funder IDs: 10.13039/501100000038 and 10.13039/501100004489).

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.

A_s =

source area (m²)

A_t =

flux channel cross-sectional area (m²)

$A$ =

characteristic length (m)

k =

thermal conductivity (W m^{– 1} K^{– 1})

Q =

heat flow rate (W)

R_s =

thermal spreading resistance (K W^{– 1})

R_total =

total thermal spreading resistance (K W^{– 1})

T_cp =

contact plane temperature (K)

T_f =

aink temperature (K)

T_i =

initial temperature (K)

T_s =

source temperature (K)

Greek Symbols

α =

thermal diffusivity, $kρCp$ (m² s^–1)

ψ =

constriction resistance parameter, $4kaR$

ϵ =

constriction ratio, $AsAt$

Dimensionless Groups

$FoA$ =

Fourier number, $αtA2$

$RA*$ =

dimensionless resistance (mean basis), $kR¯A$

$R0*$ =

dimensionless spreading resistance (centroid basis), $kR0A$

$C∞T$ =

dimensionless steady-state thermal spreading resistance for an isothermal source

$C∞q$ =

dimensionless steady-state thermal spreading resistance for an isoflux source