*Corr. Author’s Address: University of Maribor, Faculty of Mechanical Engineering, Smetanova 11, 2000 Maribor, Slovenia, matjaz.ramsak@um.si 517
Strojniški vestnik - Journal of Mechanical Engineering 68(2022)9, 517-528 Received for review: 2022-01-24
© 2022 The Authors. CC BY 4.0 Int. Licensee: SV-JME Received revised form: 2022-04-25
DOI:10.5545/sv-jme.2022.28 Original Scientific Paper Accepted for publication: 2022-06-13
“por” — 2022/8/30 — 11:38 — page 1 — #1
Strojniški vestnik - Journal of Mechanical Engineering 63(2017)3, XXX-4
© 2017 Journal of Mechanical Engineering. All rights reserved.
Review Paper — DOI: 10.5545/sv-jme.2017.4027
Received for review: 2016-11-04
Received revised form: 2017-01-14
Accepted for publication: 2017-02-12
FractalGeometryasanEffectiveHeatSink
Matjaž Ramšak
1
1
Faculty of Mechanical Engineering, University of Maribor, Slovenia
"How long is the coast of Britain?" was the question stated by Mandelbrot. Using smaller and smaller rulers the coast length limits to infinity. If this logic is
applied to the fractal heat sink geometry, infinite cooling capacity should be obtained using fractals with mathematically infinite 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 noninfinite limit. The presented fractal heat sink geometry is non-competitive to the straight fins.
Keywords: fractal heat sink, light-emitting diode and central processing unit cooling, conjugate heat transfer, laminar flow, boundary element
method, Koch snowflake
Highlights
• Infinite fractal heat transfer area leads to zero Nusselt number. Their product limits to finite value of heat flow.
• Richardson extrapolation of numerical simulation results confirm that.
• Fractal flow pattern follows fractal geometry.
• Fractal heat sink could not compete with straight fins.
• Conjugate heat transfer is simulated using boundary element method inhouse code.
0 INTRODUCTION
Theconvectiveheattransferperunittime
˙
Qbetweena
heat sink and the fluid is computed as
˙
Q=hAΔT , (1)
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).ThetestisperformedusingasequenceofKoch
snowflakefractalformationstartingfromastraightline
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
straightfinswhichisusednowadaysinmostcases. As
itisapopularslogan: thereisaroomforimprovement.
The geometry of heat sinks is subject of many
recent articles [2–6]. The fractals are not an unknown
topic in this subject. In most papers, fractal geometry
isusedtocharacterisethesurfaceroughness[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
laminarflowinamicrochannelisanalysednumerically
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–12].
Elaboratereviewarticleoffractalheatexchangers
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
fractalheatsinkhaslowerpressuredrop,moreuniform
temperaturefielddistribution,andhighercoefficientof
performancethanthatofthetraditionalhelicalchannel
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
*Corr. Author’s Address: Name of institution, Address, City, Country, xxx.yyy@xxxxxx.yyy 1

