Topology optimization of liquid-cooled microchannel heat sinks: An experimental and numerical study

International Journal of Heat and Mass Transfer 142 (2019) 118401

Contents lists available at ScienceDirect

International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt

Topology optimization of liquid-cooled microchannel heat sinks: An experimental and numerical study Shi Zeng, Poh Seng Lee ⇑ Department of Mechanical Engineering, National University of Singapore, 9 Engineering Drive 1, 117576 Singapore, Singapore

a r t i c l e

i n f o

Article history: Received 7 March 2019 Received in revised form 9 July 2019 Accepted 9 July 2019

a b s t r a c t Electronics cooling is always a critical challenge for its small profile and high heat flux. Various microchannel heat sinks have been proposed in literature to achieve reliable and energy efficient cooling effect. They are mainly different patterned fin shapes originated from researchers’ experience. In this paper, we present the process of designing liquid-cooled microchannel heat sink fin geometries based on a numerical method, topology optimization. The optimization model is developed to find designs that satisfy heat transfer requirement with low pressure drop penalty. It is implemented based on a derived accurate 2D model with minimizing pressure drop as objective and average junction temperature as a constraint. Parameter study for different velocities and junction temperature constraints is performed. The optimized structures are then post processed under geometry constraints and analyzed by 3D CFD simulation. Size optimization is implemented on conventional straight channel heat sinks, which serve as the benchmark cases. Moreover, an experimental test loop is built to validate the 3D simulation of topology optimized heat sinks and straight channel heat sinks. The performance comparison shows that the topology optimized heat sinks could save up to 50.9% pumping power under the same thermal performance requirement. Detailed CFD analysis is then engaged to examine their thermohydraulic characteristics and reveal the reasons for their superior performance. Ó 2019 Elsevier Ltd. All rights reserved.

1. Introduction Thermal management of electronics is becoming more and more challenging with the advancement of chips in miniaturization and performance. The high power density is about to make air cooling infeasible. While two-phase cooling can reach excessive heat flux, its application is limited by the complexity of fabrication, assembling and operation [1]. It is expected that single phase liquid cooling will be the dominant technology in the near future [2]. In addition, cooling strategy has to be energy efficient, especially with the popularity of data centers. Energy consumed by cooling infrastructures typically counts for 20–40% of total energy consumption of data centers [3]. It is reported that water-cooled systems could reduce energy consumption compared to aircooled systems [4]. In 1980s, Tuckerman and Pease [5] showed that microchannel heat sink was capable of dissipating a heat flux of 790 W=cm2 with a pressure drop as high as 31 psi. Since then, designing and optimizing microchannel or minichannel cooling systems has attracted a lot of interest. Enhancing heat transfer performance and decreas-

⇑ Corresponding author. E-mail addresses: [email protected] (S. Zeng), [email protected] (P.S. Lee). https://doi.org/10.1016/j.ijheatmasstransfer.2019.07.051 0017-9310/Ó 2019 Elsevier Ltd. All rights reserved.

ing pressure drop are two key factors for the success of such a system. Researchers have proposed various fin configurations in literature, such as the diamond-shaped interrupted fin [6], pin fins with various shapes [7] and the oblique fin [8]. On the other hand, based on the view of channel modification, Sui et al. [9] examined the flow phenomenon in wavy channel with rectangular cross section. Fluid mixing, Dean vortices and chaotic advection were identified as the main factors for enhancing heat transfer. Mills et al. [10] reported that sinusoidal wall wavy channels were better than straight channels when Reynolds numbers were higher than 280. Al-Neama et al. [11] found the serpentine configuration could reduce 19% of total thermal resistance while its pressure drop was up to ten times of straight rectangular channels. Besides, many experimental and numerical studies have investigated effects of introducing various kinds of interruptions into channels, like dimple [12], vortex generators [13], herringbone microstructures [2] and inverted fish scale profile [14]. More similar structures are summarized in [15]. Moreover, various novel structures have been proposed at the cooling system level. Split-flow is an effective concept to reduce pressure drop because both the flow length and fluid flow rate in a single channel are reduced compared to the single-pass arrangement [16]. The manifold microchannel heat sink could be seen as a further development of the split-flow idea. Escher et al. [17,18]

2

S. Zeng, P.S. Lee / International Journal of Heat and Mass Transfer 142 (2019) 118401

Nomenclature A c h H I k _ m P q_ € q qf =qk =qh Q R SC TO T u V v V_ L

area (m2 ) specific heat (J=kg
K) convection coefficient (W=m2
K) height (m) current (A) thermal conductivity (W=m
K) mass flow rate (kg=s) pressure (Pa) heat transfer rate (W) heat flux (W=m2 ) convexity parameter power (W) thermal resistance (°C/W) straight channel topology optimized temperature (°C) velocity vector voltage (V) velocity (m/s) volume flow rate (m3 =s) length (m)

Greek symbols D gradient (difference) a friction coefficient

designed an ultrathin heat sink based on the manifold microchannel concept and proved that it could achieve similar thermal resistance with an almost 20 times smaller pumping power compared to conventional single pass channels. Sharma et al. [19] demonstrated a hotspot-targeted microchannel design based on manifold microchannel architecture and proved it could reduce temperature variance by 57% while pumping power was less than 0.3% of the total chip power. The designs mentioned above are all originated from experience or intuition. Recently, a new numerical method called topology optimization is gaining popularity for its pure numerical based design capability. This method was first invented by Bendsøe and Kikuchi [20] in the area of solid mechanics, then extended by Borrvall and Petersson [21] in 2003 to cover fluid mechanics. It is a material distribution method, in which every discretized element of the design domain is decided to be either fluid or solid. In density based topology optimization, the discrete solid or fluid distribution is relaxed through introducing a material index, which takes values from 0 to 1. The two extremities (0 and 1) represent solid and fluid material, respectively. For values between 0 and 1, they represent intermediate material, which is allowed to exist in the process of optimization but expected to converge to either 0 or 1 in the end. The material index is updated to approach optimum based on gradient information of objective functions and constraints. The primary advantage of topology optimization over traditional numerical optimization methods, e.g., size and shape optimization, is its capability of generating new boundaries and no need for predefined geometries. Many works have been done to apply this method in coupled thermal-fluid problems. Alexandersen et al. [22] presented a work on heat sink design by natural convection. The full 3D conjugate heat transfer optimization problem with up to 330 million degrees of freedom was solved in a parallel computation framework of PETSc. It was found that branches in optimized structures evolved from simple to complex as the Grashof numbers increased. The same method was applied to design passive coolers for light-emitting diode lamps [23]. Comparing

l q c

dynamic viscosity (Pa
s) density (kg=m3 ) local design variable

Subscripts ave average b base bt base bottom surface fl fluid HS heat sink in inlet min minimum max maximum out outlet s solid ff fluid & fin layer th thermal wet wet surface fb fin base layer unfin fin base unfinned area j junction en environment sp shallow plenum dp deep plenum

with lattice-fin design, the optimized designs reduced the package temperature by 26% and material usage by 12%. In terms of forced convection problems, Yoon [24] presented a simple 2D case of cooling a hot plat. Inverse permeability, heat conductivity, density and specific heat capacity were all interpolated between fluid and solid. He highlighted the local optimum issue that different initial guess would lead to different optimized structures. Similarly, Koga et al. [25] demonstrated a complete process of developing heat sinks by topology optimization in stokes flow region, including optimization in a 2D domain, 3D computational simulation validation and experimental verification. Pressure drop minimization and heat dissipation maximization were combined into one objective function. Final channel layouts depended on weighting factors of the two competing objectives and penalty factors in property interpolation functions. Dede [26] studied a case with center top inlet and multipass branching channel network. Differently, only inverse permeability and heat conductivity were interpolated. Instead of manually setting weighting factors of two objectives, Sato et al. [27] proposed an adaptive weighting scheme in combination with topology optimization to explore different designs. The new scheme was proved to be capable of effectively escaping local Pareto fronts. Dilgen et al. [28] did 2D and 3D topology optimization in the turbulence region. Turbulence modelling led to more complicated solvers but allowed more complex and higher velocity flow, which could result in better heat sink performance comparing to optimization results confined in laminar region. In order to reduce computational cost, Zhao et al. [29] presented a ‘‘poor man’s approach” to cope with turbulence modeling in topology optimization. Though this approach could not yield a perfectly accurate velocity field, cooling channel designs from this method performed well. The same idea has been extended to solve natural convection conjugate problems recently [30]. The number of degree of freedom could be reduced by 50% using the reduced model and optimized results were found to be qualitatively similar as those from full Navier Stokes-based model.

S. Zeng, P.S. Lee / International Journal of Heat and Mass Transfer 142 (2019) 118401

Except for study of heat sinks, heat exchanger designs by topology optimization were also reported in literature. Qian [31] defined a 2D simplification of a standard shell and tube heat exchanger. Water path surrounding circular tubes was optimized with a surface tangential thermal gradient. Haertel and Nellis [32] targeted on air-side surface optimization of condensers. A simplified 2D model was first developed with fully developed internal flow assumption. Maximization of heat exchanger conductance was set as the objective function while pressure drop and air temperature change were prescribed. The topology optimized designs could outperform size optimized slot geometry by up to 71%. Recently, Haertel et al. [33] presented a work on heat sink design with a pseudo 3D model. Pressure drop and heat production rate were fixed. Minimizing heat transfer resistance was set as the objective. Higher pressure drops were found to generate more complex channel layouts. A 13.6% of thermal resistance reduction was achieved when comparing with a size optimized parallel fin heat sink. Our last work [34] was also focused on planar air heat sink design problem. A similar 2D two-layer model was developed and verified as an accurate simplification of the full 3D conjugate heat transfer model. Differently, longer design domain allowed more sectional fins to be generated in the stream wise direction through topology optimization and effective flow modulation behavior for enhancing heat transfer was observed. Experimental tests proved that the optimized structure was able to save 54.9% pumping power comparing to a traditional straight fin. The work presented in this paper is a direct extension from our previous work by changing cooling agent from air to water. The 2D model for predicting heat sink thermohydraulic characteristics is further corrected by a post processing step of assuming 1D conduction in the heat sink base. Topology optimization is set up based on the 2D model to seek for microchannel designs with high heat transfer performance and low pressure drop penalty. Then the topology optimized (TO) microchannel heat sinks are benchmarked with size optimized straight channel (SC) heat sinks by 3D numerical analysis and experimental study. Together with our previous work [34], we intend to provide solid proof of the effectiveness of topology optimization in designing micro/mini scale gas/liquid heat sinks with experimental validation. In addition, more insights about the formation process of fin structures during topology optimization and the simplified 2D model are presented. This paper is organized into the following sections: Section 2 describes the topology optimization model setup and implementation, including the development of the 2D simplified model for heat sink simulation. Topology optimized results and post processing are presented in Section 2.4. Section 3 describes the 3D numerical analysis method used for verifying performance of TO heat sinks and doing size optimization for SC heat sinks. Section 4 describes the experimental setup and experiment procedures. Performance comparison and discussion are presented in Section 5, followed by conclusions in Section 6.

