Abstract

Heterogeneous and complex electronic packages may require unique thermomechanical structures to provide optimal heat guiding. In particular, when a heat source and a heat sink are not aligned and do not allow a direct path, conventional thermal management methods providing uniform heat dissipation may not be appropriate. Here we present a topology optimization method to find thermally conductive and mechanically stable structures for optimal heat guiding under various heat source-sink arrangements. To exploit the capabilities, we consider complex heat guiding scenarios and three-dimensional (3D) serpentine structures to carry the heat with corner angles ranging from 30 deg to 90 deg. While the thermal objective function is defined to minimize the temperature gradient, the mechanical objective function is defined to maximize the stiffness with a volume constraint. Our simulations show that the optimized structures can have a thermal resistance of less than 32% and stiffness greater than 43% compared to reference structures with no topology optimization at an identical volume fraction. The significant difference in thermal resistance is attributed to a thermally dead volume near the sharp corners. As a proof-of-concept experiment, we have created 3D heat guiding structures using a selective laser melting technique and characterized their thermal properties using an infrared thermography technique. The experiment shows the thermal resistance of the thermally optimized structure is 29% less than that of the reference structure. These results present the unique capabilities of topology optimization and 3D manufacturing in enabling optimal heat guiding for heterogeneous systems and advancing the state-of-the-art in electronics packaging.

1 Introduction

Enhancing computing capabilities at the device and die levels while increasing the speed, efficiency, and reliability of electronics makes thermal management ever more challenging. Conventional cooling solutions based on higher thermal conductivity materials or heat exchangers dissipate the heat from a source to a sink only in a unidirectional manner, and they cannot guide the heat flow when the heat source and heat sink are unaligned. The earlier studies on guiding the heat flux have investigated the concept via ballistic phonon transport in holey silicon [1], through silicon via integrated thermoelectric cooling [2], thermal cloaks [35], and copper-PDMS based thermal shifters [6]. However, the existing macroscale approaches for thermal management in electronic packaging still rely on creating thermally conductive channels from a heat source to the heat sink using the through-hole vias [7], creating a cavity in the printed circuit board, and using high thermal conductivity inlays [8]. Compared to two-dimensional (2D) packaging scenarios, investigations on thermal management strategies for three-dimensional (3D) packages are still limited. For instance, common methods for resolving the thermal challenges in 3D ICs are restricted to floor plan optimization, reducing heat sink thermal resistance, lowering the power consumption of the chip, and applying thermal vias [9]. These conventional solutions are effective when a heat-generating component is aligned with a heat sink (Fig. 1(a)). However, with a limited routing space in a 3D package, alternative solutions may be required for optimizing the thermal routing of unaligned heat source and heat sink components (Fig. 1(b)). Therefore, novel heat guiding structures are required for acquiring optimal routing of heat flow under nonuniform power distribution to achieve the desired thermal objectives.

Fig. 1
(a) Heat transfer path in a conventional packaging scenario where heat-generating component and heat sink are aligned. In this case, components with different temperature tolerances can easily be impacted if positioned in the heat transfer path. (b) Heat transfer path in a nonconventional packaging scenario where a heat sink may not be aligned with a heat source. In the presence of several components with different temperature tolerances in the system, novel heat guiding structures with zigzag or S shape designs would be favorable for optimal routing of the heat flow from the heat source to the heat sink. These heat guiding structures may create an efficient heat transfer path while minimizing the thermal crosstalk with the components of different temperature tolerances. (c) Thermomechanical optimization process for complex packaging scenarios. The thermal boundary conditions include a heat flux (q″) on the heat source side while fixing the temperature on the heat sink side (Tsink). The mechanical boundary conditions include a force applied on the heat sink side and a fixed end on the heat source side.
Fig. 1
(a) Heat transfer path in a conventional packaging scenario where heat-generating component and heat sink are aligned. In this case, components with different temperature tolerances can easily be impacted if positioned in the heat transfer path. (b) Heat transfer path in a nonconventional packaging scenario where a heat sink may not be aligned with a heat source. In the presence of several components with different temperature tolerances in the system, novel heat guiding structures with zigzag or S shape designs would be favorable for optimal routing of the heat flow from the heat source to the heat sink. These heat guiding structures may create an efficient heat transfer path while minimizing the thermal crosstalk with the components of different temperature tolerances. (c) Thermomechanical optimization process for complex packaging scenarios. The thermal boundary conditions include a heat flux (q″) on the heat source side while fixing the temperature on the heat sink side (Tsink). The mechanical boundary conditions include a force applied on the heat sink side and a fixed end on the heat source side.
Close modal

Capabilities of heat guiding structures to control the heat flux effectively [4,6,1012] and to guide the heat in unaligned configurations of heat source and heat sink can lead to optimal solutions for complex packaging scenarios [13,14]. For example, 3D metamaterial solutions for heat guiding structures can potentially protect components with different temperature limits while effectively dissipating the heat to address the demand for high integration density [15,16]. There have been fundamental studies on the possibility of developing thermal heat guiding structures and manipulating heat flows in 3D structures [5]. Dede et al. established optimized 3D heat routing structures for power electronics gate drive printed circuit board (PCB) thermal management [17]. While most of these designs accommodated thermal management challenges in electronic packaging, further investigation is still required for the thermomechanical performance of the heat guiding structures in macroscale power semiconductor devices. The introduction of the topology optimization technique created a world of possibilities for the optimal design of heat guiding structures under desired objectives. Existing studies have extensively investigated topology optimization with mechanical boundary conditions for maximizing the mechanical stiffness [1824] and with thermal boundary conditions for maximizing the temperature diffusivity. Iga et al. considered thermal conduction and convection boundary loading in the topology optimization of structural designs [25]. Zhang et al. developed a numerical model of topology optimization for isotropic and anisotropic structures to achieve the least dissipation of heat transport potential capacity [26]. Lundgaard et al. presented a density based topology optimization method for thermal energy systems by coupling fluid and heat transfer models [27]. Zhou et al. applied a design-dependent convection in conjunction with conduction topology optimization for industrial applications [28]. There are several topology optimization examples for designing metamaterials in electronic systems such as topology optimization of an actively cooled system integration into downhole electronics [29], thermal-composite design optimization for thermal management of PCB electronics [13], design of a power semiconductor module using topology optimization for efficient cooling [30], topology optimization for effective heat removal in 3D packages [31], and experimental demonstration of a topology optimized heat sink for a tablet [32]. However, further investigations for coupled thermomechanical constraints are still required for electronic packaging applications.

