Topology optimization of a pseudo 3D microchannel heat sink: influence of the initial guess

Nowadays, the optimization of microchannel heat sinks involves using algorithms, e.g., topology optimization (TO), to find the optimal layout of fluid channels to enhance certain objectives, e.g., heat transfer, pressure drop. The outcome is an unintuitive material distribution in the design space, subject to constraints such as manufacturability, fluid flow requirements, and thermal performance. This approach has been shown to provide significant benefits compared to traditional design methods, leading to increased cooling performance. However, depending on the application, defining an initial (i.e., starting) guess for the optimization problem is not trivial. Initial guess for TO, i.e., starting design for the optimization process, can play a crucial role in determining the final optimized design. Its choice may affect the convergence rate, the performance of the final solution, and the computational effort, as it can be set in a variety of ways, e.g., a random distribution, a preconceived design, or a combination of these two. In this study, we propose a heuristic approach for incorporating a genetic algorithm (GA) as outer algorithm into the TO routine in order to find the initial guess that ensures the lowest thermal resistance while limiting the pumping power. The aim is to correlate initial designs with final objectives and to define new useful insights on microchannel optimization, since TO has proved to be a promising area of research for the development of advanced thermal management systems.


Introduction
Topology optimization (TO) is a current mathematical technique used in several engineering branches to define the optimal material layout within a predefined space [1].While TO has been widely used in structural design, its application in CFD problems is relatively new and presents intricate challenges, involving both purely numerical and conceptual aspects [2].When coupling TO to CFD problems, the objective is to optimize shape and topology of the solid/fluid domain to achieve desired flow characteristics and fulfil thermal targets.In this frame, the choice of design variables, the selection of appropriate constraints and boundary conditions is a well-known challenge, as well as the complex physics of fluid flow, e.g., turbulence, nonlinearity, and multiphase flow [3].
However, among these issues, one of the critical input parameters that can significantly impact the optimization process and the final design is the initial guess [4].The initial guess is an initial estimate of the optimized design that serves as a starting point for the TO algorithm.The choice of the initial guess can affect the convergence rate, as well as the topology and performance of the optimized design.In this case the problem can be defined as initial guess dependent (IGD).This means that different initial designs result in diverse topologies, although having close objective values [5].Since the initial guess can affect the convergence rate and the quality of the optimized design, a reasonable option should be a starting guess close to the optimal solution, with features of the expected design [2].Thus, the initial guess should be chosen based on prior knowledge of the design problem, which is not always feasible, especially in the area of thermofluid-dynamics where designs tend to have aerodynamic or fluid-based shapes.
Thermofluid-dynamics topology optimization problems are quite non-convex and can easily converge to only locally optimal topologies depending on the starting guess or initial design [5], [6].Thus, exploring different initial design options is mandatory to comprehensively explore the design domain and avoid trivial solutions.More specifically, regarding thermal management, microchannel heat sink optimization has proven to be IGD [4], [6].An apriorist evaluation of various configurations with different fin placements can be time-consuming, which added to the computational cost of TO makes the process lengthy and tedious.To this end, the aim of this work is to apply a heuristic approach for integrating a genetic algorithm (GA) as an outer algorithm into the TO routine: GA generates multiple pin fins' configurations that are used as initial guesses by TO to find a final design able to minimize the thermal resistance between the channel layer and the base plate layer, and the pressure drop.By doing this, starting guesses can be correlated with final objectives, defining new helpful insights on microchannel optimization.

Materials and methods
The model is based on a pseudo-3D formulation [6].A full 3D optimization is computationally expensive since topology optimization (TO) typically requires several hundred iterations before convergence to a final design.Accordingly, the pseudo 3D heat sink model approximates the original 3D one with two 2D thermally coupled problems, which are not independent themselves.Figure 1 shows a 3D sketch of the forced convection heat sink here considered.The base plate -the one concerned by the heat source -is thermally coupled to the channel layer thanks to an additional term.At the same time, this term is balanced on the channel layer side through a complementary one.Consequently, this approach -even being 2D -considers out-of-plane effects.The design domain, i.e., channel layer, is also depicted in figure 1.It consists of a rectangular domain, with an internal rectangular subdomain (Ωd) where topology takes place.The outer ones (Ωf) serve as inflow and outflow sections.Therefore, Ωd represents the domain where the initial guess has to be set.The starting design is chosen to be a microchannel heat sink with round pins.The arrangement and the number of round pins are not defined as constraints, whereas the initial volume constraint, i.e., the ratio between the area of the solid material and the total of Ωd, is fixed to a value of 0.2.Note that this constraint is only set on the initial design and not on the final one.Therefore, the aim of this work is to investigate the effect of varying pins position and size as initial guess for the TO problem.

