Reduced Order Design Optimization of Liquid Cooled Heat Sinks

The rapid densification of electronic systems as evident from the International Technology Roadmap for Semiconductors [1] has increased the demand for effective heat dissipation. Thermal management is a roadblock to further miniaturization of electronics considering safe operating condition requirements as well as performance degradation at elevated temperatures within active and passive components. Thermal and hydraulic performance along with the volume and mass of the thermal management solution at the component level as well as the auxiliary component and system level has inspired multidisciplinary, multi-objective optimization for electrothermal systems [2–5]. Effective cooling performance and simplicity have popularized the application of single-phase liquid cooled thermal management approaches for electronics cooling. These include channel flow heat sinks, as well as jet impingement cooling [6–8]. Seminal work on microchannel heat sinks by Tuckerman and Pease [9] demonstrated the potential of liquid cooling solutions for high heat flux removal. Enhancement in the heat transfer performance mapped in terms of the low thermal resistance of 0.1 (° C· cm²)/W comes at the cost of high-pressure drop approaching 2 bar, thus increasing demand for pumping power. Many researchers have investigated liquid-cooled heat sink performance for a range of design and operating conditions using experimental [10–12], analytical [13–16], and numerical approaches [17–19]. The results of these studies have guided the design process for liquid-cooled systems. Parametric modeling [20–22] of performance parameters with respect to design variables including geometric dimensions and flow boundary conditions has provided means for design optimization [23–27]. Experimental approaches have yielded the most accurate results for performance modeling but are restrictive in terms of the number of design points that can be considered due to the cost and time required in setting up the experiments for different design variations. Numerical simulations using fluid flow and conjugate heat transfer analysis enable the investigation of the necessary performance parameters. Numerical studies however need to be supported with mesh refinement and experimental validation for ensuring the accuracy of the obtained results. Even though parametric studies are a proven tool for design optimization, their scope is limited in terms of design variables available for consideration. This limitation often fails to capture the complete flexibility within the design domain. Implementation of custom channel designs and flow layouts for increasing the heat transfer performance by enhancing flow mixing [28–34] is instrumental for designing high-performance thermal management systems. Most past studies are constrained in a way that leads to optimization of hyper-specific designs. While providing utility at optimizing those designs, there is a need for additional research that allows for a broader search of the design space and consideration of different design architectures.

Inspired by exciting developments in its application to solve problems in structural mechanics, topology optimization (TO) has been adapted for designing high-performance heat sinks [35–37]. More recently, TO has been adopted for the design of liquid cooled heat sinks [38–43]. TO uses numerical simulations and sensitivity calculations for every iteration in the optimization search, which can be computationally expensive and time-consuming. Published research investigated design improvement and optimization of liquid cooled heat sinks most focuses on overall heat transfer performance at a given boundary. In practice, liquid cooled heat sinks can be subject to spatially heterogeneous temperature profiles, resulting in nonuniform heat fluxes and thermal stresses within electronics packages and reduced reliability [44,45]. Although the challenges associated with device-to-device temperature gradients have been addressed in specific studies focusing on targeted hot-spot cooling [46–48], seldom do these studies entail design optimization. Electrothermal codesign applications use reduced-order thermal performance estimation [4]. Temporal and computational resource requirements of existing thermal design optimization approaches restrict their integration with electrothermal codesign packages.

This paper develops a design optimization tool for single-phase liquid cooled heat sinks. The fluid flow layout is broken down into flow blocks that form elements divided into a coarser two-dimensional (2D) grid. First, the governing mass, momentum, and energy equations for fluid flow and heat transfer within the fluid and solid domains [48,49] are solved to obtain the performance of the design in terms of thermal resistance, pressure drop, and temperature deviation on the heated surface. Next, an optimal fluid flow layout is found using the optimal path search (OPS) algorithm based on the optimization of the thermal-hydraulic performance parameters. The flow layout geometry is then further optimized by considering variable hydraulic diameters using a gradient-based optimization solver. The layout and geometry optimization are then integrated to form an overall design optimization tool. Performance modeling coupled with our rapid design method resulted in cold plate design optimization on a time scale of 60–120 s depending on design complexity. All computations were performed on an Intel^(R) Core™ i7-10700 CPU at 2.9 GHz computer using serial computing. A case study was conducted on a SiC power module, demonstrating a 25% reduction in junction-to-coolant thermal resistance for our optimized design when compared to a commercially available aggressive single-phase cold plate solution at the identical 16 kPa pressure drop.

