arXiv:1508.04596v1 [physics.flu-dyn] 19 Aug 2015

Large scale three-dimensional topology optimisation of
heat sinks cooled by natural convection
Joe Alexandersena,∗, Ole Sigmunda , Niels Aagea
a

Department of Mechanical Engineering, Solid Mechanics
Technical University of Denmark
Nils Koppels Allé, Building 404
DK-2800, Denmark

Abstract
This work presents the application of density-based topology optimisation
to the design of three-dimensional heat sinks cooled by natural convection.
The governing equations are the steady-state incompressible Navier-Stokes
equations coupled to the thermal convection-diffusion equation through the
Bousinessq approximation. The fully coupled non-linear multiphysics system
is solved using stabilised trilinear equal-order finite elements in a parallel
framework allowing for the optimisation of large scale problems with order of
40-330 million state degrees of freedom. The flow is assumed to be laminar
and several optimised designs are presented for Grashof numbers between
103 and 106 . Interestingly, it is observed that the number of branches in
the optimised design increases with increasing Grashof numbers, which is
opposite to two-dimensional optimised designs.
Keywords: topology optimisation, heat sink design, natural convection,
large scale, multiphysics optimisation
1. Introduction
Natural convection is the phenomena where density-gradients due to temperature differences cause a fluid to move. Natural convection is therefore
a natural way to passively cool a hot object, such as electronic components,
light-emitting diode lamps or materials in food processing.
∗

Corresponding author
Email address: joealex@mek.dtu.dk (Joe Alexandersen)

Preprint submitted to arXiv.org

August 20, 2015

Topology optimisation is a material distribution method [1] used to optimise the layout of a structure in order to minimise a given performance
measure subject to design constraints and a physical model. In order to take
convective heat transfer, to an ambient fluid, into account in the design process of density-based methods, a common extension is to introduce some form
of interpolation of the convection boundaries, see e.g. [2, 3, 4]. More recently,
Dede et al. [5] used these simplified models to design and manufacture heat
sinks subject to jet impingement cooling. However, it is hard to justify the
application of a predetermined and constant convection coefficient, because
topology optimisation often leads to unanticipated designs, closed cavities
and locally varying velocities. During the optimisation process, the design
also changes significantly and the interaction with the ambient fluid changes.
Therefore, to ensure physically correct capturing of the aspects of convective
heat transfer, the full conjugate heat transfer problem must be solved.
Topology optimisation for fluid systems began with the treatment of
Stokes flow in the seminal article by Borrvall and Petersson [6] and has
since been applied to Navier-Stokes [7], as well as scalar transport problems [8] amongst others. The authors have previously presented a densitybased topology optimisation approach for two-dimensional natural convection problems [9]. Recently, Coffin and Maute presented a level-set method
for steady-state and transient natural convection problems using X-FEM
[10]. Interested readers are referred to [9] for further references on topology
optimisation in fluid dynamics and heat transfer.
Throughout this article, the flows are assumed to be steady and laminar.
The fluid is assumed to be incompressible, but buoyancy effects are taken
into account through the Boussinesq approximation, which introduces variations in the fluid density due to temperature gradients. The inclusion of a
Brinkman friction term facilitates the topology optimisation of the fluid flow.
The scope of this article is primarily to present large scale three-dimensional
results using the formulation presented in [9], as well as the computational
issues arising from solving the non-linear equation system and the subsequent
linear systems. Thus, only a brief overview of the underlying finite element
and topology optimisation formulations are given and the reader is referred
to [9] for further information.
In recent years an increasing body of work has been published on efficient
large scale topology optimisation. These works cover the use of high-level
scripting languages [11, 12], multiscale/-resolution approaches [13, 14] and
parallel programming using the message parsing interface and C/Fortran
2

[15, 16, 17, 18]. To facilitate the solution to truly large scale conjugate heat
transfer problems, the implementation in this paper is done using PETSc
[19] and the framework for topology optimisation presented in [18].
The layout of the paper is as follows: Section 2 presents the governing
equations; Section 3 presents the topology optimisation problem; Section
4 briefly discusses the finite element formulation; Section 5 discusses the
numerical implementation details; Section 6 presents scalability results for
the parallel framework and optimised designs for a test problem; Section 7
finishes with a discussion and conclusion.
2. Governing equations
The dimensionless form of the governing equations have been derived
based on the Navier-Stokes and convection-diffusion equations under the assumption of constant fluid properties, incompressible, steady flow and neglecting viscous dissipation. Furthermore, the Boussinesq approximation has
been introduced to take density-variations due to temperature-differences
into account. A domain is decomposed into two subdomains, Ω = Ωf ∪ Ωs ,
where Ωf is the fluid domain and Ωs is the solid domain. In order to facilitate
the topology optimisation of conjugate natural convective heat transfer between a solid and a surrounding fluid, the equations are posed in the unified
domain, Ω, and the subdomain behaviour is achieved through the control of
coefficients. The following dimensionless composite equations are the result.
∀x ∈ Ω :