sink is 45 mm long and 32 mm wide; the base is 4 mm high; the fin is 1.2 mm high; water inlet temperature is 25 °C; heat source heat flux is 25 W=cm2 . Different fin and channel structures will result in different thermal and hydraulic performance. Hence, the finned region is set as the design domain, where topology optimization is implemented to generate new designs. In order to reduce the computational expense, only a small spanwise strip of heat sink, 2 mm wide, is numerically studied and optimized as a representation of the full heat sink by applying symmetry boundary conditions on both sides. The reduced design domain and boundary conditions are shown in Fig. 2. Mirroring the generated structure 16 times then could get the full domain design without difference in performance. The heat sink is subjected to a constant heat flux from the bottom surface and a constant water flow rate. In terms of performance evaluation criteria, lower bottom surface average temperature and lower pressure drop to drive the water are considered favorable. 2.2. 2D thermofluid model of water-cooled heat sinks As an iterative optimization method, topology optimization gradually evolves structures in the design domain to the optimum based on the flow and thermal fields of last iteration. CFD analysis is performed at each iteration to provide the information. Even though we only study a 2 mm wide stripe of the heat sink, a full 3D CFD analysis still requires a lot of time considering that the optimization may take hundreds of iterations. Hence, a simplified 2D model is first developed to perform CFD analysis of watercooled heat sinks for the use of topology optimization. This 2D model is similar to that developed in our previous work [34] except that one more layer is added due to the large temperature gradient from the base bottom surface to the fin base surface. It is also similar to the model used in [33] except that the model accuracy is improved by utilizing a physically derived heat transfer coefficient. 2.2.1. Governing equations As shown in Fig. 3, a heat sink is divided into three layers in the height direction. The first layer is set at the middle of the fin height, representing the interaction between water and fins. The governing equations of fluid domains in this layer are written as following:

r
u¼0

ð1Þ

qfl ðu
rÞu ¼ rP þ lr2 u

ð2Þ




2.1. Description of the design problem This study is aimed to generate non-conventional water-cooled microchannel heat sinks through topology optimization and verify their performance against conventional SC heat sinks. A typical SC heat sink is shown in Fig. 1. The bottom surface is attached to a heat source, which is a high heat flux chip considered in this study. Heat is conducted upwards through the base and taken away by water pumped through the channels. Some critical dimensions and operation conditions considered in this study are: the full heat



qfl cfl u
rT ff ¼ r
kfl rT ff þ


 hfl T fb  T ff Hfin

ð3Þ

The governing equation of solid domains in this layer is only energy equation




2. Topology optimization

3



r
ks rT ff þ


 hs T fb  T ff ¼0 Hfin

ð4Þ

where T ff is the first layer (fluid & fin layer) temperature field, T fb is the second layer (fin base layer) temperature field, Hfin is the fin height, and qfl ; cfl ; kfl are the density, specific heat capacity and conductivity of liquid, respectively. hfl is the convection heat transfer coefficient between the fluid and fin base and hs is the nominal convection heat transfer coefficient between fins and the fin base. The derivation of these two coefficients will be discussed in Section 2.2.2. The second layer represents the fin base surface, where heat partition is determined, e.g., heat is either conducted into fins or

4

S. Zeng, P.S. Lee / International Journal of Heat and Mass Transfer 142 (2019) 118401

water outlet

4 mm

water inlet

= 25℃ ̈

= 25 W/cm2

Fig. 1. The full domain view of a SC heat sink with boundary conditions.

Fig. 2. 2 mm wide optimization domain and symmetry boundary conditions on both sides.

Layer 1 fluid & fin Layer 2 fin base

Layer 3 bottom

Fig. 3. Illustration of the 2D thermofluid model.

convectively transferred to water. The governing equations are only energy equations: In domains beneath fluid domains of the first layer:




r
ks rT fb




 €bt hfl T fb  T ff q  ¼0 þ Hb Hb






r
ks rT fb þ

€bt
T bt ¼ T fb þ q

ð6Þ

The third layer was not involved in the 2D model for air-cooled €bt in that case is small and the heat sink base heat sinks because q temperature was assumed uniform. It has to be mentioned that only the first and second layer need to be defined and solved in Comsol Multiphysics for numerical simulation. The temperature



€bt hs T fb  T ff q  ¼0 Hb Hb

Hb ks

ð5Þ

In domains beneath solid domains of the first layer:




The third layer is the bottom surface of the heat sink. Assuming uniform 1D vertical heat conduction from the bottom to the fin base, the bottom surface temperature could be calculated as:

€bt is the heat flux, and Hb is the heat sink base thickness. where q

ð7Þ

5

S. Zeng, P.S. Lee / International Journal of Heat and Mass Transfer 142 (2019) 118401

of the third layer is calculated manually from the simulated second layer data by Eq. (7). It is a post processing step after the numerical modeling in Comsol Multiphysics. 2.2.2. Determination of hf l and hs The temperature prediction accuracy of this 2D model is mainly dependent on the two heat transfer coefficients hfl and hs . They determined the thermal resistance from the second layer to the first layer and the heat partition from fin base to fins or fluid. hfl reflects the convection heat transfer resistance between the unfinned base surface and fluid. 3D simulation of a SC heat sink with fin width of 0.3 mm and channel width of 0.5 mm is used to determine this value. The 3D simulation method will be discussed in Section 3 in detail. Then hfl is calculated as

hfl ¼

q_ unfin Aunfin ðT unfin  T fl Þ

ð8Þ

where q_ unfin is the heat transfer rate from the fin base unfinned area to fluid, Aunfin is the fin base unfinned area, T unfin is the average wall temperature of the fin base unfinned surface, and T fl is the volumeweighted average temperature of fluid in the channel domain. hfl increases with the increase of the inlet velocity, as listed in Table 1. hs reflects the conduction heat transfer resistance between the finned base surface and the fins, so it is called nominal convection heat transfer coefficient. We use the heat conduction resistance from the fin base to the middle height of fins to approximate its value. It is calculated as

hs ¼

ks ¼ 646000 W=ðm2
KÞ 0:5Hfin

ð9Þ

2.3. Implementation of topology optimization In this section, we will introduce how to apply density based topology optimization in the design problem described in Section 2.1 with the developed 2D model in Section 2.2. The principle is to unify fluid and solid by including a design variable c to represent the phase of every discretized element in the design domain. Elements with c ¼ 1 represent fluid and elements with c ¼ 0 represent solid. The best design is then found through optimizing values of the design variables. More basic introduction about topology optimization of fluid mechanics could be seen in [35–37].



cðxÞ ¼

1; if x 2 Xf 0; if x 2 Xs

where x is the special coordinate, Xf is the fluid domain, and Xs is the solid domain. 2.3.1. Modified governing equations In topology optimization, the governing equations introduced in Section 2.2 are modified to unify fluid and solid and incorporate influence of design variables: Layer 1:

qfl ðu
ruÞ ¼ rP þ lr2 u  aðcÞu

ð10Þ






 q €bt hðcÞ T fb  T ff r
ks rT fb þ  ¼0 Hb Hb

hfl (W=ðm
KÞ)

0.05 0.1 0.15 0.2 0.25

3544.5 3984.6 4357.8 4668.4 4940.2

2

ð11Þ

ð12Þ

where au is an extra friction force term [21,38]. a is the friction force coefficient or inverse permeability. Theoretically, a should satisfy that:



aðcÞ ¼

0;

if x 2 Xf

1;

if x 2 Xs

in order to ensure zero flow velocity in solid domains while eliminating the extra fiction force in fluid domains. The thermal conductivity and nominal heat transfer coefficient should satisfy that:

 kðcÞ ¼

ks ;



if x 2 Xs

hðcÞ ¼

kfl ; if x 2 Xf

hs ;

if x 2 Xs

hfl ; if x 2 Xf

In topology optimization, c is relaxed to vary continuously from 0 to 1 to avoid a discrete optimization problem. Its value between 0 and 1 represents intermediate material, which is not physically available. Ideally, all elements will be pushed to either 0 or 1 at the end of optimization and clear boundaries between fluid and solid domains should be formed. 2.3.2. Interpolation functions The inverse permeability a, conductivity and nominal convection heat transfer coefficient should take different values for fluid and solid. They are relaxed to change continuously between these two extremities as a function of c. Rational Approximation of Material Properties (RAMP)-style interpolation functions [39] are utilized in this study:

a ¼ amin þ ðamax  amin Þ 

1c 1 þ qf c



hðcÞ ¼

ð13Þ



c kfl ð1 þ qk Þ  ks þ ks kðcÞ ¼ 1 þ qk c

ð14Þ



c hfl ð1 þ qh Þ  hs þ hs 1 þ qh c

ð15Þ

where amin is the minimum inverse permeability, which is set as 0 throughout this work. amax is the maximum inverse permeability, and qf , qk and qh are parameters that control the convexity of the interpolation functions, termed as convexity parameters. These four parameters play an important role in driving intermediate material phase to 0 and 1 during optimization and have great influence on the final results. 2.3.3. Density filter To avoid checker board problem and govern the minimum structure size, a Helmholtz-type PDE density filter is applied [31]: 