where $din, vin$⁠, and $dout, vout$ are the channel hydraulic diameter and average flow velocity at the inlet and outlet faces of the flow block, respectively, $ΔP$ represents the pressure drop within the flow block calculated using $f$⁠, $ρ$⁠, $L$⁠, $d$⁠, and $v$ which denote friction factor, fluid density, block size/length, average channel hydraulic diameter (⁠$d=(din+dout)/2$⁠) and an average of the inlet and outlet velocities (⁠$v=(vin+vout)/2$⁠), respectively.

where $Re$ is the Reynolds number based on average diameter and fluid flow velocity (Re = $ρvd/μ$⁠, where $μ$ is the dynamic viscosity of the fluid). The limits on the diameter values prevent drastic changes in flow properties within a single flow block supporting our assumption of a Reynolds number based on average values. Here $ϵ$ is the roughness of the internal surface of the channel. For our simulations, we assumed fairly smooth channels with a surface roughness of $ϵ≈$ 1.5 μm (drawn Cu tubing).

where $Nu$ is the Nusselt number based on the average hydraulic diameter (⁠$d$⁠), $Pr$ is the fluid Prandtl Number, and $kf$ is the fluid thermal conductivity. Convective heat transfer in the fluid cannot be modeled with either constant heat flux or constant wall temperature boundary conditions. Hence, an average Nusselt number for flow in a circular channel with constant wall heat flux and constant wall temperature was chosen for the laminar flow regime.

where $q″$ is the heat flux, $hamb$ is the heat transfer coefficients to ambient at corner elements, $Al$⁠, $Atop$ are the lateral and top surface areas, respectively. The block lateral surface area is assumed to be equal for all faces of the block in contact with neighboring elements or ambient. The bottom surface is assumed to be adiabatic. Here, $Tnb$ and $Tamb$ denote the neighboring element top surface and ambient temperatures. Effects of heat spreading (lateral block to block heat transfer) within the solid domain are considered within the heat conduction equation. Equations (5) and (10) result in a set of linear equations which are solved to obtain the fluid and heated surface temperature profiles for the entire cold plate.

where $min$ represents the multivariable minimization function over the entire cold plate layout. Here, the maximum top surface temperature among all element blocks $Ts,max$⁠, total pressure drop within the entire channel $ΔPtot$ and the standard deviation of the top surface temperatures $Ts,dev$ are minimized with respect to the flow layout using the minimization function. The inlet mass flow rate $m˙$ (Eq. (11)) and local heat flux $qloss″$ at individual element locations (⁠$i, j$⁠, Eq. (12)) form the optimization constraints. The multi-objective optimization of $Ts,max$⁠, $ΔP$⁠, and $Ts,dev$ is solved using three unique approaches described herein.

A genetic algorithm (GA) is a type of evolutionary search algorithm [51]. The GA is effective for discrete optimization problems due to its heuristic search method which does not depend on gradient calculation [52]. The matlab-based ga function was applied for finding the optimal discrete design variables. We started by implementing the ga function to find the optimal design variables which correlated to the discrete flow block codes for all element locations. Each design variable represented an integer code for unique flow block type and flow direction at different grid locations.

We assigned an arbitrarily high value to the performance parameters being minimized for incomplete designs with partially connected flow layouts (Fig. 2(c)). This direct implementation of the ga function for searching the optimal layout was unable to find a solution due to the limited feasible flow layouts in a large design space within a reasonable combination of the population size and number of iterations. We developed a modified objective function (Fig. 2(a)) which incorporates encoding of designs in terms of continuous variables and penalizing incomplete designs to address this.

