Multi-objective topology optimization of pin-fin heat exchangers using spectral and finite-element methods

Forced convective pin-fin heat exchangers, due to the high wet surface area per volume and the hindered thermal boundary layers, feature low thermal resistances. However, the considerable coolant pressure drop, particularly for densely packed fin arrays, imposes operational costs for pumping power supply. This paper develops a multi-objective topology optimization approach to optimize sink geometries in order to minimize thermal resistance and pressure loss, concurrently. In accordance to the pin-fin geometrical characteristics, a dedicated pseudo-3D conjugate heat transfer model is utilized, by assuming periodic flow and fin design pattern, to reasonably reduce the high cost of full-3D model optimization. For the solution of flow part, a pseudo-spectral scheme is used, which is intrinsically periodic and features a high spectral accuracy, and the finite element method for the non-periodic conjugate heat transfer model. Exact partial derivatives of the discrete solutions are obtained analytically by hand-differentiation. This task is rather convenient for the flow part, due to the simplicity of the pseudo-spectral implementation; however, the MATLAB symbolic toolbox is selectively utilized for the finite element code to ensure an error-free implementation. Finally, the sensitivities are computed directly or via a discrete adjoint method, for the flow and heat models, respectively. To examine the present approach, two approaches are used for optimization of a practical cooling task: constrained and unconstrained multi-objective formulations, where in all cases more emphasis is placed on thermal resistance minimization. A series of optimized heat sink geometries via constrained or unconstrained multi-objective optimizations are obtained to examine practical utility of the present approach. The optimized topologies demonstrated superior cooling performances at lower costs of pressure losses compared to conventional (circular) in-line and staggered fins, and confirmed the supremacy of topology over pure sizing optimization.