∂ui
∂
∂ui ∂uj
∂p
uj
− Pr
+
+
= −α(x)ui − GrP r2 egi T
(1)
∂xj
∂xj ∂xj
∂xi
∂xi
∂uj
=0
(2)
∂xj


∂
∂T
∂T
−
K(x)
= s(x)
(3)
uj
∂xj
∂xj
∂xj
where ui is the velocity field, p is the pressure field, T is the temperature field,
xi denotes the spatial coordinates, egi is the unit vector in the gravitational
direction, α(x) is the spatially-varying effective impermeability, K(x) is the
spatially-varying effective thermal conductivity, s(x) is the spatially-varying
volumetric heat source term, P r is the Prandtl number, and Gr is the Grashof
number.
3

The effective thermal conductivity, K(x), is defined as:

1 if x ∈ Ωf
K(x) = 1
if x ∈ Ωs
Ck

(4)

k

where Ck = kfs is the ratio between the fluid thermal conductivity, kf , and the
solid thermal conductivity, ks . Theoretically, the effective impermeability,
α(x), is defined as:

0 if x ∈ Ωf
α(x) =
(5)
∞ if x ∈ Ωs

in order to ensure zero velocities inside the solid domain. However, numerically this requirement must be relaxed as will be described in section 3. The
volumetric heat source term is defined as being active within a predefined
subdomain of the solid domain, ω ⊂ Ωs :

0 if x ∈
/ω
s(x) =
(6)
s0 if x ∈ ω

is the dimensionless volumetric heat generation, q is the
where s0 = ksqL
∆T
dimensional volumetric heat generation, ∆T is the reference temperature
difference and L is the reference length scale.
The Prandtl number is defined as:
ν
Pr =
(7)
Γ
where ν is the kinematic viscosity, or momentum diffusivity, and Γ is the
thermal diffusivity. It thus describes the relative spreading of viscous and
thermal effects. The Grashof number is defined as:
Gr =

gβ ∆T L3
ν2

(8)

where g is the acceleration due to gravity and β is the volumetric coefficient of
thermal expansion. It describes the ratio between the buoyancy and viscous
forces in the fluid. The Grashof number is therefore used to describe to what
extent the flow is dominated by natural convection or diffusion. For low
Gr the flow is dominated by viscous diffusion and for high Gr the flow is
dominated by natural convection. The problems in this article are assumed
to have large enough buoyancy present to exhibit natural convective effects,
but small enough Gr numbers to exhibit laminar fluid motion.
4

3. Optimisation formulation
3.1. Interpolation functions
In order to perform topology optimisation, a continuous design field, γ(x),
varying between 0 and 1 is introduced. Pure fluid is represented by γ(x) = 1
and solid by γ(x) = 0. For intermediate values between 0 and 1, the effective
conductivity is interpolated as follows:
K(γ) =

γ(Ck (1 + qf ) − 1) + 1
Ck (1 + qf γ)

(9)

and likewise the effective impermeability is interpolated using:
α (γ) = α

1−γ
1 + qα γ

(10)

The interpolation functions ensure that the end points defined in (4) and (5),
respectively, are satisfied. The effective impermeability is bounded to α in the
solid regions and this upper bound should be chosen large enough to provide
vanishing velocities, but small enough to ensure numerical stability. The
convexity factors, qf and qα , are used to control the material properties for
intermediate design values in order to promote well-defined designs without
intermediate design field values.
3.2. Optimisation problem
The topology optimisation problem is defined as:
Z
minimise:f (γ, T ) = s(x) T dV
γ∈D
Z ω
Z
subject to: g(γ) =
1 − γ dV ≤ vf
Ωd

dV

(11)

Ωd

R(γ, u, p, T ) = 0
0 ≤ γ(x) ≤ 1 ∀x ∈ Ωd

where γ is the design variable field, D is the design space, f is the thermal compliance objective functional, g is the volume constraint functional,
R(γ, u, p, T ) is the residual arising from the stabilised weak form of the governing equations, and Ωd ⊆ Ω is the design domain.
The design field is regularised using a PDE-based (partial differential
equation) density filter [20, 18] and the optimisation problem is solved using
the method of moving asymptotes (MMA) [21, 18].
5

3.3. Continuation scheme
A continuation scheme is performed on various parameters in order to
stabilise the optimisation process and to improve the optimisation results.
It is the experience of the authors that the provided continuation scheme
yields better results than starting with the end values, although this cannot
generally be proven [22, 23].
The chosen continuation strategy consists of five steps:
qf ∈ {0.881, 8.81, 88.1, 881, 881}
qα ∈ {8, 8, 8, 98, 998}

α ∈ 105 , 105 , 105 , 106 , 107

