On the Numerical Investigation of Natural-Convection Heat Sinks Across a Wide Range of Flow and Operating Conditions

Abstract

Many designs for natural-convection heat sinks and semi-empirical correlations have been proposed in the recent years, but they are only valid in a limited range of Elenbaas numbers $El$ and were mostly tested for laminar flows. To alleviate those limits, parametric studies with 2D and quasi-3D models were carried out, in ranges of Grashof numbers up to $1.55\times {10}^{11}$ and Elenbaas numbers up to $3.42\times {10}^7$. Ansys Fluent’s laminar, transition-SST, SST k-$\omega$ and k-$ϵ$ models were applied. In addition, when used in this valid range, i.e., mean Elenbaas numbers, with the simplified quasi-3D model, the transition-SST model could predict better results, overestimating the heat flux by 10 to $15\%$ compared to semi-empirical correlations. The 2D model was not deemed satisfying, regarding turbulence models. Consequently, a quasi-3D model was developed: it appeared to be an efficient trade-off between computational time and prediction accuracy, in particular for turbulence models. New grouping factors were also found, to ensure proper dimensioning of natural-convection heat sinks. They corresponded to non-dimensional parameters that dictated the physical behaviour of the heat sink with respect to the semi-empirical correlations. Typically, the ratio of the spacing to the optimal spacing predicted by Bar-Cohen’s correlation turned out to be an appropriate grouping factor with a threshold of 1, above which the fins could safely be considered as isolated, thus greatly simplifying all further calculations.

1. Introduction

As our modern societies rely ever more on electronic devices, efficient solutions for cooling must be found. This is very challenging, as those solutions must fulfil many requirements and trade-offs, such as being cheap to manufacture, efficient from a technical point of view and safe overall. Passive cooling solutions often come forward for those applications, as they do not have any mechanical parts, thus requiring minimal maintenance. In general, they can dissipate between $0.1$ and $10\mathrm{W}/{\mathrm{cm}}^2$, but this value obviously depends on many factors, such as shape, density, materials, etc. [1]. They use natural convection, as opposed to forced convection for active cooling. This value is obviously less than for the active solutions; however, they are more energy-efficient since they do not require any source of energy to work. Consequently, they have an extended lifespan, which is why they are often used for devices that will not be frequently monitored physically, such as remote radio heads. They are silent, which is a clear benefit compared to active cooling appliances, such as pumps or fans, which can be very loud. In addition, they do not require any electrical power to work, which is very beneficial for small devices that have to work remotely for long periods of time with limited power [1].

In addition, safety is a major concern when dealing with heat-dissipation devices, and temperature is by far the major source of failure [2]: it is obviously very detrimental for the heat sink itself, for the device it aims to cool and for people in the vicinity, which is why sufficient heat-dissipation capacity should be ensured.

The development of efficient, reliable and low-cost passive cooling solutions is a major concern from social, economical and environmental points of view. However, as will be discussed later, no results are known for heat sinks for a large variety of conditions, typically for $El\gtrsim {10}^5$ or $G{r}_L\gtrsim {10}^8$. Further physical interpretation based on those parameters will be provided below.

1.1. Semi-Empirical Correlations for $\overline{N{u}_S}$

This article focuses on vertical plate-fin heat sinks, i.e., an array of rectangular fins periodically placed on a baseline that is also rectangular. The definitions of the dimensions are shown in Figure 1. S is the spacing between two adjacent fins, L is the vertical length and t is the thickness of the fin, which will be fixed here.

The Nusselt number (Nu) is a dimensionless quantity used to quantify the magnitude of the convective heat transfer relative to the conductive heat transfer within a given fluid. It compares actual heat transfer occurring by convection with respect to heat transfer by pure conduction, i.e., cases in which the fluid is stationary. Heat transfer through these natural-convection heat sinks is calculated using the average Nusselt number, defined as

$$\overline{N{u}_S}={\displaystyle \frac{hS}{k}},$$

(1)

where h is the heat transfer coefficient averaged over the whole surface of the fin—hence, the overline in the notation—k is the thermal conductivity of air at the film temperature and $\Delta T$ is the difference between the temperature at the surface of the fin and the free stream temperature, i.e., $\Delta T=T_s-T_{\infty }$. The average heat transfer coefficient can be expressed in relation with dimensionless numbers (typically, Elenbaas numbers or Rayleigh numbers (Ra)), which quantify heat transfer in natural-convection devices. Rayleigh numbers are defined as the product of Grashof numbers (Gr) and Prandtl numbers (Pr). Grashof numbers (Gr) are defined as the ratio of the buoyancy force (which promotes the fluid motion) to the viscous force (which resists the fluid motion):

$$G{r}_s={\displaystyle \frac{g\beta \Delta T{S}^3}{{\nu }^2}},$$

(2)