Introduction
Alongside the trend of electronic devices generations, which is toward more compact designs and increased power densities, more advanced cooling technologies are continually demanded, to address the future industrial design requirements (Mudawar 2001). Many electronic components, such as central processing units (CPU), require a reliable heat dissipation, as well as a careful thermal Responsible Editor: Gengdong Cheng Ali Ghasemi a.ghasemi@tu-braunschweig.de 1 Institute of Aircraft Design and Lightweight Structure, Technische Universität Braunschweig, Braunschweig, Germany management, in order to maintain a prolongated lifespan while operating at the maximum capacities (Garimella et al. 2008). Such important cooling tasks have been often delivered to forced-convection heat sinks (FHS), in a wide range of applications for decades. Air and liquid-cooled sinks are two common FHS types, which the latter, due to the higher coolant heat capacity, is often preferred for intensive heat dissipation tasks. In liquid-cooled sinks, the coolant is pumped and circulated into the heat sink; hence, the cost of pumping power has to be counted for design and maintenance. As a result, the performance of such heat sinks are measured by two factors, simultaneously: heat removal capacity and required pumping power supply.
Heat sinks are designed in a variety of shapes and configurations, such as pin-fins and plate-fins. Pin-fin sinks are widely utilized as they feature high heat dissipation capacities, due to the disrupted and non-unidirectional thermal boundary layer development (Kim et al. 2008).
However, such thermal performance comes with the cost of high coolant pressure loss and increased pumping power (Ndao et al. 2009). It has been shown that the fin shapes and configurations play a crucial role in the overall performance (Abdoli et al. 2015). Accordingly, to achieve the highest performance, sink geometry has to be optimized. A series of studies have been carried out to enhance pin-fin geometries, by focusing on either sizeing optimization (Chen et al. 2005;Ndao et al. 2009) or shape optimization (Hajabdollahi et al. 2012;Tullius et al. 2012;Reddy et al. 2017). In these approaches, limitedly few geometrical aspects, such as fin diameter or spacing, are manipulated during optimization, and other aspects such as array configuration remain unchanged (see Alihosseini et al. 2020 for a review). Hence, a geometry optimization approach, capable of handling all geometrical aspects of pin-fins, is required. To this end, topology optimization techniques, which are able to concurrently optimize sizes, shapes, and configurations, can flexibly obtain optimized sinks without any limitations. Therefore, as the primary aim of this study, a dedicated topology optimization approach, considering the pin-fin geometrical aspects, is developed in this work.
Topology optimization (TO) was developed initially for structural optimization problems, using a homogenization theory (Bendsoe and Kikuchi 1988), and other techniques such as level-set methods (Sigmund and Maute 2013). TO techniques were soon applied in other applications such as flow and heat transfer problems. Development of TO for fluid flow problems was pioneered by Borrvall and Petersson (2003), who aimed to minimize power dissipation of Stokes flow using Brinkman penalization technique and introduction of a friction factor as the design variable to create solid region inside flow domain. Soon later, TO techniques were extended to steady Navier-Stokes (Gersborg-Hansen et al. 2005;Olesen et al. 2006), unsteady Navier-Stokes (Kreissl et al. 2011;Deng et al. 2011), turbulent flows (Othmer 2008;Kontoleontos et al. 2013), and fluid-structure interaction problems (Yoon 2010). Dede (2009) extended TO applications to conjugate heat transfer (CHT) problems, where heat is transferred between the solid zones (pure conduction) and the surrounding fluid (convection-diffusion). Later works suggested more utilities and potentials of TO by developing it for various multiphysics applications. The reader is referred to Dbouk (2017) and Alexandersen and Andreasen (2020) for a comprehensive review on TO studies, developments, and implementations, during the past decades with an insight to their trend and limitations, advantages, and disadvantages.
Utilizing TO for heat exchangers is computationally expensive, if solved in full-3D (Dilgen et al. 2018), particularly when it comes to compute sensitivities for every optimization step. The problem is more intensified for multi-objective analysis where a series of optimizations are required to capture a Pareto-like solution set. Pure 2D thermo-fluid models (Matsumori et al. 2013), on the other hand, are excessively simplified to properly simulate pinfin heat exchangers where heat transfer rate is a function of temperature difference between the hot base and finned area. To tackle this problem, McConnell and Pingen (2012) introduced a multi-layer (pseudo-3D) heat sink model, consisting a coupled fluid-thermal layer attached to a conductive thermal base-layer. Although their test cases were limited to a very coarse mesh (25×25), the obtained optimized multi-fin sinks are worth noticing. Later, Haertel et al. (2018) used a similar pseudo-3D model for the optimization of an air-cooled heat sink, and approved a satisfying agreement between the pseudo-3D model and full-3D simulation. Pseudo-3D CHT model has been also used for the TO of air and liquid-cooled heat exchangers (Zeng et al. 2018;Zeng et al. 2020a), and the outcomes have suggested promising performance achievements, both numerically and experimentally. Yan et al. (2019) improved the pseudo-3D model further by proposing a forth-order polynomial and linear temperature profile within finned area and base, respectively. Their study confirmed the consistency of the pseudo-3D with full-3D simulations, and the inadequacy of a single-layer model (pure-2D) to predict the temperature distributions. Zeng et al. (2020b) employed similarly a pseudo-3D model for transient TO of instantaneous cooling tasks, and reported an error of less than 1%, compared to full-3D.
In order to employ an appropriate and cost-effective CHT model for pin-fin type of heat exchangers, a similar pseudo-3D model can be used. However, the model should be adapted for the pin-fin sinks due to the geometrical aspects, and practical considerations. Pin-fin sinks, considering the manufacturing costs, are mostly designed with an array identical and repeated fins, which can simply scale up to hundreds of fins. Optimizing the whole sink geometry (as of previous studies) not only imposes significant computational burdens for optimization process, but also increases manufacturing costs. Instead, a single periodic (repeated) design geometry considerably reduces problem complexity, and is consistent with the pin-fin design convention. In addition, due to the rapid flow development of low Reynolds flows, the flow regime also forms a repeated pattern across the fins. Accordingly, a periodic flow assumption across fins is applied to significantly reduce computational loads, in favor of sensitivity calculations. For modeling the thermal part, a non-periodic computational domain is considered. But similar to previous studies, this domain is confined to a single column of fins (or design spaces) with symmetrical boundary conditions within adjacent columns, to predict thermal distributions spanning between the inlet and outlet, as well as the maximum base junction temperature. As a result, a modified pseudo-3D model, in accordance to the pin-fin sinks, is developed for this study.
Improving performance of a convective heat exchanger is essentially a bi-objective task, in the sense that the coolant pressure loss as well as the thermal resistance have to be minimization, simultaneously. Consequently, a series of designs in the form of Pareto-like solutions can be assessed as optimized heat sinks. The majority of studies, such as in Haertel et al. (2018) and Yan et al. (2019), considered minimizing thermal resistance for a fixed value of pressure loss, or vice versa, such as in Zeng et al. (2020b). Marck et al. (2013) were apparently first to develop a multi-objective topology optimization technique for CHT problems. They attained a convex Pareto-frontier with an aggregate objective function method, but reported some design limitations, such as flow blockage when thermal weighting is dominant. Subramaniam et al. (2019) also performed a similar study for flows at low to moderate Reynolds numbers using a continuous adjoint analysis. In contrast, they reported no numerical difficulties, namely flow blockage or broken path for thermal weightings up to 0.99. Also Li et al. (2019) developed a similar multi-objective TO tool for the cooling channel optimization and discussed the effects of weighting, Reynolds number, and initial guess on the Pareto-like solutions, supported by manufactured sink prototype experiments. Nevertheless, these studies developed a multi-objective optimization based on a singlelayer (pure-2D) model, which its deficiency and inaccuracy have been approved by several authors. Multi-objective topology optimization, using pseudo-3D CHT, has not been considered so far. Accordingly, it is aimed to employ the advantages of a modified pseudo-3D model to perform a rigorous multi-objective analysis, in order to provide a fully practical, flexible, and accurate approach towards optimal pin-fin sink design.
Pseudo-spectral scheme is used for the solution of the flow part of the present CHT problem. Spectral methods are powerful numerical schemes for approximation of flow, particularly on periodic domains (Canuto et al. 2012). This method is based on the Fourier transform of the discretized solutions from physical space, into Fourier space, in order to compute the spatial derivatives. Fourier spectral methods are generally more accurate than other numerical methods, such as FVM, particularly when the solution is smooth and periodic (Trefethen 2000;Naulin and Nielsen 2003). The utility of the pseudo-spectral with Brinkman penalization, as well as its advantages over FVM, namely its accuracy and flexibility, are studied by Kevlahan and Ghidaglia (2001) and Kadoch et al. (2016), respectively. The flow periodicity in the present pseudo-3D model allows using such a high accurate method, as it is consistent with the periodic nature of spectral schemes. For the solution of non-periodic thermal distributions, a standard finite-element method (FEM) is used for both convective (finned) and conductive (base) layers.
The exact partial derivatives for the sensitivity analysis are derived analytically by hand-differentiation of the discrete primal solvers; nevertheless, MATLAB symbolic toolbox has been selectively used for an error-free implementation. The discrete sensitivities are then computed directly for the flow part, and using a coupled discrete adjoint for the thermal part. Additionally, sensitivities are verified via a central-scheme finite difference method to ensure the accuracy and correct implementations. The globally convergant method of moving asymptotes (GCMMA) (Svanberg 1995) is utilized, as the optimizer. A practical heat exchanger problem with realistic physical parameters is considered for optimization with various design conditions, and is compared with the conventional circular pin-fin heat sinks (staggered and in-line). In summary, the novelties of the present work are listed as: 1. Development of the pseudo-3D CHT model with periodic flow and design assumptions, compatible with pin-fin sink geometries. 2. Utilization of a coupled pseudo-spectral/FEM numerical approach for the CHT model solutions. 3. Accurate sensitivity analysis of the flow and thermal parts, using analytical differentiation of the discrete primal solvers. 4. Performing a rigorous multi-objective topology optimization to capture a Pareto-like set of realistic pin-fin heat exchangers.
The remainder of this paper is organized as follows: in Section 2, the CHT model as well as the numerical methods are described, and in Section 3, the topology optimization formulations including the sensitivity analysis are presented. In Section 4, two classes of practical test problems are studied and quantitatively analyzed in details. In Section 5, further details of the findings are discussed, and finally, Section 6 summarizes the key features and achievements of this work.