(12a)
(12b)
(12c)

The sequence is chosen in order to alleviate premature convergence to poor
local minimum. The value of qf is slowly increased to penalise intermediate design field values with respect to conductivity. The maximum effective
permeability, α, is set relatively low during the first three steps, as this ensures better scaled sensitivities and more stable behaviour. Over the last two
steps, α is increased by two orders of magnitude in order to further decrease
the velocity magnitudes in the solid regions. The particular values of qf are
chosen by empirical inspection such as to ensure the approximate collocation
of the fluid and thermal boundaries. A more comprehensive investigation
into the modelling accuracy of the boundary layers is outside the scope of
this paper and will be investigated in future work.
It is important to note that the optimisation problem is by no means
convex and any optimised design will at best be a local minimum. The obtained design will always depend on the initial design as well as the continuation strategy. However, in the authors experience, the chosen continuation
strategy gives a good balance between convergence speed, final design performance and physicality of the modelling. The effects of the steps of the
current continuation strategy on the design distribution will be discussed in
section 6.4.
4. Finite element formulation
The governing equations are discretised using stabilised trilinear hexahedral finite elements. The design field is discretised using elementwise constant variables, in turn rendering the effective thermal conductivity and the
effective impermeability to be elementwise constant. The monolithic finite
6

element discretisation of the problem ensures continuity of the temperature
field, as well as the fluxes across fluid-solid interfaces. The particularities of
the implemented finite element formulation are as detailed in [9]. However,
simpler stabilisation parameters have been used in order to allow for consistent sensitivities to be employed, see Appendix B. The Jacobian matrix is
now fully consistent in that variations of the stabilisation parameters with
respect to design and state fields are included, in contrast to the original
work [9].
5. Numerical implementation
The discretised FEM equation is implemented in PETSc [19] based on the
topology optimisation framework presented in [18]. The PETSc framework
is used due to its parallel scalability, the availability of both linear and nonlinear solvers, preconditioners and structured mesh handling possibilities. All
components described in the following are readily available within the PETSc
framework.
5.1. Solving the non-linear system
The non-linear system of equations is solved using a damped Newton
method. The damping coefficient is updated using a polynomial L2 -norm fit,
where the coefficient is chosen as the minimiser of the polynomial fit. The
polynomial fit is built using the L2 -norm of the residual vector at the current
solution point, at 50% of the Newton step and at 100% of the Newton step.
This residual-based update scheme combined with a good initial vector (the
solution from the previous design iteration) has been observed to be very
robust throughout the optimisation process for the moderately non-linear
problems treated. To further increase the robustness of the non-linear solver,
if the Newton solver fails from the supplied initial vector, a ramping scheme
for the heat generation magnitude is applied in order to recover. Throughout,
the convergence criteria for the Newton solver is set to a reduction in the L2 norm of the residual of 10−4 relative to the initial residual.
5.2. Solving the linear systems
Due to the large scale and three-dimensional nature of the treated problems, the arising linearised systems of equations are by far the most time

7

consuming part of the Newton scheme. Therefore, to make large scale problems tractable the (unsymmetric) linear systems are solved using a fully
parallelised iterative Krylov subspace solver.
Constructing an iterative solver that is both independent of problem settings and possesses both parallel and numerical scalability, is intricate and
beyond the scope of this work, where focus is on the optimisation. However,
in order to facilitate the solution of large scale optimisation problems, an
efficient solver is required. To this end the authors use the readily available
core components in PETSc to construct a solver with focus on simplicity, as
well as reduced wall clock time. This is quite easy to obtain for solvers and
preconditioners that rely heavily on matrix-vector multiplications. To this
purpose the flexible GMRES (F-GMRES) method is used as the linear solver
combined with a geometry-based Galerkin-projection multigrid (GMG) preconditioner. Such a solver depends highly on the quality of the smoother
to guarantee fast convergence [24]. The authors have found that a simple
Jacobi-preconditioned GMRES provides a reasonable choice of smoother1 .
The convergence criterion for the Krylov solver is set to 10−5 relative to the
initial residual.
6. Results
6.1. Problem setup
The considered optimisation problem is an (academic) example of a heat
sink in a closed cubic cavity. Figure 1 shows the problem setup with dimensions and boundary conditions. The heat source (black in figure), exemplifying an electronics chip, is placed in the mid-bottom of the cavity and
modelled using a small block of solid material with volumetric heat generation, s0 = 104 . The design domain (gray in figure) is placed on top of the
heat source in order to allow the cooling fluid to pass underneath it, as well
as to allow room for wires etc. The vertical and top outer walls of the cavity
are assumed to be kept at a constant cold temperature, T = 0, while the
bottom is insulated. The height of the cavity has been used as the reference
length scale, L. Thus, the cavity dimensions are 1 × 1 × 1, the design domain dimensions are 0.75 × 0.75 × 0.75 and the heat source dimensions are
1

