# Multi-Objective Topology Optimization of Conjugate Heat Transfer Using Level Sets and Anisotropic Mesh Adaptation

Philippe Meliga, Wassim Abdel Nour, Delphine Laboureur, Damien Serret, Elie Hachem

► To cite this version:

Philippe Meliga, Wassim Abdel Nour, Delphine Laboureur, Damien Serret, Elie Hachem. Multi-Objective Topology Optimization of Conjugate Heat Transfer Using Level Sets and Anisotropic Mesh Adaptation. Fluids, 2024, 9 (5), pp.105. 10.3390/fluids9050105. hal-04750780

HAL Id: hal-04750780

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

Submitted on 23 Oct 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
Article

# Multi-objective topology optimization of conjugate heat transfer using level-sets and anisotropic mesh adaptation

Philippe Meliga<fnref n="1" />, Wassim Abdel Nour<fnref n="1" />,<fnref n="2" />, Delphine Laboureur<fnref n="3" />, Damien Serret<fnref n="2" /> and Elie Hachem <fnref n="1" />,*

**Abstract:** This study proposes a new computational framework for multi-objective topology optimization of conjugate heat transfer systems using a continuous adjoint approach. It relies on a monolithic solver for the coupled steady-state Navier-Stokes and heat equations, that combines finite elements stabilized by the variational multiscale method, level set representations of the fluid-solid interfaces and immersed modeling of heterogeneous materials (fluid/solid) to ensure that the proper amount of heat is exchanged to the ambient fluid by solid objects in arbitrary geometry. At each optimization iteration, anisotropic mesh adaptation is applied in near-wall regions automatically captured by the level-set, which reduces the computational burden related to the call of the finite element solver, compared to classical optimization schemes working on uniform grids with similar mesh refinement. Since we work here under the constraint of a given number of nodes, the approach improves the accuracy in the geometric description of all layouts, while allowing to transfer said accuracy to the state and adjoint simulation models. Both the resolution and remeshing steps are performed in a massively parallel framework allowing for the optimization of large scale systems. The developed solver is validated first by minimizing dissipation in a flow splitter device, for which the method delivers relevant optimal designs over a wide range of volume constraints and flow rate distributions over the multiple outlet orifices, while providing superior accuracy over reference data computed on isotropic meshes (in the sense that the layouts are more smooth, and the solutions are better resolved). The proposed scheme is then applied to a two-dimensional heat transfer problem, using bi-objective cost functionals combining power dissipation and thermal recoverable power. A comprehensive parametric study reveals a complex arrangement of optimal solutions on the Pareto front, with multiple branches of symmetric and asymmetric designs, some of them previously unreported. Finally, the algorithmic developments are substantiated with several three-dimensional numerical examples tackled under fixed weights for heat transfer and power dissipation, for which we show that the optimal layouts computed at low Reynolds number, that are intrinsically relevant to a broad range of microfluidic application, can also serve as smooth solutions to high-Reynolds-number engineering problems of practical interest.

**Keywords:** Topology Optimization; Fluid mechanics; Conjugate heat transfer; Level Set Method; Anisotropic mesh adaptation; Pareto front; Thermal control

## 1. Introduction

Topology optimization is the mathematical science of optimal material allocation in a volume under predefined objectives and constraints. Such an approach originates from solid mechanics [1,2], where it has matured into a powerful, reliable and increasingly available tool for engineers in the early stages of complex structural design processes [3,4], while spreading to a variety of multiphysics applications modeled after partial differential equations; see Refs. [5,6] for surveys of the evolving methods and applications. The mathematical

**Citation:** Meliga, P.; Adel Nour, W.; Laboureur, D., Serret, D.; Meliga, P.; Hachem, E. Multi-objective topology optimization of conjugate heat transfer using level-sets and anisotropic mesh adaptation.
*Fluids* **2023**, *1*, 0. https://doi.org/

Received:

Revised:

Accepted:

Published:

**Copyright:** ©2024 by the authors.
Submitted to *Fluids* for possible
open access publication under
the terms and conditions of
the Creative Commons Attri-
bution (CC BY) license (https://
creativecommons.org/licenses/by/
4.0/).

<footnote n="1">¹ Mines Paris, PSL University, Centre for material forming (CEMEF), UMR CNRS, 06904 Sophia Antipolis, France</footnote>

<footnote n="2">² TEMISTh SAS, Technocentre des Florides, 13700 Marignane, France</footnote>

<footnote n="3">³ The von Karman Institute for Fluid Dynamics, Waterloosesteenweg 72, B-1640 Sint-Genesius Rode, Belgium</footnote>

<footnote>* Correspondence: elie.hachem@minesparis.psl.eu; Tel.: +33 4 93 95 74 58</footnote>
foundation of topology optimization is rooted in iterative analysis and design update steps, often steered by gradient evaluations. Its key advantage (compared the size and shape optimization methods it has emerged from) is that it usually perform beyond expectations as the final design has unconstrained complexity and can possibly meet conflicting requirements and complex interdependencies between design parameters and system response.

In the context of fluid flow problems, topology optimisation becomes a question of what is the best path for a fluid to flow in a prescribed design domain, or equivalently where to enforce relevant boundary conditions for the flow problem. Leaving aside explicit boundary methods, that represent the fluid-solid interface by edges or faces of a body-fitted mesh, and have limited flexibility to handle complicated topological changes, the prevalent classes of methods for topology optimization are the density and the level set methods. Density methods rely on a Brinkman penalization of the solid domain, where the flow is modeled as a fictitious porous material with very low permeability [1,7,8]. They manage drastic topological changes, as the gradient information (or sensitivity) is distributed over a large part of the domain, but can lead to spurious or leaking flows if the penalization factor is not well-calibrated (since the velocity and pressure fields are computed in both the solid and fluid regions). Level set methods conversely model the solid boundaries by iso-contours of a level set function [9–11]. They lack a nucleation mechanism to create new holes, due to the sensitivities being located only at the solid-fluid interface (which is often relieved using initial designs with many holes), but easily handle complicated topological changes (e.g., merging or cancellation of holes), and allow for well defined, crisp interface representations while avoiding the intermediate material phases (grayscales) and mesh-dependent spatial oscillations of the interface geometry (staircasing) often encountered in density methods [12].

The focus here is on conjugate convective heat transfer systems, in which temperature variations occur within the fluid and solid material due to thermal interactions dominated by conduction in the solid and convection in the solid. This is a matter of great engineering interest, as many industries have embraced the ability of topology optimization to improve the performances and cut the production costs of thermal devices like heat exchangers (to regulate process temperatures and ensure that machinery, chemicals, water, gas, drugs or food remain within safe operating conditions), finned surfaces, microelectronic equipment and heat sinks, and deliver more compact designs with less mass, less frictional losses and better thermal efficiency. In this context, the early related literature can be broadly classified into two categories: pure heat conduction problems maximizing heat evacuation from singular tree-like optimal structures of high conduction material, and pure fluid flow problems minimizing the power dissipated inside the domain (alternatively minimizing drag or maximizing the outlet flow uniformity) from complex channel layouts in the diffusion and convection dominated regimes; see [13,14] for recent reviews and references therein. Since then, the topology optimization of coupled thermal-fluid problems (that combine both aspects, and thus require dual objective function strategies to increase heat transfer while keeping dissipation as low as possible) has become an active field of research. Although variants of the level set method have received attention recently [15–17], the vast majority of available studies implement a density-based monolithic approach [18–26] to overcome the fact that the fluid-solid interface is constantly changing over the optimization process, which makes using either a constant heat transfer coefficient or some specific surrogate model to model the heat transfer between the fluid and its surrounding ineffective. A variety of models have been used, ranging from oversimplified (dismissing the thermal conductivity differences between the solid and fluid regions [22] or numerically imposing a constant solid temperature [21]), to highly realistic (full coupling of flow and heat transfer under dual objective function strategies [18]).

The norm in topology optimization is to employ fixed finite element meshes with uniform (or close-to-uniform) element size, small enough that all relevant physical phenomena are reliably captured, but not so small that the cost of performing the optimization becomes unaffordable. A recent trend has been to use adaptive remeshing techniques to maintain a
competitive computational cost. Such an approach consists in generating a coarse base grid, then in adding recursively finer and finer subgrids in the regions requiring higher resolution, either until a maximum level of refinement is reached, or until the local truncation error drops below a certain tolerance (for more sophisticated implementations endowed with error estimation routines). Within the context of pure fluid flow problems, particular emphasis has been put on (but not limited to) adaptive meshing refinement (AMR) schemes, using both density [27,28] and level set methods [29,30]; see also [31] for an application to phase field methods<fnref n="1" /> and [32,33] for recent efforts applying a different remeshing scheme to a combination of level set functions and adaptive body-conforming meshes for several coupled physics applications, including convective heat transfer.

There is still ample room for progress, though, as almost all adaptive algorithms applied so far to fluid flow topology optimization support only isotropic size maps. Fluid dynamics conversely involves convection dominated phenomena for which anisotropic meshes are highly desirable [34], especially in the vicinity of the solid boundaries, where the fluid velocity exhibits steep gradients in the wall-normal direction and skin-friction plays a defining role. The premise of this study is that the ability to generate highly stretched elements in boundary layer regions can substantially increase the accuracy of the geometric representation (compared to what is often seen in topology optimization of flow problems) and naturally convey said accuracy to the numerical solution without sophisticated interpolation or discretization techniques. We note that this is all perfectly in line with the recommendations made in [14] to improve upon the current state of the art. Nonetheless, our literature review did not reveal any other study combining anisotropic mesh adaptation and topology optimization of thermal-fluid problems, besides the (pure fluid flow) density-based optimisation of Stokes flow in Ref. [28], possibly because of the notorious difficulty of finding spatial discretization schemes that meet the level of robustness required by automatic anisotropic mesh adaptation.

This research intends to fill the gap by introducing a novel numerical framework for topology optimization of conjugate heat transfer systems governed by the coupled Navier-Stokes and heat equations. The latter combines level set methods and anisotropic mesh adaptation to handle arbitrary geometries immersed in an unstructured mesh. The governing equations are solved by a variational multiscale (VMS) stabilized finite element method supporting elements of aspect ratio up to the order of 1000:1 [35]. The same numerical method is used to solve the adjoint Navier-Stokes and heat equations underlying the sensitivity analysis needed to evolve the level set function. The metric map providing both the size and the stretching of mesh elements in a very condensed information data is derived from the level set. A posteriori anisotropic error estimator is then used to minimize the interpolation error under the constraint of a prescribed number of nodes in the mesh. The latter can be adjusted over the course of optimization, meaning that the base grid can be either refined or coarsened on demand (in contrast with AMR, whose total number of mesh elements cannot be controlled, and whose mesh cannot be coarsened further than its base configuration). This is expected to achieve further speed-ups (by reducing the cost of modelling the solid material away from the interface) and also to help improve manufacturability of the optimal design, which remains an issue as most classical topology optimization methods render organic designs that can be difficult to translate into computer-aided design models.

The paper is organized as follows: the topology optimization equations are formulated in Sec. 2. The immersed, stabilized finite element numerical framework and anisotropic mesh adaptation algorithm used to perform the design update step are described in Secs. 3 and 4, respectively. The details of the implemented topology optimization algorithm are provided in Sec. 5. Finally, numerical experiments showcasing the potential of the approach to increase the recoverable thermal power while minimizing the dissipated power in two

<footnote n="1">¹ Another class of interface capturing schemes that remain less popular due to the larger computational cost and the difficulty of numerically discretizing the biharmonic phase-field equation.</footnote>
dimensional (2-D) and three-dimensional (3-D) systems are presented in section 6; with particular attention paid to highlighting the improved accuracy during all stages of the optimization. Finally, the numerical cost is discussed in section 7, where we also debate the generalization to high-Reynolds number regimes.

## 2. Gradient-based model for conjugate heat transfer topology optimization

### 2.1. Governing equations

In the following, we denote by $\Omega$ a fixed, open bounded domain in $\mathbb{R}^d$ (with $d$ the space dimension), with boundary $\partial\Omega$ oriented with inward-pointing normal vector $n$. Throughout this study, $\Omega = \Omega_f \cup \Omega_s$ is the disjoint reunion of a fluid domain $\Omega_f$ and a solid domain $\Omega_s$, separated by an interface $\Gamma = \Omega_f \cap \Omega_s$, whose position we seek to optimize with respect to a certain measure of performance. The boundary $\partial\Omega_f$ of the fluid domain is split into interface, inlet $\Gamma_i$ (defined as the combined boundary of all surfaces where fluid enters the domain), and outlet $\Gamma_o$ (the combined boundary of all surfaces where fluid leaves the domain). The boundary $\partial\Omega_s$ of the solid domain is split into interface, isothermal $\Gamma_{iso}$ (the combined boundary of all surfaces where temperature is prescribed), and adiabatic (the combined boundary of all surfaces where heat is exchanged with no gain or loss).

Mathematically, the problem is characterized by a set of physical variables determined as the solutions of partial differential equations (PDEs), themselves derived from modeling considerations and serving as constraint during the optimization. Here, the level set method is used to localize and capture the interface between the fluid and solid domains from the zero iso-value of a smooth level set function, classically the signed distance function defined as

$$ \varphi(\mathbf{x}) = \begin{cases} -\mathrm{dist}(\mathbf{x}, \Gamma) & \text{if } \mathbf{x} \in \Omega_f, \\ 0 & \text{if } \mathbf{x} \in \Gamma, \\ \mathrm{dist}(\mathbf{x}, \Gamma) & \text{if } \mathbf{x} \in \Omega_s, \end{cases} \qquad (1) $$

with the convention that $\varphi < 0$ in the fluid domain. The immerse volume method (IVM, more details in Sec. 3) is then used to extend to the whole domain $\Omega$ the coupled steady incompressible Navier-Stokes and heat equations governing the flow motion in the fluid domain $\Omega_f$. Simply put, the IVM fills $\Omega_s$ with a fictitious fluid mimicking a solid phase, but avoids introducing discontinuities at the interface. It then solves fluid-like equations in the whole domain $\Omega$, using non-homogeneous material properties adequately interpolated over a small layer around the zero level-set (whose thickness is user-defined and does not increase in size during the optimization, unlike the homogenization method or any other generalized material method) and otherwise equal to their fluid and solid values. Such an approach is simpler that the Ersatz material approach [36], that adds a Brinkman penalization term to the Navier-Stokes equations and has clear connections to density-based methods through the material distribution [30]. This yields the monolithic formulation

$$ \nabla \cdot \mathbf{u} = 0 \qquad \text{in } \Omega, \qquad (2) $$

$$ \rho \mathbf{u} \cdot \nabla \mathbf{u} = - \nabla p + \nabla \cdot (2\mu \boldsymbol{\varepsilon}(\mathbf{u})) \qquad \text{in } \Omega, \qquad (3) $$

$$ \rho c_p \mathbf{u} \cdot \nabla T = \nabla \cdot (k \nabla T) \qquad \text{in } \Omega, \qquad (4) $$

where $\mathbf{u}$ is the velocity, $p$ the pressure, $\boldsymbol{\varepsilon}(\mathbf{u}) = (\nabla \mathbf{u} + \nabla \mathbf{u}^T)/2$ the rate of deformation tensor, and $\rho$ and $\mu$ are phase dependent density, dynamic viscosity, thermal conductivity and specific heat.

We seek here to minimize a cost function $J$ that we assume can be formulated as a surface (rather than volume) integral over all or any part of inlet and/or outlet (as is most often the case in topology optimization), i.e.,

$$ J_s = \int_{\Gamma_i \cup \Gamma_o} J ds. \qquad (5) $$
This is done using the continuous adjoint method to find the cost function sensitivity to variations of a design variable $\beta$ physically representing deformations of the interface under the form of local surface normal displacements. In short, the desired sensitivity to a displacement $\beta$ at some specific point of the interface is given by

$$ \delta_{\beta} J_s = \int_{\Gamma} \beta \mu (\nabla \tilde{\mathbf{u}} \cdot \mathbf{n}) \cdot (\nabla \mathbf{u} \cdot \mathbf{n}) \, \mathrm{d}s + \int_{\Gamma} \beta k (\nabla \tilde{T} \cdot \mathbf{n}) (\nabla T \cdot \mathbf{n}) \, \mathrm{d}s, $$

where $(\tilde{\mathbf{u}}, \tilde{T})$ are the adjoint velocity and temperature solution to the adjoint Navier-Stokes and heat equations written in monolithic formulation as

$$ \nabla \cdot \tilde{\mathbf{u}} = 0 \qquad \text{in } \Omega, \qquad (6) $$

$$ -\rho \mathbf{u} \cdot \nabla \tilde{\mathbf{u}} + \rho \nabla \tilde{\mathbf{u}}^T \cdot \tilde{\mathbf{u}} = \nabla \tilde{p} + \nabla \cdot (2\mu\varepsilon(\tilde{\mathbf{u}})) \qquad \text{in } \Omega, \qquad (7) $$

