Topology Optimization Design of Micro-Channel Heat Sink by Considering the Coupling of Fluid-Solid and Heat Transfer

: To investigate the effect of the target weight coefﬁcient on the structure design of the micro-channel heat sink, an innovative method for the topology optimization design of micro-channel structures with different bifurcation angles is adopted. In this study, the improved interpolation function, density ﬁltering, and hyperbolic tangent projection methods are adopted to obtain a clear topological structure of the micro-channel heat sink. The heat transfer of the micro-channel heat sink under different bifurcation angles is compared. At the same time, the inﬂuence of the two different objective functions, heat transfer, and ﬂow energy consumption, is analyzed in the topology optimization of micro-channel heat sinks. The results show that when the bifurcation angle is 135 ◦ , both the heat transfer and the average outlet temperature of the micro-channel heat sink obtain the maximum value, and the heat transfer effect is the best. With the increase of the heat transfer weighting coefﬁcient, the distribution of solid heat sources in the main channel increases, and the reﬁnement of the branch channels also increases. On the other hand, although the heat transfer effect of the micro-channel heat sink is the best, the corresponding ﬂow energy consumption is larger.


Introduction
With the development of electronic components towards chip miniaturization and integration, the heat transfer problem caused by the high heat flux density will be a great challenge.The quality of the heat exchange problem will significantly affect the average temperature of the electronic components when they are working.Excessive average temperatures will reduce the service life, thereby reducing their thermal stability.In addition, the improvement of heat transfer performance will significantly reduce the energy consumption of the heat sink.Especially in the micro-channel heat sink, Tuekerman [1] first proposed a micro-channel water-cooled silicon integrated heat sink, which can take away 790 W/cm 2 of heat when the temperature difference between the inlet and outlet is 71 K.The micro-channel heat sink of this high-performance cooling device, which is conducive to energy conservation, low carbon, and environmental protection has become an important research field in recent years due to its huge heat transfer potential [2][3][4][5].
The traditional straight channel, V-channel, S-channel, and other high heat flux microchannel heat sinks cannot meet the heat transfer requirements of heat sinks in special environments due to its low heat transfer capacity and large flow resistance.In order to design a micro-channel structure with high heat transfer, Bejan et al. [6] proposed a configuration theory, and the designed fluid channel is a tree-like network.The microchannel structure designed by Li [7] and Peng [8], based on the configuration theory, can improve the heat transfer capacity and reduce flow resistance.Qu [9] studied the micropin-fin heat sink and the straight channel micro-channel heat sink, and found that the micro-pin-fin heat sink has a lower thermal resistance and better heat dissipation effect, but also has a higher pressure loss.Chein [10] et al. replaced the fluid in the micro-channel heat sink with nanofluids, and it was found that the performances were greatly improved when nanofluids were used as the coolants.Rahbarshahlan et al. [11] enhanced heat transfer by replacing some surfaces in the micro-channel with hydrophobic surfaces, and research shows that this method can reduce the flow rate of hydrophobic surfaces, reduce pumping power, and thus increase heat transfer.In addition, Rostamzadeh et al. [12] used a new multi-dimensional artificial characteristic-based (MACB) scheme, which is applied to different and complex geometries with heat transfer and forced convection, such as a heat exchanger, vortex tube, etc.However, this kind of structure design, based on previous traditional methods, design experience, experimental results, and other newly proposed methods has great limitations due to its small degree of freedom in design.Among the existing optimization techniques, the topology optimization method proposed by Bendsøe and Kikuchi [13] can be regarded as one of the most promising optimization tools, due to its large design freedom.These new methods could potentially integrate the topology optimization method to calculate the heat transfer topology optimization of heat sink.Borvall and Petersson [14] applied the variable density method topology optimization to the steady Stokes flow at a low Reynolds number for the first time, and established the basic model of variable density method topology optimization for fluid flow, which opened the research of scholars on fluid topology optimization.In terms of the topology optimization of thermal fluids, the variable density method has been used for the structural design of forced convection [15][16][17][18] and natural convection problems [19,20].Dede [21] used the commercial finite element analysis software COMSOL Multiphysics to optimize heat conduction and heat convection problems, and the objective function was expressed in terms of minimizing the average temperature and flow energy dissipation.Yoon [15] proposed a two-dimensional variable density method topology optimization to deal with the forced convection heat transfer problem, and the objective function was the minimum thermal compliance.Zhang et al. [22] designed the nanofluid-cooled heat sink based on topology optimization, and the research shows that reducing the diameter of nanoparticles and increasing the volume fraction of nanoparticles are conducive to improving the heat transfer performance of the heat sink.In most studies, topology optimization problems are solved with a single-objective function, such as minimizing the domain average temperature or minimizing thermal compliance [23,24].
Multi-objective optimization involves minimizing or maximizing multiple mathematical objective functions [25,26].When comparing single-objective optimization with multi-objective optimization, it is more difficult to produce a unique optimal solution; however, compromises can be made in different objective functions to make the overall objective as optimal as possible.Dong et al. [18] studied the influence of non-Newtonian fluid on heat transfer performance by taking minimizing fluid energy dissipation and maximizing recoverable heat power as objective functions under the condition of a given solid heat source.Based on the variable density topology optimization method, Koga et al. [16] took the weighting function of heat transfer and flow energy dissipation as the objective function, and studied the influence of different weighting coefficients and inlet boundary conditions on the heat sink topology.Lv and Liu [27] used the topology optimization method to design the flow channel of the micro-channel heat sink for electronic components.In the optimization, the maximization of heat transfer and the minimization of pressure drop were taken as objective functions, and the effects of different penalty factors, fluid inlet flow rates, and weighting coefficients on the topology structure of the heat sink were investigated.Ghasemi et al. [28] dealt with the heat transfer problem of a pin-fin heat exchanger based on topology optimization method, and optimized the geometric water tank with minimizing the thermal resistance and pressure loss as the objective function.The optimization results show that the optimized topology exhibits superior cooling performance with reduced pressure loss.
Based on the structure topology optimization design method of variable density, the flow equation and the convective heat transfer equation are coupled, and the optimization model of the micro-channel heat sink with an indefinite solid heat source is established.Using a new interpolation function, the performance of micro-channel heat sinks with different bifurcation angles is studied by numerical simulation.Considering the two objective functions, maximizing heat transfer and minimizing flow energy consumption, the topology optimization design of the micro-channel heat sink is carried out by using different target weight coefficients.