Due to the choice of a Krylov smoother, the multigrid preconditioner will vary slightly
with input and thus require a flexible outer Krylov method.

8

(a) Dimensions

(b) Boundary conditions

Figure 1: Illustration of the problem setup for the heat sink in a close cubic
cavity. The heat source is black and the design domain is gray.
0.1 × 0.05 × 0.1. A discussion on the definition of the temperature difference
in the Grashof number can be found in Appendix A. Initial investigations
showed that due to the symmetry of the domain and boundary conditions,
the design and state solutions remained quarter symmetric throughout the
optimisation. Thus, the computational domain is limited to a quarter of the
domain with symmetry boundary conditions. The volume fraction is set to
5%, i.e. vf = 0.05, for all examples.
6.2. Parallel performance
To demonstrate the parallel performance of the state solver, the model
optimisation problem described in section 6 is solved on a fixed mesh for different numbers of processes. All computations in this paper were performed
on a cluster, where each node is equipped with two Intel Xeon e5-2680v2 10core 2.8GHz processors. The results shown in Tables 1 and 2 are averaged
over 250 design cycles and show the performance for Gr = 103 and Gr = 106 ,
respectively. The data shows that the proposed solver scales almost linearly
in terms of speed up, and more importantly that the performance is only
slightly affected by the Grashof number.
In order to quantify the degree of numerical scalability, a second study is
performed in which the mesh resolution is varied. The study is conducted for
9

Processes
160
320
640

time [s]
53.2
28.9
14.1

scaling
1.00
0.54
0.26

Table 1: Average time taken per state solve over 250 design iterations for
Gr = 103 at a mesh resolution of 80 × 160 × 80.

Processes
160
320
640

time [s]
62.6
31.9
16.5

scaling
1.00
0.51
0.26

Table 2: Average time taken per state solve over 250 design iterations for
Gr = 106 at a mesh resolution of 80 × 160 × 80.

Mesh size
Gr = 103
80 × 160 × 80
7.5
160 × 320 × 160
10.1
320 × 640 × 320
18.4

Gr = 106
5.6
7.7
15.6

Table 3: Average iterations for the linear solver per state solve over entire
design process for Gr = 103 and Gr = 106 at varying mesh resolutions.

10

both low and high Grashof numbers and the average F-GMRES iterations
are collected in Table 3. The total number of design iterations averaged
over was 250, 500 and 1000, respectively, for the three mesh resolutions.
The data clearly shows that the computational complexity increases with
problem size, and thus that the solver is not numerically scalable. However,
since the growth in numerical effort, i.e. the number of F-GMRES iterations,
is moderate we conclude that the proposed solver is indeed applicable for
solving large scale natural convection problems.
6.3. Varying Grashof number
The problem is investigated for varying Gr under constant volumetric
heat generation, s0 = 104 , Prandtl number, P r = 1, and thermal conductivity ratio, Ck = 10−2 . The purpose of the present study is not to provide
a detailed physical example. It is rather to provide phenomenological insight into the effect of changing the governing parameter of the fluid-thermal
coupling, namely the Grashof number. However, physical interpretations of
tuning the Gr-number, while keeping the dimensionless volumetric heat generation and the P r-number constant, could be e.g. equivalent to tuning the
gravitational strength (going from microgravity towards full gravity) or the
dimensional volumetric heat generation. It is important to note that when
interpreting the results, the dimensional temperature scale will differ for the
two interpretations. While tuning the gravitational strength, the temperature scale remains the same; by varying the magnitude of the volumetric heat
generation, the temperature scale varies accordingly2 .
The computational mesh is 160 × 320 × 160 elements yielding a total
of 8, 192, 000 elements and 41, 603, 205 degrees of freedom for the quarter
domain. The design domain consists of 3, 456, 000 elements and the filter
radius is set to 2.5 times the element size.
Figure 2 shows the optimised designs for varying Gr-number with superimposed slices of the corresponding temperature fields. Due to the use of
density filtering, the interface between solid and fluid regions for the optimized designs are not exactly described but consists of a thin layer of intermediate design field values. The optimised design distributions are shown
The dimensionless volumetric heat generation is given by s0 = ksqL
∆T . Requiring that
s0 remains constant, means that the dimensional volumetric heat generation, q, and the
q
= const.. An increase
reference temperature difference, ∆T , must vary in unison, i.e. ∆T
in the Gr-number is thus achieved through the increase of ∆T .
2

11

(a) Gr = 103

(b) Gr = 104

(c) Gr = 105

(d) Gr = 106

Figure 2: Optimised designs for varying Gr-number at a mesh resolution of
160 × 320 × 160.

12

Gr
103
104
105
106

Primary
12
8
8
8

Secondary
0
16
28
48

