**3. The heat-transfer problem in the context of topology optimization**

The heat-transfer problem (as shown in **Figure 5**), in its weak form in the design domain, can be generalized as

$$\begin{aligned} \rho c \, T\_t &= k \Delta T + q\_v & \text{in } \Omega \times [0, t] \\ T &= T\_b & \text{on } \Gamma\_1 \\ -k \frac{\partial T}{\partial n} &= q & \text{on } \Gamma\_2 \\ -k \frac{\partial T}{\partial n} &= h (T - T\_\omega) & \text{on } \Gamma\_3 \\ T \vert\_{t=0} &= T\_\circ \end{aligned} \tag{1}$$

where *ρ* is the material density, *c* is the specific heat of material, *Tt* is the temperature for a particular given time in transient cases, *k* is the thermal conductivity and *qv* is an internal heat generation rate per unit volume. In general, three types of boundary conditions may exist and can be considered: a temperature condition on Γ<sup>1</sup> , a heat flux conduction on Γ<sup>2</sup> , and a convective condition on Γ<sup>3</sup> . *T*<sup>0</sup> is the initial temperature at time *t* = 0, *Tb* is a temperature imposed on Γ<sup>1</sup> , *q* is a heat flux boundary condition imposed on Γ<sup>2</sup> , *h* is the convective heat-transfer coefficient

**Figure 5.** Generalized heat-transfer problem.

method are termed 'material interpolation', 'artificial material', 'power law' and 'solid isotropic material with penalization (SIMP)' methods. Although SIMP is only one of the methods, its popularity has led for the term to be colloquially used in place of density-based methods. Density-based methods are an extension of the works on the homogenization method. This type of method has experienced much popularity in recent years in this community due to its conceptual simplicity and ease in implementation. Nearly all commercial topology optimiza-

Similarly with the homogenization method, these density-based methods operate on fixed domain of finite elements. The main difference is that, rather than a set of microstructure properties, each finite element contains only a single design variable. This variable is often

ment concerned with the physics involved, for example, the elastic modulus for structural problems or thermal conductivity for heat-transfer problems, is made as a function of the density design variable. This is usually accomplished by utilizing an interpolation function. The topology generated in **Figure 1** was based on this method. Tremendous amount of literature is available for this method and the book [13] contains much discussion on this method as well as an '99-line code' for MATLAB which pioneered the publication of codes for educational purposes in topology optimization. It has been reworked by Andreassen et al. in [35] which shortened the code as well as greatly improving its efficiency. Another rework was made by Liu et al. in [36] which provides the code's extension to 3D problems in the MATLAB environment. More recently, Aage et al. [37] has released their code which utilized Portable, Extensible Toolkit for Scientific Computation (PETSc) and can handle problem scales which

**3. The heat-transfer problem in the context of topology optimization**

The heat-transfer problem (as shown in **Figure 5**), in its weak form in the design domain, can

*ρc Tt* = *kΔT* + *qv* in Ω × [0, *t*]

*T* = *Tb* on Γ<sup>1</sup> <sup>−</sup>*<sup>k</sup>* \_\_\_ <sup>∂</sup>*<sup>T</sup>* <sup>∂</sup>*<sup>n</sup>* <sup>=</sup> *<sup>q</sup>* on <sup>Γ</sup>2

<sup>∂</sup>*<sup>n</sup>* <sup>=</sup> *<sup>h</sup>*(*<sup>T</sup>* <sup>−</sup> *<sup>T</sup>*<sup>∞</sup> ) on <sup>Γ</sup><sup>3</sup>

generation rate per unit volume. In general, three types of boundary conditions may exist and

is the initial temperature at time *t* = 0, *Tb*

−*k* \_\_\_ <sup>∂</sup>*<sup>T</sup>*

where *ρ* is the material density, *c* is the specific heat of material, *Tt*

can be considered: a temperature condition on Γ<sup>1</sup>

*q* is a heat flux boundary condition imposed on Γ<sup>2</sup>

. *T*<sup>0</sup>

*T*|*<sup>t</sup>*=0 = *T*<sup>0</sup>

particular given time in transient cases, *k* is the thermal conductivity and *qv*

, a heat flux conduction on Γ<sup>2</sup>

(1)

,

is the temperature for a

is a temperature imposed on Γ<sup>1</sup>

, *h* is the convective heat-transfer coefficient

is an internal heat

, and a convec-

. The relevant material property of each ele-

tion tools utilize a density-based method [34].

68 Heat Exchangers– Design, Experiment and Simulation

understood as the element material density, *ρ<sup>e</sup>*

are not practical in MATLAB.

be generalized as

tive condition on Γ<sup>3</sup>

on Γ<sup>3</sup> and *T∞* is a fixed reference temperature, *n* is the boundary normal vector. Special treatment is needed for methods which produce intermediate densities for problems considering convective boundary conditions since boundaries are not well defined.

Simplifying to a steady-state heat pure heat conduction case with only temperature boundary conditions and heat flux boundary conditions considered, Eq. (**1**) is reduced to

$$\begin{cases} q\_v = k \Delta T & \text{in } \Omega \\ T = T\_b & \text{on } \Gamma\_1 \\ -k \frac{\partial T}{\partial u} = q & \text{on } \Gamma\_2 \end{cases} \tag{2}$$

This form is often considered for the 'volume-to-point' problem commonly investigated in heat conduction problems. Using the virtual temperature field, *v*, the weak formulation of the heat conduction problem is given by

$$\int\_{\Omega} (\mathbf{k} \Delta \mathbf{T} - q\_v) \mathbf{v} d\mathbf{x} = \mathbf{0} \tag{3}$$

After integration by parts has been carried out and applying the heat flux boundary condition, Eq. (**3**) becomes

$$\int\_{\Omega} k \,\nabla \cdot T \cdot \nabla \cdot \nu - q\_v \nu d\mathbf{x} = \int\_{\Gamma\_v} k \,\nabla \cdot T \cdot mvd\Gamma = -\int\_{\Gamma\_v} qvd\Gamma \tag{4}$$

And the weak form can be written as

$$a(T, \nu) = \ell(\nu) \,\forall \nu \in \tilde{\mathbf{T}} \,\tilde{\mathbf{T}} = \left\{ T \in \left[ H^1(\Omega) \right] \, \middle| \, T = T\_{\flat} \text{ on } \mathbf{x} \in \Gamma\_1 \right\} \tag{5}$$

where *v* is in **T**˜ and **T˜** is a subset of a Sobolev space. The left-hand side *a*(*T, v*) represents the energy bilinear form. It is obtained from Eq. (**4**) and is given as