Pin-fin heat sink model
In this part, a liquid-cooled pin-fin heat sink device and a simplified CHT model, corresponding the physical aspects of pin-fins, are described. The developed CHT model is based on a pseudo-3D approximation, with which few assumptions are employed to reasonably and accurately simplify the full model. The motivation is that in gradient-based optimizations, the sensitivity calculation is computationally expensive, and its cost dramatically increases for 3D simulations. In addition, for multiobjective optimizations, essentially not a unique solution exists and often a series of simulations are required to obtain a set of optimized. Therefore, the main concern is to avoid significant computational costs of full-3D models, by employing a simplified, yet sufficiently accurate, CHT model for topology optimization at reasonable computational costs. Figure 1a illustrates the full-3D heat sink (detached) model, consisted of three parts: hot plate, sink, and a cover. The hot plate generates uniform volumetric heat, which exclusively is transferred to the sink, attached above. The sink is enclosed with a cover to guide the liquid coolant from the inlet to outlet, and prevent leakage. The heat is absorbed from the hot plate and transfered to the cold fluid, flowing across fins. To force heat convection, the liquid coolant is pumped into the sink, but the pump is not depicted. External surfaces are assumed adiabatic, in the sense that heat dissipation is limited to the liquid coolant, as of the focus of this study. Thus, no heat is transferred to the surrounding environment via natural convection. It should be noted that no gap space is left between the cover and top of the fins, when the parts are attached.

Model description
Computational cost of using such full-3D sink model for topology optimization is significantly high, particularly when it comes to sensitivity computations for every optimization step. To alleviate this problem, the pseudo-3D model is utilized, which has a computational cost of nearly a 2D problem. This model, as shown in Fig. 1b consists of two parts, coupled together: base and convective layer. The base layer represents the hot plate and sink base of the full-3D model in Fig. 1a, merged together. The convective layer models cross-section of a full-3D with uniform temperature and velocity in height (z) direction. Base and convective layers are thermally coupled, in the sense that heat is generated in the base and absorbed from the convective layer.
Motivated by the flow and thermal characteristics of the presented pin-fin heat sink, the aforementioned pseudo-3D model ( Fig. 1b) is further simplified. Steady laminar flows, crossing the fin arrays, quickly develop a periodic pattern, due to the repeated pattern of fins. Hence, a fully periodic flow condition is imposed to the convective layer, in which flow is modeled sufficiently on a single periodic domain, referred to as Ω (see Fig. 2). It should be noted that as a In addition, for a sink consisted of an N c × N c fin array, only a single fin column (1 × N c ) is modeled with symmetrical boundary condition with neighboring columns. Accordingly, the base layer is split in cross-flow direction, matching the split convection layer. Figure 2 also illustrates schematic and boundary conditions of the simplified base and convection layers, referred to as Ω b an Ω c , respectively.

Thermo-fluid mathematical model
Flow of the liquid coolant across fins is assumed incompressible, laminar, steady state, and mathematicaly modeled using the Navier-Stokes equations and Brinkman volume penalization technique (Arquis and Caltagirone 1984;Angot et al. 1999), on a 2D doubly periodic domain (Ω). The penalized Navier-Stokes equations is as follows where u represents the velocity vector, P the hydrodynamic pressure, ρ f the density of fluid, ν the kinematic viscosity, χ the penalization function, u ∞ the upstream flow velocity, and the penalization parameter. The Brinkman penalization technique is a member of immersed-boundary methods, which suppresses the flow velocity anywhere it is enabled, only by modification of governing equations. It forms a nearly impermeable porous medium, which implicitly creates solid zones corresponding to fins, as shown in Fig. 2. The penalization effect, used to define solid zones, is handled by spacial function χ = χ(x), ranging between 0 and 1. The fluid is conveniently blocked (penalized) where Hence, the no-slip boundary condition on the solid-fluid boundary is automatically imposed when full penalization is enabled. The penalization intensity is controlled by parameter ( → 0 + ). It has been proved that the error bound of penaliziation, compared to conventional CFD methods with explicit solid-fluid boundaries, is proportional to 1/2 (Carbou et al. 2003). However, very small values of impose system stiffness due to instability of numerical methods. Hence, a sufficiently small value of , e.g. 10 −5 , should be chosen for a sufficient accuracy level and penalization effect.
It is important to note that flow is partially blocked, if 0 < χ < 1. Such cases create a porous medium with a continuous impermeability proportional to χ/ , referred to as gray material, later in the optimization process.
The thermal distributions over the convection (finned) layer, Ω c , and the base (heated) layer, Ω b , are mathematically modeled using the convection-diffusion, and pure diffusion equations, formulated as and, where, T c and T b are temperature distributions over convection and base layers, respectively; c p is the coolant specific heat capacity,q tr the heat transfer flux rate between layers, V b the base volume, andQ gen the uniform heat generation rate in the base. z c and z b are the heights of the heat sink base and fins, respectively. It should be noted that (3) represents thermal distributions simultaneously over both solid and fluid regions of Ω c . Inside fins, velocity is zero and κ equals to k s , the solid thermal conductivity, whereas in the flow region, κ is equal to k f , the fluid thermal diffusivity.
The heat flux rate ofq tr is defined aṡ where H is a function which determines the thermal coupling magnitude between layers. H = h s where fins are pinned to the base, and H = h f where coolant flows over the base. Accordingly, h s and h f are the thermal coupling coefficients between the solid-solid (base-fin) and solid-fluid (base-coolant) contacts, respectively. The thermal boundary conditions are given as follows: T c = T ∞ at the inlet, representing the coolant initial temperature. On symmetrical or adiabatic boundaries, as specified in Fig. 2, zero normal temperature gradients are imposed, such that where n is a vector, perpendicular to the boundary. All material properties are assumed constant, and invariant to temperature variations. This assumption leads to a oneway coupling between flow and thermal equations, in the sense that the velocity solution is required for calculating temperature (in convection term), but not vice versa.