In Equation (2), g, $\beta$, $\nu$ represent acceleration due to the gravity constant, the volumetric thermal expansion coefficient and the kinematic viscosity of the fluid, respectively. Prandtl numbers (Pr) are defined as the ratio of the momentum diffusivity to the thermal diffusivity. Hence, Rayleigh numbers describe the natural-convective heat transfer in physical systems where fluid motion is governed by change in fluid density due to temperature difference, which occurs through the heat transfer to the fluid from a hot surface. Physically, Elenbaas numbers are a re-scaled version of Rayleigh numbers (Ra), with a specific focus on quantification of the natural-convective heat transfer in horizontally aligned parallel plates. Those correlations, which are often semi-empirical, are used as references to be compared with new data. The first was introduced by Elenbaas in 1942 [3], and it uses numbers that were named after him. In the following equation, $R{a}_S$ denominates Rayleigh numbers calculated with S as the typical length:

$$El=R{a}_S·{\displaystyle \frac{S}{L}}=R{a}_L·{\displaystyle \frac{S^4}{L^4}}.$$

(3)

The Elenbaas correlation shows that the heat transfer coefficient, and thus the heat flux, if the contact area is fixed, only depends on one factor, $El$:

$$\overline{N{u}_S}={\displaystyle \frac{1}{24}}El{\left(1-exp\left(-{\displaystyle \frac{35}{El}}\right)\right)}^{3/4}.$$

(4)

This correlation is only valid in a limited range, ${10}^{-1}\lesssim El\lesssim {10}^5$ [4], and it also assumes that fins are isothermal, which is not the case in most heat sink applications, as the tip of the fin will obviously be colder that its base. In addition, this correlation was developed and tested only for laminar flows, i.e., with low Grashof numbers. This correlation calculates the heat transfer coefficient for a plate-fin heat sink on the assumption that all the fins will dissipate the same amount of heat, which is not the case [5]. Further research has shown that the Elenbaas correlation has a rather unique form. Most other correlations have the following form [4,6,7,8]:

$$\overline{N{u}_S}={\left({\left(C_{1}E{l}^{\alpha }\right)}^{-\gamma }+{\left(C_{2}E{l}^{\beta }\right)}^{-\gamma }\right)}^{-1/\gamma }$$

(5)

where $C_1,C_2,\alpha ,\beta ,\gamma$ are positive constants and $\alpha >\beta$. This form of correlation shows that two limits can be taken into account: cases where the plates are isolated, corresponding to large Elenbaas numbers, and the opposite, cases where the developed regions of the flow for each fin are merged.

Other semi-empirical correlations have subsequently been proposed, and these were systematically compared with Elenbaas’. Bar-Cohen and Rohsenow, using the fact that at the inlet the plates could be considered as isolated and that boundary layers were overlapping at the outlet, proposed a new correlation in 1984 [8]:

$$\overline{N{u}_S}={\left({\displaystyle \frac{576}{E{l}^2}}+{\displaystyle \frac{2.87}{E{l}^{1/2}}}\right)}^{-1/2}$$

(6)

This correlation has a different form compared to Elenbaas’. This is because it assumes two asymptotic behaviours, and they are merged using a least square fitting method [8]:

• For small values of $El$, the fins can be considered under the fully developed limit, i.e., the boundary layers of two adjacent fins are deeply merged. This means that their typical boundary layer thickness $\delta$ is smaller or much smaller than $S/2$. In this case, the heat transfer coefficient h cannot be easily calculated; hence, the need for CFD simulations.
• For large Elenbaas numbers, S is large enough to say that the fins are isolated. In this case, the calculations are a lot easier, as correlations such as Churchill and Chu’s can be applied. This will be further explored below.

Those two asymptotic behaviours can be seen by either letting $El\to 0$ or $El\to \infty$ in Equation (6). The two asymptotes intersect at $El=34$.

Van de Pol and Tierney found a similar correlation a few years previously [9]. All these correlations aim at yielding an optimal spacing, $S_{opt}$, that maximizes the total heat flux from the plate Q, which is proportional to the total area of contact A. The optimal spacing can be derived from Equation (6):

$$S_{opt}=2.71{\left({\displaystyle \frac{R{a}_S}{S^{3}L}}\right)}^{-1/4}=2.71LR{a}_L^{-1/4}$$

(7)

Finally, a more recent correlation was introduced by Kim et al. to take into account the differences in the heat transfer coefficient among all the channels [6]. Indeed, it is an inaccurate assumption that all fins dissipate the same amount of heat: the central channels typically have higher heat transfer coefficients [5], due to a chimney effect:

$$\overline{N{u}_S}={\left({\left(0.09112E{l}^{0.6822}\right)}^{-3.5}+{\left(0.5170E{l}^{0.2813}\right)}^{-3.5}\right)}^{-1/3.5}$$

(8)