v in (m/s)

hðcÞ T fb  T ff Hfin



Layer 2:



r 2 r2 c þ c ¼ c Table 1 hfl values under different inlet velocities.



qfl cfl u
rT ff ¼ r
kðcÞrT ff þ




ð16Þ

where r controls the size of filtering, c is the original design vari

able, and c is the filtered design variable, which will replace c in all the functions of topology optimization model. The problem of density filter is that it will create a band of intermediate material phase in boundaries among solid and fluid domains. To avoid this issue, density filter is turned off in the last several iterations in order to get clear topologies of channel and fin designs.

6

S. Zeng, P.S. Lee / International Journal of Heat and Mass Transfer 142 (2019) 118401

2.3.4. Optimization scheme As stated in Section 2.1, the liquid-cooled heat sink design in this study is a dual-objective optimization problem. High heat transfer performance and low pressure drop driving the fluid are the pursued objectives. In this work, they are handled as a combined minimization and constraint problem: Minimize: DP Subject to:0  cðxÞ  1 Maximum bottom surface average temperature: T bt;av e  T con Continuity equation: (1) Modified N-S equation: equation (10) Energy equation: equation (11) & (12) where DP is the pressure drop through the heat sink and T con is the temperature constraint value. Optimization is conducted in several different inlet velocities and bottom surface average temperature constraints, which results in different designs for each case. In addition, different from our previous work, it is found that the volume constraint is not necessary for converging to a discrete fluid/solid design in this study. It is probably due to the different convex parameter and amax values used in the interpolation functions. 2.3.5. Implementation procedure In this section, we will describe the process of implementing topology optimization in the commercial software Comsol Multiphysics 5.2a. The physical model described in Section 2.2 with the modified governing equation in Section 2.3.1 is setup in Laminar Flow and Heat transfer module. The density filter is defined in the Coefficient Form PDE module and the optimization scheme is defined in the Optimization module. The design domain is set as c ¼ 0:5 uniformly as the start. Flow and temperature fields are calculated by solving the governing equations. The objective function and constraints are then calculated and their sensitivities with respect to the design variable are also calculated [39]. The design variable field is then updated by the globally convergent version of Method of Moving Asymptotes (GCMMA) based on the sensitivity information [40]. The adjusted design variable field is then taken into the next iteration of optimization. Constraints will be satisfied and the objective will be minimized gradually. In order to obtain a clearly separated 0/1 design variable field (only solid and fluid domain, no intermediate material phase), parameters, amax , qf , qk , qh and r, are carefully selected through numerical experiments and gradually changed during the optimization process. Their values used in the continuation strategy are shown in Table 2. It has to be mentioned that nouter;iter is the maximum number set in the optimization algorism, and some continuation steps did not run 30 iterations if the algorism already reached residual criteria of 0.001 in advance. At the last two continuation step, the density filTable 2 Values of optimization parameters and maximum number of outer iteration in GCMMA. Continuation step

nouter;iter

1

30

2

30

3

amax

qf

qk

qh

r (mm)

5  10

1000

10

10

0.075

5  106

500

50

50

0.075

30

107

250

100

100

0.05

4

30

107

125

125

125

0.05

5

30

107

100

150

150

0.025

6

30

107

50

200

200

0.025

7

30

107

25

200

200

\

8

30

107

10

200

200

\

6

ter is turned off to avoid intermediate material phase band in final optimized designs. 2.4. Optimization results and post processing Topology optimization of microchannel water heat sinks is performed under different inlet velocities v in and bottom surface average temperature constraints T con . The final design variable fields are shown in Fig. 4. It could be seen that the final design variable fields are clearly separated as solid and fluid regions. Intermediate material phase almost does not exist in the final results. For the same inlet velocity, the generated fin areas become smaller and channel structures become simpler as the constraint temperature increases. These designs are then extracted by a level curve of c ¼ 0:5 and analysed through the 2D model described in Section 2.2 with explicit boundaries between fluid and solid domains. Each design is simulated under various inlet velocities to generate the curve of average bottom surface temperature T bt;av e versus pumping power. The pumping power is calculated as

Q pump ¼ v in WHcha  DP

ð17Þ

where v in is the inlet velocity, W is the width of simulation domain, Hcha is the channel height and DP is the simulated pressure drop. Fig. 5 shows the curves of T bt;av e versus pumping powers for different topology optimization results. Designs that could achieve the minimum T bt;av e under the same pumping power or consume the smallest pumping power while maintaining the same T bt;av e are preferred. The curves of such designs would be those close to the origin in the graph. It could be seen that 0.1 m/s–47.43 °C is the best or quite close to the best throughout the considered flow range. Thus, it is selected for further 3D numerical analysis and experimental study. Different from manual designs in which geometry parameters could be explicitly defined, the channel width in topology optimization is not fully controlled in our algorithm and appears to be even smaller than 0.1 mm in some parts. These small features make manufacturing difficult or even impossible, thus requiring post processing. In this case, small channels are widened by cutting off parts of fins and slightly adjusting positions of some pieces of fins to meet the minimal channel width requirement. The modified structures of 0.1 m/s-47.43 °C design with minimal channel width of 0.3 mm and 0.5 mm are shown in Fig. 6. For simplicity, the two designs are called TO 0.3 and TO 0.5 in the following discussions. In addition to manual enlarging the gaps among fins, there is also a topological change in TO 0.5 as highlighted with an arrow in Fig. 6. The original large piece of fin is cut into two pieces and one part of it is moved to the symmetry line. The reason is that the channel between the fin and the symmetry line there is small (less than 0.1 mm) in the original topology optimization results. The flow resistance through this fine channel is high. Moreover, the nature of density based topology optimization is that fluid could penetrate through the ‘‘solid” areas even with large inverse permeability a, which leads to the local flow field in topology optimization as Fig. 7. It could be seen that there are streamlines going through the ‘‘solid” fin. If the original fine channel is enlarged to be 0.5 mm wide, the large flow resistance there would not exist, instead, it becomes a bypass channel. Therefore, the fin is cut into to two pieces and one part is moved to touch the symmetry line, closing the fine channel and creating a secondary channel. In such a way, it could mimic the large flow resistance and small penetration at the local spot. The treatment is imposed for TO 0.5 but not TO 0.3 because 0.3 mm channel could still enforce reasonable flow resistance. It has to be noted that length scale control methods are

S. Zeng, P.S. Lee / International Journal of Heat and Mass Transfer 142 (2019) 118401

7

(m/s) 63.43 ℃

65.43 ℃ 0.05 67.43 ℃

45.43 ℃

47.43 ℃ 0.1 49.43 ℃

39.43 ℃

41.43 ℃ 0.2 43.43 ℃

Fig. 4. Final design variable (c) field grayscale picture under different optimization conditions (black color represents solid regions and white color represents fluid regions. Water flows from left to right and arrows show the flow patterns).

available in literature [41–44]. With such extra regulators in topology optimization, it is possible to precisely control the minimum channel width and avoid the abovementioned troubles.

After post processing, the 0.3 mm and 0.5 mm designs are both applied in 3D numerical analysis, as described in Section 3, and the 0.5 mm one is manufactured and experimentally tested. 3. 3D numerical analysis Topology optimization is implemented based the developed 2D model, which is a balanced choice between computational cost and numerical modeling accuracy. To get the most accurate CFD analysis results for the TO heat sinks, 3D simulation is implemented in the commercial software ANSYS Fluent 18.2. Besides, it is also used for size optimization of SC heat sinks. Since CFD analysis of SC heat sinks is well established in literature, we will mainly focus on the simulation of topology optimized structures. 3.1. Numerical model setup

Fig. 5. Performance comparison of topology optimization results.

The computational domain and boundary conditions are shown in Fig. 8. The water inlet boundary is set as constant velocity v in , and water outlet is set as pressure outlet (P ¼ 0). The heat source €bt . Surfaces at both sides of all is set as a uniform heat flux q domains are symmetry boundaries, and the other surfaces are all adiabatic walls. Due to the complex flow channel network, a turbulence model, RNG  ke, with enhanced wall treatment for both pressure gradi-

8

S. Zeng, P.S. Lee / International Journal of Heat and Mass Transfer 142 (2019) 118401

Fig. 6. Modified structures of 0.1 m/s–47.43 °C design with minimal channel width of (a) 0.3 mm (b) 0.5 mm.

ent and thermal effects is used in all the flow conditions studied. More details about using this turbulence model in heat sink simulation could be found in [45]. For material properties, the heat sink base is copper, having a thermal conductivity of 387.6 W=ðm
KÞ. Water properties are temperature dependent, and are first taken from [46] then interpolated as 3rd degree polynomial functions. The coefficients are listed in Table 3. SIMPLEC is selected as the pressure-velocity coupling scheme with second order upwind as the discretization scheme. 106 is set as the residual convergence criteria for continuity, x-y-z veloc-

Fig. 7. The local velocity field in topology optimization around the manually cut fin, showing streamlines penetrating through the ‘‘solid” fin.

ities, k, e, x and 1010 for energy. ANSYS Mesh 18.2 is used to generate the computational grids. 15 gradually growing boundary layer mesh with first layer thickness of 0.004 mm is applied on the fin surface to ensureyþ < 0.5 such that the mesh requirement of enhanced wall treatment is sat-

Fig. 8. The computational domain and boundary conditions of 3D simulation of a TO heat sink.

Table 3 Coefficients of the polynomial functions for calculating water properties.

u

A1

A2

A3

B

qðkg=m3 Þ

6.070780382 1.38889E-05 2.77778E-09 1.86667E-09

0.016825298 0.016825298 1.29786E-05 1.95665E-06

1.38889E-05 0.008708786 0.00068891 26.71531664

314.5864138 0.909246839 0.081831866 7323.291323

kðW=ðm
KÞÞ cðJ=ðkg
KÞÞ lðkg=ðm
sÞÞ