$$ -c_p \mathbf{u} \cdot \nabla \tilde{T} = \nabla \cdot (k \nabla \tilde{T}) \qquad \text{in } \Omega, \qquad (8) $$

with $\tilde{p}$ the associated adjoint pressure. The main steps of the method are described in [1], together with the boundary conditions appended to both the state and adjoint equations; see also [37–39] for further deepening on the topic. This enables efficient design update schemes via first-order gradient descent methods, for instance the simplest steepest-descent algorithm implemented herein moves down the cost function, in the direction of the steepest slope using

$$ \beta = -\mu(\nabla \tilde{\mathbf{u}} \cdot \mathbf{n}) \cdot (\nabla \mathbf{u} \cdot \mathbf{n}) - k(\nabla \tilde{T} \cdot \mathbf{n})(\nabla T \cdot \mathbf{n}), \qquad (9) $$

up to a positive multiplicative factor to control the step taken in the gradient direction.

## 2.2. Multi-objective optimization

A classical objective in topology optimization of conjugate heat transfer is to maximize heat transfer in the domain without increasing the mechanical pumping power that need be spent to overcome friction and move the fluid through the device (nor blocking the fluid flow). This is done in some studies by maximizing heat transfer under prescribed pressure drop values [22,25], and in others by minimizing pressure drop under prescribed heat transfer performance [40]. As further explained in Sec. 6, we rather use here multi-criteria optimization and minimize the linear weighted sum of a hydraulic cost function $J_{\nu}$ associated to dissipation (to minimize) and a thermal cost function $J_{\Theta}$ measuring heat transfer efficiency (to maximize, the minimization applies to $-J_{\Theta}$). This yields

$$ J = (1 - \omega) J_{\nu} - \omega J_{\Theta}, \qquad (10) $$

where $\omega \in [0; 1]$ is the so-called thermal weigh, a scalar-valued factor weighing the priority given to each objective function ($\omega = 0$ in the pure hydraulic limit, and $\omega = 1$ in the pure thermal limit). In practice, a single point concurrently minimizing both objectives usually does not exist. The solution to such a problem thus aims at identifying the Pareto front [41], i.e., the subset of designs that best manage trade-offs between conflicting criteria, in the sense that further optimizing one cost function decreases the performance of the other one (after which the final design is selected from the Pareto optimal subset by a human decision maker based on subjective preferences).

# 3. Computational methods

A primitive pseudo-code of the procedure for solving the above topology optimization problem is provided in Alg. 1, to repeat until a maximum number of iterations or a convergence threshold has been reached. In a nutshell, this is done here using a finite element immersed numerical framework combining implicit representation of the different domains, level set description of the interface, and anisotropic remeshing capabilities. For the sake of readability, the mesh adaptation algorithm and parallel computational framework, whose implementation in the context of fluid flow topology optimization makes for the main
<caption>**Algorithm 1** Simplified update scheme</caption>

**Require:** Anisotropic mesh adapted to initial interface position

1: loop

2: Compute state

3: Compute adjoint

4: Compute cost function sensitivity

5: Set displacement in the direction of steepest slope

6: Update interface position

7: Generate anisotropic mesh adapted to new interface position

novelty of this study, are presented in the following as stand-alone sections. In the remainder of this section, we walk through each of the other steps and review the various problems involved and the numerical methods for solving them.

### 3.1. Variational multiscale modeling

We solve all equations (i.e., state, adjoint, level set) using equal order linear/linear approximations for both the velocity and pressure variables, a scheme very desirable in the context of the large-scale, three-dimensional applications considered here due to its simplicity of implementation and affordable computing cost. To this end, we solve stabilized weak forms cast in the Variational Multiscale (VMS) framework, that enhance the stability of the Galerkin method via a series of additional integrals over element interior (as a linear/linear discretization otherwise breaks the Babuska-Brezzi condition, which results in spurious node-to-node oscillations). The basic idea is to split all quantities into coarse and fine scale components, corresponding to different levels of resolution, and to approximate the effect of the fine scale (that cannot be resolved by the finite element mesh) onto the coarse scale via consistently derived residual based terms. The various coarse scale variational problems solved by VMS are provided in J, together with additional details about the implicit/explicit discretization in space.

### 3.2. Immersed volume method

The monolithic immersed volume method (IVM) is used to combine the fluid and solid phases of the problem into a single fluid with variable material properties; see [42,43] for details regarding the mathematical formulation in the context of finite element VMS methods. Simply put, this amounts to solving the state equations (2)-(4) and adjoint equations (6)-(8) on a unique mesh of the domain $\Omega$ in which the fluid and solid domains $\Omega_f$ and $\Omega_s$ are immersed, using the level set method to achieve interface tracking between the two different domains. Using the level set function (1) as criterion for anisotropic mesh adaptation (more details provided in Sec. 4) ensures that individual material properties can be distributed accurately and smoothly as possible over the smallest possible thickness around the interface (classically by linear interpolation of the fluid and solid values, using a smooth Heaviside function computed from the level set to avoid discontinuities by creating an interface transition with a thickness of a few elements). The IVM approach is especially relevant to the thermal coupling problems tackled in this research, as having composite conductivity and specific heat means that the amount of heat exchanged at the interface then proceeds solely from the individual material properties on either side of it, and removes the need for a heat transfer coefficient.

In practice, we compute the composite thermal conductivity as the harmonic mean of the solid and fluid values, i.e.,

$$ \frac{1}{k} = \frac{1}{k_f} H_{\epsilon}(\varphi) + \frac{1}{k_s} (1 - H_{\epsilon}(\varphi)), \qquad (11) $$
where $H_\epsilon$ is the smoothed Heaviside function on the fluid domain defined as

$$H_{\epsilon}(\varphi) = \begin{cases} \frac{1}{2} \left( 1 - \frac{\varphi}{\epsilon} - \frac{1}{\pi} \sin\left(\pi \frac{\varphi}{\epsilon}\right) \right) & \text{if } |\varphi| < \epsilon, \\ 0 & \text{if } |\varphi| \ge \epsilon, \end{cases} \quad (12)$$

and $\epsilon$ is a regularization parameter set to $2h_\perp$. This ensures continuity of the heat flux across the interface, as obtained from a steady, no source, one dimensional analysis of the heat flux when the conductivity varies stepwise from one medium to the next (see [44] for detailed derivation and analysis, and [45] for proof of the gain in numerical accuracy with respect to the classical arithmetic mean model). All other properties (density, dynamic viscosity and thermal heat capacity) are constant and equal to the fluid values, and we set the velocity to zero at all grid nodes located inside the solid domain $\Omega_s$. This can be seen as a hard penalty (compared to using a very high solid to fluid viscosity ratio to ensure that the velocity is zero in the solid domain) preventing the fluid from leaking across the immersed interface, that holds numerically because anisotropic mesh adaptation ensures that the interface does not intersect arbitrarily the mesh elements (it precisely aims at aligning the mesh element edges along the interface), which may otherwise compromise the accuracy of the finite element approach.

### 3.3. Interface update scheme

Once the sensitivity analysis has output a displacement $\beta$ in the direction of the steepest slope, the position of the level set is updated solving a transport equation with normal velocity $\beta n/\Delta\tau$, where $\Delta\tau$ is a pseudo-time step to convert from displacement to velocity (of no physical relevance since we are not concerned by the absolute displacement of a given point on the interface, only by its relative displacement with respect to its neighbors). This equation is posed in the whole domain $\Omega$, which is because the normal vector recovered at the interface as $n = \nabla\varphi/\|\nabla\varphi\|$ is easily extended to $\Omega$ using (1). The main problem with this approach is that the level set after transport is generally no longer a distance function, which is especially problematic when a specific remeshing strategy depending on the distance property is used at the interface (as is the case here). As a result, the distance function needs to be reinitialized, which is done here using a coupled convection-reinitialization method wherein the level set function is automatically reinitialized during the resolution of the transport equation. In practice, the signed distance function is cut off using a hyperbolic tangent filter, as defined by

$$\phi = E \tanh\left(\frac{\varphi}{E}\right), \quad (13)$$

with $E$ the cut-off thickness (so the metric property is asymptotically satisfied in the vicinity of the zero iso-value). This filtered level set is then evolved solving the auto-reinitialization equation

$$\partial_{\tau}\phi + \mathbf{a}_{\tau} \cdot \nabla\phi = S, \quad (14)$$

where we note

$$\mathbf{a}_{\tau} = \frac{\beta}{\Delta\tau}\mathbf{n} + \frac{\lambda}{\Delta\tau}\operatorname{sgn}(\phi)\frac{\nabla\phi}{\|\nabla\phi\|}, \quad S = \frac{\lambda}{\Delta\tau}\operatorname{sgn}(\phi)\left(1 - \left(\frac{\phi}{E}\right)^2\right), \quad (15)$$

and $\lambda$ is a parameter homogeneous to a length, set to the mesh size $h_\perp$ in the direction normal to the interface. Such an approach is shown in [46–48] to reduce the computational cost and to ensure a better mass conservation compared to the classical Hamilton–Jacobi method in which both steps are performed in succession). Moreover, since the filtered level set defined in (13) is bounded, Dirichlet boundary conditions $\phi = \pm E$ are easily appended
to Eq. (14) to explicitly design fluid and solid sub-regions of $\partial\Omega$ (in practice, we impose fluid at the inlet and outlet, and solid everywhere else).

# 4. Anisotropic mesh adaptation

## 4.1. Construction of an anisotropic mesh

The main idea of anisotropic, metric-based mesh adaptation is to generate a uniform mesh (with unit length edges and regular elements) in a prescribed Riemannian metric space, but anisotropic and well adapted (with highly stretched elements) in the Euclidean space. Assuming that, in the context of metric-based adaptation methods, controlling the interpolation error suffices to master the global approximation error, the objective can be formulated as finding the mesh, made up of at most $N_n$ nodes, that minimizes the linear interpolation error in the $L^1$ norm. Following the lines of [49,50], an edge-based error estimator combined to a gradient recovery procedure is used to compute, for each node, a metric tensor that prescribes a set of anisotropic directions and stretching factors along these directions, without any direct information from the elements, nor any underlying interpolation. The optimal stretching factor field is obtained by solving an optimization problem using the equi-distribution principle under the constraint of a fixed number of nodes in the mesh, after which a new mesh is generated using the procedure described in [51], based on a topological representation of the computational domain.

## 4.2. Edge error estimate

Given a mesh $\Omega_h$ of the domain $\Omega$, we denote by $\mathbf{x}^{ij}$ the edge connecting a given node $\mathbf{x}^i$ to $\mathbf{x}^j \in \Sigma(i)$, where $\Sigma(i)$ is the set of nodes connected to $\mathbf{x}^i$, and the number of such nodes is noted as $|\Sigma(i)|$. Also, given a regular analytical (scalar) function $\psi$ defined on $\Omega$, and its P1 finite element approximation $\psi_h$ computed on $\Omega_h$, we follow [49] and estimate the interpolation error along the edge $\mathbf{x}^{ij}$ as the projection along the edge of the second derivative of $\psi$. This is obtained projecting along the edge a Taylor expansion of the gradient of $\psi$ at $\mathbf{x}^j$ to give

$$ \varepsilon_{ij} = |\mathbf{g}^{ij} \cdot \mathbf{x}^{ij}|, \qquad (16) $$

where the $i$ and $j$ superscripts indicate nodal values at nodes $\mathbf{x}^i$ and $\mathbf{x}^j$, respectively, $\mathbf{g}^i = \nabla\psi(\mathbf{x}^i)$ is the exact value of the gradient at $\mathbf{x}^i$, and $\mathbf{g}^{ij} = \mathbf{g}^j - \mathbf{g}^i$ is the variation of the gradient along the edge. Although Eq. (16) involves only values of the gradient at the edge extremities and can thus be evaluated without resorting to ressource expensive Hessian reconstruction methods, this however requires the gradient of $\psi$ to be known and continuous at the nodes, which in turn requires full knowledge of $\psi$. Meanwhile, only the linear interpolate $\psi_h$ is known in practice, whose gradient is piecewise constant and discontinuous from element to element (although its projection along the edges is continuous since it depends only on the nodal values of the field).

A recovery procedure is thus used to build a continuous gradient estimator defined directly at the nodes. It is shown in [49] that a suitable error estimate preserving second-order accuracy is obtained substituting the reconstructed gradient for the exact gradient in (16), to give

$$ \varepsilon_{ij} = |\bar{\mathbf{g}}^{ij} \cdot \mathbf{x}^{ij}|, \qquad (17) $$

where $\bar{\mathbf{g}}^{ij} = \bar{\mathbf{g}}^j - \bar{\mathbf{g}}^i$ and we denote by $\bar{\mathbf{g}}^i$ the recovered gradient of $\psi_h$ at node $\mathbf{x}^i$. The latter is defined in a least-square sense as

$$ \bar{\mathbf{g}}^i = \underset{\mathbf{g} \in \mathbb{R}^d}{\operatorname{argmin}} \sum_{j \in \Sigma(i)} |(\mathbf{g} - \nabla\psi_h) \cdot \mathbf{x}^{ij}|^2, \qquad (18) $$
for which an approximate solution using the nodal values as sole input is shown in [49] to be

$$ \bar{\mathbf{g}}^i = (\mathbf{X}^i)^{-1} \cdot \sum_{j \in \Sigma(i)} (\psi_h(\mathbf{x}^j) - \psi_h(\mathbf{x}^i))\mathbf{x}^{ij}, \quad (19) $$

where $\mathbf{X}^i$ is the length distribution tensor defined as

$$ \mathbf{X}^i = \frac{1}{|\Sigma(i)|} \sum_{j \in \Sigma(i)} \mathbf{x}^{ij} \otimes \mathbf{x}^{ij}, \quad (20) $$

that gives an average representation of the distribution of the edges sharing an extremity.

### 4.3. Metric construction

In order to relate the error indicator $\epsilon_{ij}$ defined in (17) to a metric suitable for mesh adaptation purposes, we introduce the stretching factor $s_{ij}$ as the ratio between the length of the edge $\mathbf{x}^{ij}$ after and before the adaptation. The metric at node $\mathbf{x}^i$ is sought to generate unit stretched edge length in the metric space, that is,

$$ (s_{ij}\mathbf{x}^{ij})^T \cdot \mathbf{M}^i \cdot (s_{ij}\mathbf{x}^{ij}) = 1, \quad \forall j \in \Sigma(i), \quad (21) $$

for which an approximate least-square solution is shown in [49] to be

$$ \mathbf{M}^i = \left( \frac{d}{|\Sigma(i)|} \sum_{j \in \Sigma(i)} s_{ij}^2 \mathbf{x}^{ij} \otimes \mathbf{x}^{ij} \right)^{-1}, \quad (22) $$

provided the nodes in $\Sigma(i)$ form at least $d$ non co-linear edges with $\mathbf{x}^i$ (which is the case if the mesh is valid). The metric solution of (22) is ultimately computed setting a target total number of nodes $N_n$. Assuming a total error equi-distributed among all edges, the stretching factor is shown in [50] to be

$$ s_{ij} = \left( \frac{\sum_{i} N_i (1)}{N_n} \right)^{\frac{2}{d}} \epsilon_{ij}^{-1/2}, \quad (23) $$

where $N_i(1)$ is the number of nodes generated in the vicinity of node $\mathbf{x}^i$ for a unit error, given by

$$ N_i(1) = \left( \det \left( \frac{d}{|\Sigma(i)|} \sum_{j \in \Sigma(i)} \epsilon_{ij}^{1/2} \frac{\mathbf{x}^{ij}}{|\mathbf{x}^{ij}|} \otimes \frac{\mathbf{x}^{ij}}{|\mathbf{x}^{ij}|} \right) \right)^{-1/2} \quad (24) $$

### 4.4. Level set-based adaptation criteria

In order to simplify and clarify the presentation, the main steps needed for metric construction at the nodes is summarized in algorithm 2. In practice, the sole variable used for error estimation purpose is the filtered level set defined in (13), as it satisfies the metric property in a thin layer around the interface (in particular it preserves the zero iso-value of $\phi$, which is the only relevant information for mesh adaptation purposes), but avoids unnecessary adaption of the mesh further away from the interface (where the interpolation error is close-to-zero, due to $||\nabla\phi|| \sim 0$). This means that the criterion for mesh adaptation is purely geometric, i.e., the same mesh is pre-adapted around the fluid-solid interface, then used to compute all quantities needed to perform the next design update step. Nonetheless, it is worth mentioning that the approach also supports more complex adaptation criteria featuring physical quantities, thus providing the ability to dynamically adapt the mesh during the simulations. The common method to adapt a mesh to several
**Algorithm 2** Anisotropic mesh adaptation algorithm

**Require:** Anisotropic adapted mesh

1: Set number of nodes $N_n$

2: Compute $\psi_h$ on current mesh

3: **for** each node $\mathbf{x}^i$ **do**

4: Compute length distribution tensor $\mathbf{X}^i$ using (20)

5: Compute nodal recovered gradient. $\bar{\mathbf{g}}^i$ using (19)

6: **for** all edges $\mathbf{x}^{ij}$ **do**

