Fractal Geometry as an Effective Heat Sink

"How long is the coast of Britain?" was the question stated by Mandelbrot. Using smaller and smaller rulers the coast length limits to inﬁnity. If this logic is applied to the fractal heat sink geometry, inﬁnite cooling capacity should be obtained using fractals with mathematically inﬁnite surface area. The aim of this article is to check this idea using Richardson extrapolation of numerical simulation results varying the fractal element length from one to zero. As expected, the extrapolated cooling capacity has a noninﬁnite limit. The presented fractal heat sink geometry is non-competitive to the straight ﬁns.


INTRODUCTION
The convective heat transfer per unit timeQ between a heat sink and the fluid is computed aṡ where h is the heat transfer coefficient, A the heat transfer area, and ∆T the temperature difference between the solid surface and fluid free stream. If the heat transfer area tends to infinity, the heat transfer power should tend to infinity too, presuming that the heat transfer coefficient and ∆T are finite nonzero values. The aim of this research is to test this assumption using computational fluid dynamics (CFD). The test is performed using a sequence of Koch snowflake fractal formation starting from a straight line to some small finite value, as shown in the last figure in the article. Using Richardson extrapolation [1], the results are extrapolated to the infinite surface area.
If the result of this study confirms the idea, a very effective cooling device should be gained using fractal geometry.
The fractal geometry is encountered in nature in many shapes and purposes. For example, the blood veins with capillaries in lungs for heat and mass exchange. It is well known that the nature evolution solutions are superior to engineering ones in many areas. If this is true for heat sinks, then a fractal heat sink must exist to be better than simple one with the straight fins which is used nowadays in most cases. As it is a popular slogan: there is a room for improvement.
The geometry of heat sinks is subject of many recent articles [2][3][4][5][6]. The fractals are not an unknown topic in this subject. In most papers, fractal geometry is used to characterise the surface roughness [7], where the correlation between the critical heat flow and the fractal surface roughness of the pressure tubes from the cool heavy water moderator is investigated. A 3D laminar flow in a microchannel is analysed numerically in [8], where near wall swirls are obtained of size 1 µm. In the minority group of papers, fractal geometry is used to describe the material porosity, for example [9][10][11][12].
Elaborate review article of fractal heat exchangers is written in [13] where the list of similar articles is presented. Almost all of them are numerical simulations. The numerical feasibility study of fractal heat sink is published in [14], where fractal like heat sink is proposed for cooling of electronic device. The conjugate heat transfer in a fractal tree like channels network heat sink is studied numerically and experimentally [15]. A conclusion is made that a fractal heat sink has lower pressure drop, more uniform temperature field distribution, and higher coefficient of performance than that of the traditional helical channel net heat sink. In the latest review article [16] on optimization design of heat sinks the fractal geometry is mentioned only in one cited article [17]. A simple heat sink is simulated consisting of a fractal split microchannel up to second iteration formation only. The challenge of high cooling power in electronic devices and its dissipation using bionic Y-shaped
The Koch snowflake fractal geometry, again only up to the second iteration formation, is found in [19], where it is applied for micro mixer baffles geometry. The laminar flow is solved up to Re = 100, coincidentally the same as in our work. The conjugate heat transfer in this work is computed using in house code based on mixed boundary elements and subdomain boundary element method (BEM). The idea of mixed boundary elements [20] is to split the function and flux approximation using continuous interpolation polynomials for function and discontinuous for function derivative in normal direction to the boundary element. In this manner the problem of undefined normal direction on the corner flux nodal points is elegantly avoided. The main advantage of subdomain technique is sparse matrix in comparison to the classic BEM, where only the boundary of computational domain must be discretised. The subdomain technique in its limit version by [21], where each subdomain is consisted of three or four boundary elements as triangle or quadrilateral subdomain, resulted in extremely sparse system matrix like in finite element method (FEM). The interface boundary conditions between mixed elements of subdomains lead to overdetermined matrix, which is solved using fast iterative least squares method, [22]. The code has been successfully used and validated for the conjugate heat transfer Benchmark revision [23].
The paper is organised as follows. The Problem definition is stated after the Introduction. In the next section Results and Discussion, the main subsection is the last one titled The infinite heat flow idea. Prior to this, various tests are performed for the shortest simulated fractal length (1/3) 5 = 0.004 in the fifth iteration consisting of 4 5 = 1024 elements: mesh sensitivity, Reynolds dependency for isothermal solution, influence of solid/fluid thermal conductivity ratio and influence of Reynolds number value on thermal solution. The paper finishes with concluding remarks on enhanced heat sinks using fractal geometry.