uðT Þ ¼ A1 T þ A2 T 2 þ A3 T 3 þ B. where T is in the unit of Kelvin. (Valid rang 293.15 K-353.15 K).

9

S. Zeng, P.S. Lee / International Journal of Heat and Mass Transfer 142 (2019) 118401

isfied. Mesh is first generated on the top surfaces of the fluid domains, then swept down from the top to the bottom, also with fine boundary layer mesh around top and bottom surfaces of fluid domains to meet yþ criteria. As an example, the final generated mesh for TO 0.5 is shown in Fig. 9 and it has 8.5 million elements. Mesh independence study is performed by changing size parameters, including the first layer thickness of the boundary layer mesh, the face size of the top surface and the number of divisions in the sweep method, resulting in a coarser mesh of 6.4 million and a denser mesh of 10 million. The simulation results of key parameters, pressure drop and heat sink base bottom surface temperature at the highest velocity considered, are shown in Table 4. The deviation of results is very small for all three meshes, which indicates good mesh independence. For SC heat sinks, only half width of the fin and half width of the channel are simulated. The same length of inlet and outlet extension as in the TO heat sink simulation are also included. Symmetry boundary conditions are applied on all the side surfaces. Laminar model is used and the other settings are similar to the TO heat sink simulation.

3.2. Numerical data reduction

inlet pressure Pwater;in and output pressure P water;out . The equations could be expressed as

/areaweighted ¼

ð18Þ



Pn

i¼1 /i

/massweighted ¼ P n

i¼1

!

 qi ! u i A i 

!

ð19Þ

 qi ! u i A i

where / is the target parameters and n the number of cells at the respective boundaries. The volume flow rate is calculated at the water inlet boundary:

V_ water;single

strip

¼

n X !! u i A i ðsingle strip in simulationÞ

ð20Þ

i¼1

V_ water ¼ V_ water;single

strip

 # of stripes ðfull domain & experimentÞ ð21Þ

The mass flow rate is also calculated at the water inlet boundary:

_ water;single m

Area-weighted average in ANSYS Fluent is used to calculate average junction temperature T j;av e (average bottom surface temperature T bt;av e ). Mass-weighted average is used to calculate water inlet temperature T water;in , water outlet temperature T water;out , water

n 1X / jAi j A i¼1 i

strip

¼

n X

!

qi ! u i A i ðsingle strip in simulationÞ

ð22Þ

i¼1

_ water ¼ m _ water;single m

strip

 # of stripes ðfull domain & experimentÞ ð23Þ

(a)

Fins

water inlet extension

(b)

water flow channel

heat sink base

Fig. 9. (a) Top view and (b) side view of the mesh for the TO heat sink.

Table 4 Mesh independence study results.

v in ðm=sÞ

€bt (W=cm2 ) q

DP (Pa)

T av e (°C)

T min (°C)

T max (°C)

Mesh elements

0.3 0.3 0.3

25 25 25

1017.99 1018.84 1030.53

41.64 41.65 41.57

39.39 39.40 39.38

44.41 44.43 44.34

6,392,635 8,529,375 10,046,455

10

S. Zeng, P.S. Lee / International Journal of Heat and Mass Transfer 142 (2019) 118401

The pressure drop through the whole domain (including inlet and outlet extension) is calculated as

DP ¼ Pwater;in  Pwater;out

ð24Þ

The pumping power is then calculated as

_ water DP=qav e Q pump ¼ m

ð25Þ

where qav e is the mass-averaged density of water in all fluid domains of simulation. Total thermal resistance of the heat sink is calculated as

Rth ¼

T j;av e  T water;in q_ bt

ð26Þ

where q_ bt is the total heat rate supplied to the heat sink bottom surface. Finally, a geometry factor, blockage ration (volume fraction of the fin material), is defined as

BR ¼

Afin;cross Abase

ð27Þ

where Afin;cross is the total cross sectional area of fins and Abase is the total base area. This geometry factor reflects how much area is occupied and blocked by fins. 3.3. Size optimization of straight channel SC heat sinks are considered as the benchmark for the TO heat sinks. To make a fair comparison, size optimization is performed by 3D CFD analysis. Channel and fin width both vary from 0.3 mm to 0.8 mm with an interval of 0.1 mm, and inlet velocity varies from 0.05 m/s to 0.35 m/s with an interval of 0.05 m/s. Fin height is fixed as 1.2 mm and heat flux is 25 W=cm2 . 252 cases are studied and the best size parameters and operating velocities are selected based on the Pareto Front of the graph of junction temperature versus pumping power, as shown in Fig. 10. It is found that the fin width is always 0.3 mm for all the Pareto Front data points while the channel width varying from 0.8 mm to 0.3 mm. For high T j;av e and low P pump regions, the Pareto Front data points are large channel widths with low flow rates. While for low T j;av e and high P pump regions, the Pareto Front data points are all 0.3 mm channel width with different velocities. This also shows the advantage of having small channels in boosting heat transfer capability in microchannel heat sinks.

95 85 75 65 55

Pareto Front 45 35 0.0001

0.001

0.01

0.1

Fig. 10. Size optimization data and the Pareto Front of SC heat sinks.

4. Experimental study In order to verify the simulation results, a TO heat sink and a SC heat sink were manufactured by CNC machining and experimentally tested. This section will discuss the experimental work in detail. 4.1. Experiment test rig A water loop, shown in Fig. 11, was built for testing microchannel heat sinks. Deionized water is stored in a reservoir tank1 with a capacity of 3.8 L. a variable-speed gear pump2 is used to drive the fluid. After the pump, water goes through a compact liquid to liquid heat exchanger3 which is connected to a water bath4. Flow rate is measured by a flow sensor5 before the test section. At the test section, pressure drop is monitored by a differential pressure transducer6. Water inlet and outlet temperature are measured by two T-type thermocouples7 inserted in the deep plenums right before and after the heat sink. The heat sink temperature is read from 6T-type thermocouples8 at different locations of the heat sink base. Four cartridge heaters9 powered by a high DC voltage power supply unit10 supply heat to the test section. At the downstream, water goes through a filter11 and radiator12 before flowing back to the reservoir. Besides, both flow sensor and differential pressure transducer are connected to another DC power supply unit13 to provide 12 V excitation power. Data from all the sensors, e.g., thermocouples, the flow sensor and the differential pressure transducer, are taken by a data acquisition system14. Details of equipment and sensors used in the flow loop are summarized in Table 5. 4.2. Test section The test section design is the most critical part of the experimental work. It is designed to house copper heat sinks and facilitate data measurements. The exploded view of the 3D drawing of the test section is shown in Fig. 12 and its actual apparatus is shown in Fig. 13. It is mainly composed of four parts: top cover, top housing, bottom housing and the copper block. Microchannels are machined out at the top surface of the copper block with an footprint area of 45  32 mm. Its geometry is created by repeating the 2 mm wide topology optimized structure 16 times. For the SC heat sink manufactured, its fin width is 0.3 mm and its channel width is 0.5 mm. Six 0.7 mm diameter holes are drilled into the middle of the heat sink at the height level of 2 mm lower than the bottom surface of microchannels (fin base). They are distributed evenly in the streamwise direction and each side has 3 thermocouples. There is an extended part at the height level of 4 mm lower than the bottom surface of microchannels, working as a flange to fix the copper block’s position at the height direction when being inserted into the top housing. Below the flange, four 38 mm long and 6.5 mm diameter holes are drilled from the bottom for housing cartridge heaters. The 3D drawings of the copper blocks are shown in Fig. 14, focusing on their top fin geometries. The top housing is made of PEEK and designed to tightly hold the copper block. It has two deep plenums at the inlet and outlet portion to facilitate measuring water inlet and outlet temperature and pressure difference. The shallow plenums before and after the copper block are 23 mm long each, modelled as the inlet and outlet extension parts in simulations. Six 0.7 mm diameter holes are drilled at the corresponding positons of the top housing to allow thermocouples to go through. The top cover is made of Polycarbonate and works as a lid. There are sixteen 3.2 mm diameter holes for M3 screw to go though and tighten the top cover with the top

11

S. Zeng, P.S. Lee / International Journal of Heat and Mass Transfer 142 (2019) 118401

∆

Test Section

Fig. 11. Schematic diagram of the water flow loop.

Table 5 List of equipment and sensors. No.

Device

Brand

Model No.

Specification

1 2 3 4 5

Reservoir tank Gear pump Heat exchanger Water bath Flow sensor

McMaster-Carr Cole-Parmer BHE Manufacturing Ronneby Thermo Scientific McMillan

6 7

Differential pressure transducer Inlet and outlet thermocouple

Setra Omega

Portable wide mouth ASME 01 EA pressure tank 75211-70 AlfaNova 14-10H (A21, A21) Type 003-2617 002-4276 G106-5D G112-8 2301001PD2F2DAC TJ40-CPSS-116G-5

8

Heat sink thermocouple

Omega

TJ40-CPSS-020G-2

9 10 11 12 13 14

Cartridge heater Heater power supply Filter Radiator Sensor power supply DAQ system

Watlow ITECH Swagelok Thermatron ITECH Agilent

E1J35-L12 IT6726H SS-4FW-15 735SPC2A01 IT6322 34972A

3.8L – – – 50–500 mL/min 0.5–5 L/min 0–1 PSI type-T D = 1/16 in. type-T D = 0.02 in. 240 V 250 W 300 V/20 A/3000 W 15 lm – 0–30 V –

Top cover

Top housing

Copper block Bottom housing

Fig. 12. Exploded view of the 3D drawing of the test section.

housing. There are also eight 5.3 mm diameter holes for M5 screw to go through and tighten the top cover, top housing and bottom housing together. The two small holes above the inlet and outlet plenums are intended for removing trapped air inside the plenums, in case of any. A groove is cut on both top cover and top housing to hold an O-ring for sealing. The bottom housing is made of Teflon. Its standing-up part could support the copper block and ensure that it goes to the top housing. In order to prevent leakage, silicone gel and silicon sheet is applied on the contact surfaces of the stand-up part of the bottom housing and copper block and contact surfaces of the top housing and bottom housing as an addition to the tight fit between the copper block and top housing. 4.3. Experiment procedure After the test section is assembled, a widely used heat loss characterization is first performed [47]: The test section is heated at a constant heating power until the temperature measurement from the 6 thermocouples in the copper block is stale. The temperature and heating power are recorded. Then we increase the heating