Physical Model
The design model of the fluid-structure coupled with the heat transfer of the microchannel heat sink is shown in Figure 1.The red area D represents the design domain.The bifurcation angle between the inlet and outlet of the heat sink system is set to θ.The parameter L is the inlet and outlet width of the heat exchange system, which is set to 1mm, and the inlet temperature T in is set to 273.15 K.In addition, the Reynolds number and Prandtl number at the inlet are 100 and 6.97 (the standard value of water), respectively.The fluid flow at the outlet directly flows into a large water tank without resistance, so the gauge pressure is set to 0 Pa.Γ in and Γ out represent the inlet and outlet boundaries of the design system, respectively, corresponding to the Dirichlet boundary and the Neumann boundary, and other boundaries are no-slip, adiabatic boundary conditions.
Based on the structure topology optimization design method of variable density, th flow equation and the convective heat transfer equation are coupled, and the optimization model of the micro-channel heat sink with an indefinite solid heat source is established Using a new interpolation function, the performance of micro-channel heat sinks wit different bifurcation angles is studied by numerical simulation.Considering the two ob jective functions, maximizing heat transfer and minimizing flow energy consumption, th topology optimization design of the micro-channel heat sink is carried out by using dif ferent target weight coefficients.

Physical Model
The design model of the fluid-structure coupled with the heat transfer of the micro channel heat sink is shown in Figure 1.The red area D represents the design domain .Th bifurcation angle between the inlet and outlet of the heat sink system is set to θ.The pa rameter L is the inlet and outlet width of the heat exchange system, which is set to 1mm and the inlet temperature Tin is set to 273.15 K.In addition, the Reynolds number and Prandtl number at the inlet are 100 and 6.97 (the standard value of water), respectively The fluid flow at the outlet directly flows into a large water tank without resistance, so th gauge pressure is set to 0 Pa.Γin and Γout represent the inlet and outlet boundaries of th design system, respectively, corresponding to the Dirichlet boundary and the Neumann boundary, and other boundaries are no-slip, adiabatic boundary conditions.

Governing Equations
In this study, the fluid flow problem is assumed to be the viscous incompressibl steady flow, and the fluid flow model can be expressed as: The governing equation of the heat transfer field can be written as: ( )

Governing Equations
In this study, the fluid flow problem is assumed to be the viscous incompressible steady flow, and the fluid flow model can be expressed as: where u, ρ, and p represent the flow velocity, density and pressure of the fluid, respectively, µ represents the dynamic viscosity of the fluid, and F represents the body force of the fluid flow.The governing equation of the heat transfer field can be written as: where T denotes the temperature of the design domain, k denotes thermal conductivity, c p denotes constant pressure heat capacity, and Q denotes the heat production per unit volume of the solid domain.

Topology Optimization Model Based on Variable Density Method
In the variable density topological optimization model of steady Navier-Stokes flow, the optimization problem is solved by adding a body force F, including the design variable γ, to the Navier-Stokes equations.Therefore, the body force in Equation (1) can be expressed as: where α(γ) is the resistance coefficient of artificial porous media, which is controlled by the design variable γ.In the variable density method, the design variable γ represents the distribution of materials in the design domain, which varies continuously from 0 to 1, and γ = 0, γ = 1 correspond to the solid and fluid domains in the design domain, respectively.Since the design variables between 0 and 1 have no practical significance, the introduction of the interpolation function γ p eliminates more intermediate design variables (intermediate density), which in turn reduces the grayscale of the final optimization result.
where the interpolation function γ p is defined as: where q is the penalty factor, and its influence on the interpolation function is shown in Figure 2a.The smaller the q value is, the greater the degree of aggregation of design variables to both ends of 0/1 by the interpolation function γ p , which is conducive to avoiding the generation of intermediate density.And when γ ∈ (0, 0.1) and q = 0.01, its influence on the interpolation function is most obvious, that is, it reduces the function value from 1 to 0 in the smallest γ change.Therefore, in this optimization process, q value is 0.01.