To express the design in terms of continuous variables, each flow block is modeled in terms of four variables incorporating two pairs of coordinates in a 2D number line (Fig. 2(b)). The inlet and outlet flow direction for the block is determined based on the quadrant of the two variable pairs. This encoding approach correlates similar flow blocks which are not possible using the discrete representation. While obtaining the performance parameters for incomplete designs, penalty values are assigned based on different design characteristics (Fig. 2(c)). Negative penalties are given for rewarding continuous chains of flow while positive penalties are given for having flow discontinuities as well as for having inlet/outlet ports at different grid elements from the preset locations. These modifications make the objective function behavior smoother in the design space and help identify optimal designs. The modified GA implementation was able to successfully find optimal layouts for a smaller grid size of 3 × 4 blocks with approximately 300 iterations of the solver with a population size of 1,000. The same approach however was unable to find a complete design solution for moderate grid sizes of 5 × 5, highlighting the main limitation of the GA-based layout optimization approach, thus inspiring the search for alternate methods.

We were unable to perform layout optimization using the ga function because a majority of variable combinations led to incomplete or nonfeasible flow layouts. Given the limited number of feasible design layouts possible for a large design space, we considered an approach that relies on exhaustively searching for optimal performance within a set of feasible layouts (Fig. 3(a)). We used path search algorithms [53] for finding feasible fluid flow paths connecting predefined inlet and outlet locations. In the process of finding all paths, we created a fringe for storing design nodes, defined as vectors with flow block codes at each element location. At each element location starting from the inlet, we added to the fringe nodes by finding all possible flow blocks that can be assigned (Fig. 3(b)) based on the previous block. A node was terminated if the block assignment led to an already occupied location on the grid or found a boundary element other than the outlet location. Finally, a node leading to the outlet location was added to the set of feasible designs. The set of feasible designs was stored locally and is called for calculating performance given the boundary conditions for the layout optimization problem.

Using the exhaustive design space search leads to optimization which ensures the best possible design. Even though the number of possible layouts is limited, the design space grows exponentially with the number of grid elements (Fig. 3(c)). The requirements on runtime memory for maintaining the fringe with all design nodes and local memory for storing all possible layouts also increase significantly due to this exponential growth. For instance, the number of possible designs for a 3 × 4 grid size is 38 while that for a 6 × 6 grid size is 1,246,850. The time required for calculating the performance metrics for each unique heat dissipation profile for all possible designs is another detrimental factor against using this approach.

The previous two approaches were incapable of meeting the requirements for rapid layout optimization. An approach building on the merits of the last two approaches is discussed herein. The advantages offered by the heuristic search nature of the discrete optimization methods and ease of finding feasible layouts using the path search method were combined in this approach (Fig. 4(a)).

Figure 4(b) shows the flow layout formed by the shortest path starting from the inlet, connecting a series of control points leading to the outlet. We used a heuristic search method [53] to find the shortest path on a partially filled grid using the Manhattan distance between the current location and the next control point location as the heuristic. The process of connecting the shortest path through control points is repeated for a randomized population of control points. Performance parameters are calculated and compared for each of the generated designs. Progressively increasing the number of control points on the grid and randomized population search increased the possibility of finding the optimal solution using this approach for any given combination of boundary conditions. A smaller number of control points usually led to shorter designs which are desirable for having low-pressure drop. Conversely, a greater number of control points ensures that the fluid reaches all heat dissipation locations with favorable thermal performance. Several designs analyzed can be used to find not only a single optimal layout but a series of Pareto optimal layouts for the multi-objective optimization problem.

The OPS approach successfully identified the Pareto front for the bi-objective optimization of thermal resistance and pressure drop shown in Fig. 5. Optimization was carried out for grid sizes of 4 × 4, 5 × 5, and 6 × 6 for a design having the same overall dimensions with equivalent boundary conditions. The Pareto fronts were obtained using a number of control points ranging from 2 to 24 (depending on grid size) each with 500 randomized sets of control points. As seen from the results (Fig. 5), for the same design, higher grid size results in smaller block dimensions and thus overall higher pressure drop. Thermal performance improves with increasing grid size with minimum temperatures of 79.8 °C for 4 × 4, 64.6 °C for 5 × 5 and 54.5 °C for the 6 × 6 grid size. As the grid size increases, the maximum T_s,max which corresponds to the shortest path connecting the inlet and outlet port locations (simple path connecting the first row of elements on the coarse grid) increases. The increase is due to the lower portion of the total area covered by the fluid path for higher grid sizes. The overall minimum pressure drop which corresponds to these designs also decreases showing that the Pareto front improves with respect to both design variables individually with an increase in the grid size.

