# Topology optimisation of heat sinks embedded with phase-change material for minimising temperature oscillations

Mark Bjerre Müller Christensen, Joe Alexandersen

► To cite this version:

Mark Bjerre Müller Christensen, Joe Alexandersen. Topology optimisation of heat sinks embedded with phase-change material for minimising temperature oscillations. 2023. hal-04185641

HAL Id: hal-04185641

https://hal.science/hal-04185641v1

Preprint submitted on 23 Aug 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Topology optimisation of heat sinks embedded with phase-change material for minimising temperature oscillations

Mark Bjerre Müller Christensen¹ and Joe Alexandersen¹*

¹Department of Mechanical and Electrical Engineering, University of Southern Denmark, Odense, DK-5230, Denmark.

*Corresponding author(s). E-mail(s): joal@sdu.dk;

## Abstract

This paper aims at providing insight into the optimal layout of heat sinks embedded with phase-change material. Topology optimisation is applied to reduce the amplitude of temperature oscillations in electronic components resulting from a transient cyclic thermal load. To do so, a Python code is used to perform topology optimisation using FEniCS for solving the physical problem, dolfin-adjoint for automatically computing adjoint sensitivities, and an open-source version of GCMMA for solving the optimisation problem. The phase-change material heat sinks are modelled using a transient thermal diffusion problem, where the phase change is modelled with the apparent heat capacity method. The design parametrisation and material interpolation is done using a simple analytical homogenisation approach. To reduce the temperature oscillation amplitude, the temporal variance of the spatially-averaged temperature at the electronic component was used as the objective functional for the optimisation. Initially, a time period of 20 seconds from a cold state is used as input for the topology optimisation, which resulted in optimised design that reduce the objective by reducing the initial transient evolution, rather than the amplitude of the temperature oscillation as intended. By changing the objective to only consider the temperature signal from a single cycle, after reaching a quasi-steady-state, a substantial reduction in the amplitude of the temperature oscillations is observed. Compared to the designs based on 20 seconds, the amplitude was reduced with up to 78%.

**Keywords:** topology optimisation, phase-change material, temperature oscillation, heat sink, FEniCS, dolfin-adjoint
# 1 Introduction

## 1.1 Motivation

As electronics become a larger part of our everyday lives, the importance of reliability, cost and life-time of electronic components increases. A major cause for failure in electronic components is the mechanical stresses induced by the mismatch in thermal coefficient of expansion (TCE) across materials in an electronic component and an increase in temperature. In cases where the electronic component has a transient cyclic heat production, these mechanical stresses can further lead to fatigue failures (Depiver et al, 2021) if not managed correctly through adequate cooling. In practice the temperature of electronic components with high heat production are kept cool with heat sinks which are placed on top of the component and dissipate the heat produced by the component to the surrounding air, an “example” can be seen in Figure 1.

<figure><img src="image_2.png" /><figcaption>Fig. 1 Three-dimensional sketch of generic electronics cooling problem.</figcaption></figure>

This paper focuses on heat sinks with embedded phase-change material (PCM) (henceforth referred to as “PCM heat sink”) for controlling the amplitude of the temperature oscillations caused by a transient thermal load. The idea is to use the thermal storage available through the latent heat of fusion in the PCM as a thermal buffer to smoothen the temperature oscillations caused by a cyclic thermal loading.

A major issue with PCM from a thermal design aspect is that most common PCMs have a low thermal conductivity making it hard to get the heat from the electronic component into the PCM. Therefore, PCM is usually paired with a highly thermally conductive material (HCM) in the form of fins or foams which ensure the distribution of heat to the PCM (Sahoo et al, 2016). In the literature numerous simple fin designs in just as many operational modes have been investigated through extensive parametric investigations in order to determine best performing design (Sahoo et al, 2016; Kalbasi et al, 2019). Most of these studies are based on a steady thermal loading, which makes the physical problem easier to compute. Resulting in a poor understating of the best use of PCM in an PCM heat sink under a transient thermal loading is therefore. This paper aims to deepen the understanding of the optimal design of PCM
heat sinks with respect to reducing the amplitudes in thermal oscillation induced by
a cyclic thermal loading. However, instead of continuing the tradition of parametric
investigations which are costly to perform, limited by the parametrisation of the HCM
design and specific in their applications, this paper seeks to deepen the understanding
with the help of gradient-based topology optimisation (TO).

Currently, there is no commercially available software capable of performing TO
with respect to reducing the amplitudes in thermal oscillation induced by a cyclic
thermal loading. Therefore, a Python code capable of this had to be created for this
paper. The Python code uses the FEniCS package in conjunction with the dolphin-
adjoint module to perform the TO.

1.2 Literature

Several numerical and experimental studies have been conducted to improve the performance of two-dimensional plate fin-based PCM heat sinks investigating the effects of the number of fins, heat sink height, and fin thickness (Sahoo et al, 2016; Hosseinizadeh et al, 2011; Kalbasi et al, 2019). However, the studies have focused on optimising the performance with regard to the temperature or time until critical temperature is reached. The studies also suffer from a restriction to simple designs as this makes the construction of parametric studies possible. As an alternative to this restrictive approach this paper seeks to employ topology optimisation to optimise the two-dimensional plate fin-based PCM heat sinks.

For a general overview of TO for heat transfer applications, the reader is suggested
to look at the available review papers (Dbouk, 2017; Alexandersen and Andreasen,
2020; Fawaz et al, 2022) on the subject. The usage of TO for optimising the layout
PCM and HCM for heat sinks is sparsely represented in literature, whereas the use of
PCM and HCM for thermal energy storage is more common. Varanasi and Anantha-
suresh (2006) presented the first use of TO for PCM materials using a density-based
TO to minimise the difference between the initial and final temperature over a selected
part of the heat sink boundary. They modelled the PCM heat sink as a transient
thermal diffusion problem involving phase change using COMSOL Multiphysics and
MATLAB. They used an enthalpy method for modeling the phase change, a solid
isotropic material with penalisation (SIMP) method for interpolating the conductiv-
ity, and an optimality criterion method to update the design. The results showed a
35% reduction in the objective functional when the optimised design was compared to
a conventional design. However, the designs suffer from poor resolution due to the lim-
ited computational power available at the time. Ho et al (2021) optimised PCM heat
sinks for maximum cooling in lieu of convective cooling, thus, storing the energy from
the heat source in the PCM. They applied a transient thermal diffusion problem with
the modified heat capacity method, comparing performance between designs optimised
using steady and unsteady models both numerically and experimentally. Subsequently,
they added the effect of natural convection in the melted PCM (See et al, 2022) and
Iradukunda et al (2020) does indeed treat the subject of PCM-integrated heat sinks
using topology optimisation, but use a very simple steady-state heat transfer model for
the actual optimisation. Thus, it cannot be classified as true TO of PCM heat sinks.
However, they did show improved thermal buffering from topology-optimised designs,
which is the objective pursued in the present work. Similarly, Bianco et al (2023) opti-
mised using a steady-state model with heat generation in place of an actual phase
change model, but investigated the application of density-based TO to minimise the
thermal resistance under constant wall temperature. Despite using a simplified model,
the optimised design showed an up to 2 times improvement of the heat flux compared
to a conventional design.

As mentioned, the use of TO for thermal energy storage systems is more common
and still relevant to summarise, since all of the work actually makes use of transient
heat transfer analysis using models that includes the latent heat of fusion. Pizzolato
et al (2017b) was the first to use TO to improve heat transfer performance in a
latent heat thermal energy storage system (LHTES). In the initial study, only thermal
diffusion was considered, but this was later extended to include the effects of natural
convection (Pizzolato et al, 2017a, 2020). They found that optimised designs based on
pure thermal diffusion had higher heat transfer rates during the initial melting phase
compared to optimised designs based on models including convective heat transfer in
the liquid PCM. While designs based on models including convection had a higher
charge rate, if the storage unit had stored more than 80.2% of its energy storage
capacity. There exists several further examples on TO for PCM and HCM of LHTES
(Zhao et al, 2020; Yao et al, 2021; Tian et al, 2021; Laasri et al, 2022; Peremans et al,
2023; Zhang et al, 2023), the details of which will not be discussed herein.

FEniCS (Alnæs et al, 2015; Logg et al, 2012) is a widely-used collection of open-source components for automating the solution of differential equations using the finite element method. In the context of topology optimisation, it has been used by a number of authors and the following only presents a brief selection. Laurain (2018) presented a compact educational structural TO code written using FEniCS, to perform compliance minimisation using a level set method. Qian (2017) used FEniCS to perform TO to improve the printability of the optimised design on a three-dimensional printer by adding undercut and overhang angle control. Mezzadri et al (2018) used FEniCS to perform TO of self-supporting support structures for additive manufacturing. Yan et al (2022) presented a general-purpose TO platform for multidisciplinary problems with FEniCS as the multiphysics solver. Jauregui et al (2023) recently presented a reusable modular TO software, capable of using FEniCS as the simulation tool.

1.3 Contributions

This paper aims to build upon the current understanding of PCM heat sink with help of
gradient-based topology optimisation with a specific focus on reducing the temperature
oscillation in the electronic components resulting from a transient cyclic thermal load
from the electronic components. Based on a simplified model of a PCM heat sink, a
topology optimisation code using FEniCS is applied to optimise the layout of PCM
and HCM in order to reduce the amplitude of temperature oscillations. The novelty
in this research lies in the relatively complex objective functional for the optimisation,
which gives insight in the complex problem that is reducing the amplitude of the
temperature oscillation, the addition of cooling when modeling the thermal problem,
and use of the open-source program FEniCS to do this.
This paper presents the major methods and assumptions used for the simulation and optimisation of the heat sink, and the resulting designs from the optimisation under different conditions. The optimised designs are investigated and based on these investigations it is expected that some patterns will emerge, that can help further the understanding of PCM heat sinks.

## 1.4 Paper layout

The rest of the paper is organised as follows: Section 2 presents the methods and assumptions used to set up a physical model of the PCM heat sink; Section 3 presents the optimisation problem and methods used to performing TO; Section 4 describes how the methods are implemented to produce the results and presents the verification of the implementation of the methods; Section 5 presents the optimised designs of PCM heat sinks for different use cases; Section 6 discusses the result and assumptions behind the physical model and TO and their implication; Section 7 presents the main conclusions of the paper.

# 2 Physical model

As this paper aims to optimise PCM heat sinks with respect to the temperature of electronic components, the thermal problem is the most important to solve. The presence of PCM material adds a lot of complexity to the thermal problem, the melting and solidification requires a transient model, the latent heat of fusion introduces non-linearities and the liquid PCM causes natural convection heat transfer due to buoyancy effects which requires a fluid simulation to model fully.

The TO requires a lot of evaluations of the physics, Therefore, it is desired to keep the physical model as simple as possible to keep the computational time to a manageable level. Therefore, the effects of natural convection will be neglected in this study. This choice is expected to result in some errors in the physical model as a review of modeling PCM in LHTES (Jegadheeswaran and Pohekar, 2009) found that natural convection had a major effect on the melting of PCM. The effect of these errors will have to be investigated in future works. With this delimitation, it is assumed that the PCM heat sinks can be modelled as a transient thermal diffusion problem with two materials and phase change.