PROBLEM DEFINITION
The geometry of the Koch snowflake represents the first cut out of the heat sink of the high heat density source, such as light-emitting diode (LED) or central processing unit (CPU) processor, as shown in Fig. 1. The bottom of the solid wall is heated to a constant temperature. The cooling fluid is flowing into the domain with constant velocity. The zero gradient outlet boundary condition is prescribed. It is not physically adequate since the flow field is not fully developed, but the obtained flow field is as expected at the outflow region and serves well for the numerical example aim. The buoyancy effect is neglected. This could be justified by small fluid solid temperature difference or orienting the gravity in the third dimension not influencing the flow field in the 2D computational cross-section shown in Fig. 1.
The problem is nondimensionalized as follows. Length quantities are nondimensionalized by the length of the computational domain in the mean flow direction. The velocity is nondimensionalized by the inlet velocity. The Reynolds number is computed as Re = Velocity@inlet × Length/ν. The steady laminar flow of air is presumed as a cooling fluid. The Prandtl number is set to 0.71. The bottom dimensionless temperature is 1.0 and inlet temperature is zero. In this manner the temperature difference ∆T is defined to be 1 in Eq. (11). Computing the steady state solution, the fluid solid thermal diffusivity ratio is the last solution parameter investigated in next sections.
The non-dimensional form of governing equations for a 2D incompressible laminar flow are written using the nondimensional stream function vorticity formulation of Navier-Stokes equations.
where v x is the velocity in x direction computed as v x = ∂ ψ/∂ y and v y as v y = −∂ ψ/∂ x. The energy equation within the fluid region is where T is non-dimensional fluid temperature. Energy equation within the solid region is where α s and α f are diffusivities for the solid and fluid regions respectively. For details on equations derivation [23]. The mechanism of heat conduction and its background is clearly explained in Liu et al. [9]. The interface boundary conditions on the wall between solid and fluid are written as temperature equality and heat flux equality as where the n is unit normal direction to the fluid solid interface. The local Nusselt number is defined as the The average Nusselt number Nu is the integral value computed as where L is the interface wall length. The heat flowQ is computed using its definition aṡ where A is the actual interface area computed as L · 1.
The unity length is defined on the 3rd dimension. Heat transfer coefficient h is computed as where ∆T is one in this case.
In this paper only a steady solution is computed. The steady solution is solid thermal diffusivity α s independent since the time derivative is zero, see Eq. (5). In this manner the solid fluid thermal conductivity ratio k defined as k = k s /k f is the only solid material parameter, arising from interface boundary condition Eq. (7), influencing the solution.
All governing equations can be written in the same general form and solved using practically the same multidomain BEM solver applying different boundary conditions for each governing equation. The solver is explained in a detail in [24]. The validation of the developed multidomain BEM solver [23] where the benchmark solution of conjugate heat transfer of backward facing step problem is computed.

RESULTS AND DISCUSSION
The basic aim of the article is to check the assumption about the infinite heat flow for infinite length of fractal solid-fluid interface. Before this main numerical test, a few necessary tests are performed using the finest fractal geometry which is numerically the most cumbersome to solve.