Topology Optimization Model Based on Variable Density Method
In the variable density topological optimization model of steady Navier-Stoke the optimization problem is solved by adding a body force F, including the design v , to the Navier-Stokes equations.Therefore, the body force in Equation (1) can pressed as: where () is the resistance coefficient of artificial porous media, which is contro the design variable .In the variable density method, the design variable  rep the distribution of materials in the design domain, which varies continuously from and  = 0,  = 1 correspond to the solid and fluid domains in the design dom spectively.Since the design variables between 0 and 1 have no practical significan introduction of the interpolation function   eliminates more intermediate design bles (intermediate density), which in turn reduces the grayscale of the final optim result.where the interpolation function   is defined as: where q is the penalty factor, and its influence on the interpolation function is sh Figure 2a.The smaller the q value is, the greater the degree of aggregation of desig ables to both ends of 0/1 by the interpolation function   , which is conducive to av the generation of intermediate density.And when  ∈ (0, 0.1) and q = 0.01, its in on the interpolation function is most obvious, that is, it reduces the function value to 0 in the smallest  change.Therefore, in this optimization process, q value is 0.
(a) However, the use of the above-mentioned interpolation function still inevitabl to an undesired grayscale phenomenon in the final optimization result.Therefore, study, the generation of grayscale is further suppressed by further difference of th polation function   .The interpolation function of   is expressed as: However, the use of the above-mentioned interpolation function still inevitably leads to an undesired grayscale phenomenon in the final optimization result.Therefore, in this study, the generation of grayscale is further suppressed by further difference of the interpolation function γ p .The interpolation function of γ p is expressed as: where f γ p is the interpolation function of the newly added γ p .In traditional research, f γ p = γ p is taken.The comparison is shown in Figure 2b.When γ p is close to the intermediate value, the newly added interpolation function f γ p is significantly larger than the traditional method, which in turn increases the Darcy friction of the intermediate material and inhibits the appearance of the intermediate variable value.Therefore, this improved method can further reduce the gray level of the optimization domain, so as to achieve a better optimization effect.α max and α min are the minimum and maximum values of α, respectively, and represent the resistance coefficient of the liquid and solid phases, respectively.The resistance coefficient of the liquid phase is very small, which is set to 0 in this paper.The value of the resistance coefficient of the solid phase α max is determined by the Darcy number Da and the characteristic length L, which is written as: In topology optimization of the fluid-structure coupled heat transfer problem using the variable density method, we redefine material properties and heat source distribution to ensure the difference in material properties of the solid domain, the fluid domain, and a single heat source that generates heat only in solid domain.Based on this, the interpolation function relationship between material properties and design variables γ p , and the heat source are defined as: where k is the thermal conductivity; c p is the heat capacity at constant pressure.The subscripts f and s represent fluid materials and solid materials, respectively.In addition, the thermophysical parameters of the fluid and solid materials are shown in Table 1.B is the heat generation coefficient which controls the speed of the heat generation [5] and a is a constant, and the value is 167.T ref is the reference temperature, which is 298 K. Combining Equations ( 5), (7), and ( 12), it can be seen that when the design variable γ = 1 (representing the fluid domain), the amount of heat production per unit volume Q is 0, that is, no heat is produced in the fluid flow region.

Filter and Projection
In the optimization process, in order to ensure the smoothness of the design variable distribution, it needs to be filtered by the Partial Differential Equation (PDE) filter.PDE filter.The design variables are filtered using the Helmholtz partial differential equation [29,30], which is defined as: where R is the filter radius, which is selected as the minimum grid size in this paper; γ c is the design variable before filtering; and γ f is the design variable after filtering.
Helmholtz partial differential equation filtering design variables can solve problems such as grid dependence and smoothness of the distribution of design variables.However, it also increases the overall grayscale of the design domain, resulting in a large amount of gray area (intermediate density) between the solid and fluid materials in the design domain.In order to reduce the grayscale of the design area, the hyperbolic tangent projection method is used in the filtered design area, which is expressed as: where γ is the design variable after projection; γ β is the projection point, and its value is 0.5.β γ is the projection slope, and its influence on the projection is shown in Figure 2c.The larger the value of β γ , the steeper the curve is, and the better the polarization distribution effect of design variables; however, it will also lead to instability of the calculated value and poor convergence.In order to improve the stability and convergence of the optimization problem, and avoid the optimization problem from falling into a local optimal solution, we adopt the continuous approach, that is, take a small value of β γ in the initial stage of the optimization iteration, gradually increase the value of β γ , and finally make the fluid-solid boundary gradually clear [31].

Multi-Objective Optimization
The optimization of the fluid-structure interaction heat transfer needs to consider the problem of fluid target, thermal target, and solid distribution, at the same time.The problem of solid distribution is solved by adding a body force containing the design variable γ to the Navier-Stokes equations.For the fluid target and thermal target, weigh the importance of the two and adjust the proportion of the two targets to meet the actual engineering needs.Therefore, it is necessary to establish a multi-objective optimization equation first, then calculate a set of Pareto optimal solutions, and finally select the optimal solution according to the actual needs.
The larger the heat transfer in the design domain, the better the heat transfer effect.In addition, the smaller the flow energy consumption, the smaller the flow loss.Therefore, the first optimization objective considered in this study is to maximize the heat exchange in the entire design area, and the second optimization objective is to minimize the flow energy dissipation in the whole design area.The heat exchange and flow energy consumption are defined as: In the study, the simple additive weighting (SAW), that is, the linear weighted summation method, is used to combine the above two optimization objectives.In order to avoid an invalid weighting coefficient, due to the difference of the calculation results of the two functions, J 1 and J 2 need to be normalized.The normalization processing method is used to perform scaling processing on the above single-objective function, and the scaling processing method is as follows: where I takes the value of 1 and 2. J max i and J min i represent the maximum and minimum values of the objective function, respectively.Therefore, the multi-objective function weighted by the objective functions of heat exchange J 1 and flow energy consumption J 2 are expressed as follows: where J is the multi-objective function, ω 1 and ω 2 are the weighting coefficients of heat exchange and flow energy consumption, respectively.The expression of the multi-objective function can be controlled by adjusting the size of the weighting coefficient, and then different optimization structures can be obtained.However, if the values of ω 1 and ω 2 are not constrained, many repeated cases will be obtained, which is meaningless and is not conducive to the subsequent analysis of the results.For example, the working condition of ω 1 = 0.1, ω 2 = 0.2 and the working condition of ω 1 = 1, ω 2 = 2 are repeated.Therefore, we set 0 ≤ ω 1 ≤ 1, 0 ≤ ω 2 ≤ 1, and ω 1 = 1 − ω 2 , and when ω 1 > 0.5, the design domain takes the maximum heat exchange as the main objective function, that is, the objective function of heat exchange dominates, otherwise the design domain takes the minimum flow energy dissipation as the main objective function.
Based on the above model analysis, the mathematical expression of the multi-objective topology optimization of the fluid-structure coupled heat transfer system is as follows:

The Solution Process of Topology Optimization
The calculation and optimization problems are identified by the interface between the commercial finite element software COMSOL Multiphysics 5. 6 [32] and MATLAB [33].The computational fluid dynamics module (CFD) and fluid heat transfer module (ht) are used to compute related problems such as fluid flow and heat transfer, respectively.The design domain is discretized using bilinear finite elements [34].Free triangular meshes are used to generate meshes by controlling the maximum mesh size, and the dependence of calculation results on meshes is eliminated through projection and filtering.The convergence criterion of the optimization process is defined as: where, n is the number of calculation iteration steps, ε = 1 × 10 −8 .If the calculation results of the two iteration steps in the optimization process conform to Equation ( 21), or the maximum number of iteration steps is reached, then the iteration is terminated.
The flow chart of the solution process for the topology optimization of the microchannel heat sink is shown in Figure 3.The optimization process includes the following main steps: (1) Give the initial value of the design variable γ in the design domain.

Micro-Channel Heat Sink Optimization under Different θ
In practical engineering applications, different inlet and outlet an can be selected.Therefore, it is necessary to study the influence of differe positions on the heat transfer, flow energy consumption, and average o in the design domain.Taking  1 ′ as the objective function, that is, select = 0 in Equation ( 18), the topology structure with the best heat transfer pe different bifurcation angles θ is studied.
The optimal topological structure and its velocity field, temperature dissipation density field distribution under different bifurcation angles ure 4. In the optimal topology, the white and black regions represent th regions, respectively, and their corresponding design variables are 1 an However, the design variables in the grey area are between 0 and 1, wh termediate material densities, which have no practical significance and nated as much as possible in the optimal topology.It can be seen from density of the intermediate material in this study has been well eliminate the interpolation function (  ) used in this paper has played a good r