12

S. Zeng, P.S. Lee / International Journal of Heat and Mass Transfer 142 (2019) 118401

Outlet TC

Heat sink TC

Inlet TC

Water outlet

Water inlet

Outlet Pressure

Heat sink TC

Inlet Pressure

Fig. 13. Actual apparatus of the test section (TC – Thermocouples).

Fig. 14. 3D drawings of (a) topology optimized (b) straight channel heat sinks tested in experiments.

power and wait the copper block temperature to be stable again. The process is repeated up to 90 °C with an interval of 10 °C. At the same time, the environment temperature is also recorded. The heat loss is assumed to be a linear function with respect to the temperature difference of the copper block and environment. The experiment data and best fit curve are shown in Fig. 15. The best fit curve function is:

q_ loss ¼ 0:1743ðT HS;av e  T en Þ þ 0:0846

ð28Þ

14 12 2

10

= 0.9909

8

̇

6 4 2 0 0

10

20

30

40

50

60

Fig. 15. Best fit curve of the heat loss function.

70

80

Thus, the total heat dissipated through water could be calculated as:

q_ water ¼ V heater Iheater  q_ loss

ð29Þ

Next, all valves in the flow loop are opened and the pump is activated. The pump speed is adjusted to get the required flow rates. Voltage of the cartridge heater power supply is tuned to get the desired q_ water . Then water bath temperature is adjusted to ensure the water inlet temperature as 25 °C. The copper block temperature may gradually change in this process, resulting in change of heat loss q_ loss . Thus, the voltage of the cartridge heater power supply needs to be finely adjusted accordingly. Finally, stable condition is reached when: (1) q_ water deviates from the required by less than 0.1 W; (2) T water;in ¼ 25  0:2 °C; (3) all thermocouples fluctuate within 0:2 °C in 1 min and none of them are continuously increasing or decreasing; (4) pressure drop measurement fluctuates within 10 Pa and is not continuously increasing or decreasing; (5) flow rate deviates from the setting point by less than 10 mL/min and is not continuously increasing or decreasing. Then data is acquired for another 1 min with 1 Hz frequency. The flow rate is then changed to the next test condition and the process repeats. Cartridge heater power supply unit used in this paper is programmable. It is controlled by Matlab and its output voltage and current are also read through Matlab. Flow sensor, differential pressure transducer and thermocouples are all connected to an Agilent data acquisition module. Then Agilent Bench Logger 3 is

13

S. Zeng, P.S. Lee / International Journal of Heat and Mass Transfer 142 (2019) 118401

used to monitor the real-time readings and Matlab is used to take the final stable readings and post process.

A standard uncertainty analysis [49] is performed. The error for a derived parameter is defined as

4.4. Experimental data reduction

sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 2 2 @f @f @f df ¼ ð dxÞ þ ð dyÞ þ ð dzÞ þ


@x @y @z

The heat sink temperature is measured at 2 mm below the microchannel bottom surface (fin base). For temperature at 4 mm below the microchannel bottom surface, namely junction temperature, is calculated based on 1D heat conduction assumption.

T j;av e ¼ T HS;av e þ

q_ water HHS 2kHS AHS

ð37Þ

where f is the derived parameter, x, y, z, . . .are the measured parameters, and d means the absolute error. The range of uncertainties calculated for various parameters reported in this study are listed in Table 7.

ð30Þ

where T HS;av e is the arithmetic average of the 6 thermocouples inserted into the copper block. Water mass flow rate is calculated as

_ water ¼ q@T V_ water m water;in

ð31Þ

where V_ water is the volume flow rate taken from the flow sensor. The measured pressure difference from the pressure differential transducer is the combined pressure loss due to contraction from the inlet deep plenum to the shallow plenum DPc1 , friction in the inlet shallow plenum DPf 1 , contraction from the inlet shallow plenum to microchannels DP c2 , friction in the microchannels DP f 2 , expansion from microchannels to the outlet shallow plenum DPe2 , friction in the outlet shallow plenum DPf 3 and expansion from the outlet shallow plenum to the outlet deep plenum DP e1 .

DPmeasured ¼ DPc1 þ DPf 1 þ DPc2 þ DPf 2 þ DPe2 þ DPf 3 þ DPe1

ð32Þ

The inlet and outlet shallow plenums are included in the simulation domain as inlet and outlet extension parts, so only DPc1 and DPe1 are extra pressure loose in experiments comparing to simulations. According to [48], the contraction and expansion loss could be calculated as:

K 1  c1 DPc1 ¼ qfl v 2sp;in  v 2dp;in þ q v2 2 2 fl sp;in

ð33Þ

K 1  e1 DPe1 ¼ qfl v 2dp;out  v 2sp;out þ q v2 2 2 fl sp;out

ð34Þ

5. Performance comparison and discussion 5.1. 3D numerical analysis validation The experiment is designed for verifying results from 3D numerical analysis. The modified TO heat sink with 0.5 mm minimal channel width (TO 0.5) and a SC heat sink with 0.3 mm fin width and 0.5 mm channel width (SC 0.5) were manufactured and tested experimentally. Heat flux is kept as 25 W=cm2 as in the simulations while flow rates vary in the studied range. The comparison of 3D simulation and experimental results are shown in Fig. 16. The results, both junction temperature and pressure drop, agree with each other reasonably well for all cases. For the TO heat sink, maximum deviation of temperature and pressure drop between numerical and experimental results are 1.0 °C and 6.0%, respectively. For the SC heat sink, they are 0.97 °C and 5.2%, respectively. These errors are considered acceptable, thus, the 3D simulation results are taken as trustable data for further analysis. It could also be seen from Fig. 16(a) that a critical flow rate exists below which the SC heat sink exhibits better cooling capability while above which the TO heat sink is better. In terms of pressure drop, the TO heat sink is always lower than the SC heat sink, which promises the advantage of lower pumping power. The low pressure drop feature also corresponds to the topology optimization objective. 5.2. Comparison of TO and SC heat sinks

The uncertainties of various sensors used in experiments are listed in Table 6.

The experimentally validated 3D numerical analysis provides detailed data for comparing performance of TO heat sinks and SC heat sinks. As mentioned in section 2.4, we generated two topology optimized heat sink designs with different minimal channel width, 0.3 mm and 0.5 mm. Both will be compared with SC heat sinks of the same level of channel size constraints, respectively. Fig. 17 shows the comparison of TO heat sinks and SC heat sinks in terms of thermal resistance and pressure drop under various flow rates. The fin width of SC heat sinks is 0.3 mm, which is proved to be the optimized size in Section 3.3. SC 0.5 means that the channel width is 0.5 mm while SC 0.3 means 0.3 mm. It is found that SC 0.3 has the smallest thermal resistance while having much higher pressure drop penalty. TO 0.5 is better than SC 0.5 in thermal resistance under high flow rates while worse under low flow rates. Thermal resistance of TO 0.3 is always smaller than that of TO 0.5 in the flow range studied due to the smaller channel size.

Table 6 Uncertainties of the sensors.

Table 7 Range of uncertainties.

where K c1 is chosen as 1, and K e1 is calculated as

K e1 ¼ ð1 

Asp 2 Þ ¼ 0:8638 Adp

ð35Þ

where Asp and Adp are the cross section areas of shallow and deep plenum. The velocities, v sp;in , v dp;in , v sp;out and v dp;out , are calculated by dividing V_ water with corresponding plenum cross sectional areas. Then the total heat sink pressure drop is calculated as

DP ¼ DPmeasured  DPc1  DPe1

ð36Þ

4.5. Experimental uncertainties

Sensor T-type thermocouples Differential pressure transducer Flow sensor (G106-5D) Flow sensor (G112-8) Heater power supply voltage and current

Uncertainty

Parameter

Range of uncertainty

±0.5 °C 0.25% full scale 1% full scale 0.5% full scale 0.01% and 0.1%

Average junction temperature T j;av e Pressure drop DP Volume flow rate V_ water

0.51 °C  18:75 Pa 5 mL/min below 500 mL/min and  25 mL/min above 500 mL/min  0:15%

Heating power to water q_ water

14

S. Zeng, P.S. Lee / International Journal of Heat and Mass Transfer 142 (2019) 118401

80

1200 TO Num 25 TO Exp 25 SC Num 25 SC Exp 25

(a) 75 70

TO Num 25 TO Exp 25 SC Num 25 SC Exp 25

(b) 1000 800

65 60

600

55

400

50 200

45 40

0 50

150

250

350

450

550

650

50

150

250

350

450

550

650

Fig. 16. Comparison of experimental and 3D numerical results in (a) average junction temperature (b) pressure drop.

4000

0.15 TO 0.5 0.13

SC 0.5

3500

TO 0.3

3000

TO 0.5 SC 0.5 TO 0.3 SC 0.3

SC 0.3

0.11

2500 2000

0.09

1500

0.07

1000 0.05

500 0

0.03 50

250

450

650

50

250

450

650

Fig. 17. Comparison of thermal resistance and pressure drop of TO heat sinks and SC heat sinks.

