# Multi-scale topology optimisation of microchannel cooling using a homogenisation-based method

Hao Li, Pierre Jolivet, Joe Alexandersen

► To cite this version:

Hao Li, Pierre Jolivet, Joe Alexandersen. Multi-scale topology optimisation of microchannel cooling using a homogenisation-based method. 2024. hal-04607025

HAL Id: hal-04607025

https://hal.science/hal-04607025

Preprint submitted on 10 Jun 2024

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 Attribution 4.0 International License
# Multi-scale topology optimisation of microchannel cooling using a homogenisation-based method

Hao Li, Pierre Jolivet, Joe Alexandersen

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

## To cite this version:

Hao Li, Pierre Jolivet, Joe Alexandersen. Multi-scale topology optimisation of microchannel cooling using a homogenisation-based method. 2024. hal-04607025v3

HAL Id: hal-04607025

https://hal.science/hal-04607025v3

Preprint submitted on 9 Dec 2024

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Multi-scale topology optimisation of
microchannel cooling using a
homogenisation-based method

Hao Li¹*, Pierre Jolivet² and Joe Alexandersen¹

¹*Institute of Mechanical and Electrical Engineering, University of Southern Denmark, 5230, Odense, Denmark.

² Sorbonne Université, CNRS, LIP6, 75252, Paris, France.

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

Contributing authors: pierre@joliv.et; joal@sdu.dk;

**Abstract**

Microchannel cooling is often the preferred choice for compact heat sinks. However, widely adopted topology optimisation (TO) techniques, such as density-based and level-set methods, often struggle to generate very thin channel strips without imposing maximum length scale constraints. To address this limitation, multi-scale design methodologies have emerged. This paper builds upon recent advances in de-homogenisation techniques to contribute to the multi-scale design of microchannels for cooling applications. We start by selecting a single-class microstructure and employ numerical homogenisation to build an offline library. This library is then fed in online macro-scale topology optimisation, where both microstructure parameters and local orientation fields are optimised. By using a sawtooth-function-based mapping, the de-homogenised results capture fine details across different length scales through a unique homogenised design. Our findings show that the generated microchannels outperform conventional pillar arrays, offering valuable insights for heat sink designers. Additionally, imperfections observed in the de-homogenised results serve as benchmarks for future improvements, addressing concerns related to modelling accuracy, manufacturability, and overall performance enhancements.

**Keywords:** homogenisation method, topology optimisation, de-homogenisation, multi-scale design, microchannel, heat sink
# Multi-scale topology optimisation of microchannel cooling using a homogenisation-based method

Hao Li¹*, Pierre Jolivet² and Joe Alexandersen¹

¹*Institute of Mechanical and Electrical Engineering, University of Southern Denmark, 5230, Odense, Denmark.

² Sorbonne Université, CNRS, LIP6, 75252, Paris, France.

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

Contributing authors: pierre@joliv.et; joal@sdu.dk;

## Abstract

Microchannel cooling is often the preferred choice for compact heat sinks. However, widely adopted topology optimisation (TO) techniques such as density-based and level-set methods often struggle to generate very thin channel strips unless maximum length scale constraints are imposed and very fine meshes are employed. To address this limitation, multi-scale design methodologies have emerged. This paper builds upon recent advances in de-homogenisation techniques to contribute to the multi-scale design of microchannels for cooling applications. We start by selecting a single-class microstructure and employ numerical homogenisation to build an offline library. This library is then fed in online macro-scale topology optimisation, where both microstructure parameters and local orientation fields are optimised. By using a sawtooth-function-based mapping, the de-homogenised results capture fine details across different length scales through a unique homogenised design. Our findings show that the generated microchannels outperform conventional pillar arrays, offering valuable insights for heat sink designers. Additionally, imperfections observed in the de-homogenised results serve as benchmarks for future improvements, addressing concerns related to modelling accuracy, manufacturability, and overall performance enhancements.

**Keywords:** homogenisation method, topology optimisation, de-homogenisation, multi-scale design, microchannel, heat sink
# 1 Introduction

Topology optimisation (TO) is a promising numerical method that integrates sim-
ulations with mathematical optimisation algorithms to produce conceptual designs
offering very high degree of design freedom. Originating from structural optimisa-
tion, TO technique draws its roots from the seminal work of Bendsøe and Kikuchi
(1988). They introduced the concept of homogenisation, which involves represent-
ing a heterogeneous material as an equivalent homogeneous material with effective
properties. However, it was considered rather complex and always resulted in solu-
tions that could not be clearly interpreted. Therefore, after the simpler approaches
such as density-based approach (Zhou and Rozvany, 1991; Rozvany et al., 1992;
Bendsøe and Sigmund, 1999) and level-set method (Wang et al., 2003; Allaire et al.,
2004) were popularised, the homogenisation-based design approach was somewhat
lost.

On the other hand, through decades of development, TO techniques have been
applied to many other complicated multiphysics problems, one of which is the
fluid-related problems. Since the seminal work by Borrvall and Petersson (2003),
extensive research efforts have been dedicated to improving the hydrodynamic
performance by using the above-mentioned popular TO techniques (Alexander-
sen and Andreasen, 2020). In the meantime, motivated by the growing trend
towards miniaturisation in thermal-fluid devices, microchannels are often consid-
ered as the first choice for the compact heat sinks (Adham et al., 2013) and heat
exchangers (Mohammed et al., 2011). However, the authors’ experience suggests
that these popular TO techniques often fail to generate very thin channel strips
(microchannels) unless a maximum length scale constraint is imposed. Further-
more, high-resolution TO requires using very fine meshes since the design variable
(i.e., pseudo density) is defined for each pixel to represent fluid/solid phases. Con-
sequently, the forward problem becomes rather computationally intensive, and
solving the inverse design problem can be even more challenging. Therefore, multi-
scale topology optimisation techniques may be necessary to mitigate these issues.
For more comprehensive insights into multi-scale design methodologies, readers are
referred in the review papers by Wu et al. (2021) and Lee et al. (2024).

Conventional concurrent multi-scale topology optimisation is known to be computa-
tionally intensive because it requires online homogenisation of evolving microstruc-
tures at each optimisation iteration. In stark contrast to this, de-homogenisation
techniques have emerged as a promising subfield within multi-scale TO, since the
pioneering work by Pantz and Trabelsi (2008). The main idea is to reconstruct
the optimised microstructure at a specified length scale to achieve a thorough and
detailed picture of the microstructure after obtaining an optimised composite struc-
ture. During recent years, numerous research efforts have focused on developing
more efficient methods that preserve homogenised properties as well as ensure good
manufacturability (Allaire et al., 2019a,b; Groen et al., 2020).

According to the guidelines in the review paper by Lee et al. (2024), multi-scale
design methods, in a broad sense, can be further categorised into top-down and
bottom-up strategies. The top-down strategy involves first identifying the optimal
distribution of properties within the design space. Subsequently, the database of
unit cells is retrieved, and they are assembled together in a process called tiling
# 1 Introduction

Topology optimisation (TO) is a promising numerical method that integrates sim-
ulations with mathematical optimisation algorithms to produce conceptual designs
offering very high degree of design freedom. Originating from structural optimisa-
tion, TO technique draws its roots from the seminal work of Bendsøe and Kikuchi
(1988). They introduced the concept of homogenisation, which involves represent-
ing a heterogeneous material as an equivalent homogeneous material with effective
properties. However, it was considered rather complex and always resulted in solu-
tions that could not be clearly interpreted. Therefore, after the simpler approaches
such as density-based approach (Zhou and Rozvany, 1991; Rozvany et al., 1992;
Bendsøe and Sigmund, 1999) and level-set method (Wang et al., 2003; Allaire et al.,
2004) were popularised, the homogenisation-based design approach was somewhat
lost.

On the other hand, through decades of development, TO techniques have been
applied to many other complicated multiphysics problems, one of which is the
fluid-related problems. Since the seminal work by Borrvall and Petersson (2003),
extensive research efforts have been dedicated to improving the hydrodynamic
performance by using the above-mentioned popular TO techniques (Alexandersen
and Andreasen, 2020). In the meantime, motivated by the growing trend towards
miniaturisation in thermal-fluid devices, microchannels are often considered as the
first choice for the compact heat sinks (Adham et al., 2013) and heat exchang-
ers (Mohammed et al., 2011). However, the authors’ experience suggests that these
popular TO techniques often fail to generate small-scale features unless a maxi-
mum length scale constraint is imposed and very fine meshes are employed. For
example, Wu et al. (2017) used a maximum length scale constraint in the context
of structural optimisation in solid mechanics, specifically for the mean compliance
problem. They proposed a geometrical constraint on localised material volume to
generate feature-rich and bone-like structures. Similar research was conducted by
the authors (Li et al., 2023) by adopting adaptive mesh refinement to reduce the
computational cost but still required millions of mesh elements. This challenge
appears not only in structural optimisation problems but also in other multiphysics
problems, where both high-resolution meshes and proper physics representation are
crucial in encouraging fine-scale features to generate. Consequently, the forward
problem becomes rather computationally intensive, and solving the inverse design
problem can be even more challenging. Therefore, multi-scale topology optimisa-
tion techniques may be necessary to mitigate these issues. For more comprehensive
insights into multi-scale design methodologies, readers are referred in the review
papers by Wu et al. (2021) and Lee et al. (2024).

Conventional concurrent multi-scale topology optimisation is known to be computa-
tionally intensive because it requires online homogenisation of evolving microstruc-
tures at each optimisation iteration. In stark contrast to this, de-homogenisation
techniques have emerged as a promising subfield within multi-scale TO, since the
pioneering work by Pantz and Trabelsi (2008). The main idea is to reconstruct
the optimised microstructure at a specified length scale to achieve a thorough and
detailed picture of the microstructure after obtaining an optimised composite struc-
ture. During recent years, numerous research efforts have focused on developing
more efficient methods that preserve homogenised properties as well as ensure good
manufacturability (Allaire et al., 2019a,b; Groen et al., 2020).
to ensure compatibility. The bottom-up approach begins with a dataset comprising parametrised microstructures, which are then fed into the optimisation process at the macroscopic level. During this process, it is essential to ensure compatibility in geometric property mapping and, ultimately, manufacturing feasibility. This approach is more practical when dealing with problems featuring higher dimensional property space.

To date, numerous research works have successfully applied de-homogenisation techniques in the design of compliant mechanical systems including projection-based method (Geoffroy-Donders et al., 2020; Lee et al., 2021; Jensen et al., 2022), streamline-based approach (Wang et al., 2023), convolutional neural networks assisted method (Elingaard et al., 2022), and phasor-noise driven technique (Woldseth et al., 2024), among others. However, for many other more complicated multiphysics problems, de-homogenisation is still in its infancy stages.

In the context of multi-scale design of fluid-related problems, we summarise in Table 1 the state-of-the-art works, clarifying the design strategies, the selected optimisers, and the choice of microstructures in the literature. Hence, we clarify the positioning of the present work. Additionally, the categorisation of these representative works is illustrated in Fig. 1.

<table><caption>**Table 1:** The state-of-the-art works on the multi-scale design of fluid-related problems (listed in chronological order).</caption><thead><tr><th>Ref.</th><th>Scale</th><th>Design strategies</th><th>Optimiser</th><th>Unit-cell orientation</th><th>Microstructure</th><th>2D/3D</th><th>Application</th></tr></thead><tbody><tr><td>Guest and Prévost (2007)</td><td>Micro-scale</td><td></td><td>MMA</td><td>No</td><td>Class-free</td><td>2D/3D</td><td>Maximal permeability</td></tr><tr><td>Wu (2019)</td><td>Multi-scale</td><td>Concurrent</td><td>GCMMA</td><td>Yes</td><td>Class-free</td><td>2D</td><td>Minimal dissipation</td></tr><tr><td>Takezawa et al. (2019b,a)</td><td>Multi-scale</td><td>Bottom-up</td><td>MMA</td><td>No</td><td>Single-class</td><td>3D</td><td>Heat sink</td></tr><tr><td>Dede et al. (2020)</td><td>Multi-scale</td><td>Bottom-up</td><td>MMA</td><td>Yes</td><td>Single-class</td><td>2D</td><td>Flow distribution</td></tr><tr><td>Ozguc et al. (2021)</td><td>Multi-scale</td><td>Bottom-up</td><td>MMA</td><td></td><td></td><td>2D</td><td>Heat sink</td></tr><tr><td>Geng et al. (2022)</td><td>Multi-scale</td><td>Bottom-up</td><td>MMA</td><td>No</td><td>Single-class</td><td>2D</td><td>Heat sink</td></tr><tr><td>Zhou et al. (2022)</td><td>Multi-scale</td><td>Bottom-up</td><td>MMA</td><td>Yes</td><td>Single-class</td><td>2D</td><td>Microreactor</td></tr><tr><td>Hankins et al. (2023)</td><td>Multi-scale</td><td>Bottom-up</td><td>MMA</td><td>Yes</td><td>Single-class</td><td>2D</td><td>Microreactor</td></tr><tr><td>Padhy et al. (2023a)</td><td>Multi-scale</td><td>Bottom-up</td><td>Neural networks</td><td>Yes</td><td>Diversified-class</td><td>2D</td><td>Minimal dissipation</td></tr><tr><td>Padhy et al. (2023b)</td><td>Multi-scale</td><td>Bottom-up</td><td>Neural networks</td><td>Yes</td><td>Super-shaped</td><td>2D</td><td>Minimal dissipation</td></tr><tr><td>Feppon (2024b)</td><td>Multi-scale</td><td>Bottom-up</td><td>Null-space optimiser</td><td>Yes</td><td>Single-class</td><td>2D</td><td>Flow uniformity</td></tr><tr><td>This work</td><td>Multi-scale</td><td>Bottom-up</td><td>MMA</td><td>Yes</td><td>Single-class</td><td>2D</td><td>Heat sink</td></tr></tbody></table>

Guest and Prévost (2007) pioneered the idea in designing periodic porous material microstructures that maximise permeability and achieve prescribed flow symmetries in bulk materials. They computed permeability through numerical homogenisation of the base cell using finite elements, but their work was limited to the micro-scale design. Wu (2019) is known as the first who conducted the true sense of multi-scale TO of fluid mechanics problems, utilising a concurrent optimisation technique to design optimised micro pillar arrays, aiming to minimize energy loss due to local fluid drag forces. However, as mentioned earlier, the concurrent TO necessitates on-the-fly homogenisation, which requires expensive computational cost.

Ozguc et al. (2021) developed a homogenisation-based method where the so-called partial densities are physically represented as porous microstructures, enabling the design of thermofluidic structures with sub-grid features represented by "grey regions" but without providing the de-homogenised cooling channels.

Takezawa et al. (2019a,b) were among the pioneers in introducing the bottom-up design strategy to multi-scale fluid-flow designs. They employed a typical lattice structure comprising pillars, with the lattice shape parametrised by pillar diameter.
According to the guidelines in the review paper by Lee *et al.* (2024), multi-scale design methods, in a broad sense, can be further categorised into top-down and bottom-up strategies. The top-down strategy involves first identifying the optimal distribution of properties within the design space. Subsequently, the database of unit cells is retrieved, and they are assembled together in a process called tiling to ensure compatibility. The bottom-up approach begins with a dataset comprising parametrised microstructures, which are then fed into the optimisation process at the macroscopic level. During this process, it is essential to ensure compatibility in geometric property mapping and, ultimately, manufacturing feasibility. This approach is more practical when dealing with problems featuring higher dimensional property space.

To date, numerous research works have successfully applied de-homogenisation techniques in the design of compliant mechanical systems including projection-based method (Geoffroy-Donders *et al.*, 2020; Lee *et al.*, 2021; Jensen *et al.*, 2022; Liu *et al.*, 2023; Chen *et al.*, 2023a,b), streamline-based approach (Wang *et al.*, 2023), convolutional neural networks assisted method (Elingaard *et al.*, 2022), and phasor-noise driven technique (Woldseth *et al.*, 2024), among others. However, for many other more complicated multiphysics problems, de-homogenisation is still in its infancy stages.

In the context of multi-scale design of fluid-related problems, we summarise in Table 1 the state-of-the-art works, clarifying the design strategies, the selected optimisers, and the choice of microstructures in the literature. Hence, we clarify the positioning of the present work. Additionally, the categorisation of these representative works is illustrated in Fig. 1.

<table><caption>**Table 1:** The state-of-the-art works on the multi-scale design of fluid-related problems (listed in chronological order).</caption><thead><tr><th>Ref.</th><th>Scale</th><th>Design strategies</th><th>Optimiser</th><th>Unit-cell orientation</th><th>Microstructure</th><th>2D/3D</th><th>Application</th></tr></thead><tbody><tr><td>Guest and Prévost (2007)</td><td>Micro-scale</td><td>-</td><td>MMA</td><td>No</td><td>Class-free</td><td>2D/3D</td><td>Maximal permeability</td></tr><tr><td>Wu (2019)</td><td>Multi-scale</td><td>Concurrent</td><td>GCMMMA</td><td>Yes</td><td>Class-free</td><td>2D</td><td>Minimal dissipation</td></tr><tr><td>Takezawa et al. (2019a)</td><td>Multi-scale</td><td>Bottom-up</td><td>MMA</td><td>No</td><td>Single-class</td><td>3D</td><td>Heat sink</td></tr><tr><td>Dede et al. (2020)</td><td>Multi-scale</td><td>Bottom-up</td><td>MMA</td><td>Yes</td><td>Single-class</td><td>2D</td><td>Flow distribution</td></tr><tr><td>Ozguc et al. (2021)</td><td>Multi-scale</td><td>Bottom-up</td><td>MMA</td><td>-</td><td>-</td><td>2D</td><td>Heat sink</td></tr><tr><td>Geng et al. (2022)</td><td>Multi-scale</td><td>Bottom-up</td><td>MMA</td><td>No</td><td>Single-class</td><td>2D</td><td>Heat sink</td></tr><tr><td>Zhou et al. (2022)</td><td>Multi-scale</td><td>Bottom-up</td><td>MMA</td><td>Yes</td><td>Single-class</td><td>2D</td><td>Microreactor</td></tr><tr><td>Hankins et al. (2023)</td><td>Multi-scale</td><td>Bottom-up</td><td>MMA</td><td>Yes</td><td>Single-class</td><td>2D</td><td>Microreactor</td></tr><tr><td>Padhy et al. (2023a)</td><td>Multi-scale</td><td>Bottom-up</td><td>Neural networks</td><td>Yes</td><td>Diversified-class</td><td>2D</td><td>Minimal dissipation</td></tr><tr><td>Padhy et al. (2023b)</td><td>Multi-scale</td><td>Bottom-up</td><td>Neural networks</td><td>Yes</td><td>Super-shaped</td><td>2D</td><td>Minimal dissipation</td></tr><tr><td>Feppon (2021b)</td><td>Multi-scale</td><td>Bottom-up</td><td>Null-space optimiser</td><td>Yes</td><td>Single-class</td><td>2D</td><td>Flow uniformity</td></tr><tr><td>This work</td><td>Multi-scale</td><td>Bottom-up</td><td>MMA</td><td>Yes</td><td>Single-class</td><td>2D</td><td>Heat sink</td></tr></tbody></table>

Guest and Prévost (2007) pioneered the idea in designing periodic porous material microstructures that maximise permeability and achieve prescribed flow symmetries in bulk materials. They computed permeability through numerical homogenisation of the base cell using finite elements, but their work was limited to the micro-scale design. Wu (2019) is known as the first who conducted the true sense of multi-scale TO of fluid mechanics problems, utilising a concurrent optimisation technique to design optimised micro pillar arrays, aiming to minimize energy loss due to local fluid drag forces. However, as mentioned earlier, the concurrent TO necessitates on-the-fly homogenisation, which requires expensive computational cost.

Ozguc *et al.* (2021) developed a homogenisation-based method where the so-called partial densities are physically represented as porous microstructures, enabling
<figure><img src="image_2.png" /><figcaption>**Fig. 1:** Categorisation of the 11 representative existing works summarised in Table 1 and their correlations.</figcaption></figure>

Fluid flow was approximated by deriving effective properties from the Darcy–Forchheimer law, and at the macro-scale, the Brinkman–Forchheimer equation was applied. The lattice density, namely, the pillar diameter was carefully designed for coolant flow. A similar design workflow can be found in the work by Geng et al. Geng et al. (2022) who employed a lattice structure consisting of crosses. Numerical homogenisation was conducted to calculate effective permeability and thermal conductivity, followed by the construction of a surrogate model. They successfully demonstrated concurrent topology optimisation for multi-scale fluid channels, inlets, and outlets. However, the inability to rotate the unit-cell limits the degree of design freedom to some extent.

Alternatively, in the work by Dede et al. (2020), a reaction-diffusion equation driven de-homogenisation approach was proposed to build Turing pattern microchannel flow structures that fill space based on the porous media field obtained by a gradient-based optimiser. They showcased precise control over fluid flow distribution across numerous microchannels, typically numbering in the hundreds. In their subsequent works, they integrated an anisotropic porous medium rotated by the orientation tensor and showcased the effectiveness of their de-homogenisation framework in microreactor design (Zhou et al., 2022). Later in their work, they presented design cases with zones identified based on different physical objectives and made a record of 200 unique and distinctly different microreactor flow channels (Hankins et al., 2023). They also presented an experimental investigation of airflow from a Turing pattern microchannel array, with channel sizes ranging from 0.6 to 1.5 mm, emanated into an open atmosphere at very low Reynolds numbers (Dede et al., 2022).

More recently, motivated by the need for fluid-flow devices with large contact areas, Padhy et al. (2023a) introduced a neural networks-assisted bottom-up design framework. In this approach, various pre-selected classes of microstructures are
<figure><img src="image_3.png" /><figcaption>**Fig. 1:** Categorisation of the 11 representative existing works summarised in Table 1 and their correlations.</figcaption></figure>

the design of thermofluidic structures with sub-grid features represented by "grey regions" but without providing the de-homogenised cooling channels.

Takezawa et al. (2019a,b) were among the pioneers in introducing the bottom-up design strategy to multi-scale fluid-flow designs. They employed a typical lattice structure comprising pillars, with the lattice shape parametrised by pillar diameter. Fluid flow was approximated by deriving effective properties from the Darcy–Forchheimer law, and at the macro-scale, the Brinkman–Forchheimer equation was applied. The lattice density, namely, the pillar diameter was carefully designed for coolant flow. A similar design workflow can be found in the work by Geng et al. Geng et al. (2022) who employed a lattice structure consisting of crosses. Numerical homogenisation was conducted to calculate effective permeability and thermal conductivity, followed by the construction of a surrogate model. They successfully demonstrated concurrent topology optimisation for multi-scale fluid channels, inlets, and outlets. However, the inability to rotate the unit-cell limits the degree of design freedom to some extent.

Alternatively, in the work by Dede et al. (2020), a reaction-diffusion equation driven de-homogenisation approach was proposed to build Turing pattern microchannel flow structures that fill space based on the porous media field obtained by a gradient-based optimiser. They showcased precise control over fluid flow distribution across numerous microchannels, typically numbering in the hundreds. In their subsequent works, they integrated an anisotropic porous medium rotated by the orientation tensor and showcased the effectiveness of their de-homogenisation framework in microreactor design (Zhou et al., 2022). Later in their work, they presented design cases with zones identified based on different physical objectives and made a record of 200 unique and distinctly different microreactor flow channels (Hankins et al., 2023). They also presented an experimental investigation of airflow from a Turing pattern microchannel array, with channel sizes ranging
parametrised to construct an offline library using numerical homogenisation. During macro-scale topology optimisation, their algorithm permits seamless transition between different classes. Leveraging the automatic differentiation capability provided by the neural networks, manual sensitivity analysis can be eliminated. In a recent preprint (Padhy et al., 2023b), they further investigated the potential of using super-shapes as microstructures, enabling the generation of new microstructures not present in the pre-built offline library. Most recently, Feppon (2024b) derived Darcy's law asymptotically, often referred to as a “poor man’s approach” in the fluid TO community (Zhao et al., 2018; Pollini et al., 2020; Geng et al., 2022), for a periodically porous medium deformed by a diffeomorphism. He demonstrated that the homogenised permeability matrix depends not only on the local orientation but also on the local dilation of the deformed periodic medium. He built upon the existing de-homogenisation technique from compliant mechanical systems to the design of flow uniformity. Full-scale simulation results showed that the proposed designs achieved satisfactory flow uniformity but exhibited noticeable deviations from the homogenized model.

As mentioned above, multi-scale design of fluid-related problems are relatively recent topics that have been explored in only a few research works. Their potential application to microchannel cooling designs remains largely unexplored. The main contribution of this paper is to showcase that a bottom-up homogenisation-based TO and de-homogenisation techniques can be effectively applied to the context of cooling channel design considering conjugate heat transfer. As for the fluid dynamics modelling, we rely on the Darcy model for simulating fluid-flow behaviour in porous media. This model involves simplifications such as the ignorance of boundary layer effects. The thermal modelling assumes a two-dimensional setup without considering the thermal interaction between base plate and the thermofluidic design layer. Regardless of those simplification, we demonstrate that the generated microchannels can outperform conventional micro-pillar arrays through full-scale simulations with upsampled structures, offering valuable insight for heat sink designers. Moreover, the imperfections observed in those de-homogenised designs such as disconnected branches and limitations in pre-selected microstructures, triggers future research prospects and offer a benchmark for comparison.

The remainder of this paper is structured as follows: in Section 2, we outline the homogenisation problem for fluid flow and conjugate heat transfer, the unit-cell problems for computing homogenized permeability and thermal conductivity, the filtering and projection schemes used, the optimisation mathematical model, and the projection-based de-homogenisation techniques. In Section 3, we provide details of the implementation, including weak formulations for the partial differential equations (PDEs) and optimisation algorithms. In Section 4, we demonstrate several numerical examples to validate the current design procedure. Finally, we draw the conclusions, outline the limitations and future work in Section 5.

## 2 Formulation

In this section, we start by presenting the homogenized problems for fluid flow and conjugate heat transfer in Section 2.1. Subsequently, we outline the unit-cell problems utilized for computing the homogenized properties, namely the effective inverse permeability and thermal conductivity, in Section 2.2. Following this, we present the filtering and projection schemes in Section 2.3. Subsequent to that, we
from 0.6 to 1.5 mm, emanated into an open atmosphere at very low Reynolds numbers (Dede et al., 2022).