All these correlations use either CFD (Kim) or experimental data (all three correlations) only corresponding to laminar cases. Typically, Elenbaas himself claimed that his correlation is valid for $L\lesssim 0.5\mathrm{m}$ [3]. No literature, whether concerning simulations or experiments, is to be found on turbulent modelling of heat sinks. It is to be noted that the usual thresholds $G{r}_L={10}^9$–${10}^{10}$ for turbulent flow [4] may not apply here, as they are valid for an isolated vertical plate, but they still give a good order of magnitude.

In addition, transition boundary layers are present in heat sinks for natural-convective flows if appropriate conditions are met [10]. This makes this geometry particularly difficult to simulate:

• The Elenbaas correlation and other semi-empirical correlations being developed for laminar flows will underestimate $\overline{N{u}_S}$ if the flow is turbulent or under transition, as turbulence is known to enhance heat transfer.
• On the other hand, a fully turbulent model will overestimate $\overline{N{u}_S}$: the reality lies between those two extreme cases, and this study aimed at providing guidelines to conduct such simulations.

1.2. Novel Heat Sink Designs Using Shape Optimization Techniques

In recent years, many new designs have been proposed, aiming to either decrease thermal resistance or increase the heat transfer coefficient [11]. Sometimes, a variant of the heat transfer coefficient, relative to the mass of material used or its volume, for instance, is used [7]. Some designs consist of relatively simple changes to the classical plate-fin heat sink, as proposed by Abbas and Wang [12]. Offsetting a fin out of two could decrease $R_{th}$ by $60\%$ at most, providing S is small enough. It should also be pointed out that the standard and RNG k-$ϵ$ models were used. Likewise, Feng et al. introduced a design aiming to reduce the boundary layer length [13]. It consists of a cross-fin heat sink, where the central fins are perpendicular to the other ones. The heat transfer coefficient could be improved by up to $11\%$.

Chang et al. introduced a design where dimples were dug in the fins [14]. Rayleigh numbers up to $1\times {10}^8$ were tested for each of the four cases, where various depths of dimples were considered. Comparison with experimental values was also undertaken. In the best configuration, an increase of $65\%$ in $\overline{Nu}$ was claimed. Ray et al. illustrated a well-known trade-off phenomenon in heat sinks design, concerning the appropriate number of fins (which is, obviously, linked with $S_{opt}$) [15]. Ray et al. proposed a design with interrupted or branched fins, which, once again, aims at optimizing the boundary layers and increasing the heat transfer coefficient [15]. The best design was associated with a diminution in weight of $12.26\%$. Finally, McCay et al. studied many new designs and conducted parametric studies [11], including rectangular, trapezoidal, curved and a comparison with Zhang et al. [16]. The trapezoidal and curved cases displayed very good thermal performance, by both simultaneously reducing the thermal performance and increasing the heat transfer coefficient.

A summary of all the previous designs is available in

Table 1. Summary of novel natural-convection heat sink designs.

Article	Design	/C	/E	${\mathit{Gr}}_{\mathit{L}}$$\times {10}^6$	Main Results
[12], 2020	Offsetting fins	yes	yes	3∼40	$R_{th}\downarrow 60\%$
[13], 2018	Cross-fin	no	yes	N/A	$h\uparrow 11\%$
[15], 2022	Interruptions and branches	no	yes *	$0.1$∼10	Heat transfer ↑ $m\downarrow 12\%$
[14], 2017	Dimples	yes *	yes	15∼150	$\overline{Nu}\uparrow 65\%$
[17], 2014	(Inverted) Y-shape	no	yes *	N/A	$R_{th}\downarrow 30\%$
[11], 2023	Trapezoidal Curved W-shape	yes	no	3	Depends on case $R_{th}\downarrow 5\%$ $h\uparrow 5\%$
[18], 2015	Pin-fin	yes	yes	$0.1$∼4000	Plate-fin better for Q Pin-fin better for $Q/m$
[19], 2018	HT stagnation zone removed	yes	no	N/A	$R_{th}\downarrow 8\%$ $h\uparrow 15\%$
Current report, 2024	Plate-fin HS Several models	yes	yes *	$9.2$∼120,000	Presented after

, along with the typical values of $G{r}_L$ that were tested, when the heat sink was vertical, as $G{r}_L$ cannot be properly defined in other cases. Most of the time, those values were not provided by the authors themselves, so they were computed from any available data. This table was inspired by and adapted from previous work in a report submitted in January 2024. The column /C means that the results were compared with a semi-empirical correlation, and the asterisk means that the verification was conducted with data from another paper instead of a correlation. On the other hand, the column /E means that the article also contains experimental work, and the asterisk here means that the experimental data came from another paper. The main results are also presented, with upward and downward arrows $\uparrow ,\downarrow$ meaning an increase or a decrease compared to a base case. This base case was not the same for all the articles, so the comparison between those values should be made with extreme caution [11].

Among all the previously cited articles, most of them studied a limited range of Grashof or Rayleigh numbers. Typically, they ranged around $G{r}_L={10}^6$–${10}^9$. Consequently, the present study considered a much wider range of Grashof numbers to investigate more types of flows.