Numerical methods
The penalized Navier-Stokes (1-2) are solved using the powerful pseudo-spectral scheme (Kevlahan and Ghidaglia 2001;Canuto et al. 2012). This method is based on the discrete Fourier approximation of the solutions, in physical space, defined in Fourier space via where, k = (k x , k y ) is the wave vector, defined as the coordinates of an N 2 uniformly discretized Fourier space (k x , k y ∈ {−N/2 + 1, −N/2 + 2, · · · , N/2});î is the unit imaginary number; u = (u, v) is the physical velocity approximation, and correspondinglyū is the transformed velocity.
It can be seen that the required spatial partial derivatives in (1) could be calculated conveniently by exclusively differentiating the exponential term in Fourier space.
The 2D discrete Fourier transform could be performed for a given scalar field, ϕ(x), defined on an N-by-N uniformly discretized domain using transform matrices. To do so, the square matrix M ϕ is defined such that m i,j = ϕ(x i,j ), and the transform is operated via whereφ is the transformed value, and D is a complex matrix, with elements defined as It should be noted that for the transform of (8), a fast Fourier transform algorithm, such as FFTW3 (Frigo and Johnson 2005), could be used for better efficiency. By projection of (1) onto the divergence-free space, the continuity (2) is satisfied. Then, the penalized Navier-Stokes become where ω is the flow vorticity vector, and Λ = (Λ x , Λ y ) is defined equivalent to time-derivative of u, which is supposed to be zero for steady-state solutions. The procedure of pseudo-spectral is as follows: first Λ in (10) is computed via discrete Fourier method, then u is obtained using a pseudo-time-marching loop, via until convergence is achieved, i.e. Λ = 0. To compute Λ, the vorticity field is first computed by whereD is the complex conjugate transpose of D, used for the inverse transform, and • is the Hadamard matrix product operation. Then, the convective, viscous, and penalization terms are computed in Fourier space by Then, the divergence-free projector, P, is applied point-wise for each (Γ x , Γ y ) values, via As the final step, Λ is derived via an inverse transform of (14) back into the physical space It is important to note that the preceding pseudo-spectral steps are based on the basic matrix operations, which noticeably simplifies computation of the solutions and sensitivity analysis, which will be discussed later.
The same periodic flow field (Ω) is used as the design domain, according to the periodic fin patterns. The 2D design space is then defined as of a set of control points, i.e. γ i,j ∈ Ω, assigned to each domain discretization point x i,j . These control points are bounded independent variables (0 ≤ γ i,j ≤ 1), which are used to define solid-fluid parts of the domain. More precisely, the material-dependent parameters in (1), (3), and (5), namely, penalization parameter (χ), thermal coupling coefficient (H ), and thermal diffusivity (κ) are defined as functions of γ . Such functions are discussed later in part Section 4.1.2. In this manner, the sink geometry is exclusively handled via control points, in the sense that γ = 1 refers to solid material and γ = 0 to areas filled by fluid. Hence, γ 's are used as design variables in the topology optimization process.
For the solution thermal part, a standard finite-element method (FEM) is utilized. The strong forms of (3)-(4) are multiplied with appropriate test functions, integrated over Ω c and Ω b , and by performing integration by part, the matrix forms of the corresponding FEM system of equations are derived as and where, K c and K b are the global stiffness matrices corresponding to the convection-diffusion and conduction terms in (3) and (4), respectively; T c and T b are vectors of thermal solutions in finned and base domains, H is a matrix associated with heat transfer coefficient, H , and Q is the vector of heat source. R b,c denote the residuals of the system of equations, supposed to be zero once the solution is converged. The non-periodic thermal domains are discretized by 2 × N 2 × N c structured triangular elements, where N c is the number of fins in each column. This implies the periodic flow field u is mapped N c times onto Ω c , as required to construct K c . In summary, the solve procedure is as follows: for a specific geometry, the design-dependent parameters χ, κ, and H are defined; then, flow velocity solution is computed via the pseudo-spectral scheme on the periodic domain of Ω, and then mapped repeatedly onto the thermal domain Ω c ; next, the conjugate thermal solutions (17-16) are computed via FEM, with the provided flow field. Therefore, it is important to note that the solutions exclusively depend on the sink geometry, or equivalently, the design variables γ . It is then the task of topology optimization to design optimized geometries, by varying γ .

Optimization formulation
As mentioned earlier, the overall performance of a forcedconvection heat sink is evaluated by two factors: thermal resistance and pressure loss, denoted here as R and P , respectively. For a heat sink, a lower R value is demanded, since that measures the temperature rise of the hot zone, compared to a reference temperature, with a prescribed heat generation. In this work, it is defined as where,T b denotes the average temperature in the base plate, and T ∞ stands for the coolant temperature at the inlet. The secondary factor, the coolant pressure loss, P , is regarded as a cost for the cooling task, since for a fixed coolant flow rate, it determines the pumping power to be supplied externally. Hence, lower the P value, lower the required pumping power. To calculate the pressure loss, first the drag force (acting on each fin in the periodic cell) is calculated by integrating the Brinkman penalization term over the entire periodic domain (Angot et al. 1999), via Then, the total pressure loss is derived using the summation of drag forces on a column of fins, denoted as where, A i equals the inlet area. Accordingly, to enhance a heat sink performance, two scenarios are possible. In scenario 1, one may keep one factor, e.g. P , fixed and minimize the other one, R, or, try to minimize both factors simultaneously, regarded as scenario 2.
With respected to these scenarios, two class of optimization problems are considered in this work. For scenario 1, a constrained single-objective optimization problem is formulated as where P * denoted the maximum permissible pressure loss, and for scenario 2, an unconstrained multi-objective problem is formulated as where, β ∈ [0, 1] is the weight factor, and both weighted objectives are normalized with respect to their initial values. In scenario 2, the target is to minimize J , and β determines the importance of pressure loss during optimization.