In terms of pressure drop, both TO 0.5 and TO 0.3 are smaller than SC 0.5. It could also be seen that TO heat sinks are less sensitive to the channel width constraint compared to the SC heat sinks. Throughout the considered flow range, the difference in thermal resistance between SC 0.5 and SC 0.3 is around 0.015 °C/W while it is less than 0.01 °C/W for TO heat sinks. This observation is more significant for pressure drop difference, in which SC 0.3 is around 3 times as SC 0.5 while TO 0.3 is only around 1.1 times as TO 0.5. It is because that the defined minimum channel width has a global influence on SC heat sinks. On the contrary, it only influences several local spots in TO heat sinks as their channel is not uniform. The thermal and hydraulic characteristics of TO heat sinks are mostly dependent on the overall topology. In term of overall performance, considering both pressure drop penalty and heat transfer enhancement, heat sinks are compared in a junction average temperature versus pumping power graph. In Fig. 18, the SC heat sink data is taken from the Pareto Front of Fig. 10, which is the best data sets. They are composed of multiple sizes, not necessarily channel width of 0.3 mm or 0.5 mm. SC 0.5 means best data points with a geometry constraint of channel width larger than or equal to 0.5 mm and SC 0.3 means those with channel width larger or equal to 0.3 mm. Data points of TO heat sinks are on the left side of the SC heat sinks, which proves better overall performance of the TO heat sinks. For example, to cool the heat sink to the same average junction temperature of 45 °C, SC 0.3 requires 0.0065 W pumping power while TO 0.3

will only consume 0.004 W, which is 38.5% less; SC 0.5 requires 0.011 W pumping power while TO 0.5 only consume 0.0054 W, which is 50.9% less. On the other hand, if all the heat sinks operate under pumping power of 0.01 W, TO 0.3 could cool the chip to 41.4 °C while SC 0.3 is 42.8 °C; TO 0.5 could cool to 42.1 °C while

70 TO 0.5 SC ≥ 0.5 TO 0.3 SC ≥ 0.3

65 60 55 50 45 40 0

0.005

0.01

0.015

Fig. 18. Comparison of average junction temperature under various pumping powers.

S. Zeng, P.S. Lee / International Journal of Heat and Mass Transfer 142 (2019) 118401

SC 0.5 is 45.1 °C. It has to be emphasized that the geometries of TO heat sinks are fixed while the SC heat sinks compared here are composed of various channel sizes. 5.3. Flow and thermal characteristic of the TO heat sinks In this section, we will analyse in detail about the reasons for the performance difference presented. SC 0.3 shows lower thermal resistance under same flow rates. This is due to its larger heat transfer area with water. As shown in Table 8, there is a great difference between TO heat sinks and SC heat sinks in wetted surface area Awet . SC 0.5 is 1.37 times of TO 0.5 and SC 0.3 is 1.64 times of TO 0.3. To analyse their heat transfer characteristics without the huge influence of heat transfer area difference, their heat transfer coefficients are calculated:

h¼

q_ water Awet ðT wall;av e  T water;av e Þ

where T wall;av e is taken from area-weighted average of the wetted surface in 3D simulations, and T water;av e is the algebraic average of inlet and outlet water temperature. Fig. 19 shows the heat transfer coefficients of TO and SC heat sinks under various flow rates. It could be seen that the heat transfer coefficient of TO heat sinks increase almost linearly with the increase of the water flow rate. However, the improvement of heat transfer coefficient for SC heat sinks by increasing the water flow rate is limited under high velocities. This trend makes heat transfer coefficients of SC heat sinks slightly better than that of TO heat sinks in low velocity conditions while tremendously smaller in high velocity conditions. This is because of the complex microchannels in TO heat sinks which provide the heat transfer enhancement potential. As shown in Fig. 20, the velocity contour of the SC heat sink shows a simple and undisturbed flow field across the whole channel. The domain showed here is composed of half width fin and half width channel. The hydraulic boundary layer gradually devel-

Table 8 Geometry specifics of heat sinks.

Awet ðmm2 Þ BR

TO 0.5

TO 0.3

SC 0.5

SC 0.3

3828.61 25.7%

3985.52 28.9%

5248.80 37.5%

6518.40 50.0%

11400 TO 0.5 SC 0.5 TO 0.3 SC 0.3

10400 9400 8400 7400 6400 5400 4400 3400 50

250

450

650

Fig. 19. Heat transfer coefficients under various flow rates.

15

ops from the channel inlet to the outlet and the velocity at the channel center gradually increases. The flow fields of the TO heat sinks are quite different: sectional fins break hydraulic boundary layers continuously and guide the flow to turn and go through secondary channels instead of directly escaping from the main channel. The secondary channel flow diverges fluid from one streamwise main channel and injects it to the next streamwise main channel, contributing to flow mixing. From the enlarged view of Fig. 20(c) and (e), it could be found that a low flow rate region exists at the back of every piece of sectional fins, which is a vortex generated due to flow separation and local adverse pressure field. Secondary flow channels are blocked by the vortexes, and the vortex regions are larger in TO 0.5 than that in TO 0.3. Due to this effect, the increase in channel width does not necessarily lead to enhancement of secondary flow. Fig. 21 shows the temperature contours of heat sinks studied. In the SC 0.5 case, it is clear that its temperature boundary layer gradually becomes thicker along the flow stream. Differently, in the TO heat sinks, temperature boundary layers are reinitiated at the leading edge of each piece of sectional fins which greatly promotes convection heat transfer. On the other hand, there is clear temperature stratification in the SC heat sink that temperature is low in the middle of channel and high in the vicinity of the fin. Due to better flow mixing seen in the velocity contour, this stratification does not exist in the TO heat sinks, which provides lower fluid temperature in the vicinity of fins in the downstream. Comparing Fig. 21 (b) and (c), there is a region with higher temperature after each sectional fin for TO 0.5, especially at the downstream part. This is because of larger after-fin vortexes in TO 0.5 which deteriorate its heat transfer performance. In terms of pressure drop characteristic, as discussed in the last section, TO heat sinks are always lower than SC heat sinks though there are more turns, vortexes and secondary channels in TO heat sinks. Their low pressure drop is mainly due to their smaller contacting surface areas between fluid and solid in TO heat sinks, which results in smaller surface friction force. Another geometry factor that could unveil the reason is the blockage ratio (BR). Generally, larger BR means bigger flow resistance through the channels. As shown in Table 8, the blockage ratios of SC heat sinks are much larger than TO heat sinks, especially the BR of SC 0.3 is almost two times of TO 0.3. Thereforce, the extremely higher pressure drop penalty of SC 0.3 is reasonable. Lastly, as this study used similar method as that used to optimize air-cooled heat sink in [34], it is interesting to compare the characteristics of the two optimized heat sink for different coolants. TO 0.3 is used for comparison here since its modification from original optimization results is less than TO 0.5. Moreover, both are compared under similar Reynolds number of 550, which corresponds to inlet velocity of 0.149 m/s and 0.3 m/s for air and water heat sink, respectively. Their streamlines at the middle height of the channel is shown in Fig. 22. It is obvious that the water heat sink has more sectional fins with more spots of local flow perturbation while the air one prefers larger fins and less sections. A mainstream without encountering fins exists in the water heat sink, labeled as undisturbed mainstream in Fig. 22(b). This is beneficial for reducing pressure drop. However, it could deteriorate heat transfer as it would become a bypass for water if the channel width is too large. The optimized result is a trade-off between the pressure drop and heat transfer. Such undisrupted mainstream is not observed in the air heat sink. Another observation is that the air heat sink has two large unfinned regions (Region 1 and Region 2 in Fig. 22(a)), and even in the multi-stage optimization process [34], only two thin fins (Fin 1 and Fin 2 in Fig. 22(a)) were generated there and the two unfinned regions were maintained. Air gathers in these two regions then flows to the downstream. TO 0.3 obviously does not have such

16

S. Zeng, P.S. Lee / International Journal of Heat and Mass Transfer 142 (2019) 118401

(m/s) (a) (b)

(c)

(d)

(e)

Fig. 20. Velocity contours of (a) SC 0.5 (b) TO 0.5 (c) enlarged view of TO 0.5 (d) TO 0.3 (e) enlarged view of TO 0.3 (Note: the vertical scales shown are different for the SC and TO).

(℃) (a) (b) (c)

Fig. 21. Temperature contours of (a) SC 0.5 (b) TO 0.5 (c) TO 0.3 (Note: the vertical scales shown are different for the SC and TO).

Normalized Velocity

(a)

Fin 1

Region 1

10 mm

(b)

120 mm

Fin 2

Region 2

Undisturbed mainstream

2 mm 45 mm Fig. 22. Streamlines of (a) TO air heat sink (b) TO water heat sink 0.3 at the middle height of channel plane around Reynolds number of 550 (the velocity is normalized against their own locally highest velocity, showing v =v max ).

17

S. Zeng, P.S. Lee / International Journal of Heat and Mass Transfer 142 (2019) 118401

large unfinned regions. This could be because air conductivity is smaller than water, thus it needs larger space for flow mixing to make temperature uniform. 5.4. Further discussion on topology optimization 5.4.1. Accuracy of the 2D model and influence of hf l To verify the developed 2D thermofluid model described in Section 2.2, temperature and pressure drop results from this model for TO heat sinks and SC heat sinks are compared with 3D simulation results. From Fig. 23(a), it could be seen that the average junction temperature prediction from the 2D model fits with the results of 3D model reasonably well. The deviation for SC 0.5, SC 0.3 and TO 0.3 is all less than 0.5 °C. Larger deviation, 1.7 °C at the maximum, is seen at the higher flow rate conditions for TO 0.5. For pressure drop prediction, as shown in Fig. 23(b), values from the 2D model are always smaller than that of the 3D model and the discrepancy becomes larger at high flow rates. This is due to the friction force at top and bottom of channels could not be counted in 2D model and channels are taken as infinite high. However, the 2D model is able to predict the trend of pressure drop increasing with increase of flow rates. It could also provide right orders of pressure drop magnitude of different geometries (the 2D model predicts the order as SC 0.3 > SC 0.5 > TO 0.3 > TO 0.5, which is the same as that from 3D model simulation). This level of accuracy is enough for being used in topology optimization. It is also due to this discrepancy of pressure drop magnitude of the 2D model that 3D simulation is employed to verify performance of TO heat sinks. This good agreement in temperature prediction is mainly attributed to the choice of hfl as listed in Table 1. It is calculated from 3D simulation of SC 0.5 but applied for all the designs investigated here. As a comparison, Table 9 presented the hfl values calculated from 3D numerical analysis for all the designs investigated. It could be seen that hfl gradually increases with the increase of inlet velocities for all the four designs. TO heat sinks have a faster increasing speed which leads to larger discrepancy of hfl between TO and SC heat sinks under higher velocities. However, this discrepancy is still small in terms of the sensitivity of the 2D model with respect to the applied hfl values. To illustrate this point, 2D models for the four designs are evaluated at the inlet velocity of 0.2 m=s while hfl changes in the range from 1000 to 15,000 W=ðm2
KÞ. Then the average junction temperature prediction error is evaluated as the difference from the 3D numerical results, which is expressed as