Figure 6 shows the effects of different boundary conditions and design constraints on optimal layouts obtained using the OPS method. For a grid size of 4 × 4, the effect of two heat dissipation profiles including uniform and concentrated profiles was considered (Figs. 6(a) and 6(b)). For the uniform loss profile, all grid elements were assigned a fixed heat dissipation value of 200 W (⁠$q″$= 50 W/cm²). In the concentrated loss case, the four central grid elements were assigned a fixed heat dissipation value of 800 W (⁠$q″$= 200 W/cm²), with all other elements being adiabatic (⁠$q″$= 0) ensuring that the total heat dissipation in both cases is identical. The pressure drop for both uniform and concentrated heat dissipation cases is constrained at a value of $ΔP$ = 0.81 kPa for an inlet mass flow rate of 0.1 kg/s with water as coolant. Aluminum is considered a solid domain material. The fluid block height is $L$ = 2 cm and diameter values for all blocks are $d$ = 1.5 cm. The Pareto-optimal layouts identified by the OPS differ for the two loss profiles. The uniform heat dissipation profile design connects dispersed points on the grid while the concentrated loss profile design covers all of the central location. Moreover, the maximum top surface temperature among all elements with uniform loss is 130.2 °C with design 1 which is originally optimized for uniform heat dissipation (Fig. 6(a)) and 166.3 °C with design 2 which is optimized for concentrated heat dissipation (Fig. 6(b)). The maximum surface temperature in the case of a concentrated loss with design 1 is 194.2 °C while that for design 2 is 151.2 °C. Both designs perform better with the heat dissipation profile for which they are optimized. Allowing the inlet and outlet elements to change positions can lead to further improvement in the optimal designs (Figs. 6(c) and 6(d)). For this, we carried out the design optimization analysis for a grid size of 5 × 5 having flow conditions and individual block dimensions the same as before. The Pareto front shifts to the left (lower maximum surface temperature) when the inlet/outlet locations are allowed to float and are not fixed at the bottom side of the cold plate, indicating that better Pareto optimal performance is possible. Fixed inlet/outlet port locations serve as additional constraints, the removal of which acts to increase the number of possible layouts which can generate better thermal-hydraulic performance.

We carried out optimization of the cross section geometrical parameters after obtaining the optimal flow layout for given grid size and boundary condition. The optimization formulation for the cross section geometry is similar to Eq. (5) with diameter values for the flow blocks now representing design variables instead of the 2D flow layout. We note that the geometry optimization described herein is not only a part of the overall design optimization but can also independently serve as a design improvement tool for maximizing performance for a set 2D layout.

The objectives and constraints for geometry optimization are continuous functions of the design variables (diameters, $d$⁠). We can use gradient-based optimization which converges faster to local optima as compared to methods such as the GA. Here, we use the matlabfmincon function for geometry optimization. Figure 7 highlights the results of geometry optimization for a given flow layout in different cases. Figure 7(a) shows a set of Pareto optimal points for thermal-hydraulic bi-objective optimization obtained using optimization of the channel diameters. The insets in the plot show optimal designs for different combinations of thermal resistance and pressure drop objectives. On the pareto front, designs with higher weightage to pressure drop objectives have larger channel diameters, because pressure drop increases for lower channel diameters. Designs with higher weightage to thermal resistance objective have lower diameters since the heat transfer coefficient is inversely proportional to the hydraulic diameter.