For combining the objectives of maximizing the mechanical stiffness while minimizing the temperature gradient, thermomechanical topology optimization can be employed [33]. Dede [34] performed thermomechanical topology optimization to decrease component weight and increase gravimetric power density in the thermal bracket of a magnetic inductor. Takezawa et al. [35] investigated thermomechanical topology optimization in arbitrary 2D shapes considering the thermal expansion effect, and Zhu et al. [36] studied the temperature-constrained topology optimization for thermoelastic structures under a design-dependent temperature field. However, thermomechanical topology optimization in earlier studies has been mostly focused on using multiple materials [37,38], confined to 2D domains [3943], or challenging to manufacture [44].

In this study, we developed a finite element method (FEM) based on thermomechanical topology optimization for designing novel 3D homogenous heat guiding structures with optimized thermal and mechanical properties under 50% volume fraction constraint. Figure 1(c) demonstrates developing thermomechanically optimized heat guiding structures for an unaligned configuration of a heat source and a heat sink. Using thermal, mechanical, and thermomechanical objective functions results in unique topologies that satisfy the applied boundary conditions. A minimum dimension of 30 μm constraint has been applied to our design space for securing the manufacturability of 3D metallic thermomechanically optimized heat guiding structures via SLM 125 HL metal printer.

2 Topology Optimization of Three-Dimensional Heat-Guiding Structures

2.1 Density-Based Topology Optimization Algorithm Development.

Topology optimization-based finite element methods have been explored in the past to enable heat flow control in arbitrary (e.g., noncircular or nonspherical) geometries [45,46] and bifunctional cloaking [47]. In this paper, we utilized the 2D thermal conductivity design approach proposed by Dede [46] for heat flux shielding and expanded the approach to a 3D design space domain. The heat flux shielding is a suitable choice for minimizing the magnitude of the temperature gradient and heat flow within a selected region of a system. This approach allows us to achieve the optimized topology for a given volume fraction while minimizing thermal compliance. For mechanical optimization, a density model with regularization via the Helmholtz equation was used [48]. The optimization function that we solve in this study is as follows:
(1)

where the objective function J(ρ) is defined in terms of mechanical Jm and thermal Jθ objectives, where 0ω1 is the weight factor, and ω=0.5 was used for our calculations. gj,j=1,,Nc is the Jth volume constraint to subregions of the design domain (εj) [49]. νe and ρe are the area and density of element e, respectively; and νj is the upper limit for the volume constraint j. We use the equilibrium equation for the steady-state linear elastic deformation for defining the mechanical objective function—Jm—where K and u are the stiffness matrix and displacement vector, respectively. In addition, we consider steady-state conduction in a 3D domain as the thermal objective function—Jθ—where κ is the thermal conductivity matrix, and θ is the nodal temperature vector.

Sufficient penalization of the intermediate densities—ρe—by the interpolation function is required for the convergence of the solution. Solid isotropic material with penalization (SIMP) approach has been implemented in COMSOL for the topology optimization of 3D heat guiding structures. The penalization factor (p)can be set to values greater than 1 for penalizing the intermediate densities [29]; we used the value of p=3for thermal, mechanical, and thermomechanical topology optimization simulations [48]. To utilize gradient-based optimization techniques, the problem is solved by allowing intermediate values of 0ρe1 for the design parameter to exist. The design parameters are bounded from below with ρmin with the value of 0.001; therefore 0<ρminρe1 and we can derive the penalized density (ρp) as below:
(2)
We initialize the optimization process by providing an initial design space and utilize the material properties of SLM manufactured SS316 L for the calculations with Young's modulus of ∼ 150 (GPa) [50] and thermal conductivity of ∼15 (Wm−1K−1) [51]. Evaluating the sensitivity of the overall objective function with respect to the control variable ρe can be rephrased as calculating the derivative of the objective function (Jρe) which can be derived from the sensitivity analysis summarized as below [52]:
(3)

Based on the sensitivity input, the SIMP optimization approach updates the design regarding the Method of Moving Asymptotes (MMA) [53], which was initially developed for structural optimization of mechanical loads.

2.2 Thermomechanical Topology Optimization of Heat Guiding Structures.

Defining an effective thermal path design is crucial for transferring the excessive heat to the heat sink, protecting the components with different thermal tolerances from hot surroundings, and preventing heat leakage. However, the limited available space between the heat source and the heat sink makes this task challenging. Moreover, the package should be designed such that the stress in any component does not exceed the permissible extent of that material and can withstand the thermomechanical cyclic loading imposed during the system operation. Knowing that the maximum volume fraction of the utilized metal in the PCB board of packaging structures is about 50%, we used this number as a limiting design parameter. Hence, we solved the density-based topology optimization problem for 3D heat guiding structures with an objective volume fraction of 50%. Figure 2 demonstrates a series of investigated designs for possible 3D heat guiding structures where the thermal path from the heat source to the heat sink deviates from conventional packaging scenarios. The initial structures have been chosen to demonstrate the impact of varying distances from the heat source to the heat sink in nonconventional thermal management cases. With the additive manufacturing capability, multiple reference structure choices were available for evaluating the performance of the optimized heat guiding structures. Here, two reference structure designs with an identical topology of the initial structure but with a 50% volume fraction were studied. The solid reference structure mimics its initial structure but with 50% of the cross-sectional area. The infill reference structure has an identical geometry as the initial structure but with a rectilinear-infill pattern that represents 50% porosity. The rectilinear-infill pattern was chosen due to its common usage in additive manufacturing [5456] and its opportunities for complex numerical studies.