7: Compute edge recovered gradient $\bar{\mathbf{g}}^{ij}$

8: Compute edge-based error $\varepsilon_{ij}$ using (17)

9: Compute stretching factor $s_{ij}$ using (23)

10: Compute metric $\mathbf{M}^i$ using (22)

11: Generate new mesh by local improvement in the neighborhood of the nodes and edges [51]

12: Interpolate $\psi_h$ on new mesh using classical linear interpolation

variables is to combine the metrics corresponding to each individual variable using metric intersection algorithms, which is known to incur a relatively high computational cost and to have potentially non-unique, suboptimal outcome. Conversely, the present approach allows building directly a unique metric from a multi-component error vector combining level set and any relevant flow quantity of interest, as definition (17) is easily extended to account for several sources of error [52]. Indeed, if we consider $\psi = (\psi_1, \psi_2, \dots, \psi_p)$ a vector consisting of $p$ scalar variables, it comes out straightforwardly that the error is now a vector $\varepsilon_{ij} = (\varepsilon_{ij,1}, \varepsilon_{ij,2}, \dots, \varepsilon_{ij,p})$, whose $L^2$ norm can serve as simple error value for the edge from which to compute the stretching factor (23) and ultimately, the metric solution of (22). For instance, the $2d+5$ sized nodal vector field defined as

$$ \psi_h(\mathbf{x}^i) = \left( \frac{\phi_h^i}{\max_{j \in \Sigma(i)} \phi_h^j}, \frac{u_h^i_{k \in \{1 \dots d\}}}{\|\mathbf{u}_h^i\|}, \frac{\|\mathbf{u}_h^i\|}{\max_{j \in \Sigma(i)} \|\mathbf{u}_h^j\|}, \frac{T_h^i}{\max_{j \in \Sigma(i)} T_h^j}, \frac{\tilde{u}_h^i_{k \in \{1 \dots d\}}}{\|\tilde{\mathbf{u}}_h^i\|}, \frac{\|\tilde{\mathbf{u}}_h^i\|}{\max_{j \in \Sigma(i)} \|\tilde{\mathbf{u}}_h^j\|}, \frac{\tilde{T}_h^i}{\max_{j \in \Sigma(i)} \tilde{T}_h^j} \right), \quad (25) $$

can be used to combine adaptivity with respect to the norm and direction of the state and adjoint velocity vectors and with respect to the state and adjoint temperatures, in addition to the level set. Because all fields are normalized by their respective global maximum, a field much larger in magnitude cannot dominate the error estimator, meaning that the variations of all variables are fairly taken into account. This benefits problems involving more complex physics (here, heat transfer, but also turbulence, fluid-structure interaction, multiple phases, possibly in interaction with one another), all the more so in the context of topology optimization, as the difference in the spatial supports of the state and adjoint quantities (due to the non-normality of the linearized evolution operator [53]) may otherwise yield conflicting requirements in terms of the regions of the computational domain most in need of refinement.

# 5. Numerical implementation

## 5.1. Geometrical constraints

Fluid flow topology optimization is generally performed under geometrical constraints, typically, constant or upper bounded surfaces and/or volumes to avoid the two extreme cases of the solid domain clogging the entire design domain (as in pressure drop minimization problems), or disappearing altogether (as in drag minimization problems). This is usually done adding penalty terms to the Lagrangian (each of which consists of an empirical penalty parameter multiplied by a measure of violation of the constraint), whose variations with respect to the state and design variables snowballs into the derivation of the adjoint problem and of the cost function sensitivity. Here, the constraint of a constant volume of fluid
$V_{target}$ is applied a posteriori, i.e., we solve the unconstrained problem presented in Sec. 2 (in the sense that no penalty term is added to the Lagrangian, although the optimization remains subject to Navier-Stokes as state equations). Once the convective level set method presented in Sec. J.5 has updated the interface position, a first pass of anisotropic mesh adaptation is performed, after which the volume of the fluid domain is computed as

$$V_{\varphi} = \int_{\Omega} H_{\epsilon}(\varphi) \, dv, \qquad (26)$$

where $H_{\epsilon}$ is the same smoothed Heaviside function on the fluid domain defined in (12). A simple dichotomy approach is then used to optimize a constant deformation $\delta\varphi$ meant to enlarge ($\delta\varphi < 0$) or shrink ($\delta\varphi > 0$) the fluid domain, until the difference $|V_{\varphi+\delta\varphi} - V_{target}|$ between the actual and target volumes drops below a certain tolerance, at which point we cut off $\varphi + \delta\varphi$ and perform a second pass of mesh adaptation. Two points are worth mentioning: first, because each offset changes the min-max values of the truncation, the above procedure requires knowledge of the level set $\varphi$ (not just the filtered level set $\phi$). A brute force algorithm therefore performs beforehand a complete reconstruction of the distance function from the zero iso-value of $\phi$, as only the filtered level set (not the level set) is evolved during the convection-reinitialization step. Second, only small deformations are considered so that no intermediate mesh adaptation passes are required. By doing so, the total cost has been found to be essentially that of performing the second pass of mesh adaptation (not shown here for conciseness).

## 5.2. Steepest descent update rule

In practice, the displacement used to perform the update step is defined as

$$\beta = -\theta \frac{\beta_0 \chi_{\Gamma}(\mathbf{x})}{\max_{\Omega} \beta_0 \chi_{\Gamma}(\mathbf{x}) \prod_{l} \zeta(\|\mathbf{x} - \mathbf{x}_s^l\|)}, \qquad (27)$$

where $\beta_0$ is the steepest descent estimate estimated from (9), $\theta > 0$ is a descent factor controlling the step taken in the gradient direction, and $\chi_{\Gamma}$ and $\zeta$ are activation functions between 0 and 1 ensuring that the design is fittingly updated only in relevant regions of the computational domain. More details are as follows:

* $\chi_{\Gamma}$ is a binary filter returning a value of 1 only at nodes within a distance $E$ of the interface. This is because the normal vector in a level set framework is recovered as $\mathbf{n} = \nabla\phi/\|\nabla\phi\|$, so the displacement is non-zero in the whole fluid domain, even far from the interface where $\mathbf{n}$ has unit norm because $\|\nabla\phi\|$ only tends asymptotically to zero. In return, the update step can break down numerically at nodes nearly equidistant from two subparts of the interfaces (for instance the centerline of a channel).

* $\zeta$ is a smooth filter assigning 0 value to some subset $X_s^l \in \partial\Omega$ that can be either a point or a curve, and $\|\mathbf{x} - \mathbf{X}_s^l\|$ is the shortest-path distance to $X_s^l$. Such subsets are singled out prior to optimization, because the flow there may be driven to a singularity, and ill-defined velocity gradients may cause large, unphysical displacements. Such singularities can be dealt with numerically by appending fluid/solid Dirichlet boundary conditions to the level set convection-reinitialization problem. Nonetheless, they must not be included in the normalization step to avoid forcing excessively small displacements along the remaining part of the interface, and thereby considerably slowing down the convergence rate of the iterative optimization process. We use here hyperbolic tangent filters

$$\zeta(r) = \frac{1}{2} + \frac{1}{2} \tanh \left( \alpha_s \tan \left( -\frac{\pi}{2} + \frac{\pi}{2} \frac{r}{r_s + \epsilon_{s1}} + \epsilon_{s2} \right) \right), \qquad (28)$$
increasing from 0 to 1 within a distance of $2r_s$ from the singularity, with $r_s$ a transition radius such that

$$4r_s < \min_{l,m} \|\mathbf{X}_s^l - \mathbf{X}_s^m\|, \qquad (29)$$

to prevent overlaps, $\alpha_s$ a steepness parameter controlling the sharpness of the transition, and $\epsilon_{s1,2}$ small regularization parameters to avoid local discontinuities.

Ultimately, the above filtering and normalization steps ensure that the level set is updated using a displacement that is non-zero only in a thin layer of thickness $E$ about the interface, minus a certain number of singular subsets.

### 5.3. Descent factor

It follows from Eq. (27) that the descent factor $\theta$ physically represents the maximum displacement amplitude over the update region of interest. In practice, though, the actual numerical displacement (estimated from the difference between zero iso-value of the filtered level set before and after transport) has been found to be well below its theoretical value. This is because the state and adjoint velocities are forced to zero is the solid domain, so the fluid component of the displacement (driven by the velocity gradients) is also zero everywhere in the solid, except in a very narrow region about the interface. As a result, it is not possible to explicitly control the displacement achieved numerically at each iteration. A simple scheme to do so would have been to repeatedly evolve the interface with a small descent factor until the difference between the cumulated and target displacement drops below a certain tolerance, but the interface can be evolved only once per update step, as the gradient information is lost if the displacement happens to be in the direction of the solid (for the same reason mentioned above). We thus tune the descent factor manually on a case by case basis, for the achieved displacement to be slightly smaller than the cut-off thickness. This has been found to be a satisfactory trade-off between accuracy and numerical effort, as the number of iterations required for convergence remains affordable, and the position of the evolved interface is accurately tracked (displacements larger than the cut off thickness conversely move the level set into regions of the computational domain lacking the proper mesh refinement, which has been found to ultimately affect the accuracy of the interface representation).

### 5.4. Parallel resolution

The resolution of the various physical problems considered herein (e.g., Navier-Stokes, adjoint Navier-Stokes, heat equation, adjoint heat equation, level-set advection) requires to compute, store and solve large-scale linear systems (or non-linear systems that may lead to the resolution of several linear systems if an implicit discretization scheme is used). To this end, the resolution step makes a clear distinction between those large systems that need to be stored and solved, and their local contributions at the element levels: all finite element formulations are only implemented sequentially at the element level, then assembled and solved in parallel using the PETSc library [54], that offers a wide range of parallel data structures (linear and non-linear solvers as well as preconditioners) and can be run on large computing clusters. Here, only semi-implicit and explicit discretization schemes are used; see J), and the associated linear systems are sufficiently well conditioned to be solved by iterative methods. We thus use the Generalized Minimal Residual (GMRES) algorithm with block Jacobi incomplete LU preconditioning, and consider the solutions to be converged if the absolute residuals are less than $10^{-6}$.

### 5.5. General algorithm

Figure 1 shows the flowchart of the implemented topology optimization algorithm, in which anisotropic mesh adaptation is key to capture the interface with the highest precision possible. The necessary algorithmic parameters common to all examples documented in the following are given in Tab. 1. Note, as a consequence of the level set-based technique used
<figure><img src="image_2.png" /><figcaption>Figure 1. Flowchart of performance topology optimization procedure.</figcaption></figure>

to enforce the volume of fluid constraint, convergence is achieved not when the displacement
is identically zero (as would be the case using a penalized Lagrangian approach), but when
the displacement is uniform along the interface. This is not easily done on the fly, though,
so we rather iterate until a maximum number of iterations has been reached and evaluate
convergence a posteriori.

**6. Numerical benchmarks**

6.1. Preliminaries

This section assesses the efficiency of the numerical framework through a series of
topology optimization problems showcasing the accuracy to which the optimal interfaces are
captured in the simulation model. Several cases are considered in 2-D and 3-D, for which we
aim at maximizing heat transfer in the domain while minimizing the total pressure losses in
the fluid channels. A bi-objective optimization strategy is used, that consists in minimizing
the linear weighted sum (10) of two criteria inspired from [37] and [55], namely an hydraulic
cost function J<sub>ν</sub> measuring the net inward flux of total pressure through the boundaries (to
minimize) and a thermal cost function J<sub>Θ</sub> measuring the recoverable thermal power from
the domain through the inlet and outlet flow boundary conditions (to maximize). Since the
<table><caption>**Table 1.** Algorithmic parameters.</caption><tr><td>h<sub>⊥</sub> = 0.0001</td><td>Min. interface normal mesh size</td></tr><tr><td>Δt = 0.1</td><td>CFD Numerical time step</td></tr><tr><td>E = 0.005</td><td>Level set cut off thickness</td></tr><tr><td>|δφ| = [0.0005; 0.001]</td><td>Initial volume recovery offset</td></tr><tr><td>r<sub>s</sub> = 0.0125</td><td>Transition radius</td></tr><tr><td>α<sub>s</sub> = 2.1</td><td>Sharpness parameter</td></tr><tr><td>(ε<sub>s1</sub>, ε<sub>s2</sub>) = (0.0005, 0.005)</td><td>Regularization parameters</td></tr></table>

orientation of the normal **n** yields **u**·**n**|<sub>Γ<sub>i</sub></sub> > 0 and **u**·**n**|<sub>Γ<sub>o</sub></sub> < 0, this is expressed in the form of (5) as

$$J_{\nu} = p_{tot}(\mathbf{u} \cdot \mathbf{n}) = \left(p + \frac{1}{2} \rho (\mathbf{u} \cdot \mathbf{u})\right) (\mathbf{u} \cdot \mathbf{n}), \quad \text{and} \quad J_{\Theta} = \rho c_p T (\mathbf{u} \cdot \mathbf{n}). \quad (30)$$

In all cases, a reference design domain is chosen under the form of a cubic or cuboid (parallelepipedic) cavity, with cylindrical inlet and outlet at which parabolic profiles normal to the boundary are prescribed, as defined by

$$\mathbf{u}_{i,o} = u_{i,o} \left( 1 - \frac{4r^2}{e_{i,o}^2} \right) \mathbf{n}, \quad (31)$$

with *r* the distance to the inlet/outlet center, *e*<sub>*i*</sub> (resp. *e*<sub>*o*</sub>) the inlet diameter (resp. the outlet diameter) and *u*<sub>*i*</sub> (resp. *u*<sub>*o*</sub>) the inlet centerline velocity (resp. the outlet centerline velocity, adjusted for the mass flow exiting through the outlet to match exactly that entering through the inlet). For each case, the control parameters are the Reynolds (built here on inlet diameter and maximum inlet velocity) and Prandtl numbers.

The remainder of the practical implementation details are as follows:

* All design domains are initialized with solid inclusions coming in various shapes and sizes. From experience, the flow topology optimization problems tackled in the following are essentially insensitive to the initial design provided a sufficient large number of inclusions is used (additional mechanisms for seeding solid inclusions could be added to the proposed framework, but however lie outside the scope of this study).

* The admissible error on the target volume is set to 1% in 2-D, and 5% in 3-D.

* The fluid is systematically conveyed into and out of the design domain using leads of length *l*<sub>*i*</sub> (the same at all inlets) and *l*<sub>*o*</sub> (the same at all outlets) appended normal to the boundary. This is for numerical consistency, as the exact problem formulation in the literature may vary depending on the case, and it is not always clear whether such leads should be included in the design domain (which they are here, although they are not considered in the volume constraint, neither in definition of the target volume nor in the computation of the volume of fluid).

* The singular subsets excluded from the displacement normalization step are the sharp intersections between the leads and the boundary of the cavities. Note, this is not a consequence of explicitly representing the leads, as the exact same procedure has been found suitable without such appendage. In practice, since all inlets and cylindrical outlets are cylindrical, each smooth filter ζ therefore transitions from 0 to 1 over a circle of radius 2r<sub>s</sub> (in 2-D) or a torus of minor radius 2r<sub>s</sub> and major radius equal to the inlet/outlet radius (in 3-D).

* The leads are excluded from the displacement normalization step, for which we simply add to the max argument of (27) a binary filter returning a value of 0 at all nodes located inside the pipes. This is again to avoid slowing down the convergence rate of the iterative optimization process, as the maximum displacement is otherwise located
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Set-up of the two-dimensional single pipe problem with heated walls. The orange and light gray shade denote hot isothermal and adiabatic walls, respectively.</figcaption></figure>

in the leads (because the easiest way to minimize the dissipated power is to suppress the flow by having the solid entirely clogging the leads).

* Without seeking to optimize the performance, all optimization runs have been found to converge within a few hundreds iterations, which is essentially the number of steps used to fulfill the fluid volume constraint (more details in the following) while ensuring that the displacement achieved at each iteration remains below the level set cut-off thickness.

* All 3-D meshes (resp. 2-D meshes) have been checked to have an element-to-node ratio close to 5 (resp. close to 2), as should be for denses mesh made up of tetrahedral (resp. triangular) elements. In order to ease the comparison with the available literature, the mesh information is thus documented in the following in terms of its equivalent number of elements, defined as $N_{el} = 5N_n$ (resp. $N_{el} = 2N_n$).

## 6.2. Two-dimensional splitter device with two outlets

