International Journal of Heat and Mass Transfer 127 (2018) 32–40 Contents lists available at ScienceDirect International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt Efficient three-dimensional topology optimization of heat sinks in natural convection using the shape-dependent convection model Younghwan Joo, Ikjin Lee, Sung Jin Kim ⇑ Department of Mechanical Engineering, Korea Advanced Institute of Science and Technology, 291 Daehak-ro, Yuseong-gu, Daejeon 34141, Republic of Korea a r t i c l e i n f o Article history: Received 27 June 2018 Received in revised form 3 August 2018 Accepted 3 August 2018 Keywords: Three-dimensional topology optimization Shape-dependent heat transfer coefficient Heat sink Natural convection Constructal design a b s t r a c t In this study, heat sinks in natural convection are thermally optimized using the method of threedimensional topology optimization. In order to perform three-dimensional topology optimization with low computational cost, a shape-dependent convection model is proposed. This model accounts for the variation of the heat transfer coefficient depending not only on the local shape of the fins but also on the development of the thermal boundary layer. The physical validity of the proposed model is confirmed by the fin geometry of the topology-optimized design that matches the multiscale structures proposed previously by the constructal theory. For further validation, the effective heat transfer coefficient evaluated by the proposed model is compared to that obtained from numerical simulations. Because the new topology-optimized design has a complicated fin geometry, design simplification is performed to yield a more manufacturable design. The thermal performance of the topology-optimized heat sink is compared to that of the radial plate-fin heat sink optimized analytically using an existing correlation. It is found that the topology-optimized heat sink has 13% lower thermal resistance and 48% less mass than the optimized radial plate-fin heat sink. This implies that the three-dimensional topology optimization method suggested in this study can provide a heat sink design with improved thermal performance and reduced mass for various practical applications. Ó 2018 Elsevier Ltd. All rights reserved. 1. Introduction Natural convection heat sinks are widely used because of their simplicity and reliability [1,2]. Thermal engineers are interested in finding the optimum geometry of a heat sink that dissipates the maximum amount of heat under a given physical volume. Conventionally, heat sinks with a given fin shape such as plates or pins have been optimized experimentally, numerically, or analytically [3–7]. Apart from conventional methods, topology optimization, a method commonly used in structural problems , has been recently applied to the thermal optimization of heat sinks [9–13]. Since this method does not impose any constraint on the fin shape, completely new designs for heat sinks can be suggested. To optimize a heat sink using the topology optimization method, the convective heat transfer rate from the heat sink surface needs to be accurately evaluated. Alexandersen et al.  solved full conjugate heat transfer problems numerically in a 2-D computational domain to account for the local variation of the heat transfer coefficient, and they extended the method to a more complex 3-D problem . Recently, they applied their 3-D topology ⇑ Corresponding author. E-mail address: email@example.com (S.J. Kim). https://doi.org/10.1016/j.ijheatmasstransfer.2018.08.009 0017-9310/Ó 2018 Elsevier Ltd. All rights reserved. optimization method to a passive cooler for light-emitting diode lamps, and obtained new designs in the given physical domain . Their method generally requires high computational cost because the velocity field and the corresponding distribution of the heat transfer coefficient need to be calculated each time the shape of the heat sink fins changes at every iteration. The supercomputing equipment had to be utilized to cope with a significantly increased computational load caused by the 3-D full conjugate heat transfer problem. Hence, for the practical use of 3-D topology optimization, it is necessary to develop a simple and efficient method that enables 3-D topology optimization with fairly accurate predictions of the shape-dependent heat transfer coefficient. In our previous work , the shape-dependent convection model for topology optimization of heat sinks in natural convection was proposed for a 2-D computational domain. In this model, a geometric parameter that reflects the local shape of the fins was obtained, and then the heat transfer coefficient was evaluated by substituting this geometric parameter into the existing Nusselt number correlation. Since it was not necessary for this model to solve the full conjugate heat transfer problem, the computational load could be greatly reduced. Even if this model is extended to 3-D topology optimization, the computational load for this Y. Joo et al. / International Journal of Heat and Mass Transfer 127 (2018) 32–40 33 Nomenclature A cp d El f g H h k L Ne Nu Pr Q R Rth Sh Sv T1 Tb wc xi surface area [m2] specific heat [kJ/kg-K] fin diameter [m] Elenbaas number [–] volume fraction [–] standard acceleration of gravity [m/s2] fin height [m] convective heat transfer coefficient [W/m2-K] thermal conductivity [W/m-K] heat sink length [m] adjacent elements [–] Nusselt number [–] Prantl number [–] total heat input [W] outer radius of domain [m] thermal resistance [K/W] horizontal spacing [m] vertical spacing [m] ambient temperature [K] base temperature [K] channel spacing [m] spatial location of element i [–] shape-dependent convection model is expected to be manageable by a PC. However, this model developed for the 2-D computational domain cannot be directly applied to 3-D topology optimization because the variation of the heat transfer coefficient along the direction of gravity, which is caused by the development of the thermal boundary layer, was neglected in the model for the 2-D domain. Therefore, for efficient 3-D topology optimization, it is necessary to develop a new convection model in which the heat transfer coefficient can simply be evaluated depending not only on the shape of the fins but also on the development of the thermal boundary layer. In this study, heat sinks in natural convection are thermally optimized using 3-D topology optimization. To optimize the design of the fin shape in a computationally efficient way, a simple but accurate shape-dependent convection model that is applicable to 3-D topology optimization is suggested. All other details about the optimization method, including the governing equations, objective function, constraints, and optimization procedure, are the same as our previous work . In order to validate the proposed model, it is checked whether the fin geometry of the topology-optimized design matches the multiscale structures that are known to enhance thermal performance by utilizing the unheated fluid at the entrance region. For further validation, the effective heat transfer coefficient obtained by the proposed model is compared to that obtained from numerical simulations. After the validation, the complex fin geometry of the topology-optimized heat sink is simplified to improve its manufacturability while maintaining the essential features of the fin shape favorable to the thermal performance. Finally, the thermal performance of the newly obtained heat sink is compared to that of the radial platefin heat sink optimized using the existing correlation under the constraint of the same physical volume. 2. The shape-dependent convection model for 3-D topology optimization In this study, topology optimization is performed with a new shape-dependent convection model that is applicable to a 3-D Greek symbol a THERMAL diffusivity [m2/s] b volumetric thermal expansion coefficient [1/K] C boundary of design domain [–] c relative density [–] g fin efficiency [–] l dynamic viscosity [N-s/m2] m kinematic viscosity [m2/s] q density [kg/m3] X computational domain [–] Subscripts conv convection eff effective f fluid fin fin ins insulation L heat sink length q heat flux S solid sur surface computational domain. This model predicts the local heat transfer coefficient in natural convection without solving a full conjugate heat transfer problem. In this section, the shape-dependent convection model is suggested based on the existing correlation for a pin-fin heat sink . This model accounts for the variation of the heat transfer coefficient along the direction of gravity, which is caused by developing thermal boundary layers. To perform topology optimization with low computational cost, a surrogate model predicting the local heat transfer coefficient for a 2-D computational domain was suggested in our previous study . The effective channel spacing, which offers a minimum distance to neighboring structures, was introduced to reflect the local shape in the prediction of the local heat transfer coefficient. In 3-D topology optimization, the concept of the effective channel spacing is extended in order to obtain both horizontal and vertical spacing under an arbitrary shape during optimization. A staggered pin-fin array under consideration is shown in Fig. 1. The horizontal spacing, Sh, is the horizontal distance between the vertical single arrays, while the vertical spacing, Sv, is the distance between the pin-fins in a single vertical array. For pin-fin heat sinks in natural convection with a vertically oriented base plate as depicted in Fig. 1, the correlation of the heat transfer coefficient was proposed in the previous study, Ref. : 81 8 ¼ f ðS ; Sv ; dÞ ¼ ðh1:3 þ h1:3 þ h1:3 Þ1:3 þ h8 h ; fin h fin;1 fin;2 fin;3 fin;4 ð1Þ hfin;1 ¼ Sh Sv 4Sh Sv pd qf cp gbgfin ðT b T 1 Þ ; pdL 48 mf ð2Þ hfin;2 ¼ kf Sv 0:3669 0:0494 Gr1=4 L ; L d ð3Þ hfin;3 ¼ kf ; ½2:132Sh 0:4064Ra0:188 d d ð4Þ kf ; 0:85Ra0:188 d d ð5Þ 2 and hfin;4 ¼ 34 Y. Joo et al. / International Journal of Heat and Mass Transfer 127 (2018) 32–40 Fig. 2. Schematic for the calculation of the effective fin diameter and channel spacing. expressed in Eq. (7), the value of the relative density at the k-th element is checked to determine whether the condition for reaching the neighbor structures is satisfied. As shown in Fig. 3, this channel spacing is decomposed into the effective horizontal and vertical spacing as Sh;eff;j ¼ wc;j cos h; ð10Þ and Sv;eff;j ¼ wc;j sin h: Fig. 1. Configuration and geometric parameters of a pin-fin heat sink. where qf is the density of the fluid; cp is the specific heat; g is the standard acceleration of gravity; b is the volumetric thermal expansion coefficient; gfin is the fin efficiency; Tb is the base temperature; T1 is the ambient temperature; mf is the kinematic viscosity of the fluid; kf is the thermal conductivity of the fluid; Gr is the Grashof number; af is the thermal diffusivity of the fluid; and Ra is the Rayleigh number. This correlation requires the fin diameter (d), the horizontal spacing (Sh), and the vertical spacing (Sv) to predict the average heat transfer coefficient for a fin array. In order to utilize this correlation in 3-D topology optimization, three new parameters are introduced under a local shape of a structure: the effective fin diameter (deff), the effective horizontal spacing (Sh,eff), and the effective vertical spacing (Sv,eff). The design variable for topology optimization is the relative density (c) that ranges from 0 to 1. Elements with a relative density of 1 are considered as solid region, and elements with a relative density of 0 are considered as void region. For a solid-void interfacial element i, there are 26 neighbor elements named j. The effective fin diameter is obtained by calculating the maximum thickness defined as deff ¼ maxfjxk xi jjcðxi Þ cðxj Þ < 0; cðxk Þ < 0:25; j 2 Ne g; ð6Þ xk ¼ xi þ lðxj xi Þ; ð7Þ l ¼ 1; 2; :::; 1; The heat transfer coefficient of the i-th element is determined by choosing the direction where the heat transfer coefficient predicted by Eq. (1) has the minimum value h fin;i ¼ minfhfin;j hfin;j ¼ f ðSh;eff;j ; Sv;eff;j ; deff Þ; j 2 N e g ð8Þ xk ¼ xi þ lðxj xi Þ; ð9Þ l ¼ 1; 2; :::; 1: Fig. 2 presents the schematic for the calculation of the effective fin diameter and the channel spacing. The calculation for the distance to the neighbor structures starts in the direction of the j-th element. While increasing the distance by |xj-xi| step by step as ð12Þ which is a function of local geometric parameters. In 3-D topology optimization, the variation of the heat transfer coefficient along the direction of gravity needs to be considered. When there is a vertical plate that has a higher surface temperature than the ambient air, the thermal boundary layer is formed near the wall with the natural convective flow. It is known that the local heat transfer coefficient decreases along the flow direction, and the relationship between the local heat transfer coefficient and the location along the vertical direction is given by  hz / z1=4 : ð13Þ In order for the shape-dependent convection model to reflect the variation of the heat transfer coefficient in the z-direction, the variation profile expressed in Eq. (13) was combined with Eq. (1) as hfin ¼ 3 1=4 1=4 L hfin z 4 where x is the spatial location; c is the relative density of element; and Ne is the 26 adjacent elements of the i-th element. According to Eq. (6), the distances only along the directions where cðxi Þ cðxj Þ < 0 are considered in determining the maximum distance. Before obtaining the effective horizontal and vertical spacing, the channel spacing along the (xj-xi) directions, which means the distance to neighboring structures, is obtained as wc;j ¼ fjxk xi jjcðxi Þ cðxj Þ > 0; cðxk Þ > 0:25; j 2 Ne g; ð11Þ Fig. 3. Decomposition of wc into Sh,eff and Sv,eff. ð14Þ 35 Y. Joo et al. / International Journal of Heat and Mass Transfer 127 (2018) 32–40 that satisfies 1 L Z L : hfin dz ¼ h fin ð15Þ 0 The average heat transfer coefficient appearing in Eq. (14) is obtained by Eq. (12). According to Eq. (14), the local heat transfer coefficient is the function of the effective fin diameter, horizontal spacing, vertical spacing, and location along the z-axis. 3. Results and discussion In the previous section, the shape-dependent convection model accounting for the variation of the heat transfer coefficient in natural convection was proposed. In this section, the proposed model is used to perform 3-D topology optimization under a specific computational domain. The physical validity of the newly obtained design is checked by the comparison of the topology-optimized design to the multiscale structures for convection problems described in Chapter 6 of Ref. . To estimate the thermal performance of the newly suggested designs, numerical simulations are performed. Finally, a simplified design for 3-D topology optimization is developed in order to enhance its manufacturability. 3.1. Problem setup Fig. 4. Cross-section of the computational domain (x-y plane). Topology optimization with the shape-dependent convection model was performed under a 3-D rectangular computational domain. The 2-D domain in Fig. 4 is a cross-section of the 3-D domain. To compare the results of 3-D topology optimization with radial plate-fin heat sinks, the computational domain that is the same as the radial plate-fin heat sinks was chosen. The center Table 1 Properties of solid and fluid materials. Solid Fluid (Air) Void k (W/m-K) q (kg/m3) cp (J/kg-K) l (Pa-s) b (K1) 220 0.026 0.02 – 1.16 – – 1007 – – 1.85 105 – – 3.34 103 – Fig. 5. Optimization history: (a) 25th iteration, (b) 45th iteration, (c) 150th iteration, and (d) final iteration (995th iteration). 36 Y. Joo et al. / International Journal of Heat and Mass Transfer 127 (2018) 32–40 cylinder with a radius of 10 mm was set as the fixed structure. The space beyond the outer radius (R) was set as the void region, so structures are not designed on this region. The outer radius of the domain, and the fin height were set to 20 mm and 10 mm, respectively. The length of the domain was set to 40 mm. The computational domain was discretized using 60 60 40 cubic elements, and computations were performed with Intel Core i5 4core 3.5 GHz processors. The volume fraction was given as 0.084. The filter radius was set to 2. The adiabatic boundary conditions were applied at all boundary surfaces of the computational domain. Uniform heat generation of q_ = 1.59 105 W/m3 was assumed at the elements belonging to the center cylinder. Table 1 presents the properties of solid and fluid materials used in this study. A large number of fins was given as the initial design, as shown in Fig. 4. To reduce the computational time, only one quadrant was used for the computation. 3.2. Optimization history are observed in the 3-D topology-optimized heat sink, the geometric features are summarized as follows. Fins are not connected in the direction of gravity, so they look similar to pin-fins. In addition, fins are branched at their tips, and this shape seems to increase the surface area exposed to ambient air. These pin-like structures were found to be maintained at higher volume fractions or lower element sizes. The heat transfer coefficient is large at the edge of the outer cylinder where flat structures are observed. This is because the structures near the outer cylinder are directly facing the fresh air. The shape-dependent convection model accounts for this effect by using the distance calculation algorithm that yields a large value of the channel spacing when there is no material outside the outer radius of the domain. For the analysis on the geometric characteristics, the number of local fins along the z-direction was investigated, as shown in Fig. 7. According to Da Silva and Bejan , the heat transfer performance of a vertically oriented plate-fin heat sink in natural convection can be enhanced by inserting additional fins at the entrance region. Fig. 5 shows the optimization history. In the earlier stages, most fins given at the initial design were removed, and wide channel spacing was formed, except in the entrance region where z is small. At the entrance region, more fins with irregular shapes were observed. At the 45th iteration, the initially given rectangular fins were cross-cut along the z-direction, and branched fin tips were observed near the outer boundary. Near the entrance region, a staggered fin array was formed with a larger number of local fins. At the 150th iteration, most fins were divided into a larger number of fins and evenly distributed in the computational domain. These fins were branched near the outer boundary to increase the area in contact with the ambient air. In the later stages, the number of fins did not change much, and some fins were straightened along the radial direction. 3.3. Physical validity of the topology-optimized design Fig. 6 presents the configurations of the optimized radial platefin heat sink and the 3-D topology-optimized heat sink. The geometry of the radial plate-fin heat sink was optimized by using the correlation suggested by An et al. . Although complex shapes Fig. 7. Number of local fins along the z-direction. Fig. 6. Heat sink configurations: (a) radial plate-fin heat sink and (b) 3-D topology-optimized heat sink. Y. Joo et al. / International Journal of Heat and Mass Transfer 127 (2018) 32–40 These additional fins utilize unheated air at the entrance region where the thermal boundary layers have not yet formed. The same characteristics were observed in the 3-D topology optimization results. The maximum number of local fins is observed at the entrance region, and then decreases and converges to a constant value. Therefore, the 3-D topology-optimized heat sink seems to utilize unheated air near the entrance region to enhance the thermal performance. 37 fication was performed. Fig. 8 depicts the geometry of the simplified 3-D heat sink. A T-shaped fin, which reflects the branched fin shape observed in the 3-D topology-optimized heat sink, was modeled with the minimum fin thickness of 1 mm. Except for some L-shaped fins near the entrance region, these uniformly sized T-shaped fins were arranged around the center cylinder. The staggered fin array with a larger number of local fins near the entrance region was reflected. Fins were aligned at specific positions in the z-direction while maintaining the total number of fins. 3.4. Design simplification The 3-D topology-optimized heat sink shown in Fig. 6 has very complex shapes that seem to be not manufacturable. To improve manufacturability while maintaining the geometric characteristics that have advantages in heat transfer performance, design simpli- Fig. 8. Simplified heat sink based on 3-D topology-optimized heat sink. 3.5. Evaluation of the thermal performance To evaluate the thermal performance of the 3-D topologyoptimized heat sinks numerically, ICEPAK , a commercial software provided by ANSYS, Inc., was used. The full-blown 3-D fluid models are used in this simulation tool. Thus, the thermal performance evaluated by this simulation tool reflects the actual fluid field caused by the complex geometries of the 3-D topologyoptimized heat sink. To obtain a 3-D physical model of the heat sink from a 3-D relative density field, the STL writer code written by Liu and Tovar  was used with a threshold value of the relative density 0.35. This value was chosen empirically to have the volume fraction maintained while transforming the relative density field into a 3-D physical model. This means the value of 0.35 properly represents the solid-void interface, so the threshold value for the condition for reaching the neighbor structures was chosen near 0.35. The threshold value of 0.25, which was used in Section 2, was chosen in order to lower the possibility of ignoring the existence of neighbor structures in the convection model by using a value that is lower than 0.35. However, the results of topology optimization were virtually unaffected by this threshold value within the range of 0.25–0.35. After generating a mesh geometry, the geometry was refined by using MeshLab, an open-source mesh-processing tool, in order to obtained a 3-D model which is more suitable for numerical simulations. This 3-D model was directly imported to the simulation tool. Fig. 9 shows the computational domain for numerical simulations. The size of the domain is 10R 10R 10L, which is large enough to exclude the domain size effect on the performance evaluation. Numerical simulations were performed only for the Fig. 9. Computational domain for numerical simulations. 38 Y. Joo et al. / International Journal of Heat and Mass Transfer 127 (2018) 32–40 quadrant by applying symmetric boundary conditions at min-x and min-y walls. On the remaining boundaries, the pressure was set to the ambient pressure and the velocities were not predetermined. Uniform heat generation was given at the center cylinder. Grid tests were conducted for the topology-optimized heat sink. The number of grids changed from 22,855 to 513,671 nodes. The changes in the average base temperature divided by the base-toambient temperature difference did not exceed 1% after the number of grids exceeded 267,305 nodes. Thus, the grid of 267,305 nodes was used in the numerical simulations. Other details about the numerical simulations are the same as in the method suggested by Joo and Kim . The thermal resistance given by Rth ¼ Tb T1 ; Q ð16Þ where Q is the total heat input, was used as the performance index. Table 2 presents the results of the numerical simulations. The thermal resistance of the radial plate-fin heat sink in Fig. 6(a) was also evaluated numerically under the same computational domain for the comparison. The volume fraction (f) in Table 2 reflects the mass of a heat sink, which is another performance index for the comparison. The 3-D topology-optimized heat sink has 13% lower thermal resistance and 48% less mass than the optimized radial plate-fin heat sink. After the design simplification, the thermal resistance increased by 1.2%, but it is still much lower than that of the radial plate-fin heat sink. 3.6. Accuracy and efficiency of the shape-dependent convection model In order to validate the proposed shape-dependent convection model, the effective heat transfer coefficient obtained from this model was compared to that obtained from the numerical simulations, as shown in Table 2. The effective heat transfer coefficient is defined as ¼ h eff Q AðT b T 1 Þ ð17Þ where A is the surface area of a heat sink. The difference in the effective heat transfer coefficient is about 11%. Even with the simple shape-dependent convection model, the convective heat transfer from complex 3-D shapes is fairly well predicted. The spatial variation of the heat transfer coefficient was also investigated. Fig. 10 shows how the distribution of the heat transfer coefficient used in topology optimization compares to that obtained from the numerical simulation. It can be seen that the large values of the heat transfer coefficient are observed near the outer radius of the domain in both results. The heat transfer coefficient decreases gradually along the z-direction while having a large value in the entrance region. This is why the relatively large number of fins is located near the entrance region. Fig. 11 shows the temperature distributions from topology optimization and the numerical simulation. The maximum temperature is observed at the heat source, and the decrease of the temperature along the radial direction in a fin structure confirms that the heat generated within the source flows through the fin structures. Fig. 12 presents the comparison of the interface temperature distribution obtained Fig. 10. Distributions of the heat transfer coefficient: (a) the shape-dependent convection model (b) numerical simulation by ICEPAK. from topology optimization to that obtained from the numerical simulation. In the result of the numerical simulation, the lower part of the heat sink has lower temperatures than the upper part. This means that for the same geometry in this lower part, the full Table 2 Thermal resistances and effective heat transfer coefficients of the topology-optimized heat sinks and the optimized radial plate-fin heat sink at Q = 2 W. Case f (-) Rth,numerical (K/W) 2 h eff;numerical (W/m -K) 2 h eff;topology (W/m -K) Original topology-optimized HS (Fig. 6(b)) Simplified topology-optimized HS (Fig. 8) Optimized radial plate-fin HS (Fig. 6(a)) 0.084 0.093 0.163 10.44 10.57 11.76 6.82 6.96 5.19 6.13 – – Y. Joo et al. / International Journal of Heat and Mass Transfer 127 (2018) 32–40 39 Fig. 11. Distributions of the temperature at the entrance region (min-z): (a) topology optimization (b) numerical simulation by ICEPAK. conjugate heat transfer model predicted the higher cooling performance than the shape-dependent convection model. It was confirmed that the shape-dependent convection model fairly well predicts the spatial variation of the heat transfer coefficient. However, it should be noted that this convection model is not generally applicable to all types of 3-D design problems in its present form. By extending the convection model to cover the shapes of vertical fins and horizontal plates, the presented method will be able to deal with various types of 3-D design problems. The computational time required for 500 iterations was about 12 h in this study with 4 cores of processors. The shapedependent convection model took 38% of the total computation time, and this is lower than the time taken by the finite element analysis solving only the conduction equation. It was reported that the 3-D topology optimization using the full-blown fluid model took about 10 h, with 1280 cores of processors at a mesh resolution of 160 320 160 . With the same mesh resolution and the number of cores, the computational time for using the shapedependent convection model is expected to be about 1/5 of that using the full-blown fluid model. The computation time for solving the full conjugate heat transfer problem can be reduced by implementing the local meshing. On the other hand, the shapedependent convection model presented in this study is going to Fig. 12. Distributions of the interface temperature: (a) topology optimization (b) numerical simulation by ICEPAK. be combined with a more efficient topology optimization method currently being studied, so that the computational time can be even further reduced. 4. Conclusion In this study, heat sinks in natural convection were thermally optimized using the three-dimensional topology optimization method. A new shape-dependent convection model that provides accurate predictions of the heat transfer coefficient while having a significantly less computational load was developed. This model accounts for the variation of the heat transfer coefficient depending not only on the local shape of the fins but also on the development of the thermal boundary layer. The effective heat transfer coefficient obtained from the proposed model was in close 40 Y. Joo et al. / International Journal of Heat and Mass Transfer 127 (2018) 32–40 agreement with that obtained from numerical simulations, within an error of 11%. The new design obtained from three-dimensional topology optimization had the multiscale structures that enhanced thermal performance by utilizing the unheated fluid at the entrance region. Based on the original topology-optimized heat sink that had complex shapes, a more manufacturable design was suggested via design simplification. In order to compare the thermal performance of the topology-optimized designs against that of an existing optimized design, the fin geometry of a heat sink was optimized using two methods: the method proposed in this study for the former, and an analytical method based on the existing correlation for the latter. From the comparison, it was found that the topology-optimized heat sink has 13% lower thermal resistance and 48% less mass than the optimized radial plate-fin heat sink. Therefore, three-dimensional topology optimization suggested in this study is expected to provide novel designs with improved thermal performance and reduced mass for various applications. Conflict of interest None declared. Acknowledgement This work was supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIT) (No. 2012R1A3A2026427). References  T.H. Kim, D.-K. Kim, K.H. Do, Correlation for the fin Nusselt number of natural convective heat sinks with vertically oriented plate-fins, Heat Mass Transf. 49 (2013) 413–425.  B. Li, Y.-J. Baik, C. Byon, Enhanced natural convection heat transfer of a chimney-based radial heat sink, Energy Convers. Manage. 108 (2016) 422–428.  Y. Joo, S.J. Kim, Comparison of thermal performance between plate-fin and pinfin heat sinks in natural convection, Int. J. Heat Mass Transf. 83 (2015) 345– 356.  A. Bejan, A.J. Fowler, G. Stanescu, The optimal spacing between horizontal cylinders in a fixed volume cooled by natural convection, Int. J. Heat Mass Transf. 38 (1995) (1995) 2047–2055.  A. Bejan, Convection Heat Transfer, fourth ed., Wiley, 2004, pp. 180–215.  A. Bar-Cohen, W.M. Rohsenow, Thermally optimum spacing of vertical, natural convection cooled, parallel plates, J. Heat Transfer 106 (1984) 116–123.  E.M. Sparrow, S.B. Vemuri, Orientation effects on natural convection/radiation heat transfer from pin-fin arrays, Int. J. Heat Mass Transf. 29 (1986) 359–368.  M.P. Bendsøe, O. Sigmund, Topology Optimization: Theory, Methods, and Applications, Springer, 2003.  E.M. Dede, S.N. Joshi, F. Zhou, Topology optimization, additive layer manufacturing, and experimental testing of an air-cooled heat sink, J. Mech. Des. 137 (2015) 111403.  A. Iga, S. Nishiwaki, K. Izui, M. Yoshimura, Topology optimization for thermal conductors considering design-dependent effects, including heat conduction and convection, Int. J. Heat Mass Transf. 52 (2009) 2721–2732.  M. Zhou, J. Alexandersen, O. Sigmund, C.B.W. Pedersen, Industrial application of topology optimization for combined conductive and convective heat transfer problems, Struct. Multidiscip. Optim. 54 (2016) 1045–1060.  P.T. Lin, M.C.E. Manuel, J. Zhang, Y. Jaluria, H.C. Gea, Multi-objective design optimization of multiple microchannel heat transfer systems based on multiple prioritized preferences, J. Therm. Sci. Eng. Appl. 9 (2) (2017) 021011.  Y. Joo, I. Lee, S.J. Kim, Topology optimization of heat sinks in natural convection considering the effect of shape-dependent heat transfer coefficient, Int. J. Heat Mass Transf. 109 (2017) 123–133.  J. Alexandersen, N. Aage, C.S. Andreasen, O. Sigmund, Topology optimisation for natural convection problems, Int. J. Numer. Meth. Fluids 76 (2014) 699– 721.  J. Alexandersen, O. Sigmund, N. Aage, Large scale three-dimensional topology optimisation of heat sinks cooled by natural convection, Int. J. Heat Mass Transf. 100 (2016) 876–891.  J. Alexandersen, O. Sigmund, K.E. Meyer, B.S. Lazarov, Design of passive coolers for light-emitting diode lamps using topology optimization, Int. J. Heat Mass Transf. 122 (2018) 138–149.  A. Bejan, S. Lorente, Design with Constructal Theory, John Wiley and Sons, 2008, pp. 215–248.  B.H. An, H.J. Kim, D.-K. Kim, Nusselt number correlation for natural convection from vertical cylinders with vertically oriented plate fins, Exp. Therm Fluid Sci. 41 (2012) 59–66.  A.K. Da Silva, A. Bejan, Constructal multi-scale structure for maximal heat transfer density in natural convection, Int. J. Heat Fluid Flow 26 (1) (2005) 34– 44.  ANSYSÒ, ICEPAK, Release 16.1, User’s Guide, ANSYS, Inc.  K. Liu, A. Tovar, An efficient 3D topology optimization code written in Matlab, Struct. Multidiscip. Optim. 50 (6) (2014) 1175–1196.  Y. Joo, S.J. Kim, Thermal optimization of vertically oriented, internally finned tubes in natural convection, Int. J. Heat Mass Transf. 93 (2016) 991–999.