Micro-Channel Heat Sink Optimization under Different θ
In practical engineering applications, different inlet and outlet angle distributions can be selected.Therefore, it is necessary to study the influence of different inlet and outlet positions on the heat transfer, flow energy consumption, and average outlet temperature in the design domain.Taking J 1 as the objective function, that is, selecting ω 1 = 1 and ω 2 = 0 in Equation ( 18), the topology structure with the best heat transfer performance under different bifurcation angles θ is studied.
The optimal topological structure and its velocity field, temperature field, and energy dissipation density field distribution under different bifurcation angles are shown in Figure 4.In the optimal topology, the white and black regions represent the fluid and solid regions, respectively, and their corresponding design variables are 1 and 0, respectively.However, the design variables in the grey area are between 0 and 1, which represent intermediate material densities, which have no practical significance and should be eliminated as much as possible in the optimal topology.It can be seen from Figure 4 that the density of the intermediate material in this study has been well eliminated, indicating that the interpolation function f γ p used in this paper has played a good role.ure 4. In the optimal topology, the white and black regions represent the fluid and solid regions, respectively, and their corresponding design variables are 1 and 0, respectively.However, the design variables in the grey area are between 0 and 1, which represent intermediate material densities, which have no practical significance and should be eliminated as much as possible in the optimal topology.It can be seen from Figure 4 that the density of the intermediate material in this study has been well eliminated, indicating that the interpolation function (  ) used in this paper has played a good role.Comparing the optimal topological structures under four different bifurcation angles, it can be found that the flow channel structures have certain similarities: The fluid flow channel has two to three main channels distributed from the inlet to the outlet.The fluid flows into the heat sink from the inlet, and is divided by the solid heat source in the design domain when it flows through it, and finally converges in the outlet channel and flows out after the heat exchange in the design domain.The solid heat sources in the optimized design domain show a scale-like distribution, which is different from the conventional heat sink, and these scale-like distributions of solid heat sources lead to the refinement of many tiny flow channels in the heat sink.However, the difference is that with the change of the bifurcation angle, the position of the solid heat source, that acts as a diverter at the entrance of the design domain, also changes and the smaller the bifurcation angle Comparing the optimal topological structures under four different bifurcation angles, it can be found that the flow channel structures have certain similarities: The fluid flow channel has two to three main channels distributed from the inlet to the outlet.The fluid flows into the heat sink from the inlet, and is divided by the solid heat source in the design domain when it flows through it, and finally converges in the outlet channel and flows out after the heat exchange in the design domain.The solid heat sources in the optimized design domain show a scale-like distribution, which is different from the conventional heat sink, and these scale-like distributions of solid heat sources lead to the refinement of many tiny flow channels in the heat sink.However, the difference is that with the change of the bifurcation angle, the position of the solid heat source, that acts as a diverter at the entrance of the design domain, also changes and the smaller the bifurcation angle is, the closer it is to the edge of the entrance.Correspondingly, the position of the solid heat source at the outlet of the design domain, which plays the role of collecting fluid, also changes similarly.This leads to a large curvature of the fluid flowing through the inlet and outlet channels of the heat sink under a small bifurcation angle, and its flow distance is small, which is not conducive to fluid flow and heat exchange.
Comparing the velocity field and temperature field at different bifurcation angles, it is obvious that the velocity distribution and temperature distribution are similar to the above-mentioned flow channel structure.For the heat sink with θ = 45 • , a part of the flow channel connecting the inlet and outlet is short and narrow, and the overall flow channel is tapered.As a result, when the fluid flows through this part of the fine flow channel, the degree of bending changes greatly, the flow velocity increases sharply, and the velocity gradient is large, which in turn leads to high fluid flow velocity, insufficient heat exchange, and uneven outlet temperature.In addition, it can be found from Figure 4 that for the heat sink with large bifurcation angles, the inlet and outlet channels are distributed farther, the flow channels in the design domain are fully extended, the distance of fluid flow is longer, the fluid heat exchange is more sufficient, and the outlet temperature is more uniform.It is these fully extended scaly-like distributions of solid heat sources and tiny flow channels that make the outlet velocity and temperature distribution of the heat sink under large bifurcation angles more uniform, and the average temperature of the entire design domain lower, which also shows that the heat transfer performance of a micro-channel heat sink with the large bifurcation angle is better.
The energy dissipation density distribution under different bifurcation angles is shown in Figure 4d.By comparison, it can be found that the areas with high energy dissipation density are all distributed in the solid heat source side impacted by the fluid, which is also the area with large flow loss in the microchannel heat sink.On the side of a solid heat source that is not impacted by fluid, the energy dissipation density is small, and the corresponding flow loss is small.In particular, when the bifurcation angle of the microchannel heat sink is 45 • and 90 • , the fluid curvature flowing through the inlet and outlet channels of the heat sink is large, and the energy dissipation density is also greatly increased, which can be clearly shown in Figure 4d.However, the energy dissipation density distribution of the microchannel heat sink with bifurcation angles of 135 • and 180 • is relatively uniform.
Figure 5 shows the variation of the heat exchange J 1 , the flow energy consumption J 2 , the average outlet temperature T, and the thermal resistance R, with the bifurcation angle θ in the design domain.The thermal resistance R is defined as follows: where T denotes the average outlet temperature, T in denotes the inlet temperature, which is set at 273.15 K, and J 1 denotes the heat exchange of the design domain.
With the increase of the bifurcation angle θ, the J 1 , J 2 , and T of the optimized domains all show an increasing trend, and R tends to decrease first and then increase.In particular, when θ = 135 • , the heat exchange J 1 and the average outlet temperature T reach their peaks.The flow energy consumption J 2 does not increase, and the thermal resistance R is relatively small when θ > 135 • .Compared with the small bifurcation angle, the increase of J 1 and T are larger, and the increase of energy dissipation J 2 is relatively small, and compared with the large bifurcation angle, J 1 and T are larger, and J 2 is smaller.In particular, although the heat resistance R of the micro-channel heat sink at θ = 90 • is the smallest, the heat exchange and the average outlet temperature at this time are smaller.In addition, compared with the optimal topology of the heat sink with θ = 135 • and θ = 180 • in Figure 4, it can be seen that when the fluid flows through the heat sink with θ = 135 • , more fluid flows through the upper and middle flow channels, and after heat exchange with the solid heat sources on both sides of the flow channel, it gathers in other flow channels and flows out.Correspondingly, when the fluid flows through the heat sink with θ = 180 • , the fluid is divided into the upper and lower main branch channels and flows out after heat exchange with the solid heat source at the side of the central channel.However, the scale-like distribution of solid heat sources of the heat sink with θ = 135 • is more compact than that of the heat sink with θ = 180 • , and the heat transfer between fluid and heat source is more sufficient.As shown in Figure 5, the average temperature at the outlet is higher, and the heat transfer in the design domain is greater.Therefore, in this study, the bifurcation angle of 135 • with better comprehensive performance will be selected for subsequent multi-objective optimization.
The energy dissipation density distribution under different bifurcation angles is shown in Figure 4d.By comparison, it can be found that the areas with high energy dissipation density are all distributed in the solid heat source side impacted by the fluid, which is also the area with large flow loss in the microchannel heat sink.On the side of a solid heat source that is not impacted by fluid, the energy dissipation density is small, and the corresponding flow loss is small.In particular, when the bifurcation angle of the microchannel heat sink is 45° and 90°, the fluid curvature flowing through the inlet and outlet channels of the heat sink is large, and the energy dissipation density is also greatly increased, which can be clearly shown in Figure 4d.However, the energy dissipation density distribution of the microchannel heat sink with bifurcation angles of 135° and 180° is relatively uniform.
Figure 5 shows the variation of the heat exchange J1, the flow energy consumption J2, the average outlet temperature T, and the thermal resistance R, with the bifurcation angle θ in the design domain.The thermal resistance R is defined as follows: where T denotes the average outlet temperature, Tin denotes the inlet temperature, which is set at 273.15 K, and J1 denotes the heat exchange of the design domain.With the increase of the bifurcation angle θ, the J1, J2, and T of the optimized domains all show an increasing trend, and R tends to decrease first and then increase.In particular, when θ = 135°, the heat exchange J1 and the average outlet temperature T reach their peaks.The flow energy consumption J2 does not increase, and the thermal resistance R is relatively small when θ > 135°.Compared with the small bifurcation angle, the increase of J1 and T are larger, and the increase of energy dissipation J2 is relatively small, and compared with the large bifurcation angle, J1 and T are larger, and J2 is smaller.In particular, although the heat resistance R of the micro-channel heat sink at θ = 90° is the smallest, the heat exchange and the average outlet temperature at this time are smaller.In addition, compared with the optimal topology of the heat sink with θ = 135° and θ = 180° in Figure