In order to provide a first verification and characterization of the method, we consider a purely hydraulic ($\omega = 0$) problem aiming at minimizing power dissipation in a 2-D flow splitter device, a classical example of functional unit used for continuous separation and collection of particles in microchannels [56], and a class of flow that has received attention as a relevant example of topology optimization in fluid dynamics [57,58]. The design domain sketched in Fig. 2 is a square cavity of unit height, with a single inlet at the bottom and two outlets at the left size on the left and right sides. The aim is to determine the optimal design that connects the inlet to the outlets, subject to the constraint that the fluid must occupy 20% percent of the cavity, and the flow must be distributed evenly over the multiple outlet orifices for each outlet to have 1/2 of the fluid flow entering through the inlet; see Tab. 2 for provision of the remaining problem parameters. The initial design shown in Fig. 4 consists of spherical occlusions arranged for the initial volume of fluid to fill about 50% of the cavity, in violation of the volume constraint. This is because many more smaller inclusions are needed to recover the proper volume, which would either dramatically increase the surface of the interfaces that needs be captured (and thus, the number of mesh elements needed to maintain the numerical accuracy), or risk clogging the fluid path due to insufficient mesh refinement. As shown in the convergence history presented in Fig. 3, there is thus an initial transient during which the cost function, albeit low, has little physical meaning, as the constraint value is decreased up to the point where it reaches the target within the desired tolerance, and the cost function adjusts until a feasible minimum is found.
<table><caption>**Table 2.** Numerical parameters for the two-dimensional flow splitter device problem.</caption><tr><td>Ω = [0; 1]×[0; 1]</td><td>Design domain</td></tr><tr><td>d = 2</td><td>Problem dimensionality</td></tr><tr><td>V<sub>target</sub> = 0.25</td><td>Target volume of fluid</td></tr><tr><td>V<sub>φ,0</sub> = 0.50</td><td>Initial volume of fluid</td></tr><tr><td>Re = 2</td><td>Reynolds number</td></tr><tr><td>u<sub>i</sub> = 0.3</td><td>Inlet centerline velocity</td></tr><tr><td>e<sub>i</sub> = 0.125</td><td>Inlet diameter</td></tr><tr><td>l<sub>i</sub> = 0.4</td><td>Inlet leads length</td></tr><tr><td>x<sub>i1</sub> = (0.8125, -0.4)</td><td>Inlet center coordinates</td></tr><tr><td>u<sub>o1</sub> = 0.1</td><td>Outlet 1 centerline velocity</td></tr><tr><td>u<sub>o2</sub> = 0.1</td><td>Outlet 2 centerline velocity</td></tr><tr><td>e<sub>o</sub> = 0.1875</td><td>Outlet diameter</td></tr><tr><td>l<sub>o</sub> = 0.125</td><td>Outlet leads length</td></tr><tr><td>x<sub>o1</sub> = (-0.125, 0.21875)</td><td>Outlet 1 center coordinates</td></tr><tr><td>x<sub>o1</sub> = (-0.125, 0.78125)</td><td>Outlet 1 center coordinates</td></tr><tr><td>N<sub>n</sub> = 30000</td><td>Nb. mesh nodes</td></tr><tr><td>N<sub>el</sub> = 60000</td><td>Nb. mesh elements</td></tr></table>

A total of 500 iterations has been run with 30000 mesh elements. The optimization goes through several complex stages, e.g., merging or cancellation of holes, all accurately represented on anisotropic adapted meshes, as evidenced by the selected samples shown in Fig. 4. Also, all meshes exhibit the expected refinement and deformation, with coarse and regular elements away from the interface between solid and fluid (all the more so in the solid domain, where only a few ten elements are used), but fine, extremely stretched elements on either side of the interface (for the velocity to smoothly transition to zero across the boundary layer). In return, the interfaces are sharply captured, not only at optimality but during all stages of the optimization (even in the leads). Ultimately, the optimal duct for this case is a wide pipe splitting at mid length into two almost identical (in terms of diameter), thinner pipes, each connecting to an outlet. This layout stands as the better trade-off between transporting fluid the shortest way, and transporting it in the widest possible pipe, and is consistent with the results documented in Ref. [57], although the optimal shapes therein exhibit quality issues (staircase effects) in smoothly curved regions, and anisotropic mesh adaptation represents a tremendous improvement in this regards. The same trends are observed when assessing the sensitivity of the designed layout to variations in the design conditions, such as the volume of fluid allowed in the layout, or the ratio of flow rate at the outlets. This is evidenced in Fig. 5(a) showing the final ducts optimized for several volume of fluid constraints in a range from 10 to 50%. Overall, the obtained results show that decreasing the volume of fluid simply slims down the optimal channels, but increases the dissipated power, which is because most energy is dissipated by shear at low Reynolds numbers, so an optimal flow pipe is preferably as short and wide as possible (note, the analysis overlooks the mass of the layout, that obviously puts a limit on how much fluid should be allowed in a practical device). Similarly, we show in Fig. 5(b) the final ducts optimized for several flow rate ratios ranging from 1:9 (i.e., 10% at the lower outlet vs. 90% at the upper outlet) to 9:1. The results show that the larger the flow rate at a certain outlet, the thicker the channel, which is fully consistent with the previous findings. Ultimately, less power is dissipated if the flow rate is larger at the lower outlet, with the 1:9 optimal lower than its 9:1 counterpart by 25%. This is simply because a larger lower flow rate makes for a shorter path from the inlet to the outlet, and the cost of bending the fluid stream shortly after the inlet is low, given that most fluid flows in the (shorter) inner region.
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Convergence history for the two-dimensional flow splitter device presented in Fig. 2. All cost function values made non dimensional using the inlet diameter and maximum inlet velocity (equivalently, using $\rho u_i^3 e_i$ as reference cost functional value). The grey shade indicates the iterations during which the volume constraint is adjusted to its target value.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Hydraulic optimization of the two-dimensional flow splitter device presented in Fig. 2. From left to right and from top to bottom: the zero iso-value of the level set function and associated anisotropic adapted meshes are sampled over the course of optimization using the parameters given in Tab. 2. The associated volume of fluid are 49.5%, 45.4%, 40.6%, 34.7%, 24.8% and 24.9%, respectively.</figcaption></figure>

### 6.3. Two-dimensional single pipe with heated walls

The second case study is a two-dimensional conjugate heat transfer problem that has received substantial attention in the recent literature [26,32,55,59]. The design domain shown in Fig. 6 is a square cavity of unit height. It has a single inlet on the left side and a single outlet on the right side lined up in front of each other, and is discretized with 50000 elements. A cold fluid is flowing from the inlet, and is heated by the top and bottom walls, subject to a fixed (hot) temperature. All other walls (cavity and leads) are insulated from the surroundings with zero heat absorbed or released (i.e., adiabatic). The solid is set to be 10 times more diffusive than the fluid, which allows using fluid to insulate thermically inner regions from the cold inlet temperature. The aim is to determine the optimal design that connects the inlet to the outlets subject to the constraint that the fluid
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Sensitivity of the optimal design against (a) the volume of fluid allowed in the design domain, and (b) the fraction of fluid leaving the design domain through the lower outlet, for the two-dimensional flow splitter device presented in Fig. 2. All cost function values made non-dimensional using the inlet width and centerline inlet velocity (equivalently, using $ρu_i^3e_i$ as reference cost functional value).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Set-up of the two-dimensional single pipe problem with heated walls. The orange and light gray shade denote hot isothermal and adiabatic walls, respectively.</figcaption></figure>

must occupy 40% percent of the cavity (twice as mush as the straight parallel pipe fitting exactly to the inlet and outlet). All other problem parameters, including Reynolds and Prandtl numbers, are given in Tab. 3. Note, because the inlet and outlet diameters are the same, mass conservation demands the same velocity condition to be prescribed at the inlet and outlet. We do not share the view expressed in [26] that this is ill-posed, in the sense that it does force the algorithm to identify acceptable trade-offs between both hydraulic and thermal for heat and mass transfer optimization without yielding broken flow paths, dead ends or non-physical artifacts, which is the desired goal. Moreover, the argument that the optimization is limited by the fact that the sole variables left to optimization are the
<table><caption>**Table 3.** Numerical parameters for the two-dimensional single pipe problem with heated walls.</caption><tr><td>Ω = [0; 1]×[0; 1]</td><td>Design domain</td></tr><tr><td>d = 2</td><td>Problem dimensionality</td></tr><tr><td>V<sub>target</sub> = 0.4</td><td>Target volume of fluid</td></tr><tr><td>V<sub>φ,0</sub> = 0.4</td><td>Initial volume of fluid</td></tr><tr><td>Re = 4.5</td><td>Reynolds number</td></tr><tr><td>Pr = 5</td><td>Prandtl number (fluid)</td></tr><tr><td>k<sub>s</sub>/k<sub>f</sub> = 10</td><td>Solid to fluid thermal diffusivity ratio</td></tr><tr><td>T<sub>w</sub> = 10</td><td>Hot wall temperature</td></tr><tr><td>u<sub>i</sub> = 1</td><td>Inlet centerline velocity</td></tr><tr><td>T<sub>i</sub> = 0</td><td>Inlet cold temperature</td></tr><tr><td>e = 0.2</td><td>Inlet diameter</td></tr><tr><td>l<sub>i</sub> = 0.1</td><td>Inlet leads length</td></tr><tr><td>x<sub>i</sub> = (-0.1, 0.5)</td><td>Inlet center coordinates</td></tr><tr><td>u<sub>o</sub> = 1</td><td>Outlet centerline velocity</td></tr><tr><td>e<sub>o</sub> = 0.2</td><td>Outlet diameter</td></tr><tr><td>l<sub>o</sub> = 0.1</td><td>Outlet lead length</td></tr><tr><td>x<sub>o</sub> = (1.1, 0.5)</td><td>Outlet 1 center coordinates</td></tr><tr><td>N<sub>n</sub> = 25000</td><td>Nb. mesh nodes</td></tr><tr><td>N<sub>el</sub> = 50000</td><td>Nb. mesh elements</td></tr></table>

pressure (not total pressure) drop between the inlet and outlet, and the outlet temperature (which removes the need to explore, e.g., converging or narrowing channels designs aimed at increasing the fluid velocity), while true in this particular setting, does not hold if multiple inlets/outlets or different inlet/outlet diameters are used.

We show in Fig. 7 distinct optimal designs computed by increasing progressively the thermal weigh, to which we come back below for further discussion. At this stage, the point of interest is that all related optimization runs go through several complex stages all accurately represented on anisotropic adapted meshes made up of extremely stretched elements on either side of the interface. This is evidenced by the selected samples shown in Figs. 8-9, where the main difference compared to the hydraulic case in Sec. 6.2 is the finer element size used to discretize the inner solid domain (here the same as in the inner fluid domain) to accurately resolve heat conduction. This is all the more important given that the increased non-linearity of the optimization problem at large thermal weighs (where there is almost no contribution from the hydraulic cost function) yields strongly anisotropic material distributions, that require adequately capturing the formation and destruction of very fine cross-flow fluid structures whose diameter can be below 1/20 the inlet diameter.

For low thermal weighs, the optimal design is a single, straight pipe connecting the inlet to the outlet, as evidenced in Fig. 8(a). This is because the contribution of the thermal cost function is negligible, so the only objective is to minimize the fluid power dissipation, hence short and wide pipes. In return, the optimal pipe is as wide as allowed by the volume constraint, with increased cross section halfway though (compared to the inlet and outlet) to minimize shear. Increasing $\omega$ adds more priority to increase the recoverable thermal power, which opens the possibility to depart from the straight pipe even at the cost of some increase in the hydraulic objective. In this regards, our results highlight the existence of two distinct branches of solution, referred to as symmetric and asymmetric.

Symmetric designs feature a solid core forming at the center of the cavity, hence dividing the lead into a lower and an upper pipe, as shown in Fig. 8(b). This increases both the dissipated power and the recoverable thermal power, as it lengthens the distance over which the fluid travels, but moves both pipes towards the hot walls, and heats up the fluid without any temperature losses associated to heat flux conduction through the solid (as the solid thermal conductivity is 10 times that of the fluid). The asymmetric designs presented in Fig. 9(a) conversely feature a single pipe bending into either the lower or the upper half
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Optimal designs sorted by weighting $\omega$ for the two-dimensional single pipe device with heated walls presented in Fig. 6. (a) Low and intermediated thermal weighs. The dashed lines denote the results obtained at large thermal weighs, further presented by the close-ups in (b-c) for (b) symmetric and (c) asymmetric designs.</figcaption></figure>

of the domain, which is a different trade-off involving both less recoverable thermal power and less dissipated power, as the fluid is heated up at only one out of the two hot walls, but travels in a wider pipe.

Increasing the thermal weigh forces the fluid along the hot walls to expand the exchange surface. Beyond a certain threshold, the symmetric solid core splits vertically into an increasing number of subcores, as a network of fluid strips forms to act as a large thermal resistance breaking the horizontal temperature gradient to reduce the core heat conduction; see Fig. 8(c). Meanwhile, in Fig. 9(b), the asymmetric bent pipe gives way in to the more complex Z-shaped pipe successively forcing the flow along the top and bottom walls (yet another trade-off that increases both the recoverable thermal power and the dissipated power, as the heat exchange surface doubles, but the fluid travels in a thinner and longer pipe.), whose solid layout eventually fragments vertically near the Z edges to make the most of the low conductivity of the fluid; see Fig. 9(c).

Upon comparing the above design to those in Ref. [55] (the closest study to our work in the available literature), the following remarks can be made:
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Multi-objective optimization of the two-dimensional single pipe device with heated walls presented in Fig. 6. From top to bottom: the zero iso-value of the level set function and associated anisotropic adapted meshes are sampled over the course of optimization using the parameters given in Tab. 3. (a) Straight pipe solution with $\omega = 0.4$. (b) Solid core solution with $\omega = 0.7$. (c) Fragmented core solution with $\omega = 0.987$.</figcaption></figure>

* anisotropic adapted meshes dramatically improve the accuracy of all geometric representations, as most results in the recent available literature exhibit obvious staircase effects in all curved regions.
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Multi-objective optimization of the two-dimensional single pipe device with heated walls presented in Fig. 6. From top to bottom: the zero iso-value of the level set function and associated anisotropic adapted meshes are sampled over the course of optimization using the parameters given in Tab. 3. (a) Bent pipe solution with $\omega = 0.85$. (b) Z pipe solution with $\omega = 0.97$. (c) Fragmented Z pipe solution with $\omega = 0.998$.</figcaption></figure>

* the solid and fragmented core solutions are generally consistent, although they show up in [55] at surprisingly much lower values of $\omega$ (an issue already raised in [26]) and with lesser horizontal symmetry at large thermal weights (which may be because the
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Fragmented core solutions computed under various initial designs under thermal weigh $\omega = 0.987$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 11.** Pareto frontier for bi-objective topology optimization of the single pipe device with heated walls. The red symbols are the near-optimal design for which both the dissipated power and the recoverable thermal power approach their single objective optimization value.</figcaption></figure>

authors in the aforementioned study do not impose a specific target volume of fluid, but only an upper bound).

* asymmetric designs are noticeably absent from [26,55] and from other studies tackling variations of this problem. Again the explanation may lie in the constraint on the maximum volume of fluid, in the sense that for a given asymmetric design minimizing the cost function under a certain thermal weigh, a more efficient symmetric design may exist at a smaller volume of fluid. Asymmetric designs are reported in [32], for which the authors allude to the use of unstructured meshes, but we believe they are rather the consequence of different flow regimes, as the aforementioned study considers a much higher Reynolds number of 400 and a much lower Prandtl number of 0.05, and the present use of unstructured meshes does not alters the solutions symmetry.

* similar formation fluid strips to act like a heat insulation material at large $\omega$ is documented in [55]. While it is a robust mechanism, in the sense that even a non-fragmented solid cores and Z pipes computed at a slightly smaller thermal weigh end up breaking up, we have found the way the solid layout splits to be very sensitivity to the optimization path. This is evidenced in Fig. 10 showing a series of fragmented core designs generated by varying the initial design under constant thermal weigh. The number of subcores and the subcores arrangements does vary, but they yield identical cost functions (and thermal cost functions) within 1%, which suggests that fragmented solutions are actually flat minimizers.
<figure><img src="image_13.png" /><figcaption>**Figure 12.** Set-up of the three-dimensional single pipe problem with heated walls: one single hot stripe. The orange and light gray shade denote hot isothermal and adiabatic walls, respectively.</figcaption></figure>

* non-physical designs (e.g., broken flow paths, dead ends and total fluid flow blockage, all highly undesirable from manufacturing point of view) have been obtained at much higher values of $\omega > 0.999$, which is when the contribution of the hydraulic cost function becomes negligible. This may have to do with the modeling of the solid material, as a porous media does allow solutions with no fluid connection between the inlet and outlet flows, while the IVM rigorously forces the solid velocity to zero.

Additional testing has been performed in the attempt to clarify the connection between the various branches of solutions (although not the intended scope of this study). The main findings are threefold: first, the straight pipe solution exists up to $\omega \sim 0.7$, after which its centerline shifts increasingly in the upper domain (or lower domain, by vertical reflectional symmetry) and the solution evolves continuously into the bent pipe solution. Second, the solid core solution branches off the straight pipe at about $\omega \sim 0.64$ (this has been estimated using a branch-tracking technique in which the design is initialized with a solid core solution computed at a slightly larger thermal weigh), then evolves continuously into the fragmented core solution. Finally, the connection between the bent and Z pipe solutions remains uncertain: both solutions have been found to coexist over a range of thermal weighs from 0.94 to 0.98, where they yield almost identical cost functions. Meanwhile, we could not manage to have a bent pipe continuously turn into a Z, regardless of the value of $\omega$ and the number of update steps (up to several ten thousands). This raises the possibility that the Z pipe solution may branch off subcritically from the bent pipe solution, leading to hysteresis (testing this hypothesis is uneasy due to the difficulty of consistently generating Z pipe solutions unless a branch-tracking technique is used, which is why the optimization run documented in Fig. 9(b) does not start from the classical design with solid occlusions).