1.3. Turbulence Modelling for Natural-Convection Applications

Most of the studies in the previous section used Ansys Fluent with the laminar model. Unfortunately, very little data were found regarding simulations of heat sinks in free convection with turbulent models. In this section, some other natural-convection applications will be briefly presented, along with the turbulent models that were used. Less attention will be paid to the results, since these applications did not concern heat sinks. It is also to be noted that all these simulations involve heat transfer, as in the present study. Some of these geometries present a relative degree of similarity with a heat sink, although their main goal is not to dissipate heat.

Yet, many industrial applications may concern turbulent applications, as has been mentioned. For instance, in data centres, cooling solutions can be packed one on top of the other, so that the inlet of the top heat sink may be turbulent as a result of the flow motion in the bottom one. In telecommunications applications, very large heat sinks are also very common, as they sometimes only rely on this source of cooling, thus requiring very large heat sinks.

It is also to be noted that a comparison with LES (Large Eddy Simulation) was made by Ma and He in the case of the horizontal cylinder [20]. The turbulent models are listed for each application in

Table 2. Summary of the most-used turbulent models for natural-convection applications.

Article	Design	Turbulent Models
[24], 2017	Convection between heated and adiabatic wall	Spalart-Allmaras, k-$ϵ$, k-$\omega$ and derivatives
[25], 2020	Vertical hollow cylinder with isothermal wall	std k-$ϵ$, SST k-$\omega$
[26], 2017	Vertical cylinder with annular fins inside	std k-$ϵ$
[20], 2021	Horizontal cylinder with lateral thermal walls	SST k-$\omega$, k-$ϵ$, transition-SST
[27], 2012	Convection between walls at different temperatures	k-$ϵ$
Current report, 2024	Plate-fin heat sink	transition-SST, SST k-$\omega$, k-$ϵ$ and derivatives

. The transition-SST model, although it provides very good results, is seldom chosen in CFD studies [20], in particular when it comes to examples of heat sinks. There are, however, a few examples of forced convection applications where it was used and yielded very good results: for instance, for dual jets [21] or impinging jet flows [22]. The k-$ϵ$ model and its variants are by far the most-used, whereas they are known to be rather inaccurate in the near-wall region [23]: that is because it often comes as the default two-equation model for turbulence.

In summary, new correlations of the form $\overline{N{u}_S}=f\left(El\right)$ have been found, and those correlations will be compared to simulation results [3,6,8]. Nevertheless, they are only valid for a limited range of Elenbaas numbers, and they have never been tested for turbulent flows. Very little is known concerning turbulent flows around heat sinks. Still, several similar applications in free convection have been simulated using turbulent models. This study therefore aimed at comparing different flow models for natural-convection heat sinks, for several geometric models.

As a consequence, the present study aimed at investigating the behaviour of vertical natural-convection heat sinks over a wide variety of operating conditions. The ranges of parameters used are listed in

Table 3. Ranges of values for geometrical and physical parameters.

Description	Symbol	Value	Unit
Vertical length	L	$0.1$–$1.5$	$\mathrm{m}$
Half-spacing (width)	$S/2$	2–50	$\mathrm{mm}$
Temperature difference ($T_s-T_{\infty }$)	$\Delta T$	60–300	$\mathrm{K}$
Elenbaas number	$El$	$1.17$–$3.42\times {10}^7$	$-$
Grashof number	$G{r}_L$	$9.22\times {10}^6$–$1.55\times {10}^{11}$	$-$

: such a wide variety of Grashof and Elenbaas numbers had never been reached before. These are the typical ranges of operation that are commonly observed in natural-convection heat sink devices used for electronic cooling applications [28]. Two main geometrical models were developed, namely, 2D and quasi-3D, as turbulence is known to lead to high deviations in 2D [4,23,29]. Two main transition models were tested: standard k-$ϵ$, as it is by far the most commonly used, so data has to be generated for it, and transition-SST, as our case involved flow under a transition regime. Where possible, our results were compared to semi-empirical correlations, to assess their validity in different regimes.

2. Materials and Methods

2.1. Setup

All our simulations were carried out using the computational fluid dynamics (CFD) software Ansys Fluent 2023R1. Fin thickness was regarded as a constant, and it was taken as $t=2\mathrm{mm}$. The free stream temperature was fixed at $T_{\infty }=300\mathrm{K}$ for all the simulations. The parameters in

Table 4. Material properties.

Description	Symbol	Value	Unit
Dynamic viscosity	$\mu$	$1.79\times {10}^{-5}$	$\mathrm{kg}/\mathrm{m}/\mathrm{s}$
Density (Boussinesq approximation)	${\rho }_0$	$1.225$	$\mathrm{kg}/{\mathrm{m}}^3$
Thermal conductivity	k	$0.0242$	$\mathrm{W}/\mathrm{m}/\mathrm{K}$
Specific heat capacity	$C_p$	$1006.43$	$\mathrm{J}/\mathrm{kg}/\mathrm{K}$
Thermal diffusivity	a	$1.963\times {10}^{-5}$	${\mathrm{m}}^2/\mathrm{s}$
Prandtl number	$Pr$	0.744	$-$

