Sage Journals: Discover world-class research

Abstract

In this paper, two most prevalent topological optimisation approaches namely Density and Level set method are applied to a three dimensional heat sink design problem. The relative performance of the two approaches is compared in terms of design quality, robustness and computational speed. The work is original as for the first time it demonstrates the relative advantages and disadvantages for each method when applied to a practical engineering problem. It is additionally novel in that it presents the design of a convectively cooled heat sink by solving full thermo-fluid equations for two different solid-fluid material sets. Further, results are validated using a separate computational fluid dynamics study with the optimised designs are compared against a standard pin-fin-based heat sink design. The results show that the Density method demonstrates better performance in terms of robustness and computational speed, while Level-set method yields a better quality design in terms of final objective value.

Keywords

Topology optimisation density method level set method computational fluid dynamics validation heat sink design

Introduction

Thermal management is a key challenge in modern microelectronics system design due to ever increasing levels of miniaturisation, integration and operating frequency which result in significantly higher power densities.¹ The thermal design strategies adopted in order to meet these challenges include increasing utilisation of forced convection, and use of liquid cooling solutions.^2,3 Apart from effective cooling methods, an effective design method is also required to meet the thermal management challenges. Topology optimisation (TO) is an effective, optimisation-based design method, wherein optimal shapes are obtained from the given design domain meeting the specified set of constraints.⁴

The two most popular methods used in TO are the Density Method (DM) and the Level-set method (LSM).⁵ Other TO methods are namely Topology derivative method, Phase field method and Evolutionary approaches. In DM, the material density is used as a design variable, which is defined by the optimisation algorithm and its value directly indicates the phase (solid/void) of any given cell. The DM approach is most widely used in conjunction with the Method of Moving Asymptote (MMA) optimisation algorithm.^6,7

Level sets are implicit functions which are modelled to represent material interface and by advecting them in the decreasing direction of objective, optimum shapes are obtained in LSM.^8,9 The primary disadvantages of the density-based TO method are that the no slip condition is not strictly imposed on the solid walls and the interface between the solid and fluid may not be crisply captured due the presence of grey cells.⁵ In contrast, LSM provides a crisp interface with no grey region.

The DM approach has been applied to fluid flow problems by a number of researchers. Borrvall and Petersson¹⁰ pioneered TO involving Stoke’s flows. Olesen et al.¹¹ subsequently extended the study to consider full steady Navier–Stokes (NS) flows using the FEMLAB software. TO of thermo-fluidic problems started with Dede,² who optimised the liquid cooling channels of a rectangular domain with volumetric heat source. As material properties were not interpolated, the solid region created in the optimisation had zero thermal conductivity. Yoon¹² carried out the design of a heat dissipating structure subjected to forced convection and for the first time he interpolated thermal conductivity and other relevant material properties with respect to design variables. Thereby the resulting solid regions had the non-zero thermal conductivity.

Following this, many works have been published on heat sink optimisation, the most notable being by Koga et al.³ and Burger et al.¹³ The first work on natural convection cooled heat sink was carried out by Alexandersen et al.,^14,15 he optimised heat sink designs for various Grashof numbers by fully solving the thermo-fluidic governing equations. Haertel and Nellis¹⁶ optimised the air side surface of dry cooled power plant condensers by solving a steady state thermo-fluidic model with fully developed flow.

TO of fluid flow problems, starting from Stokes flow to NS flow are demonstrated using LSM. Challis and Guest¹⁷ studied TO of Stokes flow and Zhou and Li¹⁸ studied the TO of NS flow using variational LSM. Kreissl and Maute¹⁹ carried out TO of NS flows through Level set with Extended finite element method (XFEM)-based geometry mapping and he clearly brought out the advantages of LSM over DM for fluid problems, namely pressure diffusion across the solids and inefficient no-slip imposition on solid walls. Deng et al.²⁰ extended LSM to steady NS flow subjected to body forces.

Yamada et al.²¹ solved generic design dependent heat transfer problems through the LSM with Ersatz projection method and Yaji et al.²² applied LSM to a liquid cooled heat sink problem similar to the work of Koga et al.³ and Dede.² This work employed Tikhonov-based regularisation to enable qualitative control of geometric complexity. Coffin and Maute,²³ carried out TO of cooling device by approximating convective fluxes through Newton’s Law of cooling through LSM with XFEM mapping. The results obtained were better than the DM and the LSM with Ersatz mapping but the fluid flow equations are not solved and approximations used in lieu. Subsequently Coffin and Maute,²⁴ solved natural convection TO problems by solving full flow equations using an LSM with XFEM approach. Coffin also optimised a three dimensional (3D) heat sink subjected to steady low Grashof number natural convection.

