Heat Exchanger Design with Topology Optimization

Abstract

Topology optimization is proving to be a valuable design tool for physical systems, especially for structural systems. However, its application in the field of heat transfer is less evident but is constantly progressing. In this chapter, we would like to introduce topology optimization in the context of heat exchanger design to the general reader. We also provide a chronological review of available literature to see the current progress of topology optimization in the field of heat transfer and heat exchanger design. We expect that topology optimization will prove to be a valuable tool in heat exchanger design for the coming years.

1. Introduction

The need for high-performance heat-dissipating devices is highly needed in today’s rapidly changing power device and electronics markets [1, 2]. With worldwide movements on the implementation of Industry 4.0, we will see more radical changes in the way tangible products are manufactured [3]. At the same time, rapid product design cycles are becoming more of a standard rather than a demand. Thus, the need for automated design processes carried out with the use of computer as tools has never been so imperative. Computational design procedures have been more widely accepted during the past decades due to the improvements in computing technologies [4]. Together with this, rapid advancements in the algorithms and automated design procedures have flourished. Topology optimization can be viewed as one of the most promising automated design procedures, which has been an active topic of research for almost three decades.

Topology optimization is an automated, ‘best material layout’ process, which follows the governing equations of one or more physics taken into consideration under a user-defined set of conditions and limitations. Several methods and techniques are already well developed especially for the field of structural engineering. Topology optimization is slowly being used in mainstream design processes of tangible products due to the advancements in computational power of computers, the optimization methods, and techniques used in topology optimization itself.

Computational tools have been developed to aid and answer some of the engineering queries, but the main design of the structure is usually left to experienced and specialized professionals. Commonly applied modern-day topology optimization methods utilize finite element analyses (FEA) where each discretization is treated as a design variable. By choosing and varying the adequate material property related to the investigated case, we would iteratively investigate which element is helpful, thus material is ‘allocated’, and which ones are not, thus can be left as ‘void’, from the design space. We can also set areas that must be filled with material or areas where materials should not be placed. There are a number of learning materials for topology optimization, most are from one research group from Denmark. Among their developments are a free mobile app, TopOpt [5] and TopOpt3D [6], which can execute structural topology optimization and output. STL files ready for three-dimensional (3D) printing. The interface, some common definitions for structural topology optimization and an example are presented in Figure 1.

Figure 1.

TopOpt app [5] developed by DTU reflecting the essential elements of topology optimization.

The earliest work related to topology optimization can be traced back to the ingenious Australian inventor who formulated Michell’s truss theory [7] (named after inventor George Michell). The said theory dealt with the least-volume topology of trusses with a single load condition and a stress constraint. Not only was this imaginatively ingenious, it was also ahead of his time where almost nothing was known about the techniques of structural optimization. His works were ignored for almost half a decade where it was rediscovered by Cox [8, 9] and Owen [10] in the 1960s, the same time when computers were acknowledged as automation tools. It was Hemp [11] and his co-workers who had spent most of their professional lives and comprehensively studied Michell structures. Modern-day computer-aided topology optimization can probably be traced back to the works of Bendsoe and Kikuchi [12] on homogenization who had also coined and popularized the term topology optimization. For the following decades, their works had sparked the interest of many researchers and might not have necessarily had any product-related applications. In the next section, we first briefly discuss the main methods commonly used in topology optimization.

2. Topology optimization methods and learning codes

Different methods have been developed in finding solutions to the optimal layout problem. Since Bendsoe and Kikuchi’s work in 1988 [12], focus has been more on finite element (FE)-based topology optimization of continuum structures. Different methods have been developed since. The differences in the different methodologies lay in the way the design space, and consequently, the design variable is parameterized and controlled. Some methods directly define the design variables on the finite element domain, while others define a separate function from which the generated structure is interpreted. In both cases, ‘0-1’ designs or ‘void/solid’ designs are desired because they can be easily interpreted and physically realized. Here, 1 or solid means that material is allocated on the design element and 0 or void means that material is not present in the design element. In some methods, ‘grey’ or intermediate densities, which are values between 0 and 1, are encountered and observed. The following subsections outline some of the popular methods for topology optimization. We also list references at the end of each method where codes (usually written in MATLAB) are readily available for interested readers. These codes also contain more information on the mathematical background and rationale for each method. These codes, however, are usually written in the context of structural topology optimization but could be modified appropriately to solve heat-transfer problems. The detailed modifications needed are explicitly given in Appendix B of Ref. [13].