Finally, Fig. 11 recasts the obtained results into the Pareto frontier of the multi-objective optimization problem. Interestingly, it turns out that the Pareto-efficient subset consists exclusively of straight, bent and Z pipes (although the fragmented core solutions equally dominate at very high influences of the thermal objective function). Interestingly, the close-to-convex shape of the Pareto front means that a few solutions provide an acceptable trade-off by having both single cost functions close to their single objective optimization. This corresponds here to the bent pipe at $\omega = 0.97$ and the Z pipe at $\omega = 0.95$, for which the recoverable thermal power is below its single objective maximum by less than 10% (the dissipated power is 7 times as large as its single objective minimum, which is not small strictly speaking, but very reasonable given that the worst performance is actually by a factor of 100).

### 6.4. *Three-dimensional single pipe with heated walls*

The last test case stands as a three-dimensional counterpart of the two-dimensional conjugate heat transfer problem considered in Sec. 6.3. The setting inspired from [17] is depicted in Fig. 12, with detailed problem parameters given in Tab. 4. The design domain is a cuboid cavity of unit height and aspect ratio 2:1:1, with a single inlet on the left and a single outlet on the right, again lined up in front of each other. A cold fluid is flowing from
<table><caption>**Table 4.** Numerical parameters for the three-dimensional single pipe problem with heated walls.</caption><thead><tr><th></th><th></th><th></th></tr></thead><tbody><tr><td>Ω = [0; 2]×[0; 1]×[0; 1]</td><td>»</td><td>» Design domain</td></tr><tr><td>d = 3</td><td>»</td><td>» Problem dimensionality</td></tr><tr><td>V<sub>target</sub> = 0.4</td><td>»</td><td>» Target volume of fluid</td></tr><tr><td>V<sub>φ,0</sub> = 1.8</td><td>»</td><td>» Initial volume of fluid</td></tr><tr><td>Re = 12</td><td>»</td><td>» Reynolds number</td></tr><tr><td>Pr = 83.5</td><td>»</td><td>» Prandtl number</td></tr><tr><td>k<sub>s</sub>/k<sub>f</sub> = 10</td><td>»</td><td>» Thermal conductivity ratio</td></tr><tr><td>T<sub>w</sub> = 10</td><td>»</td><td>» Hot stripe temperature</td></tr><tr><td>Δx<sub>h</sub> = 0.2</td><td>»</td><td>» Hot stripe width</td></tr><tr><td>x<sub>h1</sub> = 1</td><td>0.5</td><td>» Hot stripe 1 center coordinate</td></tr><tr><td>-</td><td>-</td><td>x<sub>h2</sub> = 1.5 Hot stripe 2 center coordinate</td></tr><tr><td>u<sub>i</sub> = 1</td><td>»</td><td>» Inlet centerline velocity</td></tr><tr><td>T<sub>i</sub> = 0</td><td>»</td><td>» Inlet cold temperature</td></tr><tr><td>e = 0.2</td><td>»</td><td>» Inlet diameter</td></tr><tr><td>l<sub>i</sub> = 0.1</td><td>»</td><td>» Inlet leads length</td></tr><tr><td>x<sub>i</sub> = (-0.1, 0.5, 0.5)</td><td>»</td><td>» Inlet center coordinates</td></tr><tr><td>u<sub>o</sub> = 1</td><td>»</td><td>» Outlet centerline velocity</td></tr><tr><td>e<sub>o</sub> = 0.2</td><td>»</td><td>» Outlet diameter</td></tr><tr><td>l<sub>o</sub> = 0.1</td><td>»</td><td>» Outlet lead length</td></tr><tr><td>x<sub>o</sub> = (2.1, 0.5, 0.5)</td><td>»</td><td>» Outlet 1 center coordinates</td></tr><tr><td>N<sub>n</sub> = 1000000</td><td>»</td><td>» Nb. mesh nodes</td></tr><tr><td>N<sub>el</sub> = 5000000</td><td>»</td><td>» Nb. mesh elements</td></tr></tbody></table>

the inlet, and is heated by the cavity walls, with the difference that only a finite stripe at the middle of the cavity walls is maintained at a constant (hot) temperature, and all remaining walls (cavity and leads) are considered adiabatic. In what follows, the thermal weigh is set $\omega = 0.95$ to add more priority to increase the recoverable thermal power. Note, although the configuration has two reflectional symmetries, we do not reduce the computational cost by modeling only a quarter of the domain together with symmetry boundary conditions, which is feasible [60] but would not allow assessing the method in the context of large scale systems. The entire domain is thus discretized with 5000000 mesh elements, and we let symmetry eventually arise as a result of the optimization process.

The initialization shown in Fig. 13 corresponds to a fluid box filled with islands of solid spherical inclusions occupying about 10% of the cavity. The fluid thus fills initially 90% of the cavity, well above the 20% volume constraint, hence the same approach as in Sec. 6.2 is used, in which the constraint value decreases during the early stage of the optimization process, up to the point where it reaches the target within the desired tolerance. Again, the method is found to handle well the various topological changes occurring over the course of optimization, and all anisotropic adapted meshes exhibit extremely stretched elements regardless of the interface complexity, that allow sharply representing the fluid and solid domains (also, the edges of the hot stripe) and accurately computing the solutions during all stages of optimization. In the optimal solution shown in Fig. 13, the fluid flows in and out of the cavity through single, straight pipes. This is because the hot stripe is far from the inlet/outlet sections, so there is a good proportion of the cavity where the thermal cost function contributes little to nothing, and the best trade-off is to minimize power dissipation. Similarly to what could be observed in 2-D, a solid core forms in the stripe region. The latter divides the inlet pipe into a near-perfect symmetrical network of 8 pipes, that quickly merge themselves to deliver the fluid to the outlet via a complex 4-element
<figure><img src="image_14.png" /><figcaption>**Figure 13.** Multi-objective optimization of the three-dimensional single pipe device with heated walls (one single hot stripe) presented in Fig. 12. The zero iso-value of the level set function and associated anisotropic adapted meshes are sampled over the course of optimization using the parameters given in Tab. 4. The associated volume of fluid (from top to below) is as follows: 90%, 42.9%, 20.6%, 20.5% and 20.4%, respectively.</figcaption></figure>

comb-like arrangement (one per face of the cavity). This forcesthe fluid along the hot walls and eventually merging into a thin, square annulus shaped to the hot stripe to maximize the fluid heat up. As illustrated in Fig. 14, the optimal features thin inclusions of fluid attached to the main pipes. This is essentially reminiscent of the two-dimensional fragmentation mechanism observed at such large thermal weigh, where fluid is used to insulate thermically the inner pipes from the cold inlet temperature. The present optimal design is overall close to that documented in [17], but the pipe arrangements differ in the hot region, most likely
<figure><img src="image_15.png" /><figcaption>**Figure 14.** Optimal design of the three-dimensional single pipe device with heated walls (one single hot stripe) presented in Fig. 12, together with stream-wise and cross-wise cuts at positions shown by the grey planes. The colors hue corresponds to the temperature profile.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 15.** Set-up of the three-dimensional single pipe problem with heated walls: one single upstream hot stripe vs. two hot stripes. The orange and light gray shade denote hot isothermal and adiabatic walls, respectively.</figcaption></figure>

because the authors in the aforementioned reference optimize the thermal recoverable power by imposing an upper bound threshold for the pressure drop (which yields to a different trade-off).

Two other cases have been considered to assess the capability of designing more complex shapes by giving more importance to the thermal cost function. The associated setups depicted in Fig. 15 differ by the number and position of hot stripes, namely the first case (case 1) has one stripe shifted upstream against the inlet, and the second one (case 2) has two stripes against the inlet and outlet arranged symmetrically with respect to the middle of the cavity; see Tab. 4 for provision of other detailed problem parameters and Figs. 16-17 for illustration of the corresponding optimization runs using anisotropic adapted meshes. For case 1, the main features of the baseline optimal design hereinabove carry over, with the difference that the solid core moves upstream to follow in the footsteps of the hot stripe, hence the inlet lead immediately splits into a similar network of 8 pipes. For case 2, the presence of two separated hot spots yields a different optimal, with the 8 pipes reconfiguring into 4 wider pipes to transport fluid the shortest way downstream and avoid the cost of bending. These pipes then widen to form four quasi triangular prisms shaped to the downstream stripe (to maximize the fluid heat up) before merging to connect to the outlet. In both cases, the optimal designs shown in Fig. 18 exhibit the same thin inclusions of fluid attached to the main pipes to benefit from the insulating low-conductivity of the fluid.

# 7. Discussion

## 7.1. *Numerical cost*

It is worth noticing that the number of nodes used here is actually quite large and mostly useful during the early stage of optimization. This is because the surface of the interfaces (perimeter in 2-D) that needs be captured is initially dramatically large to the
<figure><img src="image_17.png" /><figcaption>**Figure 16.** Multi-objective optimization of the three-dimensional single pipe device with heated walls (one single upstream hot stripe) presented in Fig. 15(a). The zero iso-value of the level set function and associated anisotropic adapted meshes are sampled over the course of optimization using the parameters given in Tab. 4. The associated volume of fluid (from top to below) is as follows: 90%, 43.0%, 20.4%, 20.5% and 20.4%, respectively.</figcaption></figure>

many solid inclusions, then decreases substantially after the first dozens of iterations, as has been found computing the surface area

$$S_{\varphi} = \int_{\Omega} \delta_{\epsilon}(\varphi) \, dv, \qquad (32)$$
<figure><img src="image_18.png" /><figcaption>**Figure 17.** Multi-objective optimization of the three-dimensional single pipe device with heated walls (two hot stripes) presented in Fig. 15(b). The zero iso-value of the level set function and associated anisotropic adapted meshes are sampled over the course of optimization using the parameters given in Tab. 4. The associated volume of fluid (from top to below) is as follows: 90%, 54.4%, 31.9%, 20.6% and 20.3%, respectively.</figcaption></figure>

where $\delta_{\epsilon}$ is the Dirac function

$$ \delta_{\epsilon}(\varphi) = \begin{cases} \frac{1}{2\epsilon}\left(1 + \cos\left(\pi\frac{\varphi}{\epsilon}\right)\right) & \text{if } |\varphi| \le \epsilon, \\ 0 & \text{if } |\varphi| > \epsilon, \end{cases} \qquad (33) $$

smoothed with the same regularization parameter $\epsilon$ as the Heaviside function (12). Also, the anisotropic mesh adaptation algorithm refines the mesh in hierarchical importance of the

789

790

791
<figure><img src="image_19.png" /><figcaption>**Figure 18.** (a) Optimal design of the three-dimensional single pipe device with heated walls (one single upstream hot stripe) presented in Fig. 15(a), together with stream-wise and cross-wise cuts at positions shown by the grey planes. (b) Same as (a) for the device with heated walls (two hot stripes) presented in Fig. 15(b). The colors hue corresponds to the temperature profile.</figcaption></figure>

level set gradient. If new geometrical features appear in the solution (associated with high gradients), the mesh is automatically coarsened in regions with lower gradient and refined near the newly emerging features. If the number of nodes is large (as has been the case so far), the decrease in the interface surface area allows resolving finer, more complex patterns without degrading the accuracy in other parts of the design domain, because the coarsened regions are actually over-resolved. This shows through the progressive mesh refinement in the fluid domain in the various figures, as more and more elements become available to improve the mesh in other regions of the domain.