Although considerable work is done on heat sink design using each of the methods, no work is done yet to assess their relative performance for 3D heat sink design, which is critical for choosing a suitable method for industrial applications. In this study a 3D heat sink subjected to steady laminar forced convection is designed using both DM and LSM to allow the performance of the two methods to be directly compared. This paper is organised as follows. The next section describes motivation for this work, Section ‘Heat sink design using the density method’ describes heat sink design using DM, which includes problem formulation, numerical implementation and results. Section ‘Heat sink design using the level-set method (LSM)’ describes heat sink design using LSM. Computational fluid dynamics (CFD)-based validation of optimised heat sinks against standard pin-fin heat sink is given in Section ‘Computational fluid dynamics-based validation’. Discussion about the performance of two methods, their comparison and validation results are given in Section ‘Discussion’ and conclusions are given in the last section.

Motivation

Although the DM and LSM approaches are becoming increasingly widely adopted for TO of variety of problems, their application to industrial heat transfer problems is limited. In particular, there is little, or no, work assessing the relative benefits and restrictions of these approaches. This is particularly the case in the 3D heat sink design problem forming the focus of this work. In this study, both DM and LSM are formulated in the same environment (COMSOL 5.1)²⁵ and their performance is assessed in terms of design quality (final objective value), robustness and computational cost. Furthermore, the design quality is assessed through a CFD study with results contrasted against a standard pin-fin heat sink. The relative performance is studied at differing thermal conductivity ratios. As such, the study aims to provide an insight into the applicability and effectiveness of DM and LSM for practical engineering design problems.

Heat sink design using the density method

The DM approach to TO is the most widely adopted and the mostly researched. In this method the relevant physical parameters of the problem (e.g. Young’s modulus for structural mechanics) are modelled as function of a design variable which is modelled to take values 0 or 1. MMA is the most widely used optimiser for optimising the design variables and it has been implemented in COMSOL software.

The various steps involved in TO by DM are depicted in Figure 1. The TO is carried out in COMSOL 5.1 by combining the optimisation module with relevant physics module in a coupled manner.

Figure 1.

Flow chart illustrating DM procedure.

At the start of the TO problem some initial value for design variable is assumed in all mesh elements within the domain. Then physics of the problem is solved and the solution is used to calculate the sensitivity and the objective value. Using the sensitivity, optimiser runs and gives a new distribution of density. The difference between the new density distribution and the earlier distribution is calculated and if the difference is significant then the physics solver and optimiser loop is repeated until convergence is achieved.

Density method formulation

The governing equations for this forced convective heat sink design study are given below. An artificial frictional force term is used to differentiate the solid and fluid materials.

ρ (\nabla \cdot u) = 0

(1)

ρ (u . \nabla u) = - \nabla p + \nabla . {µ {\nabla u + {(\nabla u)}^{T}}} - α u

(2)

ρ C p (u . \nabla T) = \nabla . (k \nabla T)

(3)

where ‘α’ is the effective impermeability, and it is zero in the fluid domain and takes higher value (10⁵) in case of solid domain. Along with impermeability, thermal properties like thermal conductivity (k), specific heat capacity (Cp) and density (ρ) are varied depending on the grid cells’ design variable (γ) value. In DM, value of design variable determines whether the element is fluid (γ=0) or solid (γ=1).¹¹ The interpolation of thermal properties are carried out as per¹² and they are given in Table 1. The subscript ‘s’ stands for solid and subscript ‘f’ stands for fluid property.

α (γ) = α_{m a x} γ^{3}

(4)

Table 1.

Thermal properties interpolation formula in DM.

Name	Expression
k	(k_s−k_f)×γ³ + k_f
Cp	(Cp_s−Cp_f)×γ³+ Cp_f
ρ	(ρ_s−ρ_f)×γ³ + ρ_f

DM: density method.

Thermal properties of solid and fluid used in this study are given in Table 2. Minimisation of thermal compliance is the objective of this optimisation. The problem can be stated as below

Objective function: min \int_{Ω} k (γ) * {(\nabla T)}^{2} d Ω

(5)

Subjected to

Table 2.

Parameter values used for DM-based 3D heat sink design.

Property	Value
k_s	40 (W/(mK)) (k_f/k_s=0.001)0.4 (W/(mK)) (k_f/k_s =0.1)
ρs	8920 (kg/m³)
Cp_s	385 (J/(kg×K))
k_f	0.04 (W/(mK))
ρf	1000 (kg/m³)
Cp_f	4184 (J/(kg×K))
η_f	1.002e-3 (Pa.s)

DM: density method.

Thermo-fluidic Governing equations (1) to (3)

Volume constraint: $\int_{Ω} γ d Ω \leq 0.25 \times V$

where ‘V’ is design domain volume. An ideal heat sink has to effectively transfer the heat throughout the design domain to keep the thermal compliance at minimum.

Density method computational details

