**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**

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 optimiza-

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

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.

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

tion itself.

62 Heat Exchangers– Design, Experiment and Simulation

are presented in **Figure 1**.

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 (*μ*<sup>1</sup> , *μ*<sup>2</sup> 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

**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].

anisotropic macro-scale representation of the material can be achieved such as pure solid, pure void, composite and porous material.

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].

### **2.3. Boundary variation methods**

anisotropic macro-scale representation of the material can be achieved such as pure solid,

**Figure 2.** (a) Schematic of the treatment of design variables contained in elements showing two possible microstructures.

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

'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

pure void, composite and porous material.

64 Heat Exchangers– Design, Experiment and Simulation

(b) Result for a truss problem reflecting generated microstructures [14].

is called 'inverse homogenization'.

**2.2. 'Hard-kill' 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

**Figure 3.** Conceptual differences of different 'hard-kill' methods (a) ESO, (b) combinatorial techniques and (c) topological derivative.

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].

#### *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

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

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 phasefield 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, phasefield 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**

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

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

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].

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

and topological changes (e.g. formation, disappearance and merging of void regions).

*2.3.1. Level-set method*

*2.3.2. Phase-field method*

optimization of structural frequencies.

66 Heat Exchangers– Design, Experiment and Simulation

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, *ρ<sup>e</sup>* . 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.