We studied the effect of heat dissipation profile on optimal diameters with uniform and concentrated profiles. Once again, the uniform heat dissipation profile assumes equal heat flux on all element faces while the concentrated profile assumes heat flux only on the four central elements, with an equal total loss in both cases. In the case of uniform heat dissipation (Fig. 7(b)), the optimal diameter values are nearly equal for all element faces, resulting in low pressure drop encountered for the gradual changes in diameters. For the concentrated heat dissipation profile (Fig. 7(c)), the diameter values in the central region where heat dissipations are applied are smaller than the uniform heat dissipation case. The smaller diameters lead to increased local heat transfer coefficients and also cause higher pressure drop. To counter the elevated pressure drop while maintaining good thermal performance, the optimal design selects larger diameter values in the outer regions of the grid near the edges, with the resulting total pressure drop being similar for both the uniform and concentrated heat dissipation cases.

Figure 8 shows a design optimization framework that integrates layout optimization and cross section geometry optimization. The input for the design optimization problem is obtained as a 2D schematic of the electronic system with locations specified for heat generating devices. The layout is then transformed into a coarse grid, with the device locations translated to the grid element with specified heat dissipation. First, the optimal layout search is carried out using OPS, using the heat transfer and flow boundary conditions. A design is chosen from the Pareto-optimal points to get the desired tradeoff between thermal and hydraulic objectives. Second, geometry optimization is performed on the flow layout obtained in the previous step. For the geometry optimization step, the pressure drop with the design obtained from layout optimization is used as a constraint to minimize the thermal objective.

Figure 9 shows a case study applying the design optimization method for single-phase liquid cooling of electronics. The study targets an XM3 series half-bridge module developed by Wolfspeed (Fig. 9(a)) [54,55]. The XM3 series power modules are developed for electric vehicle (EV) chargers and traction drives. The module packaging includes discrete silicon carbide (SiC) semiconductor devices which act as the major heat generating sources with a silicon nitride (SiN) substrate and a copper (Cu) base plate for heat spreading. A region within the geometry of the module is divided into a coarse grid (Fig. 9(b)) of 8 × 4 blocks for design optimization.

First, we obtained the optimal layout for the given set of input parameters using the OPS approach (Fig. 9(c)) to minimize thermal resistance. Second, the cross section geometry optimization was used to generate the diameter values to further minimize thermal resistance with constrained pressure drop (Fig. 9(d)). We built 3D computer aided design models for the layout (Fig. 9(f)) and geometry optimized (Fig. 9(g)) designs using autodeskfusion360. Finally, we mapped the performance of the two designs using 3D computational fluid dynamics (CFD) simulations in comsolmultiphysics (version 5.6). We performed simulations for a total heat dissipation of 500 W distributed among the devices (Fig. 9(e)). This heat loss corresponds to a local device-level heat flux of 178.5 W/cm² as computed using the device base area. Considering the effect of heat spreading from device to base plate, the heat flux was reduced to an average of 62.5 W/cm² due to the larger area of the base plate. We considered water–ethylene glycol mixture (50% by volume) as the coolant for total flow rates ranging from 2 to 8 liters per minute with a bulk coolant inlet temperature of 25 °C. water and ethylene glycol (WEG) was chosen as a preferred coolant due to its use in EV and hybrid electric vehicle (HEV) onboard cooling systems due to its lower freezing points when compared to water. We chose WEG50 as our coolant for analysis in order to recreate conditions similar to the ones used in reported state-of-the-art (SOA) cold plate performance metrics which are available for comparison. The cold plate material was assumed to be aluminum alloy. A thermal interface material layer was included between the heatsink and module baseplate having a 100 μm thickness and a 5 W/(m·K) thermal conductivity.

We carried out preliminary simulations for multiple mesh densities including 0.71, 3.04, 9, and 16.01 × 10⁶ elements to check for mesh independence. Results showed <2% change in both maximum temperature and pressure drop for the three higher mesh densities. Based on this result, we chose a mesh density corresponding to 3.04 million elements considering both simulation accuracy and computational expense for further simulations. We obtained the junction-to-coolant thermal resistance which is the ratio of the difference between the maximum SiC device junction and coolant temperatures to total heat dissipation along with total cold plate pressure drop (Fig. 9(h)). The design with both optimum layout and geometry outperforms the design having only the optimum layout with an average 4.5% reduction in the junction-to-coolant thermal resistance for a given pressure drop over the range of flow rates considered. Figure 9(h) shows a comparison between the optimal cold plate and an SOA cold plate developed by Wieland Microcool (CP4012) specifically for cooling of XM3 power modules. The SOA cold plate has higher thermal resistance (0.2 °C/W) [55] for a comparable value of pressure drop when compared to the optimal design (0.15 °C/W), demonstrating a 25% reduction at approximately 16 kPa pressure drop.