Governing equations
To define the model equations, a continuous design field, γ(x) is defined in Ωd.Its values range from 0 to 1, and the upper value represents the fluid, whereas the lower one represents heat sink fin material.The reason why γ(x) takes one value over the other is specified in Section 2.2.
An incompressible, laminar, and steady-state flow is assumed throughout this case.The fluid and heat transfer problems in the thermofluid-dynamic design layer are modeled two-dimensionally, as described in the previous section.The 2D assumption is motivated by the fin height, significantly larger than the baseplate one.In addition, since the base plate's extension is much larger than its height, the original 3D thermal diffusion problem in the base plate is simplified to a 2D problem.Under these assumptions [6], the governing equations for the fluid are written as: where ρf is the fluid density, u is the fluid velocity vector, p is the pressure field, μf is the dynamic viscosity of the fluid, and F is the Brinkman friction term, which is used in TO flow problems to penalize flow through solid areas within the design domain. is the maximum inverse permeability of the porous medium and Iα(γ) a suitable function for the inverse permeability interpolation.Da is the Darcy number and Lc is a characteristic length which corresponds to the design domain width.It is worth noting that there is no Brinkman friction on the fluid outside of the design domain.
The 2D thermal convection-diffusion equations with/without a heat source are solved in the thermofluid design layer inside/outside the design domain: where Tf is the temperature field in the thermofluid design layer, cf is the specific fluid heat capacity, kf the thermal conductivity of the fluid, Ik(γ) is a function that interpolates between thermal conductivities.Δzchannel is the height of the air channel and fins, and qinter(γ) is the heat transferred from the solid base plate to the thermofluid design layer.
A 2D heat conduction problem is solved in the metal base plate: where ks is the thermal conductivity of the base plate, Ts is the temperature field in the base plate, Qsource is the prescribed heat production rate, Vbp is the base plate volume, and Δzbp is the base plate height.The base plate is assumed to produce heat at a constant rate in this study.The conductive heat transfer from the base plate into the heat sink fins, as well as the convective heat transfer from the base plate to the fluid, must be represented flexibly in the out-of-plane heat transfer between base plate and thermo-fluid design layer, qinter(γ).This is accomplished by interpolating between a high conductive heat transfer into the fins and a lower convective heat transfer to the fluid using a heat transfer coefficient: where hf is the convective heat transfer coefficient to the fluid, and Ih(γ) is a function that interpolates between hf and the heat transfer coefficient from the base plate into the heat sink fins, hs.The conductive heat transfer resistance in the z-direction in the fins is represented by the parameter hs, which is empirically determined.All walls are adiabatic, except for the inlet where the temperature Tin is fixed.The thermophysical properties of air that are used in this study are: cair = 1006 J/(kg K), kair = 0.024 W/(m K), air = 1.94 •10 -5 Pa s, and air = 1.204 kg/m 3 .
Pressure drop between inlet and outlet, Δp, is computed as difference of pressure at the inlet and outlet and a normal laminar inflow is set at the inlet boundary Γin.A symmetry boundary condition for the fluid flow is applied on bottom and upper boundaries.

Topology optimization
Density based TO is here applied [1].This method supports regularization via Helmholtz filter.The filter size is defined as two times the mesh size, which is set as 0.02 mm.Interpolation functions are used to represent intermediate values in the given design problem.An interpolation function is used in this work for the Brinkman friction term (Eq.3): where bα is a parameter that determines the interpolation's convexity.A rational approximation of material properties (RAMP)-style interpolation is used to interpolate the thermal conductivity within the design layer and the heat transfer between the heat sink base plate and the thermo-fluid design layer: where bj is the interpolation convexity parameter and Cj is defined in the respective case as kf/ks and hf/hs.
As described in equation (10), ϕ is the composition of 2 terms, ϕRth and ϕΔp, which are weighted through the terms wRth and wΔp.Minimizing both the objectives defines a multi-objective problem.To reduce the thermal resistance, the heat transfer area has to increase, while to minimize the pressure drop less flow resistance is required.In this work, the weighting factors have been considered equals.
It is worth specifying that in order to consider pressure drop as part of the objective, only inlet velocity, i.e., 0.1 m/s, has been specified as fluid boundary condition.Further, no volume constraint has been defined, apart from the one related to the initial guess.

Genetic algorithm
A genetic algorithm (GA) is run as outer algorithm into the TO routine to investigate the initial design condition that simultaneously fulfils thermal and fluid dynamic objectives.The design variables used in this work are the minimum circle radius (rmin) and the minimum distance between the circles (smin).The parameter rmin can vary from 0.06 mm to 0.13 mm, each 0.005 mm, whereas smin varies from 0.04 mm to 0.08 mm, each 0.01 mm.The reason for these ranges will be explained later.However, smin is set to be higher than the maximum mesh size.By doing this, an excessive mesh refinement between two circles -which should compromise the TO convergence -is avoided.All settings for GA tuning are taken from [7].