$$a(T, \nu) = \int\_{\Omega} k \,\,\nabla \cdot T \cdot \nabla \,\,\nu \tag{6}$$

The (*ν* ) term is called the thermal load linear form and can similarly be obtained from Eq. (**4**) and is given as

$$\mathcal{U}(\nu) = \underset{\alpha}{\int} q\_{\nu} \,\nu d\mathbf{x} - \underset{\mathbf{f}\_z}{\int} q\nu d\Gamma \tag{7}$$

This is often used for deriving the propagating velocity of the material boundaries by the material derivative theory in boundary variation methods and the homogenization method. One design objective or thermal compliance measure, *c*, that is considered as the mean temperature could be expressed as

$$c(\Omega) = \int\_{\Omega} k \cdot \nabla \cdot T \cdot \nabla \cdot \text{Tdx} \tag{8}$$

And finally the topology optimization problem of the heat conduction problem is expressed as

$$\begin{aligned} \text{minimize} & \text{provenn} \text{ or } \text{ue} \text{ near } \text{conuncu} \text{ or } \text{provenn} \text{ is expressed} \\\\ \min\_{\Omega \in \mathcal{D}} & \quad c(\Omega) \\\\ \text{s.t.} & \quad a(T, \nu) = t(\nu), \quad \text{for all } \nu \in \tilde{\mathbf{T}} \\ & \quad \int\_{\mathcal{D}} d\mathbf{x} \le V\_{\text{max}} \end{aligned} \tag{9}$$

where *D* contains the material distributed in the design domain Ω and *V*max is a volumetric constraint. In the context of density-based topology optimization, we introduce the element density, *ρ<sup>e</sup>* , and applying a discretized optimization model, for example, FEA, the heat conduction problem is defined as

$$\begin{aligned} \min\_{\boldsymbol{\rho}\_{\boldsymbol{\epsilon}}} \quad & c(\boldsymbol{\rho}\_{\boldsymbol{\epsilon}}) = \mathbf{T}^{\mathrm{T}} \mathbf{K} \mathbf{T} = \mathbf{Q}^{\mathrm{T}} \mathbf{T} \\ \text{s.t.} \quad & \mathbf{K} \mathbf{T} = \mathbf{Q} \\ & V\_{\boldsymbol{\epsilon}} = V(\boldsymbol{\rho}\_{\boldsymbol{\epsilon}}) / V\_{\boldsymbol{\epsilon}} \\ & 0 < \boldsymbol{\rho}\_{\boldsymbol{\epsilon}} \le 1 \end{aligned} \tag{10}$$

where K is the global thermal conductivity matrix, **T** is the node temperature vector and **Q** is the applied thermal load. It is to be noted that the global thermal conductivity matrix if formed from the assembly of individual element thermal conductivity matrix, **k***<sup>e</sup>* , and the material interpolation schemes is applied here as is formally given as

$$\mathbf{K}(\rho\_{\epsilon}) = \sum\_{\epsilon=1}^{N} k\_{\epsilon f}(\rho\_{\epsilon}) \,\mathbf{k}\_{\epsilon} \tag{11}$$

where *keff* is the material interpolation scheme. The objective function could then be expressed as

$$\mathbf{c}(\boldsymbol{\rho}\_{\boldsymbol{\epsilon}}) = \mathbf{Q}^{\top}\mathbf{T}, \quad \text{where } \mathbf{T} \text{ solves:}\\\sum\_{\boldsymbol{\epsilon}=\mathbf{l}}^{N} \mathbf{k}\_{\boldsymbol{\eta}}(\boldsymbol{\rho}\_{\boldsymbol{\epsilon}}) \, \mathbf{k}\_{\boldsymbol{\epsilon}}\mathbf{T} = \mathbf{Q} \tag{12}$$

One simple form of the interpolation scheme was presented when SIMP was introduced and is given as

$$k\_{\rm eff}(\rho\_{\epsilon}) = \left(k\_{\rm max} - k\_{\rm min}\right) \rho\_{\epsilon}^{\;\;\;\prime} \tag{13}$$

Gradients are usually required by the optimization algorithms needed for the update process in topology optimization. These are easily derived for the objective and constraints involving only *ρ<sup>e</sup>* . For functions that depend also on temperatures, derivative can be obtained using the chain rule. These expressions will then contain derivatives of temperature, which in turn can be obtained by taking the derivative of the equilibrium equation, **KT = Q**. The most effective method for calculating the derivatives is to use the adjoint method, where derivatives of the temperature are not calculated explicitly. For the thermal compliance problem given above, we rewrite the objective function by adding a zero function:

$$c(\rho\_{\iota}) = \mathbf{Q}^{\scriptscriptstyle{\scriptscriptstyle{\scriptscriptstyle{\scriptscriptstyle{\scriptscriptstyle{\scriptscriptstyle{\scriptscriptstyle{\scriptscriptstyle{\scriptscriptstyle{\scriptscriptstyle{\scriptscriptstyle{\scriptscriptstyle{\scriptscriptstyle{\scriptscriptstyle{\scriptscriptstyle{\scriptscriptstyle{\scriptscriptstyle{\scriptscriptstyle{\scriptscriptstyle{\scriptscriptstyle{\scriptscriptstyle{\scriptscriptstyle{\scriptscriptstyle{\scriptscriptstyle{\scriptscriptstyle{\boldsymbol{\prime}}}}}}}}}}}}}}}}}}}}}} \mathbf}}} \mathbf \leftarrow 1} \mathbf \leftarrow 1} \cdot 1 \, ^{\scriptscriptstyle{\!\!\!/} \mathbf{T} \cdots \lambda \, ^{\scriptscriptstyle{\!\!\!/} \mathbf{\dot{\mathcal{K}}} \mathbf{\dot{\mathcal{K}}} \mathbf{\!} \cdots \mathbf{\!} \mathbf{\mathcal{K}}}} \, \mathbf{\!} \mathbf{\!} \mathbf{\!} \mathbf{\!} \mathbf{\!} \mathbf{\!} \mathbf{\!} \mathbf{\!} \mathbf{\!} \mathbf{\!} \mathbf{\!} \mathbf{\!} \mathbf{\!} \mathbf{\!} \mathbf{\!} \mathbf{\!} \mathbf{\!} \mathbf{\!} \mathbf{\!} \mathbf{\!} \mathbf{\!} \mathbf{\!} \mathbf{\!} \mathbf{\!} \mathbf{\!} \mathbf{\!} \mathbf{\!} \mathbf{\!} \mathbf{\!} \mathbf{\!} \mathbf{\!} \mathbf{\!} \mathbf{\!} \mathbf{\!} \mathbf{\!} \mathbf{\!} \mathbf{\!} \mathbf{\!} \mathbf{\!} \mathbf{\!} \mathbf{\!} \mathbf{\!} \mathbf{\!} \mathbf{\!} \mathbf{\!} \mathbf{\!} \mathbf{\!} \mathbf{\!} \mathbf{\!} \mathbf{\!} \mathbf{\!} \mathbf{\!} \mathbf{\!} \mathbf{\!} \mathbf{\!} \mathbf{\!} \mathbf{\!} \mathbf{\!} \mathbf{\!} \mathbf{\!} \mathbf{\!} \mathbf{\!} \mathbf{\!} \mathbf{\!} \mathbf{\!} \mathbf{\!} \mathbf{\!} \mathbf{\!} \mathbf{\!} \mathbf{\!} \mathbf{\!} \mathbf{\!}$$

