Since the shape of the heat sink is unknown before performing optimization, the boundary between the heat sink and the air cannot be distinguished in the initial design domain. Imagining the heat dissipation material as a viscous fluid, this problem can then be described by the fluid heat transfer equation:
$$\rho {c}_{\text{P}}\varvec{u}\bullet \nabla T+\nabla \bullet \varvec{q}=Q+{Q}_{\text{P}}+{Q}_{\text{v}\text{d}}$$
1
$$\varvec{q}=-k\nabla T$$
2
where ρ is density, cP is specific heat, u is velocity, q is heat flux, k is thermal conductivity, Q is heat generation of heat source, QP is pressure work, and Qvd is viscous dissipation. \(\)The convective heat transfer term \(\rho {c}_{p}u\bullet \nabla T\), QP, and Qvd in Eq. 1 are ignored in this work for simplicity.
The variable density method is used to study the topological configuration of the structure, that is, the material density of the node is used as a variable. The density of nodes ranges from 0 to 1, with 1 being completely solid and 0 being voids filled with air. Using the variable density method, the optimization problem of material layout is reduced to binary integer variables. This paper uses Solid Isotropic Microstructures with Penalization (SIMP) interpolation, and the physical quantities in the design domain can be defined as:
$$\rho \left(\varvec{x}\right)={\rho }_{\text{s}}\left[\frac{{\rho }_{0}}{{\rho }_{\text{S}}}+\left(1-\frac{{\rho }_{0}}{{\rho }_{\text{S}}}\right){\rho }^{p}\right]k\left(\varvec{x}\right)={k}_{\text{S}}[\frac{{k}_{0}}{{k}_{\text{S}}}+(1-\frac{{k}_{0}}{{k}_{\text{S}}}\left){\rho }^{p}\right]$$
3
where x is the coordinate of the design variable. ρS and kS are the density and thermal conductivity, respectively. ρ0 and k0 are the density and thermal conductivity, respectively. ρ(x) and k(x) are the density and thermal conductivity of a node, respectively. p is the penalty factor, which penalizes the intermediate value of the material density of a node, usually between 3 and 5. The larger the p, the higher the stability of the solution; while the excessive penalty on the intermediate value will lead to less branches and thus the degraded heat dissipation performance. Thus p = 3 is adopted in this work.
The problem of this work is a standard nonlinear optimization problem (Powell, 1994):
\(\text{m}\text{i}\text{n}\text{i}\text{m}\text{i}\text{z}\text{e} f\left(x\right)={\int }_{{\Omega }}{k(\nabla T)}^{2}\text{d}{\Omega }\) | (4) |
\(\text{s}\text{u}\text{b}\text{j}\text{e}\text{c}\text{t} \text{t}\text{o} 0\le {\int }_{{\Omega }}\rho \left(\text{x}\right)\text{d}{\Omega }\le \gamma M\) | |
\(x\in \text{X}=\{{\text{x}}_{1},{\text{x}}_{2},{\text{x}}_{3},\cdots {\text{x}}_{\text{N}}\}\) | |
\(0\le \gamma \le 1\) | |
The objective function Eq. 4 in this optimization problem is to minimize the temperature variation in the design domain Ω, which is equivalent to minimizing the thermal resistance of the heat sink. The restriction is that the total amount of solid in the design domain should be below a predetermined value γM, where M is the total mass if the design domain is fully filled with the solid, and γ is the preset maximum volume fraction of the heat sink in the design domain. γ is taken as 0.4, that is, the maximum volume of solid will not exceed 40% of the design domain.
A widely used optimization method of moving asymptotes (MMA) is used in this work (Svanberg, 1993). In each iteration step, the equation system to be solved is approximated to make it closer to the target solution. The maximum number of iterations is 1000, and the residual is set to 0.001, below which is convergence.
In order to make the simulated model be manufacturable, the mesh size should no smaller than the resolution of the 3D printing, i.e. 1.2 mm. On the other hand, the complexity of the design must be limited. Otherwise, infinitesimal features will appear in the resulted model. Helmholtz-type differential equation is employed to filter too small features (Lazarov & Sigmund, 2011):
$${\theta }_{\text{f}}={R}_{\text{m}\text{i}\text{n}}^{2}{\nabla }^{2}{\theta }_{\text{f}}+{\theta }_{\text{c}}$$
5
where θc is the design variable before filtering, and θf is the design variable after filtering, Rmin is the filter radius that should to be no less than the mesh size. The filtered design variables may have a large number of intermediate values at the solid boundary. To make the boundary sharper, the hyperbolic tangent projection is used to reduce these physically insignificant regions (Datta et al., 2019):
$$\theta =\frac{\text{t}\text{a}\text{n}\text{h}\left(\beta \right({\theta }_{\text{f}}-{\theta }_{{\beta }}\left)\right)+\text{t}\text{a}\text{n}\text{h}\left(\beta {\theta }_{{\beta }}\right)}{\text{t}\text{a}\text{n}\text{h}\left(\beta \right(1-{\theta }_{{\beta }}\left)\right)+\text{t}\text{a}\text{n}\text{h}\left(\beta {\theta }_{{\beta }}\right)}$$
6
where θ is the output material volume factor, θβ is the projection point, and Β is the projection slope. This work uses θβ = 0.3, B = 6. All of the parameters used in simulation are listed in Table 1.
Table 1
Parameters used for the computation
Parameter | Symbol | Value |
Density of heat sink material | ρS | 3.132 g/cm3 |
Air density | ρ0 | 1.204 kg/m³ |
Thermal conductivity of the material | kS | 113.412 W/(m·K) |
Thermal conductivity of air | k0 | 0.024 W/(m·K) |
Penalty factor | p | 3 |
Maximum material volume fraction | γ | 0.4 |
Filter radius | Rmin | 7.625×10− 6 m |
Projection point | θβ | 0.3 |
Projection slope | β | 6 |