More recently, motivated by the need for fluid-flow devices with large contact areas, Padhy et al. (2023a) introduced a neural networks-assisted bottom-up design framework. In this approach, various pre-selected classes of microstructures are parametrised to construct an offline library using numerical homogenisation. During macro-scale topology optimisation, their algorithm permits seamless transition between different classes. Leveraging the automatic differentiation capability provided by the neural networks, manual sensitivity analysis can be eliminated. In a recent preprint (Padhy et al., 2023b), they further investigated the potential of using super-shapes as microstructures, enabling the generation of new microstructures not present in the pre-built offline library. Most recently, Feppon (2024b) derived Darcy's law asymptotically, often referred to as a “poor man's approach” in the fluid TO community (Zhao et al., 2018; Pollini et al., 2020; Geng et al., 2022), for a periodically porous medium deformed by a diffeomorphism. He demonstrated that the homogenised permeability matrix depends not only on the local orientation but also on the local dilation of the deformed periodic medium. He built upon the existing de-homogenisation technique from compliant mechanical systems to the design of flow uniformity. Full-scale simulation results showed that the proposed designs achieved satisfactory flow uniformity but exhibited noticeable deviations from the homogenised model.

As mentioned above, multi-scale design of fluid-related problems are relatively recent topics that have been explored in only a few research works. Their potential application to microchannel cooling designs remains largely unexplored. The main contribution of this paper is to showcase that a bottom-up homogenisation-based TO and de-homogenisation techniques can be effectively applied to the context of cooling channel design considering conjugate heat transfer. As for the fluid dynamics modelling, we rely on the Darcy model for simulating fluid-flow behaviour in porous media. This model involves simplifications such as the ignorance of boundary layer effects. The thermal modelling assumes a two-dimensional setup without considering the thermal interaction between base plate and the thermofluidic design layer. Regardless of those simplification, we demonstrate that the generated microchannels can outperform conventional micro-pillar arrays through full-scale simulations with upsampled structures, offering valuable insight for heat sink designers. Moreover, the imperfections observed in those de-homogenised designs such as disconnected branches and limitations in pre-selected microstructures, triggers future research prospects and offer a benchmark for comparison.

The remainder of this paper is structured as follows: in Section 2, we outline the homogenisation problem for fluid flow and conjugate heat transfer, the unit-cell problems for computing homogenised permeability and thermal conductivity, the filtering and projection schemes used, the optimisation mathematical model, and the projection-based de-homogenisation techniques. In Section 3, we provide details of the implementation, including weak formulations for the partial differential equations (PDEs) and optimisation algorithms. In Section 4, we demonstrate several numerical examples to validate the current design procedure. Finally, we draw the conclusions, outline the limitations and future work in Section 5.
formulate the optimal design problem for a thermal-fluid system in Section 2.4. Finally, in Section 2.5, we provide an overview of the de-homogenisation procedure.

<figure><img src="image_3.png" /><figcaption>Fig. 2: Schematic diagram of the homogenisation procedure.</figcaption></figure>

## 2.1 Homogenised problem

Consider a fixed and bounded domain $\Omega_{\epsilon}$ in $\mathbb{R}^d$, where $d = 2$ or $3$, comprising periodic heterogeneous media consisting of a series of microscopic unit cells with a characteristic length $\epsilon$, as illustrated in the left sub-figure of Fig. 2. Each unit cell comprises two homogeneous and isotropic material phases: the fluid phase and the solid phase, as shown in the middle sub-figure of Fig. 2. Assuming the microstructure is a square with its width denoted by $m_1$ and height denoted by $m_2$, and the unit cell is rotated by an angle $\theta$ within the global coordinate system. The representation of such heterogeneous media can be achieved through a homogeneous media with the effective material properties (also referred to as homogenised properties) as illustrated in the right sub-figure of Fig. 2. Utilizing the homogenisation procedure enables us to potentially mitigate the considerable computational expense associated with solving the governing equations on a fine mesh, transforming it into a computationally efficient model on a relatively coarse mesh with acceptable accuracy. In other words, rather than computing the detailed behaviour of fluid velocity $u_{\epsilon}$, pressure $p_{\epsilon}$, and temperature $T_{\epsilon}$ on a fine mesh, we compute their homogenized counterparts $u$, $p$, and $T$ on a coarse mesh. This approach is particularly advantageous for the conceptual design stage of multi-scale structures.

### 2.1.1 Fluid flow

A slow fluid flow (characterized by the Reynolds number $Re \ll 1$) enters domain composed of multi-scale channels via an inlet $\Gamma_{in}$ with a prescribed velocity $u_{\epsilon} = u_0$ which exhibits a parabolic profile, and the fluid flow exits the domain through the outlet boundary $\Gamma_{out}$ with a zero normal-stress boundary condition. The full-scale
## 2 Formulation

In this section, we start by presenting the homogenised problems for fluid flow and conjugate heat transfer in Section 2.1. Subsequently, we outline the unit-cell problems utilised for computing the homogenised properties, namely the effective inverse permeability and thermal conductivity, in Section 2.2. Following this, we present the filtering and projection schemes in Section 2.3. Subsequent to that, we formulate the optimal design problem for a thermal-fluid system in Section 2.4. Finally, in Section 2.5, we provide an overview of the de-homogenisation procedure.

<figure><img src="image_4.png" /><figcaption>Fig. 2: Schematic diagram of the homogenisation procedure.</figcaption></figure>

### 2.1 Homogenised problem

Consider a fixed and bounded domain $\Omega_{\epsilon}$ in $\mathbb{R}^d$, where $d = 2$ or $3$, comprising periodic heterogeneous media consisting of a series of microscopic unit cells with a characteristic length $\epsilon$, as illustrated in the left sub-figure of Fig. 2. The domain $\Omega_{\epsilon}$ is bounded by the Dirichlet boundary $\partial\Omega_{\epsilon}^D$ and Neumann boundary $\partial\Omega_{\epsilon}^N$. Each unit cell comprises two homogeneous and isotropic material phases: the fluid phase and the solid phase, as shown in the middle sub-figure of Fig. 2. Assuming the microstructure is a square with its width denoted by $m_1$ and height denoted by $m_2$, and the unit cell is rotated by an angle $\theta$ within the global coordinate system. The representation of such heterogeneous media can be achieved through a homogeneous media with the effective material properties (also referred to as homogenised properties) as illustrated in the right sub-figure of Fig. 2. Utilising the homogenisation procedure enables us to potentially mitigate the considerable computational expense associated with solving the governing equations on a fine mesh, transforming it into a computationally efficient model on a relatively coarse mesh with acceptable accuracy. In other words, rather than computing the detailed behaviour of fluid velocity $u_{\epsilon}$, pressure $p_{\epsilon}$, and temperature $T_{\epsilon}$ on a fine mesh, we compute their homogenised counterparts $u, p$, and $T$ on a coarse mesh. This approach is particularly advantageous for the conceptual design stage of multi-scale structures.
fluid-flow velocity $\mathbf{u}_\epsilon$ and pressure $p_\epsilon$ are governed by the Stokes equations in $\Omega_\epsilon$, which are as follows:

$$ egin{cases} -\Delta \mathbf{u}_\epsilon + \nabla p_\epsilon = -\alpha_\epsilon(\mathbf{x}) \mathbf{u}_\epsilon & \text{in } \Omega_\epsilon, \\ \operatorname{div}(\mathbf{u}_\epsilon) = 0 & \text{in } \Omega_\epsilon, \\ \mathbf{u}_\epsilon = \mathbf{u}_0 & \text{on } \Gamma_{\text{in}}, \\ (-p_\epsilon \mathbf{I} + (\nabla \mathbf{u}_\epsilon + \nabla \mathbf{u}_\epsilon^T)) \cdot \mathbf{n} = 0 & \text{on } \Gamma_{\text{out}}, \\ \mathbf{u}_\epsilon = 0 & \text{on } \partial\Omega_\epsilon \setminus (\Gamma_{\text{in}} \cup \Gamma_{\text{out}}), \end{cases} \quad (1) $$

where a fictitious body-force term $-\alpha_\epsilon(\mathbf{x}) \mathbf{u}_\epsilon$ is introduced to represent the presence or absence of solid phase, where the inverse permeability $\alpha_\epsilon(\mathbf{x})$ is varied with space defined as follows:

$$ \alpha_{\epsilon}(\mathbf{x}) := \begin{cases} 0 & \text{if } \mathbf{x} \in \Omega_{\epsilon,f}, \\ \alpha_{\max} & \text{if } \mathbf{x} \in \Omega_{\epsilon,s}, \end{cases} \quad (2) $$

where $\alpha_{\max}$ represents the maximum effective permeability. It is a positive value, sufficiently large to prevent fluid flow penetration and approximate solid regions. We assume that the solution for fluid velocity and pressure, denoted as $(\mathbf{u}_\epsilon, p_\epsilon)$ in the Stokes equation shown in Eq. (1) may be approximated by the effective (or homogenised) velocity and pressure $(\mathbf{u}, p)$ to the Darcy potential flow model. The mathematical proof and numerical validations supporting this asymptotic expansion are detailed in the recent preprint (Feppon, 2024a). Thus, the homogenized problem can be characterized as follows:

$$ egin{cases} \mathbf{u} = -\epsilon^2 \mathbf{X}^* (m_1, m_2, \theta) \nabla p & \text{in } D, \\ \operatorname{div}(\mathbf{u}) = 0 & \text{in } D, \\ \mathbf{u} \cdot \mathbf{n} = \mathbf{u}_0 \cdot \mathbf{n} & \text{on } \Gamma_{\text{in}}, \\ p = 0 & \text{on } \Gamma_{\text{out}}, \end{cases} \quad (3) $$

where the normalization factor $\epsilon^2$ is multiplied to render the homogenisation equation approximately independent of $\epsilon$. $\mathbf{X}^* (m_1, m_2, \theta)$ is the global permeability matrix, which can be expressed as follows:

$$ \mathbf{X}^* (m_1, m_2, \theta) = \mathbf{R}(\theta) \mathbf{X}^H (m_1, m_2) \mathbf{R}(\theta)^T, \quad (4) $$

where $\mathbf{R}$ is the well-known rotation matrix defined as follows:

$$ \mathbf{R}(\theta) = \begin{bmatrix} \cos(\theta) & -\sin(\theta) \\ \sin(\theta) & \cos(\theta) \end{bmatrix}. \quad (5) $$

$\mathbf{X}^H$ represents the permeability tensor of the undeformed unit cell $Y$, defined as:

$$ \mathbf{X}^H = \left[ \begin{array}{c|c} \mathcal{X}_{11} & \mathcal{X}_{12} \\ \hline \mathcal{X}_{21} & \mathcal{X}_{22} \end{array} \right] = \frac{1}{|Y|} \left[ \begin{array}{c|c} \int_Y \mathcal{X}_{1,y_1} \, dy_1 & \int_Y \mathcal{X}_{2,y_1} \, dy_1 \\ \hline \int_Y \mathcal{X}_{1,y_2} \, dy_2 & \int_Y \mathcal{X}_{2,y_2} \, dy_2 \end{array} \right], \quad (6) $$

where $\mathbf{X}_1 = [\chi_{1,y_1}, \chi_{1,y_2}]^T$ and $\mathbf{X}_2 = [\chi_{2,y_1}, \chi_{2,y_2}]^T$ are the solution to the unit-cell problem which will be given in Section 2.2.1. Moreover, due to the symmetric nature of the microstructure employed in this paper, $\mathbf{X}_{12} = \mathbf{X}_{21} = 0$.
### 2.1.1 Fluid flow

A slow fluid flow (characterised by the Reynolds number $Re \ll 1$) enters domain composed of multi-scale channels via an inlet $\Gamma_{in}$ with a prescribed velocity $u_{\epsilon} = u_{0}$ which exhibits a parabolic profile, and the fluid flow exits the domain through the outlet boundary $\Gamma_{out}$ with a zero normal-stress boundary condition. The full-scale non-dimensional fluid-flow velocity $u_{\epsilon}$ and pressure $p_{\epsilon}$ are governed by the Stokes equations in $\Omega_{\epsilon}$, which are as follows:

$$ egin{cases} 
-\Delta u_{\epsilon} + \nabla p_{\epsilon} = -\alpha_{\epsilon}(x) u_{\epsilon} & \text{in } \Omega_{\epsilon}, \\ 
\operatorname{div}(u_{\epsilon}) = 0 & \text{in } \Omega_{\epsilon}, \\ 
u_{\epsilon} = u_{0} & \text{on } \Gamma_{\text{in}}, \\ 
(-p_{\epsilon}\mathbf{I} + (\nabla u_{\epsilon} + \nabla u_{\epsilon}^{T})) \cdot \mathbf{n} = 0 & \text{on } \Gamma_{\text{out}}, \\ 
u_{\epsilon} = 0 & \text{on } \partial\Omega_{\epsilon} \setminus (\Gamma_{\text{in}} \cup \Gamma_{\text{out}}), 
\end{cases} 
\quad (1) $$

where a fictitious body-force term $-\alpha_{\epsilon}(x) u_{\epsilon}$ is introduced to represent the presence or absence of solid phase, where the inverse permeability $\alpha_{\epsilon}(x)$ is varied with space defined as follows:

$$ \alpha_{\epsilon}(\mathbf{x}) := \begin{cases} 0 & \text{if } \mathbf{x} \in \Omega_{\epsilon,f}, \\ \alpha_{\max} & \text{if } \mathbf{x} \in \Omega_{\epsilon,s}, \end{cases} \quad (2) $$

where $\alpha_{\max}$ represents the maximum effective permeability. It is a positive value, sufficiently large to prevent fluid flow penetration and approximate solid regions. The homogenisation of Stokes equations in porous media has been studied extensively since the late 1970s. Notably, Sánchez-Palencia (1980b), Sanchez-Palencia (1982), and Lévy (1983) demonstrated that Darcy's law can be derived as the limit of the Stokes equations in a periodic porous medium as the periodicity tends to zero. The mathematical foundations and numerical validations supporting this asymptotic expansion are presented in works such as Allaire (1991a,b) and Feppon (2024a). Following the lead of these works, we make assumption that the solution for fluid velocity and pressure, denoted as $(u_{\epsilon}, p_{\epsilon})$ in the Stokes equation shown in Eq. (1) may be approximated by the effective (or homogenised) velocity and pressure $(u, p)$ to the Darcy potential flow model. Thus, the homogenised problem can be characterised as follows:

$$ egin{cases} 
u = -\epsilon^2 \mathbf{X}^* (m_1, m_2, \theta) \nabla p & \text{in } D, \\ 
\operatorname{div}(u) = 0 & \text{in } D, \\ 
u \cdot \mathbf{n} = u_0 \cdot \mathbf{n} & \text{on } \Gamma_{\text{in}}, \\ 
p = 0 & \text{on } \Gamma_{\text{out}}, 
\end{cases} 
\quad (3) $$

where $\mathbf{X}^* (m_1, m_2, \theta)$ is the global effective permeability matrix, which is scale-invariant. Therefore, it is multiplied by a scaling factor $\epsilon^2$ to account for the effect of unit cell size (Auriault et al., 2010). $\mathbf{X}^* (m_1, m_2, \theta)$ can be expressed as follows:

$$ \mathbf{X}^* (m_1, m_2, \theta) = \mathbf{R}(\theta) \mathbf{X}^H (m_1, m_2) \mathbf{R}(\theta)^T, \quad (4) $$

where $\mathbf{R}$ is the well-known rotation matrix defined as follows:

$$ \mathbf{R}(\theta) = \begin{bmatrix} \cos(\theta) & -\sin(\theta) \\ \sin(\theta) & \cos(\theta) \end{bmatrix}. \quad (5) $$
Note that the homogenised problem given in Eq. (3) is rewritten in the form used in practise as an equivalent Laplacian equation, expressed as follows:

$$ \left\{ \begin{array}{ll} -\operatorname{div}(\boldsymbol{X}^*(m_1, m_2, \theta) \nabla p) = 0 & \text{in } D, \\ \boldsymbol{n}^T \boldsymbol{X}^*(m_1, m_2, \theta) \nabla p = -\epsilon^{-2}\boldsymbol{u}_0 \cdot \boldsymbol{n} & \text{on } \Gamma_{\text{in}}, \\ p = 0 & \text{on } \Gamma_{\text{out}}. \end{array} \right. \qquad (7) $$

## 2.1.2 Heat transfer

A cold fluid enters the domain via an inlet $\Gamma_{\text{in}}$ with a prescribed temperature $T_{\epsilon} = T_0$, and the remaining boundaries hold adiabatic boundary conditions. With the knowledge of velocity $\boldsymbol{u}_{\epsilon}$ as the solution of Eq. (1), it is then weakly coupled to the convection-diffusion equation in $\Omega_{\epsilon}$, as follows:

$$ \left\{ \begin{array}{ll} \boldsymbol{u}_{\epsilon} \cdot \nabla T_{\epsilon} - \operatorname{div}(\mathcal{K}_{\epsilon} \nabla T_{\epsilon}) = Q(\boldsymbol{x}) & \text{in } \Omega_{\epsilon}, \\ T_{\epsilon} = T_0 & \text{on } \Gamma_{\text{in}}, \\ \mathcal{K}_{\epsilon} \nabla T_{\epsilon} \cdot \boldsymbol{n} = 0 & \text{on } \partial\Omega_{\epsilon} \setminus \Gamma_{\text{in}}, \end{array} \right. \qquad (8) $$

where $\boldsymbol{u}_{\epsilon}$ is the solution to the Stokes equations Eq. (1), and the thermal conductivity and heat source is given by:

$$ \mathcal{K}_{\epsilon}(\boldsymbol{x}) = \begin{cases} \kappa_{\text{f}} & \text{if } \boldsymbol{x} \in \Omega_{\epsilon, \text{f}}, \\ \kappa_{\text{s}} & \text{if } \boldsymbol{x} \in \Omega_{\epsilon, \text{s}}, \end{cases} \quad Q(\boldsymbol{x}) = \begin{cases} 0 & \text{if } \boldsymbol{x} \in \Omega_{\epsilon, \text{f}}, \\ \dot{Q} & \text{if } \boldsymbol{x} \in \Omega_{\epsilon, \text{s}}, \end{cases} \qquad (9) $$

with the thermal conductivity ratio between fluid and solid phases $\kappa_{\text{f}}/\kappa_{\text{s}} = 1/50$. The volumetric heat source is assumed to be design-dependent with a heat generation rate $\dot{Q} = 1.0$. Applying a volumetric heat source exclusively to the solid phase in a two-dimensional conjugate heat transfer problem is justified by the observation that, in numerous practical situations like electronic devices or solid structures subjected to external heating, heat generation predominantly occurs within the solid material. Consequently, disregarding heat generation in the fluid phase is reasonable, as the primary function of the fluid is to facilitate heat transport rather than generate it. Then, the solution for temperature $T_{\epsilon}$ in the full-scale convection-diffusion equation can be approximated by the effective (or homogenised) temperature $T$ obtained by solving the following the homogenised problem (Allaire and Habibi, 2013):

$$ \left\{ \begin{array}{ll} \boldsymbol{u} \cdot \nabla T - \operatorname{div}(\boldsymbol{K}^*(m_1, m_2, \theta) \nabla T) = Q(m_1, m_2) & \text{in } D, \\ T = T_0 & \text{on } \Gamma_{\text{in}}, \\ \boldsymbol{K}^*(m_1, m_2, \theta) \nabla T \cdot \boldsymbol{n} = 0 & \text{on } \partial D \setminus \Gamma_{\text{in}}, \end{array} \right. \qquad (10) $$

where $\boldsymbol{K}^*(m_1, m_2, \theta)$ is the global thermal conductivity matrix, which can be expressed as follows:

$$ \boldsymbol{K}^*(m_1, m_2, \theta) = \boldsymbol{R}(\theta) \boldsymbol{K}^H(m_1, m_2) \boldsymbol{R}(\theta)^T, \qquad (11) $$
$\boldsymbol{\chi}^H$ represents the permeability tensor of the undeformed unit cell $Y$, defined as:

$$ \boldsymbol{\chi}^H = \left[ \frac{\chi_{11}^{\frac{1}{2}} \chi_{12}^{\frac{1}{2}}}{\chi_{21}^{\frac{1}{2}} \chi_{22}^{\frac{1}{2}}} \right] = \frac{1}{|Y|} \left[ \frac{\int_Y \chi_{1,y_1} dy_1}{\int_Y \chi_{1,y_2} dy_2} \frac{\int_Y \chi_{2,y_1} dy_1}{\int_Y \chi_{2,y_2} dy_2} \right], \quad (6) $$

where $\boldsymbol{\chi}_1 = [\chi_{1,y_1}, \chi_{1,y_2}]^T$ and $\boldsymbol{\chi}_2 = [\chi_{2,y_1}, \chi_{2,y_2}]^T$ are the solution to the unit-cell problem which will be given in Section 2.2.1. Moreover, due to the symmetric nature of the microstructure employed in this paper, $\chi_{12} = \chi_{21} = 0$.

Note that the homogenised problem given in Eq. (3) is rewritten in the form used in practise as an equivalent Laplacian equation, expressed as follows:

$$ \left\{ \begin{array}{ll} -\operatorname{div}(\boldsymbol{\chi}^* (m_1, m_2, \theta) \nabla p) = 0 & \text{in } D, \\ \boldsymbol{n}^T \boldsymbol{\chi}^* (m_1, m_2, \theta) \nabla p = -\epsilon^{-2} \boldsymbol{u}_0 \cdot \boldsymbol{n} & \text{on } \Gamma_{\text{in}}, \\ p = 0 & \text{on } \Gamma_{\text{out}}. \end{array} \right. \quad (7) $$

### 2.1.2 Heat transfer

A cold fluid enters the domain via an inlet $\Gamma_{\text{in}}$ with a prescribed non-dimensional temperature $T_{\epsilon} = T_0$, and the remaining boundaries hold adiabatic boundary conditions. With the knowledge of non-dimensional velocity $\boldsymbol{u}_{\epsilon}$ as the solution of Eq. (1), it is then weakly coupled to the convection-diffusion equation in $\Omega_{\epsilon}$, as follows:

$$ \left\{ \begin{array}{ll} \boldsymbol{u}_{\epsilon} \cdot \nabla T_{\epsilon} - \operatorname{div}(\mathcal{K}_{\epsilon} \nabla T_{\epsilon}) = Q(\boldsymbol{x}) & \text{in } \Omega_{\epsilon}, \\ T_{\epsilon} = T_0 & \text{on } \Gamma_{\text{in}}, \\ \mathcal{K}_{\epsilon} \nabla T_{\epsilon} \cdot \boldsymbol{n} = 0 & \text{on } \partial\Omega_{\epsilon} \setminus \Gamma_{\text{in}}, \end{array} \right. \quad (8) $$

where $\boldsymbol{u}_{\epsilon}$ is the solution to the Stokes equations Eq. (1), and the non-dimensional diffusive coefficient $\mathcal{K}_{\epsilon}$ and heat source are defined as follows:

$$ \mathcal{K}_{\epsilon}(\boldsymbol{x}) = \begin{cases} 1/\text{Pe}_f & \text{if } \boldsymbol{x} \in \Omega_{\epsilon,f}, \\ 1/\text{Pe}_s & \text{if } \boldsymbol{x} \in \Omega_{\epsilon,s}. \end{cases}, \quad Q(\boldsymbol{x}) = \begin{cases} 0 & \text{if } \boldsymbol{x} \in \Omega_{\epsilon,f}, \\ \dot{Q} & \text{if } \boldsymbol{x} \in \Omega_{\epsilon,s}. \end{cases} \quad (9) $$

The Peclét number is defined as follows:

$$ \mathrm{Pe} = \frac{\rho_{\epsilon} c_{p,\epsilon} UL}{\kappa_{\epsilon}}, \quad (10) $$

where $U$ is the characteristic velocity, and $L$ is the characteristic length. The Peclét number represents the ratio of convective to diffusive heat transport. A high Peclét number ($\mathrm{Pe} \gg 1$) indicates that convection dominates over diffusion, whereas a low Peclét number ($\mathrm{Pe} \ll 1$) indicates that diffusion dominates. Given that the fluid flow considered in this paper is a creeping flow, we assume that $\mathrm{Pe}_s = 1/50$, and that $\mathrm{Pe}_f = 1$.

The volumetric heat source is assumed to be design-dependent with a heat generation rate $\dot{Q} = 1.0$. Applying a volumetric heat source exclusively to the solid phase in a two-dimensional conjugate heat transfer problem is justified by the observation that, in numerous practical situations like electronic devices or solid structures subjected to external heating, heat generation predominantly occurs
where $\mathcal{K}^H$ represents the thermal conductivity tensor of the undeformed unit cell $Y$, defined as:

$$ 
\begin{aligned}
\mathcal{K}^H &= \left[ \frac{\mathcal{K}_{11}}{\mathcal{K}_{21}}; \frac{\mathcal{K}_{12}}{\mathcal{K}_{22}} \right] \\
&= \frac{1}{|Y|} \left[ \frac{\int_Y \kappa(y) (1 - \partial w_1 / \partial y_1) \, dy}{\int_Y \kappa(y) (-\partial w_1 / \partial y_2) \, dy} ; \frac{\int_Y \kappa(y) (-\partial w_2 / \partial y_1) \, dy}{\int_Y \kappa(y) (1 - \partial w_2 / \partial y_2) \, dy} \right],
\end{aligned} 
\quad (12) $$

where $\kappa(y)$ is the local thermal conductivity within the unit cell which varies with local space, and $w_1$ and $w_2$ are the solutions to the unit-cell problem which will be introduced in Section 2.2.2. Furthermore, due to the symmetric nature of the microstructure employed in this paper, the off-diagonal $\mathcal{K}_{12} = \mathcal{K}_{21} = 0$.

Given that $\mathbf{u} = -\epsilon^2 \mathbf{X}^*(m_1, m_2, \theta) \nabla p$, Eq. (10) can be rewritten as follows:

$$ 
\left\{ 
\begin{array}{ll}
-\epsilon^2 \mathbf{X}^*(m_1, m_2, \theta) \nabla p \cdot \nabla T - \operatorname{div}(\mathbf{K}^*(m_1, m_2, \theta) \nabla T) = Q(m_1, m_2) & \text{in } D, \\
T = T_0 & \text{on } \Gamma_{\text{in}}, \\
\mathbf{K}^*(m_1, m_2, \theta) \nabla T \cdot \mathbf{n} = 0 & \text{on } \partial D \setminus \Gamma_{\text{in}},
\end{array}
\right. 
\quad (13) $$

where the design-dependent volumetric heat source is given as $Q(m_1, m_2) = m_1 m_2 \dot{Q}$.

## 2.2 Unit cell problem

To determine the effective permeability matrix $\mathbf{X}^H(m_1, m_2)$ and thermal conductivity $\mathcal{K}^H(m_1, m_2)$, we parametrise the microstructure of the unit cell. In this study, the microstructure is represented as a square with side lengths $m_1$ and $m_2$, as depicted in the middle sub-figure of Fig. 2. Subsequently, the unit-cell problems are solved, and their formulations will be outlined below for clarity.

### 2.2.1 Fluid permeability

The homogenised permeability $\mathbf{X}_{ij}^H$ is determined by computing the volumetric average of the fluid velocity for the incompressible Stokes flow within the unit cell (Andreassen and Andreasen, 2014; Sanchez-Palencia, 1980). This can be represented in tensor form as follows:

$$ \mathbf{X}_{ij}^H = \frac{1}{|Y|} \int_{Y(m_1, m_2)} \mathbf{X}_j(\mathbf{y}) e_i \, dy, \quad 1 \le i, j \le 2, \quad (14) $$

where $e_i$ denotes the unit vector. $(\mathbf{X}_i(\mathbf{y}), \psi_i)$ represent the velocity and pressure, respectively to the unit-cell problem governed by the Stokes equation, where $\mathbf{X}_i(\mathbf{y})$ satisfies periodic boundary conditions. Additionally, a mean value of $\psi_i(\mathbf{y})$ is imposed to ensure solution uniqueness. The unit-cell problem can be formulated
within the solid material. Consequently, disregarding heat generation in the fluid phase is reasonable, as the primary function of the fluid is to facilitate heat transport rather than generate it. Then, the solution for temperature $T_\epsilon$ in the full-scale convection-diffusion equation can be approximated by the effective (or homogenised) temperature $T$ obtained by solving the following the homogenised problem (Allaire and Habibi, 2013):

$$ \left\{ \begin{array}{ll} \mathbf{u} \cdot \nabla T - \text{div}(\mathbf{K}^*(m_1, m_2, \theta) \nabla T) = Q(m_1, m_2) & \text{in } D, \\ T = T_0 & \text{on } \Gamma_{\text{in}}, \\ \mathbf{K}^*(m_1, m_2, \theta) \nabla T \cdot \mathbf{n} = 0 & \text{on } \partial D \setminus \Gamma_{\text{in}}, \end{array} \right. \quad (11) $$

where $\mathbf{K}^*(m_1, m_2, \theta)$ is the global thermal conductivity matrix, which can be expressed as follows:

$$ \mathbf{K}^*(m_1, m_2, \theta) = \mathbf{R}(\theta)\mathbf{K}^H(m_1, m_2)\mathbf{R}(\theta)^T, \quad (12) $$

where $\mathbf{K}^H$ represents the thermal conductivity tensor of the undeformed unit cell $Y$, defined as:

$$ \begin{aligned} \mathbf{K}^H &= \left[ \begin{array}{cc} \mathcal{K}_{11} & \mathcal{K}_{12} \\ \hline \mathcal{K}_{21} & \mathcal{K}_{22} \end{array} \right] \\ &= \frac{1}{|Y|} \left[ \begin{array}{cc} \int_Y \kappa(y) (1 - \partial w_1 / \partial y_1) \, dy & \int_Y \kappa(y) (-\partial w_2 / \partial y_1) \, dy \\ \hline \int_Y \kappa(y) (-\partial w_1 / \partial y_2) \, dy & \int_Y \kappa(y) (1 - \partial w_2 / \partial y_2) \, dy \end{array} \right], \end{aligned} \quad (13) $$

where $\kappa(y)$ is the local thermal conductivity within the unit cell which varies with local space, and $w_1$ and $w_2$ are the solutions to the unit-cell problem which will be introduced in Section 2.2.2. Furthermore, due to the symmetric nature of the microstructure employed in this paper, the off-diagonal $\mathcal{K}_{12} = \mathcal{K}_{21} = 0$.

Given that $\mathbf{u} = -\epsilon^2 \mathbf{X}^*(m_1, m_2, \theta) \nabla p$, Eq. (11) can be rewritten as follows:

$$ \left\{ \begin{array}{ll} -\epsilon^2 \mathbf{X}^*(m_1, m_2, \theta) \nabla p \cdot \nabla T - \text{div}(\mathbf{K}^*(m_1, m_2, \theta) \nabla T) = Q(m_1, m_2) & \text{in } D, \\ T = T_0 & \text{on } \Gamma_{\text{in}}, \\ \mathbf{K}^*(m_1, m_2, \theta) \nabla T \cdot \mathbf{n} = 0 & \text{on } \partial D \setminus \Gamma_{\text{in}}, \end{array} \right. \quad (14) $$

where the design-dependent volumetric heat source is given as $Q(m_1, m_2) = m_1 m_2 \dot{Q}$.

## 2.2 Unit cell problem

To determine the effective permeability matrix $\mathbf{X}^H(m_1, m_2)$ and thermal conductivity $\mathbf{K}^H(m_1, m_2)$, we parametrise the microstructure of the unit cell. In this study, the microstructure is represented as a square with side lengths $m_1$ and $m_2$, as depicted in the middle sub-figure of Fig. 2. Subsequently, the unit-cell problems are solved, and their formulations will be outlined below for clarity.

### 2.2.1 Fluid permeability

The homogenised permeability $\mathbf{X}_{ij}^H$ is determined by computing the volumetric average of the fluid velocity for the incompressible Stokes flow within the unit
as follows:

$$ egin{cases} 
\Delta \mathcal{X}_i(y) + \nabla \psi_i(y) = e_i - \alpha(y) \mathcal{X}_i(y) & \text{in } Y(m_1, m_2), \\
\text{div}(\mathcal{X}_i(y)) = 0 & \text{in } Y(m_1, m_2), \\
\int_{Y(m_1, m_2)} \psi_i(y) \, dy = 0, & \\
y \to \mathcal{X}_i(y) & \text{Y-periodic.}
\end{cases} 
\quad (15) $$

We present the solution to Eq. (15) on Fig. 3, considering the case where $m_1 = 0.4$ and $m_2 = 0.6$. It is worth noting that, for simplicity in implementation, we employ a structured mesh instead of a body-fitted one. Therefore, we adopt a fictitious body-force term within the Stokes equation to represent both the fluid and solid phases.

<figure><img src="image_4.png" /><figcaption>Fig. 3: Solution to the unit-cell problem defined in Eq. (15) for the case where $m_1 = 0.4$ and $m_2 = 0.6$.</figcaption></figure>

To construct the offline library, we solve the aforementioned unit-cell problem for all the possible parameter combinations ($m_1, m_2$) distributed across a grid discretizing the domain $(0, 1)^2$ into 100 × 100 points, eliminating 0 and 1 to avoid deterioration (a cell full of fluid and solid, respectively). Finite differences are then employed to estimate the partial derivatives $\partial \mathcal{X}_{11}/\partial m_1$, $\partial \mathcal{X}_{11}/\partial m_2$, $\partial \mathcal{X}_{22}/\partial m_1$,
cell (Sanchez-Palencia, 1980a; Andreassen and Andreasen, 2014). This can be represented in tensor form as follows:

$$ \boldsymbol{\mathcal{X}}_{ij}^{H} = \frac{1}{|Y|} \int_{Y(m_1, m_2)} \boldsymbol{\mathcal{X}}_{j}(\boldsymbol{y})e_i \, dy, \quad 1 \le i, j \le 2, \qquad (15) $$

where $e_i$ denotes the unit vector. ($\boldsymbol{\mathcal{X}}_i(\boldsymbol{y})$, $\psi_i$) represent the velocity and pressure, respectively to the unit-cell problem governed by the Stokes equation, where $\boldsymbol{\mathcal{X}}_i(\boldsymbol{y})$ satisfies periodic boundary conditions. Additionally, a mean value of $\psi_i(\boldsymbol{y})$ is imposed to ensure solution uniqueness. The unit-cell problem can be formulated as follows:

$$ \left\{ \begin{array}{ll} \Delta \boldsymbol{\mathcal{X}}_i(\boldsymbol{y}) + \nabla \psi_i(\boldsymbol{y}) = e_i - \alpha(\boldsymbol{y}) \boldsymbol{\mathcal{X}}_i(\boldsymbol{y}) & \text{in } Y(m_1, m_2), \\ \text{div}(\boldsymbol{\mathcal{X}}_i(\boldsymbol{y})) = 0 & \text{in } Y(m_1, m_2), \\ \int_{Y(m_1, m_2)} \psi_i(\boldsymbol{y}) \, dy = 0, & \\ y \to \boldsymbol{\mathcal{X}}_i(\boldsymbol{y}) & \text{Y-periodic.} \end{array} \right. \qquad (16) $$

We present the solution to Eq. (16) on Fig. 3, considering the case where $m_1 = 0.4$ and $m_2 = 0.6$. It is worth noting that, for simplicity in implementation, we employ a structured mesh instead of a body-fitted one. Therefore, we adopt a fictitious body-force term within the Stokes equation to represent both the fluid and solid phases.

To construct the offline library, we solve the aforementioned unit-cell problem for all the possible parameter combinations ($m_1, m_2$) distributed across a grid discretising the domain $(0, 1)^2$ into 100 × 100 points, eliminating 0 and 1 to avoid deterioration (a cell full of fluid and solid, respectively). Finite differences are then employed to estimate the partial derivatives $\partial\mathcal{X}_{11}/\partial m_1$, $\partial\mathcal{X}_{11}/\partial m_2$, $\partial\mathcal{X}_{22}/\partial m_1$, and $\partial\mathcal{X}_{22}/\partial m_2$. Subsequently, a linear interpolation is utilised to approximate the entire unit square. The linear interpolations of $\mathcal{X}_{11}$ and $\mathcal{X}_{22}$ are illustrated in Fig. 4.

### 2.2.2 Thermal conductivity

The homogenised thermal conductivity tensor $\mathcal{K}_{ij}^H$ in Eq. (13) can be written in tensor form as follows:

$$ \mathcal{K}_{ij}^{H} = \frac{1}{|Y|} \int_{Y(m_1, m_2)} \kappa(\boldsymbol{y}) \left( \delta_{ij} - \frac{\partial w_j(\boldsymbol{y})}{\partial y_i} \right) dy, \quad 1 \le i, j \le 2, \qquad (17) $$

where $w_i(\boldsymbol{y})$ can be obtained by solving the unit-cell problem governed by the Laplacian, with $w_i(\boldsymbol{y})$ subject to periodic boundary conditions. Additionally, a mean value of $w_i(\boldsymbol{y})$ is enforced to guarantee the uniqueness of the solution. The unit-cell problem is formulated as follows:

$$ \left\{ \begin{array}{ll} -\text{div}(\kappa(\boldsymbol{y})(e_i + \nabla w_i(\boldsymbol{y}))) = 0 & \text{in } Y(m_1, m_2), \\ \int_{Y(m_1, m_2)} w_i(\boldsymbol{y}) \, dy = 0, & \\ y \to w_i(\boldsymbol{y}) & \text{Y-periodic.} \end{array} \right. \qquad (18) $$
and $\partial \mathcal{X}_{22}/\partial m_2$. Subsequently, a linear interpolation is utilized to approximate the entire unit square. The linear interpolations of $\mathcal{X}_{11}$ and $\mathcal{X}_{22}$ are illustrated in Fig. 4.

<figure><img src="image_5.png" /><figcaption>Fig. 4: Linear interpolation of $\mathcal{X}_{11}$ and $\mathcal{X}_{22}$.</figcaption></figure>

### 2.2.2 Thermal conductivity

The homogenised thermal conductivity tensor $\mathcal{K}_{ij}^H$ in Eq. (12) can be written in tensor form as follows:

$$ \mathcal{K}_{ij}^H = \frac{1}{|Y|} \int_{Y(m_1, m_2)} \kappa(\mathbf{y}) \left( \delta_{ij} - \frac{\partial w_j(\mathbf{y})}{\partial y_i} \right) dy, \quad 1 \le i, j \le 2, \quad (16) $$

where $w_i(\mathbf{y})$ can be obtained by solving the unit-cell problem governed by the Laplacian, with $w_i(\mathbf{y})$ subject to periodic boundary conditions. Additionally, a mean value of $w_i(\mathbf{y})$ is enforced to guarantee the uniqueness of the solution. The unit-cell problem is formulated as follows:

$$ \left\{ \begin{array}{ll} -\operatorname{div}(\kappa(\mathbf{y})(e_i + \nabla w_i(\mathbf{y}))) = 0 & \text{in } Y(m_1, m_2), \\ \displaystyle\int_{Y(m_1, m_2)} w_i(\mathbf{y}) \, dy = 0, & \\ y \to w_i(\mathbf{y}) & \text{Y-periodic.} \end{array} \right. \quad (17) $$

We showcase the solution to Eq. (17) on Fig. 5, with $m_1 = 0.4$ and $m_2 = 0.6$. Following what has been discussed in Section 2.2.1, we construct an offline library for effective thermal conductivity, as can be seen in Fig. 6.

## 2.3 Filtering and projection

To mitigate the well-known checkerboard phenomenon, filtering techniques are employed (Lazarov and Sigmund, 2011; Kawamoto et al., 2011). The design variables $d = (m_1, m_2, \theta)$ are smoothened by solving the following PDE with $r$ the regularisation parameter:

$$ \left\{ \begin{array}{ll} -r^2 \Delta \tilde{d}_i + \tilde{d}_i = d_i & \text{in } \Omega, \\ \nabla \tilde{d}_i \cdot \mathbf{n} = 0 & \text{on } \partial \Omega. \end{array} \right. \quad (18) $$
<figure><img src="image_5.png" /><figcaption>(a) Characteristic function.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>(b) $\chi_1$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(c) $\psi_1$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(d) $\chi_2$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(e) $\psi_2$.</figcaption></figure>

Fig. 3: Solution to the unit-cell problem defined in Eq. (16) for the case where $m_1 = 0.4$ and $m_2 = 0.6$.

<figure><img src="image_10.png" /><figcaption>(a) $\chi_{11}(m_1, m_2)$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(b) $\chi_{22}(m_1, m_2)$.</figcaption></figure>

Fig. 4: Linear interpolation of $\chi_{11}$ and $\chi_{22}$.

We showcase the solution to Eq. (18) on Fig. 5, with $m_1 = 0.4$ and $m_2 = 0.6$. Following what has been discussed in Section 2.2.1, we construct an offline library for effective thermal conductivity, as can be seen in Fig. 6.
<figure><img src="image_6.png" /><figcaption>Fig. 5: Solution to the unit cell problem defined in Eq. (17) for the case where $m_1 = 0.4$, $m_2 = 0.6$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Fig. 6: Linear interpolation of $K_{11}$ and $K_{22}$.</figcaption></figure>

To improve the subsequent de-homogenised designs, it is crucial to avoid microstructures characterised by values of $m_1$ and $m_2$ that are close to, but not exactly equal to 0 or 1. Such values represent structures with extremely narrow solid or fluid strips, which are not ideal for manufacturing when projected on a fine-scale mesh. To address this issue, a projection technique, as described in Groen and Sigmund (2018), is employed. This projects the intermediate values of $m_1$ and $m_2$ to 0 or 1 outside the range $(\eta, 1-\eta)$, with the threshold $\eta$ and projection sharpness $\beta$. Thus, the filtered design variables $\tilde{m} = (\tilde{m}_1, \tilde{m}_2)$ are projected to yield $\hat{\tilde{m}} = (\hat{\tilde{m}}_1, \hat{\tilde{m}}_2)$ as follows:

$$ \hat{\tilde{m}}_i(\tilde{m}_i, \beta, \eta) := \tilde{m}_i(1 - H(\beta, 1-\eta, \tilde{m}_i)) H(\beta, \eta, \tilde{m}_i) + \left(\frac{\beta-1}{\beta} + \frac{\tilde{m}_i}{\beta}\right) H(\beta, 1-\eta, \tilde{m}_i), \quad (19) $$

where $H$ represents a smoothed Heaviside function that is widely used in the density-based topology optimisation, defined as follows:

$$ H(\tilde{m}_i, \beta, \eta) := \frac{\tanh(\beta\eta) + \tanh(\beta(\tilde{m}_i - \eta))}{\tanh(\beta\eta) + \tanh(\beta(1-\eta))}, \quad (20) $$

To enhance clarity for the reader, we illustrate the projection function with various settings of $\eta$ and $\beta$ in Fig. 7. In this paper, we set $\eta = 0.05$ and $\beta = 100$.
<figure><img src="image_12.png" /><figcaption>Fig. 5: Solution to the unit cell problem defined in Eq. (18) for the case where $m_1 = 0.4$, $m_2 = 0.6$.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 6: Linear interpolation of $K_{11}$ and $K_{22}$.</figcaption></figure>

## 2.3 Filtering and projection

To mitigate the well-known checkerboard phenomenon, filtering techniques are employed (Lazarov and Sigmund, 2011; Kawamoto et al., 2011). The design variables $d = (m_1, m_2, \theta)$ are smoothened by solving the following PDE with $r$ the regularisation parameter:

$$ \begin{cases} -r^2 \Delta \tilde{d}_i + \tilde{d}_i = d_i & \text{in } \Omega, \\ \nabla \tilde{d}_i \cdot \mathbf{n} = 0 & \text{on } \partial\Omega. \end{cases} \qquad (19) $$

To improve the subsequent de-homogenised designs, it is crucial to avoid microstructures characterised by values of $m_1$ and $m_2$ that are close to, but not exactly equal to 0 or 1. Such values represent structures with extremely narrow solid or fluid strips, which are not ideal for manufacturing when projected on a fine-scale mesh. To address this issue, a projection technique, as described in Groen and Sigmund (2018), is employed. This projects the intermediate values of $m_1$ and $m_2$ to 0 or 1 outside the range $(\eta, 1-\eta)$, with the threshold $\eta$ and projection sharpness $\beta$. Thus, the filtered design variables $\tilde{m} = (\tilde{m}_1, \tilde{m}_2)$ are projected to yield
<figure><img src="image_8.png" /><figcaption>Fig. 7: Projection function $\hat{m}_i(\tilde{m}_i, \beta, \eta)$, cf. Eq. (19).</figcaption></figure>

In addition, the filtered orientation angle $\tilde{\theta} \in [0, 1]$ is linearly rescaled to make it
in radians $\hat{\theta} \in [-\pi, \pi]$.

## 2.4 Optimisation problem

<figure><img src="image_9.png" /><figcaption>Fig. 8: Schematic diagram of a symmetric microchannel cooling design case.</figcaption></figure>

In the context of a two-dimensional conjugate heat transfer problem with a design-
dependent heat source applied within the solid phase, the considered objective is to
minimize the temperature in the solid regions where heat generation occurs. Lower
temperatures signify reduced heat accumulation and better cooling. Concurrently,
we aim to adhere to a maximum volume fraction constraint for the fluid phase and a
maximum allowable pressure drop. The schematic diagram of the design settings is
illustrated in Fig. 8. The optimization mathematical model is formulated as follows:
$\hat{\tilde{m}} = (\hat{\tilde{m}}_1, \hat{\tilde{m}}_2)$ as follows:

$$ \hat{\tilde{m}}_i(\tilde{m}_i, \beta, \eta) := \tilde{m}_i(1 - H(\beta, 1-\eta, \tilde{m}_i)) H(\beta, \eta, \tilde{m}_i) + \left(\frac{\beta-1}{\beta} + \frac{\tilde{m}_i}{\beta}\right) H(\beta, 1-\eta, \tilde{m}_i), \quad (20) $$

where $H$ represents a smoothed Heaviside function that is widely used in the density-based topology optimisation, defined as follows:

$$ H(\tilde{m}_i, \beta, \eta) := \frac{\tanh(\beta\eta) + \tanh(\beta(\tilde{m}_i - \eta))}{\tanh(\beta\eta) + \tanh(\beta(1-\eta))}, \quad (21) $$

To enhance clarity for the reader, we illustrate the projection function with various settings of $\eta$ and $\beta$ in Fig. 7. In this paper, we set $\eta = 0.05$ and $\beta = 100$.

<figure><img src="image_14.png" /><figcaption>Fig. 7: Projection function $\hat{\tilde{m}}_i(\tilde{m}_i, \beta, \eta)$, cf. Eq. (20).</figcaption></figure>

In addition, the filtered orientation angle $\tilde{\theta} \in [0, 1]$ is linearly rescaled to make it in radians $\tilde{\theta} \in [-\pi, \pi]$.

## 2.4 Optimisation problem

In the context of a two-dimensional conjugate heat transfer problem with a design-dependent heat source applied within the solid phase, the considered objective is to minimise the temperature in the solid regions where heat generation occurs. Lower temperatures signify reduced heat accumulation and better cooling. Concurrently, we aim to adhere to a maximum volume fraction constraint for the fluid phase and a maximum allowable pressure drop. Note that without this maximum fluid fraction constraint, the optimisation solver would favour a design dominated entirely by fluid, as the heat source is only applied in the solid, so by removing all solid, it reduces the temperatures to a minimum ("infinite heat transfer"). The schematic diagram of the design settings is illustrated in Fig. 8. The optimisation mathematical model is formulated as follows:
$$
\begin{equation}
\begin{aligned}
\min_{m_1, m_2, \theta \in D} & J(\Omega) = \int_D \hat{m}_1 \hat{m}_2 T \, d\Omega, \tag{21a} \\
\text{s.t.} \quad & \left\{
\begin{array}{@{}l@{}}
\text{Eqs. (7) and (13),} \\
G_1 = \epsilon^2 \frac{1/\text{DP}_0}{|\Gamma_{\text{in}}|} \int_{\Gamma_{\text{in}}} p \, d\Gamma - C_{\text{DP}} \le 0, \\
G_2 = \frac{\int_D 1 - \hat{m}_1 \hat{m}_2 \, d\Omega}{\int_D \, d\Omega} - V_{\text{max}} \le 0, \\
0 \le m_1(\mathbf{x}), m_2(\mathbf{x}), \theta(\mathbf{x}) \le 1 \quad \forall \mathbf{x} \in D,
\end{array}
\right.
\end{aligned}
\tag{21b}
\end{equation}
$$

where $G_1$ and $G_2$ denote the pressure drop constraint and volume constraint, respectively, $DP_0$ signifies the pressure drop at the initial iteration, $C_{DP}$ indicates the maximum permitted pressure drop fraction, and $V_{max}$ represents the maximum allowable volume fraction for the fluid phase. The sensitivity analysis is conducted using the continuous adjoint method and it is detailed in Appendix A.

## 2.5 Reconstruction of microchannels

To reconstruct the fine-scale microchannels within $\Omega_{\epsilon}$ using the optimised projected design fields $\hat{\mathbf{d}} = (\hat{m}_1, \hat{m}_2, \hat{\theta})$, which are continuously varying within the design domain $D$, we follow the lead of existing de-homogenisation techniques (Groen and Sigmund, 2018; Feppon, 2024b), based on the observation that the orientation field is free of singularities, as evidenced in Fig. 10a, where we solve a pressure drop minimisation benchmark problem as illustrated in the schematic diagram in Fig. 9. However, it is important to note that there is no guarantee of singularity-free orientation field. In case singularities are observed, one can consider employing the regularisation algorithm proposed in Allaire et al. (2019b) to achieve more regular mapped designs. In this subsection, we briefly outline the formulation employed for this procedure to make this paper self-contained.

The unit cells are arranged in a locally-periodic manner on a regular grid in Cartesian coordinates and undergo a nonlinear transformation $\Phi$ to become a curved grid in distorted coordinates. Assuming this transformation $\Phi$ is conformal, meaning that $\nabla\Phi(\phi_1, \phi_2) \simeq \mathbf{R}(\hat{\theta})$, the problem yields to find appropriate functions $\phi_1$ and $\phi_2$ such that $\nabla\phi_1 \simeq e_1(\hat{\theta})$ and $\nabla\phi_2 \simeq e_2(\hat{\theta})$, where $e_1(\hat{\theta})$ and $e_2(\hat{\theta})$ represent the local orientation vectors as illustrated in Fig. 10b, which are defined as follows:

$$ e_1(\hat{\theta}) := \begin{bmatrix} \cos(\hat{\theta}) \\ \sin(\hat{\theta}) \end{bmatrix}, \quad e_2(\hat{\theta}) := \begin{bmatrix} -\sin(\hat{\theta}) \\ \cos(\hat{\theta}) \end{bmatrix}. \qquad (22) $$

Taking $\phi_1$ as an example, the constrained least-squares minimisation problem is then formulated as follows (Groen and Sigmund, 2018):

$$
\begin{align}
\min_{\phi_1 \in D} \quad & \mathcal{I}(\phi_1(\mathbf{x})) = \frac{1}{2} \int_{\Omega} a_1(\mathbf{x}) \left\| \nabla \phi_1(\mathbf{x}) - e_1(\hat{\theta}) \right\|^2 d\Omega, \tag{23} \\
\text{s.t.} \quad & a_2(\mathbf{x}) \nabla \phi_1(\mathbf{x}) \cdot e_2(\hat{\theta}) = 0, \nonumber
\end{align}
$$
<figure><img src="image_15.png" /><figcaption>Fig. 8: Schematic diagram of a symmetric microchannel cooling design case.</figcaption></figure>

$$ \min_{m_1, m_2, \theta \in D} J(\Omega) = \int_D \hat{m}_1 \hat{m}_2 T \, d\Omega, \qquad (22a) $$

$$ \begin{array}{l} \text{s.t.} \quad \left\{ \begin{array}{l} G_1 = \epsilon^2 \frac{1/\text{DP}_0}{|\Gamma_{\text{in}}|} \int_{\Gamma_{\text{in}}} p \, d\Gamma - C_{\text{DP}} \le 0, \\ G_2 = \frac{\int_D 1 - \hat{m}_1 \hat{m}_2 \, d\Omega}{\int_D \, d\Omega} - V_{\text{max}} \le 0, \\ 0 \le m_1(\mathbf{x}), m_2(\mathbf{x}), \theta(\mathbf{x}) \le 1 \quad \forall \mathbf{x} \in D, \end{array} \right. \\ \text{Eqs. (7) and (14),} \end{array} \qquad (22b) $$

where $G_1$ and $G_2$ denote the pressure drop constraint and volume constraint, respectively, DP₀ signifies the pressure drop at the initial iteration, $C_{\text{DP}}$ indicates the maximum permitted pressure drop fraction, and $V_{\text{max}}$ represents the maximum allowable volume fraction for the fluid phase. The sensitivity analysis is conducted using the continuous adjoint or discrete adjoint method and they are detailed in Appendix A.

## 2.5 Reconstruction of microchannels

To reconstruct the fine-scale microchannels within $\Omega_{\epsilon}$ using the optimised projected design fields $\hat{d} = (\hat{m}_1, \hat{m}_2, \hat{\theta})$, which are continuously varying within the design domain $D$, we follow the lead of existing de-homogenisation techniques (Groen and Sigmund, 2018; Feppon, 2024b), based on the observation that the orientation field is free of singularities, as evidenced in Fig. 10a, where we solve a pressure drop minimisation benchmark problem as illustrated in the schematic diagram in Fig. 9. However, it is important to note that there is no guarantee of singularity-free orientation field. In case singularities are observed, one can consider employing the regularisation algorithm proposed in Allaire et al. (2019b) to achieve more regular mapped designs. In this subsection, we briefly outline the formulation employed for this procedure to make this paper self-contained.
<figure><img src="image_10.png" /><figcaption>Fig. 9: Schematic diagram illustrating a benchmark test case for minimising pressure drop, intended for demonstrating de-homogenisation.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 10: Optimised orientation field obtained through homogenisation-based topology optimisation for minimising pressure drop.</figcaption></figure>

The introduction of weighting factor $a_1$ aims to ease the projection of $\phi_1$ within $\Omega_{\epsilon,f}$ and $\Omega_{\epsilon,s}$, while maintaining regularization for lattice spacing. Additionally, $a_2$ is employed to disable strict angular enforcement in these regions (Groen and Sigmund, 2018). The weighting factors $a_1(x)$ and $a_2(x)$ are defined as follows:

$$ a_1(\boldsymbol{x}) := \begin{cases} 10^{-3} & \text{if } \boldsymbol{x} \in \Omega_f, \\ 10^{-3} & \text{if } \boldsymbol{x} \in \Omega_s, \\ 1 & \text{if } \boldsymbol{x} \in \Omega \setminus (\Omega_s \cup \Omega_f), \end{cases} \qquad a_2(\boldsymbol{x}) := \begin{cases} 0 & \text{if } \boldsymbol{x} \in \Omega_f, \\ 0 & \text{if } \boldsymbol{x} \in \Omega_s, \\ 10 & \text{if } \boldsymbol{x} \in \Omega \setminus (\Omega_s \cup \Omega_f). \end{cases} \qquad (24) $$

The domain $\Omega$ comprises three parts: a fluid domain $\Omega_f$, a solid domain $\Omega_{\epsilon,s}$, and a porous medium domain $\Omega \setminus (\Omega_s \cup \Omega_f)$, which can be determined by the correlation between the local volume fraction and the threshold value $\eta$ used in the projection scheme given in Eq. (19), as follows:

$$ \boldsymbol{x} \in \begin{cases} \Omega_f & \text{if } 1 - \hat{m}_1 \hat{m}_2 \le \eta, \\ \Omega_s & \text{if } 1 - \hat{m}_1 \hat{m}_2 \ge 1 - \eta, \\ \Omega \setminus (\Omega_s \cup \Omega_f) & \text{if } \eta < 1 - \hat{m}_1 \hat{m}_2 < 1 - \eta. \end{cases} \qquad (25) $$
The unit cells are arranged in a locally-periodic manner on a regular grid in Cartesian coordinates and undergo a nonlinear transformation $\Phi$ to become a curved grid in distorted coordinates. Assuming this transformation $\Phi$ is conformal, meaning that $\nabla\Phi(\phi_1, \phi_2) \simeq R(\hat{\theta})$, the problem yields to find appropriate functions $\phi_1$ and $\phi_2$ such that $\nabla\phi_1 \simeq e_1(\hat{\theta})$ and $\nabla\phi_2 \simeq e_2(\hat{\theta})$, where $e_1(\hat{\theta})$ and $e_2(\hat{\theta})$ represent the local orientation vectors as illustrated in Fig. 10b, which are defined as follows:

$$e_1(\hat{\theta}) := \begin{bmatrix} \cos(\hat{\theta}) \\ \sin(\hat{\theta}) \end{bmatrix}, \quad e_2(\hat{\theta}) := \begin{bmatrix} -\sin(\hat{\theta}) \\ \cos(\hat{\theta}) \end{bmatrix}. \qquad (23)$$

<figure><img src="image_16.png" /><figcaption>Fig. 9: Schematic diagram illustrating a benchmark test case for minimising pressure drop, intended for demonstrating de-homogenisation.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 10: Optimised orientation field obtained through homogenisation-based topology optimisation for minimising pressure drop.</figcaption></figure>
Then, the constrained least-squares minimisation problem defined in Eq. (23) is rewritten as an unconstrained minimisation problem as follows:

$$ \min_{\phi_1 \in \mathcal{D}} \mathcal{I}(\phi_1(\mathbf{x})) = \int_{\Omega} a_1(\mathbf{x}) \left\| \nabla \phi_1(\mathbf{x}) - e_1(\hat{\theta}) \right\|^2 d\Omega + \int_{\Omega} a_2(\mathbf{x}) \left\| \nabla \phi_1 \cdot e_2(\hat{\theta}) \right\|^2 d\Omega. \quad (26) $$

$\phi_2$ can be determined by solving the above-mentioned minimisation problem wherein $e_1(\hat{\theta})$ and $e_2(\hat{\theta})$ should be interchanged. The resulting level-set functions, corresponding to the orientation field depicted in Fig. 10a, are illustrated in Fig. 11.

<figure><img src="image_12.png" /><figcaption>Fig. 11: The level-set functions as the solution of the least-squares minimisation problem given in Eq. (26).</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Fig. 12: Optimised fields of $\hat{m}_1$ and $\hat{m}_2$ obtained through homogenisation-based topology optimisation for minimising pressure drop.</figcaption></figure>

Subsequently, given distributions of $\hat{m}_1$ and $\hat{m}_2$ illustrated in Fig. 12 alongside the transformation $\Phi$ delineated by the level-set functions $\phi_1$ and $\phi_2$, we can characterise the fluid domain on the curved grid with periodicity $\epsilon$ utilising the negative
Taking $\phi_1$ as an example, the constrained least-squares minimisation problem is then formulated as follows (Groen and Sigmund, 2018):

$$ \begin{array}{ll} \min_{\phi_1 \in \mathcal{D}} & \mathcal{I}(\phi_1(\mathbf{x})) = \frac{1}{2} \int_{\Omega} a_1(\mathbf{x}) \left\| \nabla \phi_1(\mathbf{x}) - e_1(\hat{\theta}) \right\|^2 d\Omega, \\ \text{s.t.} & a_2(\mathbf{x}) \nabla \phi_1(\mathbf{x}) \cdot e_2(\hat{\theta}) = 0, \end{array} \quad (24) $$

The introduction of weighting factor $a_1$ aims to ease the projection of $\phi_1$ within $\Omega_{\epsilon,f}$ and $\Omega_{\epsilon,s}$, while maintaining regularization for lattice spacing. Additionally, $a_2$ is employed to disable strict angular enforcement in these regions (Groen and Sigmund, 2018). The weighting factors $a_1(x)$ and $a_2(x)$ are defined as follows:

$$ a_1(\mathbf{x}) := \begin{cases} 10^{-3} & \text{if } \mathbf{x} \in \Omega_f, \\ 10^{-3} & \text{if } \mathbf{x} \in \Omega_s, \\ 1 & \text{if } \mathbf{x} \in \Omega \setminus (\Omega_s \cup \Omega_f), \end{cases} \qquad a_2(\mathbf{x}) := \begin{cases} 0 & \text{if } \mathbf{x} \in \Omega_f, \\ 0 & \text{if } \mathbf{x} \in \Omega_s, \\ 10 & \text{if } \mathbf{x} \in \Omega \setminus (\Omega_s \cup \Omega_f). \end{cases} \quad (25) $$

The domain $\Omega$ comprises three parts: a fluid domain $\Omega_f$, a solid domain $\Omega_{\epsilon,s}$, and a porous medium domain $\Omega \setminus (\Omega_s \cup \Omega_f)$, which can be determined by the correlation between the local volume fraction and the threshold value $\eta$ used in the projection scheme given in Eq. (20), as follows:

$$ \mathbf{x} \in \begin{cases} \Omega_f & \text{if } 1 - \hat{m}_1 \hat{m}_2 \le \eta, \\ \Omega_s & \text{if } 1 - \hat{m}_1 \hat{m}_2 \ge 1 - \eta, \\ \Omega \setminus (\Omega_s \cup \Omega_f) & \text{if } \eta < 1 - \hat{m}_1 \hat{m}_2 < 1 - \eta. \end{cases} \quad (26) $$

Then, the constrained least-squares minimisation problem defined in Eq. (24) is rewritten as an unconstrained minimisation problem as follows:

$$ \min_{\phi_1 \in \mathcal{D}} \mathcal{I}(\phi_1(\mathbf{x})) = \int_{\Omega} a_1(\mathbf{x}) \left\| \nabla \phi_1(\mathbf{x}) - e_1(\hat{\theta}) \right\|^2 d\Omega + \int_{\Omega} a_2(\mathbf{x}) \left\| \nabla \phi_1 \cdot e_2(\hat{\theta}) \right\|^2 d\Omega. \quad (27) $$

$\phi_2$ can be determined by solving the above-mentioned minimisation problem wherein $e_1(\hat{\theta})$ and $e_2(\hat{\theta})$ should be interchanged. The resulting level-set functions, corresponding to the orientation field depicted in Fig. 10a, are illustrated in Fig. 11.

Subsequently, given distributions of $\hat{m}_1$ and $\hat{m}_2$ illustrated in Fig. 12 alongside the transformation $\Phi$ delineated by the level-set functions $\phi_1$ and $\phi_2$, we can characterise the fluid domain on the curved grid with periodicity $\epsilon$ utilising the negative subdomain of the following pseudo-density functions (Feppon, 2024b):

$$ \gamma_{\epsilon}(\mathbf{x}) = -\max \left( \left| \text{Saw} \left( \frac{\phi_1(\mathbf{x})}{\epsilon}, T \right) \right| - \frac{\hat{m}_1(\mathbf{x})}{2}, \left| \text{Saw} \left( \frac{\phi_2(\mathbf{x})}{\epsilon}, T \right) \right| - \frac{\hat{m}_2(\mathbf{x})}{2} \right), \quad (28) $$

in which the sawtooth-function-based mapping enables the transformation of unit-cell microstructures into a continuous tiling pattern that can adapt to spatially varying properties. This is particularly important in multi-scale topology optimisation, where microstructures need to be oriented and morphed to fit the global design without compromising their physical properties (Wang et al., 2022; Groen
subdomain of the following pseudo-density functions (Feppon, 2024b):

$$ \gamma_{\epsilon}(\mathbf{x}) = -\max\left(\left|\text{Saw}\left(\frac{\phi_1(\mathbf{x})}{\epsilon}, T\right)\right| - \frac{\hat{m}_1(\mathbf{x})}{2}, \left|\text{Saw}\left(\frac{\phi_2(\mathbf{x})}{\epsilon}, T\right)\right| - \frac{\hat{m}_2(\mathbf{x})}{2}\right), \quad (27) $$

in which the sawtooth-function-based mapping ensures effective connections between oriented microstructures and maintain the original homogenised properties thanks to its conformality (Wang et al., 2022; Groen et al., 2019; Feppon, 2024b). It allows extending a function defining an implicit domain in the range $(-1/2, 1/2)^2$ through periodicity, defined as follows:

$$ \text{Saw}(t, T) := \frac{t}{T} - \left\lfloor \frac{t}{T} + 0.5 \right\rfloor, \quad (28) $$

with the periodicity $T = 1$. The de-homogenised (or upsampled) structure is illustrated in Fig. 13a, where the fluid subdomain $\Omega_{\epsilon,f}$ is defined as the negative subdomain of $\gamma_{\epsilon}$ as shown in Fig. 13b. It can be observed that “pin fins” for the minimal pressure drop problem only occur due to the filter and that the underlying problem prefers 1/0 binary structure.

<figure><img src="image_14.png" /><figcaption>Fig. 13: Interpretation of the homogenised design fields for a minimal pressure drop problem, where black regions represent solid material and white regions represent fluid.</figcaption></figure>

## 3 Implementation details

In this paper, the finite element method is used to solve the PDEs presented in Section 2. First-order elements are employed for all state variables. We use FreeFEM (Hecht, 2012; Jolivet et al., 2012), an open-source finite element software, for discretizing the PDEs, with PETSc (Balay et al., 2023) serving as the distributed linear algebra backend. It is worth noting that the homogenisation-based topology optimisation is conducted using a single-core code, while the validation of upsampled structures is executed in a distributed fashion. All the numerical experiments are conducted on an Apple MacBook Pro with a 10-core M2 processor. In the subsequent part of this section, a detailed explanation of the weak form for the PDEs utilized in this workflow will be provided.
<figure><img src="image_18.png" /><figcaption>**Fig. 11:** The level-set functions as the solution of the least-squares minimisation problem given in Eq. (27).</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Fig. 12:** Optimised fields of $\hat{m}_1$ and $\hat{m}_2$ obtained through homogenisation-based topology optimisation for minimising pressure drop.</figcaption></figure>

et al., 2019; Feppon, 2024b). The role of the sawtooth function can be summarised as follows:

1. Periodic field representation: the sawtooth function is periodic and can be used to define the spatial repetition of the microstructures. In 2D, it propagates in a specific direction determined by a unit vector. This makes it ideal for orienting microstructures while maintaining geometric compatibility between adjacent cells. By adjusting the periodicity and direction of the sawtooth function, it is possible to smoothly transition between different microstructures, ensuring that they fit together without discontinuities.

2. Conformal mapping: the sawtooth function creates a mapping that is conformal, meaning it preserves the angles between the microstructures, which is crucial for maintaining the intended material properties. This conformality ensures that when the microstructures are dehomogenised, the overall structural properties are consistent with the homogenised design.

3. Spatial adaptation: the sawtooth function is flexible in that it allows for the design of unit-cell orientations that vary across the design space. This is particularly useful in more complex design problems, where the optimal microstructures may change from one cell to another. The sawtooth function ensures that these
## 3.1 Offline library

As has been discussed in Section 2.2, the first step of this workflow is to construct the offline library by solving the unit-cell problems. Here, we define the Hilbert spaces $\mathcal{U}_{\text{per}}$ for $\mathcal{X}_i$ and $v_i$, $\mathcal{R}$ for $\psi_i$ and $q_i$, and $\mathcal{H}_{\text{per}}$ for $w_i$ and $v_i$, where $v_i$, $q_i$, and $v_i$ are the test functions, as follows:

$$ \mathcal{U}_{\text{per}} := \{v \in H_{\text{per}}^{1}(Y, \mathbb{R}^{d}); Y = [0, 1]^{2}\}, $$

$$ \mathcal{R} := \{p \in L^2(Y); \mathcal{M}_Y(p) = 0; Y = [0, 1]^2\}, $$

$$ \mathcal{H}_{\text{per}} := \{u \in H_{\text{per}}^{1}(Y); \mathcal{M}_{Y}(u) = 0; Y = [0, 1]^2\}. $$

The unit-cell problem for the Stokes flow defined in Eq. (15) can be rewritten in the weak form as:

Find $(\mathcal{X}_i, \psi_i, \lambda) \in \mathcal{U}_{\text{per}} \times \mathcal{R} \times \mathbb{R}$, such that:

$$ \begin{aligned} a(\mathcal{X}_i, v_i) + b(v_i, \psi_i) &= l(v_i) && \forall v_i \in \mathcal{U}_{\text{per}}, \\ b(\mathcal{X}_i, q_i) - \tau d(\psi_i, q_i) + \lambda m(q_i) &= 0 && \forall q_i \in \mathcal{R}, \\ m(\psi_i) &= 0, \end{aligned} $$

where we define the bilinear form $a(\mathcal{X}_i, v_i)$ and $b(v_i, \psi_i)$, linear form $l(v_i)$, pressure stabilisation term $\tau d(\psi_i, q_i)$, with the element size $h_K$ and Lagrange multiplier $\lambda$, as follows:

$$ a(\mathcal{X}_i, v_i) = \int_Y \nabla \mathcal{X}_i : \nabla v_i \, dy + \int_Y \alpha(y) \mathcal{X}_i \cdot v_i \, dy, $$

$$ b(v_i, \psi_i) = - \int_Y \nabla \cdot v_i \psi_i \, dy $$

$$ l(v_i) = \int_Y e_i \cdot v_i \, dy, $$

$$ \tau > 0, \quad d(\psi_i, q_i) = \sum_{K \in \mathcal{T}} h_K^2 \int_{Y_e} \nabla \psi_i \cdot \nabla q_i \, dy_e, $$

$$ m(\psi_i) = \int_Y \psi_i \, dy. $$

In our implementation, the discretisation of the above weak form is stored in the block matrix as follows:

$$ \mathbf{A}_{\text{Stokes}} = \begin{bmatrix} \mathbf{A} & \mathbf{B} & \mathbf{0} \\ \mathbf{B}^{\mathrm{T}} & \mathbf{S} & \mathbf{h} \\ \mathbf{0} & \mathbf{h}^{\mathrm{T}} & \mathbf{0} \end{bmatrix}, \qquad (29) $$

where $\mathbf{A}$ is the matrix discretising the bilinear form of $a(\mathcal{X}_i, v_i)$, $\mathbf{B}$ the divergence operator $b(v_i, \psi_i)$, $\mathbf{S}$ the pressure stabilisation term $-\tau d(\psi_i, q_i)$. $\mathbf{h}$ is the discretisation of the linear form associated with the Lagrange multiplier term $\lambda m(q_i)$, imposing a zero-mean value for $\psi_i$ such that the solution can be uniquely determined.
transitions are smooth, avoiding abrupt changes that could lead to undesirable effects in the performance of the structure.

More specifically, It allows extending a function defining an implicit domain in the range $(-1/2, 1/2)^2$ through periodicity, defined as follows:

$$ \mathrm{Saw}(t, T) := \frac{t}{T} - \left\lfloor \frac{t}{T} + 0.5 \right\rfloor, \qquad (29) $$

with the periodicity $T = 1$. The de-homogenised (or upsampled) structure is illustrated in Fig. 13a, where the fluid subdomain $\Omega_{\epsilon,f}$ is defined as the negative subdomain of $\gamma_{\epsilon}$ as shown in Fig. 13b. It can be observed that “pin fins” for the minimal pressure drop problem only occur due to the filter and that the underlying problem prefers 1/0 binary structure.

<figure><img src="image_20.png" /><figcaption>Fig. 13: Interpretation of the homogenised design fields for a minimal pressure drop problem, where black regions represent solid material and white regions represent fluid.</figcaption></figure>

## 3 Implementation details

In this paper, the finite element method is used to solve the PDEs presented in Section 2. First-order elements are employed for all state variables. We use FreeFEM (Hecht, 2012; Jolivet et al., 2012), an open-source finite element software, for discretising the PDEs, with PETSc (Balay et al., 2023) serving as the distributed linear algebra backend. It is worth noting that the homogenisation-based topology optimisation is conducted using a single-core code, while the validation of upsampled structures is executed in a distributed fashion. All the numerical experiments are conducted on an Apple MacBook Pro with a 10-core M2 processor. In the subsequent part of this section, a detailed explanation of the weak form for the PDEs utilised in this workflow will be provided.

### 3.1 Offline library

As has been discussed in Section 2.2, the first step of this workflow is to construct the offline library by solving the unit-cell problems. Here, we define the Hilbert
The unit cell problem for the heat transfer defined in Eq. (17) can be rewritten in the weak form as:

Find $(w_i, \lambda) \in \mathcal{H}_{\text{per}} \times \mathbb{R}$, such that:

$$ 
\begin{array}{l}
a(w_i, v_i) + \lambda m(v_i) = l(v_i) \quad \forall v_i \in \mathcal{H}_{\text{per}}, \\
m(w_i) = 0,
\end{array}
$$

where we define the bilinear form $a(w_i, v_i)$, linear form $l(v_i)$, with the Lagrange multiplier $\lambda$, as follows:

$$ 
\begin{array}{l}
a(w_i, v_i) = \displaystyle\int_Y \kappa(\mathbf{y}) \nabla w_i \cdot \nabla v_i \, dy, \\
l(v_i) = \displaystyle\int_Y \kappa(\mathbf{y}) e_i \cdot \nabla v_i \, dy, \\
m(w_i) = \displaystyle\int_Y w_i \, dy.
\end{array}
$$

The discretisation of the above weak form is stored in the block matrix as follows:

$$ \mathbf{A}_{\text{Laplacian}} = \begin{bmatrix} \mathbf{A} & \mathbf{h} \\ \mathbf{h}^T & \mathbf{0} \end{bmatrix}, \qquad (30) $$

where $\mathbf{A}$ is the matrix discretising the bilinear form $a(w_i, v_i)$. $\mathbf{h}$ is the discretisation of the linear form associated with the Lagrange multiplier term $\lambda m(v_i)$, imposing a zero-mean value for $w_i$, such that the solution can be uniquely determined.

## 3.2 Homogenisation-based topology optimisation

Now, we proceed with the homogenisation-based topology optimisation to solve the optimum design problem formulated in Section 2.4. Here, we define the Hilbert spaces $\mathcal{P}$ for the homogenised pressure $p$ and homogenised temperature $T$, their adjoint variables $p_a$ and $T_a$, as well as the test functions $q$ and $S$, as follows:

$$ \mathcal{P} := \{ \tilde{p} \in H^1(\Omega) \mid \tilde{p} = 0 \text{ on } \partial \Omega_p^D \}. $$

First, the Laplacian defined in Eq. (7) can be expressed in the following weak form: Find $p \in \mathcal{P}$ with $p = 0$ on $\Gamma_{\text{out}}$ such that:

$$ \int_D \mathbf{X}^* (\hat{\mathbf{m}}_1, \hat{\mathbf{m}}_2, \hat{\theta}) \nabla p \cdot \nabla q \, d\Omega = -\epsilon^{-2} \int_{\Gamma_{\text{in}}} (\mathbf{u}_0 \cdot \mathbf{n}) q \, d\Gamma \quad \forall q \in \mathcal{P}, p \in \mathcal{P}. \quad (31) $$

Next, the convection-diffusion equation defined in Eq. (13) can be written in weak form as:

Find $T \in \mathcal{P}$ with $T = 0$ on $\Gamma_{\text{in}}$ such that:

$$ 
\begin{align} 
& \int_D -\epsilon^2 \mathbf{X}^* (\hat{\mathbf{m}}_1, \hat{\mathbf{m}}_2, \hat{\theta}) \nabla p \cdot \nabla T S \, d\Omega + \int_D \mathbf{K}^* (\hat{\mathbf{m}}_1, \hat{\mathbf{m}}_2, \hat{\theta}) \nabla T \cdot \nabla S \, d\Omega \tag{32} \\
& = \int_D Q (\hat{\mathbf{m}}_1, \hat{\mathbf{m}}_2) S \, d\Omega \quad \forall S \in \mathcal{P}, T \in \mathcal{P}. 
\end{align} 
$$
spaces $\mathcal{U}_{\text{per}}$ for $\boldsymbol{X}_i$ and $\boldsymbol{v}_i$, $\mathcal{R}$ for $\psi_i$ and $q_i$, and $\mathcal{H}_{\text{per}}$ for $w_i$ and $v_i$, where $\boldsymbol{v}_i$, $q_i$, and $v_i$ are the test functions, as follows:

$$ \mathcal{U}_{\text{per}} := \{ \boldsymbol{v} \in H_{\text{per}}^{1}(Y, \mathbb{R}^d); Y = [0, 1]^2 \}, $$

$$ \mathcal{R} := \{ p \in L^2(Y); \mathcal{M}_Y(p) = 0; Y = [0, 1]^2 \}, $$

$$ \mathcal{H}_{\text{per}} := \{ u \in H_{\text{per}}^{1}(Y); \mathcal{M}_{Y}(u) = 0; Y = [0, 1]^2 \}. $$

The unit-cell problem for the Stokes flow defined in Eq. (16) can be rewritten in the weak form as:

Find $(\boldsymbol{X}_i, \psi_i, \lambda) \in \mathcal{U}_{\text{per}} \times \mathcal{R} \times \mathbb{R}$, such that:

$$ a(\boldsymbol{X}_i, \boldsymbol{v}_i) + b(\boldsymbol{v}_i, \psi_i) = l(\boldsymbol{v}_i) \quad \forall \boldsymbol{v}_i \in \mathcal{U}_{\text{per}}, $$

$$ b(\boldsymbol{X}_i, q_i) - \tau d(\psi_i, q_i) + \lambda m(q_i) = 0 \quad \forall q_i \in \mathcal{R}, $$

$$ m(\psi_i) = 0, $$

where we define the bilinear form $a(\boldsymbol{X}_i, \boldsymbol{v}_i)$ and $b(\boldsymbol{v}_i, \psi_i)$, linear form $l(\boldsymbol{v}_i)$, pressure stabilisation term $\tau d(\psi_i, q_i)$, with the element size $h_K$ and Lagrange multiplier $\lambda$, as follows:

$$ a(\boldsymbol{X}_i, \boldsymbol{v}_i) = \int_Y \nabla \boldsymbol{X}_i : \nabla \boldsymbol{v}_i \, dy + \int_Y \alpha(\boldsymbol{y}) \boldsymbol{X}_i \cdot \boldsymbol{v}_i \, dy, $$

$$ b(\boldsymbol{v}_i, \psi_i) = - \int_Y \nabla \cdot \boldsymbol{v}_i \psi_i \, dy $$

$$ l(\boldsymbol{v}_i) = \int_Y \boldsymbol{e}_i \cdot \boldsymbol{v}_i \, dy, $$

$$ \tau > 0, \quad d(\psi_i, q_i) = \sum_{K \in \mathcal{T}} h_K^2 \int_{Y_e} \nabla \psi_i \cdot \nabla q_i \, dy_e, $$

$$ m(\psi_i) = \int_Y \psi_i \, dy. $$

In our implementation, the discretisation of the above weak form is stored in the block matrix as follows:

$$ A_{\text{Stokes}} = \begin{bmatrix} \mathbf{A} & \mathbf{B} & \mathbf{0} \\ \mathbf{B}^T & \mathbf{S} & \mathbf{h} \\ \mathbf{0} & \mathbf{h}^T & \mathbf{0} \end{bmatrix}, \qquad (30) $$

where $\mathbf{A}$ is the matrix discretising the bilinear form of $a(\boldsymbol{X}_i, \boldsymbol{v}_i)$, $\mathbf{B}$ the divergence operator $b(\boldsymbol{v}_i, \psi_i)$, $\mathbf{S}$ the pressure stabilisation term $-\tau d(\psi_i, q_i)$. $\mathbf{h}$ is the discretisation of the linear form associated with the Lagrange multiplier term $\lambda m(q_i)$, imposing a zero-mean value for $\psi_i$ such that the solution can be uniquely determined.

The original system is indefinite because of the pressure being by default defined up to a constant. Since an indefinite linear system is generally more difficult than solving a definite linear system, we impose the zero-mean constraint. Since this can be dealt with efficiently with our approach from a discretisation and linear solver point of view, the extra work from adding a Lagrange multiplier is not a real
The weak form of the adjoint equations can be found in Appendix A.

Additionally, the weak form for the reaction-diffusion equation (PDE filter) defined in Eq. (18) is as follows:
Find $\tilde{d}_i \in \mathcal{P}$ such that:

$$ \int_D r^2 \nabla \tilde{d}_i \cdot \nabla v + \tilde{d}_i v \, d\Omega = \int_D d_i v \, d\Omega \quad \forall \delta v \in \mathcal{P}, \tilde{d}_i \in \mathcal{P}. \quad (33) $$

The homogenisation-based optimisation algorithm, summarised in Algorithm 1, is briefed as follows: it begins by initializing finite element spaces, declaring unknowns, and defining variational formulations, followed by linear interpolation of homogenised properties from an offline library described in Section 3.1. Within a main loop, the optimisation algorithm iterates until convergence or a maximum number of iterations, solving homogenisation equations to obtain velocity, pressure, and temperature fields, and computing objective and constraint functions. Adjoint equations are then solved to obtain adjoint variables, and sensitivities are computed. The design variables are updated using the method of moving asymptotes (MMA) (Svanberg, 1987).

**Algorithm 1 Homogenisation-based topology optimisation.**

**Input:** offline library: $\mathcal{X}_{11}(m_1, m_2)$, $\mathcal{X}_{22}(m_1, m_2)$, $\mathcal{K}_{11}(m_1, m_2)$, $\mathcal{K}_{22}(m_1, m_2)$, $\partial\mathcal{X}_{11}/\partial m_1$, $\partial\mathcal{X}_{11}/\partial m_2$, $\partial\mathcal{X}_{22}/\partial m_1$, $\partial\mathcal{X}_{22}/\partial m_2$, $\partial\mathcal{K}_{11}/\partial m_1$, $\partial\mathcal{K}_{11}/\partial m_2$, $\partial\mathcal{K}_{22}/\partial m_1$, $\partial\mathcal{K}_{22}/\partial m_2$

**Initialization**

create finite element spaces, declare unknowns, define variational formulations
initialise design fields $(m_1, m_2, \theta)$

linear interpolation of the homogenised properties using the offline library, cf.
Section 3.1

**for** it = 0; it < 2000; it = it + 1 **do**

solve homogenisation equation Eqs. (31) and (32) to obtain $p$ and $T$

compute $J$, $G_1$, and $G_2$, cf. Eq. (21)

**if** $\|J_{it+1} - J_{it}\| < \varepsilon, \|G_{it+1} - G_{it}\| < \varepsilon$ **then**

break

**else**

solve adjoint equations to obtain $p_a$ and $T_a$, cf. Eqs. (A.9) and (A.10)

compute sensitivities w.r.t. the projected design variables, $\delta_{\hat{m}}\mathcal{L}$, cf.

Eq. (A.7)

PDE-filter sensitivities w.r.t. the design variables, cf. Eq. (33)

update design variable $(m_1, m_2, \theta)$ using MMA

do the PDE-filtering, cf. Eq. (33), and projection, cf. Eq. (19), to get

$(\hat{m}_1, \hat{m}_2, \hat{\theta})$

linear interpolation of homogenised properties, cf. Section 3.1

**end if**

**end for**
burden. It should be noted that setting the pressure to zero in any corner would be another way to make the system definite, but a point-wise constraint is not well-defined in the finite element context (it has no variational formulation) and picking any corner implies a random choice (of a corner).

Furthermore, periodic boundary conditions are imposed within the finite element kernel which ensures that two nodes of the finite element space that represent matching periodic unknowns (on distinct faces) share the same identifier such that, once assembled, the resulting solution of the sparse linear system is indeed periodic.

The unit cell problem for the heat transfer defined in Eq. (18) can be rewritten in the weak form as:

Find $(w_i, \lambda) \in \mathcal{H}_{\text{per}} \times \mathbb{R}$, such that:

$$ 
\begin{array}{l}
a(w_i, v_i) + \lambda m(v_i) = l(v_i) \quad \forall v_i \in \mathcal{H}_{\text{per}}, \\
m(w_i) = 0,
\end{array}
$$

where we define the bilinear form $a(w_i, v_i)$, linear form $l(v_i)$, with the Lagrange multiplier $\lambda$, as follows:

$$ 
\begin{align*}
a(w_i, v_i) &= \int_Y \kappa(\mathbf{y}) \nabla w_i \cdot \nabla v_i \, dy, \\
l(v_i) &= \int_Y \kappa(\mathbf{y}) e_i \cdot \nabla v_i \, dy, \\
m(w_i) &= \int_Y w_i \, dy.
\end{align*}
$$

The discretisation of the above weak form is stored in the block matrix as follows:

$$ A_{\text{Laplacian}} = \begin{bmatrix} A & h \\ h^T & 0 \end{bmatrix}, \qquad (31) $$

where $A$ is the matrix discretising the bilinear form $a(w_i, v_i)$. $h$ is the discretisation of the linear form associated with the Lagrange multiplier term $\lambda m(v_i)$, imposing a zero-mean value for $w_i$, such that the solution can be uniquely determined.

## 3.2 Homogenisation-based topology optimisation

Now, we proceed with the homogenisation-based topology optimisation to solve the optimum design problem formulated in Section 2.4. Here, we define the Hilbert spaces $\mathcal{P}$ for the homogenised pressure $p$ and homogenised temperature $T$, their adjoint variables $p_a$ and $T_a$, as well as the test functions $q$ and $S$, as follows:

$$ \mathcal{P} := \{ \tilde{p} \in H^1(\Omega) \mid \tilde{p} = 0 \text{ on } \partial\Omega_p^D \}. $$

First, the Laplacian defined in Eq. (7) can be expressed in the following weak form: Find $p \in \mathcal{P}$ with $p = 0$ on $\Gamma_{\text{out}}$ such that:

$$ \int_D \mathbf{X}^* (\hat{\tilde{m}}_1, \hat{\tilde{m}}_2, \hat{\theta}) \nabla p \cdot \nabla q \, d\Omega = -\epsilon^{-2} \int_{\Gamma_{\text{in}}} (\mathbf{u}_0 \cdot \mathbf{n}) q \, d\Gamma \quad \forall q \in \mathcal{P}, p \in \mathcal{P}. \quad (32) $$
### 3.3 Reconstruction of fine-scale structure

The least-squares minimisation problem defined in Eq. (26) can be solved using the variational method, given as:
Find $\phi_1 \in \mathcal{P}$ such that:

$$ 
\begin{aligned}
& \int_{\Omega} a_1(\boldsymbol{x})\nabla\phi_1 \cdot \nabla v + a_2(\boldsymbol{x}) (\nabla\phi_1 \cdot \boldsymbol{e}_2(\hat{\theta})) (\nabla v \cdot \boldsymbol{e}_2(\hat{\theta})) + \tau\phi_1 v \, d\Omega \\
& = \int_{\Omega} a_1(\boldsymbol{x})\nabla v \cdot \boldsymbol{e}_1(\hat{\theta}) \, d\Omega \quad \forall v \in \mathcal{P}, \phi_1 \in \mathcal{P}.
\end{aligned}
\quad (34) $$

### 3.4 Validation for the upsampled structure

As previously noted, validating the upsampling outcomes requires conducting full-scale simulations. Here, we define the Hibert space $\mathcal{W}_{u,p}$ for $(u_\epsilon, p_\epsilon)$ and test functions $(v_\epsilon, q_\epsilon)$, and $\mathcal{P}$ for $T_\epsilon$ and test function $S_\epsilon$, as follows:

$$ 
\begin{aligned}
\mathcal{W}_{u,p} &:= \{(\tilde{u}, \tilde{p}) \in H^1(\Omega_f, \mathbb{R}^d) \times L^2(\Omega_f) \mid \tilde{u} = 0 \text{ on } \partial\Omega_f^D\}, \\
\mathcal{P} &:= \{\tilde{p} \in H^1(\Omega) \mid \tilde{p} = 0 \text{ on } \partial\Omega_p^D\}.
\end{aligned}
 $$

The weak form for the Stokes equations in Eq. (1) is as follows:

Find $(u_\epsilon, p_\epsilon) \in \mathcal{W}_{u,p}$ with $u_\epsilon = u_0$ on $\Gamma_{in}$ and $u_\epsilon = 0$ on $\partial\Omega_\epsilon \setminus (\Gamma_{in} \cup \Gamma_{out})$ such that:

$$ 
\begin{aligned}
& \int_{\Omega_{\epsilon}} \nabla u_{\epsilon} : \nabla v_{\epsilon} \, d\Omega - \int_{\Omega_{\epsilon}} \nabla \cdot v_{\epsilon} p_{\epsilon} \, d\Omega - \int_{\Omega_{\epsilon}} \nabla \cdot u_{\epsilon} q_{\epsilon} \, d\Omega \\
& + \int_{\Omega_{\epsilon}} \alpha_{\epsilon}(\boldsymbol{x}) u_{\epsilon} \cdot \boldsymbol{v}_{\epsilon} \, d\Omega - \sum_{K \in \mathcal{T}} h_K^2 \int_{\Omega_{\epsilon}} \nabla p_{\epsilon} \cdot \nabla q_{\epsilon} \, d\Omega_{\epsilon} = 0 \quad (35) \\
& \forall (v_{\epsilon}, q_{\epsilon}) \in \mathcal{W}_{u,p}, (u_{\epsilon}, p_{\epsilon}) \in \mathcal{W}_{u,p},
\end{aligned}
 $$

where the inverse permeability value in the solid phase is set to $\alpha_{max} = 10^9$.

The convection-diffusion equation in Eq. (8) can be written in weak form as follows:
Find $T_\epsilon \in \mathcal{P}$ with $T_\epsilon = 0$ on $\Gamma_{in}$ such that:

$$ \int_{\Omega_{\epsilon}} (u_{\epsilon} \cdot \nabla T_{\epsilon}) S_{\epsilon} \, d\Omega + \int_{\Omega_{\epsilon}} \mathcal{K}_{\epsilon} \nabla T_{\epsilon} \cdot \nabla S_{\epsilon} \, d\Omega = Q(\boldsymbol{x}) S_{\epsilon} \quad \forall S_{\epsilon} \in \mathcal{P}, T_{\epsilon} \in \mathcal{P}. \quad (36) $$

## 4 Numerical examples

In this section, we present two test cases: a symmetric design case in Section 4.1, and an asymmetric design case in Section 4.2. For both cases, we initialise the design field by setting $m_1 = 0.5$, $m_2 = 0.5$, and $\theta = 0.5$. We set the maximum allowed volume fraction for the fluid phase to $V_{max} = 60\%$ and the maximum allowed pressure drop coefficient to $C_{DP} = 1.0$, ensuring that the pressure drop for the optimised design does not exceed that for the initial guess. For the homogenisation-based topology optimisation, we use a 50 × 50 coarse mesh, while for the reconstruction of the de-homogenised designs, a 2000 × 2000 fine mesh is used. Based on the numerical examples presented later in this section, we verify that the narrowest microchannel spans at least a 5-element width. Increasing the mesh resolution could potentially enhance the accuracy of the simulations even further. The PDE-filter radius is set
Next, the convection-diffusion equation defined in Eq. (14) can be written in weak form as:

Find $T \in \mathcal{P}$ with $T = 0$ on $\Gamma_{in}$ such that:

$$ 
\begin{aligned}
& \int_D -\epsilon^2 \boldsymbol{\chi}^* (\hat{\tilde{m}}_1, \hat{\tilde{m}}_2, \hat{\theta}) \nabla p \cdot \nabla T S \, d\Omega + \int_D \boldsymbol{\kappa}^* (\hat{\tilde{m}}_1, \hat{\tilde{m}}_2, \hat{\theta}) \nabla T \cdot \nabla S \, d\Omega \\
& = \int_D Q (\hat{\tilde{m}}_1, \hat{\tilde{m}}_2) S \, d\Omega \quad \forall S \in \mathcal{P}, T \in \mathcal{P}.
\end{aligned}
\quad (33) $$

The weak form of the adjoint equations can be found in Appendix A.

Additionally, the weak form for the reaction-diffusion equation (PDE filter) defined in Eq. (19) is as follows:

Find $\tilde{d}_i \in \mathcal{P}$ such that:

$$ \int_D r^2 \nabla \tilde{d}_i \cdot \nabla v + \tilde{d}_i v \, d\Omega = \int_D d_i v \, d\Omega \quad \forall \delta v \in \mathcal{P}, \tilde{d}_i \in \mathcal{P}. \quad (34) $$

The homogenisation-based optimisation algorithm, summarised in Algorithm 1, is briefed as follows: it begins by initialising finite element spaces, declaring unknowns, and defining variational formulations, followed by linear interpolation of homogenised properties from an offline library described in Section 3.1. Within a main loop, the optimisation algorithm iterates until convergence or a maximum number of iterations, solving homogenisation equations to obtain velocity, pressure, and temperature fields, and computing objective and constraint functions. Adjoint equations are then solved to obtain adjoint variables, and sensitivities are computed. The design variables are updated using the method of moving asymptotes (MMA) (Svanberg, 1987).

## 3.3 Reconstruction of fine-scale structure

As discussed in Section 2.5, to reconstruct the fine-scale microchannels, first the optimised design variables $(\hat{\tilde{m}}_1, \hat{\tilde{m}}_2, \hat{\theta})$ are interpolated from the coarse mesh $T_{coarse}$, used during the homogenisation-based optimisation, to a higher resolution mesh $T_{fine}$. Then, the interpolated design variables are further smoothened by solving a reaction-diffusion equation on $T_{fine}$ with the filter radius set to the size of the fine mesh edges.

Next, the least-squares minimisation problem defined in Eq. (27) can be solved using the variational method, given as:

Find $\phi_1 \in \mathcal{P}$ such that:

$$ 
\begin{aligned}
& \int_{\Omega} a_1(\boldsymbol{x}) \nabla \phi_1 \cdot \nabla v + a_2(\boldsymbol{x}) (\nabla \phi_1 \cdot e_2(\hat{\theta})) (\nabla v \cdot e_2(\hat{\theta})) + \tau \phi_1 v \, d\Omega \\
& = \int_{\Omega} a_1(\boldsymbol{x}) \nabla v \cdot e_1(\hat{\theta}) \, d\Omega \quad \forall v \in \mathcal{P}, \phi_1 \in \mathcal{P}.
\end{aligned}
\quad (35) $$

Finally, using the distributions of $\hat{\tilde{m}}_1$ and $\hat{\tilde{m}}_2$ together with the transformation $\Phi$ described by the level-set functions $\phi_1$ and $\phi_2$, we can define the fluid domain on the curved grid with periodicity $\epsilon$ via the negative subdomain of pseudo-density functions given in Eq. (28).
to $r = 0.05$ which is $r \approx 2.4h$ where $h$ is the edge size of the mesh. The optimisation solver is terminated after 2000 iterations since suboptimal pin-fins take a long time to develop into straight connected fins.

Prior to conducting any optimisation, we first run the simulation solver to evaluate the de-homogenised initial designs across different values of $\epsilon \in \{0.05, 0.025, 0.0125\}$, as illustrated in Fig. 14 for the symmetric design case (see the design schematic in Fig. 8).

We utilise two different types of meshes in the simulation solver. Firstly, we solve the governing equations on a $2000 \times 2000$ mesh employing a Brinkmann term to represent the solid/fluid phases. The velocity, pressure, and temperature fields are illustrated in the first three columns of Fig. 15, with those obtained from the homogenised model presented in the fourth column for comparison.

Secondly, we generate a body-fitted mesh, as shown in Fig. 16. This mesh enables disjoint-reunion subdomains $\Omega_{\epsilon,f}$ and $\Omega_{\epsilon,s}$ divided by an explicit interface $\Gamma_{s,f}$. In this setup, we only need to solve the Stokes equation within the fluid domain $\Omega_{\epsilon,f}$ without incorporating any Brinkmann term. And the no-slip boundary condition $u_{\epsilon} = 0$ is imposed on the fluid-solid interface $\Gamma_{s,f}$. The velocity, pressure, and temperature fields are displayed in the first three columns of Fig. 17, while those obtained from the homogenised model are depicted in the fourth column for comparison.

We summarise the objective values obtained from the homogenised and dehomogenised structures in Table 2. The relative difference of the objective values obtained by the full-scale Stokes solution w.r.t. the homogenised Darcy solution is less than 1.0%. In the subsequent design examples, we will validate the performance of upsampled structures on body-fitted meshes.

<figure><img src="image_15.png" /><figcaption>Fig. 14: de-homogenised initial designs for various values of $\epsilon$ for the optimisation design case, where black regions represent solid material and white regions represent fluid.</figcaption></figure>

## 4.1 A symmetric test case

Now, we run Algorithm 1 to solve the symmetric test case illustrated in Fig. 8. After 2000 iterations, we obtain the optimised design fields showcased in Fig. 18. Notably, the $m_1$ field, as depicted in Fig. 18a, exhibits a uniform value of 1.0, suggesting streamline-oriented continuous fins in the design domain. Furthermore,
**Algorithm 1** Homogenisation-based topology optimisation.

**Input:** offline library: $\mathcal{X}_{11}(m_1, m_2)$, $\mathcal{X}_{22}(m_1, m_2)$, $\mathcal{K}_{11}(m_1, m_2)$, $\mathcal{K}_{22}(m_1, m_2)$, $\partial\mathcal{X}_{11}/\partial m_1$, $\partial\mathcal{X}_{11}/\partial m_2$, $\partial\mathcal{X}_{22}/\partial m_1$, $\partial\mathcal{X}_{22}/\partial m_2$, $\partial\mathcal{K}_{11}/\partial m_1$, $\partial\mathcal{K}_{11}/\partial m_2$, $\partial\mathcal{K}_{22}/\partial m_1$, $\partial\mathcal{K}_{22}/\partial m_2$

**Initialisation**

create finite element spaces, declare unknowns, define variational formulations
initialise design fields $(m_1, m_2, \theta)$

linear interpolation of the homogenised properties using the offline library, cf.
Section 3.1

for it = 0; it < 2000; it = it + 1 do

solve homogenisation equation Eqs. (32) and (33) to obtain $p$ and $T$

compute $J$, $G_1$, and $G_2$, cf. Eq. (22)

if $\|J_{it+1} - J_{it}\| < \varepsilon$, $\|G_{it+1} - G_{it}\| < \varepsilon$ then

break

else

solve adjoint equations to obtain $p_a$ and $T_a$, cf. Eqs. (A.9) and (A.10)

compute sensitivities w.r.t. the projected design variables, $\delta_{\hat{m}}\mathcal{L}$, cf.

Eq. (A.7)

PDE-filter sensitivities w.r.t. the design variables, cf. Eq. (34)

update design variable $(m_1, m_2, \theta)$ using MMA

do the PDE-filtering, cf. Eq. (34), and projection, cf. Eq. (20), to get

$(\hat{m}_1, \hat{m}_2, \hat{\theta})$

linear interpolation of homogenised properties, cf. Section 3.1

end if

end for

## 3.4 Validation for the upsampled structure

As previously noted, validating the upsampling outcomes requires conducting full-scale simulations. Here, we define the Hibert space $\mathcal{W}_{u,p}$ for $(u_\epsilon, p_\epsilon)$ and test functions $(v_\epsilon, q_\epsilon)$, and $\mathcal{P}$ for $T_\epsilon$ and test function $S_\epsilon$, as follows:

$$ \mathcal{W}_{u,p} := \{ (\tilde{u}, \tilde{p}) \in H^1(\Omega_f, \mathbb{R}^d) \times L^2(\Omega_f) \mid \tilde{u} = 0 \text{ on } \partial\Omega_f^D \}, $$

$$ \mathcal{P} := \{ \tilde{p} \in H^1(\Omega) \mid \tilde{p} = 0 \text{ on } \partial\Omega_p^D \}. $$

The weak form for the Stokes equations in Eq. (1) is as follows:

Find $(u_\epsilon, p_\epsilon) \in \mathcal{W}_{u,p}$ with $u_\epsilon = u_0$ on $\Gamma_{in}$ and $u_\epsilon = 0$ on $\partial\Omega_\epsilon \setminus (\Gamma_{in} \cup \Gamma_{out})$ such that:

$$ \int_{\Omega_{\epsilon}} \nabla u_{\epsilon} : \nabla v_{\epsilon} \, d\Omega - \int_{\Omega_{\epsilon}} \nabla \cdot v_{\epsilon} p_{\epsilon} \, d\Omega - \int_{\Omega_{\epsilon}} \nabla \cdot u_{\epsilon} q_{\epsilon} \, d\Omega \\ + \int_{\Omega_{\epsilon}} \alpha_{\epsilon}(x) u_{\epsilon} \cdot v_{\epsilon} \, d\Omega - \sum_{K \in \mathcal{T}} h_K^2 \int_{\Omega_{\epsilon}} \nabla p_{\epsilon} \cdot \nabla q_{\epsilon} \, d\Omega_{\epsilon} = 0 \quad (36) \\ \forall (v_{\epsilon}, q_{\epsilon}) \in \mathcal{W}_{u,p}, (u_{\epsilon}, p_{\epsilon}) \in \mathcal{W}_{u,p}, $$

where the inverse permeability value in the solid phase is set to $\alpha_{max} = 10^9$.

The convection-diffusion equation in Eq. (8) can be written in weak form as follows:
<figure><img src="image_16.png" /><figcaption>**Fig. 15:** The first three columns, from left to right, depict the velocity, pressure, and temperature fields obtained by solving the Stokes equations on a 2000 × 2000 fine mesh for the upsampled initial structures with various ε, cf. Fig. 14. The fourth column displays the homogenised velocity, pressure, and temperature fields acquired by solving the homogenised equation on a 50 × 50 coarse mesh. For comparison, the legend bar range is rescaled based on the homogenised model.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Fig. 16:** A snapshot of the body-fitted mesh representing the fluid subdomain Ω<sub>ε,f</sub>, derived from the initial guess where ε = 0.05.</figcaption></figure>

the optimized θ field, presented in Fig. 18c, illustrates that the fluid channel is
directed towards the far-end corners, while the channel network converges towards
Find $T_{\epsilon} \in \mathcal{P}$ with $T_{\epsilon} = 0$ on $\Gamma_{\text{in}}$ such that:

$$ \int_{\Omega_{\epsilon}} (\mathbf{u}_{\epsilon} \cdot \nabla T_{\epsilon}) S_{\epsilon} \, d\Omega + \int_{\Omega_{\epsilon}} \mathcal{K}_{\epsilon} \nabla T_{\epsilon} \cdot \nabla S_{\epsilon} \, d\Omega = Q(\mathbf{x}) S_{\epsilon} \quad \forall S_{\epsilon} \in \mathcal{P}, T_{\epsilon} \in \mathcal{P}. \quad (37) $$

# 4 Numerical examples

In this section, we present two test cases: a symmetric design case in Section 4.1, and an asymmetric design case in Section 4.2. For both cases, we initialise the design field by setting $m_1 = 0.5$, $m_2 = 0.5$, and $\theta = 0.5$. We set the maximum allowed volume fraction for the fluid phase to $V_{\text{max}} = 60\%$ and the maximum allowed pressure drop coefficient to $C_{\text{DP}} = 1.0$, ensuring that the pressure drop for the optimised design does not exceed that for the initial guess. For the homogenisation-based topology optimisation, we use a 50 × 50 coarse mesh, while for the reconstruction of the de-homogenised designs, a 2000 × 2000 fine mesh is used. Based on the numerical examples presented later in this section, we verify that the narrowest microchannel spans at least a 5-element width. Increasing the mesh resolution could potentially enhance the accuracy of the simulations even further. The PDE-filter radius is set to $r = 0.05$ which is $r \approx 2.4h$ where $h$ is the edge size of the mesh. The optimisation solver is terminated after 2000 iterations since suboptimal pin-fins take a long time to develop into straight connected fins.

When doing the homogenisation of the Stokes equation for a porous medium, we end up with a Darcy equation with an effective permeability matrix, which is scale-invariant. Therefore, the permeability matrix is multiplied by a scaling factor $\epsilon^2$ to account for the unit cell size (Auriault et al., 2010). In this way, the pressure field obtained from the full-scale Stokes equation directly corresponds to the realised physical domain and actual size. However, this would make it necessary to plot several pressure fields for the homogenised model, each associated with different unit cell sizes and different legend bar scales, making comparison difficult.

Alternatively, for the homogenised Darcy equation, one can set $\epsilon$ to 1 and scale the full-scale model $n$ times larger than the homogenised model, where the scale factor $n$ is equal to $1/\epsilon$. In this way, the pressure field obtained from the Stokes equation still directly corresponds to the realised physical domain. However, this would then be necessary to plot pressure fields for the full-scale model at different sizes and with different legend bar scales, which is also not ideal for comparison purposes.

For these reasons, $\epsilon$ is set to 1 for the homogenised model and the same computational domain size is used for both the homogenised and full-scale models. The full-scale solution of pressure is then rescaled by a multiplicative factor $\epsilon^2$. This way, only one pressure field for the homogenised model is needed, and it shares the same legend bar scale with all the full-scale solutions because they are normalised to have the same order of magnitude with the homogenised solution.

Prior to conducting any optimisation, we first run the simulation solver to evaluate the de-homogenised initial designs across different values of $\epsilon \in \{0.05, 0.025, 0.0125\}$, as illustrated in Fig. 14 for the symmetric design case (see the design schematic in Fig. 8).

We utilise two different types of meshes in the simulation solver. Firstly, we solve the governing equations on a 2000 × 2000 mesh employing a Brinkman term to
<figure><img src="image_18.png" /><figcaption>**Fig. 17:** The first three columns, from left to right, depict the velocity, pressure, and temperature fields obtained by solving the Stokes equations on body-fitted meshes for the upsampled initial structures with various $\epsilon$, cf. Fig. 14. The fourth column displays the homogenised velocity, pressure, and temperature fields acquired by solving the homogenised equation on a 50 × 50 coarse mesh. For comparison, the legend bar range is rescaled based on the homogenised model.</figcaption></figure>

<table><caption>**Table 2:** The objective values for the homogenised model and de-homogenised structures with various values of $\epsilon = \{0.05, 0.025, 0.125\}$ for the initial guess.</caption><thead><tr><th></th><th>Case</th><th>Objective value</th><th>Relative difference</th></tr></thead><tbody><tr><td>Homogenised, cf. Fig. 15i</td><td></td><td>1.89 · 10<sup>-2</sup></td><td>-</td></tr><tr><td rowspan="2">ε = 0.05</td><td>Brinkmann, cf. Fig. 15i</td><td>1.88 · 10<sup>-2</sup></td><td>0.74%</td></tr><tr><td>Body-fitted, cf. Fig. 17i</td><td>1.89 · 10<sup>-2</sup></td><td>0.29%</td></tr><tr><td rowspan="2">ε = 0.025</td><td>Brinkmann, cf. Fig. 15j</td><td>1.89 · 10<sup>-2</sup></td><td>0.02%</td></tr><tr><td>Body-fitted, cf. Fig. 17j</td><td>1.9 · 10<sup>-2</sup></td><td>0.30%</td></tr><tr><td rowspan="2">ε = 0.0125</td><td>Brinkmann, cf. Fig. 15k</td><td>1.88 · 10<sup>-2</sup></td><td>0.30%</td></tr><tr><td>Body-fitted, cf. Fig. 17k</td><td>1.88 · 10<sup>-2</sup></td><td>0.56%</td></tr></tbody></table>

the outlet. The convergence histories of the objective values, volume fraction of the fluid phase, and pressure drop constraint value are depicted in Fig. 19. The initial decrease/increase in cost functions is due to an infeasible initial guess and the MMA optimiser focuses on satisfying the volume constraint. After approximately 500 iterations, the functions largely stabilise, but the objective increases slightly
represent the solid/fluid phases. The velocity, pressure, and temperature fields are illustrated in the first three columns of Fig. 15, with those obtained from the homogenised model presented in the fourth column for comparison.

Secondly, we generate a body-fitted mesh, as shown in Fig. 16. This mesh enables disjoint-reunion subdomains $\Omega_{\epsilon,f}$ and $\Omega_{\epsilon,s}$ divided by an explicit interface $\Gamma_{s,f}$. In this setup, we only need to solve the Stokes equation within the fluid domain $\Omega_{\epsilon,f}$ without incorporating any Brinkman term. And the no-slip boundary condition $\mathbf{u}_{\epsilon} = \mathbf{0}$ is imposed on the fluid-solid interface $\Gamma_{s,f}$. The velocity, pressure, and temperature fields are displayed in the first three columns of Fig. 17, while those obtained from the homogenised model are depicted in the fourth column for comparison.

We summarise the objective values and normalised pressure drop values obtained from the homogenised and dehomogenised structures in Table 2, while the realised pressure drop values are provided in Table 3. It can be observed that smaller length scales $\epsilon$ lead to smaller pin-fins, which in turn create more boundary layers. These boundary layers contribute to increased viscous losses and, consequently, higher pressure drops. The relative difference of the objective values obtained by the full-scale Stokes solution with respect to the homogenised Darcy solution is less than 1.0%. However, there are noticeable deviations between the full-scale pressure fields and the homogenised pressure field. The homogenised model represents an averaged field and cannot fully resolve the local flow variations or capture the detailed boundary layer effects, which are more accurately captured in the full-scale simulations. In the subsequent design examples, we will validate the performance of upsampled structures on body-fitted meshes.

<figure><img src="image_21.png" /><figcaption>Fig. 14: de-homogenised initial designs for various values of $\epsilon$ for the optimisation design case, where black regions represent solid material and white regions represent fluid.</figcaption></figure>

## 4.1 A symmetric test case

Now, we run Algorithm 1 to solve the symmetric test case illustrated in Fig. 8. After 2000 iterations, we obtain the optimised design fields showcased in Fig. 18. Notably, the $m_1$ field, as depicted in Fig. 18a, exhibits a uniform value of 1.0, suggesting streamline-oriented continuous fins in the design domain. Furthermore, the optimised $\theta$ field, presented in Fig. 18c, illustrates that the fluid channel is directed towards the far-end corners, while the channel network converges towards
while the pressure drop becomes lower. Despite some slight oscillations due to its sensitivity to small changes in design fields, a trade-off is eventually reached.

<figure><img src="image_19.png" /><figcaption>Fig. 18: Optimised design fields obtained after 2000 iterations for the symmetric design case, cf. Fig. 8.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Fig. 19: Convergence histories of (a) objective value, (b) volume fraction of the fluid phase, and (c) pressure drop constraint value for the symmetric design case.</figcaption></figure>

Next, with the optimised orientation fields, we solve the least-squares minimisation problem for computing the level-set functions, as illustrated in Fig. 20. Then, we reconstruct the fine-scale microchannels using various $\epsilon$ values from the set
<figure><img src="image_22.png" /><figcaption>**Fig. 15:** The first three columns, from left to right, depict the velocity, pressure, and temperature fields obtained by solving the Stokes equations on a 2000 × 2000 fine mesh for the upsampled initial structures with various ε, cf. Fig. 14. The fourth column displays the homogenised velocity, pressure, and temperature fields acquired by solving the homogenised equation on a 50 × 50 coarse mesh. For comparison, the pressure fields obtained by the full-scale simulations are normalised by a multiplicative factor ε² and the legend bar range is rescaled based on the homogenised model.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Fig. 16:** A snapshot of the body-fitted mesh representing the fluid subdomain Ω<sub>ε,f</sub>, derived from the initial guess where ε = 0.05.</figcaption></figure>
{0.05, 0.025, 0.0125}, as depicted in Fig. 21. Following that, we validate these de-homogenised designs on the body-fitted meshes, showcasing the velocity, pressure, and temperature fields in the first three columns of Fig. 22, with corresponding results from the homogenised model presented in the fourth column for comparison. The computed objective values are summarized in Table 3.

The main observations can be summarised as follows:

1. In Fig. 21, it is evident that the de-homogenised designs feature some fluid channels connected to the wall but not connected to inlet/outlet, which is not ideal in practical applications. This phenomenon may be attributed to the absence of the no-slip boundary condition in the homogenised model, which is crucial for accurately modelling fluid motion. Consequently, in the homogenised model, we observe non-zero velocities at the boundary. This limitation warrants further investigation.

2. The thinner fins towards the right make sense from a heat transfer perspective. This might be explained by the fact that there is less solid material and, consequently, less heat generation near the outlet. Additionally, the fluid near the outlet will be hotter and less effective at cooling.

3. From a quantitative perspective, there is a notable discrepancy in the full-scale velocity magnitude compared to the homogenised one. However, such a deviation is expected as the homogenized model provides a locally averaged velocity, whereas the upsampled structure contains narrow microchannels with higher velocity magnitudes.

4. Upon computing the objective values (see Table 3), we also notice a slightly more noticeable discrepancy compared to that of the initial designs. The relative difference reaches up to 1.8%. This can be attributed to the varying microstructures within the design domain in the optimised design, unlike the initial designs which feature entirely periodic porous media.

By employing an intuitive post-processing strategy, we remove those microchannels in which the velocity magnitude is less than 0.1% of the maximum velocity magnitude within the entire domain based on the full-scale body-fitted simulation. This process effectively removes channels located at the far-end corners, as illustrated in Fig. 23. However, this may also lead to a larger discrepancy in the objective values. This is mainly due to the utilization of a design-dependent heat source. However, we argue that the proposed designs offer valuable insight for designers as they showcase a considerable improvement in heat transfer performance compared to the initial designs.

<table><caption>**Table 3:** The objective values for the homogenised design and de-homogenised designs with various values of $\epsilon = \{0.05, 0.025, 0.0125\}$, for the symmetric design case in Section 4.1.</caption><thead><tr><th>Case</th><th>Objective value</th><th>Relative difference</th></tr></thead><tbody><tr><td>Homogenised, cf. Fig. 18</td><td>2.98 &middot; 10<sup>-3</sup></td><td>-</td></tr><tr><td>&epsilon; = 0.05, cf. Fig. 21a</td><td>2.94 &middot; 10<sup>-3</sup></td><td>1.5%</td></tr><tr><td>&epsilon; = 0.025, cf. Fig. 21b</td><td>3 &middot; 10<sup>-3</sup></td><td>0.52%</td></tr><tr><td>&epsilon; = 0.0125, cf. Fig. 21c</td><td>2.93 &middot; 10<sup>-3</sup></td><td>1.8%</td></tr></tbody></table>
<figure><img src="image_24.png" /><figcaption>**Fig. 17:** The first three columns, from left to right, depict the velocity, pressure, and temperature fields obtained by solving the Stokes equations on body-fitted meshes for the upsampled initial structures with various $\epsilon$, cf. Fig. 14. The fourth column displays the homogenised velocity, pressure, and temperature fields acquired by solving the homogenised equation on a 50 × 50 coarse mesh. For comparison, the pressure fields obtained by the full-scale simulations are normalised by a multiplicative factor $\epsilon^2$ and the legend bar range is rescaled based on the homogenised model.</figcaption></figure>

<table><caption>**Table 2:** The objective values and normalised pressure drop values for the homogenised model and de-homogenised structures with various values of $\epsilon = \{0.05, 0.025, 0.125\}$ for the initial guess.</caption><thead><tr><th></th><th>Case</th><th>Objective value</th><th>Pressure drop (normalised)</th></tr></thead><tbody><tr><td>Homogenised, cf. Figs. 15 and 17 – 4th col.</td><td></td><td>1.89 · 10<sup>-2</sup></td><td>25.63</td></tr><tr><td rowspan="2">ε = 0.05</td><td>Brinkman, cf. Fig. 15 – 1st col.</td><td>1.88 · 10<sup>-2</sup></td><td>24.86</td></tr><tr><td>Body-fitted, cf. Fig. 17 – 1st col.</td><td>1.89 · 10<sup>-2</sup></td><td>26.09</td></tr><tr><td rowspan="2">ε = 0.025</td><td>Brinkman, cf. Fig. 15 – 2nd col.</td><td>1.89 · 10<sup>-2</sup></td><td>23.86</td></tr><tr><td>Body-fitted, cf. Fig. 17 – 2nd col.</td><td>1.9 · 10<sup>-2</sup></td><td>23.33</td></tr><tr><td rowspan="2">ε = 0.0125</td><td>Brinkman, cf. Fig. 15 – 3rd col.</td><td>1.88 · 10<sup>-2</sup></td><td>22.73</td></tr><tr><td>Body-fitted, cf. Fig. 17 – 3rd col.</td><td>1.88 · 10<sup>-2</sup></td><td>22.85</td></tr></tbody></table>

the outlet. The convergence histories of the objective values, volume fraction of the fluid phase, and pressure drop constraint value are depicted in Fig. 19. The initial decrease/increase in cost functions is due to an infeasible initial guess and the MMA optimiser focuses on satisfying the volume constraint. After approximately
<figure><img src="image_21.png" /><figcaption>**Fig. 20:** The level-set functions obtained by solving the least-squares minimisation problem for the symmetric design case.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**Fig. 21:** de-homogenised designs for various values of $\epsilon$ for the symmetric optimisation design case, cf. Fig. 8, where black regions represent solid material and white regions represent fluid.</figcaption></figure>

## 4.2 An asymmetric test case

Finally, we present results for an asymmetric design case illustrated in Fig. 24. After 2000 iterations, we obtain the optimised design fields as shown in Fig. 25. Similar to the observations made in the symmetric test case, within the $m_1$ field (as shown in Fig. 25a), regions with a value of 1.0 indicate a diagonally streamlined fluid channel design within the domain. Additionally, the optimized $\theta$ field (depicted in Fig. 25c) indicates that the fluid flow is directed towards the far-end corners after entering the inlet, while the channel network converges towards the outlet.

The level-set functions, obtained as the solution of the least-squares minimization problem, are displayed in Fig. 26. The upsampled microchannels with various $\epsilon$ values are illustrated in Fig. 27. The validation results on body-fitted meshes are presented in the first three columns of Fig. 28, while corresponding results from the homogenized model are provided in the fourth column for comparison. The computed objective values are summarized in Table 4. Higher deviations (relative differences) can be observed compared to the one in the symmetric test case. This may be due to the proximity of the inlet and outlet to the wall. The absence of the no-slip boundary condition in the homogenised model may exacerbate this deviation. Other than that, similar observations to those made in the symmetric
<table><caption>**Table 3:** The realised pressure drop values for the homogenised design and de-homogenised designs with various values of $\epsilon = \{0.05, 0.025, 0.0125\}$, for the initial guess, cf. Figs. 15 and 17.</caption><thead><tr><th>Case</th><th>Homogenised</th><th>Full-scale (Brinkman)</th><th>Full-scale (body-fitted)</th></tr></thead><tbody><tr><td>$\epsilon = 0.05$</td><td>$1.03 \cdot 10^4$</td><td>$1.04 \cdot 10^4$</td><td>$9.95 \cdot 10^3$</td></tr><tr><td>$\epsilon = 0.025$</td><td>$4.1 \cdot 10^4$</td><td>$3.73 \cdot 10^4$</td><td>$3.82 \cdot 10^4$</td></tr><tr><td>$\epsilon = 0.0125$</td><td>$1.64 \cdot 10^5$</td><td>$1.46 \cdot 10^5$</td><td>$1.45 \cdot 10^5$</td></tr></tbody></table>

500 iterations, the functions largely stabilise, but the objective increases slightly while the pressure drop becomes lower. Despite some slight oscillations due to its sensitivity to small changes in design fields, a trade-off is eventually reached.

Next, we further investigate the intermediate stages of the optimisation process by upscaling and presenting the corresponding snapshots of the upsampled microchannels in Fig. 20, which is called on-the-fly de-homogenisation. From the topology evolution history, it can be observed that at the initial stages e.g., $n_{it} = 10$, we observe a grid-like structure with sparse pin-fins. As the optimisation proceeds, particularly by e.g., $n_{it} = 50$, we notice that the channels start streamlining, becoming smoother and more aligned, though pin-fins still exist. By $n_{it} = 200$, the streamline structure of the channel has been significantly refined, and the pin-fins are gradually being removed, but some remain near the outlet. After 500 iterations, the microchannel design is almost completely streamlined with only minor residual structures. This long plateau in the optimisation process is likely due to the optimiser slowly eliminating these residual pin-fins, causing slight oscillations in the cost functions.

<figure><img src="image_25.png" /><figcaption>**Fig. 18:** Optimised design fields obtained after 2000 iterations for the symmetric design case, cf. Fig. 8.</figcaption></figure>

Table 4 shows a runtime breakdown for each finite element action. In addition, Fig. 21 provides a pie chart showing the percentage distribution for a clearer understanding. The total runtime for 2000 iterations was approximately 7 minutes and 7 seconds, which is manageable for two-dimensional test cases and does not pose a bottleneck. It should be also noted that the objective value itself does not affect the speed of the MMA optimiser in our case because we normalise the sensitivities with respect to the objective value. This prevents any scaling-related issues from slowing down the optimiser. The observed slow convergence is mainly due to the intended removal of residual small pin-fins, which requires quite a few iterations,
<figure><img src="image_23.png" /><figcaption>**Fig. 22:** The first three columns, from left to right, depict the velocity, pressure, and temperature fields obtained by solving the Stokes equations on a 2000 × 2000 fine mesh for the upsampled structures with various ε, cf. Fig. 21. The fourth column displays the homogenised velocity, pressure, and temperature fields acquired by solving the homogenisation equation on a 50 × 50 coarse mesh. For comparison, the legend bar range is rescaled based on the homogenised model.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>**Fig. 23:** Post-processed design achieved through the removal of minuscule channels (cf. Fig. 21) exhibiting very low velocity for the asymmetric design case.</figcaption></figure>

case outlined in Section 4.1 can be drawn from the results, highlighting both the favourable aspects and limitations of the proposed approach.
<figure><img src="image_26.png" /><figcaption>Fig. 19: Convergence histories of (a) objective value, (b) volume fraction of the fluid phase, and (c) pressure drop constraint value for the symmetric design case.</figcaption></figure>

particularly in the final stages of optimisation. As discussed, we ran the optimisation for 2000 iterations primarily to remove these residual features, which are small but take time to eliminate. In practice, the optimisation could be terminated much earlier.

<table><caption>Table 4: Time breakdown of the finite element actions for symmetric test case performed on Apple MacBook Pro with a 10-core M2 processor.</caption><thead><tr><th>Action</th><th>Runtime</th></tr></thead><tbody><tr><td>Solve Darcy equation</td><td>47 s</td></tr><tr><td>Solve energy equation</td><td>78 s</td></tr><tr><td>Solve adjoint energy equations</td><td>13 s</td></tr><tr><td>Solve adjoint Darcy equations</td><td>71 s</td></tr><tr><td>Assemble sensitivity</td><td>35 s</td></tr><tr><td>Update design variable using MMA</td><td>18 s</td></tr><tr><td>PDE filtering</td><td>72 s</td></tr><tr><td>Offline library interpolation</td><td>93 s</td></tr><tr><td>Total</td><td>7 min 7 s</td></tr></tbody></table>
<figure><img src="image_25.png" /><figcaption>**Fig. 24:** Schematic diagram of an asymmetric microchannel cooling design case.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>**Fig. 25:** Optimised design fields obtained after 2000 iterations for the asymmetric design case, cf. Fig. 24.</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>**Fig. 26:** The level-set functions obtained by solving the least-squares minimisation problem for the asymmetric design case.</figcaption></figure>

# 5 Conclusions

To conclude, the present work leverages recent advancements in de-homogenisation techniques to contribute to the multi-scale microchannel cooling designs. Beginning with the selection of a single-class microstructure—a square hole in this case—we employ numerical homogenisation to construct an offline library. This library feeds
<figure><img src="image_27.png" /><figcaption>**Fig. 20:** Snapshots of the microchannels upsampled at intermediate stages of the optimisation process for the symmetric design case. Black regions indicate solid material, while light blue regions represent the fluid. From left to right, microchannels are shown for different values of the unit cell size: $\epsilon = 0.05, 0.025, 0.0125$.</figcaption></figure>

Next, with the optimised orientation fields, we solve the least-squares minimisation problem for computing the level-set functions, as illustrated in Fig. 22. Then,
<figure><img src="image_28.png" /><figcaption>Fig. 27: de-homogenised designs for various values of $\epsilon$ for the asymmetric optimisation design case, cf. Fig. 24, where black regions represent solid material and white regions represent fluid.</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>Fig. 28: The first three columns, from left to right, depict the velocity, pressure, and temperature fields obtained by solving the Stokes equations on a 2000 × 2000 fine mesh for the upsampled structures with various $\epsilon$, cf. Fig. 27. The fourth column displays the homogenised velocity, pressure, and temperature fields acquired by solving the homogenisation equation on a 50 × 50 coarse mesh. For comparison, the legend bar range is rescaled based on the homogenised model.</figcaption></figure>

into online macro-scale topology optimisation, where both microstructure parameters and local orientation fields are optimised. Post-processing of the homogenisation results in the convergence towards optimal micro-fluid channels, capturing
<figure><img src="image_28.png" /><figcaption>**Fig. 21:** Pie chart illustrating the percentage distribution of runtime breakdown allocated to the finite element actions.</figcaption></figure>

we reconstruct the fine-scale microchannels using various $\epsilon$ values from the set $\{0.05, 0.025, 0.0125\}$, as depicted in Fig. 23. Following that, we validate these de-homogenised designs on the body-fitted meshes, showcasing the velocity, pressure, and temperature fields in the first three columns of Fig. 24, with corresponding results from the homogenised model presented in the fourth column for comparison. The computed objective values and normalised pressure drop values are summarised in Table 5.

The main observations can be summarised as follows:

1. In Fig. 23, it is evident that the de-homogenised designs feature some fluid channels connected to the wall but not connected to inlet/outlet, which is not ideal in practical applications. This phenomenon may be attributed to the absence of the no-slip boundary condition in the homogenised model, which is crucial for accurately modelling fluid motion. Consequently, in the homogenised model, we observe non-zero velocities at the boundary. This limitation warrants further investigation.

2. The thinner fins towards the right make sense from a heat transfer perspective. This might be explained by the fact that there is less solid material and, consequently, less heat generation near the outlet. Additionally, the fluid near the outlet will be hotter and less effective at cooling.

3. From a quantitative perspective, there is a notable discrepancy in the full-scale velocity and pressure drop magnitude compared to the homogenised one. However, such a deviation is expected as the homogenised model provides a locally averaged velocity, whereas the upsampled structure contains narrow microchannels with higher velocity magnitudes.

4. Upon computing the objective values (see Table 5), we also notice a slightly more noticeable discrepancy compared to that of the initial designs. The relative difference reaches up to 1.8%. This can be attributed to the varying microstructures within the design domain in the optimised design, unlike the initial designs which feature entirely periodic porous media.

By employing an intuitive post-processing strategy, we remove those microchannels in which the velocity magnitude is less than 0.1% of the maximum velocity magnitude within the entire domain based on the full-scale body-fitted simulation. This
<table><caption>**Table 4:** The objective values for the homogenised design and de-homogenised designs with various values of $\epsilon = \{0.05, 0.025, 0.0125\}$, for the asymmetric design case in Section 4.2.</caption><thead><tr><td>Case</td><td>Objective value</td><td>Relative difference</td></tr></thead><tbody><tr><td>Homogenised, cf. Fig. 25</td><td>4.02 · 10<sup>-3</sup></td><td>-</td></tr><tr><td>ε = 0.05, cf. Fig. 27a</td><td>4.22 · 10<sup>-3</sup></td><td>4.8%</td></tr><tr><td>ε = 0.025, cf. Fig. 27b</td><td>4.26 · 10<sup>-3</sup></td><td>5.9%</td></tr><tr><td>ε = 0.0125, cf. Fig. 27c</td><td>4.15 · 10<sup>-3</sup></td><td>3.1%</td></tr></tbody></table>

<figure><img src="image_30.png" /><figcaption>**Fig. 29:** Post-processed design achieved through the removal of minuscule channels (cf. Fig. 27) exhibiting very low velocity for the asymmetric design case.</figcaption></figure>

fine details through a unique computation for all. The de-homogenised outcomes showcase streamlined designs that are different from periodic arrangements of micro-pillar arrays exhibiting locally varying shapes. Through full-scale simulations with upsampled structures, we demonstrate that these generated microchannel cooling designs can outperform conventional micro-pillar arrays, providing valuable insights for heat sink designers. Additionally, imperfections observed in these de-homogenised designs serve as benchmarks for future improvements, addressing issues such as modelling accuracy, manufacturability of de-homogenised designs, and overall performance enhancements.

We note that although only $m_1 = 1$ solutions are observed in the presented examples, the methodology is capable of producing pin-fins if the underlying physics requires it. This renders the methodology ready for problems where breaking the boundary layer is captured in the physics and beneficial to the cooling performance. But the current model does not capture this due to the instantaneous development of the boundary layer in Stokes flow.

Below are the primary limitations of the current work, along with prospects aimed at addressing these challenging issues:

1. The homogenised model for the fluid flow essentially overlooks the boundary layers near the wall and disregards the transition to the Stokes regime in regions with large portion of fluid phase. Introducing a high-order homogenisation of Stokes could potentially enhance the accuracy of the homogenized model by integrating the Stokes system, the Brinkman model, and Darcy's law (Feppon, 2021).
process effectively removes channels located at the far-end corners, as illustrated in Fig. 25. However, this may also lead to a larger discrepancy in the objective values. This is mainly due to the utilisation of a design-dependent heat source. However, we argue that the proposed designs offer valuable insight for designers as they showcase a considerable improvement in heat transfer performance compared to the initial designs.

<figure><img src="image_29.png" /><figcaption>**Fig. 22:** The level-set functions obtained by solving the least-squares minimisation problem for the symmetric design case.</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>**Fig. 23:** de-homogenised designs for various values of $\epsilon$ for the symmetric optimisation design case, cf. Fig. 8, where black regions represent solid material and white regions represent fluid.</figcaption></figure>

<table><caption>**Table 5:** The objective values and normalised pressure drop values for the homogenised design and de-homogenised designs with various values of $\epsilon = \{0.05, 0.025, 0.0125\}$, for the symmetric design case in Section 4.1.</caption><thead><tr><th>Case</th><th>Objective value</th><th>Pressure drop (normalised)</th></tr></thead><tbody><tr><td>Homogenised, cf. Fig. 18</td><td>2.98 · 10<sup>-3</sup></td><td>25.61</td></tr><tr><td>ε = 0.05, cf. Fig. 23a</td><td>2.94 · 10<sup>-3</sup></td><td>18.56</td></tr><tr><td>ε = 0.025, cf. Fig. 23b</td><td>3 · 10<sup>-3</sup></td><td>17.69</td></tr><tr><td>ε = 0.0125, cf. Fig. 23c</td><td>2.93 · 10<sup>-3</sup></td><td>16.94</td></tr></tbody></table>
2. The de-homogenised designs feature some fluid channels which are not connected to those belong to inlet or outlet. The straightforward post-processing strategy demonstrated in this paper is not satisfying to resolve this issue. A potential solution is to employ a phasor noise de-homogenisation technique where an image morphology-based method was proposed to close these branches, ensuring structural connectivity while minimising the deviations in prescribed orientations, periodicity, and homogenised properties (Woldseth et al., 2024).

3. In this paper, we limit ourselves to use the simple square hole as the pre-selected microstructure. More diversified parametrised microstructure such as bar groups (Wang et al., 2022) and super-shapes (Padhy et al., 2023b) may provide a broader spectrum of properties, which is a focus of our future work.

**Supplementary Information**

No supplementary information is used in this paper.

**Acknowledgments**

The first and last authors gratefully acknowledge financial support from the European Union through the Marie Skłodowska-Curie Actions Postdoc Fellowship (no. 101106842). The first author also extends gratitude to Dr. Minghao Yu for his contribution to integrate the MMA algorithm into the workflow and for fruitful discussions.

**Declarations**

**Conflict of interest**

The authors declare that they have no competing financial interests for this paper.

**Replication of results**

A detailed implementation and algorithm of the proposed method have been presented in Section 3 and one can follow them and reproduce the results. In case of further queries, please contact the corresponding author(s).

**Author contributions**

Hao Li (hli@sdu.dk): Conceptualization, Methodology, Software, Validation, Writing – original draft. Pierre Jolivet (pierre@joliv.et): Software, Validation, Supervision, Writing – Review & Editing. Joe Alexandersen (joal@sdu.dk): Conceptualization, Methodology, Supervision, Writing – Review & Editing.

**References**

Ahmed Mohammed Adham, Normah Mohd-Ghazali and Robiah Ahmad. Thermal and hydrodynamic analysis of microchannel heat sinks: A review. *Renewable and Sustainable Energy Reviews*, 21:614–622, 2013.

Joe Alexandersen and Casper Schousboe Andreasen. A review of topology optimisation for fluid-based problems. *Fluids*, 5(1):29, 2020.
<figure><img src="image_31.png" /><figcaption>**Fig. 24:** The first three columns, from left to right, depict the velocity, pressure, and temperature fields obtained by solving the Stokes equations on a 2000 × 2000 fine mesh for the upsampled structures with various $\epsilon$, cf. Fig. 23. The fourth column displays the homogenised velocity, pressure, and temperature fields acquired by solving the homogenisation equation on a 50 × 50 coarse mesh. For comparison, the pressure fields obtained by the full-scale simulations are normalised by a multiplicative factor $\epsilon^2$ and the legend bar range is rescaled based on the homogenised model.</figcaption></figure>

<figure><img src="image_32.png" /><figcaption>**Fig. 25:** Post-processed design achieved through the removal of minuscule channels (cf. Fig. 23) exhibiting very low velocity for the asymmetric design case.</figcaption></figure>
Grégoire Allaire and Zakaria Habibi. Homogenization of a conductive, convective, and radiative heat transfer problem in a heterogeneous domain. *SIAM Journal on Mathematical Analysis*, 45(3):1136–1178, 2013.

Grégoire Allaire, François Jouve and Anca-Maria Toader. Structural optimization using sensitivity analysis and a level-set method. *Journal of Computational Physics*, 194(1):363–393, 2004.

Grégoire Allaire, Lorenzo Cavallina, Nobuhito Miyake, Tomoyuki Oka and Toshiaki Yachimura. The homogenization method for topology optimization of structures: old and new. *Interdisciplinary Information Sciences*, 25(2):75–146, 2019.

Grégoire Allaire, Perle Geoffroy-Donders and Olivier Pantz. Topology optimization of modulated and oriented periodic microstructures by the homogenization method. *Computers & Mathematics with Applications*, 78(7):2197–2229, 2019.

Erik Andreassen and Casper Schousboe Andreasen. How to determine composite material properties using numerical homogenization. *Computational Materials Science*, 83:488–495, 2014.

Satish Balay, Shrirang Abhyankar, Mark F. Adams, Steven Benson, Jed Brown, Peter Brune, Kris Buschelman, Emil Constantinescu, Lisandro Dalcin, Alp Dener, Victor Eijkhout, Jacob Faibussowitsch, William D. Gropp, Václav Hapla, Tobin Isaac, Pierre Jolivet, Dmitry Karpeev, Dinesh Kaushik, Matthew G. Knepley, Fande Kong, Scott Kruger, Dave A. May, Lois Curfman McInnes, Richard Tran Mills, Lawrence Mitchell, Todd Munson, Jose E. Roman, Karl Rupp, Patrick Sanan, Jason Sarich, Barry F. Smith, Stefano Zampini, Hong Zhang, Hong Zhang, and Junchao Zhang. PETSc/TAO users manual. Technical Report ANL-21/39 - Revision 3.19, Argonne National Laboratory, 2023.

Martin Philip Bendsøe and Noboru Kikuchi. Generating optimal topologies in structural design using a homogenization method. *Computer Methods in Applied Mechanics and Engineering*, 71(2):197–224, 1988.

Martin P Bendsøe and Ole Sigmund. Material interpolation schemes in topology optimization. *Archive of Applied Mechanics*, 69(9-10):635–654, 1999.

Thomas Borrvall and Joakim Petersson. Topology optimization of fluids in stokes flow. *International journal for numerical methods in fluids*, 41(1):77–107, 2003.

Ercan M Dede, Yuqing Zhou and Tsuyoshi Nomura. Inverse design of microchannel fluid flow networks using turing pattern dehomogenization. *Structural and Multidisciplinary Optimization*, 62(4):2203–2210, 2020.

Ercan M Dede, Yuqing Zhou, Tomoki Tambo, Feng Zhou, Danny J Lohan and Tsuyoshi Nomura. Measurement of low reynolds number flow emanating from a turing pattern microchannel array using a modified bernoulli equation technique. *Experimental Thermal and Fluid Science*, 139:110722, 2022.

Martin Ohrt Elingaard, Niels Aage, Jakob Andreas Bærentzen and Ole Sigmund. De-homogenization using convolutional neural networks. *Computer Methods in Applied Mechanics and Engineering*, 388:114197, 2022.

Florian Feppon. High order homogenization of the stokes system in a periodic porous medium. *SIAM Journal on Mathematical Analysis*, 53(3):2890–2924, 2021.

Florian Feppon. Asymptotic expansions of Stokes flows in finite periodic channels. working paper or preprint, February 2024.

Florian Feppon. Multiscale topology optimization of modulated fluid microchannels based on asymptotic homogenization. *Computer Methods in Applied Mechanics and Engineering*, 419:116646, 2024.
## 4.2 An asymmetric test case

Finally, we present results for an asymmetric design case illustrated in Fig. 26. After 2000 iterations, we obtain the optimised design fields as shown in Fig. 27. Similar to the observations made in the symmetric test case, within the $m_1$ field (as shown in Fig. 27a), regions with a value of 1.0 indicate a diagonally streamlined fluid channel design within the domain. Additionally, the optimised $\theta$ field (depicted in Fig. 27c) indicates that the fluid flow is directed towards the far-end corners after entering the inlet, while the channel network converges towards the outlet.

The level-set functions, obtained as the solution of the least-squares minimisation problem, are displayed in Fig. 28. The upsampled microchannels with various $\epsilon$ values are illustrated in Fig. 29. The validation results on body-fitted meshes are presented in the first three columns of Fig. 30, while corresponding results from the homogenised model are provided in the fourth column for comparison. The computed objective values and normalised pressure drop values are summarised in Table 6. Higher deviations (relative differences) can be observed compared to the one in the symmetric test case. This may be due to the proximity of the inlet and outlet to the wall. The absence of the no-slip boundary condition in the homogenised model may exacerbate this deviation. Other than that, similar observations to those made in the symmetric case outlined in Section 4.1 can be drawn from the results, highlighting both the favourable aspects and limitations of the proposed approach.

<figure><img src="image_33.png" /><figcaption>**Fig. 26:** Schematic diagram of an asymmetric microchannel cooling design case.</figcaption></figure>

<table><caption>**Table 6:** The objective values and normalised pressure drop values for the homogenised design and de-homogenised designs with various values of $\epsilon = \{0.05, 0.025, 0.0125\}$, for the asymmetric design case in Section 4.2.</caption><thead><tr><td>Case</td><td>Objective value</td><td>Pressure drop (normalised)</td></tr></thead><tbody><tr><td>Homogenised, cf. Fig. 27</td><td>4.02 · 10<sup>-3</sup></td><td>34.15</td></tr><tr><td>ε = 0.05, cf. Fig. 29a</td><td>4.22 · 10<sup>-3</sup></td><td>25.81</td></tr><tr><td>ε = 0.025, cf. Fig. 29b</td><td>4.26 · 10<sup>-3</sup></td><td>23.38</td></tr><tr><td>ε = 0.0125, cf. Fig. 29c</td><td>4.15 · 10<sup>-3</sup></td><td>22.35</td></tr></tbody></table>
Da Geng, Chuang Wei, Yichang Liu and Mingdong Zhou. Concurrent topology optimization of multi-scale cooling channels with inlets and outlets. *Structural and Multidisciplinary Optimization*, 65(8):234–2022.

Perle Geoffroy-Donders, Grégoire Allaire and Olivier Pantz. 3-d topology optimization of modulated and oriented periodic microstructures by the homogenization method. *Journal of Computational Physics*, 401:108994–2020.

Jeroen P Groen and Ole Sigmund. Homogenization-based topology optimization for high-resolution manufacturable microstructures. *International Journal for Numerical Methods in Engineering*, 113(8):1148–1163, 2018.

Jeroen P Groen, Jun Wu and Ole Sigmund. Homogenization-based stiffness optimization and projection of 2d coated structures with orthotropic infill. *Computer Methods in Applied Mechanics and Engineering*, 349:722–742, 2019.

Jeroen P Groen, Florian C Stutz, Niels Aage, Jakob A Bærentzen and Ole Sigmund. De-homogenization of optimal multi-scale 3d topologies. *Computer Methods in Applied Mechanics and Engineering*, 364:112979, 2020.

James K Guest and Jean H Prévost. Design of maximum permeability material structures. *Computer Methods in Applied Mechanics and Engineering*, 196(4-6):1006–1017, 2007.

Sarah N Hankins, Yuqing Zhou, Danny J Lohan and Ercan M Dede. Generative design of large-scale fluid flow structures via steady-state diffusion-based dehomogenization. *Scientific Reports*, 13(1):14344, 2023.

Frédéric Hecht. New development in FreeFem++. *Journal of Numerical Mathematics*, 20(3-4):251–266, 2012.

Peter Dørffler Ladegaard Jensen, Ole Sigmund and Jeroen P Groen. De-homogenization of optimal 2d topologies for multiple loading cases. *Computer Methods in Applied Mechanics and Engineering*, 399:115426, 2022.

Pierre Jolivet, Victorita Dolean, Frédéric Hecht, Frédéric Nataf, Christophe Prud'homme and Nicole Spillane. High-Performance Domain Decomposition Methods on Massively Parallel Architectures with FreeFem++. *Journal of Numerical Mathematics*, 20(4):287–302, 2012.

Atsushi Kawamoto, Tadayoshi Matsumori, Shintaro Yamasaki, Tsuyoshi Nomura, Tsuguo Kondoh and Shinji Nishiwaki. Heaviside projection based topology optimization by a pde-filtered scalar function. *Structural and Multidisciplinary Optimization*, 44(1):19–24, 2011.

Boyan Stefanov Lazarov and Ole Sigmund. Filters in topology optimization based on helmholtz-type differential equations. *International Journal for Numerical Methods in Engineering*, 86(6):765–781, 2011.

Jaewook Lee, Chiyoung Kwon, Jeonghoon Yoo, Seungjae Min, Tsuyoshi Nomura and Ercan M Dede. Design of spatially-varying orthotropic infill structures using multiscale topology optimization and explicit de-homogenization. *Additive Manufacturing*, 40:101920, 2021.

Doksoo Lee, Wei Chen, Liwei Wang, Yu-Chin Chan and Wei Chen. Data-driven design for metamaterials and multiscale systems: A review. *Advanced Materials*, 36(8):2305254, 2024.

HA Mohammed, G Bhaskaran, NH Shuaib and Rahman Saidur. Heat transfer and fluid flow characteristics in microchannels heat exchanger using nanofluids: a review. *Renewable and Sustainable Energy Reviews*, 15(3):1502–1512, 2011.

Serdar Ozguc, Liang Pan and Justin A Weibel. Topology optimization of microchannel heat sinks using a homogenization approach. *International Journal of Heat and Mass Transfer*, 169:120896, 2021.
<figure><img src="image_34.png" /><figcaption>Fig. 27: Optimised design fields obtained after 2000 iterations for the asymmetric design case, cf. Fig. 26.</figcaption></figure>

<figure><img src="image_35.png" /><figcaption>Fig. 28: The level-set functions obtained by solving the least-squares minimisation problem for the asymmetric design case.</figcaption></figure>

<figure><img src="image_36.png" /><figcaption>Fig. 29: de-homogenised designs for various values of $\epsilon$ for the asymmetric optimisation design case, cf. Fig. 26, where black regions represent solid material and white regions represent fluid.</figcaption></figure>

## 5 Conclusions

To conclude, the present work leverages recent advancements in de-homogenisation techniques to contribute to the multi-scale microchannel cooling designs. Beginning with the selection of a single-class microstructure—a square hole in this case—we employ numerical homogenisation to construct an offline library. This library feeds into online macro-scale topology optimisation, where both microstructure parameters and local orientation fields are optimised. Post-processing of the homogenisation results in the convergence towards optimal micro-fluid channels, capturing
Rahul Kumar Padhy, Aaditya Chandrasekhar and Krishnan Suresh. Fluto: Graded multi-scale topology optimization of large contact area fluid-flow devices using neural networks. *Engineering with Computers*, pages 1-17, 2023.

Rahul Kumar Padhy, Krishnan Suresh and Aaditya Chandrasekhar. Tomas: Topology optimization of multiscale fluid devices using variational autoencoders and super-shapes. *arXiv preprint arXiv:2309.08435*, 2023.

Olivier Pantz and Karim Trabelsi. A post-treatment of the homogenization method for shape optimization. *SIAM Journal on Control and Optimization*, 47(3):1380-1398, 2008.

Nicolò Pollini, Ole Sigmund, Casper Schousboe Andreasen and Joe Alexandersen. A "poor man's" approach for high-resolution three-dimensional topology design for natural convection problems. *Advances in Engineering Software*, 140:102736, 2020.

George IN Rozvany, Ming Zhou and Torben Birker. Generalized shape optimization without homogenization. *Structural Optimization*, 4(3-4):250-252, 1992.

Enrique Sanchez-Palencia. Fluid flow in porous media. *Non-homogeneous media and vibration theory*, pages 129-157, 1980.

Krister Svanberg. The method of moving asymptotes—a new method for structural optimization. *International journal for numerical methods in engineering*, 24(2):359-373, 1987.

Akihiro Takezawa, Xiaopeng Zhang, Masaki Kato and Mitsuru Kitamura. Method to optimize an additively-manufactured functionally-graded lattice structure for effective liquid cooling. *Additive Manufacturing*, 28:285-298, 2019.

Akihiro Takezawa, Xiaopeng Zhang and Mitsuru Kitamura. Optimization of an additively manufactured functionally graded lattice structure with liquid cooling considering structural performances. *International Journal of Heat and Mass Transfer*, 143:118564, 2019.

Michael Yu Wang, Xiaoming Wang and Dongming Guo. A level set method for structural topology optimization. *Computer Methods in Applied Mechanics and Engineering*, 192(1-2):227-246, 2003.

Liwei Wang, Zhao Liu, Daicong Da, Yu-Chin Chan, Wei Chen and Ping Zhu. Generalized de-homogenization via sawtooth-function-based mapping and its demonstration on data-driven frequency response optimization. *Computer Methods in Applied Mechanics and Engineering*, 395:114967, 2022.

Junpeng Wang, Rüdiger Westermann and Jun Wu. A streamline-guided dehomogenization approach for structural design. *Journal of Mechanical Design*, 145(2), 2023.

Rebekka V Woldseth, J Andreas Bærentzen and Ole Sigmund. Phasor noise for dehomogenisation in 2d multiscale topology optimisation. *Computer Methods in Applied Mechanics and Engineering*, 418:116551, 2024.

Jun Wu, Ole Sigmund and Jeroen P Groen. Topology optimization of multi-scale structures: a review. *Structural and Multidisciplinary Optimization*, 63(3):1455-1480, 2021.

Tong Wu. *Topology optimization of multiscale structures coupling fluid, thermal and mechanical analysis*. PhD thesis, Purdue University, 2019.

Xi Zhao, Mingdong Zhou, Ole Sigmund and Casper Schousboe Andreasen. A "poor man's approach" to topology optimization of cooling channels based on a darcy flow model. *International Journal of Heat and Mass Transfer*, 116:1108-1123, 2018.
<figure><img src="image_37.png" /><figcaption>**Fig. 30:** The first three columns, from left to right, depict the velocity, pressure, and temperature fields obtained by solving the Stokes equations on a 2000 × 2000 fine mesh for the upsampled structures with various $ϵ$, cf. Fig. 29. The fourth column displays the homogenised velocity, pressure, and temperature fields acquired by solving the homogenisation equation on a 50 × 50 coarse mesh. For comparison, the pressure fields obtained by the full-scale simulations are normalised by a multiplicative factor $ϵ^2$ and the legend bar range is rescaled based on the homogenised model.</figcaption></figure>

<figure><img src="image_38.png" /><figcaption>**Fig. 31:** Post-processed design achieved through the removal of minuscule channels (cf. Fig. 29) exhibiting very low velocity for the asymmetric design case.</figcaption></figure>

fine details through a unique computation for all. The de-homogenised outcomes showcase streamlined designs that are different from periodic arrangements of
Ming Zhou and George IN Rozvany. The coc algorithm, part ii: Topological, geometrical and generalized shape optimization. *Computer Methods in Applied Mechanics and Engineering*, 89(1-3):309–336, 1991.

Yuqing Zhou, Danny J Lohan, Feng Zhou, Tsuyoshi Nomura and Ercan M Dede. Inverse design of microreactor flow fields through anisotropic porous media optimization and dehomogenization. *Chemical Engineering Journal*, 435:134587, 2022.

## A Sensitivity analysis

The sensitivity analysis is conducted to derive the sensitivity of objective function $J$, and pressure drop constraint function $G_1$ w.r.t. the projected design variable $\hat{\mathbf{d}} = (\hat{\mathbf{m}}_1, \hat{\mathbf{m}}_2, \hat{\mathbf{\theta}})$. For conciseness, throughout the remainder of this appendix, we utilize $\mathbf{X}^*$ and $\mathbf{K}^*$ for brevity.

Now, the sensitivity analysis is performed. First, the Lagrange function $\mathcal{L}$ is constructed as follows:

$$ \mathcal{L} := J + \langle p_a, -\operatorname{div}(\mathbf{X}^*\nabla p) \rangle + \langle T_a, -\epsilon^2 \mathbf{X}^*\nabla p \cdot \nabla T - \operatorname{div}(\mathbf{K}^*\nabla T) - Q \rangle, \quad (A.1) $$

where $p_a$ and $T_a$ are the adjoint pressure and adjoint temperature, respectively. According to the Karush-Kuhn-Tucker (KKT) conditions for PDE constrained optimisation problems,

$$ \left\{ \begin{array}{l} \delta_T \mathcal{L} = 0, \\ \delta_p \mathcal{L} = 0. \end{array} \right. \qquad (A.2) $$

The variation of the Lagrange function w.r.t. the state variables $T$ and $p$ are derived as follows:

$$ \delta_T \mathcal{L} = \delta_T J + \int_D -\epsilon^2 \mathbf{X}^* \nabla p \cdot \nabla \delta T T_a \, d\Omega + \int_D \mathbf{K}^* \nabla \delta T \cdot \nabla T_a \, d\Omega = 0 \quad \forall \delta T \in \mathcal{P}, T_a \in \mathcal{P}, \qquad (A.3a) $$

$$ \delta_p \mathcal{L} = \delta_p J + \int_D \mathbf{X}^* \nabla p_a \cdot \nabla \delta p \, d\Omega - \int_D \epsilon^2 \mathbf{X}^* \nabla \delta p \cdot \nabla T T_a = 0 \quad \forall \delta p \in \mathcal{P}, p_a \in \mathcal{P}. \qquad (A.3b) $$

As for the optimum design problem given in Eq. (21), the variation of the cost functions w.r.t. the state variables are derived as follows:

$$ \left\{ \begin{array}{l} \delta_T J = \int_D \hat{\mathbf{m}}_1 \hat{\mathbf{m}}_2 \delta T \, d\Omega, \\ \delta_p J = 0, \end{array} \right. \qquad (A.4a) $$

$$ \left\{ \begin{array}{l} \delta_T G_1 = 0, \\ \delta_p G_1 = \epsilon^2 \frac{1/\text{DP}_0}{|\Gamma_{\text{in}}|} \int_{\Gamma_{\text{in}}} \delta p \, d\Gamma. \end{array} \right. \qquad (A.4b) $$

Substituting Eqs. (A.3) and (A.4) into Eq. (A.2), the adjoint system can be obtained, see the corresponding weak forms Eqs. (A.9) and (A.10).
micro-pillar arrays exhibiting locally varying shapes. Through full-scale simulations with upsampled structures, we demonstrate that these generated microchannel cooling designs can outperform conventional micro-pillar arrays, providing valuable insights for heat sink designers. Additionally, imperfections observed in these de-homogenised designs serve as benchmarks for future improvements, addressing issues such as modelling accuracy, manufacturability of de-homogenised designs, and overall performance enhancements.

We note that although only $m_1 = 1$ solutions are observed in the presented examples, the methodology is capable of producing pin-fins if the underlying physics requires it. This renders the methodology ready for problems where breaking the boundary layer is captured in the physics and beneficial to the cooling performance. But the current model does not capture this due to the instantaneous development of the boundary layer in Stokes flow.

Below are the primary limitations of the current work, along with prospects aimed at addressing these challenging issues:

1. The homogenised model for the fluid flow essentially overlooks the boundary layers near the wall and disregards the transition to the Stokes regime in regions with large portion of fluid phase. Introducing a high-order homogenisation of Stokes could potentially enhance the accuracy of the homogenised model by integrating the Stokes system, the Brinkman model, and Darcy's law (Feppon, 2021).

2. The de-homogenised designs feature some fluid channels which are not connected to those belong to inlet or outlet. The straightforward post-processing strategy demonstrated in this paper is not satisfying to resolve this issue. A potential solution is to employ a phasor noise de-homogenisation technique where an image morphology-based method was proposed to close these branches, ensuring structural connectivity while minimising the deviations in prescribed orientations, periodicity, and homogenised properties (Woldseth et al., 2024).

3. In this paper, we limit ourselves to use the simple square hole as the pre-selected microstructure. More diversified parametrised microstructure such as bar groups (Wang et al., 2022) and super-shapes (Padhy et al., 2023b) may provide a broader spectrum of properties, which is a focus of our future work.

**Supplementary Information**

No supplementary information is used in this paper.

**Acknowledgments**

The first and last authors gratefully acknowledge financial support from the European Union through the Marie Skłodowska-Curie Actions Postdoc Fellowship (no. 101106842). The first author also extends gratitude to Dr. Minghao Yu for his contribution to integrate the MMA algorithm into the workflow and for fruitful discussions.
Finally, the functional derivative of the Lagrange function w.r.t. design variables, $\delta_d \mathcal{L}$, can be derived using the chain rule as follows:

$$ \left\{ \begin{array}{l} \frac{\delta \mathcal{L}}{\delta \tilde{m}_1} = \frac{\delta \mathcal{L}}{\delta \hat{m}_1} \frac{\partial \hat{m}_1}{\partial \tilde{m}_1} \frac{\partial \tilde{m}_1}{\partial \tilde{m}_1}, \\ \\ \frac{\delta \mathcal{L}}{\delta \tilde{m}_2} = \frac{\delta \mathcal{L}}{\delta \hat{m}_2} \frac{\partial \hat{m}_2}{\partial \tilde{m}_2} \frac{\partial \tilde{m}_2}{\partial \tilde{m}_2}, \\ \\ \frac{\delta \mathcal{L}}{\delta \hat{\theta}} = \frac{\delta \mathcal{L}}{\delta \hat{\theta}} \frac{\partial \hat{\theta}}{\partial \tilde{\theta}} \frac{\partial \tilde{\theta}}{\partial \hat{\theta}}, \end{array} \right. \qquad (A.5) $$

where the derivative of the projected variable w.r.t. the $(\tilde{m}_1, \tilde{m}_2)$ is expressed as follows:

$$ \frac{\partial \hat{m}_i}{\partial \tilde{m}_i} := \frac{\left(1 - \frac{\tanh(\beta(1-\eta)) + \tanh(\beta(\tilde{m}_i - 1 + \eta))}{\tanh(\beta\eta) + \tanh(\beta(1-\eta))}\right) (\tanh(\beta\eta) + \tanh(\beta(\tilde{m}_i - \eta)))}{\tanh(\beta\eta) + \tanh(\beta(1-\eta))} \\ - \frac{\tilde{m}_i \beta (1 - \tanh(\beta(\tilde{m}_i - 1 + \eta))^2) (\tanh(\beta\eta) + \tanh(\beta(\tilde{m}_i - \eta)))}{(\tanh(\beta\eta) + \tanh(\beta(1-\eta)))^2} \\ + \frac{\tilde{m}_i \left(1 - \frac{\tanh(\beta(1-\eta)) + \tanh(\beta(\tilde{m}_i - 1 + \eta))}{\tanh(\beta\eta) + \tanh(\beta(1-\eta))}\right) \beta (1 - \tanh(\beta(\tilde{m}_i - \eta))^2)}{\tanh(\beta\eta) + \tanh(\beta(1-\eta))} \\ + \frac{\tanh(\beta(1-\eta)) + \tanh(\beta(\tilde{m}_i - 1 + \eta))}{\beta(\tanh(\beta\eta) + \tanh(\beta(1-\eta)))} \\ + \frac{\left(\frac{\beta-1}{\beta} + \frac{\tilde{m}_i}{\beta}\right) \beta (1 - \tanh(\beta(\tilde{m}_i - 1 + \eta))^2)}{\tanh(\beta\eta) + \tanh(\beta(1-\eta))}. \qquad (A.6) $$

Therefore, $\delta_{\hat{d}_i} \mathcal{L}$ reads:

$$ \delta_{\hat{d}_i} \mathcal{L} = \frac{\partial J}{\partial \hat{d}_i} + \frac{\partial \mathcal{X}^*}{\partial \hat{d}_i} \nabla p \cdot \nabla p_a - \epsilon^2 \frac{\partial \mathcal{X}^*}{\partial \hat{d}_i} \nabla p \cdot \nabla T T_a + \epsilon^2 \frac{\partial \mathcal{K}^*}{\partial \hat{d}_i} \nabla T \cdot \nabla T_a + \frac{\partial Q(\hat{m}_1, \hat{m}_2)}{\partial \hat{d}_i} \quad \text{in D,} \qquad (A.7) $$

where $\frac{\partial \mathcal{X}^*}{\partial \hat{d}_i}$ can be expanded as follows:

$$ \frac{\partial \mathcal{X}^*}{\partial \hat{m}_1} := \left[ \frac{\partial \mathcal{X}_{11}}{\partial \hat{m}_1} \cos^2(\hat{\theta}) + \frac{\partial \mathcal{X}_{22}}{\partial \hat{m}_1} \sin^2(\hat{\theta}) \left( \frac{\partial \mathcal{X}_{11}}{\partial \hat{m}_1} - \frac{\partial \mathcal{X}_{22}}{\partial \hat{m}_1} \right) \sin(\hat{\theta}) \cos(\hat{\theta})}{\left( \frac{\partial \mathcal{X}_{11}}{\partial \hat{m}_1} - \frac{\partial \mathcal{X}_{22}}{\partial \hat{m}_1} \right) \sin(\hat{\theta}) \cos(\hat{\theta})} + \frac{\partial \mathcal{X}_{11}}{\partial \hat{m}_1} \sin^2(\hat{\theta}) \right], \qquad (A.8a) $$

$$ \frac{\partial \mathcal{X}^*}{\partial \hat{m}_2} := \left[ \frac{\partial \mathcal{X}_{11}}{\partial \hat{m}_2} \cos^2(\hat{\theta}) + \frac{\partial \mathcal{X}_{22}}{\partial \hat{m}_2} \sin^2(\hat{\theta}) \left( \frac{\partial \mathcal{X}_{11}}{\partial \hat{m}_2} - \frac{\partial \mathcal{X}_{22}}{\partial \hat{m}_2} \right) \sin(\hat{\theta}) \cos(\hat{\theta})}{\left( \frac{\partial \mathcal{X}_{11}}{\partial \hat{m}_2} - \frac{\partial \mathcal{X}_{22}}{\partial \hat{m}_2} \right) \sin(\hat{\theta}) \cos(\hat{\theta})} + \frac{\partial \mathcal{X}_{11}}{\partial \hat{m}_2} \sin^2(\hat{\theta}) \right], \qquad (A.8b) $$
# Declarations

## Conflict of interest

The authors declare that they have no competing financial interests for this paper.

## Replication of results

A detailed implementation and algorithm of the proposed method have been presented in Section 3 and one can follow them and reproduce the results. In case of further queries, please contact the corresponding author(s).

## Author contributions

Hao Li (hli@sdu.dk): Conceptualisation, Methodology, Software, Validation, Writing – original draft. Pierre Jolivet (pierre@joliv.et): Software, Validation, Supervision, Writing – Review & Editing. Joe Alexandersen (joal@sdu.dk): Conceptualisation, Methodology, Supervision, Writing – Review & Editing.

# References

Ahmed Mohammed Adham, Normah Mohd-Ghazali and Robiah Ahmad. Thermal and hydrodynamic analysis of microchannel heat sinks: A review. *Renewable and Sustainable Energy Reviews*, 21:614–622, 2013.

Joe Alexandersen and Casper Schousboe Andreasen. A review of topology optimisation for fluid-based problems. *Fluids*, 5(1):29, 2020.

Grégoire Allaire and Zakaria Habibi. Homogenization of a conductive, convective, and radiative heat transfer problem in a heterogeneous domain. *SIAM Journal on Mathematical Analysis*, 45(3):1136–1178, 2013.

Grégoire Allaire, François Jouve and Anca-Maria Toader. Structural optimization using sensitivity analysis and a level-set method. *Journal of Computational Physics*, 194(1):363–393, 2004.

Grégoire Allaire, Lorenzo Cavallina, Nobuhito Miyake, Tomoyuki Oka and Toshiaki Yachimura. The homogenization method for topology optimization of structures: old and new. *Interdisciplinary Information Sciences*, 25(2):75–146, 2019.

Grégoire Allaire, Perle Geoffroy-Donders and Olivier Pantz. Topology optimization of modulated and oriented periodic microstructures by the homogenization method. *Computers & Mathematics with Applications*, 78(7):2197–2229, 2019.

Grégoire Allaire. Homogenization of the navier-stokes equations in open sets perforated with tiny holes i. abstract framework, a volume distribution of holes. *Archive for Rational Mechanics and Analysis*, 113:209–259, 1991.

Grégoire Allaire. Homogenization of the navier-stokes equations in open sets perforated with tiny holes ii: Non-critical sizes of the holes for a volume distribution and a surface distribution of holes. *Archive for rational mechanics and analysis*, 113:261–298, 1991.

Erik Andreassen and Casper Schousboe Andreasen. How to determine composite material properties using numerical homogenization. *Computational Materials Science*, 83:488–495, 2014.

Jean-Louis Auriault, Claude Boutin and Christian Geindreau. *Homogenization of coupled phenomena in heterogenous media*, volume 149. John Wiley & Sons, 2010.
$$ \frac{\partial \mathbf{X}^*}{\partial \hat{\theta}} := \begin{bmatrix} \sin(\hat{\theta}) \cos(\hat{\theta}) (2\mathcal{X}_{22} - 2\mathcal{X}_{11}) & (\mathcal{X}_{22} - \mathcal{X}_{11}) (\sin^2(\hat{\theta}) - \cos^2(\hat{\theta})) \\ (\mathcal{X}_{22} - \mathcal{X}_{11}) (\sin^2(\hat{\theta}) - \cos^2(\hat{\theta})) & \sin(\hat{\theta}) \cos(\hat{\theta}) (2\mathcal{X}_{11} - 2\mathcal{X}_{22}) \end{bmatrix}. \quad (A.8c) $$

Likewise, $\frac{\partial \mathcal{K}^*}{\partial \hat{d}_i}$ can be expanded using the same formulation as described above, and will be elaborated on here.

The weak form of the adjoint equations derived above are given here. As for the objective function, the adjoint equations are as follows:
Find $T_a \in \mathcal{P}$ with $T_a = 0$ on $\Gamma_{in}$ such that:

$$ \begin{aligned} & \int_D - (\epsilon^2 \mathbf{X}^* \nabla p \cdot \nabla \delta T) T_a \, d\Omega + \int_D \mathcal{K}^* \nabla \delta T \cdot \nabla T_a \, d\Omega \\ & = - \int_D \hat{m}_1, \hat{m}_2 \delta T \, d\Omega \quad \forall \delta T \in \mathcal{P}, T_a \in \mathcal{P}, \end{aligned} \quad (A.9a) $$

Find $p_a \in \mathcal{P}$ with $p_a = 0$ on $\Gamma_{out}$ such that:

$$ \int_D \mathbf{X}^* \nabla p_a \cdot \nabla \delta p \, d\Omega = \int_D (\epsilon^2 \mathbf{X}^* \nabla \delta p \cdot \nabla T) T_a \, d\Omega \quad \forall \delta p \in \mathcal{P}, p_a \in \mathcal{P}. \quad (A.9b) $$

As for the pressure drop constraint, $T_a = 0$ in $D$. Thus, the adjoint system reads:
Find $p_a \in \mathcal{P}$ with $p_a = 0$ on $\Gamma_{out}$ such that:

$$ \int_D \mathbf{X}^* \nabla p_a \cdot \nabla \delta p \, d\Omega = -\epsilon^2 \frac{1/\mathrm{DP}_0}{|\Gamma_{in}|} \int_{\Gamma_{in}} \delta p \, d\Gamma \quad \forall \delta p \in \mathcal{P}, p_a \in \mathcal{P}. \quad (A.10) $$
Satish Balay, Shrirang Abhyankar, Mark F. Adams, Steven Benson, Jed Brown, Peter Brune, Kris Buschelman, Emil Constantinescu, Lisandro Dalcin, Alp Dener, Victor Eijkhout, Jacob Faibussowitsch, William D. Gropp, Václav Hapla, Tobin Isaac, Pierre Jolivet, Dmitry Karpeev, Dinesh Kaushik, Matthew G. Knepley, Fande Kong, Scott Kruger, Dave A. May, Lois Curfman McInnes, Richard Tran Mills, Lawrence Mitchell, Todd Munson, Jose E. Roman, Karl Rupp, Patrick Sanan, Jason Sarich, Barry F. Smith, Stefano Zampini, Hong Zhang, Hong Zhang, and Junchao Zhang. PETSc/TAO users manual. Technical Report ANL-21/39 - Revision 3.19, Argonne National Laboratory, 2023.

Martin Philip Bendsøe and Noboru Kikuchi. Generating optimal topologies in structural design using a homogenization method. *Computer Methods in Applied Mechanics and Engineering*, 71(2):197-224, 1988.

Martin P Bendsøe and Ole Sigmund. Material interpolation schemes in topology optimization. *Archive of Applied Mechanics*, 69(9-10):635-654, 1999.

Thomas Borrvall and Joakim Petersson. Topology optimization of fluids in stokes flow. *International journal for numerical methods in fluids*, 41(1):77-107, 2003.

Lianxiong Chen, Yifei Fan, Ran Zhang, Xihua Chu and Hui Liu. Scale-separated vcut level set method for designing the graded stochastic voronoi cellular structure. *Structural and Multidisciplinary Optimization*, 66(3):54, 2023.

Lianxiong Chen, Yu Pan, Xihua Chu, Hui Liu and Xinzhong Wang. Multiscale design and experimental verification of voronoi graded stochastic lattice structures for the natural frequency maximization problem. *Acta Mechanica Sinica*, 39(8):422445, 2023.

Ercan M Dede, Yuqing Zhou and Tsuyoshi Nomura. Inverse design of microchannel fluid flow networks using turing pattern dehomogenization. *Structural and Multidisciplinary Optimization*, 62(4):2203-2210, 2020.

Ercan M Dede, Yuqing Zhou, Tomoki Tambo, Feng Zhou, Danny J Lohan and Tsuyoshi Nomura. Measurement of low reynolds number flow emanating from a turing pattern microchannel array using a modified bernoulli equation technique. *Experimental Thermal and Fluid Science*, 139:110722, 2022.

Martin Ohrt Elingaard, Niels Aage, Jakob Andreas Bærentzen and Ole Sigmund. De-homogenization using convolutional neural networks. *Computer Methods in Applied Mechanics and Engineering*, 388:114197, 2022.

Florian Feppon. High order homogenization of the stokes system in a periodic porous medium. *SIAM Journal on Mathematical Analysis*, 53(3):2890-2924, 2021.

Florian Feppon. Asymptotic expansions of Stokes flows in finite periodic channels. working paper or preprint, February 2024.

Florian Feppon. Multiscale topology optimization of modulated fluid microchannels based on asymptotic homogenization. *Computer Methods in Applied Mechanics and Engineering*, 419:116646, 2024.

Da Geng, Chuang Wei, Yichang Liu and Mingdong Zhou. Concurrent topology optimization of multi-scale cooling channels with inlets and outlets. *Structural and Multidisciplinary Optimization*, 65(8):234, 2022.

Perle Geoffroy-Donders, Grégoire Allaire and Olivier Pantz. 3-d topology optimization of modulated and oriented periodic microstructures by the homogenization method. *Journal of Computational Physics*, 401:108994, 2020.

Jeroen P Groen and Ole Sigmund. Homogenization-based topology optimization for high-resolution manufacturable microstructures. *International Journal for Numerical Methods in Engineering*, 113(8):1148-1163, 2018.
Jeroen P Groen, Jun Wu and Ole Sigmund. Homogenization-based stiffness optimization and projection of 2d coated structures with orthotropic infill. *Computer Methods in Applied Mechanics and Engineering*, 349:722-742, 2019.

Jeroen P Groen, Florian C Stutz, Niels Aage, Jakob A Bærentzen and Ole Sigmund. De-homogenization of optimal multi-scale 3d topologies. *Computer Methods in Applied Mechanics and Engineering*, 364:112979, 2020.

James K Guest and Jean H Prévost. Design of maximum permeability material structures. *Computer Methods in Applied Mechanics and Engineering*, 196(4-6):1006-1017, 2007.

Sarah N Hankins, Yuqing Zhou, Danny J Lohan and Ercan M Dede. Generative design of large-scale fluid flow structures via steady-state diffusion-based dehomogenization. *Scientific Reports*, 13(1):14344, 2023.

Frédéric Hecht. New development in FreeFem++. *Journal of Numerical Mathematics*, 20(3-4):251-266, 2012.

Peter Dørffler Ladegaard Jensen, Ole Sigmund and Jeroen P Groen. De-homogenization of optimal 2d topologies for multiple loading cases. *Computer Methods in Applied Mechanics and Engineering*, 399:115426, 2022.

Pierre Jolivet, Victorita Dolean, Frédéric Hecht, Frédéric Nataf, Christophe Prud'homme and Nicole Spillane. High-Performance Domain Decomposition Methods on Massively Parallel Architectures with FreeFem++. *Journal of Numerical Mathematics*, 20(4):287-302, 2012.

Atsushi Kawamoto, Tadayoshi Matsumori, Shintaro Yamasaki, Tsuyoshi Nomura, Tsuguo Kondoh and Shinji Nishiwaki. Heaviside projection based topology optimization by a pde-filtered scalar function. *Structural and Multidisciplinary Optimization*, 44(1):19-24, 2011.

Boyan Stefanov Lazarov and Ole Sigmund. Filters in topology optimization based on helmholtz-type differential equations. *International Journal for Numerical Methods in Engineering*, 86(6):765-781, 2011.

Jaewook Lee, Chiyoung Kwon, Jeonghoon Yoo, Seungjae Min, Tsuyoshi Nomura and Ercan M Dede. Design of spatially-varying orthotropic infill structures using multiscale topology optimization and explicit de-homogenization. *Additive Manufacturing*, 40:101920, 2021.

Doksoo Lee, Wei Chen, Liwei Wang, Yu-Chin Chan and Wei Chen. Data-driven design for metamaterials and multiscale systems: A review. *Advanced Materials*, 36(8):2305254, 2024.

Thérese Lévy. Fluid flow through an array of fixed particles. *International Journal of Engineering Science*, 21(1):11-23, 1983.

Hao Li, Minghao Yu, Pierre Jolivet, Joe Alexandersen, Tsuguo Kondoh, Tiannan Hu, Kozo Furuta, Kazuhiro Izui and Shinji Nishiwaki. Reaction-diffusion equation driven topology optimization of high-resolution and feature-rich structures using unstructured meshes. *Advances in Engineering Software*, 180:103457, 2023.

Hui Liu, Lianxiong Chen, Yi Jiang, Dezhou Zhu, Yufa Zhou and Xinzhong Wang. Multiscale optimization of additively manufactured graded non-stochastic and stochastic lattice structures. *Composite Structures*, 305:116546, 2023.

HA Mohammed, G Bhaskaran, NH Shuaib and Rahman Saidur. Heat transfer and fluid flow characteristics in microchannels heat exchanger using nanofluids: a review. *Renewable and Sustainable Energy Reviews*, 15(3):1502-1512, 2011.

Serdar Ozguc, Liang Pan and Justin A Weibel. Topology optimization of microchannel heat sinks using a homogenization approach. *International Journal of Engineering Science and Technology*, 2019:100001, 2019.
of Heat and Mass Transfer, 169:120896, 2021.

Rahul Kumar Padhy, Aaditya Chandrasekhar and Krishnan Suresh. Fluto: Graded multi-scale topology optimization of large contact area fluid-flow devices using neural networks. *Engineering with Computers*, pages 1–17, 2023.

Rahul Kumar Padhy, Krishnan Suresh and Aaditya Chandrasekhar. Tomas: Topology optimization of multiscale fluid devices using variational autoencoders and super-shapes. *arXiv preprint arXiv:2309.08435*, 2023.

Olivier Pantz and Karim Trabelsi. A post-treatment of the homogenization method for shape optimization. *SIAM Journal on Control and Optimization*, 47(3):1380–1398, 2008.

Nicolò Pollini, Ole Sigmund, Casper Schousboe Andreasen and Joe Alexandersen. A “poor man’s” approach for high-resolution three-dimensional topology design for natural convection problems. *Advances in Engineering Software*, 140:102736, 2020.

George IN Rozvany, Ming Zhou and Torben Birker. Generalized shape optimization without homogenization. *Structural Optimization*, 4(3-4):250–252, 1992.

Enrique Sanchez-Palencia. Fluid flow in porous media. *Non-homogeneous media and vibration theory*, pages 129–157, 1980.

Enrique Sánchez-Palencia. Non-homogeneous media and vibration theory. *Lecture Note in Physics*, Springer-Verlag, 320:57–65, 1980.

E Sanchez-Palencia. On the asymptotics of the fluid flow past an array of fixed obstacles. *International Journal of Engineering Science*, 20(12):1291–1301, 1982.

Krister Svanberg. The method of moving asymptotes—a new method for structural optimization. *International journal for numerical methods in engineering*, 24(2):359–373, 1987.

Akihiro Takezawa, Xiaopeng Zhang, Masaki Kato and Mitsuru Kitamura. Method to optimize an additively-manufactured functionally-graded lattice structure for effective liquid cooling. *Additive Manufacturing*, 28:285–298, 2019.

Akihiro Takezawa, Xiaopeng Zhang and Mitsuru Kitamura. Optimization of an additively manufactured functionally graded lattice structure with liquid cooling considering structural performances. *International Journal of Heat and Mass Transfer*, 143:118564, 2019.

Michael Yu Wang, Xiaoming Wang and Dongming Guo. A level set method for structural topology optimization. *Computer Methods in Applied Mechanics and Engineering*, 192(1-2):227–246, 2003.

Liwei Wang, Zhao Liu, Daicong Da, Yu-Chin Chan, Wei Chen and Ping Zhu. Generalized de-homogenization via sawtooth-function-based mapping and its demonstration on data-driven frequency response optimization. *Computer Methods in Applied Mechanics and Engineering*, 395:114967, 2022.

Junpeng Wang, Rüdiger Westermann and Jun Wu. A streamline-guided dehomogenization approach for structural design. *Journal of Mechanical Design*, 145(2), 2023.

Rebekka V Woldseth, J Andreas Bærentzen and Ole Sigmund. Phasor noise for dehomogenisation in 2d multiscale topology optimisation. *Computer Methods in Applied Mechanics and Engineering*, 418:116551, 2024.

Jun Wu, Niels Aage, Rüdiger Westermann and Ole Sigmund. Infill optimization for additive manufacturing—approaching bone-like porous structures. *IEEE Transactions on Visualization and Computer Graphics*, 24(2):1127–1140, 2017.
Jun Wu, Ole Sigmund and Jeroen P Groen. Topology optimization of multi-scale structures: a review. *Structural and Multidisciplinary Optimization*, 63(3):1455–1480, 2021.

Tong Wu. *Topology optimization of multiscale structures coupling fluid, thermal and mechanical analysis*. PhD thesis, Purdue University, 2019.

Xi Zhao, Mingdong Zhou, Ole Sigmund and Casper Schousboe Andreasen. A “poor man’s approach” to topology optimization of cooling channels based on a darcy flow model. *International Journal of Heat and Mass Transfer*, 116:1108–1123, 2018.

Ming Zhou and George IN Rozvany. The coc algorithm, part ii: Topological, geometrical and generalized shape optimization. *Computer Methods in Applied Mechanics and Engineering*, 89(1-3):309–336, 1991.

Yuqing Zhou, Danny J Lohan, Feng Zhou, Tsuyoshi Nomura and Ercan M Dede. Inverse design of microreactor flow fields through anisotropic porous media optimization and dehomogenization. *Chemical Engineering Journal*, 435:134587, 2022.

# A Sensitivity analysis

## A.1 Continuous adjoint method

The sensitivity analysis is conducted to derive the sensitivity of objective function $J$, and pressure drop constraint function $G_1$ w.r.t. the projected design variable $\hat{d} = (\hat{m}_1, \hat{m}_2, \hat{\theta})$. For conciseness, throughout the remainder of this appendix, we utilise $\mathbf{X}^*$ and $\mathbf{K}^*$ for brevity.

Now, the sensitivity analysis is performed. First, the Lagrange function $\mathcal{L}$ is constructed as follows:

$$ \mathcal{L} := J + \langle p_a, -\operatorname{div}(\mathbf{X}^*\nabla p) \rangle + \langle T_a, -\epsilon^2 \mathbf{X}^*\nabla p \cdot \nabla T - \operatorname{div}(\mathbf{K}^*\nabla T) - Q \rangle, \quad (A.1) $$

where $p_a$ and $T_a$ are the adjoint pressure and adjoint temperature, respectively. According to the Karush-Kuhn-Tucker (KKT) conditions for PDE constrained optimisation problems,

$$ \left\{ \begin{array}{l} \delta_T \mathcal{L} = 0, \\ \delta_p \mathcal{L} = 0. \end{array} \right. \qquad (A.2) $$

The variation of the Lagrange function w.r.t. the state variables $T$ and $p$ are derived as follows:

$$ \delta_T \mathcal{L} = \delta_T J + \int_D -\epsilon^2 \mathbf{X}^* \nabla p \cdot \nabla \delta T T_a \, d\Omega + \int_D \mathbf{K}^* \nabla \delta T \cdot \nabla T_a \, d\Omega = 0 \quad \forall \delta T \in \mathcal{P}, T_a \in \mathcal{P}, \qquad (A.3a) $$

$$ \delta_p \mathcal{L} = \delta_p J + \int_D \mathbf{X}^* \nabla p_a \cdot \nabla \delta p \, d\Omega - \int_D \epsilon^2 \mathbf{X}^* \nabla \delta p \cdot \nabla T T_a = 0 \quad \forall \delta p \in \mathcal{P}, p_a \in \mathcal{P}. \qquad (A.3b) $$
As for the optimum design problem given in Eq. (22), the variation of the cost functions w.r.t. the state variables are derived as follows:

$$ \left\{ \begin{aligned} \delta_T J &= \int_D \hat{m}_1 \hat{m}_2 \delta T \, d\Omega, \\ \delta_p J &= 0, \end{aligned} \right. \qquad (A.4a) $$

$$ \left\{ \begin{aligned} \delta_T G_1 &= 0, \\ \delta_p G_1 &= \epsilon^2 \frac{1}{|\Gamma_{in}|} \int_{\Gamma_{in}} \delta p \, d\Gamma. \end{aligned} \right. \qquad (A.4b) $$

Substituting Eqs. (A.3) and (A.4) into Eq. (A.2), the adjoint system can be obtained, see the corresponding weak forms Eqs. (A.9) and (A.10).

Finally, the functional derivative of the Lagrange function w.r.t. design variables, $\delta_d \mathcal{L}$, can be derived using the chain rule as follows:

$$ \left\{ \begin{aligned} \frac{\delta \mathcal{L}}{\delta m_1} &= \frac{\delta \mathcal{L}}{\delta \hat{m}_1} \frac{\partial \hat{m}_1}{\partial \tilde{m}_1} \frac{\delta \tilde{m}_1}{\delta m_1}, \\ \frac{\delta \mathcal{L}}{\delta m_2} &= \frac{\delta \mathcal{L}}{\delta \hat{m}_2} \frac{\partial \hat{m}_2}{\partial \tilde{m}_2} \frac{\delta \tilde{m}_2}{\delta m_2}, \\ \frac{\delta \mathcal{L}}{\delta \theta} &= \frac{\delta \mathcal{L}}{\delta \hat{\theta}} \frac{\partial \hat{\theta}}{\partial \tilde{\theta}} \frac{\delta \tilde{\theta}}{\delta \theta}, \end{aligned} \right. \qquad (A.5) $$

where the derivative of the projected variable w.r.t. the $(\tilde{m}_1, \tilde{m}_2)$ is expressed as follows:

$$ \frac{\partial \hat{m}_i}{\partial \tilde{m}_i} := \frac{\left(1 - \frac{\tanh(\beta(1-\eta)) + \tanh(\beta(\tilde{m}_i - 1 + \eta))}{\tanh(\beta\eta) + \tanh(\beta(1-\eta))}\right) (\tanh(\beta\eta) + \tanh(\beta(\tilde{m}_i - \eta)))}{\tanh(\beta\eta) + \tanh(\beta(1-\eta))} \\ - \frac{\tilde{m}_i \beta (1 - \tanh(\beta(\tilde{m}_i - 1 + \eta))^2) (\tanh(\beta\eta) + \tanh(\beta(\tilde{m}_i - \eta)))}{(\tanh(\beta\eta) + \tanh(\beta(1-\eta)))^2} \\ + \frac{\tilde{m}_i \left(1 - \frac{\tanh(\beta(1-\eta)) + \tanh(\beta(\tilde{m}_i - 1 + \eta))}{\tanh(\beta\eta) + \tanh(\beta(1-\eta))}\right) \beta (1 - \tanh(\beta(\tilde{m}_i - \eta))^2)}{\tanh(\beta\eta) + \tanh(\beta(1-\eta))} \\ + \frac{\tanh(\beta(1-\eta)) + \tanh(\beta(\tilde{m}_i - 1 + \eta))}{\beta(\tanh(\beta\eta) + \tanh(\beta(1-\eta)))} \\ + \frac{\left(\frac{\beta-1}{\beta} + \frac{\tilde{m}_i}{\beta}\right) \beta (1 - \tanh(\beta(\tilde{m}_i - 1 + \eta))^2)}{\tanh(\beta\eta) + \tanh(\beta(1-\eta))}. \qquad (A.6) $$

Therefore, $\delta_{\hat{d}_i} \mathcal{L}$ reads:

$$ \delta_{\hat{d}_i} \mathcal{L} = \frac{\partial J}{\partial \hat{d}_i} + \frac{\partial \boldsymbol{\chi}^*}{\partial \hat{d}_i} \nabla p \cdot \nabla p_a - \epsilon^2 \frac{\partial \boldsymbol{\chi}^*}{\partial \hat{d}_i} \nabla p \cdot \nabla T T_a + \epsilon^2 \frac{\partial \boldsymbol{\mathcal{K}}^*}{\partial \hat{d}_i} \nabla T \cdot \nabla T_a + \frac{\partial Q(\hat{m}_1, \hat{m}_2)}{\partial \hat{d}_i} \quad \text{in D,} \qquad (A.7) $$