Surface area
0.887
0.853
0.834
0.846

Table 4: The number of primary and secondary branches, as well as the
surface area for the optimised designs of figure 2.
thresholded at γ = 0.05, which is the approximate location of the computational interface. The general characteristics of all the designs are similar, i.e.
all designs are “thermal trees” with conductive branches moving the heat
away from the heat source. However, it can clearly be seen that the design
changes considerably with increasing convection-dominance (increasing Gr).
For increasing Gr-number the conductive branches contract, resulting in a
smaller spatial extent of the overall heat sink. This intuitively makes sense as
the problem goes from one of conduction/diffusion at Gr = 103 to convection
at Gr = 106 . When diffusion dominates, the goal for the branches essentially
becomes to conduct the heat directly towards the cold walls. As convection
begins to matter, the fluid movement aids the transfer of heat away from the
heat sink and the branches do not need to be as long. Instead, the design
forms higher vertical interfaces in order to increase surface area perpendicular to the flow direction and thus better transfer of heat to the fluid. At the
same time the complexity of the designs increases as the importance of convection increases. This can be seen by studying the number of primary and
secondary branches. Primary branches are defined as those connected to the
heat source directly and secondary branches are those connected to primary
branches. Table 4 shows the number of primary and secondary branches of
the optimised designs shown in figure 2. The number of primary branches is
largest for the diffusion-dominated case, but more or less constant thereafter.
However, it can be seen that the number of secondary branches significantly
increases as the Gr-number is increased. Table 4 also shows that the total
surface area3 is decreasing for the three lower Gr-numbers and then increases
slightly.
3

The surface area is computed using an isosurface at the selected threshold design field
value.

13

(a) Gr = 103

(b) Gr = 106

Figure 3: Temperature distribution in optimised designs for Gr = 103 and
Gr = 106 at a mesh resolution of 160 × 320 × 160 - view from below.
It is interesting to note, that the trend of increasing complexity with
increasing Gr-number is the reverse of what was observed for two-dimensional
problems [9]. There, the complexity of the design decreased as the Gr-number
increased. This difference is likely due to the fact that going into threedimensions allows the fluid to move around and through the design making
it more a question of “topology”, in contrast to in two-dimensions where it is a
question of surface shape. Physically, in two-dimensions additional branches
disturb the flow and thus the heat transfer; in three-dimensions, vertical
branches improve the heat transfer through an increased vertical surface
area. Figure 3 shows the optimised designs for Gr = 103 and Gr = 106
from below. The radial extent of the designs are emphasised from these
views. It is also seen that the branches for Gr = 106 are positioned to ensure
that the structure is open from below, with the branches forming vertical
14

Analysis Gr
103
104
105
106

3

10
2.064962
1.933460
1.488278
1.134963

Optimisation Gr
104
105
2.066856 2.241125
1.880946 1.994538
1.450276 1.404511
1.121979 1.061452

106
2.362470
2.111726
1.439089
1.025340

Table 5: Crosscheck objective function values for the designs shown in figure
2.
Gr
103
104
105
106

Time

Non-linear:
avg. (max)
9:56
1.9 (2)
10:25
2.0 (3)
10:28
2.1 (10)
10:35
2.1 (7)

Linear:
contin - avg. (max)
7.6,8.4,7.8,8.1,18.4 - 10.1 (25)
8.3,8.6,8.3,8.6,22.7 - 11.3 (29)
8.4,8.6,8.7,8.2,15.7 - 9.9 (34)
7.3,7.4,7.5,8.0,8.4 - 7.7 (14)

Table 6: Computational time, average non-linear iterations and linear iterations for the optimised designs of figure 2.
walls as discussed above.
Table 5 shows a crosscheck of the objective functions for the optimised
designs. It can be seen that all designs optimised for certain flow conditions
outperform the other designs for the specified Gr-number.
The optimisations was run for 500 design iterations for each optimised
design. Table 6 shows the computational time, average non-linear iterations
and linear iterations for the optimised designs of figure 2 using 1280 cores.
It can be seen that the computational time only weakly increases as the Grnumber is increased and remains close to 1.2 minutes per design iteration
on average. It is interesting to note, that Gr = 106 seems easier to solve
than the others as it exhibits lower average number of linear iterations than
all others, as well as a lower maximum number of non-linear iterations than
Gr = 105 . Furthermore, it can be seen that due to the Newton method
starting from a good initial vector, only 2 non-linear iterations are needed
for most of the design iterations independent of Gr-number.

15