were also considered for the working fluid.

A pressure-based solver was used, and flow and energy equations were solved. The Boussinesq approximation was enabled for the fluid.

For the laminar and k-$ϵ$ simulations, the flow was initialized with a velocity of $1\mathrm{m}/\mathrm{s}$ in the +y direction, all other parameters being set at zero. Without a proper initialization of the flow fluid, the transition-SST model is very long in converging, even with under-relaxation factors close to 1, so the flow fluid was first initialized with a converged corresponding laminar case.

The following convergence criteria were enforced. Residuals were monitored, for continuity, velocities and energy, and also for turbulent parameters if such a model was used. But the most useful convergence criterion was to monitor the area-weighted average of the heat flux, and enforce that its residuals on the last ten iterations were smaller than $1\times {10}^{-5}$.

A mesh convergence study was conducted, using a total number of elements ranging from 340 k to 1.25 M. This led to a Richardson extrapolation grid convergence ratio of $r_g=0.63$ (being the ratio of the difference between the medium and fine grid solutions and the coarse and medium ones), thus comprising between 0 and 1, which meant that monotonic convergence was reached. A grid convergence index of $GCI=1.1\mathrm{W}/{\mathrm{m}}^2$ or $GCI=0.25\%$, as defined by Roache [30], which is small enough to consider the mesh convergence study validated.

A domain size sensitivity study was also carried out. It aimed to choose the size of the fluid domain encapsulating the heat sink. It had to be as small as possible to limit computational cost, yet large enough that the results would be independent of the size. For such applications, it is recommended to use a length of $3L$ above and $2L$ below [23]. After careful consideration for our case, this was also deemed sufficient. In particular, the length below had a negligible influence as long as it was not too small.

The pressure-velocity coupling was chosen as SIMPLE after various tests, although it had a rather small influence. The boundary conditions’ turbulent parameters were also chosen accordingly to match theoretical predictions, namely, the turbulent Prandtl number $P{r}_t$, the turbulent intensity and the length scale.

2.2. Definition of Geometrical Models and Boundary Conditions

Three models were defined, of increasing complexity, but mostly two were used because of computational limitations. They are all represented in Figure 1. It is to be noted that the fluid domain encapsulating this heat sink is not drawn here for clarity:

• The full model is in grey. It simulates the whole heat sink—thus, the conduction in the fins—and it is the most accurate, as it does not consider that all the fins dissipate an identical amount of heat. Nevertheless, it is computationally very expensive even with small domain sizes, and thus was not adapted for our study, as large values for L could not be properly simulated in reasonable time. For that model, inflation layers were used, whereas the edge sizing method was preferred for the others, as it meshes the model faster, and it leads to a structured mesh.
• An intermediate solution was proposed in blue, still in Figure 1. It will be referred to, hereafter, as quasi-3D, although strictly speaking it was a 3D model from the point of view of the simulation. Periodic (or symmetry) boundary conditions were applied on sides normal to the z direction. As a consequence, conduction in the fins was neglected, so the fin efficiency was $100\%$. It is important to note that the Elenbaas and Bar-Cohen correlations were developed by calculating the heat flux between the two isothermal plates. Hence, the fin efficiency would always be $100\%$. The quasi-3D models were developed to mimic this scenario, so that the fin efficiency would be $100\%$. The interest resided in the fact that three-dimensional models are known to predict much better turbulent phenomena [4,23,29]. The model encapsulated half a fin and half a channel, thus resulting in a width of $S/2$, as the channel had a total width of $S-t$.
• The simplest model is represented in red. It simply consisted of a 2D slice, and the computational complexity was decreased a lot compared to the last two models. It can nevertheless be assumed that its performance on turbulent simulations will be rather poor.

Very few instances in the literature have compared these different models for natural-convection applications. One goal of this report was to see what level of accuracy simpler models can give.

The boundary conditions for the 2D model were chosen as in Figure 2. The wall boundary conditions consisted of a no-slip isothermal wall. Imposed temperature was preferred over fixed heat flux, as this might lead to divergence, especially for the 2D case, given that the slice had no thickness. As already mentioned, in the quasi-3D case the planes normal to the z direction were defined as periodic boundary conditions. Sometimes, for technical reasons that arose from the creation of periodic boundary conditions by Fluent, they were replaced by symmetry boundary conditions. Multiple tests were performed, and this replacement had a negligible influence. For example, with $L=1\mathrm{m}$, $S/2=6\mathrm{mm},\Delta T=60\mathrm{K}$, the difference between the symmetry and the periodic boundary conditions in the z direction was $0.8\%$. For the pressure inlet and outlet boundaries, a gauge pressure of $0\mathrm{Pa}$ was chosen.