Micro-Channel Heat Sink Optimization under Different ω 1
In order to study the performance characteristics under different weighting coefficients ω 1 , micro-channel heat sinks under different ω 1 were selected for analysis.Figure 6 shows the optimal topology, velocity field, temperature field and energy dissipation density field distribution under different ω 1 .Comparing the optimal topological structures under different ω 1 , it can be found that when ω 1 = 0, the flow channel of the micro-channel is the main channel connecting the inlet and outlet, and no branches appear.After the fluid flows into the design domain from the inlet, it flows directly through the main channel and flows out from the outlet.This is because the objective function J 1 of heat exchange under this working condition does not contribute to the multi-objective function J, and the optimization objective degenerates into a single-objective function that minimizes the flow energy dissipation J 2 .With the increase of ω 1 , the proportions of J 1 and J 2 in the multiobjective function J increase and decrease, respectively, which means that the objective function of heat exchange gradually dominates.At this time, the main channel connecting the inlet and outlet of the micro-channel is gradually shunted by the solid heat source, and the area of the main channel is reduced.The branch channels appear preferentially below the main channel, and gradually expand to the entire optimization domain, and the branch channels are gradually further refined by the scaly-like distribution of solid heat sources.In particular, when ω 1 = 1, the flow energy dissipation J 2 does not contribute to the multi-objective function J, and the optimization objective degenerates into a single-objective function that minimizes heat exchange J 1 .The fluid at the entrance of the optimization domain is completely diverted by the solid heat source, and the main channel no longer exists, but is refined into several tributaries by a large number of scaly-like distribution of solid heat sources.Compared with the working condition of ω 1 < 1, the distribution of the solid heat source is more refined, and the small tributaries are more diverse.We comprehensively analyzed the optimal topology and energy dissipation density, and found that the area with a large energy dissipation density was also distributed on the solid heat source side impacted by the fluid.The smaller the weighting coefficients ω 1 is, the smaller the energy dissipation density of the whole heat sink domain is.This change trend also matches the flow energy dissipation in the optimization goal.
Energies 2022, 15, x FOR PEER REVIEW 13 of 19 flow energy dissipation J2.With the increase of ω1, the proportions of J1 and J2 in the multiobjective function J increase and decrease, respectively, which means that the objective function of heat exchange gradually dominates.At this time, the main channel connecting the inlet and outlet of the micro-channel is gradually shunted by the solid heat source, and the area of the main channel is reduced.The branch channels appear preferentially below the main channel, and gradually expand to the entire optimization domain, and the branch channels are gradually further refined by the scaly-like distribution of solid heat sources.In particular, when ω1 = 1, the flow energy dissipation J2 does not contribute to the multi-objective function J, and the optimization objective degenerates into a singleobjective function that minimizes heat exchange J1.The fluid at the entrance of the optimization domain is completely diverted by the solid heat source, and the main channel no longer exists, but is refined into several tributaries by a large number of scaly-like distribution of solid heat sources.Compared with the working condition of ω1 < 1, the distribution of the solid heat source is more refined, and the small tributaries are more diverse.We comprehensively analyzed the optimal topology and energy dissipation density, and found that the area with a large energy dissipation density was also distributed on the solid heat source side impacted by the fluid.The smaller the weighting coefficients ω1 is, the smaller the energy dissipation density of the whole heat sink domain is.This change trend also matches the flow energy dissipation in the optimization goal.In addition to the above optimal topological structures under different ω1, we also compared and analyzed the corresponding velocity field and temperature field, which was similar to the above flow channel distribution.As there is no branching effect, the flow is the smoothest, and the flow energy dissipation is also the smallest; the contact between the fluid and the solid wall is only the wall part of the main channel, so the fluidsolid contact area is the smallest and the heat transfer effect is the worst, which corresponds to the objective function after degeneration.With the increase of ω1, in addition to flowing through the main channel, the fluid flow in the design domain will also be divided by the branch flow channels.The larger the ω1, the smaller the area occupied by the main channel, the more scaly-like distribution of solid heat sources, the more branch flow channels, and the more obvious the diversion effect.The above features also mean that the fluid-solid contact area is increased, thus the fluid flow energy dissipation is increased, and the corresponding heat transfer effect is also improved.When ω1 = 1, the fluid-solid contact area reaches the maximum value, and the fluid flow becomes the most complicated.The fluid flows out from the three branch flow channels and several small channels, and the energy dissipation increases greatly.At the same time, the fluid takes away the heat of the optimization domain to the greatest extent, and the heat transfer effect is greatly increased.We can also find from Figure 6 that when ω1 is large, the temperature distribution of the entire design domain is more even, with no large temperature gradient, which is more reasonable.
Figure 7 shows the optimal solutions, Pareto frontier, and optimal configuration under different weighting factors, where the single-objective functions J1 and J2 are the abscissa and ordinate, respectively.The Pareto frontier of the micro-channel heat sinks in this paper all have the characteristics of an upward convex curve, which is consistent with the convex function shape described by Timothy et al. [35].Its topological structure change has obvious regularity with the increase of ω1: the branch channels are preferentially split under the main runner, and the main channel is slowly replaced by the upper, middle and lower branch channels.The distribution of the solid heat source tends to be uniform and presents a scaly-like distribution, the aspect ratio of a single solid heat source gradually decreases, and the shape of the solid heat source evolves from a large strip to a small square (approximately).Using the linear logarithmic transform function (Log3P1) to perform nonlinear curve fitting on the Pareto frontier, the fitting equation is obtained using equation (23).Based on this equation, engineering designers can select the corresponding performance parameters according to their actual engineering needs, and the optimal design has been achieved.