For this paper, the PCM heat sink is modelled as a simplified two-dimensional model consisting of a unit square with a heat flux applied at the boundary $\partial\Omega_{HS}$ simulating the heat produced in the electronic component and a cooling flux applied at the boundary $\partial\Omega_C$. The remaining boundaries are considered adiabatic. A sketch of the simplified two-dimensional model of the PCM heat sink can be seen in Figure 2.

## 2.1 The transient thermal diffusion problem

The PCM heat sink is modelled using the following thermal diffusion problem:

$$ \rho_m c_p \frac{\partial T}{\partial t} - k_T \frac{\partial^2 T}{\partial x_i \partial x_i} = 0 \qquad (1) $$
<figure><img src="image_3.png" /><figcaption>Fig. 2 two-dimensional sketch of the simplified PCM heat sink used to define the design domain for the topology optimisation.</figcaption></figure>

where: $\rho_m c_p \frac{\partial T}{\partial t}$ is the thermal storage term; $-k_T \frac{\partial^2 T}{\partial x_i^2}$ is the diffusion term; T is the continuous temperature field; $x_i$ is a vector containing all the spatial dimension; $k_T$ is the thermal conductivity; $\rho_m$ is the mass density; $c_p$ is the specific heat; t is time.

To get FEniCS to solve the thermal diffusion problem, the weak form has to be derived. Using a backward difference scheme to handle the time derivative, the weak form is derived to the following after integration-by-parts:

$$R = \int_{\Omega} \rho_m c_p v \frac{T^{(n)} - T^{(n-1)}}{\Delta t} dA + \int_{\Omega} k_T \frac{\partial v}{\partial x_i} \frac{\partial T^{(n)}}{\partial x_i} dA + \int_{\partial \Omega} v \left(-k_T \frac{\partial T^{(n)}}{\partial n}\right) ds = 0 \quad (2)$$

where: $\Delta t$ is the size of the time step used for the backward difference scheme; $T^{(n)}$ is the continuous temperature field at the current time step n; $T^{(n-1)}$ is the continuous temperature field at the previous time step; v is test function; $\partial \Omega$ is the boundary of the design domain; $-k_T \frac{\partial T}{\partial n}$ is the normal flux over the boundary. Note that the weak form is provided in the semi-discrete form, being continuous in the spatial dimension and discrete over time, in order to be consistent with the notation used in the code.

## 2.2 Material interpolation