It is to be noted that a simple models comparison could be conducted with a 3D model. Obviously, this had a better correlation with semi-empirical formulas, typically less than $\pm 3\%$ with Elenbaas’.

3. Results and Discussion

3.1. Parametric Study in 2D

In the first place, a parametric study was carried out in 2D with the k-$ϵ$ and laminar models. It used more than 60 design points, which were different combinations of the following parameters:

• $L=0.1\mathrm{m},0.5\mathrm{m},1\mathrm{m},1.5\mathrm{m}$;
• $S/2=2\mathrm{mm},6\mathrm{mm},10\mathrm{mm},15\mathrm{mm},20\mathrm{mm},50\mathrm{mm}$;
• $\Delta T=60\mathrm{K},100\mathrm{K},300\mathrm{K}$.

This allowed for a large range of Elenbaas and Grashof numbers—more specifically, $El\in [1,3.42\times {10}^7]$ and $G{r}_L\in [9.20\times {10}^6,1.55\times {10}^{11}]$. Many points were also added in the intermediate region $10\lesssim El\lesssim 100$, which also coincided with $S/S_{opt}\approx 1$. The heat transfer coefficient was calculated by dividing the total heat transfer rate (from the fin surface to the fluid) by the product of half the fin area and the temperature difference ($\Delta T$). The detailed calculation procedures are available in our recent publication (refer to Equation (2) in McCay et al. [11] for details).

This ensured that a wide variety of flows was covered, from fully laminar to mostly turbulent. The results for the k-$ϵ$ model are presented below, but many interesting conclusions can already be drawn from the results for the Fluent laminar model. It is to be noted it has been applied whatever the $G{r}_L$ value was, so very likely some turbulent situations were modelled with the laminar model. This is illustrated in Figure 3. Particular attention was paid to having many values for $G{r}_L$ for all the ranges of Elenbaas numbers. This was particularly relevant, as it allowed for uncorrelating those variables.

The Elenbaas semi-empirical correlation was valid up to $El\approx {10}^7$ at least, whereas the current upper limit that had previously been tested was ${10}^5$ [4].

However, at a high Grashof number, the laminar model only dissipated heat through viscous effects, since the Kolmogorov cascade was not considered. Turbulence, because of eddies, led to more vigorous fluid mixing in the boundary layers, thus increasing the heat transfer. This phenomenon was not captured by the laminar model. The boundary layer thickness and profile also largely depended on the model, and the laminar model could not predict a turbulent boundary layer.

Interestingly, the main region of difference between the current data points and the correlation was in the range where it was supposed to be the most precise, i.e., the intermediate region $10\lesssim El\lesssim 100$. The number of data points in that region seemed to corroborate the fact that the semi-empirical correlations might be inaccurate in that domain, and the deviation was all the more important for $G{r}_L\gtrsim {10}^9$. This was because the correlation was not used in its valid range for most points in that range. In other words, the laminar model was forcing the modelling of turbulent flows as laminar. This proved that the semi-empirical correlations were not valid for turbulent flows, as stated earlier.

Also, this discrepancy could be the result of inaccurate predictions regarding asymptotic behaviours. Typically, the Bar-Cohen correlation assumes two asymptotic behaviours, as mentioned earlier, but, unfortunately, the region where the relation is supposed to be the most accurate, typically $10\lesssim El\lesssim 100$, is also where the two asymptotic behaviours intersect, thus leading to high deviations, as the merging process, corresponding to the exponent $-1/2$, has no physical meaning but is only a statistical modelling.

3.2. New Grouping Factors: Examining Physical Parameters Governing Heat Sink Performance

One of the major goals of this work was also to find grouping factors. Grouping factors are physical quantities that dictate the overall behaviour of the system, with respect to semi-empirical correlations. If possible, those quantities should be dimensionless, so that they can be easily transposed to any geometry. In that sense, they could also be named scaling factors. However, further research would be needed to prove that they are indeed independent of the size and shape of the heat sink. The Elenbaas number is a clear example of a grouping factor for this problem, as all the correlations from Section 1.1 are of the form $\overline{N{u}_S}=f\left(El\right)$.

The first strategy to find other grouping factors was to group data points by value of external parameters that were not directly related to the Elenbaas number itself, as depicted in Figure 3: $G{r}_L$ was not a grouping factor. Conversely, $S/S_{opt}$ stood out as an appropriate one. As a reminder, $S_{opt}$ was here defined as the optimal grouping factor according to Bar-Cohen’s correlation, as in Equation (7). $S/S_{opt}$ was, consequently, an intricate function of the physical and geometrical parameters of the problem, i.e., $L,S$ and $\Delta T$, as $R{a}_L$ itself depended on L and $\Delta T$. Figure 4 clearly shows that for $S/S_{opt}$ slightly under 1, the current data points would switch between both asymptotic behaviours. $S/S_{opt}$ stands out as the best way to normalize S, and the following physical interpretations can be given, which are backed up by the numerical simulations:

• $S/S_{opt}\lesssim 1$: the boundary layers of two adjacent fins are deeply merging, thus corresponding to low $El$ behaviour,
• $S/S_{opt}\gtrsim 1$: the fins can, to some extent, be considered as isolated—that is, high $El$ behaviour.

3.3. Comparison with Churchill and Chu’s Correlation

Another means of validation of the data would be to compare it with Churchill and Chu’s correlation. This correlation is valid for natural convection for vertical plates, so it does not properly apply to an array of fins. However, it is valid for a much wider range of flows compared to the previous correlations, namely, ${10}^4\lesssim R{a}_L\lesssim {10}^{13}$, and, in particular, can account for turbulent flows. Obviously, S is not a parameter for this correlation, which is why it will be plotted in a different manner than for Figure 3 and Figure 4, and it reads [4,31].

$$\overline{N{u}_L}={\left(0.825+{\displaystyle \frac{0.387R{a}_L^{1/6}}{{[1+{(0.492/Pr)}^{9/16}]}^{8/27}}}\right)}^2$$

(9)

As a consequence, the data points from the 2D Fluent laminar model can be plotted in a graph $(R{a}_L,\overline{N{u}_L})$, and they can be compared with Churchill and Chu’s relation, as in Figure 5a. This graph is interesting, in the sense that it validates our model for large values of S, typically larger than $S_{opt}$: the values have a small deviation from the correlation. It also means that Churchill and Chu’s correlation is applicable for large spacings to find the total heat transfer coefficient from $\overline{N{u}_L}$, and that although the other correlations are valid in those cases, they are not necessary [31]. Finally, it also confirms that $S/S_{opt}$ is a grouping factor for natural-convection heat sink applications.

In order to precisely find the threshold for $S/S_{opt}$, a more detailed study was performed with values of $S/S_{opt}$ around 1, as mentioned above. Most values for $S/S_{opt}$ ranged between $0.8$ and $1.2$; hence, the rescaling for the values of $S/S_{opt}$ represented in Figure 5b. This shows that the threshold between both behaviours occurs at exactly $S/S_{opt}=1$, and also that this transition is less sharp for low Rayleigh numbers.

As mentioned above, the 2D model is not really adapted for turbulence modelling. Indeed, turbulence is inherently a 3D model [4,23], which is why the transition-SST led to very poor performance on the 2D model. The k-$ϵ$ models yielded somewhat better results, but were still not acceptable when compared to the correlations. This is why a similar parametric study was conducted with the quasi-3D model.

3.4. Parametric Study in Quasi-3D

For this section, the simulations were carried out using the quasi-3D model, also referred to as the periodic model. Once again, special attention was paid to having all kinds of flows, i.e., a wide range of $G{r}_L$ in every Elenbaas number regions.

A comparison was made of the Fluent laminar, transition-SST and standard k-$ϵ$ models. The results are presented in Figure 6. Once again, good correlation was observed between the theoretical correlations and the current CFD data, which validated the models.

The k-$ϵ$ was retained for study, as it is the most widely used turbulent model in the literature, as shown in Section 1.3. In addition, under certain conditions, it is possible to optimize the meshing strategy to keep $y^{+}>30$, thus significantly increasing computational time, which is impossible with SST models.

The transition-SST model had fewer points than the other two. This was because it had to be initialized from a converged laminar case, thus making direct parametric studies harder to conduct and extremely long to converge without proper initialization.

Turbulence is known to enhance heat transfer, as it can be theoretically predicted [4], or by numerical simulations [11,29]. Therefore, it seemed logical that the value predicted by the transition-SST model would lie between the laminar and the other turbulence models, as it is known to accurately take into account the transition zone between both regimes. Furthermore, the k-$ϵ$ model is known to have relatively poor performance in the near-wall region, i.e., in the boundary layers. This was why its value was the most extreme and the furthest from the predictions of the semi-empirical correlations.

However, this order is not true for high $El$, typically greater than ${10}^4$. As a consequence, the transition-SST model should not be used in that range. It can also be confirmed for a wider range of Elenbaas numbers that the k-$ϵ$ model overestimates the heat flux with respect to the laminar model: this is a well-known consequence of turbulence models. Overall, all those conclusions are consistent with what was observed so far, and also in the literature [11,29]. In summary:

• The laminar model provided good estimations of the heat flux, whatever the operating conditions: for instance, it was not dependant on the value of $G{r}_L$; 2D simulations were sufficient for the laminar model, but it was the only model for which this was true.
• The transition-SST model was deemed to be the most accurate for $El\lesssim {10}^4$, which indeed seems consistent, as it was especially developed for this kind of application.
• The k-$ϵ$ model tended to overestimate the heat flux over the whole range of Elenbaas numbers. Consequently, it should be used with precaution, although it gives a good estimation of the heat flux in a shorter delay as the mesh can be coarser, thanks to the usage of wall functions, notably Enhanced Wall Treatment (EWT). More precisely, the k-$ϵ$ model can, under certain conditions, be used with a meshing strategy allowing $y^{+}>30$ for the first cell, unlike many other models, typically the SST models [23]. For this reason, calculations could be much faster, as the number of elements in the mesh went down from ${10}^5$–${10}^6$ to ${10}^4$–${10}^5$ in order of magnitude. Also, the deviation for average values of the y velocity from the laminar model was the highest, which might mean poorer prediction regarding the flow equations.