(
) 9466.9 1005.5 log 4.377 10 Figure 8 shows the heat exchange J1, flow energy dissipation J2, the average outlet temperature T and the thermal resistance R under different weighting coefficients ω1.According to the comparative analysis of Figure 7, it can be concluded that for the optimal In addition to the above optimal topological structures under different ω 1 , we also compared and analyzed the corresponding velocity field and temperature field, which was similar to the above flow channel distribution.As there is no branching effect, the flow is the smoothest, and the flow energy dissipation is also the smallest; the contact between the fluid and the solid wall is only the wall part of the main channel, so the fluid-solid contact area is the smallest and the heat transfer effect is the worst, which corresponds to the objective function after degeneration.With the increase of ω 1 , in addition to flowing through the main channel, the fluid flow in the design domain will also be divided by the branch flow channels.The larger the ω 1 , the smaller the area occupied by the main channel, the more scaly-like distribution of solid heat sources, the more branch flow channels, and the more obvious the diversion effect.The above features also mean that the fluid-solid contact area is increased, thus the fluid flow energy dissipation is increased, and the corresponding heat transfer effect is also improved.When ω 1 = 1, the fluid-solid contact area reaches the maximum value, and the fluid flow becomes the most complicated.The fluid flows out from the three branch flow channels and several small channels, and the energy dissipation increases greatly.At the same time, the fluid takes away the heat of the optimization domain to the greatest extent, and the heat transfer effect is greatly increased.We can also find from Figure 6 that when ω 1 is large, the temperature distribution of the entire design domain is more even, with no large temperature gradient, which is more reasonable.
Figure 7 shows the optimal solutions, Pareto frontier, and optimal configuration under different weighting factors, where the single-objective functions J 1 and J 2 are the abscissa and ordinate, respectively.The Pareto frontier of the micro-channel heat sinks in this paper all have the characteristics of an upward convex curve, which is consistent with the convex function shape described by Timothy et al. [35].Its topological structure change has obvious regularity with the increase of ω 1 : the branch channels are preferentially split under the main runner, and the main channel is slowly replaced by the upper, middle and lower branch channels.The distribution of the solid heat source tends to be uniform and presents a scaly-like distribution, the aspect ratio of a single solid heat source gradually decreases, and the shape of the solid heat source evolves from a large strip to a small square (approximately).Using the linear logarithmic transform function (Log3P1) to perform nonlinear curve fitting on the Pareto frontier, the fitting equation is obtained using Equation (23).Based on this equation, engineering designers can select the corresponding performance parameters according to their actual engineering needs, and the optimal design has been achieved.J 1 = 9466.9+ 1005.5 × log J 2 − 4.377 × 10 −5 (23) optimal topology with ω1 = 0.8, the fluid at the inlet is completely shunted by the solid heat source, and the main flow channel is reduced to no longer exist.Therefore, the increase of J1 is larger than that of ω1 < 0.8, and the J2 is also greatly increased.When ω1 = 0.8, the average outlet temperature reaches the maximum value, the flow energy consumption is smaller than that of ω1 > 0.8, and the change of heat transfer is not large.Therefore, when choosing the best topology design, we can pay special attention to the case where the main channel is reduced until it no longer exists.Figure 8 shows the heat exchange J 1 , flow energy dissipation J 2 , the average outlet temperature T and the thermal resistance R under different weighting coefficients ω 1 .According to the comparative analysis of Figure 7, it can be concluded that for the optimal topology with ω 1 < 0.3, the fluid directly flows through the main channel connecting the inlet and outlet of the design domain, with no branches and small channels, so J 1 , J 2 , T, and R do not change much.For the optimal topology with ω 1 = 0.3, there are branch flow channels in the lower part of the main channel; however, most of the fluid still flows out through the main channel, and J 1 and J 2 show little change.Compared to the working condition of ω 1 = 0.2, J 1 increases by 670 W/m, J 2 increases by 0.11 Mw/m, T increases by 3 K, and R decreases by 0.58 × 10 −3 K•m/W.For the optimal topology with ω 1 > 0.3, the branch flow channels near the main runner are gradually refined and increased, so J 1 , J 2 , and T also gradually increase, and correspondingly, R gradually decreases.In particular, for the optimal topology with ω 1 = 0.8, the fluid at the inlet is completely shunted by the solid heat source, and the main flow channel is reduced to no longer exist.Therefore, the increase of J 1 is larger than that of ω 1 < 0.8, and the J 2 is also greatly increased.When ω 1 = 0.8, the average outlet temperature reaches the maximum value, the flow energy consumption is smaller than that of ω 1 > 0.8, and the change of heat transfer is not large.Therefore, when choosing the best topology design, we can pay special attention to the case where the main channel is reduced until it no longer exists.For the optimal topology with ω1 > 0.3, the branch flow channels near the main runner are gradually refined and increased, so J1, J2, and T also gradually increase, and correspondingly, R gradually decreases.In particular, for the optimal topology with ω1 = 0.8, the fluid at the inlet is completely shunted by the solid heat source, and the main flow channel is reduced to no longer exist.Therefore, the increase of J1 is larger than that of ω1 < 0.8, and the J2 is also greatly increased.When ω1 = 0.8, the average outlet temperature reaches the maximum value, the flow energy consumption is smaller than that of ω1 > 0.8, and the change of heat transfer is not large.Therefore, when choosing the best topology design, we can pay special attention to the case where the main channel is reduced until it no longer exists.

Mesh Independence Test
We take ω1 = 1 and θ = 135° as an example to verify the mesh independence.The optimal structure under three different mesh sizes is shown in the Figure 9, and the mesh size selected in this paper is the medium mesh.The mesh number, heat exchange J1, flow energy consumption J2, average outlet temperature T, and its relative former difference of the micro-channel heat sink are shown in Table 2.In addition, the cell dimensions of the coarse mesh, medium mesh, and fine mesh are L/50, L/60, and L/70, respectively.The mesh independence test shows that the optimal topology under different mesh sizes is consistent on the whole, while there are some differences in the distribution of scale-like solid heat sources in some regions.The relative former difference of parameters under different mesh sizes is very small, which also shows that the final optimization design result is independent of the number of meshes.Therefore, the medium mesh size selected in this paper is reasonable and reliable.

Mesh Independence Test
We take ω 1 = 1 and θ = 135 • as an example to verify the mesh independence.The optimal structure under three different mesh sizes is shown in the Figure 9, and the mesh size selected in this paper is the medium mesh.The mesh number, heat exchange J 1 , flow energy consumption J 2 , average outlet temperature T, and its relative former difference of the micro-channel heat sink are shown in Table 2.In addition, the cell dimensions of the coarse mesh, medium mesh, and fine mesh are L/50, L/60, and L/70, respectively.The mesh independence test shows that the optimal topology under different mesh sizes is consistent on the whole, while there are some differences in the distribution of scale-like solid heat sources in some regions.The relative former difference of parameters under different mesh sizes is very small, which also shows that the final optimization design result is independent of the number of meshes.Therefore, the medium mesh size selected in this paper is reasonable and reliable.