Fig. 2
Boundary conditions and the resultant thermal, mechanical, and thermomechanical topology optimization on the following designs. (a) Initial structure with parameters of ⊝ = 30, L = 10 mm, W = 5 mm, t = 2 mm, and D = 4 mm. (b) Initial structure with parameters of ⊝ = 45, L = 10 mm, W = 7 mm, t = 2 mm, and D = 4 mm. (c) Initial structure with parameters of ⊝ = 60, L = 10 mm, W = 12 mm, t = 2 mm, and D = 4 mm. (d) Initial structure with parameters of L = 10 mm, W = 10 mm, t = 2 mm, and D = 4 mm. e) Initial structure with parameters of L1 = 30 mm, L2 = 12 mm, W = 15 mm, t1 = 4 mm, t2 = 2.5 mm, and D = 15 mm. The solid reference structures for designs (a)–(d) have identical dimensions as their representing initial structures but 2 mm thickness in the y-direction. The solid reference structure for design (e) has identical dimensions as its initial structure but with W/√2 = D/√2 = 10.6 mm. The infill reference structure for all designs (a)–(e) is identical to their initial design but with a rectilinear infill pattern that creates a 50% volume fraction.
Fig. 2
Boundary conditions and the resultant thermal, mechanical, and thermomechanical topology optimization on the following designs. (a) Initial structure with parameters of ⊝ = 30, L = 10 mm, W = 5 mm, t = 2 mm, and D = 4 mm. (b) Initial structure with parameters of ⊝ = 45, L = 10 mm, W = 7 mm, t = 2 mm, and D = 4 mm. (c) Initial structure with parameters of ⊝ = 60, L = 10 mm, W = 12 mm, t = 2 mm, and D = 4 mm. (d) Initial structure with parameters of L = 10 mm, W = 10 mm, t = 2 mm, and D = 4 mm. e) Initial structure with parameters of L1 = 30 mm, L2 = 12 mm, W = 15 mm, t1 = 4 mm, t2 = 2.5 mm, and D = 15 mm. The solid reference structures for designs (a)–(d) have identical dimensions as their representing initial structures but 2 mm thickness in the y-direction. The solid reference structure for design (e) has identical dimensions as its initial structure but with W/√2 = D/√2 = 10.6 mm. The infill reference structure for all designs (a)–(e) is identical to their initial design but with a rectilinear infill pattern that creates a 50% volume fraction.
Close modal

The initial design domains with 100% volume fraction provide the smallest thermal resistance for the heat transfer. Therefore, we normalize the thermal resistance of both the optimized topology and the reference structures to the thermal resistance of their initial structure. The topology optimization algorithm creates a path between the heat source and the heat sink with the smallest thermal resistance while preserving the requirements of the density-based objective function. A lower thermally resistive path guarantees a lower maximum temperature at the heat source.

By combining thermal and mechanical requirements, we used Comsol's optimization module and the Globally Convergent of the Method of Moving Asymptotes (GCMMA) [57] as the optimization method with thermal and mechanical boundary conditions. The method of moving asymptotes was chosen for two main reasons. First, the method was proven to work very well on a large variety of topology optimization problems [58]. Second, the method is very popular for parallel computations because of the separable nature of approximations. The thermal boundary conditions consist of applying a heat flux with a representative value of 250 Wcm−2 [59] on the heat source side and a constant temperature of 27 °C on the heat sink side. The mechanical boundary conditions include a 100 N compression force on the heat sink side (conversion of the maximum stress around through-silicon vias on the device layer [60]) and a fixed end on the heat source side. At each structure, the boundary conditions are applied to a small window of 2 mm × 2 mm placed at the beam ends. The boundary conditions for the investigated design domains were applied to an initial structure with a 100% volume fraction (Figs. 2(a)–2(e)).

2.3 Thermomechanical Optimization Results.

Thermal, mechanical, and thermomechanical topology optimization results are shown in Figs. 2(a)2(e) for the investigated designs. The solid reference structure of each design for Figs. 2(a)2(d) mimics the initial structure but with 2 mm thickness in the y-direction. For Fig. 2(e), W and D values in the solid reference structure were reduced to 10.6 mm, and other parameters remained as demonstrated in the initial structure to reach 50% volume fraction in the reference structure. The infill reference structure has the same thickness as the initial structure in the y-direction (D), but the implemented porosity is responsible for creating the 50% volume fraction. By applying the thermal boundary conditions demonstrated for each design in Figs. 2(a)2(e) and examining the temperature gradient across the structure, the thermal resistance for each of the designs was calculated using Rstructure=ΔTq. For evaluating the thermal performance of each design, thermal resistance was also calculated in their corresponding reference structures. Table 1 summarizes the normalized thermal resistance (Rstructure/Rinitial) for both thermal and thermomechanical optimization results of each investigated design at which Rinitial is the thermal resistance of the initial structure. The initial structures have the largest cross-sectional area (smallest thermal resistance) when compared to the optimized and the reference structures; therefore, normalizing the thermal resistance values of the optimized and reference structures to the thermal resistance of the initial structure results in values greater than 1. Since the functions that were applied to the optimization process do not result in the same topology, the thermomechanically optimized structures have inherent tradeoffs in both thermal and mechanical performances. This tradeoff is demonstrated as a higher normalized thermal resistance in thermomechanically optimized structures when compared to the thermally optimized structures for all the studied designs. As simulation results demonstrate, the normalized thermal resistance ratios of the infill reference structures are ∼50% greater in all designs except for the serpentine design (e) with a ∼33% increase when compared to the solid reference structures. This increase in thermal resistance is due to a smaller cross-sectional area of the infill reference structure—2.6 mm2 in designs (a)–(d) and 86.8 mm2 in the serpentine design (e)—compared to the solid reference structure—4 mm2 in designs (a)–(d) and 112.4 mm2 in the serpentine design (e). The solid reference structure is selected for future comparisons in this study.

Table 1

Comparison between the normalized thermal resistance (Rstructure/Rinitial) of the optimized topology structure and the reference structures of the studied heat guides shown in Fig. 2 

DesignThermally optimized structureThermomechanically optimized structureSolid reference structureInfill reference structure
a1.311.431.892.85
b1.341.491.872.83
c1.631.781.932.99
d1.551.821.942.89
e1.291.372.012.67
DesignThermally optimized structureThermomechanically optimized structureSolid reference structureInfill reference structure
a1.311.431.892.85
b1.341.491.872.83
c1.631.781.932.99
d1.551.821.942.89
e1.291.372.012.67