Sensitivity analysis
Accurate sensitivity analysis is a crucial step for the gradient-based topology optimization. The goal is to compute the derivatives of the aforementioned objective/constraint functions, R/P , to track their changes with respect to the design variables during optimization process. For the flow part, total derivative of the drag force, dF Ω /dγ , is derived by differentiation of (19), as where, γ , χ, u x , and u y are vectors for γ , χ, u, and v values, respectively. The partial derivatives with respect to F Ω are calculated simply by differentiation of (19), and dχ/dγ is calculated, according to the choice of interpolation function, which is discussed later.
To obtain du x,y /dχ , (11) is differentiated with respect to χ , as ⎡ where, I denotes the identity matrix. To derive the partial derivative matrices of Λ x,y with respect to χ and u x,y , the chain rule of differentiation is applied to the pseudo-spectral scheme steps, i.e. (12-15). For example: and, where the exact partial derivatives are obtained analytically by hand-differentiation. Now, to compute du x,y /dχ, (24) is reformulated in a succinct form, by where, X 0 is zero, since the initial velocity is fixed for the forward solve process. To reduce computational costs, the matrices A and B are only calculated with the converged flow solution; therefore, their subscript n is dropped. As a result, X n+1 can be formulated as where, by assuming n + 1 = 2 k and applying an algebraic factorization, it becomes: where, in practice, X converges quickly (after 10 iterations), as n increases exponentially with k steps. It should be noted that this technique relies only on basic matrix operations, and its main computational costs, which are B products, are independent of number of design variables (please refer to Ghasemi and Elham 2020 for more details).
For the sensitivity of the thermal part, a coupled adjoint analysis is applied to the conjugate system of (17-16), via a Lagrangian approach, defined as where λ c and λ b are Lagrange multiplier vectors. The sensitivity of the Lagrangian becomes where zero terms, such as ∂R/∂γ , have been already ignored. All the partial derivatives in (31) are analytically derived via hand-differentiation of the FEM code; however, the MATLAB symbolic toolbox has been selectively employed to avoid implementation errors. The total derivatives of du x,y /dγ are provided from (24), as discussed earlier. However, the total derivatives of dT b /dγ and dT c /dγ are principally expensive to compute. To avoid this issue, adjoint problem of (31) can be rearranged as: where the last two terms can be canceled by proper values of Lagrange multipliers. Therefore, the following system of equations has to be solved: ⎡ and plugged into (32), to obtain sensitivity of the thermal part. The optimization is performed using GCMMA ((Svanberg 1995)), which is designed for large-scale nonlinear gradient-based optimization problems, particularly with large number of design variables. The summary of the flow and thermal solutions, sensitivity calculations, and optimization process is summarized in flowchart (Fig. 3).

Numerical optimizations
A water-cooled heat sink, made of aluminum, with the physical parameters listed in Table 1 is considered for optimization. The problem parameters such as generated heat, as well as the sink geometries, are purposely defined to be similar to a high-power CPU cooling task. All physical parameters are assumed invariant with respect to thermal changes and constant during optimization. As mentioned earlier, a single fin column (Fig. 2) with periodical boundary conditions between fin columns in row-wise direction is assumed for lower computational costs. Accordingly, the generated heat in this portion of the base plate is used for calculation of thermal resistance, while one may scale it for the overall thermal resistance of the full sink, using since Q gen of a single column is 1/N c of the total Q gen . The single − column subscript of thermal resistance is dropped for sake of brevity. In the following, some numerical examples as well as some implementations considerations are discussed. The Phoenix computing facilities of TU Braunschweig have been used for executing the optimization examples.

Implementation notes
Before proceeding to numerical examples, two important aspects, pertaining to the present approach, are discussed: first, the pseudo-3D solver validation, and second, the choice of interpolation parameters, which the latter is less addressed in the literature.