Mesh Independence Test
We take ω1 = 1 and θ = 135° as an example to verify the mesh independence.The optimal structure under three different mesh sizes is shown in the Figure 9, and the mesh size selected in this paper is the medium mesh.The mesh number, heat exchange J1, flow energy consumption J2, average outlet temperature T, and its relative former difference of the micro-channel heat sink are shown in Table 2.In addition, the cell dimensions of the coarse mesh, medium mesh, and fine mesh are L/50, L/60, and L/70, respectively.The mesh independence test shows that the optimal topology under different mesh sizes is consistent on the whole, while there are some differences in the distribution of scale-like solid heat sources in some regions.The relative former difference of parameters under different mesh sizes is very small, which also shows that the final optimization design result is independent of the number of meshes.Therefore, the medium mesh size selected in this paper is reasonable and reliable.

Conclusions
The method of topology optimization design is used to study the micro-channel heat sink by considering the coupling of fluid-solid and heat transfer.Based on the two objective functions of achieving high heat transfer capacity and low flow energy dissipation, a multiobjective function is proposed to optimize the structure topology of the heat sink, and the influences of the inlet and outlet bifurcation angles and heat transfer weighting coefficients on the optimal topology of the radiator are revealed.The main conclusions are as follows.

1.
In the numerical implementation of topology optimization design, the improved interpolation function, PDE filter, and hyperbolic tangent projection effectively eliminate the gray area of the topology structure and tiny solid particles, which improve the smoothness of the flow channels and avoid checkerboard results.

2.
The optimal topology of the micro-channel heat sink has obvious regularity with the change of θ.With the increase of θ, the overall performance parameters show a monotonically increasing trend.Among them, the heat exchange and the average outlet temperature reach the maximum value when θ = 135 • , and the overall flow energy dissipation is relatively small.Therefore, in the selection of the heat sink, different bifurcation angles can be considered based on the calculation results, rather than the intuitive or empirical selection of θ = 180 • .

3.
Based on the multi-objective topology optimization design, micro-channel heat sinks with different heat transfer weighting coefficients are obtained.When minimizing energy dissipation dominates, the heat sink is occupied by the main channel connecting the inlet and outlet, and the fluid mainly flows out through the main flow channel with few branch flow channels.Its flow energy dissipation is smaller and the heat exchange effect is worse.When the maximum heat exchange is dominant, the main channel of the heat sink is filled with scaly-like distribution of solid heat sources, and the main channel gradually shrinks and is replaced by the upper, middle, and lower branch channels.In addition, the branch flow channels are further refined into many tiny flow channels, which makes the heat exchange effect of the heat sink better and the flow energy consumption larger.

Figure 1 .
Figure 1.Design model, dimensions, and boundary conditions of the micro-channel heat sink.
, and p represent the flow velocity, density and pressure of the fluid, respec tively, μ represents the dynamic viscosity of the fluid, and F represents the body force o the fluid flow.

Figure 1 .
Figure 1.Design model, dimensions, and boundary conditions of the micro-channel heat sink.

Figure 2 .
Figure 2. Interpolation function and hyperbolic tangent projection.(a) Effect of penalty factor on interpolation function; (b) comparison of two interpolation functions; and (c) effect of projection slope on hyperbolic tangent function.

( 2 )
Solve equations at various current parameters.(3) Calculate the objective function and constraints.(4) Calculate the design sensitivity under the current concomitant variables.(5) Update the design variable γ using the Global Convergence Method of Moving Asymptotes (GCMMA).(6) Loop the above (2)-(5) steps until the iteration termination condition is met.

Figure 3 .
Figure 3.The flow chart of the optimization design.

Figure 3 .
Figure 3.The flow chart of the optimization design.

Figure 4 .
Figure 4. Variation of optimal topology, velocity field, temperature field, and energy dissipation density field distribution with θ.(a) Optimal topology; (b) Velocity field; (c) Temperature field; and (d) Energy dissipation density field.

Figure 4 .
Figure 4. Variation of optimal topology, velocity field, temperature field, and energy dissipation density field distribution with θ.(a) Optimal topology; (b) Velocity field; (c) Temperature field; and (d) Energy dissipation density field.

Figure 5 .
Figure 5. Variation of performance parameters under different θ values.(a) The heat exchange J1 and the flow energy consumption J2, and (b) The average outlet temperature T and the thermal resistance R.

Figure 5 .
Figure 5. Variation of performance parameters under different θ values.(a) The heat exchange J 1 and the flow energy consumption J 2, and (b) The average outlet temperature T and the thermal resistance R.

Figure 6 .
Figure 6.Variation of optimal topology, velocity field, temperature field and energy dissipation density field distribution with ω 1 .(a) Optimal topology; (b) Velocity field; (c) Temperature field; and (d) Energy dissipation density field.

Energies 2022 ,
15, x FOR PEER REVIEW 15 of 19 topology with ω1 < 0.3, the fluid directly flows through the main channel connecting the inlet and outlet of the design domain, with no branches and small channels, so J1, J2, T, and R do not change much.For the optimal topology with ω1 = 0.3, there are branch flow channels in the lower part of the main channel; however, most of the fluid still flows out through the main channel, and J1 and J2 show little change.Compared to the working condition of ω1 = 0.2, J1 increases by 670 W/m, J2 increases by 0.11 mW/m, T increases by 3 K, and R decreases by 0.58 × 10 −3 K•m/W.

Figure 8 .
Figure 8. Variation of performance parameters under different ω 1 values.(a) The heat exchange J 1 and the flow energy consumption J 2 and (b) The average outlet temperature T and the thermal resistance R.

Figure 8 .
Figure 8. Variation of performance parameters under different ω1 values.(a) The heat exchange J1 and the flow energy consumption J2 and (b) The average outlet temperature T and the thermal resistance R.

Table 1 .
Physical properties of materials.

Table 2 .
Mesh details for the mesh independence test.

Table 2 .
Mesh details for the mesh independence test.

Table 2 .
Mesh details for the mesh independence test.