Mesh Sensitivity Study
The aim of this test is to choose the appropriate mesh density and numerical solution accuracy estimation using standard procedure described in [1]. Four mesh densities were used: coarse, medium, fine and finest, see  numerical solution accuracy is the worst among all, being 0.50 %, which is to be expected since the interface vorticity is the most nonlinear and therefore difficult to solve. The thermal solutions are less mesh sensitive, resulting in maximal 0.02 % error for average temperatures and 0.16 % error for the Nusselt number. Obviously and intuitively, the temperature profile over the interface is less mesh dependent than the vorticity one. Based on the CPU consumption and basic aim, the fine mesh was chosen as a default mesh for all further computations resulting in less than 1 % numerical solution accuracy. Similar numerical accuracy is published in the work [23] where using the same BEM code the benchmark conjugate heat transfer problem was solved.

Reynolds Number Dependency of Isothermal Solution
The contour plots are presented in Figs. 3 and 3 for stream function and vorticity field at Re = 1 and Re = 100. Close to the Koch snowflake boundary, the streamlines in Figs. 3 are almost symmetrical with respect to the upstream and downstream sides for Re = 1. For a higher Re number the symmetry feature is altered, since the recirculation zone appears.
In contrast to Re = 1, in general, the bright blue colour representing positive stream value change to dark blue colour representing negative stream value on the downstream side. The dark blue zones represent the clockwise and light blue counter clockwise rotation vortices. The symmetry issues are more evident in the smallest cavities. This feature is also visible on vorticity contour plots as a dark and bright red colour, representing the zero-value vorticity contour.
The flow self-similarity is discussed in the next paragraph. Streamlines are shown in Fig. 4, where two successive zooms are enlarged on the right-hand side of the full-scale plot. Zooms 1 and 2 were selected in such a manner that the geometric self-similarity is evident. The flow pattern self-similarity is obvious too, especially for Re = 1. The flow pattern is discussed in many aspects in our latest work completely dedicated to the isothermal solution [25].

Influence of Solid Fluid Conductivity Ratio k
The real solid fluid conductivity ratio k values are 16000, 9000 and 2300 for Cu, Al and steel heat sink material, respectively. The test values are chosen to be in the range from 1 to 10 5

Influence of Reynolds Number on the Thermal Solution
For this test, the fluid solid conductivity ratio k value is fixed to 10 in order to obtain temperature changes in the solid. In this academic case the heat sink is made using isolative material such as wood or plastic. Increasing the Re number, the heat convection prevails over diffusion, resulting in nearly linear growth of the Nusselt number values, see Tab. 3 and Fig. 6. One should expect that increasing Re and consequently the cooling heat rate, the outlet temperature should increase too. Wrong, the outlet temperature decreases since the mass flow  [1] increases. The proper thermodynamic conclusion in this case would be: the enthalpy of outflow, computed as mass flow and temperature product, increases and exergy decreases. Increasing Re number the interface temperature decreases too, indicating higher temperature gradients and consequently, higher heat flux in the solid domain.