Comparison of pseudo-3D with full-3D simulations
The accuracy of the flow solutions obtained from pseudospectral method, particularly the drag force acting on the fins, or equivalently the pressure losses, has been previously verified by the authors, for incompressible laminar flows (Ghasemi and Elham 2020). The thermal solutions, however, are eccentrically required to be analyzed and verified before proceeding to optimization. In this regards, the Ansys FLU-ENT, which is a finite-volume based commercial software, is used to model heat sinks in full-3D scale. An in-line cylindrical pin-fin configuration with four different cross-section diameters of 4, 3.5, 3, and 2.5 (mm) are considered for comparison and solved via both Ansys and pseudo-3D. For the best accuracy of the full-3D model, body-fitted meshes with nearly 3.3 million cells are generated and utilized. It should be noted that it takes, on average, 20 min for the full-3D model to converge (on a 4-core, 4.2 GHz processor), whereas it takes less than 15 s for the pseudo-3D model to converge with total of 0.148 million degrees of freedom. Figure 4 shows the temperature distributions in full-3D, pseudo-3D, as well as the temperature profiles along x direction, in the base plate. It is clear that for the case of D = 4 mm, the pseudo-3D solution perfectly matches the full-3D. As the fin diameter reduces, e.g. D = 2.5 mm, the pseudo-3D tends to slightly overestimate the base temperature near the entrance. The reason is that the flow field in the pseudo-3D model is periodic and therefore fully developed across the entire domain. However, in practice, the flow regime is developing near the inlet, and after a certain distance becomes developed. Flow regime develops more rapidly, for the case in which the fin diameter is relatively large, due to the small gap between fin columns. On the other hand, it takes longer for the flow to fully develop, when fin diameter is relatively small, in the sense that the developing distance is larger. It is also known that the heat transfer rate of the developed laminar flow is rather smaller than the developing flow. Therefore, the pseudo-3D may locally, near the inlet, demonstrate a limitedly lower heat transfer rate, due to the fully developed flow profile, Fig. 4 a Comparison of pseudo-3D and full-3D thermal distribution in the base plate. The curves illustrate averaged temperature profiles in y-z plane, along x direction. F3D and P3D stand for full and pseudo-3D, respectively. b Full-3D thermal distributions using Ansys Fluent. c Full-3D temperature at 50% fin-height cross-section. d T c in convection layer, and e T b in base layer although for the rest of the domain, it is well consistent and satisfyingly accurate. This effect can be observed in the entrance of the finned area in Fig. 4c, where the coolant is slightly heated in the space between the inlet and the first fin. Table 2 compares the calculated thermal resistance, obtained by both approaches. Considering the relatively low computational cost of pseudo-3D compared to full-3D, its error is satisfyingly small and negligible. Such consistency and accuracy level confirms the utility of the present approach for analysis and optimization of the pin-fin sinks.
The full-3D results are also used to quantify the thermal coupling parameters of h s and h f , as required for the coupling factor H , in (5). In the literature, these coupling factors are defined either heuristically (Haertel et al. 2018), or using an analogy of full-3D solvers (Zeng et al. 2018), such that where A f and A s are non-finned and finned areas over the base, andh values can be obtained from either full-3D simulations or experiments. However, the problem is that the obtained values are valid only for a unique geometry. In topology optimization, the main difficulty is that the geometry changes in every optimization iterations, which accordingly, the heat transfer coefficients changes continually. In addition, for the intermediate geometries during optimization, where the design resembles a porous medium, an explicit solid-fluid boundary does not exist to calculate the corresponding solid-fluid convection rate.
To tackle this problem, the full-3D results from several geometries are correlated for an estimation of the convection heat transfer coefficients for optimization, as listed in Table 1. It is important to note that h s /h f ≥ 100, which means that the heat is predominantly absorbed via fins, and therefore the value of h s determines the base average temperature, or equivalently the thermal resistance. However, in practice the optimization process is almost insensitive to the h s value. It is because the optimizer attempts to minimize R/R initial ratio, where it is independent of h s itself. For the accurate values of h f and h s , after optimization, one may consider a full-3D simulation of the optimized fin geometry to obtain the accurate temperature profile in the base.

Interpolation parameters tuning
In density-based topology optimization approaches, it is very important to control behavior of the gray materials (intermediate designs). For example, in structural topology optimization, the elasticity of the intermediate designs are intentionally reduced via a material penalization approach (SIMP), to achieve a full solid-void design with negligible grayness (Bendsøe and Sigmund 2003).
In the present problem, three design dependent parameters, namely χ, H , κ, require a proper interpolation strategy. Motivated by the interpolation function suggested by Borrvall and Petersson (2003), the design-dependent parameters are defined as: and, where I a , I b , and I c are fixed parameters to control the convexity of the interpolation functions for the intermediate design variables during topology optimization process. The interpolation parameters are defined generally heuristically and are less discussed in the literature, although, it has been turned out that they play an important role for the   Optimized fin topologies, for constrained optimization, in a periodic setting. Imaginary red lines specify the periodic design domains optimization process and final results. In the present work, the process in which these parameters are set is included. Figure 5a demonstrates the variation of the penalization parameter χ as a function of a design variable γ for various I a values. It can be seen that the convexity of the curve is controlled by I a , i.e. smaller the I a value, more convexity of χ, and vice versa. The same behavior applies to the variation of H and κ.
To analyze the effect of the interpolation parameters, a cylindrical fin with various grayness is simulated. More precisely, the set of design variables are scaled, between 0 and 1, via γ gray = scaling-factor × γ f in , to study and control the behavior of the intermediate (gray) designs. Figure 5b shows how pressure loss changes in the intermediate designs for various interpolation parameters of I a . It is found that an upward convex variation works better, since pressure loss in any case is the considered as a cost for the design which should be reduced. In this manner, I a = 10 produces the best results since lower values generate large pressure losses initially, and consequently cause oscillations in the optimization process.
In contrast, for thermal resistance, a downward convexity works better, since in all optimization cases lower R is desired (Fig. 5c), and such convexity promotes the optimizer to create/expand the design in empty parts of the domain. In practice, it is observed that I b = 10 is an appropriate value. For interpolation of thermal conductivities, κ, as shown in Fig. 5d, thermal resistance is rather insensitive to I c . However, it has been found that a relatively large value between 10 3 -10 4 is more appropriate, while lower values lead to designs with gray solid-fluid boundaries. The reason is that the convection heat transfer rate, which is proportional to the fluid thermal diffusivity (k f ), is boosted in a gray boundary. Such gray boundaries permit fluid to flow through for a convective process, where a linear interpolation between k f and k s yields a large conductivity that artificially increases the convection rate. To avoid this problem, I c has to be increased such that the thermal conductivity of the gray materials are closer to k f unless where γ → 1.

Sensitivity verification
A crucial step for gradient-based optimization is verification of sensitivities. This step is initially carried out to ensure the correct code implementation as well as validation of sensitivity analysis accuracy. For this purpose, the sensitivity of the multi-objective function J (in (22)) is approximated numerically via finite difference method (FDM) and compared with the analytical sensitivities. The central-scheme FDM sensitivities are derived by: where δ is the perturbation step size. As a test case, a problem with a coarse discretization is considered and for each design variable the relative error between analytic and FDM is computed by: For verification purpose, the same example in part Section 4.1.1 is used with a coarser discretization (11 2 design variables). Figure 6 illustrates the verification results for three values of weighting indexes (β), obtained by a proper value of δ. It can be observed that sensitivities computed via different approaches are perfectly matching for the entire design domain, which confirms the correctness of the implementation and analysis. In addition, the gradients for various β values are satisfyingly accurate, as required for the multi-objective optimization approach in the present study.