6.4. High resolution design
The problem is now investigated with a computational mesh of 320 ×
640 × 320 elements yielding a total of 65, 536, 000 elements and 330, 246, 405
degrees of freedom for the quarter domain. The design domain consists of
27, 648, 000 elements and the filter radius is set to 2.5 times the element size,
i.e. in absolute measures half the size of before. The optimisation is run for
1000 design iterations and the computational time was 107 and 108 hours,
respectively, for Gr = 103 and Gr = 106 , using 2560 cores. This yields an
average of 6.4 and 6.5 minutes per design iteration, respectively.
Figure 4 shows the optimised design for Gr = 106 with the fine mesh
resolution and small length scale at various steps of the continuation strategy. The complexity of the design can be seen to be significantly higher than
for the design with a larger length scale, figure 2d. The subfigures are the
final iterations of the 1st, 2nd, 3rd and 5th (final) continuation steps. It
can be seen that the complexity of the design decreases during the optimisation process once the overall topology has been settled (iteration 400 and
onwards). This is due to the harder and harder penalisation of intermediate
design field values, with respect to conductivity, which are present at the
interface between solid and fluid. Therefore, smaller features are progressively removed as the surface area is more heavily penalised. The reason
for going to such high penalisation of the conductivity is to ensure the approximate collocation of the fluid and thermal boundaries. However, if one
starts directly with these physically-relevant parameters, particularly poor
local minima have been observed.
Figure 5 shows the final optimised designs for Gr = 103 and Gr = 106
with the fine mesh resolution and small length scale. The complexities of
both designs can be seen to be significantly higher than the previous, figures
2a and 2d, due to the smaller length scale. The increase in complexity with
increasing Grashof number is not as apparent at this resolution and length
scale, however, it is still argued that the number of primary and secondary
members follows the trend from the previous study, section 6.3. Also, the
shift from long conducting branches to shorter members, with significant
vertical surface area, is still obvious.
7. Discussion and conclusion
In this paper we have applied topology optimisation to the design of
three-dimensional heat sinks using a fully coupled non-linear thermofluidic
16

(a) Iteration 200

(b) Iteration 400

(c) Iteration 600

(d) Iteration 1000

Figure 4: Optimised designs for Gr = 106 at a mesh resolution of 320 ×
640 × 320. Please note that the freely hanging material in (a) is due to only
elements below the threshold, γ = 0.05, are shown - the design is connected
by intermediate design field values.

17

(a) Gr = 103

(b) Gr = 106

Figure 5: Optimised designs for Gr = 103 and Gr = 106 at a mesh resolution
of 320 × 640 × 320. The outline of the outer walls are shown in black.
model. In contrast to previous works, that considered simplified convection
models, the presented methodology is able to recover interesting physical
effects and insights, and avoids problems with the formation of non-physical
internal cavities, length-scale effects and artificial convection assumptions.
The implementation of the code in a PETSc framework suitable for large
scale parallel computations allows for running examples with more than 300
million degrees of freedom and almost 30 million design variables on regular
grids.
The example considered in the paper was primarily of academic nature.
Nevertheless, some interesting insight is obtained, showing that optimised
structures go from exhibiting simple branches that conduct heat towards the
cold outer boundaries for diffusion-dominated problems, towards complex
and compact multi-branched structures that maximize the convection heat
transfer for higher Grashof numbers.
Current and future work includes applications to real life problems (see
preliminary work in [25]), irregular meshes, multiple orientations, as well
as the extension to transient problems and detailed investigations into the
modelling accuracy of the boundary layer.

18

Computational Gr
103
104
105
106
103
106

Mesh
Figure Actual Gr
160 × 320 × 160
2a
1.69 × 103
160 × 320 × 160
2b
1.55 × 104
160 × 320 × 160
2c
1.16 × 105
160 × 320 × 160
2d
8.60 × 105
320 × 640 × 320
5a
1.42 × 103
320 × 640 × 320
5b
7.16 × 105

Table A.7: Computational and actual Gr-numbers for the designs presented
in this paper.
8. Acknowledgements
This work was funded by Villum Fonden through the NextTop project,
as well as by Innovation Fund Denmark through the HyperCool project. The
first author was also partially funded by the EU FP7-MC-IAPP programme
LaScISO.
Appendix A. Computational versus actual Grashof number
The Grashof numbers used above are all based on an a priori defined reference temperature difference. However, the actual temperature difference
observed between the heat source and the walls are not known a priori due
to the fact that the problem only has a single known temperature (Dirichlet
boundary condition) and a volumetric heat source. This is why the maximum temperature observed for the designs, see figure 2, is not equal to 1.
Therefore, after the optimisation, one can define an a posteriori Grashof
number based on the actual computed temperature difference. The a priori
version can be termed the computational Grashof number, the one that goes
into the dimensionless governing equations; and the a posteriori version can
be called the actual, or optimised, Grashof number for the given optimised
design. The actual Grashof number can be useful for experimental studies, as
well as for future comparisons. The actual Grashof numbers for the designs
shown in this paper are listed in Table A.7. The values for the fine meshes
are lower due to the lower thermal compliance, equivalent to the temperature
of the heat source, allowed by the smaller design length scale.

19