2.1. Homogenization method

Pioneering work by Bendsøe and Kikuchi [12] posed a structural layout problem within the context of homogenization theory. In their method, now known as the homogenization method, they treat each element as porous material whose microstructures can be modelled and controlled. By tuning these microstructures, macro-scale material properties are realized which are best suited for the stress experienced from each element. The periodic microstructures are defined for each discretized unit cell in the finite element domain. In their work, they had demonstrated two potential microstructures that could be generated on each unit cell. The first one being a perforated microstructure in the form of a square cell with a rectangular void with three control parameters (μ₁, μ₂ and θ). The second was a layered microstructure with two isotropic constituents with the same control parameters. These two types of microstructure definition are visualized in Figure 2 (a). Under the assumption of infinitesimally small periodic unit cells and the adequate microstructure definition, it was deemed that any anisotropic macro-scale representation of the material can be achieved such as pure solid, pure void, composite and porous material.

Figure 2.

(a) Schematic of the treatment of design variables contained in elements showing two possible microstructures. (b) Result for a truss problem reflecting generated microstructures [14].

A single set of variables corresponding to each microstructure can be used for each design element or can be extended to a sub-mesh to generate finer structures. Topology optimization, in this sense, becomes a problem to determine the optimal combination of these design variables, which corresponds to the optimal macro-scale distribution of properties which minimize a given objective function. This approach was investigated in the 1990s but has received less attention in the recent years due to the emergence of more efficient methods. Nevertheless, it gave the fundamental concepts and ideas in the other methods. Additionally, some methods that will be later mentioned apply alternative formulations to alleviate the common numerical issues found in explicit topology parameterization. Nowadays, it has found its application in finding ways of how to realize high-performing microstructures and is called ‘inverse homogenization’.

2.2. ‘Hard-kill’ methods

‘Hard-kill’ methods are a generalization of methods that explicitly treat each element as material or void. Unlike other methods, they do not relax the ‘0-1’ problem on topology optimization. These methods gradually remove (in some cases add) elements that represent absence (or presence) of a material into the design domain explicitly for each iteration step. A few of these methods utilize combinatorial techniques such as genetic algorithms and simulated annealing, to name a few. Another ‘hard-kill’ method that is based on sensitivity information is known as the concept of using topological derivatives (or topological sensitivity) [15]. The concept of topological derivatives is that undesired computational nodes are explicitly removed. The most well-known ‘hard-kill’ method in topology optimization is the evolutionary structural optimization (ESO) [16] and, more recently, the bi-directional evolutionary structural optimization (BESO) [17]. BESO is differentiated from ESO in a way that ESO only allows for the removal of elements while BESO allows for both the addition and removal of elements that represent the presence or absence of a material based on an ‘optimization criterion’, which is evaluated in each small domain or element. This is analogous to slowly evolving the shape of a structure towards the desired optimum result by removing (or adding) the elements that do not contribute to the improvement of the desired objective function. The choice of material to be removed (or added) is based on heuristic criteria, which is based on sensitivity information of the iteration steps. As a result of these heuristic features, the technicality of this method is often questioned for a robust theoretical basis does not exist [18, 19]. One of the most attractive features of these hard-kill methods is its simplicity with which they can be utilized with commercial finite element packages. It is claimed that the integration of algorithms based on hard-kill methods with finite element analysis (FEA) solvers requires only minor modification in the pre- or post-processing steps [19]. Also, structures generated are free from intermediate or ‘grey’ material representations due to the nature of its solution method of explicitly removing (or adding) material in the finite element system. More recently in [20], BESO has been relaxed to prevent the concerns given in [18, 19] and was termed ‘soft-kill’ method. An attempt to visually present the conceptual differences between the different ‘hard-kill’ methods is presented in Figure 3. In Figure 3 (a), elements are essentially removed from the FEA routine as executed in the original ESO. In Figure 3 (b), a void element is essentially allowed to ‘roam’ on neighbouring elements until such time it finds an optimal location and this is common for combinatorial techniques. In Figure 3 (c), a node is essentially removed and creates an area of void elements, and this is the concept behind the topological derivatives. A MATLAB code of the relaxed BESO implementation is also given in Ref. [20].

Figure 3.