The optimization process results in a heat guiding structure with maximum stiffness of 6.08 (N/m) in design (e), while its solid reference counterpart showed a stiffness of 4.24 (N/m). By normalizing the stiffness values to the stiffness of the initial structure of design (e) (8.95 (N/m), the optimized heat guiding structure demonstrates a 43% improvement in the normalized stiffness. Using this optimization process, various electronic packaging designs with numerous constraints can be studied.

3 Manufacturing and Characterization of Three-Dimensional Heat-Guiding Structures

To experimentally validate the concept of heat guiding for packaging applications and characterize their thermal properties, the heat guiding structure of design (e) was benchmarked for manufacturing using a selective laser melting (SLM) technique. During the 3D manufacturing process, the structures were fabricated corresponding to their input CAD files by selectively melting and consolidating thin layers of SS316 L metal powder using a scanner laser beam. The process parameters were optimized via a broad sweep of experiments at many parameter sets, which is a common procedure for fabricating porosity-free samples [61,62]. The structures were printed at 400 W laser power, 230 mm/s scan speed, 60 μm hatch spacing, and 30 μm layer thickness. Figure 3 demonstrates SS316 L 3D printed samples based on the design (e). As demonstrated in Fig. 2(e), the dimensions of this design with 15 mm width and 30 mm height were slightly larger than designs (a–d). The reason behind this selection was the bottom-up manufacturability of the structures. Mechanical supports were implemented within the structural gaps for the 3D printing process that were later removed manually. A reference characterization sample was fabricated for characterizing the intrinsic thermal properties of the 3D printed SS316 L (Fig. 3(a)). Figure 3(b) demonstrates the fabricated initial structure followed by its corresponding solid reference structure with 50% volume fraction (Fig. 3(c)) and the topology optimized heat guiding structure (Fig. 3(d)) after postprocessing.

Fig. 3
Demonstration of the design (e) 3D printed steel structures after postprocessing and emissivity coating. (a) A reference characterization sample was fabricated for characterizing the properties of additively manufactured steel along with (b) the initial structure (ν=100%), (c) the reference structure with 50% volume of the initial structure (ν=50%), and (d) the optimized heat guiding structure (ν=50%).
Fig. 3
Demonstration of the design (e) 3D printed steel structures after postprocessing and emissivity coating. (a) A reference characterization sample was fabricated for characterizing the properties of additively manufactured steel along with (b) the initial structure (ν=100%), (c) the reference structure with 50% volume of the initial structure (ν=50%), and (d) the optimized heat guiding structure (ν=50%).
Close modal

The postprocessed samples were coated with a high emissivity coating spray (RUST-OLEUM high heat primer) and were thermally characterized using the infrared thermography method [63]. The samples were placed in a JANIS VPF-800 vacuum chamber in a vacuum level below 10−5 Torr (high vacuum range). The surface emissivity of the applied coating was calibrated as 0.810 ± 0.013. The estimated error in the calibrated emissivity results in temperature reading with ± 0.5 °C inaccuracy. Using the calibrated emissivity value, the thermal conductivity of the reference characterization sample was measured as 14.18 ± 1.13 Wm−1K−1 through utilizing the infrared thermography methodology discussed in Farzinazar et al. [63]. This value is in good agreement with previously reported thermal conductivity values for additively manufactured stainless steel 316 L (SS316L) samples [51]. Iceberg Thermal DRIFTIce 0.5 mm thermal pad with thermal conductivity of 13 Wm−1K−1 (provided by the vendor) was used as the Thermal Interface Material (TIM) between the fabricated samples and the heat source/heat sink.

A constant heat flux of 5.1 Wcm−2 was applied to both reference design and optimized heat guiding structures which have a similar volume fraction of 50%. A reference quartz disk with known thermal conductivity of 1.38 Wm−1K−1, a diameter of 25 mm, and a thickness of 6.33 mm was used for calibrating the incoming heat flux from the heater (Fig. 4(a)). The actual experimental setup was oriented such that the heater was placed on the top and the heat sink at the bottom. However, since in most packaging scenarios the heat sink is placed on top of the package, the experimental setup image was flipped (Fig. 4(a)) for a direct comparison with a practical application design. The IR images show the temperature contour across the samples and the quartz disk as a result of the applied heat flux (Fig. 4(c)). The measurement of incoming heat flux and temperature gradient in both samples resulted in thermal resistance of 171.1 ± 12.8 KW−1 and 121.1 ± 9.3 KW−1 in the reference and optimized heat guiding structures, respectively.

Fig. 4
(a) Demonstration of the IR thermography experimental setup for the reference structure (top) and the optimized heat guiding structure (bottom). (b) Numerical analysis of the reference sample (top) and the optimized heat guiding structure (bottom) with the measured boundary conditions of heat flux of 5.1 Wcm−2 on the top and a constant temperature of 33 °C in the copper sample holder. The inset figures represent the direction of heat flow using the red arrows. While the reference structure restrictions cause the heat flux concentration in the narrow regions, the optimized heat guiding structure has an optimized distribution of material to prevent the concentration of heat flux. c) IR images demonstrating the temperature contour in the reference structure (top) and the optimized heat guiding structure (bottom). The scale bar in inset images represents 5 mm.
Fig. 4
(a) Demonstration of the IR thermography experimental setup for the reference structure (top) and the optimized heat guiding structure (bottom). (b) Numerical analysis of the reference sample (top) and the optimized heat guiding structure (bottom) with the measured boundary conditions of heat flux of 5.1 Wcm−2 on the top and a constant temperature of 33 °C in the copper sample holder. The inset figures represent the direction of heat flow using the red arrows. While the reference structure restrictions cause the heat flux concentration in the narrow regions, the optimized heat guiding structure has an optimized distribution of material to prevent the concentration of heat flux. c) IR images demonstrating the temperature contour in the reference structure (top) and the optimized heat guiding structure (bottom). The scale bar in inset images represents 5 mm.
Close modal

As demonstrated in Fig. 4(c), the maximum temperature in the reference structure is 71.3 ± 0.8 °C whereas, in the optimized heat guiding structure with a similar volume fraction, this value reaches 61.1 ± 0.6 °C. Numerical simulations with the measured heat flux value of 5.1 Wcm−2 and a constant temperature of 33 °C in the copper anchored heat sink was performed for both reference and the optimized heat guiding structures (Fig. 4(b)). The thermal contact resistance is a crucial parameter for determining the maximum and minimum temperatures across the fabricated samples. We characterized the thermal contact resistance on the sample/TIM interfaces on the heat source and the heat sink sides using the calibrated heat flux values and the temperature gradient across the sample/TIM interface. On the heat sink side, the values were 1.5 ± 0.3 × 10−4 km2/W and 1.3 ± 0.2 × 10−4 km2/W in the reference and the optimized structure, , respectively. On the heat source side, the thermal contact resistance was estimated as 2.8 ± 0.6 × 10−4 km2/W and 2.6 ± 0.5 × 10−4 km2/W in the reference and the optimized structure, , respectively. The thermal contact resistances are smaller on the heat sink side. The main reasons are (1) a higher apparent contact pressure (due to the weight of the system) and (2) a smoother surface quality (due to better polishing). On the heat sink interface, the corresponding estimated thermal contact resistances were used in the simulations. On the heat source side, the contact resistance values were adjusted in the simulation using a parametric study to 2.6 × 10−4 km2/W in the reference structure and 2.2 × 10−4 km2/W in the optimized structure for meeting the measured maximum temperature values. The difference between the estimated average thermal contact resistance values and the simulation results can be attributed to the uncertainties in temperature measurement, heat flux calibration, and the resolution of the Infrared camera—25 μm/pixel—for thermal interface investigations. It is worth mentioning that the reported thermal contact values are consistent with the thermal contact resistance values of elastomer-like gap pads reported in the literature [64].