We developed a rapid design optimization tool for single phase liquid cooled heat sinks such as cold plates and microchannel coolers used in electronics cooling. Using our method, optimized solutions were typically obtained within a time scale of 60–120 s depending on the design complexity and grid size. Here, we used serial computing for performing all computations with the aim of developing a simple software tool for reduced-order optimization. Parallel computing can be employed in future work for several functions including layout search on randomized populations [56], which will further reduce the time required for finding optimal designs. We used a modeling approach that predicts thermal-hydraulic performance on a coarse grid resulting in low computational effort. Although the key assumptions made to develop our simplified modeling approach can result in a lower-fidelity estimation of thermal-hydraulic parameters, the comparative performance of designs under investigation can be expected to be the same as that obtained using high-fidelity simulation methods. This enables the use of reduced-order modeling for design optimization. Despite the nonoptimal prediction fidelity, modeling accuracy can be further improved with the incorporation of correlations developed for flow and temperature parameters using high fidelity CFD simulations.

Flow layout can have a significant impact on cooling performance [28–34]. A layout optimization submodule was developed to find flow route extracting maximum thermal performance while maintaining low-pressure drop. Discrete search and optimization techniques were explored to develop a method for finding the flow layout that optimizes the thermal-hydraulic performance. An approach relying on the advantages of the path-finding algorithm and a population-based search was realized in the form of the OPS approach. Combining the heuristic nature of GA with the assurance of finding a complete layout with path search algorithms, our OPS approach enables a rapid layout search. Even though the OPS approach is able to capture layouts performing better than alternatives for most cases, it does not guarantee an optimal solution as evident from the missed Pareto optimal points while generating the Pareto front. The chances of finding an optimum could be improved by increasing the number of search points (granularity of the fundamental design elements) at the cost of increased time requirement to conduct the search and estimate thermal-hydraulic performance.

Our work focuses on designs with nonsplitting flow layouts (serial channels). Dedicated correlations to differentiate the flow behavior and resulting thermal performance changes in elbows as compared to straight blocks are required. Furthermore, additional development is needed to establish a performance modeling approach that can include flow branching. In addition, the method should be modified to accommodate nonuniform grids to allow local refinement within the cold plate geometry. Layout search and geometry optimization methods also need to be modified to consider parallel flow architectures. Similar to the simple flow blocks used in this work, correlations can be developed for blocks including performance enhancing aspects such as internal fins and turbulence inducing features. Advanced manufacturing methods including additive manufacturing can make such designs realizable [57–62]. In the future, it would be interesting to utilize our tool to conduct additional case studies for practical cooling applications with differing heat dissipation profiles and compare the results with the well-established cooling methods currently used [63]. Furthermore, future work is needed to obtain rigorous design optimization performance benchmarks against state-of-the-art design methods and design optimization software packages (ansys, comsol) in the context of required computational resources and timescale for the identification of an optimal solution.

Two-phase heat sinks for cooling electronics are a topic of high interest. Flow boiling in channel flow offers potential for achieving higher heat flux removal capabilities when compared to single-phase liquid flow while also providing isothermal cooling [64–66]. The coarse meshing approach for performance prediction presented in this work can be translated for heat sinks undergoing two-phase flow without impacting the computational speeds. This can be achieved with reliable correlations which are available for predicting macroscale flow and thermal behavior. The reliability of such applications of the tool needs to be validated and improved with the help of high-fidelity simulations [67–69] and experimental approaches [70–72].