The same methodology as in Section 3.2 was applied with the quasi-3D results, with similar results. $S/S_{opt}$ was also found to be an appropriate factor with a very clear threshold at 1, whatever the Fluent model. This acted again as a validation of the model, and confirmed once again $S/S_{opt}$ as an important parameter regarding the simulation strategy.

The turbulent models themselves should be accurately validated, to fit them to general simulations in natural convection, without even considering the aforementioned correlations. For instance, LES (Large Eddy Simulation) models could be used for this purpose.

3.5. Influence of the Number of Cells in the z Direction

The k-$ϵ$ model provided significantly different results when comparing the 2D with the quasi-3D models. As a consequence, a study was conducted on the number of cells in the z direction—only for this model, as the results were more blatant. As a reminder, the z direction is the direction normal to the periodic boundary condition planes. A higher number of cells in the z direction should mean better prediction, but also increases the number of cells in a linear way as the mesh is fully structured.

Three design points were chosen at low, medium and high Elenbaas numbers: cases C1, C2 and C3, and the number of cells varied from 1 to 14.

C1.

$L=0.1\mathrm{m}$, $S/2=2\mathrm{mm}$, $\Delta T=60\mathrm{K}$, $El=17$, $G{r}_L=9.20\times {10}^6$

C2.

$L=1.5\mathrm{m}$, $S/2=6\mathrm{mm}$, $\Delta T=60\mathrm{K}$, $El=95$, $G{r}_L=3.10\times {10}^{10}$

C3.

$L=1.5\mathrm{m}$, $S/2=50\mathrm{mm}$, $\Delta T=300\mathrm{K}$, $El=2.28\times {10}^6$, $G{r}_L=1.55\times {10}^{11}$

If that number equalled zero then it meant that a 2D result was considered. As usual, the average values of $v_y$ and $q^{″}$ are presented on two separate scales in Figure 7. For a small number of cells, typically less than four, there were significant variations in the results provided by the models. The predictions deviated more from the correlations for low Elenbaas numbers (C1), as shown above: thus, it is recommended to have more than 10 cells in the z direction for those cases. This also explains why the 2D model is not well adapted for k-$ϵ$ simulations and should be linked with the previous subsection.

There was also a higher discrepancy between the 2D and 3D results for the high Grashof numbers, and this has also been verified in the literature [29] and with other design points in this study. Typically, the difference between 2D and 3D results regarding the average velocity $v_y$ is given in each caption (Figure 7). This error increased with $G{r}_L$, corroborating the fact that the 2D model cannot capture highly turbulent effects. For such applications, the quasi-3D model is recommended.

4. Conclusions

This work focused on providing an efficient and accurate strategy to simulate buoyancy-driven flows around heat sinks. Plate-fin heat sinks were mostly studied, as they are the only ones for which several semi-empirical correlations were developed.

The 2D model overestimated the heat flux by approximately 10 to $20\%$ compared to the Elenbaas correlation, but was obviously less computationally demanding. It is recommended to use only the laminar model on the 2D model, as any turbulent model will result in high deviation from the correlation. The laminar model also proved wrong for high Grashof numbers, as could intuitively be predicted: turbulent models must absolutely be used in such cases. It was also found that the Elenbaas correlation is in reality valid up to $El\approx {10}^7$ at least. Elenbaas numbers up to $3.42\times {10}^7$ and Grashof numbers up to $1.55\times {10}^{11}$ were simulated.

The quasi-3D model provided the best results, and proved to be a very good trade-off between quality and computational time. With the Fluent laminar model, the results were very close to the 2D case, but this model had a smaller deviation for the transition-SST and k-$ϵ$ models. It was found that the k-$ϵ$ model provided relatively good results with a coarser mesh for any Elenbaas number. For low and medium Elenbaas numbers, the transition-SST model provided the most consistent results, but it should be avoided for high Elenbaas numbers.

Further research could be carried out to validate the models against more results. This could be by experimental means, but it seems that an LES simulation on natural-convection heat sinks, even small ones, has never been conducted. Other turbulent models could also be taken into consideration. In addition, turbulent models were only applied to vertical heat sinks, as the theoretical correlations are mostly valid for those, but other geometries, such as those mentioned in Section 1.2, could also easily be simulated. The effect of turbulent inlet conditions also merits further investigation, since in practical application natural-convection heat sinks are often used in harsh environments or near other components (e.g., in exposed outdoor locations, such as remote radio heads, or in high-vibration automotive or aviation applications).