2.3.1.
Coupling GA with TO.The fitness function code generates a figure of randomly positioned circles within the rectangular domain.The code takes as input the dimensions of the rectangular domain, the minimum radius and distance between the circles, and a vector of target ratios between the total area of the circles and the area of the rectangle.For each target ratio, the code generates circles randomly until the total area of the circles reaches or exceeds the target ratio times the area of the rectangle.The code uses a while loop to generate circles one at a time.For each new circle, the code generates a random radius and position, and checks whether it intersects with any of the existing circles or exceeds the boundary limits of the rectangle.If the new circle is valid, it is added to the figure and the total area of the circles is updated.The loop continues until the total area of the circles reaches or exceeds the target ratio times the area of the rectangle, or the maximum number of iterations, i.e., 100, is reached.The code also sets up a figure window with a fixed size and saves the final image as a PNG file with a name that includes the target ratio.Note that the code assumes that the aspect ratio of the rectangular domain is 1.7:1 -equal to the design domain one -and that the circles are black with a solid fill.By doing this, the image can be imported in the COMSOL file and processed as initial guess.The method of moving asymptotes (MMA) is used to perform TO.

Results and discussions
This section reports the validation of the presented model under a Δp equal to 3 Pa and minimizing the thermal resistance Rth, and the results on the investigation of the initial guess by optimization using GA + TO.Each GA simulation lasts around 40 min, using a continuation scheme for the TO algorithm similar to the reference one [6].The number of overall iterations has been reduced from 1160 (reference) to 400, without losing accuracy (as proved in the following).

Model validation
Before running simulations about the optimization problem of equation ( 10), the TO reliability is tested replacing the optimization problem solved in [6].Thus, the model validation is assessed comparing qualitatively (figure 3) the final design, velocity magnitude field, temperature field of The main difference that can be spotted in the final design is the presence of material on the vertical boundaries, since in the present work there is no prescribed density condition on them.

Initial guess influence
170 Individuals computed through the GA + TO routine are depicted in figure 4.Among these, 5 relevant individuals have been highlighted: the optimal solution (ϕ = 61.95) and four worst solutions.In detail, in order to visualize a hypothetical dependency between the value of the objective function and the initial design, the plot of the initial design is shown for each of these individuals.By doing so, it can immediately be seen that the presence of a pin with a larger radius leads to a reduction in performance, while the presence of several pins with a smaller radius results in an optimal solution.To confirm what is presented through the plot of the objective function, the values of the design variables are grouped according to population and the frequency with which they occurred (figure 4).It shows that the GA preferred smaller radii and longer distances.Nevertheless, with regards to rmin, a non-negligible number of individuals can be noted, i.e., 26 designs with rmin = 0.11 mm, to which correspond to individuals with intermediate ϕ.At this point, it is useful to check the advantage of adopting a proper initial design (Optimal solution), as opposed to one of those "worst" highlighted above (Case 4).The comparison (figure 5) is made in terms of initial guess, final design, channel layer temperature field and velocity magnitude field.Starting the optimization with more pins with lower radii allows to reach a better final solution.The increase in the heat transfer allows to reduce the thermal resistance (table 1), and more aerodynamically shaped fins ensure a lower pressure drop than in Case 4, i.e., 12.46 Pa against 15.66 Pa.However, pressure drop values are consistently higher than the one used in the validation step, i.e., 3 Pa, which explains why the thermal resistance is significantly lower.Moreover, as listed in table 1, the maximum temperature at both channel and base plate layers is lower, because of the reduced thermal resistance.It is worth noting that the resulting solid fractioncomputed on the design domain Ωd differs from the optimal solution (31%) to Case 4 (48%).This is an interesting outcome because even if the Case 4 design has more solid material available, its displacement does not imply improved thermal diffusion.Lastly, the difference between thermal resistance values is slight.This confirms the need for searching a proper initial guess before running the TO.Indeed, running TO without exploring the field of possible initial solutions would risk experiencing a local minimum, whose value is really close to the global one.

Conclusions
This work investigates the effects of using different initial designs in the use of topology optimization (TO) for a microchannel heat sink.The CFD model is based on a pseudo-3D approach, which allows to obtain unintuitive design without resorting to 3D simulations.The initial design is made of multiple pin fins, and the starting solid fraction in constrained to 0.2.To explore different pin dimensions (varying the radius) and arrangements (varying the distance), an outer algorithm, i.e., a GA has been implemented.The aim of the TO is to find a final design able to minimize -with the same magnitude -the thermal resistance between the channel layer and the base plate layer, and the pressure drop.
Results confirm that both dimensions and displacements of pins have an influence on both objective functions and final design.In detail, starting with bigger pins leads to local minima, while using more pins with lower radii enables to reach a design with higher heat transfer and more aerodynamically shaped fins.

Figure 1 .
Figure 1.Sketch of 3D heat sink design and its reduction to a pseudo 3D model.

Figure 2 .
Figure 2. Parameters affecting the round pins size and displacement (rmin, smin) used as design variables for the GA.

Figure 3 .
Figure 3. Model validation results: design domain; velocity field, design plate temperature; baseplate temperature field.

Figure 4 .
Figure 4. left: GA population's objective function (ϕ) history: detection of the optimal solution and four individuals with worst ϕ, and related initial designs.Right: Design variables distribution.

Figure 5 .
Figure 5. Optimal solution and "case 4": initial guess, final design, temperature channel layer and velocity modulus field.

Table 1 .
Model results comparison.