Conceptual differences of different ‘hard-kill’ methods (a) ESO, (b) combinatorial techniques and (c) topological derivative.

2.3. Boundary variation methods

Boundary variation methods are among the most recent and noteworthy contributions that lead to advancements in structural topology optimization. Boundary variation methods have originated in shape optimization techniques and had been recently introduced to structural topology optimization. They are differentiated from the other methods from the fact that structure domain and boundaries are represented based on implicit functions rather than an explicit parameterization of the design domain. In most methods, the design variable in the domain is given explicitly as values from 0 to 1 where 0 would represent the absence of material and 1 represents the presence of material. For boundary variation methods, the structural boundaries are implicitly defined as the contour line of a field which is a function of the design variable. Boundary variation methods are currently dominated by two methods: the level-set method and the phase-field methods. Both of these methods produce results in the design domain with crisp and smooth edges that require little post-processing effort to realize the relevant structural features. Additionally, these methods are fundamentally different from shape optimization techniques because they allow both the movement of the structural boundary and topological changes (e.g. formation, disappearance and merging of void regions).

2.3.1. Level-set method

Level sets for moving interface problems in physics were first developed by Osher and Sethian [21], with the fundamental goal of tracking the motion of curves and surfaces. This method has been applied in a wide variety of research areas [22, 23] including topology optimization. The level-set method was first applied to topology optimization in the early 2000s by Sethian and Wiegmann [24], where it was used to represent the free boundary of a structure for linearly elastic problems in structural design. In another direction, Osher and Santosa [25], in about the same time, combined level sets with a shape sensitivity analysis framework for the optimization of structural frequencies.

In the level-set method, the boundaries of the structure are represented on the zero-level curve (or contour) of the scalar function Φ which is consequently called the level-set function. Topological functions such as the changes in the boundary, merging of boundaries and formation of new voids are performed on the level-set function. The geometric boundary shape is modified by controlling the motion of the level set according to the physical problem and optimization conditions [26]. It is worth noting that most level-set formulations still rely on finite elements despite the smooth boundary representation. Thus, boundaries are still represented by discretized mesh which leads to some unsmooth results. Alternative techniques such as the extended finite elements (XFEMs) [27] have been utilized to represent the geometry in the analysis of the model which produces superior, smooth and continuous boundaries.

The level-set method does not exhibit intermediate material densities since the presence or absence of material on the domain is determined at the zero-level-set function. However, current level-set methods are known for their dependency on the initial design and locations of the level-set functions. This drawback poses a severe problem in the acceptability of solutions of level-set functions but new developments have been made to address and improve this deficiency [28]. Also, at some cases the level-set method might require re-initialization during the process when the level-set function becomes too flat or too steep. This adds computational complexity and additional tuning parameters to the algorithms which is undesirable especially for implementation with commercially available software. A visualization of this concept is presented in Figure 4 (a). A MATLAB code for the level-set method is available in Ref. [29].

Figure 4.

Conceptual difference between (a) level-set method and (b) phase-field method.

2.3.2. Phase-field method

The phase-field method originates from theories developed to track and represent phase transition and phase interface phenomena in surface dynamics [30]. This method has been utilized for solid-liquid transitions, diffusion, solidification, crack propagation, multiphase flow and eventually in topology optimization [31]. In the application of these theories, a phase-field function is specified over the design domain that is composed of two phases (e.g. A and B), which are represented by two variables as a function of the phase-field function. The boundary region between phases is a continuously varying region of thin finite thickness.

In topology optimization utilizing the phase-field method [31–33], this interface region defines the structural boundary, thus separating material from void, and is modified via a dynamic evolution of the phase-field function. The primary difference between the level-set and phase-field methods is mainly due to the fact that in the phase-field method, the interface between the boundaries of the two distinct phases is not tracked throughout optimization. Whereas in the level-set method, the boundary is tracked as it moves during the optimization process. In other words, the governing equations of phase transition are solved over the complete design domain without the initial information of the phase interface location. Consequently, phase-field methods do not require the re-initialization step as do level-set functions. Its conceptual difference with the level-set method is presented in Figure 4 (b). A MATLAB code for the phase-field method is available for download by visiting the website of Ref. [31].

2.4. Density-based methods

Currently, the most widely used methods for structural topology optimization are explicit parameterizations that are broadly classified as density-based methods. Variations of this 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 optimization tools utilize a density-based method [34].

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 understood as the element material density, ρ_e. The relevant material property of each element 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 are not practical in MATLAB.

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