The Infinite Heat Flow Idea
The infinite heat flow idea is tested by forming Koch snowflake fractal geometry starting from the flat heat sink geometry denoted by l = 1, where l is the fractal element length. The next geometry iteration is obtained by dividing each fractal element by factor 3 and creating 4 new elements using the Koch snowflake formation procedure, see Fig. 8. In this manner, the Koch boundary length is increased by factor 4/3 at each formation iteration giving an infinite limit. The final boundary length in this research is (4/3) 5 = 4.214 using 5th iteration formation and fractal element size l = (1/3) 5 = 0.004. Using Richardson extrapolation [1], the numerical result values limit is computed using l as mesh size parameter which tends to zero. The accuracy estimation is also a part of the extrapolation procedure results. The infinite heat flow assumption is tested using Re = 100 and solid fluid conductivity ratio k = 10 4 which matches the Aluminium heat sink material approximately.
In the Figs. 8, 9 and 10 the contour plots of stream function, vorticity and temperature are plotted respectively for each fractal element length l. It is interesting and natural to expect that the flow fields have the same fractal nature as the geometry has in the Koch snowflake formation procedure. This is more evident in the upwind side of the heat sink in comparison to the downwind side where the recirculation zone is present.
The quantitative results are presented in Tab. 4 and graphically in Fig. 7. Observing the heat flowQ dynamics the smooth response is obtained generally. The only exception is in the first iteration between l = 1.0 and l = 0.333. The first question would be how it is possible, that the cooling heat flow is slightly lower for 0.8 % using a single rib (l = 0.333) comparing to the flat surface (l = 1.0). While the surface area of ribbed surface A is increased, the heat transfer coefficient h is significantly decreased resulting in slightly lower heat flow which is their productQ = h · A · ∆T . In this case the cooling rib is more of a fluid flow and thermal obstacle than cooling enhanced as expected intuitively. Additionally, the recirculation zone performs thermal isolation increasing the thermal boundary layer thickness in downwind area comparing to the flat surface case, see Fig. 10. The next objectivity of discussion is to verify that the five iterations of the Koch snowflake formation are enough for the testing aim. Comparing results in the last two figures, namely l = 0.012 and l = 0.004, the macro flow solution does not change any more significantly. Decreasing l further approaching the roughness size value, the flow would change only very close to the boundary until the viscous forces would damp the smallest swirls in cavities limiting to the hydraulic smooth surface flow. If the dimension of the heat sink would be 1 cm, then the shortest fractal element length l = 0.004 would be 40 µm, which is equivalent to the N11 Roughness Grade Number which is obtained using sand cast or hot roll manufacturing of heat sink, [26]. Finally, the minimal edge dimension of l = 40 µm is still big enough to be in the continuum mechanics having a Knudsen number value of 600.

CONCLUSIONS
The infinite cooling capacity idea is very naive. The numerical experiment annulated the idea of course. Decreasing the fractal length l to zero the main conclusions are: • The area of interface surface converges to infinity.
• Nu and h converged to zero setting the convective heat flow to zero (bearing in mind the extrapolation error).
• TheQ is increasing monotonically to a specific finite value: heat diffusion. In this manner the resulting heat transfer Eq. (11) could be written aṡ as stated in Eq. (10). Since the Nusselt number represents the ratio between convection and diffusion, setting the Nu → 0 annulated heat convection leaving the diffusion the only heat transfer mechanism in the solid fluid interface, as it is the fact. The fact is also, that each real surface has a roughness, might be in the fractal manner or not, and that the heat convection from solid surface to fluid is achieved by heat diffusion only at the fluid solid interface. The numerical experiment in this article confirms this fact.
The fractal geometry heat sink as an effective heat sink? No. This kind of fractal heat sink is non-competitive to the simple straight fin heat sink. This is clearly revealed by almost constant fractal fin temperature for aluminium -air configuration in Tab. 2 and 4 indicating the optimal fin should be slenderer generating larger temperature changes. In this manner the resulting heat transfer Eq. (11) could be written aṡ as stated in Eq. (10). Since the Nusselt number represents the ratio between convection and diffusion, setting the Nu → 0 annulated heat convection leaving the diffusion the only heat transfer mechanism in the solid fluid interface, as it is the fact. The fact is also, that each real surface has a roughness, might be in the fractal manner or not, and that the heat convection from solid surface to fluid is achieved by heat diffusion only at the fluid solid interface. The numerical experiment in this article confirms this fact.
The fractal geometry heat sink as an effective heat sink? No. This kind of fractal heat sink is non-competitive to the simple straight fin heat sink. This is clearly revealed by almost constant fractal fin temperature for aluminium -air configuration in Tab. 2 and 4 indicating the optimal fin should be slenderer generating larger temperature changes.  [9] Liu, J., Xue, Y., Zhang, Q., Wang, H., and Wang, S., Coupled thermo-hydro-mechanical modelling for geothermal doublet system with 3d fractal fracture, Applied Thermal Engineering 200 (2022) 117716.
Fractal Geometry as an Effective Heat Sink