4 Discussion

Advanced 2.5D and 3D packages could suffer from nontrivial thermal management challenges such as thermal crosstalk and high hot spot heat fluxes that are in the order of 1 kWcm−2. The associated packaging structures may benefit from novel capabilities to guide heat in unaligned configurations of heat source and heat sink through topology optimization. As demonstrated in Table 1, the thermomechanical topology optimization of heat guiding structures results in a degradation in performance as a tradeoff for solving for mechanical objective function in conjunction with the thermal objective function. However, both thermal topology optimization and thermomechanical topology optimization result in structures with a smaller normalized thermal resistance (lower average temperature) compared to their reference structure counterpart at an identical volume fraction. Lower thermal resistance in the topology optimized structures has been confirmed by our experimental and numerical investigations on the design (e) heat guiding structure. On average, the optimized heat guiding structure demonstrated a 50 KW−1 smaller thermal resistance compared to its reference structure counterpart. This difference in thermal resistance can be translated into a 180 °C higher temperature in the reference structure when the heat flux varies from 10 to 100 Wcm−2. In addition to the thermal-mechanical co-optimization, future research and development in multidisciplinary optimization may provide joint efforts to improve thermal management techniques at multiple levels. The optimized complex heat guiding designs can be achieved through 3D printing techniques, which require more studies to improve their process flow, compatibility, and scalability. Moreover, investigations on the implementation methods are also needed to address the chip integration and to provide package-level solutions.

5 Summary and Conclusions

Recent advancements in electronic packaging require customized thermal management solutions to enable guiding the excessive heat in unaligned configurations of heat source and heat sink along with satisfying the mechanical requirements. We presented a density-based topology optimization method to find thermally conductive and mechanically stable structures for optimal heat guiding in nonconventional packaging scenarios. We investigated serpentine designs with different effective lengths from a heat source to a heat sink that is potentially suited for guiding the heat from the system level hotspot to the heat sink while securing the possible adjacent components with different temperature tolerances. Using density-based thermomechanical topology optimization and assuming 50% volume fraction constraint, we fabricated a prototype of the thermally optimized heat guiding structure—design (e). Our experimental results demonstrated 29% lower thermal resistance compared to its reference counterpart with an identical volume fraction. This finding suggests that topology-optimized heat guiding structures can control the heat flux and can provide a higher thermal performance compared to their reference structure counterparts at an identical volume fraction constraint. Future electronic packaging may benefit from the combined efforts across the advanced optimization process and manufacturing techniques for optimal thermal management solutions. Consequently, the resultant heat guiding structures from topology optimization and 3D manufacturing can potentially be used as a stepping stone toward thermal computing.

Acknowledgment

The authors gratefully acknowledge the support from UCI's IDMI facility—RapidTech, in our additive manufacturing efforts.

Funding Data

  • Samsung Electronics (Funder ID: 10.13039/100004358).

  • National Science Foundation Grant CMMI-1902685 issued through the Mission Directorate (Funder ID: 10.13039/100000001).

Nomenclature

f =

load vector

K =

stiffness matrix

p =

penalization factor

q¨ =

heat flux

u =

displacement vector

θ =

nodal temperature vector

κ =

thermal conductivity matrix

ν =

volume fraction

ρ =

density

ω =

weight factor

ϵ =

design domain

References

