Thermal Performance Optimization of Perforated Fins for Flat Plate Heat Sinks Using CFD Approach

The rapid evolution of modern technology has placed increasing demands on researchers to develop innovative heat dissipation systems that can effectively manage the significant heat output of high-tech electronic systems, such as high-speed microprocessors [1]. These systems necessitate highly efficient thermal management solutions to maintain optimal operation and prevent severe thermal damage.

Conventional heat sinks, equipped with plate fins, have been extensively utilized in the electrical and electronic applications industry [2]. Their wide use is primarily attributed to their ease of manufacture, cost-effectiveness, and their ability to enhance the heat transfer coefficient (h) between the system and its surroundings. However, such conventional heat removal methodologies have proven to be inadequate in handling the heat dissipation demands of advanced high-tech systems [3, 4]. Their wide use is primarily attributed to their ease of manufacture, cost-effectiveness, and their ability to enhance the heat transfer coefficient (h) between the system and its surroundings. However, such conventional heat removal methodologies have proven to be inadequate in handling the heat dissipation demands of advanced high-tech systems [5, 6]. This has prompted various efforts to improve the mechanical design of plate fins using different techniques to optimize their thermal performance [7-10]. A notable example of these adopted strategies is the incorporation of perforations into solid fins. The impact of perforations on the thermal performance of flat-plate heat sinks has been the subject of extensive numerical and experimental investigations [3, 11].

Several studies have been reviewed by Ahmed et al. [3] and Kumar et al. [11], focusing on strategies to enhance the thermal performance of utilized heat sinks. These studies primarily examined the influence of fin shapes and orientations, perforation geometry, and the use of microchannels on the heat dissipation rate of the cooling systems. Notably, most of the studies reviewed by Ahmed et al. [3] and Kumar et al. [11], and a majority of related research, have not explored the combined effect of varying perforation shapes and positions on the thermal performance of a heat sink equipped with flat-plate fins. Ibrahim et al. [12] observed significant increases in the heat transfer from the heat source to the fin tip when different shapes of perforations were applied to the tested solid fin. Their findings revealed that the greatest ΔT, 51.29%, was achieved with circular perforations, followed by 45.57% and 42.28% for rectangular and triangular perforations, respectively, in comparison to solid fins. Moreover, the impact of perforation surface area on h and temperature distribution across the tested fins was experimentically investigated by Awasarmol and Pise [13] and Al-Doori [14]. Their findings suggest a decrease in ΔT and an increase in h as the perforation surface area increased. The obtained results were largely attributed to the enhanced heat conduction within the perforated fins and the turbulent airflow near the fins' surfaces.

In this context, computational fluid dynamics (CFD) modelling with three-dimensional (3D) geometry has been increasingly employed as an effective tool for visualizing aerodynamic flow and temperature distribution within and near solid and perforated fins [15, 16]. Feng et al. [17] introduced a novel heat sink (HS) equipped with a crossed flat-plate fin array, tested under conditions of natural convection and radiation heat transfer. A numerically simulated model, validated experimentally, was utilized to investigate the system's thermal performance. The authors reported superior thermal performance of the proposed HS compared with conventional parallel fin HS. Further, a CFD simulation model using the ANSYS Fluent commercial code was carried out to evaluate the thermal performance (average Nu) of an HS equipped with wavy fins, comparing the results with a flat-plate HS [18]. The numerical outcomes were primarily analyzed and discussed in terms of wave number and amplitude of the fins. The authors concluded that heat transfer from the fins' surface to the turbulent air flow increased with the use of wavy fins, and thermal resistance was reduced as wave number and amplitude increased. Hussain et al. [19] numerically investigated the impact of curvature profile at the junction of the HS and fin-base, as well as airflow direction, on the thermal resistance of a flat-plate HS using the ANSYS Fluent code. The proposed CFD model exhibited superior thermal performance compared to the conventional HS.