Appendix B. Stabilisation parameters
The stabilised weak form of the equations are given in [9]. The equations
are stabilised using the pressure-stabilising Petrov-Galerkin (PSPG) [26, 27]
and streamline-upwind Petrov-Galerkin (SUPG) methods [28], for more information please see [9].
The stabilisation parameters for SUPG and PSPG are defined as one and
the same using the following approximate min-function:
τSU = τP S = τ =



1
1
1
+ 2+ 2
2
τ1
τ2
τ3

− 21

(B.1)

The limit factors are given by:
4he
||ue ||2
he 2
τ2 =
4P r
1
τ3 =
αe

τ1 =

(B.2a)
(B.2b)
(B.2c)

where he is a characteristic element size (for cubes the element edge length)
and ue is the element vector of velocity degrees of freedom. The first limit
factor, τ1 , has been simplified based on evaluation at element centroids under the assumption of a single Gauss-point, yielding a constant stabilisation
factor within each element.
In order to define a consistent Jacobian matrix, and thus a consistent
adjoint problem, the derivatives of the stabilisation factors are needed with
respect to the velocity field. This can be found to be:
 2  2 !−1
−1 T
∂τ
τ1
τ1
ue T ue
= −τ 1 +
ue
(B.3)
+
∂ue
τ2
τ3
Furthermore, to define consistent design sensitivities, the derivatives of
the stabilisation factors with respect to the design field is needed. This can
be found to be:
!−1
 2  2
∂τ
τ
τ3
∂αe
τ3
=−
(B.4)
+
+1
∂γe
τ3
τ1
τ2
∂γe
20

Including the derivatives in the definition of consistent adjoint sensitivities
has been observed to provide a one to two order of magnitude improvement
in accuracy of sensitivities with respect to a finite difference approximation.
Significant differences have not been observed in optimisation behaviour or
in final designs. However, this cannot be guaranteed in general and it is
therefore best to ensure consistency, as also highlighted by the similar issue
of frozen turbulence [29].
Similar definitions and derivations are carried out for the thermal SUPG
stabilisation.
References
[1] M. P. Bendsøe, O. Sigmund, Topology Optimization: Theory, Methods
and Applications, Springer, 2003, ISBN: 3-540-42992-1.
[2] L. Yin, G. Ananthasuresh, A novel topology design scheme for
the multi-physics problems of electro-thermally actuated compliant micromechanisms, Sensors and Actuators 97-98 (2002) 599–609.
doi:10.1016/S0924-4247(01)00853-6.
[3] T. Bruns,
Topology optimization of convection-dominated,
steady-state heat transfer problems,
International Journal
of Heat and Mass Transfer 50 (15-16) (2007) 2859–2873.
doi:10.1016/j.ijheatmasstransfer.2007.01.039.
[4] A. Iga, S. Nishiwaki, K. Izui, M. Yoshimura, Topology optimization for thermal conductors considering design-dependent effects, including heat conduction and convection, International Journal of Heat and Mass Transfer 52 (11-12) (2009) 2721–2732.
doi:10.1016/j.ijheatmasstransfer.2008.12.013.
[5] E. Dede, S. N. Joshi, F. Zhou, Topology optimization, additive layer manufacturing, and experimental testing of an air-cooled
heat sink, ASME Journal of Mechanical Design (available online).
doi:10.1115/1.4030989.
[6] T. Borrvall, J. Petersson, Topology optimization of fluids in Stokes flow,
International Journal for Numerical Methods in Fluids 41 (1) (2003) 77–
107. doi:10.1002/fld.426.
21

[7] A. Gersborg-Hansen, O. Sigmund, R. Haber, Topology optimization of
channel flow problems, Structural Multidisciplinary Optimization 30 (3)
(2005) 181–192. doi:10.1007/s00158-004-0508-7.
[8] C. S. Andreasen, A. R. Gersborg, O. Sigmund, Topology optimization
of microfluidic mixers, International Journal for Numerical Methods in
Fluids 61 (5) (2009) 498–513. doi:10.1002/fld.1964.
[9] J. Alexandersen, N. Aage, C. S. Andreasen, O. Sigmund, Topology optimisation of natural convection problems, International Journal for Numerical Methods in Fluids 76 (10) (2014) 699–721.
doi:10.1002/fld.3954.
[10] P. Coffin, K. Maute, A level-set method for steady-state and transient
natural convection problems, submitted.
[11] E. Andreassen, A. Clausen, M. Schevenels, B. S. Lazarov, O. Sigmund, Efficient topology optimization in MATLAB using 88 lines of
code, Structural and Multidisciplinary Optimization 43 (2011) 1–16.
doi:10.1007/s00158-010-0594-7.
[12] O. Amir, N. Aage, B. S. Lazarov, On multigrid-CG for efficient topology
optimization, Structural and Multidisciplinary Optimization (2013) 1–
15doi:10.1007/s00158-013-1015-5.
[13] J. Alexandersen, B. S. Lazarov, Topology optimisation of manufacturable microstructural details without length scale separation using a spectral coarse basis preconditioner, Computer Methods in Applied Mechanics and Engineering 290 (2015) 156–182.
doi:10.1016/j.cma.2015.02.028.
[14] T. H. Nguyen, G. H. Paulino, J. Song, C. H. Le, Improving multiresolution topology optimization via multiple discretizations Tam, International Journal for Numerical Methods in Engineering 92 (June) (2012)
507–530. doi:10.1002/nme.
[15] T. Borrvall, J. Petersson, Large-scale topology optimization in 3D using
parallel computing, Computer Methods in Applied Mechanics and Engineering 190 (2001) 6201–6229. doi:10.1016/S0045-7825(01)00216-X.