Figure 5.

Generalized heat-transfer problem.

$\begin{array}{ll}\rho c{T}_t=k\Delta T+q_v\hfill & \text{in }\text{Ω}\times \text{[}0,t\text{]}\hfill \\ T=T_b\hfill & \text{on }{\text{Γ}}_1\hfill \\ -k\frac{\partial T}{\partial n}=q\hfill & \text{on }{\text{Γ}}_2\hfill \\ -k\frac{\partial T}{\partial n}=h(T-T_{\infty })\hfill & \text{on }{\text{Γ}}_3\hfill \\ {T|}_{t=0}=T_0\hfill & \hfill \end{array}$E1

where ρ is the material density, c is the specific heat of material, T_t is the temperature for a particular given time in transient cases, k is the thermal conductivity and q_v 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 Γ₁, a heat flux conduction on Γ₂, and a convective condition on Γ₃. T₀ is the initial temperature at time t = 0, T_b is a temperature imposed on Γ₁, q is a heat flux boundary condition imposed on Γ₂, h is the convective heat-transfer coefficient on Γ₃ 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{array}{ll}{q}_v=k\Delta T\hfill & \text{in }\text{Ω}\hfill \\ T=T_b\hfill & \text{on }{\text{Γ}}_1\hfill \\ -k\frac{\partial T}{\partial n}=q\hfill & \text{on }{\text{Γ}}_2\hfill \end{array}$E2

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

$\underset{\text{Ω}}{\int }(k\Delta T-q_v)\nu dx=0$E3

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

$\underset{\text{Ω}}{\int }k\nabla T\cdot \nabla \nu -q_v\nu dx=\underset{{\text{Γ}}_1}{\int }k\nabla T\cdot n\nu d\text{Γ}=-\underset{{\text{Γ}}_2}{\int }qvd\text{Γ}$E4

And the weak form can be written as

$a(T,\nu )=\mathcal{l}\left(\nu \right)\forall \nu \in \tilde{\mathbf{T}}\widetilde{\mathbf{ T}}=\{T\in [H^1\left(\text{Ω}\right)\left]|T=T_b\text{on}x\in {\text{Γ}}_1\right\}$E5

where v is in ${\displaystyle \tilde{T}}$ and ${\displaystyle \tilde{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 )=\underset{\text{Ω}}{\int }k\nabla T\cdot \nabla \nu$E6

The $\mathcal{l}(\nu )$ term is called the thermal load linear form and can similarly be obtained from Eq. (4) and is given as

$l\left(\nu \right)=\underset{\text{Ω}}{\int }{\mathrm{ }q}_v\nu dx-\underset{{\text{Γ}}_2}{\int }qvd\text{Γ}$E7

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\left(\text{Ω}\right)=\underset{\text{Ω}}{\int }k\nabla T\cdot \nabla Tdx$E8

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

$\begin{array}{lll}{\displaystyle \underset{\text{Ω}\subseteq D}{\min }}\hfill & c(\text{Ω})\hfill & \hfill \\ \text{s}\text{.t}\text{.}\hfill & a(T,\nu )=\mathcal{l}(\nu ),\hfill & \text{for all }\nu \in {\displaystyle \tilde{T}}\hfill \\ \hfill & {\displaystyle {\int }_{D}dx\le V_{\max }}\hfill & \hfill \end{array}$E9

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, ρ_e, and applying a discretized optimization model, for example, FEA, the heat conduction problem is defined as

$\begin{array}{ll}{\displaystyle \underset{{\rho }_e}{\min }}\hfill & c({\rho }_e)=T^{T}KT=Q^{T}T\hfill \\ \text{s}\text{.t}\text{.}\hfill & KT=Q\hfill \\ \hfill & V_f=V({\rho }_e)/V_o\hfill \\ \hfill & 0<{\rho }_e\le 1\hfill \end{array}$E10

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_e, and the material interpolation schemes is applied here as is formally given as

$K({\rho }_e)={\displaystyle {\displaystyle \sum _{e=1}^N}{k}_{eff}({\rho }_e)k_e}$E11

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

$c({\rho }_e)=Q^{T}T,\text{where }T\text{solves: }{\displaystyle {\displaystyle \sum _{e=1}^N}{k}_{eff}({\rho }_e)k_e}T=Q$E12

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

$k_{eff}({\rho }_e)=\left(k_{\max }-k_{\min }\right){\rho }_e{}^p$E13

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 ρ_e. 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 }_e)=Q^{T}T-{\lambda }^T(KT-Q)$E14

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 }_e}=(Q^T-{\lambda }^{T}K)\frac{\partial T}{\partial {\rho }_e}-{\lambda }^T\frac{\partial K}{\partial {\rho }_e}T$E15