Strojniški vestnik - Journal of Mechanical Engineering 68(2022)9, 517-528
518 Ramšak, M.
“por” — 2022/8/30 — 11:38 — page 2 — #2
Strojniški vestnik - Journal of Mechanical Engineering 63(2017)3, XXX-4
fractals is the subject of very recent work of He et
al.[18].
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,
coincidentallythesameasinourwork.
Table 1. Results sensitivity on mesh density for Re= 1 and solid fluid
heat conductivity ratiok=1
Mesh data
Mesh name coarse medium fine finest
Number of Nodes [*1000] 47 92 231 457
Nodes on Koch shortest element 3357
Nodes on Outlet 90 180 226 300
Results for isothermal computation
Number of iterations 940 880 955 940
CPU [h] (serial run) 14 36 100 351
AVG (Interface vorticity) -2.872 -2.779 -2.682 -2.662
Error to coarser mesh [%] - 3.35 3.62 0.75
Num. acc. estimation [1] [%] - - 1.43 0.50
Results for thermal computation
AVG (Interface temp.) 0.6536 0.6533 0.6529 0.6528
Error to coarser mesh [%] - 0.05 0.06 0.02
Num. acc. estimation [1] [%] - - 0.04 0.02
Nusselt number 0.5173 0.5206 0.5211 0.5218
Error to coarser mesh [%] - 0.62 0.10 0.13
Num. acc. estimation [1] [%] - - 0.08 0.16
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
functionderivativeinnormaldirectiontotheboundary
element. In this manner the problem of undefined
normal direction on the corner flux nodal points is
elegantly avoided. The main advantage of subdomain
techniqueissparsematrixincomparisontotheclassic
BEM, where only the boundary of computational
domainmustbediscretised. Thesubdomaintechnique
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
iterativeleastsquaresmethod,[22]. Thecodehasbeen
successfully used and validated for the conjugate heat
transferBenchmarkrevision[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
remarksonenhancedheatsinksusingfractalgeometry.
1 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
domainwithconstantvelocity. Thezerogradientoutlet
boundary condition is prescribed. It is not physically
adequatesincetheflowfieldisnotfullydeveloped,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-sectionshowninFig. 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
temperatureis1.0andinlettemperatureiszero. Inthis
manner the temperature differenceΔT is defined to be
1inEq. (11). Computingthesteadystatesolution,the
fluid solid thermal diffusivity ratio is the last solution
parameterinvestigatedinnextsections.
Thenon-dimensionalformofgoverningequations
for a 2D incompressible laminar flow are written
using the nondimensional stream function vorticity
formulation of Navier-Stokes equations. Stream
functionequationψ is
∂
2
ψ
∂x
2
+
∂
2
ψ
∂y
2
=−ω . (2)
2 Author’s Name Surname, Co-author’s Name Surname

Strojniški vestnik - Journal of Mechanical Engineering 68(2022)9, 517-528
519 Fractal Geometry as an Effective Heat Sink
“por” — 2022/8/30 — 11:38 — page 3 — #3
Strojniški vestnik - Journal of Mechanical Engineering 63(2017)3, XXX-4
Inlet: v=(1,0),
ψ =y,
T =0.
Outlet:
∂ψ
∂n
=0,
∂T
∂n
=0.
Free flow: v=(1,0),
ψ =y ,
∂ψ
∂n
=1,
T =0.
Interface:
 v=(0,0),
ψ =0,
∂ψ
∂n
=0,
T
s
=T
f
,
k
s
∂T
s
∂n
=k
f
∂T
f
∂n
.
Lower wall: T =1.


x
y
Fluid: Re=
1×1
ν
,
Pr=
ν
α
=0.71,
k=
k
s
k
f
.
Left/Right
solid wall:
∂T
∂n
=0.
 
1




0.1
0.9
Fig. 1. The geometry of computational domain and boundary conditions
Table 2. Results sensitivity on solid fluid conductivity ratiok for Re=1
Solid fluid conductivity ratiok 1 10 100 1000 10
4
10
5
AVG (Interface temperature) 0.6529 0.9175 0.9900 0.9990 0.9999 1.0000
Change to prior - lowerk [%] - 28.84 7.32 0.90 0.09 0.01
AVG (Outlet temperature) 0.4030 0.5236 0.5526 0.5562 0.5565 0.5565
Change to prior - lowerk [%] - 23.03 5.25 0.65 0.05 0.00
AVG (Nusselt number) 0.5211 1.1167 1.3659 1.3994 1.4033 1.4036
Change to prior - lowerk [%] - 53.33 18.25 2.39 0.28 0.02
Heat flow
˙
Q [W] 0.0309 0.0663 0.0811 0.0831 0.0833 0.0833
Change to prior - lower k [%] - 53.33 18.25 2.39 0.28 0.02
Heat transfer coef. h [W/m
2
K] 0.0073 0.0157 0.0192 0.0197 0.0198 0.0198
Change to prior - lower k [%] - 53.33 18.25 2.39 0.28 0.02
Table 3. Results sensitivity on Reynolds number fork=10
Reynolds numberRe 1 10 100
AVG (Interface vorticity) -2.682 -2.777 -3.629
AVG (Interface temperature) 0.9175 0.9054 0.8471
AVG (Outlet temperature) 0.5236 0.4766 0.3331
Nusselt number 1.1167 1.2231 1.7366
Heat flow
˙
Q [W] 0.0663 0.0726 0.1031
Heat transfer coef. h [W/m
2
K] 0.0157 0.0172 0.0245
Vorticity equationω is
∂ω
∂t
+
∂(v
x
ω)
∂x
+
∂(v
y
ω)
∂y
=
1
Re
∇
2
ω , (3)
wherev
x
isthevelocityinxdirectioncomputedasv
x
=
∂ψ/∂y and v
y
as v
y
=−∂ψ/∂x. The energy equation
within the fluid region is
∂T
∂t
+
∂(v
x
T)
∂x
+
∂(v
y
T)
∂y
=
1
RePr
∇
2
T , (4)
where T is non-dimensional fluid temperature.
Energy equation within the solid region is
∂T
∂t
=

α
s
α
f

1
RePr
∇
2
T , (5)
where α
s
and α
f
are diffusivities for the solid and
fluid regions respectively. For details on equations
derivation [23]. The mechanism of heat conduction
anditsbackgroundisclearlyexplainedinLiuetal.[9].
The interface boundary conditions on the wall
between solid and fluid are written as temperature
equality
T
f
=T
s
, (6)
and heat flux equality as
k
f

∂T
∂n

f,interface
=−k
s

∂T
∂n

s,interface
, (7)
where the n is unit normal direction to the fluid solid
interface. The local Nusselt number is defined as the
Fractal Geometry as an Effective Heat Sink 3

Strojniški vestnik - Journal of Mechanical Engineering 68(2022)9, 517-528
520 Ramšak, M.
“por” — 2022/8/30 — 11:38 — page 4 — #4
Strojniški vestnik - Journal of Mechanical Engineering 63(2017)3, XXX-4
Coarse N = 47,000
Finest N = 457,000
Fig. 2. Meshes used; Increasing the mesh density the number of outlet nodes are 90, 180, 226 and 300 while on the shortest fractal element of the Koch
boundary are 3, 3, 5 and 7 nodes
Table 4. Average values at the Koch fluid solid interface for fractal
geometry sequence; the fractal element length is denoted using l and
the complete interface length using L; the Richardson extrapolation [1]
to infinite boundary is presented in the line l = 0 including extrapolation
accuracy estimation interval. In the last line, the best guestimate values
are stated
lLω
Koch
T
Koch
Nu
˙
Qh
1.000 1.000 -18.58 0.9999 8.158 0.1149 0.1149
0.333 1.333 -15.52 0.9998 6.071 0.1140 0.0855
0.111 1.778 -8.958 0.9998 4.980 0.1247 0.0701
0.037 2.370 -6.397 0.9998 3.954 0.1320 0.0557
0.012 3.160 -4.578 0.9998 3.051 0.1358 0.0430
0.004 4.214 -3.629 0.9998 2.332 0.1384 0.0328
0.0 ∞ -2.518 0.9998 -0.899 0.1446 -0.0127
acc. est.± 1.0 0.0000 1.6 0.0081 0.0219
guestimate -2.518 0.9998 0 0.1446 0
temperaturederivativeatthefluidsideofthefluidsolid
interface as
Nu=−
∂T
∂n




f,interface
. (8)
The average Nusselt number Nu is the integral value
computed as
Nu=
1
L

L
0
Nudl , (9)
where L is the interface wall length. The heat flow
˙
Q
is computed using its definition as
˙
Q=−Ak
f

∂T
∂n

f
=(L·1)k
f
Nu, (10)
where A is the actual interface area computed as L·1.
Theunitylengthisdefinedonthe3rddimension. Heat
transfer coefficient h is computed as
h=
˙
Q
AΔT
(11)
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
4 Author’s Name Surname, Co-author’s Name Surname

Strojniški vestnik - Journal of Mechanical Engineering 68(2022)9, 517-528
521 Fractal Geometry as an Effective Heat Sink
“por” — 2022/8/30 — 11:38 — page 5 — #5
Strojniški vestnik - Journal of Mechanical Engineering 63(2017)3, XXX-4
Re=1 Re=100
Re=1 Re=100
Fig. 3. The Reynolds number dependency of the stream function and vorticity contour plots
parameter, arising from interface boundary condition
Eq. (7), influencing the solution.
Allgoverningequationscanbewritteninthesame
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.
2 RESULTSANDDISCUSSION
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.
2.1 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 Fig. 2. In Tab. 1 the results are shown for the
selected case Re= 1 and solid fluid conductivity ratio
k=1.
Three integral values are selected as numerical
solution accuracy indicators: average (AVG) value
of vorticity, temperature and Nusselt value on the
solid fluid interface. The AVG (Interface vorticity)
Fractal Geometry as an Effective Heat Sink 5

Strojniški vestnik - Journal of Mechanical Engineering 68(2022)9, 517-528
522 Ramšak, M.
“por” — 2022/8/30 — 11:38 — page 6 — #6
Strojniški vestnik - Journal of Mechanical Engineering 63(2017)3, XXX-4
Re=1
Re=100
Fig. 4. Geometry and flow (stream function) self similarity at Reynolds numbers 1 and 100
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
meshforallfurthercomputationsresultinginlessthan
1 % numerical solution accuracy. Similar numerical
accuracyispublishedinthework[23]whereusingthe
sameBEMcodethebenchmarkconjugateheattransfer
problemwassolved.
2.2 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
darkbluecolourrepresentingnegativestreamvalueon
the downstream side. The dark blue zones represent
theclockwiseandlightbluecounterclockwiserotation
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,
representingthezero-valuevorticitycontour.
The flow self-similarity is discussed in the next
paragraph. StreamlinesareshowninFig. 4,wheretwo
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. Theflowpatternself-similarityisobvioustoo,
especially for Re=1. The flow pattern is discussed in
many aspects in our latest work completely dedicated
totheisothermalsolution[25].
2.3 Influence of Solid Fluid Conductivity Ratiok
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
. The Nusselt
6 Author’s Name Surname, Co-author’s Name Surname

Strojniški vestnik - Journal of Mechanical Engineering 68(2022)9, 517-528
523 Fractal Geometry as an Effective Heat Sink
“por” — 2022/8/30 — 11:38 — page 7 — #7
Strojniški vestnik - Journal of Mechanical Engineering 63(2017)3, XXX-4
k=1 k=10 k=100
Fig. 5. Influence of solid fluid conductivity ratiok on temperature field forRe=1
Re=1 Re=10 Re=100
Fig. 6. Influence of Reynolds number on temperature field fork=10
values are not changing significantly (only 2.39 %)
when increasing k over 100, see Tab. 2 and Fig.
5, indicating no significant changes. Increasing k,
solid gradients decrease obtaining an almost constant
solid temperature. Replacing the steel with Al the
cooling heat flow increases for 0.3 % in this problem
configuration. The neglecting improvement of 0.02 %
isobtainedreplacingAlwithCu.
Interesting temperature contours are obtained for
k = 1, equalling the heat conductivity in fluid and
solid, see Fig. 5. In this case (Re = 1), the
soliddomaininfluenceontemperaturenearlyvanishes,
since the conduction equals convection, obtaining
largetemperaturegradientsinthesoliddomain.
The Biot number (Bi = Length ·h/k
s
) analysis
follows. The characteristic Length is defined in this
case as Length = A/L where the A is finite heat sink
cross section approximately 0.1 ·1 and L the length
of the fractal cooling surface area. Increasing L to
the infinity the characteristic length and Bi limits to
zero. Additionally, increasing k and consequently k
s
values the Bi number also limits to zero value. Both
parameters clearly indicating very low Bi values and
consequently the uniform solid temperature as already
mentioned. Using results in Tab. 2 the maximum
Bi number value is Bi = 0.0002 using L = 4.214,
Length = 0.0237, h = 0.0073 and k = 1, confirming
that heat conduction in solid prevails heat convection
tofluid.
2.4 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
Fractal Geometry as an Effective Heat Sink 7

Strojniški vestnik - Journal of Mechanical Engineering 68(2022)9, 517-528
524 Ramšak, M.
“por” — 2022/8/30 — 11:38 — page 8 — #8
Strojniški vestnik - Journal of Mechanical Engineering 63(2017)3, XXX-4
−20
−15
−10
−5
0
0.001 0.01 0.11
ω
Koch
Fractal element length (l)
Simulations results
Extrapolated value tol =0
0.11
0.115
0.12
0.125
0.13
0.135
0.14
0.145
0.15
0.001 0.01 0.11
˙
Q
Fractal element length (l)
Simulation results
Extrapolated value tol =0
−4
−2
0
2
4
6
8
10
0.001 0.01 0.11
Nu
Fractal element length (l)
Simulation results
Extrapolated value tol =0
guestimate
−0.04
−0.02
0
0.02
0.04
0.06
0.08
0.1
0.12
0.001 0.01 0.11
h
Fractal element length (l)
Simulation results
Extrapolated value tol =0
guestimate
Fig. 7. Graphical presentation of tabulated results in Tab. 4; the shaded area is accuracy estimation of extrapolated result value computed using [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.
2.5 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
andcreating4newelementsusingtheKochsnowflake
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
finalboundarylengthinthisresearchis(4/3)
5
=4.214
using 5th iteration formation and fractal element size
l =( 1/3)
5
= 0.004. Using Richardson extrapolation
[1],thenumericalresultvalueslimitiscomputedusing
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 flow
˙
Q
dynamics the smooth response is obtained generally.
The only exception is in the first iteration between
l=1.0andl=0.333. Thefirstquestionwouldbehow
itispossible,thatthecoolingheatflowisslightlylower
for 0.8 % using a single rib (l = 0.333) comparing
to the flat surface (l = 1.0). While the surface area
8 Author’s Name Surname, Co-author’s Name Surname

Strojniški vestnik - Journal of Mechanical Engineering 68(2022)9, 517-528
525 Fractal Geometry as an Effective Heat Sink
“por” — 2022/8/30 — 11:38 — page 9 — #9
Strojniški vestnik - Journal of Mechanical Engineering 63(2017)3, XXX-4
l=1.0 l=0.333 l=0.111
l=0.037 l=0.012 l=0.004
Fig. 8. Stream function contours for the Koch snowflake formation procedure denoted by fractal element lengthl; the flow field have the same fractal nature
as the geometry
of ribbed surface A is increased, the heat transfer
coefficient h is significantly decreased resulting in
slightly lower heat flow which is their product
˙
Q=
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
zoneperformsthermalisolationincreasingthethermal
boundarylayerthicknessindownwindareacomparing
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 onlyvery
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
equivalenttotheN11RoughnessGradeNumberwhich
isobtainedusingsandcastorhotrollmanufacturingof
heat sink, [26]. Finally, the minimal edge dimension
ofl=40μmisstillbigenoughtobeinthecontinuum
mechanics having a Knudsen number value of 600.
3 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:
• Theareaofinterfacesurfaceconvergestoinfinity.
• Nu and h converged to zero setting the
convective heat flow to zero (bearing in mind the
extrapolation error).
• The
˙
Q is increasing monotonically to a specific
finite value: heat diffusion.
Fractal Geometry as an Effective Heat Sink 9

Strojniški vestnik - Journal of Mechanical Engineering 68(2022)9, 517-528
526 Ramšak, M.
5  REFERENCES
[1] AIAA (1994). Editorial policy statement on numerical accuracy 
and experimental uncertainty, AIAA Journal, vol. 32, no. 1., p. 
3-8, DOI:10.2514/3.48281.
[2] Močnik, U., Blagojevič, B., Muhič, S. (2020). Numerical 
analysis with experimental validation of single-phase fluid flow 
in a dimple pattern heat exchanger channel. Strojniški vestnik 
- Journal of Mechanical Engineering, vol. 66,  no. 9, p. 544-
553, DOI:10.5545/sv-jme.2020.6776.
[3] Yadav, S., Pandey, K.M. (2018). A parametric thermal 
analysis of triangular fins for improved heat transfer in 
forced convection. Strojniški vestnik - Journal of Mechanical 
Engineering, vol. 64, no. 6, p. 401-411, DOI:10.5545/sv-
jme.2017.5085.
[4] Al-Kouz, W.G., Kiwan, S., Alkhalidi, A., Sari, M., Alshare, A. 
(2018). Numerical study of heat transfer enhancement for 
“por” — 2022/8/30 — 11:38 — page 10 — #10
Strojniški vestnik - Journal of Mechanical Engineering 63(2017)3, XXX-4
l =1.0 l =0.333 l =0.111
l =0.037 l =0.012 l =0.004
Fig. 9. Vorticity function contours for the Koch snowflake formation procedure denoted by fractal element lengthl
In this manner the resulting heat transfer Eq. (11)
could be written as
˙
Q=lim
l→0

(h→0)(A→∞)(ΔT =1)

(12)
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
fractalmannerornot,andthattheheatconvectionfrom
solidsurfacetofluidisachievedbyheatdiffusiononly
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.
4 ACKNOWLEDGEMENTS
Authors acknowledge the financial support from the
Slovenian Research Agency (research core funding
No. P2-0196).
5 REFERENCES
[1] AIAA, Editorial policy statement on numerical accuracy
and experimental uncertainty, AIAA Journal 32(1)
(1994) 3.
[2] Moˇ cnik, U., Blagojeviˇ c, B., and Muhiˇ c, S., Numerical
analysis with experimental validation of single-phase
fluid flow in a dimple pattern heat exchanger channel,
Strojniški vestnik - Journal of Mechanical Engineering
66 (2020) 544.
[3] Yadav, S. and Pandey, K., A parametric thermal
analysis of triangular fins for improved heat transfer
in forced convection, Strojniški vestnik - Journal of
Mechanical Engineering 64 (2018) 401.
[4] Al-Kouz, W., Kiwan, S., Alkhalidi, A., Sari, M.,
and Alshare, A., Numerical study of heat transfer
enhancement for low-pressure flows in a square cavity
10 Author’s Name Surname, Co-author’s Name Surname
“por” — 2022/8/30 — 11:38 — page 10 — #10
Strojniški vestnik - Journal of Mechanical Engineering 63(2017)3, XXX-4
l =1.0 l =0.333 l =0.111
l =0.037 l =0.012 l =0.004
Fig. 9. Vorticity function contours for the Koch snowflake formation procedure denoted by fractal element lengthl
In this manner the resulting heat transfer Eq. (11)
could be written as
˙
Q=lim
l→0

(h→0)(A→∞)(ΔT =1)

(12)
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
fractalmannerornot,andthattheheatconvectionfrom
solidsurfacetofluidisachievedbyheatdiffusiononly
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.
4 ACKNOWLEDGEMENTS
Authors acknowledge the financial support from the
Slovenian Research Agency (research core funding
No. P2-0196).
5 REFERENCES
[1] AIAA, Editorial policy statement on numerical accuracy
and experimental uncertainty, AIAA Journal 32(1)
(1994) 3.
[2] Moˇ cnik, U., Blagojeviˇ c, B., and Muhiˇ c, S., Numerical
analysis with experimental validation of single-phase
fluid flow in a dimple pattern heat exchanger channel,
Strojniški vestnik - Journal of Mechanical Engineering
66 (2020) 544.
[3] Yadav, S. and Pandey, K., A parametric thermal
analysis of triangular fins for improved heat transfer
in forced convection, Strojniški vestnik - Journal of
Mechanical Engineering 64 (2018) 401.
[4] Al-Kouz, W., Kiwan, S., Alkhalidi, A., Sari, M.,
and Alshare, A., Numerical study of heat transfer
enhancement for low-pressure flows in a square cavity
10 Author’s Name Surname, Co-author’s Name Surname

Strojniški vestnik - Journal of Mechanical Engineering 68(2022)9, 517-528
527 Fractal Geometry as an Effective Heat Sink
low-pressure flows in a square cavity with two fins attached 
to the hot wall using Al
2
O
3
-air nanofluid. Strojniški vestnik - 
Journal of Mechanical Engineering, vol. 64, no. 1, p. 26-36, 
DOI:10.5545/sv-jme.2017.4989.
[5] Čarnogurská, M., Příhoda, M., Lázár, M., Jasminská, N., 
Gallik, R., Kubik, M. (2016). Measuring selected parameters 
of polypropylene fibre heat exchangers. Strojniški vestnik - 
Journal of Mechanical Engineering, vol. 62, no. 6, p. 381-388, 
DOI:10.5545/sv-jme.2015.3202.
[6] Yayla, S. (2013). Flow characteristic of staggered multiple 
slotted tubes in the passage of a fin tube heat exchanger. 
Strojniški vestnik - Journal of Mechanical Engineering, vol. 59, 
no. 7-8, p. 462-472, DOI:10.5545/sv-jme.2012.902.
[7] Fong, R.W.L., McRae, G.A., Coleman, C.E., Nitheanandan, T., 
Sanderson, D.B. (2001). Correlation between the critical heat 
flux and the fractal surface roughness of zirconium alloy tubes. 
Journal of Enhanced Heat Transfer, vol. 8, no. 2, p. 137-146, 
DOI:10.1615/JEnhHeatTransf.v8.i2.60.
[8] Chen, Y., Zhang, C., Shi, M., Peterson, G.P. (2009). Role of 
surface roughness characterized by fractal geometry on 
laminar flow in microchannels. Physical Review E, vol. 80, art. 
ID 026301, DOI:10.1103/PhysRevE.80.026301.
[9] Liu, J., Xue, Y., Zhang, Q., Wang, H., Wang, S. (2022). 
Coupled thermo-hydro-mechanical modelling for geothermal 
doublet system with 3D fractal fracture. Applied Thermal 
Engineering, vol. 200, art. ID 117716, DOI:10.1016/j.
applthermaleng.2021.117716.
[10] Cai, W., Chen, W., Wang, F. (2017). Three-dimensional 
hausdorff derivative diffusion model for isotropic/anisotropic 
fractal porous media. Thermal Science, vol. 22,  no. suppl. 1, 
p. s1-s6, DOI:10.2298/TSCI170630265C.
[11] Gao, Y., Gao, F., Dong, G., Yan, W., Yang, X. (2019). The 
mechanical properties and fractal characteristics of the 
coal under temperature-gas-confining pressure. Thermal 
Science, vol. 23, no. suppl. 3, p. 789-798, DOI:10.2298/
TSCI180610094G.
[12] Wang, Q.-L., Li, Z.-B., Kong, H.-Y., He, J.-H. (2015). Fractal 
analysis of polar bear hairs. Thermal Science, vol. 19, no. 
suppl. 1, p. 143-144, DOI:10.2298/TSCI15S1S43W.
[13] Huang, Z., Hwang, Y., Aute, V., Radermacher, R. (2016). Review 
of Fractal Heat Exchangers. International Refrigeration and Air 
Conditioning Conference. Paper 1725.
[14] Ikegami, S., Umetani, K., Hiraki, E., Sakai, S., Higashino, S. 
(2018). Feasibility study of fractal-fin heat sink for improving 
cooling performance of switching power converters. IEEE 
International Telecommunications Energy Conference, p. 1-8, 
DOI:10.1109/INTLEC.2018.8612377.
“por” — 2022/8/30 — 11:38 — page 11 — #11
Strojniški vestnik - Journal of Mechanical Engineering 63(2017)3, XXX-4
l=1.0 l=0.333 l=0.111
l=0.037 l=0.012 l=0.004
Fig. 10. Temperature contours for the Koch snowflake formation procedure denoted by fractal element length l
with two fins attached to the hot wall using al2o3-air
nanofluid, Strojniški vestnik - Journal of Mechanical
Engineering 64 (2018) 26.
[5]
ˇ
Carnogurská,M.etal., Measuringselectedparameters
of polypropylene fibre heat exchangers, Strojniški
vestnik - Journal of Mechanical Engineering 62 (2016)
381.
[6] Yayla, S., Flow characteristic of staggered multiple
slotted tubes in the passage of a fin tube heat
exchanger, Strojniški vestnik - Journal of Mechanical
Engineering 59 (2013) 462.
[7] Fong, R., McRae, G., Coleman, C., Nitheanandan,
T., and Sanderson, D., Correlation between the
critical heat flux and the fractal surface roughness of
zirconium alloy tubes, JOURNAL OF ENHANCED
HEAT TRANSFER 8 (2001) 137+.
[8] Chen, Y., Zhang, C., Shi, M., and Peterson, G. P., Role
ofsurfaceroughnesscharacterizedbyfractalgeometry
on laminar flow in microchannels, Phys. Rev. E 80
(2009) 026301.
[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 Engineering200 (2022) 117716.
[10] Cai, W., Chen, W., and Wang, F., Three-dimensional
hausdorff derivative diffusion model for
isotropic/anisotropic fractal porous media, Thermal
Science 22 (2017) 265.
[11] Gao, Y., Gao, F., Dong, G., Yan, W., and Yang, X.,
The mechanical properties and fractal characteristics
of the coal under temperature-gas-confining pressure,
Thermal Science (2019) 94.
[12] Wang, Q.-L., Li, Z.-B., Kong, H.-Y., and He, J.-H.,
Fractal analysis of polar bear hairs, Thermal Science
19 (2015) 143.
[13] Zhiwei, H., Hwang, Y., Aute, V., and Radermacher, R.,
Reviewoffractalheatexchangers, in Review of Fractal
Heat Exchangers, 2016.
[14] Ikegami, S., Umetani, K., Hiraki, E., Sakai, S.,
and Higashino, S., Feasibility study of fractal-fin
heat sink for improving cooling performance of
switching power converters, 2018 IEEE International
Telecommunications Energy Conference (INTELEC)
(2018) 1.
[15] Xia,C.,Fu,J.,Lai,J.,Yao,X.,andChen,Z., Conjugate
heat transfer in fractal tree-like channels network heat
sink for high-speed motorized spindle cooling, Applied
Thermal Engineering90 (2015) 1032 .
Fractal Geometry as an Effective Heat Sink 11

Strojniški vestnik - Journal of Mechanical Engineering 68(2022)9, 517-528
528 Ramšak, M.
[15] Xia, C., Fu, J., Lai, J., Yao, X., Chen, Z. (2015). Conjugate 
heat transfer in fractal tree-like channels network heat 
sink for high-speed motorized spindle cooling. Applied 
Thermal Engineering, vol. 90, p. 1032-1042, DOI:10.1016/j.
applthermaleng.2015.07.024.
[16] Ahmed, H.E., Salman, B.H., Kherbeet, A.Sh., Ahmed, M.I. 
(2018). Optimization of thermal design of heat sinks: A review. 
International Journal of Heat and Mass Transfer, vol. 118, p. 
129-153, DOI:10.1016/j.ijheatmasstransfer.2017.10.099.
[17] Xu, S., Wang, W., Fang, K., Wong, C.-N. (2015). Heat transfer 
performance of a fractal silicon microchannel heat sink 
subjected to pulsation flow. International Journal of Heat 
and Mass Transfer, vol. 81, p. 33-40, DOI:10.1016/j.
ijheatmasstransfer.2014.10.002.
[18] He, Z., Yan, Y., Zhao, T., Zhang, L., Zhang, Z. (2021). Multi-
objective optimization and multi-factors analysis of the 
thermal/hydraulic performance of the bionic y-shaped fractal 
heat sink. Applied Thermal Engineering, vol. 195, art. ID 
117157, DOI:10.1016/j.applthermaleng.2021.117157.
[19] Zhang, S., Chen, X., Wu, Z., Zheng, Y. (2019). Numerical study 
on stagger koch fractal baffles micromixer. International 
Journal of Heat and Mass Transfer, vol. 133, p. 1065-1073, 
DOI:10.1016/j.ijheatmasstransfer.2019.01.009.
[20] Ramšak, M., Škerget, L. (1999). Mixed boundary elements for 
laminar flows. International Journal for Numerical Methods in 
Fluids, vol. 31, no. 5, p. 861-877, DOI:10.1002/(SICI)1097-
0363(19991115)31:5<861::AID-FLD900>3.0.CO;2-S.
[21] Žagar, I., Škerget, L. (1995). The numerical simulation of 
non-linear separation columns by boundary-domain integral 
formulation. Computers & Chemical Engineering, vol. 19, 
Suppl. 1, p. 785-790, DOI:10.1016/0098-1354(95)87130-6.
[22] Ramšak, M., Škerget, L. (2000). Iterative least squares 
methods for solving over-determined matrix for mixed 
boundary elements. ZAMM Berlin, vol. 80, no. 3, p. 657.
[23] Ramšak, M. (2015). Conjugate heat transfer of backward-
facing step flow: A benchmark problem revisited. International 
Journal of Heat and Mass Transfer, vol. 84, p. 791-799, 
DOI:10.1016/j.ijheatmasstransfer.2015.01.067.
[24] Ramšak, M., Škerget, L. (2014). A highly efficient multidomain 
bem for multimillion subdomains. Engineering Analysis 
with Boundary Elements, vol. 43, p. 76-85, DOI:10.1016/j.
enganabound.2014.03.009.
[25] Ramšak, M. (2019). Multidomain bem for laminar flow 
in complex fractal geometry. Engineering Analysis with 
Boundary Elements, vol. 101, p. 310-317, DOI:10.1016/j.
enganabound.2019.01.003.
[26] DeGarmo, E.P., Black, J.T., Kohser, R.A. (1979). Materials and 
processes in manufacturing, Macmillan, New York.