Rapid design optimization creates the potential for integration with electrothermal codesign tools already available [2–5]. Traditional approaches mainly rely on detailed thermal analysis based on numerical methods requiring significant time and computational effort. These approaches limit the number of design iterations available for identifying the optimal packaging of electronics. These also require a level of user expertise in setting up the numerical simulations to ensure accurate results. The methods demonstrated in this work allow for the rapid design of compact electronics packaging systems with high performance cooling solutions. The designs generated using this method can also be used as starting points for high fidelity design optimization using TO.

We demonstrate that the design of single-phase liquid cooled channel heat sinks can be optimized using limited time and computing resources. Our approach combines modeling on a coarse grid using existing performance correlations for obtaining pressure drop and temperature profiles, flow layout optimization through multiple path searches using randomized populations and gradient-based optimization for cross section geometries. The reduced-order nature of this method offers significant benefits in terms of computational and temporal requirements when compared to traditional design optimization methods using CFD simulations coupled with numerical optimization. We show the ability of the developed method to address heterogeneous loss profiles commonly seen in electronics packaging. We conducted a case study demonstrating the application of the developed design optimization method for a commercial XM3 power electronics module. The results indicate a potential reduction of 25% in junction-to-coolant thermal resistance compared to an aggressive state-of-the-art cooling solution. The design tool enables the versatile incorporation of enhancements such as consideration of parallel flow designs, two-phase flows, internal flow features (fins), and nonuniform grids. The rapid solution and low computational requirement of the developed tool enable potential design application in electrothermal codesign for electronics packaging.

The authors thank Professor James Allison, Dr. Satya Ravi Teja Peddada, Cary Laird for their support with the channel area optimization model. The authors also thank Professors Mehdi Asheghi and Kenneth Goodson of Stanford University for useful discussions regarding thermal-hydraulic considerations and modeling. The authors also gratefully acknowledge funding for this work in part from the Power Optimization of Electro-Thermal Systems (POETS) National Science Foundation Engineering Research Center with Cooperative Agreement No. EEC-1449548, and the International Institute for Carbon Neutral Energy Research (WPI-I2CNER), sponsored by the Japanese Ministry of Education, Culture, Sports, Science, and Technology.

• POETS: National Science Foundation Engineering Research Center (Cooperative Agreement No. EEC-1449548; Funder ID: 10.13039/100000149).

• I2CNER: Japanese Ministry of Education, Culture, Sports, Science and Technology (Funder ID: 10.13039/501100007068).

$A$ =

area (m²)

$C$ =

specific heat (kJ kg⁻¹ K⁻¹)

CFD =

computational fluid dynamics

$d$ =

channel hydraulic diameter (m)

EV/HEV =

electric vehicles/hybrid electric vehicles

$f$ =

friction factor

GA =

genetic algorithm

$h$ =

heat transfer coefficient (W m⁻² K⁻¹)

$k$ =

thermal conductivity (W m⁻¹ K⁻¹)

$L$ =

fluid block size/length (m)

$m˙$ =

mass flow rate (kg s⁻¹)

$Nu$ =

Nusselt number

OPS =

optimal path search

$P$ =

pressure (Pa)

$Pr$ =

Prandtl number

$q”$ =

heat flux (W m⁻²)

$Re$ =

Reynolds number

$S$ =

conduction shape factor (m)

SOA =

state-of-the-art

$T$ =

temperature (⁠$°C$⁠)

TO =

topology optimization

$v$ =

average flow velocity (m s⁻¹)

WEG50 =

water and ethylene glycol mixture, 50% by volume

Greek Symbols

$ΔP$ =

pressure drop (Pa)

$ϵ$ =

roughness (m)

$Θ$ =

thermal resistance (⁠$°C$ W⁻¹)

$μ$ =

dynamic viscosity (N s m⁻²)

$ρ$ =

density (kg m⁻³)

Subscripts

$c$ =

cross-section

$eff$ =

effective

$f$ =

fluid

$in$ =

inlet

$l$ =

lateral

$out$ =

outlet

$s$ =

solid

$s,dev$ =

surface standard deviation

$s,max$ =

surface maximum

$top$ =

top surface