Scenario 1: constrained optimization
In this example, the constrained optimization formulation of (21) is considered, where the constraint is the maximum permissible pressure loss. This is important from practical point of view, since in many applications the pump power supply is fixed or limited. Here, P * is set to 0.1109 (Pa), which is the pressure loss of an in-line cylindrical pin-fin with 25% of fin volume occupancy (D = 2.82 mm). A set of optimization cases are studied, namely, P ≤ nP * , where n ∈ {1, 2, 3, 4, 5, 6}. The goal is then to minimize thermal resistances with respect to the limit of P .
No volume constraint is imposed here. From Fig. 5, it is clear that by formation of a fin, R decreases while P increases. This means that the existence of a fin is in favor of the objective function and the optimizer tends to use more solid material in the design domain. On the other hand, the trade-off of using occupying the sink with fins is a high pressure loss. Therefore, volume constraint is not required to limit the solid-volume occupancy.
The optimized fin topologies for the six cases are illustrated in Fig. 7, and their corresponding thermal resistances are listed in Table 3. For the first case with the lowest permissible pressure loss, the optimized topology consists of three fins with low frontal cross-section area and high aspect ratios. It has been attempted to reduce drag by keeping a large gap between fins to guide the flow mainstream. Accordingly, as plotted in Fig. 8, the first case also features relatively the minimum of the maximum flow velocities. From thermal point of view, this case has the maximum thermal resistance since the fin wet surface is relatively small. As a results the base temperature increases, due to the limited heat transfer rate of finned area (see Fig. 10), and consequently, fins temperatures are higher, as demonstrated in Fig. 9a. It is interesting that the optimized design has a staggered configuration, since the staggered arrays due to the interrupted thermal boundary layer and better thermal mixture produce higher heat transfer rates.
By relaxing of the pressure loss constraint, the number of fin elements is increased up to 10 unique fins per periodic domain, in order to increase the total wet surface area and to better mix the cold coolant with the hot fins. This certainly increases the viscous drag force since the fluidsolid contact as well as the average velocity increases (see Fig. 8f). In addition to the enhanced heat transfer, the overall fin volume is increased and is spread almost everywhere in the domain. The optimized designs attempt to not only remove heat from base, but also manage the temperature on the base plate to avoid local temperature rises. As a result, the design of case 6, features 38.2% less thermal resistance, or approximately 16 (K) less average base temperature.
An important aspect of the topology optimization, which essentially has to be investigated, is the effect of initial  The base temperature, averaged in y-direction, plotted along x-direction for the constrained optimization cases design before optimization process. It is ideal that the final optimization solutions are insensitive of the initial state; however, due to the complexity of the multi-physics problems, the final solutions can be triggered by that. To investigate this effect, the optimization problem of case 2 (P ≤ 2P * ) has been studied with six different initial designs, as shown in Fig. 11b. Designs (a)-(e) are different in shape or topology, but with equivalent solidoccupancy (25%). Design (f), in contrast, has a very low solid-occupancy (0.83%). Figure 11a shows the optimization history for each individual initial designs. It is clear that the optimization process of design (f), which is a very thin cylindrical fin, achieves the minimum thermal resistance. And surprisingly, this case has the best convergence rate, and in practice, it has been found to be the most reliable case to start the optimization with. The reason is that the design of (f) can be considered as a minimal design, in the sense that the majority of the domain (more than 99%) are left empty for the optimizer to perform the design. In all other cases, the optimizer is more or less affected by the initial flow and thermal solutions. It has been also noticed that starting from fully empty domain leads to high oscillations in optimization due to the extreme variation of the objective function. Therefore, in the present paper, the initial design of (f) has been used for all cases, to begin the optimization process with.

Scenario 2: unconstrained multi-objective optimization
For this set of examples, the unconstrained optimization formulation of (22) is considered. Similar to the previous examples, no volume constraint is imposed. It is because Fig. 11 Studying the effect of initial design on convergence, in order to find the best case to start the optimization from one side, the optimizer needs a minimal design to cool down the hot base plate, and from other side, full solid design blocks the flow (no convective cooling) and dramatically increases pressure. Therefore, volume constraint is unnecessary for this case. Particularly in this multi-objective optimization problem, β determines the importance of the pressure loss for optimization. Several optimization cases are considered here, in which the objective function weights are chosen differently, i.e. β ∈ {0.001 0.01 0.025 0.05 0.075 0.1 0.25}. The rest of the parameters are chosen similar to the previous example.
The optimized fin topologies, corresponding velocity fields, and temperature distributions are shown in Figs. 12, 13, and 14, respectively. For β = 0.25, the sink topology consists of only two hydrofoil-shape elements per periodic cell, to keep the drag as low as possible. Fin columns are well distanced such that through the remained gap space, a straight jet-like flow stream is formed, to better guide the flow and keep the flow average velocity enough low (see Fig. 13a). Such a 45 • staggered alignment better spreads fins over the entire hot base for an improved cooling effect and reduced thermal resistance; nevertheless, the thermal boundary becomes fully developed (see Fig. 14a). In contrast, as β reduces, number of fin elements per design cell increases, in order to maximize solid-fluid contact surface and mixture of the cold coolant with the hot film layer around fins. Moreover, increased number of elements increases heat transfer between fins and the base, leading to a better overall cooling performance and thermal management. As expected, the trade-off of increased heat transfer performance is the increased pressure loss. From Fig. 13, it can be seen that the more number of elements, the less gap space between fins and consequently the higher flow velocity. The increased flow velocity and fluid-solid contact leads to higher viscous drag forces. Figure 15 displays the base temperature profile, thermal resistance, and pressure loss for various weights. The R/P profile indicates that at a very low cost of pressure loss (P = 0.13 P a), the thermal resistance is still acceptable for many applications, in the sense that the hottest base temperature is only 43.4 K above ambient temperature. It is 57.3 K for P = 0.067 P a which can be considered for cases where power supply is highly limited. On the other side, the pressure loss increases considerably to achieve the lowest thermal resistance. At the cost of P = 2.06 P a, the base plate is enormously cooled down such that the hottest point is interestingly 27.8 K above ambient temperature, for such a heavy cooling task.
It is worth noting to the geometry of case of β = 0.001, which generates the maximum cooling effect. As illustrated in Fig. 16, in every periodic cell, the optimized fin array utilizes two different techniques to enhance the heat transfer rate. At the zones, denoted as A, three nozzlelike channels are created, where beside the prolongated solid-fluid contact surfaces, the flow is accelerated up to 3.5 times that the upstream velocity. Such increased flow velocities, regardless of the high pressure losses, enhance the convective heat transfer rate. Once the flow is accelerated, it is guided through the next zones, denoted as B. In B-zones, coolant crosses three sub-domains containing three small fins each. The heat convection in these zones are highly boosted because of the high flow velocity, and certainly the increased wet contact surfaces. This can be observed from the negligible temperature difference between the fins and surrounding coolant. It is also noticeable that regardless of the low pressure loss index, the small fins feature a hydrofoil shape to further minimize Fig. 15 Analysis of thermal resistance, pressure loss, and base temperature for various multi-objective weights β drag. Such a design can be used for applications where the power supply is abundantly provided, and a high cooling performance is desired.