Figure 19 presents detailed timing results obtained by averaging (and normalizing to achieve unit average time per iteration) dedicated update steps of the various multi-objective conjugate heat transfer cases performed on 64 cores (200 steps in 2-D, 100 steps in 3-D). In 2-D, the cost of an iteration is dominated by that of computing the state solution (which takes about 10 Navier-Stokes iterations representing 55% of the total cost), and otherwise by that of adapting the mesh (about 20% of the total cost). Using the same number of processors, the cost of a 3-D iteration is larger than its 2-D counterpart by roughly three orders of magnitude, the cost of which is essentially that of the two passes of mesh adaptation (about a cumulative 70% of the total cost, although the cost of the first pass is much larger since (i) the volume constraint is not applied at each design step, only when the difference between the actual and target volumes exceeds the 5% tolerance, and (ii) less elements and nodes need to be moved and migrated across processors. Meanwhile, the cost of both geometrically reinitializing the signed distance function level set and of optimizing the volume constraint offset is very affordable, as it represents less than 2% in total, with 4-5 dichotomy iterations needed to reach the desired accuracy. All three 3-D cases yield almost identical timing results, the only difference being in the cost of the volume constraint step, as the frequency at which consecutive corrections are applied depends on the geometrical specifics of the layout under consideration. This gives hope that the same conclusions may carry over to other multiple inlet/outlet duct flow problems of same dimensionality, tackled with comparable parameters. The close similarity carries over to the absolute run times per iteration shown in Fig. 20, meaning that the total run times in Tab. 5 are essentially driven by the number of design steps needed to converge (that in turn is driven by the
<figure><img src="image_20.png" /><figcaption>**Figure 19.** (a) Computational cost of the implemented algorithm, as obtained averaging 200 update steps of the 2-D single pipe problem with heated walls presented in Fig. 2. (b) Same as (a) for 100 update steps of the 3-D single pipe problem with one centered hot stripe presented in Fig. 12. (c-d) Same as (b) for the 3-D single pipe problem with (c) one upstream hot stripe and (d) two hot stripes presented in Fig. 15. All simulation parameters are those provided in Tabs. 2 and 4. The LS and LSF labels stand for level set (LS) and filtered level set (LSF), respectively.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Figure 20.** Average run time per iteration for the various 3-D cases documented in Fig. 19.</figcaption></figure>

number of steps needed to recover the proper volume of fluid). The reported run times, while large in a vacuum, are actually much lower than those that required to converge on a fixed uniform grid with similar mesh refinement. To give a taste, discretizing the single inlet/single outlet case with a uniform element size of $5 \times 10^{-3}$ would require about 140M elements, even though this would not suffice to match the interface value achieved herein. It is also worth emphasizing in this regards that we did not seek to optimize efficiency, neither by adjusting the initial design (we actually used numerous inclusions on purpose to showcase the ability of the method to support complex topological changes), nor by fine tuning the descent factor (the only requirement being that the displacement achieved at each step must be below the cut-off thickness of the level set for the evolved interface to remain accurately tracked).

## 7.2. *High-Reynolds-number flows*

We keep in mind that all cases reported herein are low-Reynolds-number, making the approach very relevant for a broad range of microfluidic applications, for instance the optimal design of microchannel heat sinks. Conversely, the significance of the computed optima as solutions to actual industrial heat exchanger problems is questionable, given the high Reynolds number at play. It lies out of the scope to tackle such cases right away, as the adjoint-based sensitivity analysis of high-Reynolds-number flow has its challenges, including the need to evaluate all sensitivities from time-dependent adjoint solutions to be solved repeatedly backwards in time, which is close to intractable in 3-D without sophisticated integration, interpolation and/or checkpointing schemes. Relevance is thus assessed a
<table><caption>**Table 5.** Run times for the various 3-D cases documented in Fig. 19.</caption><thead><tr><th></th><th></th><th></th><th></th></tr></thead><tbody><tr><td>64</td><td>»</td><td>»</td><td>Nb. cores</td></tr><tr><td>270h</td><td>340h</td><td>240h</td><td>Run time</td></tr><tr><td>800</td><td>1000</td><td>700</td><td>Nb. design steps</td></tr><tr><td>200</td><td>200</td><td>200</td><td>Nb. steps to target volume</td></tr></tbody></table>

<table><caption>**Table 6.** Performance of the optimal layouts under various high-Reynolds-number settings. All cost function values made non-dimensional using the inlet diameter and maximum inlet velocity (equivalently, using $ρu_i^3e_i^2$ as reference cost functional value).</caption><thead><tr><th></th><th></th><th></th><th></th></tr></thead><tbody><tr><td>1.74</td><td>1.62</td><td>1.59</td><td>Re = 12</td></tr><tr><td>2.23</td><td>1.94</td><td>1.89</td><td>Re = 1.2 × 10<sup>4</sup></td></tr></tbody></table>

posteriori by performing high-Reynolds-number simulations of several shapes generated over the course of optimization (the only requirement being that the volume of fluid constraint be satisfied). Results are reported here for the one single upstream hot stripe case presented in Fig. 15, for which 800 time-steps have been carried out at $Re = 1.2 \times 10^4$, after which the cost function has been averaged over an additional 200 time-steps. The results pertaining to the last three shapes in Fig. 16 are reported in Tab. 6, together with their baseline counterparts obtained at $Re = 12$. We also provide in Fig. 21 3-D instantaneous streamlines colored by the temperature, to emphasize that the flow undergoes a increasingly complex swirling motion after crossing the solide core, that keeps being accurately resolved by virtue of the anisotropic mesh adaptation procedure. While the cost function increases with $Re$ (which has been found to be because the optimal heat transfer performance decreased), it is interesting to note that the optimal shape keeps performing best even at such high Reynolds number, and that its superiority with respect to the other two shapes is more pronounced (the improvement in the cost function being respectively by 8.5% and 1.8% at $Re = 12$, but 15.2% and 2.5% at $Re = 1.2 \times 10^4$). This suggests that the present approach can deliver relevant robust and practical solutions to real-life heat exchanger applications.

## 8. Conclusion

The present study proves feasible to perform topology optimization of conjugate heat transfer systems using anisotropic meshes adapted under the constraint of a fixed number of nodes. The proposed approach combines a level set method to represent the boundary of the fluid domain by the zero iso-value of a signed distance function, and stabilized formulations of the state, adjoint, and level set transport equations cast in the Variational Multiscale (VMS) framework. The method has been shown to allow for drastic topology changes during the optimization process. Nonetheless, the main advantage over existing methods is the ability to capture all interfaces to a very high degree of accuracy using adapted meshes whose anisotropy matches that of the numerical solutions. The approach also considerably decreases the cost of improving the numerical precision, as the number of nodes needs be increased only in the anisotropy direction, hence only 2 times as many nodes are required to improve the resolution by a factor of 2, as opposed to 4 and 8 times in classical 2-D and 3-D isotropic calculations. This gives hope that the method can ease the transition
<figure><img src="image_22.png" /><figcaption>**Figure 21.** (a) Representative steady 3-D streamlines colored by the temperature computed at $Re = 12$ for the last three shapes in Fig. 16 by the magnitude of velocity. (b) Same as (a) for the instantaneous streamlines computed at $Re = 1.2 \times 10^4$.</figcaption></figure>

to manufacturable CAD models closely resembling the optimal topology, with important applicability to the engineering of low mass, high efficiency thermal devices, including heat exchangers, heat sinks or cold plates.

The method has been tested on several examples of bi-objective minimization weighing power dissipation and recoverable thermal power, including a large scale three-dimensional example involving several dozen million state degrees of freedom. The obtained optimal designs agree well with the existing literature, which assesses the relevance of the present implementation for engineering various thermal fluidic devices, and that even at non-small Reynolds numbers, as has been assessed a posteriori from dedicated simulations comparing the performance of various shapes generated over the optimization process. It is worth emphasizing in this regards that we did not seek to optimize efficiency, neither by adjusting the initial design, nor by fine tuning the descent factor (the only requirement being that the displacement achieved at each step must be below the cut-off thickness of the level set for the evolved interface to remain accurately tracked). Future work should aim at improving the numerical efficiency, for instance by adding nucleation mechanisms to alleviate the need for initial design with holes, as this makes it difficult to fulfill the proper volume constraint from the outset, and requires substantial mesh refinement to avoid clogging the fluid path in the early stage of optimization. Other research directions include application to complex physics more representative of real life situations (e.g., multiphase flows, fluid-structure interactions), as well as assessment of multi-component adaptation criteria taking into account the difference in the spatial supports of the state and adjoint solutions, to further improve the accuracy of the gradient evaluations.

**Author Contributions:** Conceptualization, E.H.; methodology, E.H. and P.M.; software and validation, E.H., P.M. and W.A.N.; formal analysis, E.H. and P.M.; investigation, W.A.N.; resources, E.H.; data curation, W.A.N.; writing—original draft preparation, P.M.; writing—review and editing, E.H., P.M., D.L., D.S. and W.A.N.; visualization, E.H., P.M. and W.A.N.; supervision, E.H., P.M. and D.S.; project administration, D.L., E.H. and D.S.; funding acquisition, D.L. and D.S.. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work is part of the PANTTHER project, that has received funding from the Clean Sky2 Joint Undertaking (JU) under grant agreement No 886698. The JU receives support from the European Union's Horizon 2020 research and innovation program and the Clean Sky 2 JU members
other than the Union. It reflects only the authors' view and the JU is not responsible for any use that may be made of the information it contains.

**Institutional Review Board Statement:** Not applicable

**Informed Consent Statement:** Not applicable

**Data Availability Statement:** The data presented in this study are available on reasonable request from the corresponding author. The data are not publicly available due to confidential reasons.

**Conflicts of Interest:** The authors declare no conflict of interest.

# Appendix I Adjoint-based sensitivity analysis

## Appendix I.1 General setting

This first appendix presents the main steps of the adjoint method used to derive the sensitivity of the cost function sensitivity to a design variable $\beta$ representing deformations of the interface under the form of local surface normal displacements. We recall that the fluid is characterized by its density $\rho_f$, dynamic viscosity $\mu_f$, thermal heat capacity $c_{pf}$ and thermal conductivity $k_f$. Conversely, the solid is only characterized by a thermal conductivity $k_s$. The flow motion in the fluid domain $\Omega_f$ is modeled after the coupled steady incompressible Navier-Stokes and heat equations, hence

$$ \nabla \cdot \mathbf{u} = 0 \qquad \text{in } \Omega_f, \qquad (A34) $$

$$ \rho_f \mathbf{u} \cdot \nabla \mathbf{u} = - \nabla p + \nabla \cdot (2\mu_f \varepsilon(\mathbf{u})) \qquad \text{in } \Omega_f, \qquad (A35) $$

$$ \rho_f c_{pf} \mathbf{u} \cdot \nabla T_f = \nabla \cdot (k_f \nabla T_f) \qquad \text{in } \Omega_f, \qquad (A36) $$

with open flow boundary conditions

$$ \mathbf{u} = \mathbf{u}_i \quad \text{on } \Gamma_i, \qquad (A37) $$

$$ T_f = T_i \quad \text{on } \Gamma_i, \qquad (A38) $$

$$ \mathbf{u} = \mathbf{u}_o \quad \text{on } \Gamma_o, \qquad (A39) $$

$$ \nabla T_f \cdot \mathbf{n} = 0 \quad \text{on } \Gamma_o, \qquad (A40) $$

consisting of a prescribed velocity and temperature at the inlet, and a prescribed velocity and zero heat flux at the outlet (with outflow velocity adjusted for the total amount of mass flow exiting through the outlet to match exactly that entering through the inlet). Finally, we assume zero velocity at the interface, together with a (Robin) convective heat flux condition, i.e.,

$$ \mathbf{u} = 0 \quad \text{on } \Gamma, \qquad (A41) $$

$$ k_f \nabla T_f \cdot \mathbf{n} = \eta (T_f - T^*) \quad \text{on } \Gamma, \qquad (A42) $$

where $T^*$ is a reference temperature and $\eta$ is a heat transfer coefficient driving the budget of heat-flux and temperature variance between the fluid and solid domains (as it is widely recognized that neither isothermal nor isoflux boundary conditions can realistically mimic actual heat transfer in practical applications, especially when the thermal diffusivity of the solid and the fluid are of the same order of magnitude [61]).

The problem of minimizing the cost function subject to coupled Navier-Stokes and heat equations as state equations is tackled using the continuous adjoint method [37–39]. One first forms the Lagrangian

$$ \mathcal{L} = \int_{\Gamma_i \cup \Gamma_o} J \, ds - \int_{\Omega_f} \tilde{p} \nabla \cdot \mathbf{u} \, dv - \int_{\Omega_f} \tilde{\mathbf{u}} \cdot (\rho_f \mathbf{u} \cdot \nabla \mathbf{u} + \nabla p - \nabla \cdot (2\mu_f \varepsilon(\mathbf{u})) \, dv \\ - \int_{\Omega_f} \tilde{T}_f \cdot (\rho_f c_{pf} \mathbf{u} \cdot \nabla T_f - \nabla \cdot (k_f \nabla T_f)) \, dv, \qquad (A43) $$
featuring the adjoint pressure $\tilde{p}$ as the Lagrange multiplier for the continuity equation (A34), the adjoint velocity $\tilde{u}$ as the Lagrange multiplier for the momentum equations (A35), and the adjoint temperature $\tilde{T}_f$ as the Lagrange multiplier for the heat equation (A36). One then seeks to decompose the variation of $\mathcal{L}$ due to a change in the interface position into individual variations with respect to the adjoint, state and design variables. The variation with respect to the adjoint variables

$$ \delta(\tilde{u}, \tilde{p}, \tilde{T}_f)\mathcal{L} = - \int_{\Omega_f} \delta\tilde{p} \nabla \cdot \mathbf{u} \, dv - \int_{\Omega_f} \delta\tilde{u} \cdot (\rho_f \mathbf{u} \cdot \nabla \mathbf{u} + \nabla p - \nabla \cdot (2\mu_f \varepsilon(\mathbf{u}))) \, dv \\ - \int_{\Omega_f} \delta\tilde{T}_f \cdot (\rho_f c_{p_f} \mathbf{u} \cdot \nabla T_f - \nabla \cdot (k_f \nabla T_f)) \, dv, \quad (A44) $$

is trivially zero as long as $(\mathbf{u}, p, T_f)$ is solution to the above coupled Navier-Stokes and heat equations, in which case $\mathcal{L} = J$. After integrating by parts, the variation with respect to the state variables is

$$ \begin{aligned} \delta(\mathbf{u}, p, T_f)\mathcal{L} = & \int_{\Omega_f} (\nabla \cdot \tilde{\mathbf{u}})\delta p \, dv \\ & + \int_{\Omega_f} (-\rho_f \mathbf{u} \cdot \nabla \tilde{\mathbf{u}} + \rho_f \nabla \mathbf{u}^T \cdot \tilde{\mathbf{u}} - \nabla \tilde{p} - \nabla \cdot (2\mu_f \varepsilon(\tilde{\mathbf{u}})) - \rho_f c_{p_f} T_f \nabla \tilde{T}_f) \cdot \delta \mathbf{u} \, dv \\ & + \int_{\Omega_f} (\rho_f c_{p_f} \mathbf{u} \cdot \nabla \tilde{T}_f + \nabla \cdot (k_f \nabla \tilde{T}_f)) \delta T_f \, dv \\ & + \int_{\Gamma_i \cup \Gamma_o} \partial_{\mathbf{u}} J \cdot \delta \mathbf{u} \, ds + \int_{\partial \Omega_f} (\tilde{p} \mathbf{n} + 2\mu_f \varepsilon(\tilde{\mathbf{u}}) \cdot \mathbf{n} + \rho_f (\mathbf{u} \cdot \mathbf{n}) \tilde{\mathbf{u}} + \rho_f c_{p_f} T_f \tilde{T}_f \mathbf{n}) \cdot \delta \mathbf{u} \, ds \\ & - \int_{\Gamma_i \cup \Gamma_o} \partial_p J_{\Gamma} \mathbf{n} \cdot (-\delta p \mathbf{n} + 2\mu_f \varepsilon(\delta \mathbf{u}) \cdot \mathbf{n}) \, ds - \int_{\partial \Omega_f} \tilde{\mathbf{u}} \cdot (-\delta p \mathbf{n} + 2\mu_f \varepsilon(\delta \mathbf{u}) \cdot \mathbf{n}) \, ds \\ & + \int_{\Gamma_i \cup \Gamma_o} \partial_{T_f} J \delta T_f \, ds + \int_{\partial \Omega_f} ((k_f \nabla \tilde{T}_f \cdot \mathbf{n} + \rho_f c_{p_f} (\mathbf{u} \cdot \mathbf{n}) \tilde{T}_f) \delta T_f - k_f \tilde{T}_f \nabla \delta T_f \cdot \mathbf{n}) \, ds, \end{aligned} \quad (A45) $$

on behalf of the viscous stress being purely tangential in incompressible flows. At this stage, adjoint equations and boundary conditions are designed to ensure $\delta(\mathbf{u}, p, T_f)\mathcal{L} = 0$, which requires the domain and boundary integrals to vanish individually in (A45). Keeping in mind that we work here under the assumption of a fixed interface (since the design variable is constant) and that the cost function does not depend on the quantities on the wall, we obtain the linear, homogeneous problem

$$ \begin{aligned} & \nabla \cdot \tilde{\mathbf{u}} = 0 && \text{in } \Omega_f, \quad (A46) \\ & -\rho_f \mathbf{u} \cdot \nabla \tilde{\mathbf{u}} + \rho_f \nabla \mathbf{u}^T \cdot \tilde{\mathbf{u}} = \nabla \tilde{p} + \nabla \cdot (2\mu_f \varepsilon(\tilde{\mathbf{u}})) + \rho c_{p_f} T_f \nabla \tilde{T}_f && \text{in } \Omega_f, \quad (A47) \\ & -c_{p_f} \mathbf{u} \cdot \nabla \tilde{T}_f = \nabla \cdot (k_f \nabla \tilde{T}_f) && \text{in } \Omega_f, \quad (A48) \end{aligned} $$

driven by the non-homogeneous boundary conditions

$$ \begin{aligned} \tilde{\mathbf{u}} &= - \partial_p J_{\Gamma \mathbf{n}} && \text{on } \Gamma_i, \quad (A49) \\ \tilde{T}_f &= 0 && \text{on } \Gamma_i, \quad (A50) \\ \tilde{\mathbf{u}} &= - \partial_p J_{\Gamma \mathbf{n}} && \text{on } \Gamma_o, \quad (A51) \\ k_f \nabla \tilde{T}_f \cdot \mathbf{n} + \rho_f c_{p_f} (\mathbf{u} \cdot \mathbf{n}) \tilde{T}_f &= - \partial_{T_f} J_{\Gamma} && \text{on } \Gamma_o, \quad (A52) \\ \tilde{\mathbf{u}} &= 0 && \text{on } \Gamma, \quad (A53) \\ k_f \nabla \tilde{T}_f \cdot \mathbf{n} &= \eta \tilde{T}_f && \text{on } \Gamma, \quad (A54) \end{aligned} $$

associated to (A37)-(A42), including an interface adjoint convective heat flux condition using the same heat transfer coefficient $\eta$. The key difference between the state and adjoint equations lies in the minus sign in front of the convective term of Eqs. (A47)-(A48), to reflect
that adjoint information is convected upstream, not downstream, due to the non-normality of the linearized evolution operator [62]. Expressing the interface perturbations after [63] as

954

955

$$ \delta \mathbf{u} = \beta \nabla \mathbf{u} \cdot \mathbf{n}, \quad \delta T_f = \beta \nabla T_f \cdot \mathbf{n}, \qquad (A55) $$

the variation with respect to the design variable (now encompassing the domain deformation) is ultimately computed as

956

957

$$ \delta_{\beta} J_s \equiv \delta_{\beta} \mathcal{L} = \beta \int_{\Gamma} (\tilde{\mathbf{p}}\mathbf{n} + 2\mu_f \varepsilon(\tilde{\mathbf{u}}) \cdot \mathbf{n} + \rho_f c_p f_f \tilde{T}_f \mathbf{n}) \cdot (\nabla \mathbf{u} \cdot \mathbf{n}) \, ds \\ + \beta \int_{\Gamma} k_f (\nabla \tilde{T}_f \cdot \mathbf{n})(\nabla T_f \cdot \mathbf{n}) \, ds - \beta \int_{\Gamma} \tilde{T}_f (k_f \nabla (\nabla T_f \cdot \mathbf{n}) \cdot \mathbf{n}) \, ds. \qquad (A56) $$

which reduces to

958

$$ \delta_{\beta} J_s = \beta \int_{\Gamma} \mu_f (\nabla \tilde{\mathbf{u}} \cdot \mathbf{n}) \cdot (\nabla \mathbf{u} \cdot \mathbf{n}) \, ds \\ + \beta \int_{\Gamma} k_f (\nabla \tilde{T}_f \cdot \mathbf{n})(\nabla T_f \cdot \mathbf{n}) \, ds - \beta \int_{\Gamma} \tilde{T}_f (k_f \nabla (\nabla T_f \cdot \mathbf{n}) \cdot \mathbf{n}) \, ds, \qquad (A57) $$

due to the incompressibility of the state and adjoint solutions [37]. The simplest steepest-descent algorithm implemented herein therefore moves down the cost function, in the direction of the steepest slope using

959

960

961

$$ \beta = -\mu_f (\nabla \tilde{\mathbf{u}} \cdot \mathbf{n}) \cdot (\nabla \mathbf{u} \cdot \mathbf{n}) - k_f (\nabla \tilde{T}_f \cdot \mathbf{n})(\nabla T_f \cdot \mathbf{n}) + \tilde{T} (k_f \nabla (\nabla T_f \cdot \mathbf{n}) \cdot \mathbf{n}), \qquad (A58) $$

up to a positive multiplicative factor to control the step taken in the gradient direction.

962

## Appendix I.2 Extension to the IVM-VMS resolution framework

963

In practice, the heat transfer coefficient η ensuring that the fluid and solid exchange the proper amount of heat remains an unknown. Computing said coefficient is no small task, as it often requires solving an inverse problem to assimilate relevant experimental data, which in turn requires such data to be available. Such a lack of availability is generally acknowledged to be a limiting issue for practical applications, especially for topology optimization where varying the shape, amount, and distribution of the solid domain is integral to the optimization process itself. The immersed volume method underlying this research combines both the fluid and solid phases into a single fluid with variable material properties. It thus solves Navier-Stokes and heat equations

964

965

966

967

968

969

970

971

972

$$ \nabla \cdot \mathbf{u} = 0 \qquad \text{in } \Omega, \qquad (A59) $$

$$ \rho \mathbf{u} \cdot \nabla \mathbf{u} = - \nabla p + \nabla \cdot (2\mu \varepsilon(\mathbf{u})) \qquad \text{in } \Omega, \qquad (A60) $$

$$ \rho c_p \mathbf{u} \cdot \nabla T = \nabla \cdot (k \nabla T) \qquad \text{in } \Omega, \qquad (A61) $$

identical to (A34)-(A36), but with variable density, viscosity, thermal heat capacity and thermal conductivity, adequately interpolated over a small layer around the interface, and otherwise equal to their fluid and solid values. This allows dropping altogether the interface thermal condition (and thus alleviates the need for a heat transfer coefficient) because the

973

974

975

976
amount of heat exchanged at the interface is entirely determined by the individual material properties heat on either side of it, hence the associated boundary conditions

$$ \mathbf{u} = \mathbf{u}_i \quad \text{on } \Gamma_i, \qquad (\text{A62}) $$

$$ T = T_i \quad \text{on } \Gamma_i, \qquad (\text{A63}) $$

$$ \mathbf{u} = \mathbf{u}_o \quad \text{on } \Gamma_o, \qquad (\text{A64}) $$

$$ \nabla T \cdot \mathbf{n} = 0 \quad \text{on } \Gamma_o, \qquad (\text{A65}) $$

$$ \mathbf{u} = \mathbf{0} \quad \text{on } \Gamma_w^i \cup \Gamma_w^a, \qquad (\text{A66}) $$

$$ T = T_w \quad \text{on } \Gamma_w^i, \qquad (\text{A67}) $$

$$ \nabla T \cdot \mathbf{n} = 0 \quad \text{on } \Gamma_w^a. \qquad (\text{A68}) $$

Provided the velocity is zero in the solid domain (either because a very high value of the solid-to-fluid viscosity ratio is used, or, as is the case here, because the constraint is hard-coded; see Sec. 2) and the no-slip interface condition is satisfied, the convective term vanishes in (A61), that reduces to the pure conduction equation for the solid, together with prescribed temperature and zero heat flux conditions at the solid isothermal and adiabatic walls, respectively.

The exact same approach is applied to the adjoint equations, by solving adjoint Navier-Stokes and heat equations

$$ \nabla \cdot \tilde{\mathbf{u}} = 0 \qquad \text{in } \Omega, \qquad (\text{A69}) $$

$$ -\rho \mathbf{u} \cdot \nabla \tilde{\mathbf{u}} + \rho \nabla \mathbf{u}^T \cdot \tilde{\mathbf{u}} = \nabla \tilde{p} + \nabla \cdot (2\mu\varepsilon(\tilde{\mathbf{u}})) + \rho c_p T \nabla \tilde{T} \quad \text{in } \Omega, \qquad (\text{A70}) $$

$$ -c_p \mathbf{u} \cdot \nabla \tilde{T} = \nabla \cdot (k \nabla \tilde{T}) \qquad \text{in } \Omega, \qquad (\text{A71}) $$

identical to (A46)-(A48), but with variable density, viscosity, thermal heat capacity and thermal conductivity, together with boundary conditions

$$ \tilde{\mathbf{u}} = - \partial_p J_{\Gamma \mathbf{n}} \quad \text{on } \Gamma_i, \qquad (\text{A72}) $$

$$ \tilde{T} = 0 \quad \text{on } \Gamma_i, \qquad (\text{A73}) $$

$$ \tilde{\mathbf{u}} = - \partial_p J_{\Gamma \mathbf{n}} \quad \text{on } \Gamma_o, \qquad (\text{A74}) $$

$$ k \nabla \tilde{T} \cdot \mathbf{n} + \rho c_p (\mathbf{u} \cdot \mathbf{n}) \tilde{T} = - \partial_T J_{\Gamma} \quad \text{on } \Gamma_o, \qquad (\text{A75}) $$

$$ \tilde{\mathbf{u}} = \mathbf{0} \quad \text{on } \Gamma_w^i \cup \Gamma_w^a, \qquad (\text{A76}) $$

$$ \tilde{T} = 0 \quad \text{on } \Gamma_w^i, \qquad (\text{A77}) $$

$$ \nabla \tilde{T} \cdot \mathbf{n} = 0 \quad \text{on } \Gamma_w^a. \qquad (\text{A78}) $$

Assuming the adjoint velocity is zero in the solid domain (again this is hard-coded in the present context of extremely stretched, anisotropic mesh elements), Eq. (A71) reduces to the pure conduction equation for the solid with zero adjoint temperature at the solid isothermal walls and zero adjoint heat flux at the solid adiabatic walls, just as what would be obtained adding the solid conduction equation to the Lagrangian (A43) and evaluating the variation with respect to the state variables. Finally, we compute the steepest-descent displacement as

$$ \beta = -\mu(\nabla \tilde{\mathbf{u}} \cdot \mathbf{n}) \cdot (\nabla \mathbf{u} \cdot \mathbf{n}) - k(\nabla \tilde{T} \cdot \mathbf{n})(\nabla T \cdot \mathbf{n}) + \tilde{T}(k\nabla(\nabla T \cdot \mathbf{n}) \cdot \mathbf{n}). \qquad (\text{A79}) $$

Since $\nabla(\mathbf{n} \cdot \mathbf{n}) = 0$ due to the normal vector having unit norm, it can be shown that

$$ \nabla(\nabla T \cdot \mathbf{n}) \cdot \mathbf{n} = (\nabla(\nabla T) \cdot \mathbf{n}) \cdot \mathbf{n} + (\nabla \wedge \mathbf{n}) \cdot (\mathbf{n} \wedge \nabla T), \qquad (\text{A80}) $$

and thus

$$ \nabla(\nabla T \cdot \mathbf{n}) \cdot \mathbf{n} = (\nabla \wedge \mathbf{n}) \cdot (\mathbf{n} \wedge \nabla T), \qquad (\text{A81}) $$
because the second derivatives vanish due to the use of P1 linear finite element approximations. Also, since the normal vector in a level set framework is computed as $\mathbf{n} = \nabla\phi/\|\nabla\phi\|$, we have

$$ \nabla \left( \frac{\phi}{\|\nabla\phi\|} \right) = \mathbf{n} - \frac{\phi}{\|\nabla\phi\|^3} H(\phi) \cdot \nabla\phi, \qquad (A82) $$

where $H$ is the Hessian of $\phi$. Since $\phi = 0$ by design on the interface, we thus have

$$ \nabla \wedge \mathbf{n} = \nabla \wedge \nabla \left( \frac{\phi}{\|\nabla\phi\|} \right) = \mathbf{0}. \qquad (A83) $$

It follows that

$$ \nabla(\nabla T \cdot \mathbf{n}) \cdot \mathbf{n} = 0, \qquad (A84) $$

and Eq. (A79) reduces to

$$ \beta = -\mu(\nabla \tilde{\mathbf{u}} \cdot \mathbf{n}) \cdot (\nabla \mathbf{u} \cdot \mathbf{n}) - k(\nabla \tilde{T} \cdot \mathbf{n})(\nabla T \cdot \mathbf{n}). \qquad (A85) $$

The last step in the process is to specify the state derivatives of the cost function, that, for the present linear weighted sum of dissipated power (hydraulic component) and recoverable thermal power (thermal component)

$$ J = (1 - \omega)(p + \frac{1}{2}\rho(\mathbf{u} \cdot \mathbf{u}))(\mathbf{u} \cdot \mathbf{n}) - \omega\rho c_p T (\mathbf{u} \cdot \mathbf{n}), \qquad (A86) $$

are given by

$$ \partial_p J = (1 - \omega)(\mathbf{u} \cdot \mathbf{n}), \quad \partial_{\mathbf{u}} J = (1 - \omega)(p_{tot} \mathbf{n} + \rho(\mathbf{u} \cdot \mathbf{n})\mathbf{u}) - \omega\rho c_p T \mathbf{n}, \quad \partial_T J = -\omega\rho c_p (\mathbf{u} \cdot \mathbf{n}). \qquad (A87) $$

## Appendix J Coarse scale VMS variational problems

This second appendix is devoted to the stabilized finite element numerical framework used to compute all solutions of interest on anisotropic adapted meshes and to perform the design update steps. For the sake of simplicity in the notations (and as long as it does not lead to ambiguity), we omit in what follows the distinction between all continuous variables (e.g., domains, solutions, operators) and their discrete finite element counterparts, as well as the dependency of all variables on the iteration of the optimization process. In practice, the state equations are solved sequentially, i.e., we solve first the Navier–Stokes equations, then use the resulting velocity to solve the heat equation. Due to the reversal in space-time directionality (and thus in causality), the adjoint equations are also solved sequentially but in reverse order, i.e., we solve first the adjoint heat equation, then use the resulting adjoint temperature to solve the adjoint Navier–Stokes equations.

### Appendix J.1 Navier-Stokes equations

In practice, the state solution is computed by time-stepping the unsteady Navier–Stokes equations with large time steps to accelerate convergence towards a steady state (the stopping criterion being here for two consecutive time steps to differ by less than $10^{-6}$ in $L^\infty$ norm). In order to deal with the time-dependency and non-linearity of the momentum equation, the transport time of the time scale is assumed much smaller than that of the coarse scale. In return, the fine scale contribution to the transport velocity is neglected, and the fine scale is not tracked in time (although it is driven by the coarse-scale, time-dependent residuals and therefore does vary in time in a quasi-static manner). In-depth technical and mathematical details together with extensive discussions regarding the relevance of the
approximations can be found in [64]. Ultimately, the coarse scale variational problem to be solved is formulated as

$$ \int_{\Omega} (\rho \partial_t \mathbf{u} + \rho \mathbf{u} \cdot \nabla \mathbf{u}) \cdot \mathbf{w} \, dv + \int_{\Omega} 2\mu \varepsilon(\mathbf{u}) : \varepsilon(\mathbf{w}) \, dv - \int_{\Omega} p(\nabla \cdot \mathbf{w}) \, dv + \int_{\Omega} (\nabla \cdot \mathbf{u}) q \, dv \\ - \sum_{k=1}^{N_e} \int_{\Omega_k} \tau_1 \mathbf{r}_1 \cdot (\rho \mathbf{u} \cdot \nabla \mathbf{w}) \, dv - \sum_{k=1}^{N_e} \int_{\Omega_k} \tau_1 \mathbf{r}_1 \cdot \nabla q \, dv - \sum_{k=1}^{N_e} \int_{\Omega_k} \tau_2 r_2 (\nabla \cdot \mathbf{w}) \, dv = 0, \quad (A88) $$

where we have considered a discretization of $\Omega$ into $N_e$ non-overlapping elements (triangles or tetrahedrons), $\Omega_k$ is the domain occupied by the $k$th element, and $\mathbf{r}_1$ and $r_2$ are the momentum and continuity residuals

$$ -\mathbf{r}_1 = \rho \partial_t \mathbf{u} + \rho \mathbf{u} \cdot \nabla \mathbf{u} + \nabla p, \quad -r_2 = \nabla \cdot \mathbf{u}, \quad (A89) $$

whose second derivatives vanish since we use linear interpolation functions. Finally, $\tau_1$ and $\tau_2$ are the ad-hoc stabilization coefficients defined in [43,65], computed on each element from a characteristic size $h$ (here, the element diameter in the direction of the velocity, to support using anisotropic meshes with highly stretched elements [66]), and velocity $u$ (here, the average $L^2$ norm of the nodal element velocities). Equation (A88) is discretized with a first-order-accurate time-integration scheme combining semi-implicit treatment of the convection term, implicit treatment of the viscous, pressure and divergence terms, and explicit treatment of the stabilization coefficients.

## Appendix J.2 Heat equation

The coarse scale variational problem for the heat equation reads

$$ \int_{\Omega} (\rho c_p \mathbf{u} \cdot \nabla T) s \, dv + \int_{\Omega} k \nabla T \cdot \nabla s \, dv - \sum_{k=1}^{N_e} \int_{\Omega_k} \tau_3 r_3 \rho c_p \mathbf{u} \cdot \nabla s \, dv - \sum_{k=1}^{N_e} \int_{\Omega_k} \tau_4 r_4 \rho c_p \mathbf{u}_{\parallel} \cdot \nabla s \, dv = 0, \quad (A90) $$

where $\mathbf{u}_{\parallel}$ is the (normalized) velocity projected along the direction of the temperature gradient defined as

$$ \mathbf{u}_{\parallel} = \frac{\mathbf{u} \cdot \nabla T}{\|\nabla T\|^2} \nabla T, \quad (A91) $$

$r_3$ is the heat equation residual

$$ -r_3 = \rho c_p \mathbf{u} \cdot \nabla T, \quad (A92) $$

and $\tau_{3,4}$ are mesh-dependent stabilization parameters acting both in the direction of the solution and of its gradient, that proceed from the stabilization of the ubiquitous convection-diffusion-reaction equation [67,68], whose definition is given in [69,70]. Equation (A90) is solved with implicit treatment of the convection term and conduction terms (as the convection velocity is taken as a given) and explicit treatment of the stabilization coefficients.

## Appendix J.3 Adjoint heat equation

Since the adjoint heat equation (8) is formally identical to its state counterpart (save for the change in the sign of the convection velocity), its coarse scale variational problem deduces straightforwardly as

$$ \begin{aligned} & -\int_{\Omega} (\rho c_p \mathbf{u} \cdot \nabla T) s \, dv + \int_{\Omega} k \nabla T \cdot \nabla s \, dv + \sum_{k=1}^{N_e} \int_{\Omega_k} \tau_3 r_3 \rho c_p \mathbf{u} \cdot \nabla s \, dv + \sum_{k=1}^{N_e} \int_{\Omega_k} \tau_4 r_4 \rho c_p \mathbf{u}_{\parallel} \cdot \nabla s \, dv \\ & - \int_{\Gamma_o} \rho (\mathbf{u} \cdot \mathbf{n}) \tilde{T} s \, ds = \int_{\Gamma_o} \partial_T J s \, ds, \end{aligned} \quad (A93) $$
and features the same residual $r_3$ and stabilization coefficients $\tau_{3,4}$ as in J.2. Note, Eq. (A93) also includes boundary terms evaluated at the outlet, which is because the integration by part of the conductive term unveils a boundary term

$$ \int_{\partial\Omega} (k\nabla T \cdot \mathbf{n}) s \, ds = -\int_{\Gamma_o} (\rho(\mathbf{u} \cdot \mathbf{n})\tilde{T} + \partial_T J)s \, ds, \qquad (A94) $$

due to the adjoint outflow thermal condition (A52). In practice, though, preliminary tests have assessed that the adjoint thermal power term $\rho c_p(\mathbf{u} \cdot \mathbf{n})\tilde{T}$ consistently dominates by at least four orders of magnitude over the heat flux term $k\nabla\tilde{T} \cdot \mathbf{n}$. We thus end up simplifying the numerical implementation using the approximate condition

$$ \tilde{T} = -\omega \quad \text{on } \Gamma_o, \qquad (A95) $$

and solve the scale variational problem (without boundary term)

$$ -\int_{\Omega} (\rho c_p \mathbf{u} \cdot \nabla T) s \, dv + \int_{\Omega} k \nabla T \cdot \nabla s \, dv + \sum_{k=1}^{N_e} \int_{\Omega_k} \tau_3 r_3 \rho c_p \mathbf{u} \cdot \nabla s \, dv + \sum_{k=1}^{N_e} \int_{\Omega_k} \tau_4 r_4 \rho c_p \mathbf{u}_\| \cdot \nabla s \, dv = 0, \qquad (A96) $$

with implicit treatment of the convection term and conduction terms (as the convection velocity is taken as a given) and explicit treatment of the stabilization coefficients.

## Appendix J.4 Adjoint Navier-Stokes equations

Application of the stabilized formulation, as described above, to the adjoint Navier-Stokes equations yields the following coarse scale variational problem

$$ \int_{\Omega} (-\rho \mathbf{u} \cdot \nabla \tilde{\mathbf{u}} + \rho \nabla \mathbf{u}^T \cdot \tilde{\mathbf{u}}) \cdot \mathbf{w} \, dv + \int_{\Omega} 2\mu \varepsilon(\tilde{\mathbf{u}}) : \varepsilon(\mathbf{w}) \, dv + \int_{\Omega} \tilde{p}(\nabla \cdot \mathbf{w}) \, dv - \int_{\Omega} \rho c_p T \nabla \tilde{T} \cdot \mathbf{w} \, dv + \int_{\Omega} (\nabla \cdot \tilde{\mathbf{u}}) q \, dv - \sum_{k=1}^{N_e} \int_{\Omega_k} \tau_1 \tilde{\mathbf{r}}_1 \cdot (-\rho \mathbf{u} \cdot \nabla \mathbf{w}) \, dv - \sum_{k=1}^{N_e} \int_{\Omega_k} \tau_1 \tilde{\mathbf{r}}_1 \cdot \nabla q \, dv - \sum_{k=1}^{N_e} \int_{\Omega_k} \tau_2 \tilde{\mathbf{r}}_2 (\nabla \cdot \mathbf{w}) \, dv = 0. \qquad (A97) $$

The associated momentum and continuity residuals read

$$ -\tilde{\mathbf{r}}_1 = -\rho \mathbf{u} \cdot \nabla \tilde{\mathbf{u}} + \rho \nabla \mathbf{u}^T \cdot \tilde{\mathbf{u}} - \nabla \tilde{p}, \quad -\tilde{\mathbf{r}}_2 = \nabla \cdot \tilde{\mathbf{u}}, \qquad (A98) $$

and the stabilization coefficients $\tau_{1,2}$ are the same as those in J.1. This implicitly amounts to neglecting the additional stabilization stemming from the $\rho\nabla\mathbf{u}^T \cdot \tilde{\mathbf{u}}$ term describing the production of adjoint perturbations, that has been found to have no effect on the numerical results, as the problems considered herein are in the convection (not reaction) dominated limit.

## Appendix J.5 Interface update scheme using the convective level set method

The auto-reinitialization level set problem (14) is solved with an SUPG method, whose stabilization proceeds from that of the ubiquitous convection-diffusion-reaction equation [67,68]. The associated variational problem is formulated as

$$ \int_{\Omega} (\partial_{\tau} \phi + \mathbf{a}_{\tau} \cdot \nabla \phi) \xi \, dv - \int_{\Omega_k} \tau_5 r_5 \mathbf{a}_{\tau} \cdot \nabla \xi \, dv = \int_{\Omega} S \xi \, dv, \qquad (A99) $$

with residual

$$ -r_5 = \partial_{\tau} \phi + \mathbf{a}_{\tau} \cdot \nabla \phi - S, \qquad (A100) $$
and stabilization coefficient defined in the convection-dominated limit in [65]. It is easily checked that all terms scale as $1/\Delta\tau$, so we can set $\Delta\tau = 1$ without any loss of generality because the solution is ultimately independent on the pseudo-time step value. Equation (A99) is solved with semi-implicit treatment of the convection term (as the convection velocity $a_\tau$ depends on main unknown $\phi$) and explicit treatment of the source term and stabilization coefficients.

References

1. Bendsøe, M.P.; Kikuchi, N. Generating optimal topologies in structural design using a homogenization method. *Comput. Methods Appl. Mech. Engrg.* **1988**, *71*, 197-224.

2. Bendsøe, M.P.; Sigmund, O. *Topology optimization: theory, methods, and applications; Springer Science & Business Media*, 2003.

3. Yang, R.J.; Chahande, A.I. Automotive applications of topology optimization. *Struct. Opt.* **1995**, *9*, 245-249.

4. Zhu, J.H.; Zhang, W.H.; Xia, L. Topology optimization in aircraft and aerospace structures design. *Arch. Comput. Method. Engrg.* **2016**, *23*, 595-622.

5. Sigmund, O.; Maute, K. Topology optimization approaches. *Struct. Multidiscipl. Optim.* **2013**, *48*, 1031-1055.

6. Deaton, J.D.; Grandhi, R.V. A survey of structural and multidisciplinary continuum topology optimization: post 2000. *Struct. Multidiscipl. Optim.* **2014**, *49*, 1-38.

7. Suzuki, K.; Kikuchi, N. A homogenization method for shape and topology optimization. *Comput. Methods Appl. Mech. Engrg.* **1991**, *93*, 291-318.

8. Allaire, G.; Bonnetier, E.; Francfort, G.; Jouve, F. Shape optimization by the homogenization method. *Numer. Math.* **1997**, *76*, 27-68.

9. Sethian, J.A.; Wiegmann, A. Structural boundary design via level set and immersed interface methods. *J. Comput. Phys.* **2000**, *163*, 489-528.

10. Wang, M.Y.; Wang, X.; Guo, D. A level set method for structural topology optimization. *Comput. Methods Appl. Mech. Engrg.* **2003**, *192*, 227-246.

11. Allaire, G.; Jouve, F.; Toader, A.M. Structural optimization using sensitivity analysis and a level-set method. *J. Comput. Phys.* **2004**, *194*, 363-393.

12. Van Dijk, N.P.; Maute, K.; Langelaar, M.; Van Keulen, F. Level-set methods for structural topology optimization: a review. *Struct. Multidiscipl. Optim.* **2013**, *48*, 437-472.

13. Dbouk, T. A review about the engineering design of optimal heat transfer systems using topology optimization. *Appl. Therm. Eng.* **2017**, *112*, 841-854.

14. Alexandersen, J.; Andreasen, C.S. A review of topology optimization for fluid-based problems. *Fluids* **2020**, *5*, 29.

15. Yaji, K.; Yamada, T.; Kubo, S.; Izui, K.; Nishiwaki, S. A topology optimization method for a coupled thermal-fluid problem using level set boundary expressions. *Int. J. Heat Mass Transfer* **2015**, *81*, 878-888.

16. Coffin, P.; Maute, K. Level set topology optimization of cooling and heating devices using a simplified convection model. *Struct. Multidiscipl. Optim.* **2016**, *53*, 985-1003.

17. Feppon, F.; Allaire, G.; Dapogny, C.; Jolivet, P. Topology optimization of thermal fluid-structure systems using body-fitted meshes and parallel computing. *J. Comput. Phys.* **2020**, *417*, 1095-1141.

18. Dede, E.M. Multiphysics topology optimization of heat transfer and fluid flow systems. In Proceedings of the Procs. of the COMSOL Users Conference, 2009.

19. Yoon, G.H. Topological design of heat dissipating structure with forced convective heat transfer. *J. Mech. Sci. Technol.* **2010**, *24*, 1225-1233.

20. Marck, G. Optimisation topologique des transferts de chaleur et de masse: application aux échangeurs de chaleur. PhD thesis, Ecole Nationale Supérieure des Mines de Paris, 2012.

21. Kontoleontos, E.A.; Papoutsis-Kiachagias, E.M.; Zymaris, A.S.; Papadimitriou, D.I.; Giannakoglou, K.C. Adjoint-based constrained topology optimization for viscous flows, including heat transfer. *Eng. Optim.* **2013**, *45*, 941-961.

22. Matsumori, T.; Kondoh, T.; Kawamoto, A.; Nomura, T. Topology optimization for fluid-thermal interaction problems under constant input power. *Struct. Multidiscipl. Optim.* **2013**, *47*, 571-581.

23. Alexandersen, J.; Aage, N.; Andreasen, C.; Sigmund, O. Topology optimisation for natural convection problems. *Int. J. Numer. Meth. Fl.* **2014**, *76*, 699-721.

24. Qian, X.; Dede, E.M. Topology optimization of a coupled thermal-fluid system under a tangential thermal gradient constraint. *Struct. Multidiscipl. Optim.* **2016**, *54*, 531-551.

25. Dugast, F.; Favennec, Y.; Josset, C.; Fan, Y.; Luo, L. Topology optimization of thermal fluid flows with an adjoint Lattice Boltzmann method. *J. Comput. Phys.* **2018**, *365*, 376-404.

26. Subramaniam, V.; Dbouk, T.; Harion, J.L. Topology optimization of conjugate heat transfer systems: A competition between heat transfer enhancement and pressure drop reduction. *Int. J. Heat Fluid Flow* **2019**, *75*, 165-184.

27. Duan, X.B.; Li, F.F.; Qin, X.Q. Adaptive mesh method for topology optimization of fluid flow. *Appl. Math. Lett.* **2015**, *44*, 40-44.
28. Jensen, K.E. Topology optimization of Stokes flow on dynamic meshes using simple optimizers. *Comp. Fluids* **2018**, 174, 66–77.

29. Duan, X.B.; Qin, X.Q. Optimality criteria coupled adaptive mesh method for optimal shape design of Stokes flow. *Math. Methods Appl. Sci.* **2016**, 39, 3910–3920.

30. Duan, X.B.; Li, F.F.; Qin, X.Q. Topology optimization of incompressible Navier-Stokes problem by level set based adaptive mesh method. *Comput. Math. Appl.* **2016**, 72, 1131–1141.

31. Garcke, H.; Hecht, C.; Hinze, M.; Kahle, C. Numerical approximation of phase field based shape and topology optimization for fluids. *SIAM J. Sci. Comput.* **2015**, 37, A1846–A1871.

32. Feppon, F.; Allaire, G.; Bordeu, F.; Cortial, J.; Dapogny, C. Shape optimization of a coupled thermal fluid-structure problem in a level set mesh evolution framework. *SeMA* **2019**, 76, 413–458.

33. Feppon, F.; Allaire, G.; Dapogny, C.; Jolivet, P. Body-fitted topology optimization of 2D and 3D fluid-to-fluid heat exchangers. *Comput. Methods Appl. Mech. Engrg.* **2021**, 376, 113638.

34. Alauzet, F.; Loseille, A. A decade of progress on anisotropic mesh adaptation for computational fluid dynamics. *Comput. Aided Des.* **2016**, 72, 13–39.

35. Sari, J.; Cremonesi, F.; Khalloufi, M.; Cauneau, F.; Meliga, P.; Mesri, Y.; Hachem, E. Anisotropic adaptive stabilized finite element solver for RANS models. *Int. J. Numer. Meth. Fl.* **2018**, 86, 717–736.

36. Dambrine, M.; Kateb, D. On the ersatz material approximation in level-set methods. *ESAIM Contr. Optim. Ca.* **2010**, 16, 618–634.

37. Othmer, C. A continuous adjoint formulation for the computation of topological and surface sensitivities of ducted flows. *Int. J. Numer. Meth. Fl.* **2008**, 58, 861–877.

38. Hinterberger, C.; Olesen, M. Automatic geometry optimization of exhaust systems based on sensitivities computed by a continuous adjoint CFD method in OpenFOAM. Technical Paper 2010-01-1278, SAE, 2010.

39. Hinterberger, C.; Olesen, M. Industrial application of continuous adjoint flow solvers for the optimization of automotive exhaust systems. In *Proceedings of the Procs. of the ECCOMAS Thematic Conference on CFD & Optimization Methods and Applications, Antalya, Turkey, ECCOMAS 2011.*, 2011, pp. 1–17.

40. Zeng, S.; Kanargi, B.; Lee, P.S. Experimental and numerical investigation of a mini channel forced air heat sink designed by topology optimization. *Int. J. Heat Mass Transfer* **2018**, 121, 663–679.

41. Athan, T.W.; Papalambros, P.Y. A note on weighted criteria methods for compromise solutions in multi-objective optimization. *Eng. Optim.* **1996**, 27, 155–176.

42. Hachem, E.; Digonnet, H.; Massoni, E.; Coupez, T. Immersed volume method for solving natural convection, conduction and radiation of a hat-shaped disk inside a 3D enclosure. *Int. J. Numer. Method H.* **2012**, 22, 718–741.

43. Hachem, E.; Feghali, S.; Codina, R.; Coupez, T. Immersed stress method for fluid-structure interaction using anisotropic mesh adaptation. *Int. J. Numer. Meth. Eng.* **2013**, 94, 805–825.

44. Patankar, S.V. *Numerical heat transfer and fluid flow*; Taylor & Francis, 1980.

45. Patankar, S.V. A numerical method for conduction in composite materials, flow in irregular geometries and conjugate heat transfer. In *Proceedings of the Procs. of the 6th International Heat Transfer Conference, 1978*, pp. 297–302.

46. Ville, L.; Silva, L.; Coupez, T. Convected level set method for the numerical simulation of fluid buckling. *Int. J. Numer. Meth. Fl.* **2011**, 66, 324–344.

47. Coupez, T.; Silva, L.; Hachem, E. Implicit boundary and adaptive anisotropic meshing. In *New challenges in grid generation and adaptivity for scientific computing; Perotto, S.; Formaggia, L., Eds.; Springer, 2015*; pp. 1–18.

48. Bonito, A.; Guermond, J.L.; Lee, S. Numerical simulations of bouncing jets. *Int. J. Numer. Meth. Fl.* **2016**, 80, 53–75.

49. Coupez, T. Metric construction by length distribution tensor and edge based error for anisotropic adaptive meshing. *J. Comput. Phys.* **2011**, 230, 2391–2405.

50. Jannoun, G.; Hachem, E.; Veysset, J.; Coupez, T. Anisotropic meshing with time-stepping control for unsteady convection-dominated problems. *Appl. Math. Model.* **2015**, 39, 1899–1916.

51. Coupez, T. Génération de maillage et adaptation de maillage par optimisation locale. *Rev. Eur. Elem. Finis* **2000**, 9, 403–423.

52. Coupez, T.; Jannoun, G.; Veysset, J.; Hachem, E. Edge-based anisotropic mesh adaptation for CFD applications. In *Proceedings of the Procs. of the 21st International Meshing Roundtable; Jiao, X.; Weill, J.C., Eds. Springer, 2013*, pp. 567–583.

53. Meliga, P.; Chomaz, J.M.; Sipp, D. Unsteadiness in the wake of disks and spheres: Instability, receptivity and control using direct and adjoint global stability analyses. *J. Fluid Struct.* **2009**, 25, 601–616.

54. Balay, S.; Abhyankar, S.; Adams, M.F.; Brown, J.; Brune, P.; Buschelman, K.; Dalcin, L.; Dener, A.; Eijkhout, V.; Gropp, W.; et al. PETSc Users Manual (Rev. 3.13). Technical report ANL-95/11-Rev. 3.13, Argonne National Lab., 2020.

55. Marck, G.; Nemer, M.; Harion, J.L. Topology optimization of heat and mass transfer problems: laminar flow. *Numer. Heat Tr. B-Fund.* **2013**, 63, 508–539.

56. Takagi, J.; Yamada, M.; Yasuda, M.; M., S. Continuous particle separation in a microchannel having asymmetrically arranged multiple branches. *Lab Chip* **2005**, 5, 778–784.

57. Liu, Z.; Gao, Q.; Zhang, P.; Xuan, M.; Wu, Y. Topology optimization of fluid channels with flow rate equality constraints. *Struct. Multidiscipl. Optim.* **2011**, 44, 31–37.
58. Zhou, T.; Liu, T.; Deng, Y.; Chen, L.; Qian, S.; Liu, Z. Design of microfluidic channel networks with specified output flow rates using the CFD-based optimization method. *Microfluid. Nanofluid.* **2017**, *21*, 1–8.

59. Gargantini, G. Design of bi-fluid heat exchangers using topology optimization. PhD thesis, Master Thesis Politecnico di Milano, 2020.

60. Villanueva, C.H.; Maute, K. CutFEM topology optimization of 3D laminar incompressible flow problems. *Comput. Methods Appl. Mech. Engrg.* **2017**, *320*, 444–473.

61. Flageul, C.; Benhamadouche, S.; Lamballais, É.; Laurence, D. DNS of turbulent channel flow with conjugate heat transfer: Effect of thermal boundary conditions on the second moments and budgets. *Int. J. Heat Fluid Flow* **2015**, *55*, 34–44.

62. Chomaz, J.M. Global instabilities in spatially developing flows: Non-normality and nonlinearity. *Annu. Rev. Fluid Mech.* **2005**, *37*, 357–392.

63. Soto, O.; Löhner, R. On the computation of flow sensitivities from boundary integrals. *AIAA-2004-0112 2004*.

64. Hachem, E.; Rivaux, B.; Kloczko, T.; Digonnet, H.; Coupez, T. Stabilized finite element method for incompressible flows with high Reynolds number. *J. Comput. Phys.* **2010**, *229*, 8643–8665.

65. Codina, R. Stabilized finite element approximation of transient incompressible flows using orthogonal subscales. *Comput. Methods Appl. Mech. Engrg.* **2002**, *191*, 4295–4321.

66. Tezduyar, T.; Osawa, Y. Finite element stabilization parameters computed from element matrices and vectors. *Comput. Methods Appl. Mech. Engrg.* **2000**, *190*, 411–430.

67. Codina, R. Comparison of some finite element methods for solving the diffusion-convection-reaction equation. *Comput. Methods Appl. Mech. Engrg.* **1998**, *156*, 185–210.

68. Badia, S.; Codina, R. Analysis of a stabilized finite element approximation of the transient convection-diffusion equation using an ALE framework. *SIAM J. Numer. Anal.* **2006**, *44*, 2159–2197.

69. Brooks, A.; Hughes, T.J.R. Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations. *Comput. Methods Appl. Mech. Engrg.* **1982**, *32*, 199–259.

70. Galeão, A.C.; Do Carmo, E.G.D. A consistent approximate upwind Petrov-Galerkin method for convection-dominated problems. *Comput. Methods Appl. Mech. Engrg.* **1988**, *68*, 83–95.

**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.