where *λ* is called the Lagrangian multiplier which is an arbitrary, but fixed real vector. We then obtain the derivative as

$$\frac{\partial c}{\partial \rho\_{\epsilon}} = (\mathbf{Q}^T - \lambda^T \mathbf{K}) \frac{\partial \mathbf{T}}{\partial \rho\_{\epsilon}} - \lambda^T \frac{\partial \mathbf{K}}{\partial \rho\_{\epsilon}} \mathbf{T} \tag{15}$$

which can be re-written as

After integration by parts has been carried out and applying the heat flux boundary condi-

Γ1

where *v* is in **T**˜ and **T˜** is a subset of a Sobolev space. The left-hand side *a*(*T, v*) represents the

Ω

Ω

Ω

*<sup>D</sup> dx* ≤ *V*max

The (*ν* ) term is called the thermal load linear form and can similarly be obtained from Eq. (**4**)

*qv νdx* − ∫ Γ2

This is often used for deriving the propagating velocity of the material boundaries by the material derivative theory in boundary variation methods and the homogenization method. One design objective or thermal compliance measure, *c*, that is considered as the mean tem-

And finally the topology optimization problem of the heat conduction problem is expressed

where *D* contains the material distributed in the design domain Ω and *V*max is a volumetric constraint. In the context of density-based topology optimization, we introduce the element

, and applying a discretized optimization model, for example, FEA, the heat con-