22

[16] A. Evgrafov, C. J. Rupp, K. Maute, M. L. Dunn, Large-scale parallel topology optimization using a dual-primal substructuring solver,
Structural and Multidisciplinary Optimization 36 (2008) 329–345.
doi:10.1007/s00158-007-0190-7.
[17] N. Aage, B. Lazarov, Parallel framework for topology optimization using
the method of moving asymptotes, Structural and Multidisciplinary Optimization 47 (4) (2013) 493–505. doi:10.1007/s00158-012-0869-2.
[18] N. Aage, E. Andreassen, B. S. Lazarov, Topology optimization using
PETSc: An easy-to-use, fully parallel, open source topology optimization framework, Structural and Multidisciplinary Optimization 51 (3)
(2014) 565–572. doi:10.1007/s00158-014-1157-0.
[19] S. Balay, S. Abhyankar, M. F. Adams, J. Brown, P. Brune, K. Buschelman, V. Eijkhout, W. D. Gropp, D. Kaushik, M. G. Knepley, L. C.
McInnes, K. Rupp, B. F. Smith, H. Zhang, PETSc Users Manual, Tech.
Rep. ANL-95/11 - Revision 3.5, Argonne National Laboratory (2014).
URL http://www.mcs.anl.gov/petsc
[20] B. S. Lazarov, O. Sigmund, Filters in topology optimization based on
Helmholtz-type differential equations, International Journal for Numerical Methods in Engineering 86 (2011) 765–781. doi:10.1002/nme.3072.
[21] K. Svanberg, The method of moving asymptotes - a new method for
structural optimization, International Journal for Numerical Methods
in Engineering 24 (2) (1987) 359–373. doi:10.1002/nme.1620240207.
[22] M. Stolpe, K. Svanberg, On the trajectories of penalization methods
for topology optimization, Structural Multidisciplinary Optimization 21
(2001) 128–139.
[23] S. Rojas-Labanda, M. Stolpe, Benchmarking optimization solvers for
structural topology optimization, Structural and Multidisciplinary Optimization (available online). doi:10.1007/s00158-015-1250-z.
[24] P. S. Vassilevski, Multilevel Block Factorization Preconditioners,
Springer, 2007. doi:10.1007/978-0-387-71564-3.
[25] J.
Alexandersen,
O.
Sigmund,
N.
Aage,
Topology optimisation of passive coolers for light-emitting diode lamps,
23

in: 11th World Congress on Structural and Multidisciplinary Optimization, 2015. doi:10.13140/RG.2.1.3906.5446.
URL www.aeromech.usyd.edu.au/WCSMO2015/papers/1264_paper.pdf
[26] T. J. Hughes, L. P. France, M. Balestra, A new finite element formulation for computational fluid dynamics V - circumventing the
Babuska-Brezzi condition: a stable Petrov-Galerkin formulation of the
Stokes problem accomodating equal-order interpolations, Computer
Methods in Applied Mechanics and Engineering 59 (1) (1986) 85–99.
doi:10.1016/0045-7825(86)90025-3.
[27] T. E. Tezduyar, S. Mittal, S. Ray, R. Shih, Incompressible flow computations with stabilized bilinear and linear equalorder-interpolation velocity-pressure elements, Computer Methods
in Applied Mechanics and Engineering 95 (2) (1992) 221–242.
doi:10.1016/0045-7825(92)90141-6.
[28] A. N. Brooks, T. J. Hughes, Streamline Upwind/Petrov-Galerkin
formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations, Computer Methods in Applied Mechanics and Engineering 32 (1-3) (1982) 199–259.
doi:10.1016/0045-7825(82)90071-8.
[29] A. S. Zymaris, D. I. Papadimitriou, K. C. Giannakoglou, C. Othmer,
Continuous adjoint approach to the spalart-allmaras turbulence model
for incompressible flows, Computers and Fluids 38 (2009) 1528–1538.
doi:doi:10.1016/j.compfluid.2008.12.006.

24