Employing the RNG-based k-epsilon turbulent flow model and Navier-Stokes equations, major airflow parameters of the 3D fluid flow and convective heat transfer from solid and perforated fin arrays mounted on a flat-plate heat source were predicted [20]. The numerical investigation indicated superior heat transfer when linear holes were used compared to solid fins. This reduction in weight also improved thermal performance, potentially beneficial for mobile micro and high-tech applications. A CFD model was implemented on a pin fin geometry matrix to estimate optimal thermal design in terms of perforation forms and diameters [21]. The numerical investigation revealed that manipulation of dimensions, such as perforation height and diameter, are required for optimal design. The study concluded that the Nusselt number (Nu) increased and average temperature decreased compared to the original design. Al-Damook et al. [22] numerically analyzed the thermal performance of a pin fin HS with multiple circular perforations. An increase in Nu value by approximately 11% was observed when the number of perforations was increased compared to the solid fin. This also enhanced the temperature distribution over the fins body, leading to increased heat transfer between the surface and the surrounding.

There exists an algorithm proposed by Balachandar et al. [23] that targets optimized design for superior thermal performance. This algorithm is accompanied by a validated computational fluid dynamics (CFD) model, which is employed to numerically analyze the impact of elliptical perforations on the heat dissipation rate of solid fins within the utilized heat sink (HS) relative to the surrounding environment. The result demonstrated improved performance of perforated fins compared to solid fins. Furthermore, the post-processing results suggested that the dimensions and positions of the ellipse's perforations could potentially affect the heat transfer rate from the source to the surrounding environment. Al-Essa and Al-Hussien [24] investigated the effects of the orientation and geometry of two square perforations on the thermal performance of horizontal fins. To numerically resolve the presented problem, the variational approach finite element technique was adopted. The authors concluded that the thermal resistance of the solid fin decreased with the utilization of square perforations. Moreover, the heat dissipation rate was enhanced when the perforations were oriented horizontally, compared to the inclined fin.

The majority of the reviewed works, which encompass both numerical and experimental investigations, focus on optimizing the thermal performance of different types of HSs. Numerous numerical investigations have highlighted the use of multidimensional CFD modeling tools. However, few numerical studies have been conducted that investigate the combined effect of the perforations' positions and shapes on the heat transfer by conduction and natural convection on a flat plate HS. Specifically, the numerical simulation of conduction heat transfer from the heat source to the HS fins (both solid and perforated) and the natural convection heat transfer from the fins to the ambient environment, present challenges in terms of design and numerical calculation. This work proposes a developed CFD model for natural convection and conduction heat transfer to numerically investigate the combined influence of the perforations’ shapes and their positions on the thermal performance of the flat plate HS. Additionally, the effect of the generated turbulent flow in the vicinity of the fins surfaces on the thermal performance of the proposed HS is examined, and the resulting data are discussed and analyzed. Validation of the proposed CFD model is achieved through additional experimental work.

3.1 Geometry

The proposed flat-plate heat sink geometry was designed by using a SolidWorks software program and numerically examined using ANSYS Fluent code version 2019 R1. The dimensions of the proposed model were set based on the tested geometry where the required validation for the proposed numerical CFD model was performed. Three perforation zones (positions) were chosen, namely as a top, middle, and bottom, as shown in Figure 2. Besides, three perforations shapes were selected, namely circle, square, and triangle, to be tested at each perforation position in a testing matrix. The perforations were arranged in a regular matrix made of two rows and multi-columns. In order to make comparable testing conditions, the cut surface area of 78.5 mm² was set to be the same for all the used perforations.

2a.png	2b.png	2c.png
2a1.png	2a2.png	2a3.png
(a) Top circular perforation	(d) Top square perforation	(g) Top triangular perforation
1.png	2.png	3.png
4.png	5.png	6.png
(b) Mid. circular perforation	(e) Mid. square perforation	(h) Mid. triangular perforation
7.png	8.png	9.png
10.png	11.png	12.png
(c) Bot. circular perforation	(f) Bot. square perforation	(i) Bot. triangular perforation

Figure 2. Computational meshes and geometries (a-i) for the used perforated HSs

3.2 Computational mesh