*c*(*ρ<sup>e</sup>* ) = **T***<sup>T</sup>* **KT** = **QT T** s.t. **KT** <sup>=</sup> **<sup>Q</sup>** *Vf* <sup>=</sup> *<sup>V</sup>*(*ρ<sup>e</sup>* ) /*Vo*

0 < *ρ<sup>e</sup>* ≤ 1

s.t. *<sup>a</sup>*(*T*, *<sup>ν</sup>* ) <sup>=</sup> (*<sup>ν</sup>* ) , for all *<sup>ν</sup>* <sup>∈</sup> **<sup>T</sup>**˜

*k* ∇ *T* ⋅ *nνd*Γ= − ∫

(Ω ) ]

Γ2

*qvd*Γ (4)


*k* ∇ *T* ⋅ ∇ *ν* (6)

*k* ∇ *T* ⋅ ∇ *Tdx* (8)

(9)

(10)

*qvd*Γ (7)

*k* ∇ *T* ⋅ ∇ *ν* − *qv νdx* = ∫

tion, Eq. (**3**) becomes

∫

and is given as

as

density, *ρ<sup>e</sup>*

Ω

*a*(*T*, *ν* ) = ∫

*l*(*ν* ) = ∫

*c*(Ω ) = ∫

min <sup>Ω</sup>⊆*<sup>D</sup> <sup>c</sup>*(<sup>Ω</sup> )

∫

min*<sup>ρ</sup><sup>e</sup>*

perature could be expressed as

duction problem is defined as

*a*(*T*, *ν* ) = (*ν* ) ∀*ν* ∈ **T**˜ **T**˜ = {*T* ∈ [ *H*<sup>1</sup>

energy bilinear form. It is obtained from Eq. (**4**) and is given as

And the weak form can be written as

70 Heat Exchangers– Design, Experiment and Simulation

$$\frac{\partial c}{\partial \rho\_{\epsilon}} = -\lambda^{\top} \frac{\partial \mathbf{K}}{\partial \rho\_{\epsilon}} \mathbf{T} \tag{16}$$

When λ satisfies the adjoint equation:

$$\mathbf{Q}^{\mathrm{T}} - \lambda^{\mathrm{T}} \mathbf{K} = \mathbf{0} \tag{17}$$

This equation is in the form of an equilibrium equation and for thermal compliance we see that we obtain directly that λ = **T**. Moreover, the form of the stiffness matrix means that the derivatives of the thermal compliance *c*(*ρ<sup>e</sup>* ) for the main problem in Eq. (**10**), considering the SIMP interpolation as presented in Eq. (**13**), are

$$\frac{\partial c}{\partial \rho\_{\epsilon}} = -p(k\_{\text{max}} - k\_{\text{min}}) \, \rho\_{\epsilon} \, ^{p-1} \mathbf{T}^{\top} \mathbf{K} \mathbf{T} \tag{18}$$

Thus, the derivative for the thermal compliance problem becomes easier to compute. It is also worth noting that the derivative is 'localized' to the element level; however, there is an effect from other design variables hidden in the temperature, **T**. The sensitivity is negative for all elements, so intuitively, additional material in any element decreases compliance, and makes the overall objective go lower. Using this sensitivity information, the material is redistributed and the process is repeated until a convergence criterion for the topology optimization process is attained. Each of the paper in the following chapter discusses the complete topology optimization process with more depth and varies depending on the method they have utilized.

### **4. Chronology**

Interests in topology optimization can be represented by the recent amount of publications and citations over the years as presented in **Figure 6(a)** and **(b)**respectively. Although this figure might not accurately represent the exact number of papers, we can still see that the contribution of the papers related to 'heat' is roughly around 1/20th of the total contributions for topology optimization. It has also been increasing especially within the past decade. The amount of papers that are directly related to heat exchangers is arguably much less in number. The following subsections present a number of papers related to the interest of this chapter in its chronological order. For the completeness of the review, some papers at the end of each year are cited but no further elaborations are made due to access restrictions.

**Figure 6.** Some scholarly metrics [38] for topology optimization: (a) publication count and (b) citation count. Highlighted in red are search results for 'topology optimization' + 'heat'. Note that results for this search are reflected on the red axis at the right-hand side and is scaled 1:10.

### **4.1. Prior to 2005**

This equation is in the form of an equilibrium equation and for thermal compliance we see that we obtain directly that λ = **T**. Moreover, the form of the stiffness matrix means that the

= −*p*(*k*max − *k*min ) *ρ<sup>e</sup>*

Thus, the derivative for the thermal compliance problem becomes easier to compute. It is also worth noting that the derivative is 'localized' to the element level; however, there is an effect from other design variables hidden in the temperature, **T**. The sensitivity is negative for all elements, so intuitively, additional material in any element decreases compliance, and makes the overall objective go lower. Using this sensitivity information, the material is redistributed and the process is repeated until a convergence criterion for the topology optimization process is attained. Each of the paper in the following chapter discusses the complete topology optimization process with more depth and varies depending on the method they have utilized.

Interests in topology optimization can be represented by the recent amount of publications and citations over the years as presented in **Figure 6(a)** and **(b)**respectively. Although this figure might not accurately represent the exact number of papers, we can still see that the contribution of the papers related to 'heat' is roughly around 1/20th of the total contributions for topology optimization. It has also been increasing especially within the past decade. The amount of papers that are directly related to heat exchangers is arguably much less in number. The following subsections present a number of papers related to the interest of this chapter in its chronological order. For the completeness of the review, some papers at the end

of each year are cited but no further elaborations are made due to access restrictions.

**Figure 6.** Some scholarly metrics [38] for topology optimization: (a) publication count and (b) citation count. Highlighted in red are search results for 'topology optimization' + 'heat'. Note that results for this search are reflected on the red axis

) for the main problem in Eq. (**10**), considering the

*<sup>p</sup>*−<sup>1</sup> **T***<sup>T</sup>* **KT** (18)

derivatives of the thermal compliance *c*(*ρ<sup>e</sup>*

72 Heat Exchangers– Design, Experiment and Simulation

\_\_\_ <sup>∂</sup>*<sup>c</sup>*

**4. Chronology**

at the right-hand side and is scaled 1:10.

SIMP interpolation as presented in Eq. (**13**), are

∂*ρ<sup>e</sup>*

Rodrigues and Fernandez [39, 40] and Jog [41] utilized topology optimization for designing thermoelastic structures. Heat transfer was treated as one of the involved physics and as an extension to structural mechanics problems. It is worth noting that this was the beginning of the consideration for heat-transfer applications for topology optimization. However, in this chapter, we restrict ourselves to papers that focus more on heat conduction topology optimization (and a few convection cases) that is more directly related to cooling applications, such as the case for heat exchanger design.

Bejan [42, 43] introduced constructal theory in the context of heat transfer. Although it is not directly categorized as topology optimization due to restrictions on size and orientation of each building block, it has provided interesting discussions and has formulated a fundamental problem for the heat-transfer community. The problem is now known as 'volume-to-point (VP) problem' or 'access problem' and discusses the need to layout a fixed amount of material in a heat-generating domain (such as a CPU).

Xie et al. [44] used ESO explicitly for conduction problems. Several generalized claims were given regarding topology optimization, which might not necessarily be true on other methods. The paper is recognized as the first topology optimization paper presented directly solving pure conduction problems. In this paper, an element's rejection is based on the integral of different thermal parameters, more specifically, integral of the temperature surrounding the element. They have highlighted the simplicity of the ESO method to generate novel structures and had considered anisotropic cases in one of the examples. He had also multiple loading cases and had presented two ways to introduce the loading cases, which generated distinct designs.

Turteltaub [45] used SIMP for finding optimal material properties for transient heat conduction problems. Although the generated final designs were rich in intermediate densities due to the lack of penalization, this paper had first offered the possibility to extend topology optimization for transient problems. It was also mentioned that in the heat-transfer problems, special care should be given especially for convective boundary conditions. Though he did not use any explicit boundary-tracking scheme, it was already recognized that difficulty in convective boundaries are present.

Haslinger et al. [46] applied the original homogenization method for conducting structures. Although the paper had focused more on convergence analysis and approximation strategies, it has utilized rank-two laminated structures to demonstrate the optimal heat conductor configurations for its test problems. The effective conductivity of the rank-one laminates was assumed to consist of harmonic and arithmetic means. Numerical minimization was performed by a subroutine from the Numerical Algorithms Group (NAG) Numerical Library which implemented a sequential quadratic programming (SQP) approach.

Cheng et al. [47] introduced bionic optimization strategy for constructing better performing conductive paths. This was directly addressing and comparing results with Bejan's original work. There is not much detail regarding their implementation but it can be viewed as a heuristic 'hard-kill' method. In the same year, Novotny et al. [15] introduced the concept of topological derivative. Although it is viewed as a 'hard-kill' method since it explicitly creates holes in the design domain, the concept is very far from ESO since it utilizes concepts from shape sensitivity analysis to evaluate the topological derivatives. Several theorems were presented in how it can be utilized for the design of conducting structures and has considered Robin boundary conditions in the formation of new holes. In a seemingly unrelated development, Borvall and Petersson [48] introduced the use of topology optimization for fluidic systems governed by Stokes flow. This field of topology optimization has its own unique developments and only a few papers which are relevant to the context of this chapter are mentioned. In this year, Bendsoe and Sigmund [49] also published their book on topology optimization which had some mention of heat-transfer topology as well as instructions for converting the learning code to conduction heat-transfer topology optimization. Guo et al. [50] presented the least dissipation principle.

Xie [51] presented some changes in the ESO method for heat conduction applications which includes some revision for the criterion for rejection through some sensitivity measure. It is also mentionable that this paper had contained a good compilation of literature for shape sensitivity analysis in the field of heat transfer. It was not explicitly stated but the methods implemented were not as aggressive to the original ESO paper where degeneration was considered. This is more properly termed nowadays as 'soft-kill' ESO. Also, the design variable was constructed in terms of the element's thermal conductivity. New interesting problems are given in the context of proper insulation design. Alberto and Sigmund [52] also published on multiphysics problems governed by Poisson's equation, which includes conduction heat transfer. Ha et al. [53] presented non-linear heat conduction problems. Moon et al. [54] presented reliability-based topology optimization considering convection heat transfer.

In this transitory stage, we can see that most of the existing methods are directly being migrated from structural topology optimization into the context of heat transfer. Here, we see ESO, homogenization method and SIMP which is more complex and harder to understand compared to papers in the next years. It is also worth noting that SIMP has already considered transient problems. The topological derivative is also introduced first in the context of heattransfer problems which will later be a very powerful addition to level-set methods. It can also be said that in this year, fluid flow topology optimization has just started.

#### **4.2. 2005–2010**

Yoon and Kim [55] introduced an element connectivity parameterization (ECP) to alleviate problems in applying SIMP to multiphysics problems. A more specific problem of temperature undershooting was emphasized as a numerical instability when applying SIMP to include heat convection formulations on the generated structure boundaries (termed as 'side convection' in the paper). These undershootings in temperatures were deemed to be impossible and infeasible solutions which needed to be strongly addressed. Thus, their paper had given special attention to heat transfer utilizing the zero-length heat conductors as element connectivity measures in ECP. Good results were obtained using the method which was again extended to heat-dissipating structures and electro-thermal actuators. In the same year, Ha and Cho [56] introduced the level-set approach explicitly for heat conduction problems. Their paper contains a detailed yet understandable introduction for level-set methods in the context of heat-conducting structures. It is also worth noting that due to the nature of the level-set method of clear and well-defined boundaries, convection heat transfer was already considered. However, it was not directly applied to the evolving boundaries. Also, it was well reported in this pioneering paper for level-set method for heat-transfer topology optimization that the generated structures were highly dependent on the initial distribution of holes in the design domain since this implementation cannot create new holes during the optimization process. It was also mentioned that density-based method (SIMP) yielded better results for most cases, in terms of the number of iterations needed to achieve the converged results. Thermal compliance values for both methods were in very close agreements.

topological derivative. Although it is viewed as a 'hard-kill' method since it explicitly creates holes in the design domain, the concept is very far from ESO since it utilizes concepts from shape sensitivity analysis to evaluate the topological derivatives. Several theorems were presented in how it can be utilized for the design of conducting structures and has considered Robin boundary conditions in the formation of new holes. In a seemingly unrelated development, Borvall and Petersson [48] introduced the use of topology optimization for fluidic systems governed by Stokes flow. This field of topology optimization has its own unique developments and only a few papers which are relevant to the context of this chapter are mentioned. In this year, Bendsoe and Sigmund [49] also published their book on topology optimization which had some mention of heat-transfer topology as well as instructions for converting the learning code to conduction heat-transfer topology optimization. Guo et al.

Xie [51] presented some changes in the ESO method for heat conduction applications which includes some revision for the criterion for rejection through some sensitivity measure. It is also mentionable that this paper had contained a good compilation of literature for shape sensitivity analysis in the field of heat transfer. It was not explicitly stated but the methods implemented were not as aggressive to the original ESO paper where degeneration was considered. This is more properly termed nowadays as 'soft-kill' ESO. Also, the design variable was constructed in terms of the element's thermal conductivity. New interesting problems are given in the context of proper insulation design. Alberto and Sigmund [52] also published on multiphysics problems governed by Poisson's equation, which includes conduction heat transfer. Ha et al. [53] presented non-linear heat conduction problems. Moon et al. [54] pre-

sented reliability-based topology optimization considering convection heat transfer.

be said that in this year, fluid flow topology optimization has just started.

In this transitory stage, we can see that most of the existing methods are directly being migrated from structural topology optimization into the context of heat transfer. Here, we see ESO, homogenization method and SIMP which is more complex and harder to understand compared to papers in the next years. It is also worth noting that SIMP has already considered transient problems. The topological derivative is also introduced first in the context of heattransfer problems which will later be a very powerful addition to level-set methods. It can also

Yoon and Kim [55] introduced an element connectivity parameterization (ECP) to alleviate problems in applying SIMP to multiphysics problems. A more specific problem of temperature undershooting was emphasized as a numerical instability when applying SIMP to include heat convection formulations on the generated structure boundaries (termed as 'side convection' in the paper). These undershootings in temperatures were deemed to be impossible and infeasible solutions which needed to be strongly addressed. Thus, their paper had given special attention to heat transfer utilizing the zero-length heat conductors as element connectivity measures in ECP. Good results were obtained using the method which was again extended to heat-dissipating structures and electro-thermal actuators. In the same year, Ha and Cho [56] introduced the level-set approach explicitly for heat conduction problems.

[50] presented the least dissipation principle.

74 Heat Exchangers– Design, Experiment and Simulation

**4.2. 2005–2010**

Gersborg-Hansen et al. [57] introduced the use of the finite volume method (FVM) for heat conduction problems. It is worth noting that all previous papers were utilizing finite element methods and formulations. Their justification for utilizing FVM was made in the context of guaranteeing element-wise conservation of the physical quantities and to give access to FVM users to topology optimization. Element interface heat fluxes were calculated using the value of thermal conductivities based on the arithmetic and harmonic means of the surrounding nodes. The SIMP method was utilized for their implementations. Two unique compliance measures were investigated. It was mentioned that the results from FEM and FVM were qualitatively similar and the designs suffered high-mesh dependence when the compliance measure for arithmetic average was used even though the penalty value, *p*, in SIMP was increased up to 5 using the continuation approach. Using the harmonic average in the FVM formulation also reduces checkerboard formation up in their test cases. Donoso [58] revisited the VP problem in 3D space and used the optimality criteria (OC) method to find the solution.

Zhuang et al. [59] utilized the concept of topological derivative in conjunction with the levelset method. The topological derivative was used to create new holes during the topology optimization process and thus eliminating the dependence on the initial hole distribution. A fixed cutting ratio was set for the topological derivative for generating new holes. Multiple load cases were also one of the highlights of their paper and highly consistent results and convergence curves were presented. Xu et al. [60], on the other hand, presented a combinatorial approach for optimizing the heat conduction paths. In their paper, they tried to solve the volume-topoint problem using simulated annealing and genetic algorithm (hard-kill approaches). Their paper had clearly presented their implementation scheme and had made comparisons with the results of bionic optimization. The optimal results were generalized as all high conductivity materials are continuous, no holes are present. For cases in which the thermal conductivity ratio is relatively small, shapes are thick and short surrounding the heat sink. When the thermal conductivity ratio is increased, the shape becomes more slender. Mathieu-Potvin and Gosselin [61] developed an evolutionary algorithm which tries to solve the VP problem. Their evolutionary algorithm aimed to minimize the hotspot temperature by displacing elements. Displacements were either by swapping of a heat-generating cell with a void cell and swapping a heat-generating cell with a conductive cell based on heat flux or by element-averaged heat flux and temperature. It is worth noting that in their implementation, an extended domain was utilized and during the evolution process, the cell elements can rearrange themselves in the extended domain. Due to the nature of the algorithm, exact repeatability of results is most unlikely but measures were adopted to find approximate performances and determine the algorithm's robustness. Results were also compared to constructal theory in terms of the temperature objective function, *kϕ*, dimensionless distance measure, uniformity distribution measure and fractal dimension. Good discussions were made regarding each of these performance measures. Also in the same year, Bruns [62] clarified and resolved the problems presented by Yoon and Kim in their 2005 paper for convection-dominated heat-transfer problems in the context of density-based topology optimization. He has discussed the necessary techniques to prevent the 'undershoot' in temperatures mentioned by ensuring that the convection term contributions are treated as lumped matrix. Side convection terms are weighed by a density difference interpolation scheme and half of the total contribution is associated with two elements connected along the same edge. He has also used the SINH (pronounced as 'cinch') [63] method. He has concluded that poor convection modelling can greatly influence the design process. Kim et al. [64] reconsidered the printed circuit board (PCB) cooling problem but had included mechanical constraints. Zhuang et al. [65] minimized the quadratic mean temperature using the level-set method. Yoo and Kim [66] considered three-dimensional cooling fins using the ECP method.

He and Liu [67] used the bi-directional evolutionary structural optimization (BESO). BESO is differentiated from ESO since BESO allows element addition which is not allowed in ESO. Using a uniform heat distribution problem, he compared the results with SIMP-based solutions. Special attention was given to the lack of intermediate elements in the ESO results, thus, easier manufacturability. Gao et al. [68] published another BESO paper and have considered both design-independent and -dependent loading. Design-dependent loading in this paper was defined as heat loads that vary whether or not material is present. In other words, without the presence of material heat cannot be generated. One case was presented to elaborate the difference and the effect of this assumption. Zhang and Liu [69] mentioned a new method for designing heat-conducting paths based on SIMP. This is related to a later publication mentioned in 2011. Yamasaki et al. [70] presented level-set method for both vibration and heat conduction problems.

Iga et al. [71] introduced convection and heat generation design-dependent effects. He has used a different homogenization approach (termed as the homogenization design method in their paper) and defined a hat function in which the convection boundary conditions are easily applied. The hat function serves as the boundaries between the solid and the void regions. Interest in this paper is given for the utilization of a surrogate model for several fin models for including a better representation of the convection condition. They have also utilized sequential linear programming (SLP) for the update process during topology optimization. They have presented several examples which exhibit the adverse effects if an inappropriate convection modelling is used. Marczak and Anflor [72] introduced the boundary element method (BEM) as an alternative to FEA and FVM. In their paper, topological derivative was used as the means to generate the optimal topologies. BEM is differentiated from FEA and FVM since it does not directly compute based on cells or elements. BEM is considered as 'mesh-free' methods. Although nodes are still present inside the modelled domain, they are treated more as 'recovery points'. Their examples and results were compared to the first ESO paper by Xie et al. in 1999. Dede [73] presented the use of COMSOL Multiphysics coupled with MATLAB for multiphysics topology optimization of heat and flow systems. Kim et al. [74] considered non-linear heat conduction and had designed structures based on the level set with topological derivatives. Pingen and Meyer [75] presented topology optimization for thermal transport.

Yoon [76] considered a sequential computational procedure to design heat-dissipating structures that considers forced convective heat transfer. A staggered approach was used where the flow field was solved first. Artificial damping force was introduced to the Navier-Stokes equation, which was similar to techniques used in immersed boundary methods (IBMs). A total of four material properties were interpolated in his implementation. He had utilized density-based approach and SIMP interpolations for the material properties. Kim et al. [77] compared results for different sensitivity analyses formulations. They have reported the computational time for the finite difference method and two different design sensitivity analyses (DSA). It was reported that a large difference in terms of performance was present between the direct differentiation method and the adjoint variable method (a factor of about 142). The SIMP method was utilized but was not mentioned explicitly. A 3D example was also provided in one of their examples which considered a single convection boundary condition. Dede [78] performed investigations on topology-optimized designs for impinging jets. Single-jet geometry was investigated from coupled thermal-fluidic simulations in a commercial software package. SIMP-based topology optimization was then performed in MATLAB with MMA. The result from the single impinging jet was made as basis for a textured surface geometry for a 3D slot jet. It is worth noting that the two-dimensional (2D) model was made under the assumptions of laminar flow and the 3D multi-jet structure is expected to fall within the turbulent regime. Zhuang et al. [79] utilized level-set method for the design of multi-material heat-conducting structures.

It can be said that the interest in heat-transfer topology optimization started in this time period. The papers presented in this time period were mostly dedicated and developed for heat transfer and design of heat exchangers. Investigations to include convection heat transfer as well as other design-dependent effects are also evident. Level-set method that is combined with the concept of topological derivative can be treated as state of the art during this time period. Also, 'mesh-less' topology optimization was introduced. SIMP has remained as a key method and its integration for FVM users has been mostly utilized.

#### **4.3. 2011–2015**

unlikely but measures were adopted to find approximate performances and determine the algorithm's robustness. Results were also compared to constructal theory in terms of the temperature objective function, *kϕ*, dimensionless distance measure, uniformity distribution measure and fractal dimension. Good discussions were made regarding each of these performance measures. Also in the same year, Bruns [62] clarified and resolved the problems presented by Yoon and Kim in their 2005 paper for convection-dominated heat-transfer problems in the context of density-based topology optimization. He has discussed the necessary techniques to prevent the 'undershoot' in temperatures mentioned by ensuring that the convection term contributions are treated as lumped matrix. Side convection terms are weighed by a density difference interpolation scheme and half of the total contribution is associated with two elements connected along the same edge. He has also used the SINH (pronounced as 'cinch') [63] method. He has concluded that poor convection modelling can greatly influence the design process. Kim et al. [64] reconsidered the printed circuit board (PCB) cooling problem but had included mechanical constraints. Zhuang et al. [65] minimized the quadratic mean temperature using the level-set method. Yoo and Kim [66] considered three-dimensional cooling fins

He and Liu [67] used the bi-directional evolutionary structural optimization (BESO). BESO is differentiated from ESO since BESO allows element addition which is not allowed in ESO. Using a uniform heat distribution problem, he compared the results with SIMP-based solutions. Special attention was given to the lack of intermediate elements in the ESO results, thus, easier manufacturability. Gao et al. [68] published another BESO paper and have considered both design-independent and -dependent loading. Design-dependent loading in this paper was defined as heat loads that vary whether or not material is present. In other words, without the presence of material heat cannot be generated. One case was presented to elaborate the difference and the effect of this assumption. Zhang and Liu [69] mentioned a new method for designing heat-conducting paths based on SIMP. This is related to a later publication mentioned in 2011. Yamasaki et al. [70] presented level-set method for both vibration and heat

Iga et al. [71] introduced convection and heat generation design-dependent effects. He has used a different homogenization approach (termed as the homogenization design method in their paper) and defined a hat function in which the convection boundary conditions are easily applied. The hat function serves as the boundaries between the solid and the void regions. Interest in this paper is given for the utilization of a surrogate model for several fin models for including a better representation of the convection condition. They have also utilized sequential linear programming (SLP) for the update process during topology optimization. They have presented several examples which exhibit the adverse effects if an inappropriate convection modelling is used. Marczak and Anflor [72] introduced the boundary element method (BEM) as an alternative to FEA and FVM. In their paper, topological derivative was used as the means to generate the optimal topologies. BEM is differentiated from FEA and FVM since it does not directly compute based on cells or elements. BEM is considered as 'mesh-free' methods. Although nodes are still present inside the modelled domain, they are treated more as 'recovery points'. Their examples and results were compared to the first ESO paper by Xie et al. in 1999. Dede [73] presented the use of COMSOL Multiphysics coupled

using the ECP method.

76 Heat Exchangers– Design, Experiment and Simulation

conduction problems.

Yamada et al. [80] utilized the level-set method to include design-dependent effects such as convection boundary loading. A fictitious interface energy term was introduced for the design-dependent boundary conditions. Three-dimensional examples were given which clearly demonstrates clear and smooth optimal configurations. A regularization parameter was also utilized to tune the complexity of the optimal results. Convection loading was based on a fixed value. Zhang et al. [81] emphasized on the objective functional in topology optimization. It was highlighted that the cooling problem, as given by Bejan, needs to minimize the maximum temperature but most problems minimize the heat dissipation efficiency (termed as dissipation of heat potential capacity, DHTPC). A one-dimensional problem was revisited and a new objective of minimizing the geometric average temperature is presented in the context of a topology optimization problem. It was not explicitly mentioned what method was used and it is hypothesized that ESO was used due to the chosen problem and the welldefined boundaries of the optimal results (they have mentioned feasible direction method). Li et al. [82] had used the rational approximation of material property (RAMP) material interpolation scheme, OC based on density approach and a density filter was explored. Papoutsis-Kiachagias et al. [83] presented a constrained topology optimization for laminar and turbulent flows, including heat transfer.

Marck et al. [84] performed multi-objective optimization (MOO) using the SIMP method. The MOO was carried out with the two separate goals of minimizing the average temperature and minimizing the variance in the temperature. A very detailed and elaborate description of the FVM-based topology optimization was given. Tests regarding the mesh dependence, sensitivity and density filters as well as the heat transfer in the domain of the VP problem were carried out. Interestingly, this paper had obtained results which had discontinuity in the structure. Dede [85] optimized and designed multi-pass-branching microchannels with topology optimization as a tool. Gregersen et al. [86] considered finite volume-based topology optimization of coupled fluid dynamic and thermal conduction systems. Lee [87] completed his dissertation for topology optimization of convective cooling systems.

Koga et al. [88] demonstrated the complete product development cycle of a topology-optimized water-cooled heat exchanger. A fully coupled problem was solved using finite element method with some modifications to avoid numerical instabilities. A weighted logarithmic multi-objective function was used which contained a function to represent the power dissipation for the fluid flow and the heat dissipation for the heat-transfer problem. SINH was used in their implementation together with a weighed density filter. The heat exchanger was manufactured through electrical discharge machining and precision CNC milling. The experiments have matched well with the numerical simulations. It is worth noting that although the heat exchanger is a three-dimensional device, 2D modelling was employed for the topology optimization process. Burger et al. [89] explored the 3D solution for the volume to point (called volume to surface in this case) utilizing SIMP with method of moving asymptotes (MMAs) implementation. Full and partial Dirichlet boundary conditions were considered. In the partial Dirichlet boundary condition, only a square surface was given a fixed temperature. In the full Dirichlet boundary condition consideration, a volume of non-designable domain was set and the temperature conditions were set at the surfaces of a small volume before the allowable design domain. Different conductivity ratios varying volumetric constraints were explored as well as multiple boundary condition locations. Tree-like structures with four main branches extending to the corners of the design domain were the dominant optimal design features. Zhuang et al. [90] utilized triangular meshes on a transient heat conduction problem. Level-set method with topological derivative was used for the topology optimization process. Radial basis functions were used for defining the boundaries. A narrow band algorithm on the triangular mesh further improves the numerical efficiency. Dirker and Meyer [91] have performed performance tests for SIMP with MMA in an FVM setting. The VP problem was considered. It is worth noting that this work did not utilize any filtering techniques. A total of seven implementation cases were investigated. Six of the cases used predefined penalization parameters ranging from 1 to 5 in 0.5 intervals. Two volumetric constraints as well as three conductivity ratios were considered. Marck et al. [92] discussed topology optimization for heat and mass transfer problems in great detail for laminar flows. Jing et al. [93] has used BEM and level-set method for 2D heat conduction problems. Matsumori et al. [94] published fluidthermal interaction problems under constant input power. Kontoleontos et al. [95] published an adjoint-based constrained topology optimization for viscous flows, including heat transfer.

and a new objective of minimizing the geometric average temperature is presented in the context of a topology optimization problem. It was not explicitly mentioned what method was used and it is hypothesized that ESO was used due to the chosen problem and the welldefined boundaries of the optimal results (they have mentioned feasible direction method). Li et al. [82] had used the rational approximation of material property (RAMP) material interpolation scheme, OC based on density approach and a density filter was explored. Papoutsis-Kiachagias et al. [83] presented a constrained topology optimization for laminar and turbulent

Marck et al. [84] performed multi-objective optimization (MOO) using the SIMP method. The MOO was carried out with the two separate goals of minimizing the average temperature and minimizing the variance in the temperature. A very detailed and elaborate description of the FVM-based topology optimization was given. Tests regarding the mesh dependence, sensitivity and density filters as well as the heat transfer in the domain of the VP problem were carried out. Interestingly, this paper had obtained results which had discontinuity in the structure. Dede [85] optimized and designed multi-pass-branching microchannels with topology optimization as a tool. Gregersen et al. [86] considered finite volume-based topology optimization of coupled fluid dynamic and thermal conduction systems. Lee [87] completed

Koga et al. [88] demonstrated the complete product development cycle of a topology-optimized water-cooled heat exchanger. A fully coupled problem was solved using finite element method with some modifications to avoid numerical instabilities. A weighted logarithmic multi-objective function was used which contained a function to represent the power dissipation for the fluid flow and the heat dissipation for the heat-transfer problem. SINH was used in their implementation together with a weighed density filter. The heat exchanger was manufactured through electrical discharge machining and precision CNC milling. The experiments have matched well with the numerical simulations. It is worth noting that although the heat exchanger is a three-dimensional device, 2D modelling was employed for the topology optimization process. Burger et al. [89] explored the 3D solution for the volume to point (called volume to surface in this case) utilizing SIMP with method of moving asymptotes (MMAs) implementation. Full and partial Dirichlet boundary conditions were considered. In the partial Dirichlet boundary condition, only a square surface was given a fixed temperature. In the full Dirichlet boundary condition consideration, a volume of non-designable domain was set and the temperature conditions were set at the surfaces of a small volume before the allowable design domain. Different conductivity ratios varying volumetric constraints were explored as well as multiple boundary condition locations. Tree-like structures with four main branches extending to the corners of the design domain were the dominant optimal design features. Zhuang et al. [90] utilized triangular meshes on a transient heat conduction problem. Level-set method with topological derivative was used for the topology optimization process. Radial basis functions were used for defining the boundaries. A narrow band algorithm on the triangular mesh further improves the numerical efficiency. Dirker and Meyer [91] have performed performance tests for SIMP with MMA in an FVM setting. The VP problem was considered. It is worth noting that this work did not utilize any filtering techniques. A total of seven implementation cases were investigated. Six of the cases used predefined penalization

his dissertation for topology optimization of convective cooling systems.

flows, including heat transfer.

78 Heat Exchangers– Design, Experiment and Simulation

Zhuang and Xiong [96] proposed a new compliance measure for transient heat conduction problems. They have suggested that the peak values of the given compliance during the time iterations are to be minimized. SIMP with MMA utilized for this study. The equivalent static load-based topology optimization for transient problems was deemed to be more practical and computationally efficient. Cheng and Chen [97] introduced a non-constrained formulation with a volume-of-solid (VOS) function to represent the bounds of the domain. This work is interesting since the objective function was defined as the heat-transfer index (*Q*˙/*m*). Oevelen and Baelsmans [98] demonstrated solutions to a conjugate heat-transfer problem using a twolayer-reduced model to represent a full-3D solution. A test case considered Stokes flow and a highly branching flow network was obtained. SIMP interpolation scheme was used. They have acknowledged the artificial flow through a solid network if the penalization for the flow equations is not sufficient. They have further explored the effects of target temperature and bottom-layer thickness. Dede et al. [99] utilized topology optimization in the design and justification of novel structures that can shield, focus or reverse heat flux over a target domain. Anisotropic material constituents utilizing two-phase material microstructure descriptions for non-symmetric inclusions embedded in a matrix medium were manipulated to obtain the desired performance of the structures. Results were compared with a test case and extensions to arbitrary geometries were explored. Alexandersen et al. [100] made tremendous efforts for investigating heat topology optimization considering buoyancy forces. In this situation, the strongly coupled physics phenomena are very hard to model and thus making topology optimization for this kind of systems even more cumbersome. They have utilized SIMP interpolations for some of the material properties and density-based methods were sought for their implementations. They have demonstrated that effects of buoyancy affect the generated design significantly and have presented a natural convection heat exchanger as well as a buoyancy-driven micro-pump. In their paper, elements of large-scale simulations are already evident and they mentioned the difficulties they have encountered as well as their proposed solutions to overcome them. Lee [101] presented a multi-material heat conduction problem using a multiphase level set. Jing et al. [102] presented the topological sensitivity of the objective function on morphing boundaries.

Yaji et al. [103] utilized the level-set method to obtain the optimal design for a fully coupled thermo-fluidics problem. Tikhonov-based regularization scheme enabled the qualitative control for the geometric complexity of the generated structures. An optimization algorithm together with a smoothed Heaviside function was needed for the stabilization of the numerical computations. In this paper, both 2D and 3D examples were demonstrated with smooth and well-defined boundaries. Zhuang and Xiong [104] introduced additional temperature constraints on a defined region in the design domain. Their work had considered transient problems and had utilized the equivalent temperature field as a more effective means to solve the time-dependent finite element problem. In addition, they have utilized three materials in some of their examples using SIMP method. Jing et al. [105] utilized the BEM for the implementation of level-set method and considers design-dependent boundary conditions. The level-set method is used to represent the structural boundary and the boundary mesh for the BEM analysis is constructed on the iso-surface of the level-set function. Topological derivative is also utilized to make new holes. Cheng and Chen [106] utilized their volume-of-solid method for the topological design of the laminated metallic composite materials arranged in two predefined configurations. Similar to the previous paper, they have presented two new very interesting objective functions (*Q*˙/*V* and *Q*˙/*USD*). Dede et al. [107] recently demonstrated a complete product cycle development for developing a forced air-cooled heat-sink-made additive manufacturing (AM). They have applied SIMP-based topology optimization and had utilized a modified hat function to define the heat convection loading surface for their problem. A parabolic distribution of the heat-transfer coefficient was assumed in relation to the forced air cooling. Two-dimensional models are first tested and compared to some common heatsink geometries found in the market. A quarter of a 3D model was then implemented and volume reconstruction was also mentioned to obtain a water-tight design suitable for additive manufacturing. Experiments were then conducted and the topology-optimized structures are compared with the commercially available design. Results showed that the designed heat sink performed better compared to other heat sinks but due to the inferior material properties and porous structure of the AM-produced design, it was not performing as to its numerical design specifications. Alexandersen et al. [108] recently published the culmination of their buoyancy flow works by implementing a large-scale three-dimensional model of designed heat sinks. A total of 16.38 million design elements with 83.08 million degrees of freedom were solved in one of their examples for a passive heat-sink cooler for light-emitting diodes (LEDs). Lohan et al. [109] presented generative design algorithms for heat conduction. A dissertation study utilizing boundary element method was recently finished by Jing [109]. Dede [110] designed and fabricated a multi-device single-phase-branching microchannel cold plate.

In this time period, highlight is given to product design cycles and actual realization of topology-optimized designs. It is also evident that trends are going for incorporation of fluid flow either directly (through coupled analysis of both the fluid and heat-transfer domains) or indirectly (through convection boundary conditions). Interests for transient problems have also re-emerged with techniques such as the equivalent temperature field being utilized to reduce the burden of the finite analysis for the governing equations of the system. Level-set method is also evolving rapidly by utilizing other techniques such as topological derivative and BEM to make up for their weak points. Density-based methods, especially SIMP, are still staple with most of the works for 3D modelling and thermo-fluidic systems. Massive implementations with millions of DOFs are also slowly being realized, mostly utilizing density-based methods. As an additional foresight, it can be mentioned that none of the above works have considered radiation effects, though some problem formulations can accommodate radiation by utilizing the convection form of radiation. In the future, this work could be sought but would pose the problem for the discretized method of properly identifying cavities and formations inside the evolving domain. View factor computation is also one complication which would be very expensive to perform since radiating boundaries would change in each iteration.