The computational domain considered in this study is shown in Figure 2. The computational domain considered is one quadrant of the total domain, making use of symmetry boundary condition on the two sides to reduce computational costs. Note, due to limitation of COMSOL, symmetry heat flux boundary condition is applied based on thermal conductivity of fluid material (k_f) instead of with design material thermal conductivity (k(γ)). The design domain is a cube of length 0.1 m while the computational domain is of size 0.7 × 0.7 × 0.3 m. A constant heat flux is applied on the front corner of the bottom wall (10,000 W/m² at area of 0.01 × 0.01 m²). The upper surface of the computational domain is defined as an inlet and assigned a velocity of 4e-5 m/s while a pressure outlet condition is assigned on the two side walls, which are adjacent to symmetry condition. The lower surface is defined as zero flux except the heat flux boundary region. The flow is simulated for a Reynolds number of 8 at which the Prandtl number corresponds to 104.6. The volume fraction of solid material is constrained at 25%. The design domain is discretised with 35 × 35 × 35 hexahedral cells giving a total mesh size of 147,000 elements. Due to computational resource limitations further mesh refinement is not performed in design domain, but from the experience of pervious two dimensional (2D) optimisation studies we believe mesh independence is achieved. More details regarding setting up the density-based topological optimisation of 2D heat sink are provided in Santhanakrishnan et al.²⁶

Figure 2.

Computational domain of 3D heat sink design study.

Since gradient-based optimiser is used in this study, in order to find global optimum, the optimisation study has been carried out using different γ initial values, namely 0.10, 0.25, 0.45 and 0.55. Linear discretisation is used for both velocity and pressure along with streamwise diffusion stabilisation. Although the use of linear elements doesn’t fulfil the Babuska-Brezzi condition, usage of streamwise diffusion helps to circumvent this. Temperature and optimisation variable γ are also discretised linearly. The governing equations are solved in segregated manner with the linear system of equations solved using a GMRES solver. The optimisation is assumed converged if the change in objective value between consecutive iterations is less than 0.01.

Density method results

High conductivity solid case (k_f/k_s=0.001)

DM simulation for k_f/k_s=0.001 is carried out at initial gamma values of 0.10, 0.25, 0.45 and 0.55. Simulation for gamma initial values higher than 0.55 failed because higher gamma indicates higher solid volume in the domain and that leads to low fluid permeability and flow stability problems. Each optimisation study resulted in a slightly different optimal shape, indicating the presence of many local minima in the problem domain. By comparing objective value and maximum temperature in the domain of different solutions, better optimal solution is identified and that corresponds to γ=0.25 solution (Figure 3).

Figure 3.

DM optimised heat sink for k_f/k_s=0.001.

It has to be noted that the solution contains some grey regions, hence in the figure γ of 0.6 is used as a threshold value. The green square surface at the bottom of the heat sink shape indicates the region where heat flux is applied. The full view (four quadrants) of this heat sink is shown in Figure 4 and temperature distribution within the design domain is shown in Figure 5.

Figure 4.

Full view of optimised heat sink (isometric, top view) by DM for k_f/k_s=0.001.

Figure 5.

Temperature contour in the design domain for optimised heat sink by DM method for k_f/k_s=0.001 (top view & bottom view).

Low conductivity solid case (k_f/k_s=0.1)

DM TO simulations for k_f/k_s=0.1 is also carried out at different γ initialisations. In heat transfer involving liquid metal cooling, this kind of conductivity ratio is possible. For example, in copper metal and gallium liquid cooling, this conductivity ratio is possible. The optimised shape nearly remains same for different γ initialisation runs. Gamma initial value of 0.55 yields minimum objective among the tested values. To obtain the optimal shape γ threshold value of 0.9 is used. The optimised heat sink and its full view are given in Figure 6 and corresponding temperature distribution in the design domain is shown in Figure 7. A convergence plot of the objective value for k_f/k_s=0.001 and k_f/k_s=0.1 are shown in Figure 8. The computational time given in Table 3, relates to ten real cores/twenty hyper-threaded cores on a Dual Xeon CPU cluster node.

Figure 6.

DM optimised heat sink for k_f/k_s=0.1 and its full symmetrised view.

Figure 7.

Temperature contour in the design domain for optimised heat sink by DM for k_f/k_s=0.1 (top & bottom view).

Figure 8.

Convergence of Objective value for k_f/k_s=0.001 and k_f/k_s=0.1.

Table 3.

Summary of DM results.

k_f/k_s	Thermal compliance (WK)	Maximum temperature (k)	Cumulative number of model evaluations	Computational time
0.001	4.146	378	216	33 h 7min
0.1	149.5	532	135	19 h 49 min

DM: density method.

Table 4.

Thermal properties interpolation formula in LSM.