78

TO 0.5 3D TO 0.5 2D SC 0.5 3D SC 0.5 2D TO 0.3 3D TO 0.3 2D SC 0.3 3D SC 0.3 2D

(a)

73 68 63

Table 9 hfl values of different heat sink designs under different inlet velocities.

v in (m/s)

hfl (W=ðm2
KÞ)

0.05 0.1 0.15 0.2 0.25

%T j;error ¼

SC 0.5

SC 0.3

TO 0.5

TO 0.3

3544.5 3984.6 4357.8 4668.4 4940.2

4363.9 4645.0 4944.9 5230.8 5495.1

3454.7 4299.4 5195.3 6175.1 7157.7

3761.7 4700.8 5569.2 6554.7 7574.0

ðT j;av e;2D  T water;in Þ  ðT j;av e;3D  T water;in Þ T j;av e;3D  T water;in

As shown in Fig. 24, the maximum absolute error for all four designs is 3.6 °C in the large range of hfl studied, which is 16.9% of relative error. In the real hfl variation range of the four designs under 0.2 m=s (4668.4 to 6554.7 W=ðm2
KÞ), the maximum absolute error is even smaller as 1.6 °C (7% of relative error). Thus, as long as the adopted hfl is reasonably close to the accurate value, the 2D model presents satisfying accuracy. For example, to achieve an error smaller than 9% for all four designs, hfl could be chosen from 4000 to 10,000 W=ðm2
KÞ. Using the values for SC heat sinks at similar length scales is sufficient to fulfill requirement for the accuracy of hfl .

Fig. 24. Average junction temperature prediction errors of the 2D model when using different hfl values.

3000

TO 0.5 3D TO 0.5 2D SC 0.5 3D SC 0.5 2D TO 0.3 3D TO 0.3 2D SC 0.3 3D SC 0.3 2D

(b)

2500 2000

58

ð38Þ

1500

53

1000

48 500

43 38

0 50

150

250

350

450

550

650

50

150

250

350

450

Fig. 23. Comparison of average junction temperature (a) and pressure drop (b) prediction in 2D and 3D simulation.

550

650

18

S. Zeng, P.S. Lee / International Journal of Heat and Mass Transfer 142 (2019) 118401

To further investigate the influence of hfl on the topology optimization results, case 0.1 m/s–47.43 °C is implemented under different hfl values. As shown in Fig. 25, the results from hfl ¼ 2000 W=ðm2
KÞ to 6000 W=ðm2
KÞ are similar. It is not until the large deviation of hfl ¼ 500 W=ðm2
KÞ and 20,000 W=ðm2
KÞ that obvious difference in topologies is observed. This means that the topology optimization is not sensitive to the selected hfl value and using the value calculated from the straight channel is adequate. It has to be mentioned that full 3D topology optimization is still attractive even the 2D optimization is proved to be successful as it may lead to novel structures especially in the fin height direction. 5.4.2. Formation of sectional fins in topology optimization The TO heat sinks are composed of multiple pieces of fins with irregular layouts. During optimization, the fins are not formed piece by piece. Instead, larger pieces of fins are generated in the first several steps then secondary channels grow out to cut them into sectional fins. Table 10 shows the change of design variable fields of case 0.1 m/s–47.43 °C at different steps of topology optimization with parameter continuation strategy in Table 2. At the first two steps, the penalization of fluid penetrating through solid or intermediate area is small (low a and high qf ), so large area of the design domain is occupied by the intermediate and solid area to take advantage of non-physical high heat transfer performance in those penetration area. As the increase of Brinkman penaliza-

tion, intermediate material will be greatly reduced at steps 3 and 4. But the blockage topology begins to form. Though it is not totally blocking the flow section, some dead end channels are formed, so that all the fluid inside the dead end is forced to penetrate through the solid areas. In step 5, sectional channels are gradually formed to act as an escape path for the penetration fluid to decrease the pressure drop caused by penetration. However, the dead-endkind of blockage is still remained because of its non-physical contribution to heat transfer and optimization is stuck in this topology. In the following 3 steps, the small channels become slightly larger and the general topology remains unchanged. The last two steps are conducted mainly for getting sharp and clear boundaries between fluid and solid domains. The forced penetration in the optimization results will boost its heat transfer performance, which will be deteriorate with explicit boundaries in validation simulation. This issue is also reported in [22,50] and it is a main issue in density based topology optimization, which calls for further refinement. It also has to be mentioned that the topology evolution process described here is dependent on the values of parameters used in continuation strategy (Table 2). The present parameter set results in the best TO heat sink performance among hundreds of trials. It is different from the continuation strategy used in our previous air-cooled heat sink optimization, in which qh and qk were first increased and then decreased because of the severe blockage problem during its optimization.

Fig. 25. Optimization results of 0.1 m/s–47.43 °C under various hfl values.

Table 10 Design variable fields at different steps of topology optimization.

Continuation Step 1 2 3 4 5 6 7 8

Design variable field

Color legend

S. Zeng, P.S. Lee / International Journal of Heat and Mass Transfer 142 (2019) 118401

6. Conclusion In this paper, we explored the ability of density based topology optimization method in designing planar water-cooled heat sinks. Non-conventional channel layouts and fin structures with good performance are generated through topology optimization. 3D numerical simulation is implemented to compare their performance with SC heat sinks and study their flow and thermal characteristics. Experiments are designed and performed to validate numerical results. The following conclusions could be drawn: 1. Topology optimization is a powerful tool and good supplementation for designing microchannel heat sinks. The generated heat sink structures are non-conventional and show better performance in comparison with size optimized SC heat sinks. Quantitatively speaking, over 50% of pumping power saving could be reached with the same cooling requirements. 2. The layouts and shapes of sectional fins in TO heat sinks create more complex flow phenomena, e.g., flow mixing, secondary flow and boundary initiation, contributing to the enhancement of convection heat transfer. On the other hand, lower blockage ratio results in lower pressure drop comparing to SC heat sinks. Both factors guarantee better overall performance of TO heat sinks. 3. The 3D numerical simulation method presented in this paper is proved to be accurate in predicting performance of microchannel TO heat sinks. The pressure and temperature prediction errors comparing to experimental data are less than 6% and 1 °C, respectively. 4. The original structures from topology optimization contain very small channels that are out of the capability of CNC machining. In such cases, post processing of manually enlarging channels is required. It is possible to avoid such inconvenience by involving length scale control techniques reported in literature into our current topology optimization scheme.

Declaration of Competing Interest Author declares that there is no conflicts of interest. References [1] J.B. Marcinichen, J.A. Olivier, V. de Oliveira, J.R. Thome, A review of on-chip micro-evaporation: experimental evaluation of liquid pumping and vapor compression driven cooling systems and control, Appl. Energy 92 (2012) 147– 161, https://doi.org/10.1016/j.apenergy.2011.10.030. [2] J. Marschewski, R. Brechbühler, S. Jung, P. Ruch, B. Michel, D. Poulikakos, Significant heat transfer enhancement in microchannels with herringboneinspired microstructures, Int. J. Heat Mass Transf. 95 (2016) 755–764, https:// doi.org/10.1016/j.ijheatmasstransfer.2015.12.039. [3] A. Capozzoli, M. Chinnici, M. Perino, G. Serale, Review on performance metrics for energy efficiency in data center: The role of thermal management, Lect. Notes Comput. Sci. (Including Subser. Lect. Notes Artif. Intell. Lect. Notes Bioinformatics). 8945 (2015) 135–151, https://doi.org/10.1007/978-3-31915786-3_9. [4] S. Alkharabsheh, B. Ramakrishnan, B. Sammakia, Pressure drop analysis of direct liquid cooled (DLC) rack, in: 16th IEEE Intersoc. Conf. Therm. Thermomechanical Phenom. Electron. Syst., 2017, pp. 815–823. [5] D.B. Tuckerman, R.F.W. Pease, High-performance heat sinking for VLSI, IEEE Electron Dev. Lett. 2 (1981) 126–129, https://doi.org/10.1109/EDL.1981.25367. [6] T. Kishimoto, S. Sasaki, Cooling characteristics of diamond-shaped interrupted cooling fin for high-power LSI devices, Electron. Lett. 23 (1987) 456–457, https://doi.org/10.1049/el:19870328. [7] C.A. Rubio-Jimenez, S.G. Kandlikar, A. Hernandez-guerrero, Numerical Analysis of Novel Micro Pin Fin Heat Sink With Variable Fin Density, 2 (2012) pp. 825– 833. [8] Y.J. Lee, P.S. Lee, S.K. Chou, Enhanced thermal transport in microchannel using oblique fins, J. Heat Transf. 134 (2012) 101901, https://doi.org/10.1115/ 1.4006843.

19