which can be re-written as

$\frac{\partial c}{\partial {\rho }_e}=-{\lambda }^T\frac{\partial K}{\partial {\rho }_e}T$E16

When λ satisfies the adjoint equation:

$Q^T-{\lambda }^{T}K=0$E17

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(ρ_e) for the main problem in Eq. (10), considering the SIMP interpolation as presented in Eq. (13), are

$\frac{\partial c}{\partial {\rho }_e}=-p(k_{\max }-k_{\min }){\rho }_e{}^{p-1}{T}^{T}KT$E18

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

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 heat-transfer 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.

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 level-set 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-to-point 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

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 well-defined 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 fluid-thermal interaction problems under constant input power. Kontoleontos et al. [95] published an adjoint-based constrained topology optimization for viscous flows, including heat transfer.

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 (${\displaystyle \dot{Q}}/m$). Oevelen and Baelsmans [98] demonstrated solutions to a conjugate heat-transfer problem using a two-layer-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 (${\displaystyle \dot{Q}}/V$ and ${\displaystyle \dot{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 heat-sink 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.

5. Conclusion

In this chapter, we have re-introduced topology optimization with special focus on the progress of heat exchanger design over the past two decades. We have first given an overview of its historical background in terms of structural topology optimization. We have then conceptually introduced the different methods developed over the years in topology optimization. Learning references for each of the methods mentioned, together with MATLAB codes, were cited and is expected to help those who are interested in further learning and investigating topology optimization. A chronological review highlighting the different progress over the years related to heat exchanger design was also given.

Novel heat-transfer structures are still being realized to further drive design performance to its limits. Topology optimization, as a physics-based and automated layout optimization method, will indeed serve as a valuable design tool for heat-transfer systems. Heat exchanger designs arising from topology optimization has now been realized and continuous efforts are still being made to further improve both methods and implementation. Topology optimization is expected to play a bigger role in the coming years for heat exchanger design.

Nomenclature Abbreviations

2D/3D

Two-dimensional/three-dimensional

BEM

Boundary element method

BESO

Bi-directional evolutionary structural optimization

DHTPC

Dissipation of heat-transfer potential capacity

DOF

Degrees of freedom

DSA

Design sensitivity analysis

ECP

Element connectivity parameterization

ESO

Evolutionary structural optimization

FEA

Finite element analysis

FEM

Finite element method

FVM

Finite volume method

IBM

Immersed boundary method

MMA

Method of moving asymptotes

MOO

Multi-objective optimization

OC

Optimality criteria

RAMP

Rational approximation of material properties

SIMP

Solid isotropic material with penalization

SLP

Sequential linear programming

SQP

Sequential quadratic programming

TO

Topology optimization

VP

Volume-to-point problem

Symbols and variables (in the order of appearance)

μ

Dimension control parameter in homogenization method

θ

Orientation control parameter in homogenization method

ρ

Density variable for density-based topology optimization, material density

c

Standard variable for the design objective or compliance, material-specific heat

T

Temperature

k

Thermal conductivity

q

Heat flux

n

Boundary normal vector

h

Convection heat-transfer rate

Ω

Domain

t

Time

Γ

Boundaries

v

Virtual temperature field

a

Energy bi-linear form

l

Thermal load linear form

x

Design variable

V

Volume

T˜

Test temperature vector

H

Sobolev space

T

Temperature vector in heat-transfer TO

K

Global stiffness matrix for FEA

Q

Applied thermal load in heat-transfer TO

k

Global stiffness matrix for FEA

λ

Langrangian multiplier

Subscripts

0, 1, 2, .., i

Standard discrete numerical counter

e

Element in discretization

mat

Material

min

Minimize/minimummax

t

Time, to imply transient case in derivation

v

Per unit volume, in derivation

b

Imposed boundary condition, w/temperature in derivation

eff

Effective, used with thermal conductivity, k

p

Penalty parameter for SIMP