Property	Symbol	Expression
Thermal conductivity	k	(k_s−k_f)×H + k_f
Specific heat capacity	Cp	(Cp_s−Cp_f)×H + Cp_f
Density	ρ	(ρ_s−ρ_f)×H + ρ_f

LSM: level-set method.

Heat sink design using the level-set method (LSM)

Level set based TO of 3D heat sink is carried out in COMSOL 5.1 following the works of Liu et al.²⁷ and Deng et al.²⁰ The shape and size of the computational domain are same as that used in DM-based heat sink design. The solution of physics and advection of the level set are considered in a coupled manner within the COMSOL.

Level-set method formulation

In this problem, positive Signed Distance Function (SDF) ( $ψ$ ) is considered to represent solid and negative SDF is considered to represent fluid (Figure 9). This is enforced by the Ersatz projection approach,⁹ using a Heaviside function.

ψ = {\begin{cases} = 0 \forall x \in \partial Ω (B o u n d a r y) \\ > 0 \forall x \in Ω^{+} (S o l i d r e g i o n) \\ < 0 \forall x \in Ω^{-} (F l u i d r e g i o n) \end{cases}

(6)

Figure 9.

Design domain and level set function.

The governing equations are the same as the ones used in the DM approach. Brinkman’s porosity term (α) is used to differentiate solid and liquid as modelled below

α = (α_{max} - α_{min}) \times H + α_{min}

(7)

where H is Heaviside function, which is equal to unity when level-set function (LSF) is positive, equal to zero when LSF is negative and it has smooth transition between the two levels in order to enable differentiability. The derivative of Heaviside function is δ function whose expression is also given below

H (ψ) = \frac{1}{2} + \frac{15}{16} (\frac{ψ}{h}) - \frac{5}{8} {(\frac{ψ}{h})}^{3} + \frac{3}{16} {(\frac{ψ}{h})}^{5}

(8)

δ (ψ) = \frac{15}{16 h} {(1 - {(\frac{ψ}{h})}^{2})}^{2}

(9)

where α_max =1e4 and α_min =0.01

At any point within the design domain, the thermal properties k, Cp and ρ are computed based on the values of Ψ and H as given in Table 4.

Objective function and constraints are exactly same as in DM (equation (5)). Hamilton Jacobi (HJ) equation is marched in time to convect the level set function in the decreasing direction of objective value. This is done by taking the velocity of convection equal to sum of shape sensitivity, Lagrange multiplier and volume constraint terms. Thermal compliance minimisation is a self adjoint problem for Stokes flow velocities, whose shape sensitivity is given below

HJ equation : \frac{\partial ψ}{\partial t} = V_{n} | \nabla ψ |

(10)

Shape sensitivity

V n = [k * ({(\nabla T)}^{2}) + λ + Λ (V o l u m e d i f f e r e n c e)]

(11)

where λ is the Lagrangian multiplier calculated through

λ = - \frac{\int_{Ω}^{} [k * ({(\nabla T)}^{2})] δ^{2} (ψ) | \nabla ψ | d Ω}{\int_{Ω}^{} δ^{2} (ψ) | \nabla ψ | d Ω}

(12)

Λ is the volume penalty factor, which needs to be suitably selected to ensure the volume constraint is met. This is achieved by trial and error method and the suitable value for this problem found to be −50. It should be noted that the Lagrange multiplier only preserves the area or it assumes that the initial level set distribution satisfies the area constraint. More details on procedure of LSM modelling can be found in the reference.²⁷

Note re-initialisation of level set function is not carried out in the present study. Although optimised shape will be less accurate in detail without re-initialisation, the overall shape of optimal can still be evaluated in this approach. This can be used to obtain quick first estimate of the topology optimised shape.

Level-set method computational details

The computational domain used for the study is identical to that used in the DM study (Figure 2). The design domain is discretised with 43 × 43 × 43 mesh cells. As the final optimum shape depends on the initial level set distribution, two different level set initialisations are tried, namely uniform sphere distribution (A) and cube distribution (B) and they are given in Figure 10.

Figure 10.

Initial Level set distributions uniform sphere (A) and cube distributions (B).

Reynolds number and Material properties used for this simulation are same as in the DM. COMSOL automatically selects suitable time step size for time marching HJ equation depending on the stability of the numerical system. Note that in this study new hole nucleation is enabled by extending the velocity field throughout the domain, rather than only around the zero level set boundaries.

Level-set method results

High conductivity solid case:

The run time of 3D coupled LSM optimisations are generally higher. The optimisation study with cube-like initialisation ran for 23 days on the same cluster node used to conduct the DM study. The final optimised shape for high conductivity solid with spherical LS initialisation is shown in Figure 11 and with cubic LS distribution is shown in Figure 12. Although both the shapes have similar order of thermal compliance value, spherical LS had more grey areas. Hence the shape obtained through cubical LS initialisation is the best optimised shape.

Figure 11.

Optimised shape for k_f/k_s=0.001 through LSM with initialisation A and the convergence plot.

Figure 12.

Optimised shape for k_f/k_s=0.001 through LSM with initialisation B.

The LSM optimised shape overall resembles like the optimised shape generated using the DM approach (Figure 4) with webs connecting the three outer edges of design domain with heat flux boundary. But the LSM shape differs in finer details with DM shape; LSM shape has many gaps or holes between the radial arms and further there are two smaller web like structures created in the LSM design, which were absent in DM design. The temperature distribution in the design domain for this case is shown in Figure 13. The same temperature scale which is used to plot DM result (Figure 5) is used, and the figure shows temperature is uniformly distributed throughout the domain and hence thermal compliance and maximum temperature are lower than the DM. Presence of grey cells will also contribute to uniform temperature distribution. To determine the magnitude of this effect, CFD validation has been performed on the optimised shapes.

Figure 13.

Temperature contour on design domain for optimised shape for k_f/k_s=0.001 through LSM with initialisation B.

Low conductivity solid case:

For the low conductivity solid case, optimisation run could not achieve the volume constraint of 25%. Volume decreased upto 43% then the simulation stagnated indicating a presence of local minimum there. The optimised shape obtained is given in Figure 14. It should be noted that, in both DM and LSM, optimised shape of high conductivity solid case has more grey cells than the low conductivity solid case.

Figure 14.

Optimised shape for k_f/k_s=0.1 through LSM with initialisation B.

Table 5 summarises the optimisation results of DM and LSM methods. For k_f/k_s = 0.001, LS yields a much lower objective value than the DM design, this indicates LSM has given better shape than DM, but amount of grey cells present in both methods, needs to be taken into account before concluding that.

Table 5.

Comparison of DM and LSM results (Volume fraction 0.25).

	Thermal compliance		Max temperature (K)
k_f/k_s	LSM	DM	LSM	DM
0.001	2.045	4.146	376	378
0.1	163.17	149.5	546.54	532

DM: density method; LSM: level-set method.

Computational fluid dynamics-based validation

It is necessary to validate the optimised shape obtained because (1) gradient-based optimisers which are prone to initialisation effect are used in DM, (2) LSM is sensitive to initialisation and (3) result of both the methods have grey cells and the threshold parameter for demarcating the solid from fluid regions is chosen by visual judgement rather than by scientific support. For the purpose of validation, we now compare the cooling effectiveness of DM and LSM optimised heat sinks with a standard heat sink through a CFD study.

The ‘standard’ heat sink is designed based on an article by Yang et al.²⁸ He has found an optimum pin-fin heat sink cooled by air impingement by Taguchi method. But since, the Reynolds number in the present optimisation study is very much lower than the Yang’s study, a uniform inter-fin spacing is selected in this study.

Validation of the high conductivity solid case (k_f/k_s=0.001)

From DM result to extract the heat sink shape γ threshold of 0.5 is used. The resulting heat sink shape has a material volume of 19%. Hence, standard heat sink is also designed to have material volume of 19%. In order to compare it with equivalent TO result, additional TO runs are conducted with volume constraint of 19% in both DM and LSM.

The DM heat sink geometry used for CFD study is relatively simplified to enable meshing. Figure 4 denotes the actual optimised shape, but in this unattached regions are removed and very thin plates attached with branches are removed to carry out the meshing and the simplified geometry is shown in Figure 16. A tetrahedral mesh is generated in the design domain, which has two different material domains, namely solid region created through optimisation (k = 40 W/m/k) and fluid created through optimisation (k = 0.04 W/m/k). In total 1.3 million tetrahedral elements were used to discretise the entire computational domain.

The pin-fin heat sink is designed in such a manner that it occupies 19% of the domain volume in order to correlate to the DM heat sink study. Each fin has square cross section of side 0.00703 m and a height of 0.1 m including the fin base of height 0.01m. Fin base size is 0.1 × 0.1 m. The inter-fin space is kept uniform at 0.0125 m and the domain meshed with 1.3 million tetrahedral cells to correlate to the DM CFD study. The conjugate heat transfer physics module is used to carry out the CFD study.

The LSM optimised shape is more complex geometry with many small surfaces and gaps. CFD simulation on the actual TO geometry is nearly impossible in COMSOL 5.1, with the current geometry import features. Hence a simplified geometry (optimised shape obtained with LS initialisation A) has been analysed. However, mesh and CFD setup are very similar to DM and pin-fin design.

From the CFD study, the thermal compliance of the design domain is computed for all the three heat sinks and they are compared (Table 6) against the value obtained during TO. Temperature contours of standard heat sink, DM and LSM shape are presented in Figures 15 to 17, respectively. The standard heat sink and DM designed heat sink have a more uniform temperature distribution than the LSM heat sink. The following points can be observed from the table.

Thermal compliance and maximum temperature of standard pin-fin heat sink and DM heat sink are of similar order but LSM is slightly higher.

Thermal compliance obtained in CFD study is higher than the compliance obtained during the TO study.

Maximum Temperature obtained through CFD study is lower than the temperature obtained during the optimisation study.

Figure 15.

Temperature distribution from CFD study on Standard heat sink for k_f/k_s=0.001.

Figure 16.

Temperature distribution from CFD study on DM heat sink for k_f/k_s=0.001.

Figure 17.

Temperature distribution from CFD study on LSM heat sink for k_f/k_s=0.001.

Table 6.

Validation of TO results for k_f/k_s=0.001 for volume fraction 19%.

	DM results	LSM results	CFD result of standard heat sink	CFD result of DM shape	CFD result of LSM shape
Thermal compliance (kgm²K/s³)	6.518	2.05	9.498	10.10	13.86
Maximum temperature (K)	383.9	378.58	313.89	315.89	314.68

TO: topology optimisation; DM: density method; LSM: level-set method; CFD: computational fluid dynamics.

The reason for the CFD study of the DM design reporting a higher objective value is due to the grey regions formed during the TO process. These grey cells have a relatively high conductivity and distribute heat effectively within the domain, thereby reducing the objective value during TO. Alexandersen et al.¹⁵ observed 20% difference in objective value between optimisation result and CFD result. The difference observed in the present study (54.9%) is higher than the value reported by Alexandersen, mostly because of the lower alpha max (1e5) value used during this simulation and use of ‘Solid Isotropic Material with Penalisation’ (SIMP) instead of ‘Rational Approximation of Material Properties’ (RAMP) penalisation method.

The surface area exposed to convective cooling is calculated for the three heat sinks as obtained in TO and the simplified shape used for CFD validation. Table 7 shows the comparison.

Table 7.

Surface area of different heat sinks for volume fraction 19%.

	Surface area (m²) (volume fraction 19%)
Heat sink	TO shape	CFD shape
Standard pin-fin	–	0.06525
Density based heat sink	0.05308	0.04427
LS based heat sink	0.0820	0.03907

CFD: computational fluid dynamics.

Heat sink with higher exposed surface area will have better cooling performance. The LSM optimised heat sink has 26% higher surface area than standard heat sink, hence this should have better performance than standard heat sink. However, geometry simplification done to enable the CFD meshing, lead to poor performance of LSM heat sink. Hence if CFD simulation is carried out on actual TO shape, they will certainly perform better than standard heat sink.

Figure 18 shows temperature gradient on the heat sink surfaces obtained through CFD study. Heat transfer is higher on the outer top edges of cubic domain when subjected to fluid injection from top. The DM heat sink has more solid on the outer edges, so as to decrease the temperature gradient, thereby minimising the thermal compliance. The LSM heat sink also has more solid on outer edge (Figure 12), but the simplified LSM shape used for CFD study does not have much solid near the top edges leading to poor performance.

Figure 18.

Temperature gradient comparison between Standard heat sink, DM & LSM heat sinks for k_f/k_s=0.001.

Validation of low conductivity solid case (k_f/k_s=0.1)

The DM heat sink geometry with gamma threshold value of 0.925 had a material volume of 25%, hence CFD simulation is conducted using this shape. It has to be noted that the DM result for the present case, has fewer grey cells than k_f/k_s=0.001 case. A 720,000 element tetrahedral mesh was generated over the computational domain.

The standard heat sink geometry is different from the one used for k_f/k_s=0.001 case. Its fin height is 0.08 including the fin base and the total material volume is 25% of design domain. The fin material has thermal conductivity of k_s=0.4W/m/K. From the CFD study, the thermal compliance of design domain is computed for both the heat sinks and they are compared (Table 8) against the value obtained during TO. Temperature contours of the DM heat sink and standard heat sink are given in Figure 19, which shows that the DM heat sink has effectively distributed the heat thereby decreasing the maximum temperature and thermal compliance contrary to the behaviour of the standard heat sink. The following points can be observed from the table;

Table 8.

Validation of DM TO results for k_f/k_s=0.1 for volume fraction 25%.

	DM results	CFD result of standard heat sink	CFD result of DM shape
Thermal compliance (kgm²K/s³)	149.5	248.12	157.74
Maximum temperature (K)	532	555.81	455.16

TO: topology optimisation; DM: density method; CFD: computational fluid dynamics.

Figure 19.

CFD-based temperature contours of Standard and DM heat sinks for k_f/k_s=0.1

The DM optimised heat sink performs much better than the standard heat sink. Optimised heat sink lowers the maximum temperature in the design domain by 100 K compared to standard heat sink.

As in the previous case, the objective value computed through CFD simulation of the DM heat sink is higher than the value obtained during TO. Maximum temperature computed through CFD simulation is lower than the temperature obtained during TO.

As like high conductivity solid case, the difference observed in thermal compliance between standard and DM heat sink can be explained from the layout of solid material in relation to temperature gradient. In this case, high temperature gradients are observed only near the heat source and the DM heat sink has more solid material near the heat source so as to decrease the temperature gradient value. Eventually, the thermal compliance decreases for the DM heat sink but since the standard heat sink has uniform material distribution throughout the design domain it could not decrease the thermal gradient and hence the thermal compliance.

Discussion

Based on the application of DM and LSM to the 3D heat sink design problem it can be observed that DM provided optimal solutions for both the material set problems considered. LSM provided a slightly better solution for the high conductivity ratio case but for low conductivity case the obtained solution is a poor local optimum which has not met the volume constraint. As such, DM should be considered to be more robust than LSM as it provides optimum designs for all considered problems.

Furthermore, the LSM run time (23 days) is significantly longer than DM (∼1.5 days). This is offset by the DM requirement for multiple differing γ initialisation runs to identify the global optimal solution whereas in LSM 1 or 2 intuitive initialisation runs can provide superior optimal solutions. The higher computational cost of LSM is due to coupled solution of thermo-fluidic and HJ equations.

Grey regions are observed in both the methods. Presence of grey cells will artificially understate the thermal compliance during TO. So to decide on best optimal shape, grey cell information has to be taken in to account. The volume of the grey regions is also dependent on the nature of the optimisation problem and its proximity to a convex optimisation problem. In LSM study, absence of re-initialisation lead to more grey cells, but it simplified the computation considerably.

Optimised shapes obtained through both methods have some disconnected regions. In a practical heat sink, those disconnected regions are not viable and hence they are removed during CFD validation. Implementing thin feature control mechanism in the optimisation method, will prevent their formation.

Alexandersen et al.¹⁴ interpolated the thermal properties using a convex natured RAMP relation and also followed continuation approach to reduce the grey regions and to reach global optimum. If a continuation approach is followed for α, then its value can be gradually increased to higher values (1e7) thereby enhancing the accuracy in modelling of solids. The LSM TO for k_f/k_s=0.001 case is to some extent similar to the 3D LSM TO carried out by Coffin and Maute.²⁴ The primary differences between this study are that Coffin and Maute work utilises a simple Newton’s law of cooling model for calculation of heat transfer, while here complete thermo-fluidic equations are solved.

Conclusions

TO of a 3D heat sink cooled by laminar forced convection is conducted using both the DM and the LSM. The density-based optimisation is carried out with MMA optimiser and Level set optimisation with Ersatz material mapping. Complete thermo-fluidic equations and HJ equations are solved in a fully coupled manner in Level set method. Two different types of heat sink materials were considered; one with high thermal conductivity and other with low thermal conductivity solid. The objective of the optimisation study was to minimise thermal compliance.

The optimised shapes obtained for the high conductivity solid resemble a web connecting the outer edges of cube with the heat source. The optimised shape is conceptually similar in both the methods but with some differences in the finer details. The LSM shape has multiple gaps keeping the overall shape same, thereby it has much higher surface area and lower thermal compliance than DM shape for the same material volume. For the low conductivity solid case, DM gave good results whereas LSM seems to reach a local minima, for the two different initialisations tried. Results indicate that DM demonstrates better performance in terms of robustness and computational speed, while LSM yields a better quality design.

The optimised shapes were validated through comparison of their CFD performance against the CFD result of standard pin-fin heat sink. The CFD validation of k_f/k_s=0.1 case shows that, TO heat sink preforms better than standard pin-fin heat sink. For k_f/k_s=0.001 case, optimised shapes are performing on equal level to standard pin-fin heat sink. Since the TO heat sink shapes are not directly amenable for CFD mesh generation, some geometry simplifications were required. The simplification reduced the surface area of the optimised shapes considerably. If simulations are performed on actual TO shapes (which has higher surface area), they will have superior performance than the standard pin-fin design.

Footnotes

Declaration of Conflicting Interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: The first author would like to express gratitude for financial support from the University of Greenwich for sponsoring this work as part of a PhD study.

References

Tummala

RR.

Packaging: past, present and future. In: Proceedings of the 6th international conference on electronic packaging technology, Shenzhen, China, 30 August- 2 September, 2005.

Dede

EM.

Multiphysics topology optimization of heat transfer and fluid flow systems. In: Proceedings of the COMSOL conference, Boston, October, 2009

3 Koga

Lopes

ECC

Nova

HFV

et al . Development of heat sink device by using topology optimization. Int J Heat Mass Transf 2013; 64: 759–772.

Bendsoe

Kikuchi

Generating optimal topologies in structural design using a homogenization method. Comput Methods Appl Mech Eng 1988; 71: 197–224.

Sigmund

Maute

Topology optimization approaches: a comparative review. Struct Multidisc Optim 2013; 48: 1031–1055.

Svanberg

Method of moving asymptotes – a new method for structural optimization. Int J Numer Meth Eng 1987; 24: 359–373.

Svanberg

A class of globally convergent optimization methods based on conservative convex separable approximations.

SIAM J Optim 2002; 12: 555–573.

Wang

Guo

A level set method for structural topology optimization. Comput Methods Appl Mech Eng 2003; 192: 227–246.

Allaire

Jouve

Toader

A-M.

Structural optimization using sensitivity analysis and a level-set method. J Comput Phys 2004; 194: 363–393.

10.

Borrvall

Petersson

Topology optimization of fluids in Stokes flow. Int J Numer Meth Fluids 2003; 41: 77–107.

11.

Olesen

Okkels

Bruus

A high-level programming-language implementation of topology optimization applied to steady-state Navier–Stokes flow. Int J Numer Meth Eng 2006; 65: 957–1001.

12.

Yoon

GH.

Topological design of heat dissipating structure with forced convective heat transfer. J Mech Sci Technol 2010; 24: 1225–1233.

13.

Burger

Dirker

Meyer

JP.

Three dimensional conductive heat transfer topology optimization in a cubic domain for the volume to surface problem. Int J Heat Mass Transf 2013; 67: 214–224.

14.

Alexandersen

Aage

Andreasen

et al . Topology optimization for natural convection problem. Int J Numer Methods Fluids 2013; 00: 1–23.

15.

Alexandersen

Sigmund

Aage

Large scale three dimensional topology optimization of heat sinks cooled by natural convection. Int J Heat Mass Transf 2016; 100: 976–891.

16.

Haertel

JHK

Nellis

GF.

A fully developed flow thermofluid model for topology optimization of 3D printed air cooled heat exchangers. Appl Therm Eng 2017; 119: 10–24.

17.

Challis

Guest

JK.

Level set topology optimization of fluids in Stokes flow. Int J Numer Meth Eng 2009; 79: 1284–1308.

18.

Zhou

A variational level set method for the topology optimization of steady-state Navier–Stokes flow. J Comput Phys 2008; 227: 10178–10195.

19.

Kreissl

Maute

Level set based fluid topology optimization using the extended finite element method. Struct Multidisc Optim 2012; 46: 311–326.

20.

Deng

Liu

et al . Topology optimization of steady Navier–Stokes flow with body force. Comput Methods Appl Mech Eng 2013; 255: 306–321.

21.

Yamada

Izui

Nishiwaki

A level set based topology optimization method for maximizing thermal diffusivity in problems including design dependent effects. J Mech Des 2011; 133.

22.

Yaji

Yamada

Kubo

et al . A topology optimization method for a coupled thermal-fluid problem using level set boundary expressions. Int J Heat Mass Transf 2015; 81: 878–888.

23.

Coffin

Maute

Level set topology optimization of cooling and heating devices using a simplified convection model. Struct Multidisc Optim 2016; 53: 985–1003.

24.

Coffin

Maute

A level-set method for steady-state and transient natural convection problems. Struct Multidisc Optim 2016; 53: 1047–1067.

25.

www.comsol.com. Accessed between October 2016 to February 2017.

26.

Santhanakrishnan

Tilford

Bailey

On the application of topology optimisation techniques to the thermal management of microelectronics systems. In: Proc. of 17th EurosimE conference, 18-20 April 2016, Montepellier, France.

27.

27. Liu

Korvink

Huang

Structure topology optimization: fully coupled level set method via FEMLAB. Struct Multidisc Optim 2005; 29: 407–217.

28.

Yang

Peng

Hsu

HT.

Numerical optimization of pin-fin heat sink with forced cooling. Int J Electr Comput Energ Electron Commun Eng 2013; 7.

Performance assessment of density and level-set topology optimisation methods for three dimensional heat sink design

Abstract

Keywords

Introduction

Motivation

Heat sink design using the density method

Density method formulation

Density method computational details

Density method results

Heat sink design using the level-set method (LSM)

Level-set method formulation

Level-set method computational details

Level-set method results

Computational fluid dynamics-based validation

Validation of the high conductivity solid case (kf/ks=0.001)

Validation of low conductivity solid case (kf/ks=0.1)

Discussion

Conclusions

Footnotes

Declaration of Conflicting Interests

Funding

References

Validation of the high conductivity solid case (k_f/k_s=0.001)

Validation of low conductivity solid case (k_f/k_s=0.1)