As the PCM heat sink consists of both PCM and HCM, the physical model has to be able to handle two materials. To make the subsequent implementation of TO easier the two materials were implemented with an interpolation scheme depending on the material density variable field, $\rho$, representing the percentage of HCM, such that $\rho = 0$
is PCM, $\rho = 1$ is HCM, and $\rho \in ]0; 1[$ is a mixture. The material properties important for the physical model are $\rho_m c_p$ and $k_T$.

As seen in Section 1.2, the SIMP method is a common method for material interpolation when performing TO. It introduces a power-law interpolation that artificially decreases the material properties per volume for intermediate densities which usually incentives TO to create optimised designs without intermediate material (henceforth referred to as “discrete designs”) to make the best use of the material. However, during the initial testing it was found that the optimised designs based on the SIMP method had a tendency to contain significant amounts of intermediate material making them unphysical. Forcing the design to become more discrete by further penalising intermediate materials or using a projection scheme led to a drastic decrease in performance. Therefore, it was necessary to implement an interpolation scheme that has a physical meaning for intermediate density variables. To do this, a simple analytical homogenisation was chosen.

Homogenisation is a method that computes effective macroscopic material properties of a microstructure consisting of a locally periodically repeating pattern of a base cell, see Figure 3.

<figure><img src="image_4.png" /><figcaption>Fig. 3 Homogenisation is based on the assumption that intermediate material can be seen as a locally regular periodically repetition of a base cell.</figcaption></figure>

A microstructural design variable is then connected to the macroscopic effective properties to allow for optimisation. The chosen microstructure is a square base cell with a square of PCM in the centre with a frame of HCM around it. The macroscopic design variable is defined as the relative material density of HCM, $\rho = \frac{|\Omega_{HCM}|}{|\Omega_{HCM} \cup \Omega_{PCM}|}$, which in turn is coupled to the microscopic design variable, defined as the width of the HCM frame dictated by the volume of HCM. This definition is further described in A, where the simple analytical homogenisation applied herein is also described.

The thermal conductivity $k_T$ is computed as the total thermal conductivity based on an one-dimensional steady-state heat conduction analysis of the base cell:

$$k_T(\rho) = \frac{1}{\frac{1-\sqrt{1-\rho}}{k_{THCM}} + \frac{\sqrt{1-\rho}}{(1-\sqrt{1-\rho})k_{THCM}+(\sqrt{1-\rho})k_{TPCM}}} \quad (3)$$
The volumetric heat capacity, defined as the mass density multiplied by the specific heat capacity, $\rho_m c_p$, is computed as a volumetric average:

$$ \rho_m c_p(\rho) = (\rho)\rho_m HCM c_p HCM + (1 - \rho)\rho_m PCM c_p PCM \quad (4) $$

## 2.3 Latent heat of fusion

To model the latent heat of fusion, the apparent heat capacity method is used. The method introduces a temperature-dependent apparent heat capacity, $c_{pPCM} \Rightarrow c_{pPCM}(T^{(n)})$, where the value of the heat capacity is increased in the phase change temperature range to incorporate the latent heat of fusion into the thermal storage term. Assuming the change in heat capacity between the liquid and solid state of PCM is negligible, $c_{pPCM}(T^{(n)})$ can be modelled as the following piece-wise function:

$$ c_{pPCM}(T^{(n)}) = \begin{cases} c_{pPCM} & \text{if } T^{(n)} < T_{melt} - \frac{\Delta T_{melt}}{2} \\ c_{pPCM} + \frac{L_{heat}}{\Delta T_{melt}} & \text{if } T_{melt} - \frac{\Delta T_{melt}}{2} \le T^{(n)} \le T_{melt} + \frac{\Delta T_{melt}}{2} \\ c_{pPCM} & \text{if } T^{(n)} < T_{melt} + \frac{\Delta T_{melt}}{2} \end{cases} \quad (5) $$

where: $T_{melt}$ is the melting temperature; $\Delta T_{melt}$ is the phase change temperature range; $L_{heat}$ is the latent heat of fusion. Deriving the sensitivities in later steps requires a continuous function, therefore, the piece-wise function is approximated with smooth Heaviside step functions:

$$ c_{pPCM}(T^{(n)}) = c_{pPCM} + \frac{L_{heat}}{\Delta T_{melt}} \left( \frac{1}{1 + e^{-2k_H \left( T^{(n)} - \frac{\Delta T_{melt}}{2} \right)}} - \frac{1}{1 + e^{-2k_H \left( T^{(n)} + \frac{\Delta T_{melt}}{2} \right)}} \right) \quad (6) $$

where: $k_H$ is a factor that corresponds to the sharpness of the transition. For clarity of what the Heaviside function does, the temperature dependent apparent heat capacity, $c_{pPCM}(T^{(n)})$ is plotted using the piece-wise function and the smooth Heaviside step function, see Figure 4.

## 2.4 Boundary conditions

The fluxes used to model the heat source and cooling are modelled as Neumann boundary conditions. For the heat source, a uniformly distributed time-dependent heat rate is used to simulate the transient thermal loading from the electrical component:

$$ -k_T \left. \frac{\partial T}{\partial n} \right|_{\partial \Omega_{HS}} = -q_{HS}(t) = -\frac{P_{elec}}{A} (1 + \sin(2\pi\omega t)) \quad (7) $$

where: $-k_T \left. \frac{\partial T}{\partial n} \right|_{\partial \Omega_{HS}}$ is the heat flux normal to the boundary $\partial \Omega_{HS}$; $P_{elec}$ is the average heat transfer rate produced by the electronic component; $A$ is the area of the electronic component; $\omega$ is the oscillation frequency of the heat rate produced by the electronic component.
<figure><img src="image_5.png" /><figcaption>Fig. 4 The temperature dependent apparent heat capacity, $c_{pPCM}(T^{(n)})$ plotted at different temperatures with the piece-wise function and the smooth Heaviside step function.</figcaption></figure>

The cooling is modelled as a convection cooling rate using Newton’s law of cooling:

$$-k_T \frac{\partial T}{\partial n} \Big|_{\partial \Omega_C} = q_C(T_s) = h_{conv}(T_s - T_\infty) \quad (8)$$

where: $-k_T \frac{\partial T}{\partial n}|_{\partial \Omega_C}$ is the heat flux normal to the boundary $\partial \Omega_C$; $h_{conv}$ is the convection heat transfer coefficient; $T_s$ is the surface temperature; $T_\infty$ is the temperature of the surroundings. Note that the term for the heat source is negative while the term for the cooling is positive. This is due to the fact that FEniCS considers the normal to be positive in the direction away from the body.

## 2.5 Weak form and approximations

Inserting Equation (3), (4), (6), (7), and (8) into (2) results in the following weak form, which can be used to model the physics behind the PCM heat sink:

$$R = \int_{\Omega} \rho_m c_p (\rho, T^{(n)}) v \frac{T^{(n)} - T^{(n-1)}}{\Delta t} dA + \int_{\Omega} k_T(\rho) \frac{\partial v}{\partial x_i} \frac{\partial T^{(n)}}{\partial x_i} dA - \int_{\partial \Omega_{HS}} v q_{HS}(t) ds + \int_{\partial \Omega_C} v h_{conv}(T^{(n)} - T_\infty) ds = 0 \quad (9)$$

This weak form is a non-linear function due to the implementation of the latent heat of fusion, where the heat capacity is dependent on the current temperature $T^{(n)}$, which means it requires a non-linear solver to solve the model. A non-linear solver takes extra
iterative solver steps to deal with the non-linearity, which increases the computational cost and time of the solving the model.

To keep the computational time as short as possible, the non-linear model is approximated by a time-lagging model by making the $\rho_m c_p$ dependent on the temperature from the previous time step $T^{(n-1)}$ rather than the current time step $T^{(n)}$. This removes the need for the extra iterative solver steps from the non-linear solver at the cost of introducing an error into the model. The final weak form used to modeling the PCM heat sink in this paper thereby becomes:

$$ R = \int_{\Omega} \rho_m c_p(\rho, T^{(n-1)}) v \frac{T^{(n)} - T^{(n-1)}}{\Delta t} dA + \int_{\Omega} k_T(\rho) \frac{\partial v}{\partial x_i} \frac{\partial T^{(n)}}{\partial x_i} dA - \int_{\partial \Omega_{HS}} v q_{HS}(t) ds + \int_{\partial \Omega_{C}} v h_{conv}(T^{(n)} - T_{\infty}) ds = 0 \quad (10) $$

# 3 Optimisation problem

## 3.1 Topology optimisation

TO aims to find the most optimal distribution of two materials to minimise a chosen objective functional. In density-based topology optimisation, the material distribution is described with a spatially-varying material density variable field, $\rho$. To avoid simple solutions and keep different cases comparable, a volume constraint is added to the optimisation problem, restricting the use of HCM. The optimisation problem used in this paper can be seen in the following equation:

$$ \begin{array}{ll} \min_{\rho} & f_0(T(\rho), \rho) \\ \text{s.t.} & h_0(\rho) = \frac{\int_{\Omega} \rho(x_i) dA}{\int_{\Omega} \Phi dA} - 1 \\ & 0 \le \rho(x_i) \le 1 \end{array} \quad (11) $$

where: $f_0$ is the objective functional; $\rho$ is the continuous material density variable field; $T(\rho)$ is the temperature field that satisfies the physical model; $h_0(\rho)$ is the volume constraint; $\Omega$ is the continuous design domain; $\Phi$ is the maximum allowable volume fraction of HCM.

## 3.2 Objective functional

The objective of the optimisation is to minimise the amplitude of the temperature oscillation at the electronic component, that is introduced by the cyclic transient heat rate from the electronic component. This is quantified with the $\phi$, which is the temporal average variance of the spatial average temperature at $\partial \Omega_{HS}$, as this is where the PCM heat sink is connected to the electronic component (henceforth refereed to as “the variance”). In order to get the TO to reduce the amplitude of the temperature oscillation, the variance, $\phi$, is chosen as the objective functional $f_0$ for the TO. Thereby
the objective function is defined as:

$$f_0 = \phi = \frac{1}{N_t} \sum_{n=1}^{N_t} \left( T_{elec}^{(n)} - \bar{T}_{elec} \right)^2 \qquad (12)$$

where: $N_t$ is the final number of time steps; $n$ is the time step; $T_{elec}^{(n)}$ is the spatial average temperature over the heat source boundary, $\partial\Omega_{HS}$, at each time step (henceforth refereed to as “the temperature at the heat source”); $\bar{T}_{elec}$ is the temporal mean of $T_{elec}^{(n)}$. The temperature at the heat source $T_{elec}^{(n)}$ is defined as:

$$T_{elec}^{(n)} = \frac{\int_{\partial\Omega_{HS}} T^{(n)} ds}{\int_{\partial\Omega_{HS}} L_z ds} \qquad (13)$$

where: $L_z$ is the out-of-plane length of the heat sink; $T^{(n)}$ in the temperature field at time step $n$. The temporal mean of $T_{elec}^{(n)}$ is computed with:

$$\bar{T}_{elec} = \frac{1}{N_t} \sum_{n=1}^{N_t} T_{elec}^{(n)} \qquad (14)$$

## 3.3 Filtering

A direct implementation of the density-based TO method to a heat transfer problem can result in designs that are mesh-dependent and contains unphysical checkerboard patterns (Sigmund and Petersson, 1998). To alleviate these problems different filtering strategies have been proposed (Lazarov et al, 2016). For this paper, a filtering method based on a partial differential equation (PDE) (Lazarov and Sigmund, 2011) is used, as it takes advantage of the ease of implementing and solving PDEs in FEniCS. The filtered density field $\tilde{\rho}$ is found by solving the following PDE (Lazarov and Sigmund, 2011):

$$-r^2 \frac{\partial^2 \tilde{\rho}}{\partial x_i \partial x_i} + \tilde{\rho} = \rho \qquad (15)$$

where: $r$ is a filter parameter linked to the amount of smoothing applied by the filter. The PDE is combined with a homogeneous Neumann boundary condition, to ensure that volume is conserved throughout the filtering process. For the optimisation, the filtered density variable $\tilde{\rho}$ is used as the design input for the physical model of the PCM heat sink.

## 3.4 The adjoint method

In order to perform gradient-based topology optimisation, the sensitivities of the constraint and objective functionals with respect to the material density variable field, $\rho$, have to be computed. As the topology optimisation is characterised by having a large number of design variables and a small number of constraints, the adjoint method is particularly efficient for computing the sensitivities as it only requires one additional problem per functional to compute it.
In short the adjoint method computes the sensitivities using the following equation:

$$ \frac{df_j}{d\rho} = \frac{\partial f_j}{\partial \rho} - \lambda^T \frac{\partial R}{\partial \rho} \qquad (16) $$

where: $f_j$ are functionals of interest; $\rho$ is a vector containing the density variables; $R$ is a vector containing the residuals; $\lambda$ is a vector containing the adjoint variables. The adjoint variables are computed by solving the adjoint system of equations:

$$ \left( \frac{\partial R}{\partial u} \right)^T \lambda = \left( \frac{\partial f_j}{\partial u} \right)^T \qquad (17) $$

where: $u$ is a vector containing the state variables.

For time-dependent systems, the adjoint system also becomes time-dependent and, thus, another full time series solution must be solved. Furthermore, the full time history of the state variables and the system Jacobian must be saved (or recomputed) for the adjoint system. Although this is “just” a doubling of the computational cost, this is much more significant for time-dependent problems, since their computational cost is so high already from the time stepping procedure. For further details on adjoint sensitivity analysis of time-dependent systems, please see the literature on the subject (Michaleris et al, 1994; Choi and Kim, 2005). To circumvent the cumbersome work of deriving and implementing the full time-dependent adjoint system, the automatic adjoint capabilities of the dolfin-adjoint module (Mitusch et al, 2019) is used in combination with FEniCS.

# 4 Implementation

## 4.1 Packages

In order to implement the methods described in the previous section a Python code was created. Python was chosen due to its open source nature and vast number of modules that can automate many of the computational steps required for TO. For this paper, the FEniCS (legacy FEniCS version 2019.1.0) module is used for solving the physical model, the dolfin-adjoint (dolfin-adjoint version 2019.1) module is used to compute the adjoint sensitivities. The optimisation problem is solved using the Globally Convergent Method of Moving Asymptotes (GCMMA) (Svanberg, 2002) and was implemented into the Python code through the code provided by Deetman (2020).

## 4.2 The code

In order to perform TO, the process is broken down into the following computational steps:

1. Defining the initial design as a uniform distribution of material.

2. Solve the physical problem with FEniCS.

3. Compute the objective and constraint sensitivities with the adjoint method.

4. Update the design using GCMMA.
5. Check for convergence. If the convergence criteria are not satisfied, update the geometry and go back to step 2, otherwise stop the optimisation loop.

6. Save the final optimised design.

First the physical model is defined. The material properties, model parameters, geometry, discretisation schemes and boundary conditions are defined using FEniCS. The geometry is meshed with a mapped mesh using triangular elements in a cross pattern, see Figure 5. The mesh uses triangular elements as these are the only ones compatible with the FEniCS version used for this paper. The temperature fields, $T^{(n)}$ and $T^{(n-1)}$ are approximated with continuous first-order elements while the material density variable field, $\rho$, is approximated using piece-wise constant discontinuous elements.

<figure><img src="image_6.png" /><figcaption>Fig. 5 Plot of the mapped mesh using triangular element used for this paper.</figcaption></figure>

The solver for the physical problem is set up as a custom function, `forward()`, that takes in the filtered material design variable field, $\tilde{\rho}$, as an input. In order to solve the physical model, an iterative time loop is set up, solving the physical model for each time step until the final time step is reached, $n = N_t$. The temperature at the heat source $T_{elec}^{(n)}$ is stored for each time step and used to compute the variance $\phi$.

The computation of the sensitivities is automated with `dolfin-adjoint`. For that the objective and constraint functionals have to be defined in a `ReducedFunctional` class, which stores all the operations performed with FEniCS and defines the control parameter. With this the dependencies are stored and both the functionals and the sensitivities of the functionals can be computed based on the current material design variable field, $\rho$.

In order to update the design, `GCMMA` requires the current material design variable field, $\rho$, the resulting objective and constraint functionals together with the sensitivities of the functionals, which pair up nicely with the capabilities of the `ReducedFunctional` class from `dolfin-adjoint`. The `GCMMA` is set to run a maximum of
two inner iterations to keep the computational time to a reasonable level since transient
simulations are very costly. To check whether the optimisation has converged the code
uses the change in the objective functional, $f_0$, and the measure of non-discreetness,
$M_{nd}$ (described in detail in Section 5.4.1), over the optimisation iteration. The opti-
misation is considered converged if the absolute changes in both $f_0$ and $M_{nd}$ is below
$10^{-3}$ for three consecutive optimisation iterations. If the convergence criteria is not
met, the optimisation is allowed to run a maximum of 300 optimisation iterations.

The TO is set up as a custom function, optimisation(), which takes in the initial
material density variable field, ρ, the initial filtered material density variable field
ρ̃, Reduced functionals of the objective and constraint functionals, lists of historical
values for convergence plots, and the α value. The operations of optimisation() are
enveloped in with a "with pyadjoint.stop_annotating() as _:" command which stops the
dolfin-adjoint from storing the operations, which can otherwise lead to a large memory
consumption. A minimal working example of the code is available on Github - see the
"Replication of results" section at the end of this paper.

4.3 Verification

4.3.1 Comparison with COMSOL

To verify the implementation of the physical model in FEniCS, a 20 second simulation
of a benchmark design was conducted and compared to results from the same model
simulated with the commercial finite element program COMSOL Multiphysics (COM-
SOL, 2023) using a non-linear solver. Figure 6 shows the design and results from the
comparison. It was found that there is a good agreement between the results with a
maximum relative error of 1.451%, whereby the implementation of physical model in
FEniCS is considered verified.

<figure><img src="image_7.png" /><figcaption>**Fig. 6** Verification of the implementation of the physical problem with FEniCS by comparison with COMSOL Multiphysics. The results are based on a 4 × 100 × 100 mesh using a time step size of 0.02 s. Left: Design used as benchmark for the comparison. Middle: The temperature at the heat source, *T*<sub>elec</sub>, over time. Right: The relative error of the temperature at the heat source, *T*<sub>elec</sub>, over time.</figcaption></figure>

The accuracy of the physical model largely depends on the choice of time step size,
due to the highly non-linear smooth Heaviside step function used for the apparent
heat capacity method and by the approximation of the non-linear physical model with
a time-lagging one. However, the time step size also has a great impact on the compu-
tational time. Therefore, the effects of time step size on the error and computational
time is investigated by comparing the results from FEniCS with varying time step
sizes with the results from COMSOL using a much smaller time step size and a non-
linear solver. For the FEniCS simulation both the time-lagging and non-linear solver
are used to also show the effects of using time-lagging physical model. The results of
this comparison is the maximum relative error in T<sub>elec</sub> plotted against the time step
size and computational time, see Figure 7.

<figure><img src="image_8.png" /><figcaption>**Fig. 7** Investigation of maximum relative error in the temperature at the heat source $T_{Elec}$ for different time step sizes when comparing FEniCS simulations with a COMSOL simulation using a non-linear solver. The results are based on a COMSOL simulation using a time step size of $2 \cdot 10^{-3}$. Left: The error plotted against time step size. Right: The error plotted against computational time. The computational time is based on the average simulation time when the physical model is run on i7-1165G7 @ 2.80GHz with 4 GB.</figcaption></figure>

From Figure 7 it can be seen that the maximum relative error convergences with
a decrease in time step size, as there is an linear relationship between the time step
size and maximum relative error in the log-log plot. It can be seen that the maximum
relative error is similar when using a time-lagging and non-linear physical model when
using the same time step size, with the non-linear physical model a slightly higher
error. This is unexpected as the non-linear solver is expected to more precise due
to the additional iterative solver steps. However looking at the convergence rate it
can be seen that the non-liner solver has a higher convergence rate, so the error is
expected to be lower as the time step size is increased. Looking at the computational
time for solving the physical model, the non-linear model is much more computational
expensive compared to the time-lagging model. In order to limit the computational
time, an maximum relative error of 5% is deemed acceptable whereby the time-lagging
physical model with a time step size of 0.04 s is chosen for the simulation of the
physical model when evaluating the designs for the TO.
4.3.2 Finite difference check

To verify the accuracy of the adjoint sensitivities, a finite difference check is performed using a central difference scheme. Figure 8 shows the relative error between the finite difference approximation and the adjoint sensitivity of the objective functional at element number 2400, $\frac{\partial f_0}{\partial \rho_{2400}}$, at different perturbation sizes. From Figure 8, it can be seen that the convergence rate for $\varepsilon \ge 10^{-3}$ agrees with the expected convergence rate of a central difference scheme of $\approx O(\varepsilon^2)$ (Haftka and Gürdal, 1992). For perturbation size $\varepsilon < 10^{-3}$, the rounding errors become significant which increases the relative error. This relation between perturbation size and error is well known in finite difference approximations and is called “step-size dilemma” (Haftka and Gürdal, 1992). Conducting finite difference checks for additional element numbers showed that the tendencies seen in $\frac{\partial f_0}{\partial \rho_{2400}}$ are representative for all elements. As the convergence rate for $\varepsilon \ge 10^{-3}$ matches the expected convergence rate for the central difference scheme and the relative error is small, the implementation of the adjoint method is considered verified.

<figure><img src="image_9.png" /><figcaption>Fig. 8 Finite difference check, showing the relative error between the adjoint sensitivity computed with dolfin-adjoint and central difference scheme at different perturbation sizes at different stages of the optimisation. The expected convergence rate of a central difference scheme $O(\varepsilon^2)$ is plotted for reference. The optimisation is based on a 4 × 40 × 40 mesh.</figcaption></figure>

4.4 Data

This section presents the data used for the production of the results, unless stated otherwise. The dimensions used for the model are presented in Table 1, the material
properties are presented in Table 2, and the model parameters are presented in Table 3.

<table><caption>**Table 1 Geometry and discretisation**</caption><thead><tr><th>Description</th><th>Symbol</th><th>Value</th><th>Units</th></tr></thead><tbody><tr><td>Domain size</td><td>L<sub>x</sub> × L<sub>y</sub> × L<sub>z</sub></td><td>1 × 1 × 1</td><td>[m]</td></tr><tr><td>Mesh size</td><td>n<sub>epsq</sub> × n<sub>x</sub> × n<sub>y</sub></td><td>4 × 100 × 100</td><td>[-]</td></tr><tr><td>Final time</td><td>t<sub>fin</sub></td><td>20</td><td>[s]</td></tr><tr><td>Number of time steps</td><td>N<sub>t</sub></td><td>500</td><td>[-]</td></tr><tr><td>Time step size</td><td>Δt</td><td>0.04</td><td>[s]</td></tr></tbody></table>

<table><caption>**Table 2 Material properties**</caption><thead><tr><th>Description</th><th>Symbol</th><th>Value</th><th>Units</th></tr></thead><tbody><tr><td>Thermal conductivity HCM</td><td>k<sub>THCM</sub></td><td>10</td><td>[W/(m·K)]</td></tr><tr><td>Density HCM</td><td>ρ<sub>mHCM</sub></td><td>1</td><td>[kg/m<sup>3</sup>]</td></tr><tr><td>Heat capacity HCM</td><td>c<sub>pHCM</sub></td><td>1</td><td>[J/(kg·K)]</td></tr><tr><td>Thermal conductivity PCM</td><td>k<sub>TPCM</sub></td><td>0.01</td><td>[W/(m·K)]</td></tr><tr><td>Density PCM</td><td>ρ<sub>mPCM</sub></td><td>1</td><td>[kg/m<sup>3</sup>]</td></tr><tr><td>Heat capacity PCM</td><td>c<sub>pPCM</sub></td><td>1</td><td>[J/(kg·K)]</td></tr><tr><td>Melting temperature</td><td>T<sub>melt</sub></td><td>0.5</td><td>[K]</td></tr><tr><td>Melting temperature range</td><td>ΔT<sub>melt</sub></td><td>0.5</td><td>[K]</td></tr><tr><td>Latent heat of fusion</td><td>L<sub>heat</sub></td><td>10</td><td>[J/(kg)]</td></tr></tbody></table>

<table><caption>**Table 3 Model parameters**</caption><thead><tr><th>Description</th><th>Symbol</th><th>Value</th><th>Units</th></tr></thead><tbody><tr><td>Heat transfer coefficient</td><td>h<sub>conv</sub></td><td>5</td><td>[W/(m<sup>2</sup>·K)]</td></tr><tr><td>Initial temperature</td><td>T<sub>ini</sub></td><td>0</td><td>[K]</td></tr><tr><td>Surrounding temperature</td><td>T<sub>∞</sub></td><td>0</td><td>[K]</td></tr><tr><td>Average heat rate from electronic component</td><td>P<sub>elec</sub></td><td>1</td><td>[W]</td></tr><tr><td>Heat rate oscillation frequency</td><td>ω</td><td>1</td><td>[Hz]</td></tr><tr><td>Maximum volume fraction of HCM</td><td>Φ</td><td>0.3</td><td>[-]</td></tr><tr><td>Steepness factor for smooth Heaviside step function</td><td>k<sub>H</sub></td><td>25</td><td>[-]</td></tr><tr><td>Filter parameter</td><td>r</td><td>0.01</td><td>[-]</td></tr></tbody></table>
# 5 Results

This section presents and discusses the results obtained from running the TO under different conditions. All designs are visualised with grey-scale plots of the filtered material density variable field, $\tilde{\rho}$, where white indicates PCM, black indicates HCM, and grey indicates a mix of the two materials. The colourbar for all the designs can be seen in Figure 9.

<figure><img src="image_10.png" /><figcaption>Fig. 9 Colourbar for all the designs presented in this paper.</figcaption></figure>

To test the effect of including phase change, three different cases are considered for the TO. The three cases considered are: Case 1 - Without phase change and with $c_{pPCM} = 1$; Case 2 - Without phase change and with $c_{pPCM} = 21$; Case 3 - With phase change and with $c_{pPCM} = 1$. Note that the two materials are still referred to as HCM and PCM, even though no phase change occurs to keep the descriptions consistent. The phase change is designed so that the effective heat capacity of the PCM is equal to the high $c_{pPCM}$ when phase change is occurring and low $c_{pPCM}$ else. The $c_{pPCM}$ used for the cases are shown in Figure 10. These cases are chosen to see how effective the TO can make use of the phase-change process itself and to determine what favourable design features can be attributed to the inclusion of phase change and what is due to the increase in thermal storage.

<figure><img src="image_11.png" /><figcaption>Fig. 10 Effect heat capacity for the three cases investigated: Case 1 - No phase change with low $c_{pPCM}$; Case 2 - No phase change with high $c_{pPCM}$; Case 3 - With phase change.</figcaption></figure>
<table><caption>**Table 4** Results from cross-check designs based on full time history, showing normalised $\phi$ for each design tried in three different cases: Case 1 - No phase change and low $c_{pPCM}$; Case 2 - No phase change and high $c_{pPCM}$; Case 3 - With phase change.</caption><thead><tr><th rowspan="2">Tested at</th><th colspan="3">Optimised for</th></tr><tr><th>Case 1</th><th>Case 2</th><th>Case 3</th></tr></thead><tbody><tr><td>Case 1</td><td>1</td><td>2.96</td><td>2.63</td></tr><tr><td>Case 2</td><td>2.06</td><td>1</td><td>1.29</td></tr><tr><td>Case 3</td><td>1.59</td><td>0.99</td><td>1</td></tr></tbody></table>

of the major fins. The optimised designs for Case 2 and Case 3 consists of shorter fins close to the heat source reaching towards the cooling side, followed by a section of intermediate material, and then an expanding structure connecting the section of intermediate material to the cooling side. The designs generally have less fins reaching into the PCM than when compared to the design based on Case 1. The optimised design for Case 3 has more PCM layered in-between the HCM when compared to the optimised design for Case 2.

From the plots of $T_{elec}$ over time in Figure 11, it can be seen that all optimised design result in a reduced amplitude in the temperature oscillation. For Case 2 and Case 3, the optimised designs also reduces the time until a quasi-steady-state is reached. This makes sense as both the amplitude and the initial transient have an effect on the variance $\phi$ as the TO considers the entire 20 seconds simulated as input.

It can also be seen that all three optimised designs improved the performance of the heat sink by reducing the variance, $\phi$. Comparing the $\phi$ values to those of the initial design, the optimised design based on Case 1 shows a 50% reduction, the optimised design based on Case 2 shows a 52% reduction, and the optimised design based on Case 3 shows a 41% reduction. Overall, Case 3 has the best performance of the three cases which is a bit surprising as Case 2 has more thermal storage to smoothen the temperature oscillations. The plots of the accumulated $\phi$ over time in Figure 11, shows that the improvements in $\phi$ occurs at different points of the time history. It can be seen that the optimised design based on Case 1 reduces $\phi$ gradually over time, while the designs based on Case 2 and Case 3 reduce the variance $\phi$ drastically in the first 5 seconds of the simulated time, after which it is reduced further gradually over time. This further supports that reducing the time until quasi-steady-state has an effect on the performance for Case 2 and Case 3.

To verify the optimised designs from the TO, a cross-check is performed, where the variance $\phi$ of the designs based on the three cases are tested. If the TO is set up correctly the designs should perform best in the cases they were optimised for and worse in other cases, compared to the design that were optimised for that case. The result from the cross-check are shown in Table 4. It can be seen that the majority of the optimised designs perform best in the cases they were optimised for, indicating that the TO generally is effective at optimising the design in most cases. However, when testing at Case 3, the optimised design for Case 2 performs ever so slightly better than the optimised design for Case 3. This is likely a result of the complicated objective functional or the severe non-linearity introduced with the apparent heat
<figure><img src="image_12.png" /><figcaption>Fig. 12 The thermal energy field and heat flux in the y direction at t = 19.24s, t = 19.48s, and t = 19.76s. Based on the initial design under the conditions of Case 1</figcaption></figure>

capacity method, which makes it difficult for the TO to optimise the designs without getting stuck in a local minimum.

### 5.1.1 Thermal energy and heat flux fields

To further understand the optimised designs, the time history of the thermal energy field and heat flux in the y-direction are investigated. The thermal energy field is defined as:

$$E_{term}^{(n)} = \rho_m c_p (\tilde{\rho}) T^{(n)} + (1 - \tilde{\rho}) L_{heat} f_{melt}^{(n)} \quad (18)$$

where $f_{melt}^{(n)}$ is the melt fraction. The melt fraction is a piece-wise function that varies linearly from zero, when no material has change phase to one when all material has changed phase, over the phase change temperature range $\Delta T_{melt}$. For this paper it is approximated with smooth Heaviside step functions:

$$f_{melt}^{(n)} = \frac{T^{(n)} - \left(T_{melt} - \frac{\Delta T_{melt}}{2}\right)}{\Delta T_{melt}} \left( \frac{1}{1 + e^{-2k_H \left(T^{(n)} - \frac{\Delta T_{melt}}{2}\right)}} \right) - \frac{T^{(n)} - \left(T_{melt} + \frac{\Delta T_{melt}}{2}\right)}{\Delta T_{melt}} \left( \frac{1}{1 + e^{-2k_H \left(T^{(n)} + \frac{\Delta T_{melt}}{2}\right)}} \right) \quad (19)$$

The results of all the simulations are presented as animated GIFs in the supplementary material. In order to illustrate the transient process in this paper, plots of the thermal energy field and heat flux are presented at t = 19.24s, t = 19.48s, and t = 19.76s. The chosen times correspond to the time steps closest to the highest, medium, and lowest values of the oscillating heat input rate $q_{HS}(t)$ for the last oscillation period simulated.

First, the fields for the initial design is investigated. The overall characteristics of the physics are very similar for all three cases, so only Case 1 is shown in Figure 12. It can be seen that both the thermal energy and the heat flux in the y-direction are smooth. The stored thermal energy is highest when $q_{HS}(t)$ is at its medium level and
<figure><img src="image_13.png" /><figcaption>**Fig. 13** The thermal energy field and heat flux in the y direction at high $q_{HS}(t)$, medium $q_{HS}(t)$, and low $q_{HS}(t)$, for the optimised design based on Case 1</figcaption></figure>

most of the stored energy is located close to the heat source. The thermal energy is relatively low when $q_{HS}(t)$ is at its lowest level. From the flux in the y-direction, it can be seen that the heat is transferred from the heat source to the cooling side when $q_{HS}(t)$ is at its high and medium level. When $q_{HS}(t)$ is at its lowest level only very limited heat transfer occurs. This indicates that the initial design of the heat sink transports most of the thermal energy from heat source to the cooling side, which leads to the large variance seen in Figure 11.

For the optimised design based on Case 1, the thermal energy and the heat flux in the y-direction can be seen in Figure 13. It can be seen that most of the thermal energy is stored in the PCM in the bottom half of the PCM heat sink when $q_{HS}(t)$ is at its high and medium level. When $q_{HS}(t)$ is at its low level, there is a lot less thermal energy in the heat sink. It can be seen that the thermal energy is lead out into all the corners of the heat sink resulting in a fairly even energy distribution, which probably is a result of the PCM and HCM having the same heat capacity. By looking at the heat flux in the y-direction, the heat transfer into the PCM in the corners of the heat sink can be observed in the high positive heat flux in the smaller fins reaching into the PCM when $q_{HS}(t)$ is at its high and medium level. When $q_{HS}(t)$ is at its low level, a negative heat flux can be seen in the smaller fins, indicating that the thermal energy is drawn out of the PCM.

For the optimised design based on Case 2, the thermal energy and the heat flux in the y-direction can be seen in Figure 14. The design looks very different from the design based on Case 1. From the thermal energy, it can be seen that energy is stored very locally in the PCM closest to the heat source when $q_{HS}(t)$ is at its high and medium level and is expended when $q_{HS}(t)$ is at its low level. This design seems to keep a lot of the PCM isolated behind PCM resulting in it staying resulting in it a relatively low amount of energy being stored in it throughout the presented time frames. The section with intermediate material also seems to limit heat transfer away from the lower half of the PCM heat sink, as the thermal energy in the upper half of the heat sink is a lot lower than close to the heat sink. Looking at the heat flux in the y-direction, it can be seen that the highest heat flux, when $q_{HS}(t)$ is high, is located
<figure><img src="image_14.png" /><figcaption>**Fig. 14** The thermal energy field and heat flux in the y direction at high $q_{HS}(t)$, medium $q_{HS}(t)$, and low $q_{HS}(t)$, for the optimised design based on Case 2</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. 15** The thermal energy field and heat flux in the y direction at high $q_{HS}(t)$, medium $q_{HS}(t)$, and low $q_{HS}(t)$, for the optimised design based on Case 3</figcaption></figure>

close to the heat source. When $q_{HS}(t)$ is medium and low, it is located just above the section with intermediate material. Furthermore, when $q_{HS}(t)$ is low, a negative heat flux can be seen close to the heat source indicating heat transfer back to the heat source. It is the interpretation of the authors that the optimised design is aimed at reducing the transient time as much as possible. As the shorter fins close to the heat source are surrounded by the low thermally-conductive PCM, the optimised design keeps the heat from being stored in the PCM further away which reduces the thermal storage capacity available to the heat sink leading to a shorter transient period.

For the optimised design based on Case 3, the thermal energy and the heat flux in the y-direction can be seen in Figure 15. In general the thermal energy looks very similar to the case 2. Thermal energy is stored very locally in the PCM closest to the heat source when $q_{HS}(t)$ is at its high and medium level and is expended when $q_{HS}(t)$ is low. This design also keeps a lot of the PCM isolated behind PCM resulting in it a relatively low amount of energy being stored in it throughout the presented time frames. The section with intermediate material also seems to limit heat transfer away
from the lower half of the PCM heat sink, as the thermal energy in upper half of the
heat sink is a lot lower than close to the heat sink. When looking at the heat flux in the
y-direction, it can be seen that the flux is lot lower when $q_{HS}(t)$ is low, when compared
to the Case 2 design and that there is no negative heat flux close to the heat source.

In general it can be seen that the TO is able to produce optimised designs that
reduce the variance in $T_{elec}$ compared to the initial designs. It was found that apart
from the amplitude of the temperature oscillation, the transient time until the quasi-
steady-state is reached also has an effect on the variance and that the optimised designs
exploit this. Especially for the designs based on Case 2 and Case 3 see a reduction in
contributions to $\phi$ in the first 5 seconds of the simulated time.

The investigation of thermal energy and heat flux in the y-direction shows that
the optimised designs based on Case 1 tries to spread out the thermal energy into the
PCM when the $q_{HS}(t)$ is at its high and medium level, while the designs based on
Case 2 and Case 3 keep the thermal energy very close to the heat source. For Case 2
and Case 3, the heat transfer to the cooling side is further restricted with a section
consisting of intermediate material, which limits the heat flux in the y-direction.

## 5.2 Considering only the quasi-steady-state

There are three ways the TO can reduce the variance $\phi$: by reducing the amplitude in the temperature oscillation; by reducing the transient period to move the mean closer to the center of the oscillation; or by reducing the mean which also reduces the effect of the initial transient period. In terms of preventing thermal fatigue, the amplitude of the temperature oscillation is the most important factor. Therefore, the transient part is removed by only considering the physics once the quasi-steady-state is reached - a point which will depend on the current design. It would be a cumbersome exercise to derive and implement the adjoint problem by hand, since this introduces design-dependent simulation times and initial conditions for the integration. However, this process is fully automated with dolfin-adjoint and, thus, relatively easy to implement. This is expected to result in optimised designs that are better at reducing the amplitude of the temperature oscillation compared to the optimised designs based on the entire 20 seconds temperature history.

For this TO, the simulation is run until the maximum relative change in the tem-
perature oscillation period is below 1%. The $T_{elec}$ signal from the last period is then
used to calculate the mean and variance. In order to make comparisons, the TO was
performed on all the three cases. The TO takes 10 hours for Case 1, 40 hours for Case
2, and 34 hours for Case 3 to produce the optimised results on 1 core of Intel Xeon
Gold 6130 @ 2.776GHz. This is a big reduction for Case 1, as the simulation reaches
the the quasi-steady state a lot earlier than 20 seconds, due to it low thermal stor-
age. Where as the computational time for Case 2 and Case 3 are comparable to the 20
second simulations, as their simulations reach steady state around the 20 second mark.

5.2.1 Optimised designs and performance

The resulting optimised designs and their performances compared to the initial design
and 20 second designs are presented in Figure 16. The temperature at the heat source,
<figure><img src="image_16.png" /><figcaption>**Fig. 16** Optimised designs based on the quasi-steady-state compared to the 20 second designs and the initial design for the three cases. The $\phi$ values are computed using $T_{elec}$ from the last period only. The plots shows $T_{elec}$ over time at the last two periods of the simulated time for the three cases.</figcaption></figure>

$T_{elec}$, is plotted for the last two periods after the quasi-steady state is reached. The variance, $\phi$ is computed based on the last period for all cases, in order to make the results comparable. It can be seen that the optimised design based on Case 1 and its $T_{elec}$ over time look very similar to those based on the 20 seconds simulation, which makes sense as Case 1 only has a short transient period. The optimised designs based on Case 2 and Case 3, however, differ a lot from the 20 second designs. The design based on Case 2 consists of one major fin and two less prominent fins of intermediate material connecting the heat source to the cooling side, shorter fins that is reaching from the heat source to the PCM in the lower half of the heat sink, and intermediate material which is located close to the heat source. The optimised design based on Case 3 consist of three major fin structures reaching from the heat source towards the side
<table><caption>**Table 5** Results from cross-check designs based on the quasi-steady-state, showing normalised $\phi$ for each design tried in three different cases: Case 1 - No phase change and low $c_{pPCM}$; Case 2 - No phase change and high $c_{pPCM}$; Case 3 - With phase change.</caption><thead><tr><th rowspan="2">Tested at</th><th colspan="3">Optimised for</th></tr><tr><th>Case 1</th><th>Case 2</th><th>Case 3</th></tr></thead><tbody><tr><td>Case 1</td><td>1</td><td>1.39</td><td>1.26</td></tr><tr><td>Case 2</td><td>1.44</td><td>1</td><td>1.01</td></tr><tr><td>Case 3</td><td>1.32</td><td>1.36</td><td>1</td></tr></tbody></table>

with cooling, smaller fins reaching from the heat source into the PCM in the lower half of the heat sink and intermediate material at the ends of the fins leaving more PCM close to the heat source when compared to the Case 2 design.

The optimised designs based on Case 2 and 3 show a reduction in the amplitude of the oscillation of $T_{elec}$ and an increase in the average $T_{elec}$ when compared to the 20 second designs. Comparing the variance, $\phi$, to the 20 second designs, the optimised design based on Case 2 shows a 78% reduction and the optimised design based on Case 3 shows a 32% reduction. This is expected, as the TO can only reduce the $\phi$ by reducing the amplitude of the oscillation of $T_{elec}$ when considering the quasi-steady-state.

It can be seen that when considering only the quasi-steady-state, the temperature variance, $\phi$, is lower in Case 2 than for Case 3. This is contrary to what was found when looking at the 20 second designs, but makes sense from a physics perspective. When reducing the transient part has a big impact on the variance, then Case 3 must be better as it has a lower thermal storage capacity and therefore can reach the quasi-steady-state earlier. While in this case where the transient part is neglected, the higher thermal storage capacity of Case 2 allows it to store more heat to smoothen out the temperature oscillations.

To verify the optimised designs from the TO, a cross-check is performed, where the variance $\phi$ of the designs based on the three cases are tested. If the TO is set up correctly the designs should perform best in the cases they were optimised for and worse in other cases, compared to the design that were optimised for that case. The result from the cross-check are shown in Table 5. It can be seen that all the optimised designs are performing best in the cases they were optimised for, indicating that the TO is effective at optimising the design for this problem. However, it can be seen that Case 2 has a very similar performance to Case 3 when tested at Case 2. As it, in this case, is Case 2 that is close to performing worse than the other designs, it is postulated that the complicated objective function has a larger effect on the TO reaching a local minimum than the non-linearity on the apparent heat capacity method.

### 5.2.2 Thermal energy and heat flux fields

The optimised design are further investigated by looking at the thermal energy and flux in the y-direction. The results of this investigation are shown in Figures 17, 18, and 19.
<figure><img src="image_17.png" /><figcaption>**Fig. 17** The thermal energy field and heat flux in the y direction at high $q_H S(t)$, medium $q_H S(t)$, and low $q_H S(t)$ for the optimised design based on the quasi-steady-state under Case 1</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Fig. 18** The thermal energy field and heat flux in the y direction at high $q_H S(t)$, average $q_H S(t)$, and low $q_H S(t)$ for the optimised design based on the quasi-steady-state under Case 2</figcaption></figure>

For the optimised design based on Case 1, the thermal energy and the heat flux in the y-direction can be seen in Figure 17. It can be seen that the thermal energy and the heat flux in the y-direction are very similar to what was seen in Figure 13 for the 20 seconds design. This makes sense as the designs are very similar.

For the optimised design based on Case 2, the thermal energy and the heat flux in the y-direction can be seen in Figure 18. Looking at the thermal energy, it can be seen that the thermal energy is stored in much more evenly in the PCM, when compared to the 20 second design. Most of the thermal energy is stored in the PCM in the sides. The differences in the thermal energy over the presented time steps are fairly small compared to the general level of thermal energy in the heat sink, making the hard to detect in the presented times frames. However, a higher thermal energy can be seen in PCM just outside the intermediate material close to the heat sink when $q_H S(t)$ is at its medium level and low level than when $q_H S(t)$ is at its high level. The small change in thermal energy indicates a smaller temperature oscillation, and the reason for seeing the high thermal energy when $q_H S(t)$ is at its medium and low level can be explained with the slow thermal response due to the high thermal storage and low
<figure><img src="image_19.png" /><figcaption>**Fig. 19** The thermal energy field and heat flux in the y direction at high $q_H S(t)$, average $q_H S(t)$, and low $q_H S(t)$ for the optimised design based on the quasi-steady-state under Case 3</figcaption></figure>

heat conductivity in the intermediate material around the heat source. From the flux in the y-direction, it can be seen that the flux in the upper half of the heat sink stays approximately the same over time, indicating that the heat transferred when $q_{HS}(t)$ is low must come from stored energy in the bottom half of the PCM heat sink.

For the optimised design based on Case 3, the thermal energy and the heat flux in the y-direction can be seen in Figure 19. Looking at the thermal energy, it can be seen that the thermal energy is stored in much more evenly in the PCM, when compared to the 20 second design. The majority of the stored thermal energy is located in the PCM at the sides and between the major fins close to the heat source. The fact that the PCM is placed closer to the heat source, than in the Case 2 design, makes sense as the PCM needs to enter the phase change temperature range before it obtains its increased thermal storage - whereas for Case 2, the full thermal storage potential is achieved at any temperature. From the flux in the y-direction, it can be seen that the flux in the upper half of the heat sink stays approximately the same over time, again indicating that the heat transferred when $q_{HS}(t)$ is low must come from stored energy in the bottom half of the PCM heat sink.

From the investigation of the optimised designs based on the quasi-steady-state, it is found that the amplitude of the temperature oscillation can be better reduced by only considering the quasi-steady-state as input for the TO. This has the greatest impact on designs based on Case 2 and Case 3, which utilise a lot more of the PCM for storing thermal energy compared to the 20 seconds designs resulting in a more even distribution of thermal energy and a final lower variance, $\phi$. It is found that reaching the quasi-steady-state is computationally faster for Case 1 and Case 3, but took longer time for Case 2, when comparing the computation times from the TO of the 20 seconds designs.

## 5.3 Varying time step sizes

Reaching the quasi-steady-state can take a very long simulation time, especially if the thermal storage is high, which increases the computational time. Therefore, the possibility of using varying time step sizes is investigated.
For this investigation, the simulation initially uses large time steps, scaled by the thermal diffusivity of the PCM, until the heat transfer rate over the cooling side has surpassed 95% of the average heat transfer from the heat source. At that point, the time step is reduced to the time steps size used in the other simulations and the simulation is run until the maximum relative change in the temperature oscillation period is below 1%. The initial large time step size is computed by the following expression:

$$ \Delta t_{large} = L_y^2 \frac{\rho_{PCM} c_{PCM}}{k_{PCM}} \qquad (20) $$

which characterises the diffusion time scale. As for the previous quasi-steady-state designs, $T_{elec}$ from the last period is then used to calculate the variance, $\phi$. In order to make comparisons, the TO is performed on the same three cases from earlier sections. The TO using the varying time step size takes 6 hours for Case 1, 10 hours for Case 2, and 13 hours for Case 3 to produce the optimised results on a single core of an Intel Xeon Gold 6130 @ 2.776GHz. This is equivalent to a 40% to 75% reduction in computational time compared to using a constant time step size.

Figure 20 shows the resulting optimised designs and their performances, compared to the optimised designs based on quasi-stead-state from Section 5.2, which use a constant time step size. It can be seen that the optimised designs generally look similar, but differ in the amount of fins, placement of major fins, and the amount of intermediate material. From the temperature signal over time, it can be seen that for Case 1 and Case 2 there is a good agreement in $T_{elec}$ between the two methods for time stepping. However, for Case 3 there is a rather large difference, where the design based on varying time step size results in a higher average and a larger amplitude of the oscillation. This error is probably caused by the initial large time step size introducing a substantial error due to ignoring the full non-linearity introduced to the model with the apparent heat capacity method, as described in Section 2.5. This size of the error could be reduced by decreasing the large time step size, but this increases the computational time again and makes the varying time step size approach less effective.

## 5.4 Effect of forcing discrete designs

Thus far, a composite material has been considered, where intermediate material density variables represent a porous HCM skeleton embedded with PCM. This causes the manufacture to be more expensive, as this requires a precise control of the density of a HCM foam essentially. Therefore, the effects of forcing discrete designs was investigated to see whether discrete designs are possible without losing too much performance.

### 5.4.1 Explicit penalisation

As briefly discussed in Section 2.2, classical material penalisation approaches, such as SIMP, cannot be applied successfully to the presented design problem. Thus, to push the optimised designs towards being discrete, the objective functional was modified
<figure><img src="image_20.png" /><figcaption>Fig. 20 Optimised designs based on the quasi-steady-state with adaptive time step sizes compared to optimised designs based on the quasi-steady-state with a constant time step size for the three cases. The plots shows $T_{elec}$ over time at the last two periods of the simulated time for the three cases.</figcaption></figure>

using the following explicit penalty on the design variables:

$$f_0 = \phi + \alpha \int_{\Omega} \rho(1 - \rho)dA \qquad (21)$$

where $\alpha$ is the explicit penalisation factor, used to control the level of explicit penalisation. The integral $\int_{\Omega} \rho(1 - \rho)dA$ has a minimum when $\rho$ is equal to zero or one, whereby the non-discrete designs are penalised. It should be noted that the explicit penalty is applied on the mathematical design field, $\rho$, not the physical density field, $\tilde{\rho}$. This approach is prone to converge to poor local minima if the $\alpha$ value is set to
<figure><img src="image_21.png" /><figcaption>Fig. 21 Overview of optimised designs for $t_{fin} = 20s$ of three different cases with varying amount of explicit penalisation.</figcaption></figure>

high. Therefore, a continuation approach is used, where the optimisation is allowed to converge before the explicit penalisation factor, $\alpha$, is gradually increased.

The non-discreteness of a design is quantified with the measure of non-discreteness, $M_{nd}$, defined by Sigmund (2007) as the following functional:

$$M_{nd} = 100 \frac{4 \int_{\Omega} \rho(1 - \rho)dA}{\int_{\Omega} 1dA} \qquad (22)$$

This normalises the integral from (21) to be zero for an entirely discrete design and 100 for an entirely non-discrete design with $\rho = 0.5$.

### 5.4.2 Optimised designs and performance

In order to investigate the effect of removing the intermediate material, explicit penalisation is added optimised designs based on the quasi-steady-state for all three cases, presented in Section 5.2. All optimisations reached the defined convergence criteria. The resulting optimised designs can be seen in Figure 21. In general, it can be seen that the optimised designs becomes more discrete and simple as the explicit penalisation factor $\alpha$ is increased. From the optimised designs, it can be seen that the smaller fins reaching into the PCM disappear and the HCM is moved to the major fin structures as explicit penalisation is increased. This makes sense as the filter in combination with a high penalisation of intermediate material effectively penalises the boundary between PCM and HCM. It is found that for higher explicit penalisation factors, $\alpha > 10^0$, further penalisation does not cause any major design changes.
<figure><img src="image_22.png" /><figcaption>Fig. 22 Variance in $T_{elec}$, $\phi$, plotted against the measure of non discreteness $M_{nd}$ for the three cases.</figcaption></figure>

<table><caption>**Table 6** Results from cross-check, showing normalised $\phi$ for each design tried in 3 different cases: Case 1 - No phase change and low $c_{pPCM}$; Case 2 - No phase change and high $c_{pPCM}$; Case 3 - With phase change. The designs are optimised using the quasi-steady-state.</caption><thead><tr><th rowspan="3">Tested at</th><th colspan="9">Optimised for</th></tr><tr><th colspan="3">α = 0</th><th colspan="3">α = 10<sup>-1</sup></th><th colspan="3">α = 10<sup>1</sup></th></tr><tr><th>Case 1</th><th>Case 2</th><th>Case 3</th><th>Case 1</th><th>Case 2</th><th>Case 3</th><th>Case 1</th><th>Case 2</th><th>Case 3</th></tr></thead><tbody><tr><td>Case 1</td><td>1</td><td>1.39</td><td>1.26</td><td>1</td><td>2.02</td><td>1.61</td><td>1</td><td>1.48</td><td>1.05</td></tr><tr><td>Case 2</td><td>1.44</td><td>1</td><td>1.01</td><td>1.01</td><td>1</td><td>0.93</td><td>0.70</td><td>1</td><td>0.71</td></tr><tr><td>Case 3</td><td>1.32</td><td>1.36</td><td>1</td><td>1.02</td><td>1.35</td><td>1</td><td>1.12</td><td>1.52</td><td>1</td></tr></tbody></table>

To investigate the relation between the performance and discreteness, the variance, $\phi$, was plotted against the measure of non-discreteness, $M_{nd}$, for the three cases with different levels of explicit penalisation. The result can be seen in Figure 22 and shows that $\phi$ generally increases as $M_{nd}$ decreases. However, it can be seen that the optimised designs generally can become more discrete without worsening the performance too much. By allowing a 10% increase in the variance, $\phi$, from the optimised design without penalisation, the measure of non-discreteness, $M_{nd}$, can be reduced to $M_{nd} = 22$ with $\alpha = 10^{-1}$ for Case 1, $M_{nd} = 20$ with $\alpha = 10^{-2}$ for Case 2, and $M_{nd} = 22$ with $\alpha = 10^{-2}$ for Case 3.

To verify the optimised designs, a cross-check is performed, where the performance of the designs with varying explicit penalisation factors are tested. To attribute significance to design characteristics, the designs should perform best in the cases they were optimised for and worse in other cases, compared to the design that were optimised for that case. For the cross-check $\alpha = [0, 10^{-1}, 10^{1}]$ were considered. Table 6 shows that the majority of the optimised designs are performing best in the cases they were optimised for, indicating that the TO is effectively optimising the design in most
cases. However, the optimised designs for Case 2 when $\alpha = 10^{-1}$ and $\alpha = 10^1$ are not the best performing designs when tested at Case 2. When $\alpha = 10^{-1}$, the design optimised for Case 3 performs slightly better than Case 2 and when $\alpha = 10^1$ both the designs optimised for Case 1 and Case 3 perform significantly better than Case 2. This indicates that the explicit penalisation gets the optimised design based on Case 2 stuck at a local minimum. A cause for this could be the initial design based on Case 2 only contains one major fin instead of two major fins found in the optimised designs based on Case 1 and Case 3. As the explicit penalisation removes the smaller fins close to the heat source the contact boundary between the HCM and PCM is reduced, resulting in a reduction in the thermal storage capacity. As the two major fin designs have a larger contact boundary between the HCM and PCM than a single major fin, it makes sense that the optimised designs based on Case 1 and Case 3 perform better than optimised design based on Case 1, a the penalisation is increased.

# 6 Discussion

From the presented results, it can clearly be seen that it is possible to generate designs that reduce the temperature oscillation by applying topology optimisation to the layout of PCM and HCM. The results show that a PCM heat sink has the best performance in Case 2, indicating that actual PCM with phase change is not as effective at storing heat as a material with a constant higher heat capacity. However, in practice not a lot of materials have a heat capacity to match the thermal storage capacity equivalent to a material in phase change. So the quest for the most optimal layout of PCM is still relevant. Furthermore, the effects of natural convection in the melted PCM is neglected in this paper, which results in a lower effective heat transfer through the PCM compare to the real world. This could further affect the ranking of PCM versus solid material with high heat capacity.

The physical model, material properties, and so on, used to model the heat sink in this paper, are purely academic and the designs produced are therefore not directly applicable to any real life problems. The reason for this was to investigate the effect of certain parameters in a systematic way. For real life application, the PCM will most likely have a smaller phase change temperature range compared to the general temperature ranges and have a significantly higher latent heat of fusion, making the problem even more non-linear. For a real life problem, the transient model is therefore expected to require significantly smaller time steps or solving the full non-linear system every time step. Both will increase the computational time significantly.

The computational times reported in this paper are relatively high (up to 40 hours) for the relatively small models with only 40000 degrees of freedom. Moving forwards, this can be sped up by running the TO code in parallel, but the limited size of the spatial discretisation limits the scalability. Furthermore, since the physical model requires a transient solver to model the phase change, the time stepping part of the code will inherently be serial and, thus, the speed up from parallel computing will be limited to the scalability of the solver for each time step. The approach of using varying time step size shows some promise in reducing the computational time, but also shows an
increase in error for the case the phase change. This error could probably be reduced
with a more sophisticated method of adapting the time step size.

7 Conclusion

This paper aimed to improve the understanding of the optimal of PCM and HCM
layout in heat sinks embedded with PCM, through the use of TO to reduce the temper-
ature oscillation in an electronic component resulting from a transient cyclic thermal
load. To do this, a Python code capable of TO was made using FEniCS for solving
the physical problem, dolfin-adjoint for performing the adjoint sensitivity analysis,
and an open-source implementation of GCMMA for solving the optimisation prob-
lem. For the physical model, the PCM heat sink was modelled as a transient thermal
diffusion problem neglecting the natural convection occurring in the liquid melted
PCM. The phase change was modelled with the apparent heat capacity method and
material interpolation was handled with an analytical homogenisation approach. To
quantify the temperature oscillation, the temporal variance of the spatially-averaged
temperature at the electrical component, was used as the objective functional for the
optimisation.

The TO was performed for three different cases: Case 1 - No phase change with low $c_{pPCM}$; Case 2 - No phase change with high $c_{pPCM}$; Case 3 - With phase change. When using the entire time history from a 20 seconds simulation for computing $\phi$, it was found that for all three cases the optimised designs performed better than the initial design. Looking at the temperature at the electrical component over time, it could be seen that optimised designs based on Case 2 and Case 3 reduced the variance by reducing the amplitude of the temperature oscillations as well as reducing the transient time until the quasi-steady-state is reached. This made it hard to attribute design features to a reduction in temperature oscillations alone.

Therefore, a second batch was performed using only the temperatures of one period
after the quasi-steady-state was reached. For Case 1, it had a limited effect as the 20
second design had reached the quasi-steady-state early in the simulated time. How-
ever, for Case 2 and Case 3, the optimised designs showed a substantial decrease in
the temperature oscillations when compared to the 20 second design. Comparing the
variance to the 20 second designs, the optimised design based on Case 2 showed a 78%
reduction and the optimised design based on Case 3 showed a 32% reduction. Overall,
it was found that Case 2 had the lowest variance indicating that the available thermal
storage is important for the reduction of the amplitude of the temperature oscillations.

In an attempt to reduce the computational time, a simple varying time step size
approach was investigated. It was found to reduce the computational times signifi-
cantly for Case 2 and Case 3, with a reduction of 40 and 75% in overall computational
time, respectively. Comparing the temperatures over time with the results from using
a constant time step size, showed that the varying time step size introduced a small
error for Case 1 and Case 2, while introducing a large error in Case 3. It would be
interesting to investigate more sophisticated adaptive time stepping methods to see if
the errors can be reduced without reintroducing too much computational time.
Lastly, the possibility of removing intermediate material in order to improve man- ufacturability and production costs was investigated. This was done by implementing an explicit penalisation on the intermediate material. It was found that the explicit penalisation reduced the amount of intermediate material and removed smaller fin structures, which in turn lead to a increase in the variance. By allowing a 10% increase in the variance, from the optimised design without penalisation, the measure of non- discreteness, $M_{nd}$, can be reduced to $M_{nd} = 22$ with $\alpha = 10^{-1}$ for Case 1, $M_{nd} = 20$ with $\alpha = 10^{-2}$ for Case 2, and $M_{nd} = 22$ with $\alpha = 10^{-2}$ for Case 3.

Hereby, it can be concluded that gradient-based TO can be used to optimise the PCM and HCM layout in a heat sink embedded with PCM in order to reduce the temperature oscillation induced by a cyclic thermal load. The biggest reduction in the temperature oscillation was found using the temperatures from the last period once the quasi-steady-state has been reached as input for calculating the variance. The results indicate that the thermal storage capacity has an important influence on reducing the variance.

For future work, the performance of the optimised designs consisting of varying composite microstructure would need to be verified using dehomogenised designs and validated experimentally to verify the physical model. Furthermore, the influence of geometry, material, and boundary conditions should be investigated in the context of a real life problem with real physical material properties.

## Statements and Declarations

**Competing Interest.** The authors declare that they have no conflict of interest.

**CRediT authorship contribution statement.** Mark BM Christensen: Methodology, software, verification, investigation, computations, writing - original draft, writing - review & editing, visualisation. **Joe Alexandersen:** Conceptualisation, methodology, investigation, resources, supervision, writing - review & editing, project administration, funding acquisition.

**Acknowledgments.** Part of the computations were performed on the UCloud interactive HPC system, which is managed by the eScience Center at the University of Southern Denmark.

**Replication of results.** A minimal working example is available on GitHub: https://github.com/sdu-multiphysics/heatSinkPCM/.

## Appendix A Derivation of $k_T(\rho)$ with homogenisation

The base cell can be viewed as a composite wall. Assuming one-direction heat transfer, the thermal conductivity $k_T$ of the base cell can be derived as the total thermal resistance of the composite wall:

$$k_{Total} = \frac{L_{total}}{R_{total}A_{total}} \qquad (A1)$$
where: $k_{Total}$ its the total thermal conductivity through the composite wall; $R_{total}$ is the total thermal resistance; $L_{total}$ is the total length of the composite wall in the direction of the heat transfer; $A_{total}$ is the total cross-sectional area of the wall; $this$ thickness of the wall,

In order to compute the total thermal resistance the base cell is broken down into three distinct section and their thermal resistance are computed, see Figure A1. The

<figure><img src="image_23.png" /><figcaption>Fig. A1 By assuming one-dimensional heat transfer the base cell can be broken down into sections with their own thermal resistance</figcaption></figure>

thermal resistance depends on the material properties and dimensions the section:

$$R = \frac{L}{kA} \qquad (A2)$$

where: $L$ is the length in the direction of the heat transfer; $k$ is the thermal conductivity; $A$ is the cross sectional area through which the heat transfer occurs. With this the which the thermal resistances of the three sections are derived:

$$R_1 = \frac{\frac{a}{2}}{k_{THCM}1ths} \qquad (A3)$$

$$R_2 = \frac{1-a}{k_{THCM}(\frac{a}{2})ths} \qquad (A4)$$

$$R_3 = \frac{1}{k_{TPCM}ths} \qquad (A5)$$

where: $this$ is the thickness of the base cell (out of plane dimension); $a$ is a parameter describing the amount of HCM.

The middle section of the composite wall consists of three sections in parallel, the equivalent resistance is expressed as:

$$\frac{1}{R_{mid}} = \frac{1}{R_2} + \frac{1}{R_3} + \frac{1}{R_2} \qquad (A6)$$
leading to:

$$R_{mid} = \frac{R_2 R_3}{2R_3 + R_2} \qquad (A7)$$

Now the total thermal resistance can considered as sections in series:

$$R_{total} = R_1 + R_{mid} + R_1 \qquad (A8)$$

The total thermal conductivity becomes:

$$k_{Total} = \frac{1}{(2R_1 + R_{mid})1ths} \qquad (A9)$$

which can be reduced to the following term where $k_{Total}$ is dependent of a:

$$k_{Total}(a) = \frac{1}{\frac{a}{k_{THCM}} + \frac{1-a}{ak_{THCM}+(1-a)k_{TPCM}}} \qquad (A10)$$

The relation between a and the ρ can be expressed by the following function:

$$\rho = 1 - (1 - a)^2 \qquad (A11)$$

from which follows

$$a = 1 - \sqrt{1 - \rho} \qquad (A12)$$

Combining (A12) into (A10), $k_T$ can be expressed in terms of ρ:

$$k_T(\rho) = \frac{1}{\frac{1-\sqrt{1-\rho}}{k_{THCM}} + \frac{\sqrt{1-\rho}}{(1-\sqrt{1-\rho})k_{THCM}+(\sqrt{1-\rho})k_{TPCM}}} \qquad (A13)$$

## References

Alexandersen J, Andreasen CS (2020) A Review of Topology Optimisation for Fluid-Based Problems. Fluids 5(1):29. https://doi.org/10.3390/fluids5010029, URL https://www.mdpi.com/2311-5521/5/1/29

Alnæs M, Blechta J, Hake J, et al (2015) The FEniCS Project Version 1.5. <p>Archive of Numerical Software Vol 3:<strong>Starting Point and Frequency: </strong>Year: 2013</p>. https://doi.org/10.11588/ANS.2015.100.20553, URL http://journals.ub.uni-heidelberg.de/index.php/ans/article/view/20553, publisher: University Library Heidelberg Version Number: 1.0.0

Bianco N, Fragnito A, Iasiello M, et al (2023) Design of PCM-based heat sinks through topology optimization. Journal of Physics: Conference Series 2509(1):012001. https://doi.org/10.1088/1742-6596/2509/1/012001, URL https://iopscience.iop.org/article/10.1088/1742-6596/2509/1/012001
Choi KK, Kim NH (2005) Structural Sensitivity Analysis and Optimization 1. Mechanical Engineering Series, Springer New York, New York, NY, https://doi.org/10.1007/b138709, URL http://link.springer.com/10.1007/b138709

COMSOL (2023) COMSOL 2023. URL www.comsol.com

Dbouk T (2017) A review about the engineering design of optimal heat transfer systems using topology optimization. Applied Thermal Engineering 112:841–854. https://doi.org/10.1016/j.applthermaleng.2016.10.134, URL https://linkinghub.elsevier.com/retrieve/pii/S135943111632645X

Deetman A (2020) GCMMA-MMA-Python. URL https://github.com/arjendeetman/GCMMA-MMA-Python

Depiver Ja, Mallik S, Harmanto D (2021) Solder joint failures under thermo-mechanical loading conditions – A review. Advances in Materials and Processing Technologies 7(1):1–26. https://doi.org/10.1080/2374068X.2020.1751514, URL https://www.tandfonline.com/doi/full/10.1080/2374068X.2020.1751514

Fawaz A, Hua Y, Le Corre S, et al (2022) Topology optimization of heat exchangers: A review. Energy 252:124053. https://doi.org/10.1016/j.energy.2022.124053, URL https://linkinghub.elsevier.com/retrieve/pii/S0360544222009562

Haftka RT, Gürdal Z (1992) Elements of structural optimization, 3rd edn. No. v. 11 in Solid mechanics and its applications, Kluwer Academic Publishers, Dordrecht ; Boston

Ho J, See Y, Leong K, et al (2021) An experimental investigation of a PCM-based heat sink enhanced with a topology-optimized tree-like structure. Energy Conversion and Management 245:114608. https://doi.org/10.1016/j.enconman.2021.114608, URL https://linkinghub.elsevier.com/retrieve/pii/S0196890421007846

Hosseinizadeh S, Tan F, Moosania S (2011) Experimental and numerical studies on performance of PCM-based heat sink with different configurations of internal fins. Applied Thermal Engineering 31(17-18):3827–3838. https://doi.org/10.1016/j.applthermaleng.2011.07.031, URL https://linkinghub.elsevier.com/retrieve/pii/S1359431111003905

Iradukunda AC, Vargas A, Huitink D, et al (2020) Transient thermal performance using phase change material integrated topology optimized heat sinks. Applied Thermal Engineering 179:115723. https://doi.org/10.1016/j.applthermaleng.2020.115723, URL https://linkinghub.elsevier.com/retrieve/pii/S1359431120332051

Jauregui CM, Hyun J, Neofytou A, et al (2023) Avoiding reinventing the wheel: reusable open-source topology optimization software. Structural and Multidisciplinary Optimization 66(6):145. https://doi.org/10.1007/s00158-023-03589-7
Jegadheeswaran S, Pohekar SD (2009) Performance enhancement in latent heat thermal storage system: A review. Renewable and Sustainable Energy Reviews 13(9):2225-2244. https://doi.org/10.1016/j.rser.2009.06.024, URL https://linkinghub.elsevier.com/retrieve/pii/S1364032109001221

Kalbasi R, Afrand M, Alsarraf J, et al (2019) Studies on optimum fins number in PCM-based heat sinks. Energy 171:1088-1099. https://doi.org/10.1016/j.energy.2019.01.070, URL https://linkinghub.elsevier.com/retrieve/pii/S0360544219300726

Laasri IA, Elmaazouzi Z, Outzourhit A, et al (2022) Investigation of different topology-optimized fin structures in a cylindrical latent heat thermal energy storage unit. Thermal Science and Engineering Progress 33:101372. https://doi.org/10.1016/j.tsep.2022.101372, URL https://linkinghub.elsevier.com/retrieve/pii/S2451904922001792

Laurain A (2018) A level set-based structural optimization code using FEniCS. Structural and Multidisciplinary Optimization 58(3):1311-1334. https://doi.org/10.1007/s00158-018-1950-2, URL http://link.springer.com/10.1007/s00158-018-1950-2

Lazarov BS, Sigmund O (2011) Filters in topology optimization based on Helmholtz-type differential equations. International Journal for Numerical Methods in Engineering 86(6):765-781. https://doi.org/10.1002/nme.3072, URL https://onlinelibrary.wiley.com/doi/10.1002/nme.3072

Lazarov BS, Wang F, Sigmund O (2016) Length scale and manufacturability in density-based topology optimization. Archive of Applied Mechanics 86(1-2):189-218. https://doi.org/10.1007/s00419-015-1106-4, URL http://link.springer.com/10.1007/s00419-015-1106-4

Logg A, Mardal KA, Wells G (eds) (2012) Automated Solution of Differential Equations by the Finite Element Method: The FEniCS Book, Lecture Notes in Computational Science and Engineering, vol 84. Springer Berlin Heidelberg, Berlin, Heidelberg, https://doi.org/10.1007/978-3-642-23099-8, URL https://link.springer.com/10.1007/978-3-642-23099-8

Mezzadri F, Bouriakov V, Qian X (2018) Topology optimization of self-supporting support structures for additive manufacturing. Additive Manufacturing 21:666-682. https://doi.org/10.1016/j.addma.2018.04.016, URL https://linkinghub.elsevier.com/retrieve/pii/S2214860418301519

Michaleris P, Tortorelli DA, Vidal CA (1994) Tangent operators and design sensitivity formulations for transient non-linear coupled problems with applications to elastoplasticity. International Journal for Numerical Methods in Engineering 37(14):2471-2499. https://doi.org/10.1002/nme.1620371408, URL https://onlinelibrary.wiley.com/doi/10.1002/nme.1620371408
Mitusch S, Funke S, Dokken J (2019) dolfin-adjoint 2018.1: automated adjoints for FEniCS and Firedrake. Journal of Open Source Software 4(38):1292. https://doi.org/10.21105/joss.01292, URL http://joss.theoj.org/papers/10.21105/joss.01292

Peremans B, Blommaert M, Baelmans M (2023) Topology optimization of a rectangular phase change material module. Journal of Energy Storage 70:107891. https://doi.org/10.1016/j.est.2023.107891, URL https://linkinghub.elsevier.com/retrieve/pii/S2352152X23012884

Pizzolato A, Sharma A, Maute K, et al (2017a) Design of effective fins for fast PCM melting and solidification in shell-and-tube latent heat thermal energy storage through topology optimization. Applied Energy 208:210-227. https://doi.org/10.1016/j.apenergy.2017.10.050, URL https://linkinghub.elsevier.com/retrieve/pii/S0306261917314642

Pizzolato A, Sharma A, Maute K, et al (2017b) Topology optimization for heat transfer enhancement in Latent Heat Thermal Energy Storage. International Journal of Heat and Mass Transfer 113:875-888. https://doi.org/10.1016/j.ijheatmasstransfer.2017.05.098, URL https://linkinghub.elsevier.com/retrieve/pii/S0017931017303034

Pizzolato A, Sharma A, Ge R, et al (2020) Maximization of performance in multi-tube latent heat storage - Optimization of fins topology, effect of materials selection and flow arrangements. Energy 203:114797. https://doi.org/10.1016/j.energy.2019.02.155, URL https://linkinghub.elsevier.com/retrieve/pii/S0360544219303548

Qian X (2017) Undercut and overhang angle control in topology optimization: A density gradient based integral approach. International Journal for Numerical Methods in Engineering 111(3):247-272. https://doi.org/10.1002/nme.5461, URL https://onlinelibrary.wiley.com/doi/10.1002/nme.5461

Sahoo SK, Das MK, Rath P (2016) Application of TCE-PCM based heat sinks for cooling of electronic components: A review. Renewable and Sustainable Energy Reviews 59:550-582. https://doi.org/10.1016/j.rser.2015.12.238, URL https://linkinghub.elsevier.com/retrieve/pii/S1364032115016214

See Y, Ho J, Leong K, et al (2022) Experimental investigation of a topology-optimized phase change heat sink optimized for natural convection. Applied Energy 314:118984. https://doi.org/10.1016/j.apenergy.2022.118984, URL https://linkinghub.elsevier.com/retrieve/pii/S0306261922003932

Sigmund O (2007) Morphology-based black and white filters for topology optimization. Structural and Multidisciplinary Optimization 33(4-5):401-424. https://doi.org/10.1007/s00158-006-0087-x, URL http://link.springer.com/10.1007/s00158-006-0087-x

Sigmund O, Petersson J (1998) Numerical instabilities in topology optimization: A survey on procedures dealing with checkerboards, mesh-dependencies and local minima.
Structural Optimization 16(1):68-75. https://doi.org/10.1007/BF01214002, URL
http://link.springer.com/10.1007/BF01214002

Svanberg K (2002) A Class of Globally Convergent Optimization Methods Based on Conservative Convex Separable Approximations. SIAM Journal on Optimization 12(2):555-573. https://doi.org/10.1137/S1052623499362822, URL http://epubs.siam.org/doi/10.1137/S1052623499362822

Tian Y, Liu X, Xu Q, et al (2021) Bionic topology optimization of fins for rapid latent heat thermal energy storage. Applied Thermal Engineering 194:117104. https://doi.org/10.1016/j.applthermaleng.2021.117104, URL https://linkinghub.elsevier.com/retrieve/pii/S1359431121005469

Varanasi S, Ananthasuresh GK (2006) Analysis and Topology Optimization of Heat Sinks with a Phase-Change Material on COMSOL Multiphysics™ Platform. COMSOL, Bangalore, URL https://www.researchgate.net/profile/G-Ananthasuresh-2/publication/267679913_Analysis_and_Topology_Optimization_of_Heat_Sinks_with_a_Phase_Change_Material_on_COMSOL_Multiphysics_Platform/links/548eca3f0cf214269f247d97_Analysis-and-Topology-Optimization-of-Heat-Sinks-with-a-Phase-Change-Material-on-COMSOL-Multiphysics.pdf

Yan J, Xiang R, Kamensky D, et al (2022) Topology optimization with automated derivative computation for multidisciplinary design problems. Structural and Multidisciplinary Optimization 65(5):151. https://doi.org/10.1007/s00158-022-03168-2

Yao Q, Zhao C, Zhao Y, et al (2021) Topology optimization for heat transfer enhancement in latent heat storage. International Journal of Thermal Sciences 159:106578. https://doi.org/10.1016/j.ijthermalsci.2020.106578, URL https://linkinghub.elsevier.com/retrieve/pii/S1290072920310309

Zhang X, Yang X, Zhang Y, et al (2023) Phase change heat transfer enhancement based on topology optimization of fin structure. International Journal of Heat and Mass Transfer 214:124402. https://doi.org/10.1016/j.ijheatmasstransfer.2023.124402, URL https://linkinghub.elsevier.com/retrieve/pii/S0017931023005483

Zhao M, Tian Y, Hu M, et al (2020) Topology optimization of fins for energy storage tank with phase change material. Numerical Heat Transfer, Part A: Applications 77(3):284-301. https://doi.org/10.1080/10407782.2019.1690338, URL https://www.tandfonline.com/doi/full/10.1080/10407782.2019.1690338