Discussion of results
From practical point of view, it is beneficial to compare the results of the present topology optimization approach, with the conventional sizing optimization approach. In sizing optimization, the sink design is optimized with respect to some parameters such as fin diameter or cross-section length; however, shape, alignment, and topology remains unchanged. In this regards, a series of (cylindrical) pin-fin heat sinks are simulated with various solid volume fractions, ranging from 0.05 to 0.6, with an in-line and staggered configurations, to compare with the optimized topologies, obtained in this study. In Fig. 17, such results are put together on a multi-objective (R/P ) basis.
The most important point is that all the pin-fin heat sinks designed by topology optimization outperform the conventional designs in terms of thermal resistance and pressure loss, and are laid on a Pareto-like curve. To better clarify the difference, one can observe that the optimized sinks produce much less pressure loss than conventional sinks with equivalent thermal resistances. For instance, the sink of β = 0.01 generates 11.8 times less pressure loss that the staggered cylindrical fins to achieve to thermal resistance of R = 0.8 [K/W]. Moreover, the optimized sinks achieved low thermal resistances such that conventional sinks can either achieve with an increased cost of pressure loss, or never can reach to. As an example, the sink of β = 0.001 has 18% less thermal resistance (4.2 [K] lower base temperature average) compared to the staggered cylindrical fins with equivalent pressure loss.
In conventional pin-fins, particularly for in-line arrays (see Fig. 4c), due to the flow circulation in the wake, a portion of coolant is trapped, which deteriorates the Fig. 17 Multi-objective analysis and Pareto-like solutions convection rate. In the optimized sinks, fin elements mostly feature a high-aspect ratio hydrofoil, not only to reduce drag, but also to eliminate wake. This causes an improved convention rate, as coolant is never trapped.
In sizing optimization, the design is limited to a fixed topology, whereas it was shown that the topology and number of repeated fin elements changes, depending on the permissible pressure loss, in topology optimization approach. Additionally, sizing optimization has limitations in design since the heat sinks may experience flow blockage if the fins volume ratios are increased. However, it was shown previously that in topology optimization, flow blockage is automatically avoided.
It is also noticeable that the unconstrained optimization results , although negligible, are slightly better than constrained cases. In constrained mode, the optimized faces some limitations, where in unconstrained mode, no limitations exist on design in any case. The constrained approach is important for practical applications, while the unconstrained approach is valuable for multi-objective analysis of the heat sink, since there in no direct control on pressure loss with β.

Conclusions
This paper develops a numerical tool, using a topology optimization approach, particularly for the design and optimization of forced-convection liquid-cooled pin-fin heat sinks. To avoid the significant computational cost of a full-3D model, a pseudo-3D conjugate heat transfer model with periodic flow condition was utilized. Flow and thermal solutions were numerically computed via pseudo-spectral scheme and finite element method, respectively. A multiobjective optimization approach was employed, as the performance of forced-convection heat sinks is measured simultaneously by thermal resistance as well as coolant pressure loss. Two classes of constrained and unconstrained optimization formulations were studied, and a series of optimized heat sinks were designed, accordingly.
Hence, the following conclusions could be drawn: 1. It was shown that the present pseudo-3D model has less heat transfer rate, locally near the inlet, due to the periodic flow boundary condition. But, overall, its thermal distribution in the main portion of the simulated domain as well as the predicted thermal resistances are satisfyingly close and consistent with the full-3D results. 2. In all test cases, a cylindrical fin (with 1/120 solid volume ratio) was used as the initial design, as it demonstrated the best convergence rate and overall performance. 3. The layout of the optimized structures is connected with the pressure loss, in the sense that with the larger permissible pressure loss, the number of fin elements per periodic cell increases, to better mix the thermal boundary layers and increase the fin-coolant contact surface. 4. In most of the cases, the optimized fins feature hydrofoil shapes, to better reduce the drag force. Fin elements are spread evenly with a staggered configuration to produce a uniform cooling effect over the base. In addition, such staggered format has better convection rate due to the frequent interruption of thermal boundary layer. 5. In the multi-objective analysis, the optimized sinks expectedly featured lower thermal resistance and higher pressure loss, as the weight of pressure loss (β) was reduced. For the case of β = 0.001, a high performance sink was obtained, where the flow was periodically accelerated by nozzle-like channels and injected towards groups of small fins to boost the convective heat transfer rate. 6. The optimized sinks showed considerably higher performances, compared to the conventional in-line and staggered designs, by producing a Pareto-like behavior. In sizing optimization, design is limited due to flow blockage, whereas, in topology optimization, flow blockage is avoided, due to the significant rise of pressure loss.
Acknowledgments The authors would like to thank Karthikeyan Kurinjivanan for his help to prepare Ansys simulations.
Funding Open Access funding enabled and organized by Projekt DEAL.

Conflict of interests
The authors declare that they have no conflict of interest.

Replication of results
The prototype MATLAB script files will be available upon request.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.