Multi-objective Bayesian topology optimization of a lattice-structured heat sink in natural convection

Additive manufacturing (AM) has an affinity with topology optimization to think of various designs with complex structures. Hence, this paper aims to optimize the design of a lattice-structured heat sink, which can be manufactured by AM. The design objectives are to maximize the thermal performance of convective heat transfer in natural convection simulated by computational fluid dynamics (CFD) and to minimize the material cost required for AM process at the same time. The lattice structure is represented as a node/edge system via graph theory with a moderate number of design variables. Bayesian optimization, which employs the non-dominated sorting genetic algorithm II and the Kriging surrogate model, is conducted to search for better designs with the minimum CFD cost. The present topology optimization successfully finds better lattice-structured heat sink designs than a reference fin-structured design regarding thermal performance and material cost. Also, several optimized lattice-structured designs outperform reference pin-fin-structured designs regarding thermal performance though the pin-fin structure is still advantageous for a material cost-oriented design. This paper also discusses the flow mechanism observed in the heat sink to explain how the optimized heat sink structure satisfies the competing design objectives simultaneously.


Introduction
In recent years, additive manufacturing (AM) using a threedimensional printer has attracted attention and has been utilized in various fields such as mechanical engineering, aerospace engineering, biomechanical medical engineering, and architecture (Gardan and Schneider 2015;Gao et al. 2015). According to Pérez-Pérez et al. (2018), the main advantages of the additive fabrication concept used in AM are as follows.
1. It is possible to develop and manufacture complex structures, which has been difficult or impossible by conventional numerical control manufacturing such as cutting, drilling, and milling.
2. It is effective to save materials since they are not processed by removal or deformation as in the conventional methods. 3. It is unnecessary to change the tools for modeling since the structure is modeled by hardening material powder (e.g., resin and metal) with a heat source. 4. It is easy and quick to start the process of modeling once the shape data is ready.
This paper wishes to take AM process's advantages to design a heat sink, which transfers thermal energy from higher-temperature devices (e.g., CPU and GPU) to lower-temperature fluid mediums (e.g., air and liquid coolant). An issue about heating is getting critical as CPU/GPU performance grows up (Garimella et al. 2008). The AM process is expected to sophisticate the heat sink structure design and keep up with further development in high-performance computing.
The AM process has an affinity with topology optimization, which has been gaining popularity in the recent field of engineering design, as reviewed in Liu et al. (2018). Topology optimization optimizes material layout within a given design domain. It enables designers to think of various 1 Page 2 of 15 structures to be employed for a new product, such that the topological features (e.g., continuity and connectivity) can be variant. Such structures can hardly be considered by classical sizing optimization and shape optimization due to insufficient degrees of variety.
Topology optimization needs to choose a method to describe material distribution ( ) where is the position vector in the design domain. ( ) = 0 and ( ) = 1 correspond to the absence and presence of material, respectively. Topology optimization began with the homogenization method (Bendsøe and Kikuchi 1988), which models material layouts with periodic microstructures. This method allows grayscale material ( 0 < ( ) < 1 ), which leads to a vague outline of the structure. Then, the homogenization method was subsequently refined into the density method (Bendsøe 1989), which defines material density in the range of 0 ≤ ( ) ≤ 1 while penalizing the intermediate values of ( ) to approach a black ( ( ) = 0 ) or white ( ( ) = 1 ) design. When introducing penalization, however, an unrealistic material layout with checkerboard patterns appears. Such patterns are commonly removed through regularization (Sigmund and Petersson 1998). For grayscale-free topology optimization, a level-set function method (Sethian and Wiegmann 2000) has been proposed. The level-set function defines the absence and presence of material as ( ) < 0 and ( ) > 0 , respectively, and represents the outline of structure with the zero contour ( ) = 0.
In topology optimization, the material distribution ( ) is parameterized with n design variables = x 1 , x 2 , … , x n . The density method discretizes the design domain into n elements, which are located at (1) , (2) , … , (n) , and assigns material density to the ith element as the ith design variable, i.e., x i = ( (i) ) ( i = 1, 2, … , n ). Similarly, the level-set method collocates n control points (located at (1) , (2) , … , (n) ) in the design domain, and assigns the value of the level-set function to the ith control point as the ith design variable (i.e., x i = ( (i) ) ). The level-set function is then constructed in the whole design domain by interpolating the n control points. are continuous variables in both methods to allow for gradient-based optimizers because a discrete problem is tough to solve. In addition, n should be huge for high-resolution topology representation, and thus topology optimization usually results in many-variable optimization. Thus, a gradient-based optimizer, which is often assisted with adjoint-based sensitivity analysis, has been used for many-variable topology optimization. However, this optimizer tends to get stuck to a local optimal solution rather than the global optimal solution. In particular, the obtained optimal solution strongly depends on initial settings in a thermal-fluid problem where the solution is evaluated by solving nonlinear governing equations (Yaji et al. 2015).
Many works of literature studied the topology optimization of heat sink design. The design objective was to minimize the thermal compliance (Alexandersen et al. 2014(Alexandersen et al. , 2016Joo et al. 2017;Lohan et al. 2017), to maximize the total potential energy (Iga et al. 2009), to minimize the thermal resistance (Dede et al. 2015;Haertel et al. 2018), to minimize the temperature (Martínez-Maradiaga et al. 2019;Lundgren et al. 2019), etc. In this way, heat sink design was usually optimized in thermal performance (i.e., single-objective optimization). Besides, most of the works considered a volume constraint on the amount of material use (i.e., material cost) specified by practitioners. From the viewpoint of AM process, it is essential to study a tradeoff between the thermal performance and the material cost. The single-objective optimization explores a single optimal design corresponding to a different allowable volume, which means that a gradient-based optimizer needs to be run many times for the trade-off study while changing the allowable volume. Also, the allowable volume may not be deterministic and may change depending on practitioners' preference or AM machine capability.
Instead, for an efficient trade-off study, reducing the material cost also needs to be considered as an additional design objective rather than a constraint. As in Joo et al. (2017), minimizing the thermal compliance sometimes reduces the material cost simultaneously; however, the thermal performance and the material cost are the design objectives that usually compete with each other (i.e., multi-objective optimization). The multi-objective optimization explores not single, but many optimal designs (so-called non-dominated solutions or Pareto-optimal solutions), which exist between competing design objectives and are helpful for the tradeoff study, without specifying any constraint limit. Dbouk (2017) presented a review of topology optimization methods to design optimal heat transfer systems. He claimed that topology optimization is not yet a robust numerical design technique for finding optimal designs of thermal systems; e.g., optimization algorithms and multi-objective functions libraries need to be studied still more. Besides, Alexandersen and Andreasen (2020) provided an overview of the literature for topology optimization in fluid-based problems to suggest a future perspective on more complex applications (e.g., transient and turbulent flow). Such applications will be highly nonlinear optimization problems and are not easy to solve using gradient-based optimizers, which have been considered promising to solve the conventional single-objective topology optimization problems.
Therefore, this paper aims to perform multi-objective topology optimization of a heat sink designed in natural convection, which considers a balance between thermal performance improvement and material cost reduction. A gradientfree population-based optimizer, such as the non-dominated sorting genetic algorithm II (NSGA-II) (Deb et al. 2002), is employed for the present multi-objective optimization. The population-based optimizer explores the optimal solution by evolving a population with many solutions, whereas the gradient-based optimizer mathematically explores the optimal solution by starting from a single initial solution given by the user. Therefore, the population-based optimizer has the efficient capability to explore the set of non-dominated solutions at the same time though these solutions are not necessarily guaranteed global optimal due to the nature of random exploration. On the other hand, gradient-based optimizers have already been applied to multi-objective topology optimization in thermo-fluid problems. However, the gradient-based optimizers are almost limited to use with varying weights for density-method-based topology optimization (Marck et al. 2012;Li et al. 2019;Dong and Liu 2020) and with adaptive weights for level-set-method-based topology optimization (Sato et al. 2017). These gradientbased optimizers are still restricted to convex Pareto-optimal fronts, while the population-based optimizer does not have such the restriction.
This paper also considers a lattice-structured heat sink parameterized with a moderate number of design variables (currently, n = 68 ) based on graph theory (Bender and Williamson 2010). This theory is similar to the ground structure method proposed by Dorn et al. (1964) for optimizing a truss layout with a variable truss cross-sectional area. This paper employs the graph theory to discretely structure the heat sink with a fixed truss cross-sectional area in a design domain. This paper also employs commercial software for computational fluid dynamics (CFD) simulation to evaluate the thermal performance for convenience. Therefore, we cannot compute gradients based on an adjoint method due to the limitation of the commercial software and need to use a gradient-free population-based optimizer instead.
The population-based optimizer is often computationally expensive because it needs to evaluate objective functions for many solutions in a population. It is more critical in real-world design optimization where evaluation is computationally expensive even for one solution, such as in the present study that performs CFD. Hence, this paper employs a surrogate model to assist the NSGA-II for efficient multi-objective topology optimization with a less computational cost. The surrogate model is constructed as a regression or interpolation function of the sample dataset, in which the true objective function is evaluated at several locations in the design variable space. The surrogate model can instantaneously estimate the objective function value at an arbitrary unsampled location in the design space. This study employs the Kriging surrogate model (Sacks et al. 1989) based on Bayesian statistics because it is suitable for approximating multi-modal functions and quantifying the model uncertainty.
The organization of this paper is as follows. Section 2 defines the optimization problem to design the lattice topology of a heat sink. Then, Sect. 3 describes a numerical method to implement multi-objective Bayesian topology optimization. Section 4 shows and discusses optimized designs from a scientific point of view. Finally, Sect. 5 concludes this paper.

Design domain
One of the present reference designs is a commercially versatile heat sink (named "21F50") shown in Fig. 1. This is a fin-structured heat sink and conventionally manufactured by extrusion. The size is 21 mm height (15 mm fin + 6 mm base) × 50 mm width × 50 mm length. The upper fin part (15 mm × 50 mm × 50 mm, shaded in red in Fig. 1a) is set as the design domain while the lower base part (6 mm × 50 mm × 50 mm, shaded in violet in Fig. 1a) is fixed.
Besides, a pin-fin heat sink shown in Fig. 1b is another reference design. This design has the same height, width, and length as 21F50 and has 7 × 7 pins (diameter = 2.6 (mm)), consistent with the node/edge system via the graph theory employed in this paper explained in Sect. 2.3 later. The pin-fin heat sink is expected to perform better than 21F50 in natural convection when placed flat. Figure 2 shows how to represent a lattice-structured heat sink to be optimized. The lattice is replaced with the edge network constructed with 3 × 7 × 7 = 147 nodes collocated with an 8.3 mm interval in the design domain, which are denoted by the black dots in Fig. 2a.

Topology representation
Considering that a heat sink is placed flat, which will be explained in Sect. 2.4, the design domain can be reduced to 1/8 (hatched in red in Fig. 2b) by considering flow symmetry. In the reduced design domain, 3 × 10 = 30 nodes are divided into two layers to be optimized separately: lower layer denoted by the green spheres and upper layer denoted by the blue circles in Fig. 2a. Each layer consists of 2 × 10 = 20 nodes and, as denoted by the black lines in Fig. 2c, linking all these nodes constructs 20 C 2 = 190 edges. However, most of the edges cannot be manufactured in an AM process if they form an angle of higher than 45° with the build direction. Thus, 34 manufacturable edges (green lines in Fig. 2d) remains in each layer.
Consequently, 34 × 2 = 68 edges are considered in the two layers in total. Each edge is then replaced with a solid cylinder with a diameter = 2.6 (mm). The lattice structure is obtained by cutting off the protruding cylinders from the design domain and attaching them to the base part, as shown in Fig. 2e.

Design variables via graph theory
As stated in Sect. 2.2, the lattice structure is represented with a node/edge system. Therefore, the present topology optimization parameterizes the node/edge system through graph theory (Bender and Williamson 2010). In general, this theory quantifies the presence or absence of each edge linking the jth and kth nodes using an adjacent matrix with the (j, k) entry a jk . Note that this is a symmetric matrix satisfying a jk = a kj and the diagonal entries a jj do not represent edges. For the present lattice-structured heat sink, the adjacent matrix is sized 30 × 30 and contains 68 valid entries corresponding to manufacturable edges while the others correspond to non-manufacturable edges. The present topology optimization treats each effective entry of the adjacent matrix as the ith design variable Yoshimura et al. (2019), who used the graph theory to optimize the groove structure of a flow micromixer. The present design variables need to be continuous for a surrogate model used in Bayesian optimization. Another reason to bias the cut-off of x i toward the lower value ( = 0.2) comes from our intention to consider design candidates as various as possible during the present optimization. A higher cut-off (e.g., > 0.5) is expected to increase the possibility of an empty lattice structure or a lattice-free structure, which is useless for design optimization. Consequently, the present topology optimization considers 68 design variables in total.

Thermal performance
The first objective function f 1 ( ) to be considered here is to maximize a heat sink's thermal performance in natural convection. The total heat transfer rate at the interface between the heat sink and the surrounding air, Q, is the present quantity of f 1 ( ).
Q is evaluated by CFD simulation using commercial software, ANSYS Fluent 2019 R1.2 (ANSYS, Inc 2018) with taking into account convective heat transfer on the surface of a heat sink. In this study, the radiative heat transfer is not considered. The steady-state incompressible Navier-Stokes equations (continuity, momentum conservation, and energy conservation) for the laminar air (density = 1.225 (kg/m 3 ), specific heat capacity c p = 1006.43 (J/kg K), thermal conductivity k = 0.0242 (W/m K), viscosity =1.7894 × 10 −5 (kg/m s)) with gravity | | = 9.81 (m/s 2 ) are given by Eqs. 1-3 and solved in the simulation domain shown in Fig. 3. where is the velocity vector, p is the pressure, and T is the temperature to be calculated. In Eq. 2, the Boussinesq approximation is applied to the natural convection flow, with the thermal expansion coefficient = 0.00333 (1/K) and the operating temperature T 0 = 288.16 (K). The heat sink is made of an aluminum alloy, AlSi10Mg (density = 2719 (kg/m 3 ), specific heat capacity c p = 871 (J/kg K), thermal conductivity k = 190.5 (W/m K)), which is a typical material used for metal-based AM, and placed flat at the bottom of the simulation domain. The temperature field in the solid heat sink is simulated by solving Eq. 3 with = . The numerical schemes chosen in Fluent are listed in Table 1.
The simulation domain is discretized into mesh grids with the minimum element size of 0.2 mm by the cut-cell method. It enables automatic mesh generation for any complex  structure considered in topology optimization. Also, the mesh grids perfectly fit wall surfaces. Hence, the present CFD is expected to yield a reasonable solution satisfying the flow conservation laws.
Regarding the boundary conditions, the bottom surface of the heat sink (colored in red in Fig. 3) is heated at constant temperature = 323.15 (K). The outer boundaries of the simulation domain (colored in blue in Fig. 3) are set as the thermally insulating wall (density = 50 (kg/m 3 ), specific heat capacity c p = 800 (J/kg K), thermal conductivity k = 0.09 (W/m K)), where convective heat transfer occurs between the interior (temperature T in to be calculated here) and the exterior (ambient temperature T ex = 298.15 (K), heat transfer coefficient h = 5 (W/m 2 K)) of the domain through Newton's law of cooling as where ∇ ⟂ T in stands for the temperature gradient in the direction normal to the insulating wall on the interior side. Thermal conditions are coupled at the interface between the heat sink (solid) and the air (fluid).
The present CFD simulation is parallelized with 16 CPUs in the supercomputer system AFI-NITY owned by the Institute of Fluid Science, Tohoku University, Japan. It takes 25-40 min to evaluate Q for each design.

Material cost
The second objective function f 2 ( ) is to minimize the heat sink's material cost. This study considers the material volume occupying the design domain, V, as the cost measure. For more straightforward evaluation, the present topology optimization calculates the total length of the present edges (i.e., the green lines in Fig. 2d), L, instead of V during the optimization process. Note that L can be analytically evaluated by neither CFD nor surrogate model (Sect. 3.2) as soon as the values of the design variables are determined. We confirm in Fig. 4 that L is almost linearly correlated with V though L is sometimes overestimated due to the overlap of the cylinders, which means that L can be substituted for V.

Constraint function
Heat conduction must occur over the whole lattice structure. It means that, as illustrated in Fig. 5, all the present edges must be linked to the base part either directly or indirectly, i.e., there is no floating edge. Therefore, members that are not directly supported underneath are accepted if they are attached to bars that do reach the base plate. To ensure the connectivity, the present topology optimization defines the constraint function as g 1 ( ) = 0 and g 1 ( ) = 1 for feasible and infeasible solutions, respectively. The connectivity can be analytically checked by sweeping through all entries in the adjacent matrix.

Bayesian optimization method
One of the present objective functions, f 1 ( ) = Q , is evaluated by expensive CFD simulation. Hence, Bayesian optimization is performed by replacing the true function f 1 ( ) with the estimate f 1 ( ) that can be inexpensively calculated by the Kriging surrogate model. The flowchart of the present optimization is shown in Fig. 6, and the details will be explained in the following sections.

Generate initial samples
This step generates a dataset of N initial samples in the design variable space, each of which corresponds to a different combination of n design variables = x 1 , x 2 , … , x n (currently, a different lattice structure represented with n = 68 ). This study employs the Latin hypercube sampling (LHS) (McKay et al. 1979), which generates samples orthogonal to each other in the design variable space to comprehend the whole design space even with a smaller sample size than the Monte Carlo sampling. Then, true objective functions are evaluated (currently, f 1 ( ) = Q is evaluated by CFD) for each initial sample.
In the present topology optimization, 100 LHS points are generated in the design space, and 98 of them are available for the initial samples associated with the true f 1 ( ) values

Construct the Kriging surrogate model
Using the initial samples, this step constructs a surrogate model that approximates an unknown objective function with a known algebraic function. The Kriging surrogate model (Sacks et al. 1989) employed here is a stochastic process that realizes N given samples as illustrated in Fig. 7. For an objective function f ( ) to be maximized, the stochastic process yields the mean f ( ) and the standard deviation ŝ( ) to model the estimate of f ( ) and the uncertainty (i.e., error bar), respectively. Note that ŝ( ) = 0 at a sampled location. Model uncertainty ŝ( ) is not available in other surrogate models (e.g., polynomial regression, radial basis function (Powell 1987)) and is useful for efficient surrogate-based global optimization, which will be explained in Sect. 3.3.
In the present topology optimization, only f 1 ( ) = Q is replaced with the Kriging surrogate model during the optimization process while f 2 ( ) = L is analytically evaluated without the surrogate model.

Explore non-dominated solutions
This step explores non-dominated solutions on the surrogate model. Now the NSGA-II (Deb et al. 2002), which is the most prevalent genetic algorithm worldwide, is employed as the optimizer. The NSGA-II performs non-dominated sorting and crowding distance sorting to search for the solutions that seem located close to and distributed widely over the true Pareto-optimal front in multi-objective optimization.
The NSGA-II is now assisted with the Kriging surrogate model as follows. For an objective function f ( ) to be maximized, the Kriging surrogate model can calculate the expected improvement EI[f ( )] , which is illustrated in Fig. 7 and defined as where f ref is a reference value, and F ∼ Norm[f ( ),ŝ 2 ( )] . EI[f ( )] indicates how much better a solution will get than the reference solution under the model uncertainty ŝ( ) .
Instead of maximizing f ( ) , therefore, maximizing EI[f ( )] is expected to lead to an optimal solution on the surrogate model. In addition, this solution is suitable as an additional sample, which will be explained in Sect. 3.5.
In the present optimization, the NSGA-II is executed with the population size of 200, the number of generations of 200, the mutation rate of 1/n (currently, n = 68 ), the  distribution index for crossover of 30, and the distribution index for mutation of 20. Assisted with the Kriging surrogate model, the NSGA-II finds the non-dominated solutions that maximize EI[f 1 ( )] ( f 1 ( ) = Q , and f ref is set to Q evaluated for the reference design 21F50 (Fig. 1a)) and minimize f 2 ( ) = L subject to g 1 ( ) ≤ 0.

Choose representative solutions
This step chooses several representatives from many nondominated solutions explored by the NSGA-II. First, the extreme (i.e., EI[f 1 ( )]-maximum and f 2 ( )-minimum) solutions are chosen as depicted by the closed circles in Fig. 8a. In addition, this study performs K-means clustering (Jain et al. 1999) to choose K solutions, which correspond to the centroids of K clusters (the closed circles in Fig. 8b). The present topology optimization sets K = 3 , hence this step chooses at most K + 2 = 5 representatives.

Generate and add new samples
For corresponding to each representative solution, the true objective function f 1 ( ) = Q is evaluated by CFD. The representative solutions associated with f 1 ( ) are then used as new samples added to the current sample dataset (i.e., N ← N + K + 2 ). Finally, the optimization process goes back to the previous step to update the Kriging surrogate model (Sect. 3.2), including the additional samples evaluated by new CFD simulations, and is iterated until the additional samples show convergence in the objective function space. Note that "update" corresponds to the outer loop for the Kriging surrogate model as written in Fig. 6, while "generation" corresponds to the inner loop for the NSGA-II optimization.
As mentioned above, optimization assisted with the Kriging surrogate model proceeds in a data-driven manner, which expects additional samples to reach global optima and improve the surrogate model accuracy at the same time. This data-driven optimization is often called Bayesian optimization.
Regarding computational time, the present Bayesian optimization takes around 3 h for Kriging surrogate model construction + 5 min for NSGA-II optimization + 25-40 min for new sample generation every update. On the other hand, if CFD evaluates all solutions without the surrogate model, the NSGA-II optimization will take 100 h, even with 100% parallel processing. Thus, using the surrogate model can reduce the total computational time. Figure 9 plots the samples generated through the present topology optimization in the objective function space. Here note that the vertical axis represents V instead of f 2 ( ) = L . The green plots are the 98 initial samples generated by the LHS. The orange plots are 69 additional samples obtained by two-objective Bayesian optimization (maximizing EI[f 1 ( )] ( f 1 ( ) = Q ) and minimizing f 2 ( ) = L ) through 15 updates (i.e., CFD is conducted for 98 + 69 = 167 samples in total in the two-objective Bayesian optimization). For comparison, single-objective Bayesian optimization (maximizing EI[f 1 ( )] only) is also performed, and 18 additional samples generated through 18 updates are plotted in blue (i.e., CFD is conducted for 98 + 18 = 116 samples in total in the singleobjective Bayesian optimization). Among all these samples, the non-dominated ones are marked with red dots. Moreover, several reference data, 21F50 (Fig. 1a), pin-fin ( 2.6 ) (Fig. 1b) and other variants with different pin diameters ( 1.4 , 2.0 , 3.2 , 3.8 , and 4.4 ), and full lattice (Fig. 2e) are plotted.

Results and discussion
In Fig. 9, the two-objective optimization can obtain better additional samples than the initial samples in terms of both Q and V. For example, the sample with maximum Q,   Fig. 9 Samples plotted in the objective function space which is labeled "11-005 (2-obj)" meaning the fifth sample generated at the eleventh update in the two-objective optimization, achieves 26% improvement in Q and 61% reduction in V compared to . The single-objective optimization leads to the Q-maximum solution, 14-001 (1-obj), which is relatively comparable to 11-005 (2-obj) regarding Q but worse regarding V. It also indicates that the full lattice design is not worthy of being adopted as a heat sink due to the low thermal performance and high material cost. Figure 9 also shows that the two-objective optimization obtains the non-dominated set of the lattice-structured designs in a trade-off relationship between Q maximization and V minimization. Another reference design, pin-fin ( 2.6 ), is almost located in the non-dominated lattice-structured designs with the same edge diameter ( 2.6 ). The family of pin-fin designs makes a balance between Q and V with different pin diameters in the range = 1.4 -3.8 (mm). Pinfin ( 3.8 ) is the best Q performer in this family; however, 11-005 (2obj) still outperforms pin-fin ( 3.8 ) regarding Q by 8% while keeping the similar material cost V ≃ 6.0 × 10 −6 (m 3 ). Comparing the Pareto-optimal fronts between the lattice-structured designs and the pin-fin designs indicates that, as written in Fig. 9, a lattice structure is advantageous for a performance-oriented heat sink design, while a pin-fin structure is advantageous for a cost-oriented design.
Besides, Fig. 10 shows the history of the hypervolume indicator (Zitzler and Thiele 1998) of the samples calculated every update in the current two-objective Bayesian optimization, where a zero hypervolume corresponds to 21F50 in the Q and V space. It indicates that 15 updates significantly increase the hypervolume and lead to convergence after 11 updates; i.e., the present Bayesian optimization improves both Q and V of the heat sink designs simultaneously.
Next, this paper visualizes the structure and the flow fields for two reference designs: 21F50 and pin-fin ( 2.6 ) in Figs. 11 and 12, respectively, and three heat sinks designed by the two-objective optimization: 11-005 (2-obj), 11-004 (2-obj), and 10-001 (2-obj) in Figs. 13, 14, and 15, respectively. These figures visualize the temperature fields with colored contours and the velocity fields with white arrows.
In 21F50 (Fig. 11), the surrounding cool air flows from two sides into the gaps between the fins. After the air is heated at the core of the heat sink, it flows out upward. The pin-fin ( 2.6 ) design (Fig. 12), on the other hand, can take the surrounding cool air from all four sides into the heat sink core. This design is more beneficial to improve the thermal performance in natural convection compared to 21F50. Similarly, in the optimized designs (Figs. 13, 14, 15), the surrounding cool air flows from all four sides into the heat sink core due to the symmetric structure (Fig. 2b). It means that thermal performance is not affected by the orientation of a heat sink, i.e., the present optimized designs can keep the performance robust.
Comparison of the optimized designs (Figs. 13,14,15) suggests the following common features in the lattice structures.
1. In the lower layer (corresponding to the green shade in Fig. 2e), intakes exist at all the sides. 2. In the upper layer (corresponding to the blue shade in Fig. 2e), the lattice is dense around the core.
The intakes are essential for getting a sufficient amount of cool air from the surroundings into the heat sink core. The upward natural convection drives the flow at the core. The intakes are also favorable to reduce the lattice volume. With limited lattice volume, it is more effective to enhance heat transfer between the cool air and the heated lattice more locally around the core, where the cool air converges and passes. Note that, in 21F50 (Fig. 11), cool air cannot reach the core, which indicates that the heat sink does not work around the core. In addition, the optimized designs have a wider intake than the pin-fin ( 2.6 ) design (Fig. 12), i.e., this feature is the primary advantage of the present latticestructured design over the conventional pin-fin design to achieve better thermal performance in natural convection. Thus, considering the above features is essential to create a balanced design of a lattice-structured heat sink between the thermal performance and the material cost.

Conclusions
This paper developed a multi-objective Bayesian topology design optimization method to maximize the thermal performance (total heat transfer rate) and minimize the material cost (volume) of a lattice-structured heat sink, which can be manufactured by AM, in natural convection. The lattice structure was represented as a node/edge system via graph theory with a moderate number of design variables. The NSGA-II was employed as the present multi-objective global optimizer. Since the thermal performance needs to be evaluated by expensive CFD simulation, the Kriging surrogate model was also employed together with the NSGA-II. The Kriging-surrogate-based optimization is useful to explore non-dominated solutions under the model uncertainty with a less computational cost. The present multi-objective Bayesian optimization successfully found new lattice-structured heat sink designs better than the reference fin-structured design in terms of both the thermal performance and the material cost. The present optimization also indicated that the optimized lattice-structured designs outperform the reference pin-fin designs regarding thermal performance, which are conventionally considered efficient in natural convection, though the pin-fin designs are still advantageous for material cost reduction. Moreover, as useful design knowledge, the present optimization suggested characteristic features in a lattice structure to improve the thermal performance traded with limited lattice volume and revealed a characteristic flow mechanism to explain these features. Acknowledgements The authors are grateful for using the supercomputer system "AFI-NITY" owned by the Institute of Fluid Science, Tohoku University, Japan. They also would like to appreciate the collaboration with NTT Data Engineering Systems Co., Ltd., Japan, on AM technology.

Conflict of interest
On behalf of all authors, the corresponding author states that there is no conflict of interest.

Replication of results
Essential details for replication of results have been described in the manuscript. Except for ANSYS Fluent used for CFD, the C codes are made open sources and available upon request.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.