Figure 2 shows the adopted geometries and computational meshes of the perforated flat-plate HS with circular perforated fins. The HS enclosure was also included, as shown in Figure 3. A pre-processing meshing tool in the ANSYS Fluent was used to generate the computational mesh for the proposed CFD model comprising the HS and enclosure. Four mesh resolutions were tested to tackle the best mesh accuracy, starting from element numbers of 350e3 to 500e3 with an increment of 50e3 elements, as shown in Table 1(a). The mesh size of 400e3 elements was chosen because it gives the best possible simulation results and minimum calculation run time. As shown in Figure 3, the mesh consists of tetrahedral grids. However, at the sharp edges, where the triangular perforation was used, hexahedral grids were adopted for mesh stability. In general, the element size was set to 3.0 mm; meanwhile, the cell size at the sharp edges was set to be 2.0 mm. Consequently, the total elements’ number comprised in the selected mesh was 400e3 elements. Table 1(a) shows the mean value of ΔT between the HS base and enclosure. Moreover, the Nusselt number (Nu) values of the convective heat transfer were used to ensure the best elements’ number for the proposed CFD model in terms of calculation's time cost and computer's hardware complexity. It can be noticed from Table 1(a) that the difference in temperature between the HS surfaces and the enclosure (T_fin-T_∞) and Nu values are stable after the element number is increased from 350e3 to 400e3 and go on reaching 500e3. Accordingly, the number of elements of 400e3 was chosen for conducting the CFD simulation for this study. Table 1(b) shows that a Tetrahedral element type with a growth rate of 1.2 are used. The maximum and minimum element sizes are 0.001 m and 0.00085 m respectively.

3a.jpg	3b.jpg
(a) Enclosure HS mesh domain and resolution	(b) Perforated fin’s mesh domain and resolution

Figure 3. Enclosure computational mesh and a zoomed in mesh sample for a circular-perforated HSs

Table 1. Mesh division statements (a) and mesh details (b)

3.3 Boundary and initial conditions

The experimental testing conditions of this study were used to set up the boundary and initialization conditions of the proposed CFD model. An enclosure of air covered the model to create a contact region between the heat sink surfaces and air inside the enclosure, which represent the interface region to enrol the convection heat transfer and buoyancy effect phenomenon. At the beginning of the calculation, the enclosure air temperature and pressure were set at 298 K and 101.325 kPa, respectively, as listed in Table 2. Besides, the major air properties were selected and set up at the standard ambient air conditions. The metal properties of the aluminium HS were set, as shown in Table 3. Equal initial surface temperatures were set for all the HS surfaces to be at 300 K. The initial calculation pressure and temperature were set to be similar to the enclosure initial conditions (T_o=298 K and P_o=101.325 kPa) while the airflow was set at a static condition (0.0 m/sec) inside the enclosure and near the fins' surfaces, as shown in Table 2. This is because selecting these options to make the numerical calculations more stable, and the results are closer to the real-time outcomes [25]. The bottom side of the model (heat sink base) was set up to be a fixed wall, and five supplied heat fluxes were used. The supplied heat fluxes of 1.0, 1.5, 2.0, and 2.5 kW/m² were adopted for the experimental validation and numerical investigation. The numerical analysis was considered by depending on the SIMPLE algorithm method to couple the velocity and pressure. In addition, the second-order upwind scheme method was used to increase the accuracy of the discretizing of the convective terms in the governed equations. The convergence limits were set at 10e-4 for the momentum equations and 10e-6 for energy equations. Besides, the convergence criteria for residuals was set at 10e-8. It worth to be mentioned that the Boussinesq approximation was used to simulate the buoyancy effect of the ambient air caused by the natural convection heat transfer from the HS surfaces to enclosure air.

Table 2. CFD simulation initial and boundary conditions

Table 3. Properties of air and aluminum

3.4 Numerical modelling: Governing equations

The investigated CFD model involved three governing stages of heat transfer from the heat source to the surrounding. Firstly, direct supplies of heat flux were applied at the bottom surface of the proposed HS base. Five values of the subjected heat flux were chosen, as shown in Table 2. Secondly, the supplied heat would be transferred by conduction inside the HS body, starting from the HS base toward the fins. Finally, the heat would be transferred by convection from the HS surfaces (base and fins) to the surrounded atmosphere (HS enclosure). Besides, the heat transfer by convection from the HS surfaces would be carried out by the working fluid (air) via the buoyancy effect. The turbulent flow nearby the surfaces and at the vicinity of perforations could also play an important role in improving the convection heat transfer. Therefore, in order to mathematically govern the three main stages of the aforementioned heat transfer methods, the governing equations from Eq. (1) to Eq. (12) were used. These governing equations can be described and summarized as follows:

1. Continuity equation can be expressed, as shown in Eq. (1).