1.
Lee
,
J.
,
Lim
,
J.
, and
Yang
,
P.
,
2015
, “
Ballistic Phonon Transport in Holey Silicon
,” ACS Publications, Washington, DC.
2.
Ren
,
Z.
,
Yu
,
Z.
,
Kim
,
J. C.
, and
Lee
,
J.
,
2019
, “
TSV-Integrated Thermoelectric Cooling by Holey Silicon for Hot Spot Thermal Management
,”
Nanotechnology
,
30
(
3
), p.
035201
.10.1088/1361-6528/aaea3a
3.
Nguyen
,
D. M.
,
Xu
,
H.
,
Zhang
,
Y.
, and
Zhang
,
B.
,
2015
, “
Active Thermal Cloak
,”
Appl. Phys. Lett.
,
107
(
12
), p.
121901
.10.1063/1.4930989
4.
Narayana
,
S.
, and
Sato
,
Y.
,
2012
, “
Heat Flux Manipulation With Engineered Thermal Materials
,”
Phys. Rev. Lett
,
108
(
21
), p. 214303.10.1103/PhysRevLett.108.214303
5.
Xu
,
H.
,
Shi
,
X.
,
Gao
,
F.
,
Sun
,
H.
, and
Zhang
,
B.
,
2014
, “
Ultrathin Three-Dimensional Thermal Cloak
,”
Phys. Rev. Lett
,
112
(
5
), p. 054301.10.1103/PhysRevLett.112.054301
6.
Park
,
G.
,
Kang
,
S.
,
Lee
,
H.
, and
Choi
,
W.
,
2017
, “
Tunable Multifunctional Thermal Metamaterials: Manipulation of Local Heat Flux Via Assembly of Unit-Cell Thermal Shifters
,”
Sci. Rep.
,
7
(
41000
), p.
41000
.10.1038/srep41000
7.
Langer
,
G.
,
Leitgeb
,
M.
,
Nicolics
,
J.
,
Unger
,
M.
,
Hoschopf
,
H.
, and
Wenzl
,
F. P.
,
2014
, “
Advanced Thermal Management Solutions on PCBs for High Power Applications
,”
J. Microelectron. Electron. Packag.
,
11
(
3
), pp.
104
114
.10.4071/imaps.422
8.
Saums
,
D. L.
, and
Hay
,
R. A.
,
2016
, “
Developments for Copper-Graphite Composite Thermal Cores for PCBs for High-Reliability RF Systems Expansion-Matched Substrates
,”
Addit. Pap. Present
, (
HiTEC
), pp.
73
78
.10.4071/2016-HIT EC-73
9.
Sapatnekar
,
S. S.
,
2009
, “
Addressing Thermal and Power Delivery Bottlenecks in 3D Circuits
,”
2009 Asia and South Pacific Design Automation Conference
, Yokohama, Japan, Jan. 19–22, pp.
423
428
.
10.
Bandaru
,
P. R.
,
Vemuri
,
K. P.
,
Canbazoglu
,
F. M.
, and
Kapadia
,
R. S.
,
2015
, “
Layered Thermal Metamaterials for the Directing and Harvesting of Conductive Heat
,”
AIP Adv.
,
5
(
5
), p.
053403
.10.1063/1.4916220
11.
Schittny
,
R.
,
Kadic
,
M.
,
Guenneau
,
S.
, and
Wegener
,
M.
,
2013
, “
Experiments on Transformation Thermodynamics: Molding the Flow of Heat
,”
Phys. Rev. Lett.
,
110
(
19
), p. 195901.10.1103/PhysRevLett.110.195901
12.
Dede
,
E. M.
,
Nomura
,
T.
,
Schmalenberg
,
P.
, and
Lee
,
J. S.
,
2013
, “
Heat Flux Cloaking, Focusing, and Reversal in Ultra-Thin Composites Considering Conduction-Convection Effects Considering Conduction-Convection Effects
,”
Appl. Phys. Lett.
,
103
(
6
), p.
063501
.10.1063/1.4816775
13.
Dede
,
E. M.
,
Zhou
,
F.
,
Schmalenberg
,
P.
, and
Nomura
,
T.
,
2018
, “
Thermal Metamaterials for Heat Flow Control in Electronics
,”
ASME J. Electron. Packag.
,
140
(
1
), p.
010904
.10.1115/1.4039020
14.
Kim
,
J. C.
,
Ren
,
Z.
,
Yuksel
,
A.
,
Dede
,
E. M.
,
Bandaru
,
P. R.
,
Oh
,
D.
, and
Lee
,
J.
,
2021
, “
Recent Advances in Thermal Metamaterials and Their Future Applications for Electronics Packaging
,”
ASME J. Electron. Packag.
,
143
(
1
), p.
010801
.10.1115/1.4047414
15.
J. H
,
L.
,
2019
, “
Recent Advances and Trends in Fan-Out Wafer/Panel-Level Packaging
,”
ASME J. Electron. Packag.
,
141
(
4
), p.
040801
.10.1115/1.4043341
16.
Lau
,
J. H.
,
2019
, “
Recent Advances and Trends in Heterogeneous Integrations
,”
J. Microelectron. Electron. Packag
,
16
(
2
), pp.
45
77
.10.4071/imaps.780287
17.
Dede
,
E. M.
,
Liu
,
Y.
,
Joshi
,
S. N.
,
Zhou
,
F.
,
Lohan
,
D. J.
, and
Shin
,
J. W.
,
2019
, “
Optimal Design of Three- Dimensional Heat Flow Structures for Power Electronics Applications
,”
ASME J. Thermal Sci. Eng. Appl.
,
11
(
2
), p.
021011
.10.1115/1.4041440
18.
Xie
,
Y. M.
, and
Steven
,
G. P.
,
1993
, “
A Simple Evolutionary Procedure for Structural Optimization
,”
Comput. Struct.
,
49
(
5
), pp.
885
896
.10.1016/0045-7949(93)90035-C
19.
Bendsoe
,
M. P.
,
1988
, “
Generating Optimal Topologies in Structural Design Using a Homogenization Method
,”
Comput. Methods Appl. Mech. Eng.
,
71
, pp.
197
224
.10.1016/0045-7825(88)90086-2
20.
Sethian
,
J. A.
, and
Wiegmann
,
A.
,
2000
, “
Structural Boundary Design Via Level Set and Immersed Interface Methods
,”
J. Comput. Phys.
,
163
(
2
), pp.
489
528
.10.1006/jcph.2000.6581
21.
Wang
,
M. Y.
,
Wang
,
X.
, and
Guo
,
D.
,
2003
, “
A Level Set Method for Structural Topology Optimization
,”
Comput. Methods Appl. Mech. Eng.
,
192
(
1–2
), pp.
227
246
.10.1016/S0045-7825(02)00559-5
22.
Xie
,
Y. M.
,
Yang
,
X.
,
Shen
,
J.
,
Yan
,
X.
,
Ghaedizadeh
,
A.
,
Rong
,
J.
,
Huang
,
X.
, and
Zhou
,
S.
,
2014
, “
Designing Orthotropic Materials for Negative or Zero Compressibility
,”
Int. J. Solids Struct.
,
51
(
23–24
), pp.
4038
4051
.10.1016/j.ijsolstr.2014.07.024
23.
Sigmund
,
O.
,
1995
, “
Tailoring Materials With Prescribed Elastic Properties
,”
Mech. Mater.
,
20
(
4
), pp.
351
368
.10.1016/0167-6636(94)00069-7
24.
Wang
,
F.
,
Sigmund
,
O.
, and
Jensen
,
J. S.
,
2014
, “
Design of Materials With Prescribed Nonlinear Properties
,”
J. Mech. Phys. Solids
,
69
, pp.
156
174
.10.1016/j.jmps.2014.05.003
25.
Iga
,
A.
,
Nishiwaki
,
S.
,
Izui
,
K.
, and
Yoshimura
,
M.
,
2009
, “
Topology Optimization for Thermal Conductors Considering Design-Dependent Effects, Including Heat Conduction and Convection
,”
Int. J. Heat Mass Transfer
,
52
(
11–12
), pp.
2721
2732
.10.1016/j.ijheatmasstransfer.2008.12.013
26.
Zhang
,
J.
,
Wang
,
S.
,
Zhou
,
G.
,
Gong
,
S.
, and
Yin
,
S.
,
2020
, “
Topology Optimization of Thermal Structure for Isotropic and Anisotropic Materials Using the Element-Free Galerkin Method
,”
Eng. Optim.
,
52
(
7
), pp.
1097
1118
.10.1080/0305215X.2019.1636979
27.
Lundgaard
,
C.
,
Engelbrecht
,
K.
, and
Sigmund
,
O.
,
2019
, “
A Density-Based Topology Optimization Methodology for Thermal Energy Storage Systems
,”
Struct. Multidiscip. Optim.
,
60
(
6
), pp.
2189
2204
.10.1007/s00158-019-02375-8
28.
Zhou
,
M.
,
Alexandersen
,
J.
,
Sigmund
,
O.
, and
Pedersen
,
C. B. W.
,
2016
, “
Industrial Application of Topology Optimization for Combined Conductive and Convective Heat Transfer Problems
,”
Struct. Multidiscip. Optim.
,
54
(
4
), pp.
1045
1060
.10.1007/s00158-016-1433-2
29.
Haertel
,
K.
,
Hendrik
,
J.
, and
Stefanov
,
B.
,
2015
, “
Topology Optimization of an Actively Cooled Electronics Section for Downhole Tools
,”
Proceedings COMSOL Conference
, Grenoble, France, Oct.
14
16
.https://www.comsol.com/paper/topology-optimization-of-an-actively-cooled-electronics-section-for-downhole-too-28731
30.
Matsumori
,
T.
,
Kawamoto
,
A.
, and
Kondoh
,
T.
,
2019
, “
Topology Optimization for Thermal Stress Reduction in Power Semiconductor Module
,”
Struct. Multidiscip. Optim.
,
60
(
6
), pp.
2615
2620
.10.1007/s00158-019-02341-4
31.
Chen
,
C.
,
Weng
,
Y.
, and
Subbarayan
,
G.
,
2016
, “
Topology Optimization for Efficient Heat Removal in 3D Packages
,”
IEEE ITHERM Conference
, Las Vegas, NV, May 31–June 3.10.1109/ITHERM.2016.7517556
32.
Martínez-Maradiaga
,
D.
,
Damonte
,
A.
,
Manzo
,
A.
,
Haertel
,
J. H. K.
, and
Engelbrecht
,
K.
,
2019
, “
Design and Testing of Topology Optimized Heat Sinks for a Tablet
,”
Int. J. Heat Mass Transfer
,
142
(
118429
), p.
118429
.10.1016/j.ijheatmasstransfer.2019.07.079
33.
Chen
,
Y.
,
Zhou
,
S.
, and
Li
,
Q.
,
2010
, “
Multiobjective Topology Optimization for Finite Periodic Structures
,”
Comput. Struct.
,
88
(
11–12
), pp.
806
811
.10.1016/j.compstruc.2009.10.003
34.
Dede
,
E. M.
,
2021
, “
Simulation-Based Design and Optimization for Future Mobility Electronics Systems
,” 22nd International Conference on Thermal, Mechanical and Multi-Physics Simulation and Experiments in Microelectronics and Microsystems (
EuroSimE
), St. Julian, Malta, Apr.
19
21
.10.1109/EuroSimE52062.2021.9410864
35.
Takezawa
,
A.
,
Yoon
,
G. H.
,
Jeong
,
S. H.
,
Kobashi
,
M.
, and
Kitamura
,
M.
,
2014
, “
Structural Topology Optimization With Strength and Heat Conduction Constraints
,”
Comput. Methods Appl. Mech. Eng.
,
276
, pp.
341
361
.10.1016/j.cma.2014.04.003
36.
Zhu
,
X.
,
Zhao
,
C.
,
Wang
,
X.
,
Zhou
,
Y.
,
Hu
,
P.
, and
Ma
,
Z. D.
,
2019
, “
Temperature-Constrained Topology Optimization of Thermo-Mechanical Coupled Problems
,”
Eng. Optim.
,
51
(
10
), pp.
1687
1709
.10.1080/0305215X.2018.1554065
37.
Giraldo-Londoño
,
O.
,
Mirabella
,
L.
,
Dalloro
,
L.
, and
Paulino
,
G. H.
,
2020
, “
Multi-Material Thermomechanical Topology Optimization With Applications to Additive Manufacturing: Design of Main Composite Part and Its Support Structure
,”
Comput. Methods Appl. Mech. Eng.
,
363
, p.
112812
.10.1016/j.cma.2019.112812
38.
Wu
,
C.
,
Fang
,
J.
, and
Li
,
Q.
,
2019
, “
Multi-Material Topology Optimization for Thermal Buckling Criteria
,”
Comput. Methods Appl. Mech. Eng.
,
346
, pp.
1136
1155
.10.1016/j.cma.2018.08.015
39.
Du
,
Y.
,
Luo
,
Z.
,
Tian
,
Q.
, and
Chen
,
L.
,
2009
, “
Topology Optimization for Thermo-Mechanical Compliant Actuators Using Mesh-Free Methods
,”
Eng. Optim.
,
41
(
8
), pp.
753
772
.10.1080/03052150902834989
40.
Wu
,
T.
,
Najmon
,
J. C.
, and
Tovar
,
A.
,
2019
, “
Thermomechanical Topology Optimization of Lattice Heat Transfer Structure Including Natural Convection and Design Dependent Heat Source
,”
ASME
Paper No. DETC2019-97628.10.1115/DETC2019-97628
41.
Bruggi
,
M.
, and
Taliercio
,
A.
,
2013
, “
Design of Masonry Blocks With Enhanced Thermomechanical Performances by Topology Optimization
,”
Constr. Build. Mater.
,
48
, pp.
424
433
.10.1016/j.conbuildmat.2013.07.023
42.
Li
,
L.
,
Du
,
Z.
, and
Kim
,
H. A.
,
2020
, “
Design of Architected Materials for Thermoelastic Macrostructures Using Level Set Method
,”
JOM
,
72
(
4
), pp.
1734
1744
.10.1007/s11837-020-04046-2
43.
De Kruijf
,
N.
,
Zhou
,
S.
,
Li
,
Q.
, and
Mai
,
Y.-W.
,
2007
, “
Topological Design of Structures and Composite Materials With Multiobjectives
,”
Int. J. Solids Struct.
,
44
(
22–23
), pp.
7092
7109
.10.1016/j.ijsolstr.2007.03.028
44.
Li
,
L.
, and
Kim
,
H. A.
,
2020
, “
Multiscale Topology Optimization of Thermoelastic Structures Using Level Set Method
,”
AIAA
Paper No. 2020-0890.10.2514/6.2020-0890
45.
Dede
,
E. M.
,
2010
, “
Simulation and Optimization of Heat Flow Via Anisotropic Material Thermal Conductivity
,”
Comput. Mater. Sci.
,
50
(
2
), pp.
510
515
.10.1016/j.commatsci.2010.09.012
46.
Dede
,
E. M.
,
Nomura
,
T.
, and
Lee
,
J.
,
2014
, “
Thermal-Composite Design Optimization for Heat Flux Shielding, Focusing, and Reversal
,”
Struct. Multidiscip. Optim.
,
49
(
1
), pp.
59
68
.10.1007/s00158-013-0963-0
47.
Fujii
,
G.
, and
Akimoto
,
Y.
,
2019
, “
Optimizing the Structural Topology of Bifunctional Invisible Cloak Manipulating Heat Flux and Direct Current
,”
Appl. Phys. Lett.
,
115
(
17
), p.
174101
.10.1063/1.5123908
48.
Lazarov
,
B. S.
, and
Sigmund
,
O.
,
2011
, “
Filters in Topology Optimization Based on Helmholtz-Type Differential Equations
,”
Int. J. Numer. Methods Eng.
,
86
(
6
), pp.
765
781
.10.1002/nme.3072
49.
Sanders
,
E. D.
,
Aguiló
,
M. A.
, and
Paulino
,
G. H.
,
2018
, “
Multi-Material Continuum Topology Optimization With Arbitrary Volume and Mass Constraints
,”
Comput. Methods Appl. Mech. Eng.
,
340
, pp.
798
823
.10.1016/j.cma.2018.01.032
50.
Zhang
,
B.
,
Dembinski
,
L.
, and
Coddet
,
C.
,
2013
, “
The Study of the Laser Parameters and Environment Variables Effect on Mechanical Properties of High Compact Parts Elaborated by Selective Laser Melting 316 L Powder
,”
Mater. Sci. Eng. A
,
584
, pp.
21
31
.10.1016/j.msea.2013.06.055
51.
Yakout
,
M.
,
Elbestawi
,
M. A.
,
Veldhuis
,
S. C.
, and
Nangle-Smith
,
S.
,
2020
, “
Influence of Thermal Properties on Residual Stresses in SLM of Aerospace Alloys
,”
Rapid Prototyp.
,
26
(
1
), pp.
213
222
.10.1108/RPJ-03-2019-0065
52.
Li
,
D.
,
Zhang
,
X.
,
Guan
,
Y.
, and
Zhan
,
J.
,
2010
, “
Topology Optimization of Thermo-Mechanical Continuum Structure
,”
IEEE/ASME International Conference on Advanced Intelligent Mechatronics
, Montreal, QC, Canada, July 6–9, pp.
403
408
.10.1109/AIM.2010.5695845
53.
Svanberg
,
K.
,
1987
, “
The Method of Moving Asymptotes-a New Method for Structural Optimization
,”
Int. J. Numer. Methods Eng.
,
24
(
2
), pp.
359
373
.10.1002/nme.1620240207
54.
Comotti
,
C.
,
Regazzoni
,
D.
,
Rizzi
,
C.
, and
Vitali
,
A.
,
2017
, “
Additive Manufacturing to Advance Functional Design: An Application in the Medical Field
,”
ASME J. Comput. Inf. Sci. Eng.
,
17
(
3
), p.
031006
.10.1115/1.4033994
55.
Ivorra-Martinez
,
J.
,
Quiles-Carrillo
,
L.
,
Lascano
,
D.
,
Ferrandiz
,
S.
, and
Boronat
,
T.
,
2020
, “
Effect of Infill Parameters on Mechanical Properties in Additive Manufacturing
,”
Dynamics
,
95
(
1
), pp.
412
417
.10.6036/9674
56.
Steuben
,
J. C.
,
Iliopoulos
,
A. P.
, and
Michopoulos
,
J. G.
,
2016
, “
Implicit Slicing for Functionally Tailored Additive Manufacturing
,”
Comput. Aided Des.
,
77
, pp.
107
119
.10.1016/j.cad.2016.04.003
57.
Svanberg
,
K.
,
2002
, “
A Class of Globally Convergent Optimization Methods Based on Conservative Convex Separable Approximations
,”
Soc. Ind. Appl. Math.
,
12
(
2
), pp.
555
573
.10.1137/S1052623499362822
58.
Bendsøe
,
M. P.
, and
Sigmund
,
O.
,
2004
,
Topology Optimization Theory, Methods and Applications
, Springer, Berlin.
59.
Sohel Murshed
,
S. M.
, and
Nieto de Castro
,
C. A.
,
2017
, “
A Critical Review of Traditional and Emerging Techniques and Fluids for Electronics Cooling
,”
Renew. Sustain. Energy Rev.
,
78
(
May
), pp.
821
833
.10.1016/j.rser.2017.04.112
60.
Jung
,
M.
,
Pan
,
D. Z.
, and
Lim
,
S. K.
,
2012
, “
Chip/Package Co-Analysis of Thermo-Mechanical Stress and Reliability in TSV-Based 3D ICs
,”
Proceedings - Design Automation Conference
, San Francisco, CA, June 3–7, pp.
317
326
.10.1145/2228360.2228419
61.
Rashid
,
R.
,
Masood
,
S. H.
,
Ruan
,
D.
,
Palanisamy
,
S.
,
Rashid
,
R. A. R.
, and
Brandt
,
M.
,
2017
, “
Effect of Scan Strategy on Density and Metallurgical Properties of 17-4PH Parts Printed by Selective Laser Melting (SLM)
,”
J. Mater. Process. Tech.
,
249
, pp.
502
511
.10.1016/j.jmatprotec.2017.06.023
62.
Kasperovich
,
G.
,
Haubrich
,
J.
,
Gussone
,
J.
, and
Requena
,
G.
,
2016
, “
Correlation Between Porosity and Processing Parameters in TiAl6V4 Produced by Selective Laser Melting
,”
JMADE
,
105
, pp.
160
170
.10.1016/j.matdes.2016.05.070
63.
Farzinazar
,
S.
,
Schaedler
,
T.
,
Valdevit
,
L.
, and
Lee
,
J.
,
2019
, “
Thermal Transport in Hollow Metallic Microlattices
,”
APL Mater.
,
7
(
10
), p.
101108
.10.1063/1.5114955
64.
Swamy
,
M. C. K.
, and,
Satyanarayan
,
2019
, “
A Review of the Performance and Characterization of Conventional and Promising Thermal Interface Materials for Electronic Package Applications
,”
J. Electron. Mater.
,
48
(
12
), pp.
7623
7634
.10.1007/s11664-019-07623-7