[9] Y. Sui, C.J. Teo, P.S. Lee, Y.T. Chew, C. Shu, Fluid flow and heat transfer in wavy microchannels, Int. J. Heat Mass Transf. 53 (2010) 2760–2772, https://doi.org/ 10.1016/j.ijheatmasstransfer.2010.02.022. [10] Z.G. Mills, A. Warey, A. Alexeev, Heat transfer enhancement and thermalhydraulic performance in laminar flows through asymmetric wavy walled channels, Int. J. Heat Mass Transf. 97 (2016) 450–460, https://doi.org/10.1016/ j.ijheatmasstransfer.2016.02.013. [11] A.F. Al-Neama, N. Kapur, J. Summers, H.M. Thompson, An experimental and numerical investigation of the use of liquid flow in serpentine microchannels for microelectronics cooling, Appl. Therm. Eng. 116 (2017) 709–723, https:// doi.org/10.1016/j.applthermaleng.2017.02.001. [12] P. Li, Y. Luo, D. Zhang, Y. Xie, Flow and heat transfer characteristics and optimization study on the water-cooled microchannel heat sinks with dimple and pin-fin, Int. J. Heat Mass Transf. 119 (2018) 152–162, https://doi.org/ 10.1016/j.ijheatmasstransfer.2017.11.112. [13] M.T. Al-Asadi, A. Al-damook, M.C.T. Wilson, Assessment of vortex generator shapes and pin fin perforations for enhancing water-based heat sink performance, Int. Commun. Heat Mass Transf. 91 (2018) 1–10, https://doi. org/10.1016/j.icheatmasstransfer.2017.11.002. [14] A.L. Goh, K.T. Ooi, Nature-inspired inverted fish scale microscale passages for enhanced heat transfer, Int. J. Therm. Sci. 106 (2016) 18–31, https://doi.org/ 10.1016/j.ijthermalsci.2016.03.010. [15] I.A. Ghani, N.A.C. Sidik, N. Kamaruzaman, Hydrothermal performance of microchannel heat sink: the effect of channel design, Int. J. Heat Mass Transf. 107 (2017) 21–44, https://doi.org/10.1016/j.ijheatmasstransfer.2016.11.031. [16] S.G. Kandlikar, H.R. Upadhye, Extending the heat flux limit with enhanced microchannels in direct single-phase cooling of computer chips, in: Semicond. Therm. Meas. Manag. IEEE Twenty First Annu. IEEE Symp. 2005, 2005, https:// doi.org/10.1109/STHERM.2005.1412152. [17] W. Escher, B. Michel, D. Poulikakos, A novel high performance, ultra thin heat sink for electronics, Int. J. Heat Fluid Flow. 31 (2010) 586–598, https://doi.org/ 10.1016/j.ijheatfluidflow.2010.03.001. [18] W. Escher, T. Brunschwiler, B. Michel, D. Poulikakos, Experimental Investigation of an ultrathin manifold microchannel heat sink for liquidcooled chips, J. Heat Transf. 132 (2010) 081402, https://doi.org/10.1115/ 1.4001306. [19] C.S. Sharma, G. Schlottig, T. Brunschwiler, M.K. Tiwari, B. Michel, D. Poulikakos, A novel method of energy efficient hotspot-targeted embedded liquid cooling for electronics: an experimental study, Int. J. Heat Mass Transf. 88 (2015) 684– 694, https://doi.org/10.1016/j.ijheatmasstransfer.2015.04.047. [20] M.P. Bendsøe, N. Kikuchi, Generating optimal topologies in structural design using a homogenization method, Comput. Methods Appl. Mech. Eng. 71 (1988) 197–224, https://doi.org/10.1016/0045-7825(88)90086-2. [21] T. Borrvall, J. Petersson, Topology optimization of fluids in Stokes flow, Int. J. Numer. Methods Fluids 41 (2003) 77–107, https://doi.org/10.1002/fld.426. [22] 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, https://doi.org/10.1016/j. ijheatmasstransfer.2016.05.013. [23] J. Alexandersen, O. Sigmund, K.E. Meyer, B.S. Lazarov, Design of passive coolers for light-emitting diode lamps using topology optimisation, Int. J. Heat Mass Transf. 122 (2018) 138–149, https://doi.org/10.1016/j. ijheatmasstransfer.2018.01.103. [24] G.H. Yoon, Topological design of heat dissipating structure with forced convective heat transfer, J. Mech. Sci. Technol. 24 (2010) 1225–1233, https:// doi.org/10.1007/s12206-010-0328-1. [25] A.A. Koga, E.C.C. Lopes, H.F. Villa Nova, C.R.D. Lima, E.C.N. Silva, Development of heat sink device by using topology optimization, Int. J. Heat Mass Transf. 64 (2013) 759–772, https://doi.org/10.1016/j.ijheatmasstransfer.2013.05.007. [26] E.M. Dede, Optimization and design of a multipass branching microchannel heat sink for electronics cooling, J. Electron. Packag. 134 (2012) 041001, https://doi.org/10.1115/1.4007159. [27] Y. Sato, K. Yaji, K. Izui, T. Yamada, S. Nishiwaki, An optimum design method for a thermal-fluid device incorporating multiobjective topology optimization with an adaptive weighting scheme, J. Mech. Des. 140 (2018), https://doi.org/ 10.1115/1.4038209 031402. [28] S.B. Dilgen, C.B. Dilgen, D.R. Fuhrman, O. Sigmund, B.S. Lazarov, Density based topology optimization of turbulent flow heat transfer systems, Struct. Multidiscip. Optim. (2018), https://doi.org/10.1007/s00158-018-1967-6. [29] X. Zhao, M. Zhou, O. Sigmund, C.S. Andreasen, A ‘‘poor man’s approach” to topology optimization of cooling channels based on a Darcy flow model, Int. J. Heat Mass Transf. 116 (2018) 1108–1123, https://doi.org/10.1016/j. ijheatmasstransfer.2017.09.090. [30] J. Asmussen, J. Alexandersen, O. Sigmund, C.S. Andreasen, A ‘‘poor man’s” approach to topology optimization of natural convection problems, Struct. Multidiscip. Optim. 59 (2019) 1105–1124, https://doi.org/10.1007/s00158019-02215-9. [31] X. Qian, E.M. Dede, Topology optimization of a coupled thermal-fluid system under a tangential thermal gradient constraint, Struct. Multidiscip. Optim. 54 (2016) 531–551, https://doi.org/10.1007/s00158-016-1421-6. [32] J.H.K. Haertel, G.F. Nellis, A fully developed flow thermofluid model for topology optimization of 3D-printed air-cooled heat exchangers, Appl. Therm. Eng. 119 (2017) 10–24, https://doi.org/10.1016/j.applthermaleng.2017.03.030.

20

S. Zeng, P.S. Lee / International Journal of Heat and Mass Transfer 142 (2019) 118401

[33] J.H.K. Haertel, K. Engelbrecht, B.S. Lazarov, O. Sigmund, Topology optimization of a pseudo 3D thermofluid heat sink model, Int. J. Heat Mass Transf. 121 (2018) 1073–1088, https://doi.org/10.1016/j.ijheatmasstransfer.2018.01.078. [34] S. Zeng, B. Kanargi, P.S. Lee, Experimental and numerical investigation of a mini channel forced air heat sink designed by topology optimization, Int. J. Heat Mass Transf. 121 (2018) 663–679, https://doi.org/10.1016/j. ijheatmasstransfer.2018.01.039. [35] K. Lee, Topology Optimization of Convective Cooling System Designs, University of Michigan, MI, 2012. [36] J. Alexandersen, Efficient topology optimisation of multiscale and multiphysics problems, 2016. doi:10.13140/RG.2.2.15890.45761. [37] T. Van Oevelen, Optimal Heat Sink Design for Liquid Cooling of Electronics, 2014. https://mail.google.com/mail/u/0/?shva=1%5Cnpapers3://publication/ uuid/4DFCDED7-E75B-461E-911A-B466327B3A4D. [38] A. Gersborg-Hansen, O. Sigmund, R.B. Haber, Topology optimization of channel flow problems, Struct. Multidiscip. Optim. 30 (2005) 181–192, https://doi.org/ 10.1007/s00158-004-0508-7. [39] J. Alexandersen, Topology Optimisation for Coupled Convection Problems, Technical University of Denmark, Lyngby, 2013. [40] K. Svanberg, MMA and GCMMA – Fortran versions March 2013 Krister Svanberg Considered optimization problem, (2013) pp. 1–23. [41] J.V. Carstensen, J.K. Guest, Projection-based two-phase minimum and maximum length scale control in topology optimization, Struct. Multidiscip. Optim. 58 (2018) 1845–1860, https://doi.org/10.1007/s00158-018-2066-4. [42] F. Wang, B.S. Lazarov, O. Sigmund, On projection methods, convergence and robust formulations in topology optimization, Struct. Multidiscip. Optim. 43 (2011) 767–784, https://doi.org/10.1007/s00158-010-0602-y.

[43] M. Zhou, B.S. Lazarov, F. Wang, O. Sigmund, Minimum length scale in topology optimization by geometric constraints, Comput. Methods Appl. Mech. Eng. 293 (2015) 266–282, https://doi.org/10.1016/j.cma.2015.05.003. [44] B.S. Lazarov, F. Wang, O. Sigmund, Length scale and manufacturability in density-based topology optimization, Arch. Appl. Mech. 86 (2016) 189–218, https://doi.org/10.1007/s00419-015-1106-4. [45] O.B. Kanargi, P.S. Lee, C. Yap, A numerical and experimental investigation of heat transfer and fluid flow characteristics of a cross-connected alternating converging-diverging channel heat sink, Int. J. Heat Mass Transf. 106 (2017) 449–464, https://doi.org/10.1016/j.ijheatmasstransfer.2016.08.057. [46] D.R. Lide, Handbook of Chemistry and physics, 84th ed.,., CRC Press, 2004, 10.1136/oem.53.7.504. [47] M. Law, P.S. Lee, A comparative study of experimental flow boiling heat transfer and pressure characteristics in straight- and oblique-finned microchannels, Int. J. Heat Mass Transf. 85 (2015) 797–810, https://doi.org/ 10.1016/j.ijheatmasstransfer.2015.01.137. [48] R.D. Blevins, Applied Fluid Dynamics Handbook, Van Nostrand Reinhold, Company Inc., 1984. [49] J.R. Taylor, An Introduction to error analysis: the study of uncertainties in physical measurements, 2nd ed., University Science Books, 1997. [50] T. Matsumori, T. Kondoh, A. Kawamoto, T. Nomura, Topology optimization for fluid-thermal interaction problems under constant input power, Struct. Multidiscip. Optim. 47 (2013) 571–581, https://doi.org/10.1007/s00158-0130887-8.