$\frac{\partial u}{\partial x} i+\frac{\partial v}{\partial y} j+\frac{\partial w}{\partial z} k=0$                 (1)

2. The Momentum equations for the three-dimension coordinates are defined in Eqns. (2)-(4).

$-\frac{\partial p}{\partial x}+\mu\left(\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}+\frac{\partial^2 u}{\partial z^2}\right)=\rho\left(u \frac{\partial u}{\partial x}+v \frac{\partial u}{\partial y}+w \frac{\partial u}{\partial z}\right)$                     (2)

$-\frac{\partial p}{\partial y}+\mu\left(\frac{\partial^2 v}{\partial x^2}+\frac{\partial^2 v}{\partial y^2}+\frac{\partial^2 v}{\partial z^2}\right)=\rho\left(u \frac{\partial v}{\partial x}+v \frac{\partial v}{\partial y}+w \frac{\partial v}{\partial z}\right)$                       (3)

$-\frac{\partial p}{\partial z}+\mu\left(\frac{\partial^2 w}{\partial x^2}+\frac{\partial^2 w}{\partial y^2}+\frac{\partial^2 w}{\partial z^2}\right)=\rho\left(u \frac{\partial w}{\partial x}+v \frac{\partial w}{\partial y}+w \frac{\partial w}{\partial z}\right)$                       (4)

3. Energy equation for the fluid phase flow can be represented by Eq. (5).

$u \frac{\partial T}{\partial x}+v \frac{\partial T}{\partial y}+w \frac{\partial T}{\partial z}=\frac{\mathrm{k}}{\rho c}\left(\frac{\partial^2 T}{\partial x^2}+\frac{\partial^2 T}{\partial y^2}+\frac{\partial^2 T}{\partial z^2}\right)$                  (5)

4. Energy equation for the solid phase state can be represented by Eq. (6).

$0=\left(\frac{\partial^2 T}{\partial x^2}+\frac{\partial^2 T}{\partial y^2}+\frac{\partial^2 T}{\partial z^2}\right)$                          (6)

The entire flow of the air through the perforations and near the wall of the HS fins was assumed as a turbulent flow. Therefore, Re-Normalisation Group k-epsilon (RNG k-ε) turbulence model was used to numerically investigate and analyse the turbulence flow inside the enclosure [25]. It is worth to be mentioned that the Boussinesq approximation was used for the proposed model to simulate the buoyancy effect of the moving air caused by the density variation throughout the convection heat transfer. The governing equations of the turbulent flow can be defined in Eqns. (7)-(10), and for Boussinesq approximation is defined in Eqns. (11) and (12).

A- RNG k-ε mathematical model.

$\frac{\partial(\rho k)}{\partial t}+\operatorname{div}(\rho k U)=\operatorname{div}\left[a_k \mu_{e f f} \operatorname{grad} k\right]+\tau_{i j} \cdot S_{i j}-\rho \varepsilon$                    (7)

$\frac{\partial(\rho k)}{\partial t}+\operatorname{div}(\rho \varepsilon U) =\operatorname{div}\left[a_{\varepsilon} \mu_{e f f} \operatorname{grad} \varepsilon\right] +C_{1 \varepsilon} \frac{\varepsilon}{k} \tau_{i j} \cdot S_{i j}-C_{2 \varepsilon} \rho \frac{\varepsilon^2}{k}$                    (8)

$\tau_{i j}=2 \mu_t S_{i j}-\frac{2}{3} \rho k \sigma_{i j}$                  (9)

$\mu_{e f f}=\mu+\mu_t, \mu_t=\rho C_\mu \frac{k^2}{\varepsilon}$                     (10)

where: $C_\mu=0.0845, a_k=a_{\varepsilon}=1.39, C_{1 \varepsilon}=1.42, C_{2 \varepsilon}=1.68$

B- Boussinesq approximation.

$\overline{u_i^{\prime} u_j^{\prime}}=2 v_t S_{i j}-\frac{2}{3} k \sigma_{i j}$                   (11)

$-\overline{u_j^{\prime} T_a^{\prime}}=\frac{C p v_t}{P r_t} \cdot \frac{\partial T_a}{\partial x_j}$ and $\mathrm{S}_{\mathrm{ij}}=\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i} / 2$                   (12)

where: S_ij is the mean strain rate tensor.